**1. Introduction**

Discontinuous fiber reinforced polymer composites in compression molding have been applied successful in a variety of industries due to their cost-e fficient processing of complex parts and their grea<sup>t</sup> mechanical properties to weight ratio. Their availability in a variety of polymer resins (thermoplastics and thermosets) as well as their reinforcement fibers (glass, carbon, natural, etc.) satisfy a wide range of requirements for industrial application, e.g., within the automotive and aviation industry. The design of components with fiber reinforced composites faces the challenges of anisotropic and fluctuating properties inside a component due to changes in the local microstructure during processing [1]. These changes are caused by complex fiber behavior such as fiber orientation, attrition and accumulation during the polymer flow. Simulative models are generally applied to predict these changes and improve component design, while the prediction of fiber orientation and attrition has been the main focus in the past [1–3].

With longer reinforcement fibers as traditionally found in compression molding of sheet molding compound (SMC) or long fiber thermoplastics (LFT), increasing fiber interactions lead to fiber matrix phase separation during polymer flow and hence to fluctuations in fiber content inside complex components [1–4]. Early experiments by Schmachtenberg et al. [5] showed increasing fiber-matrix separation during the compression molding with thermoset sheet molding compounds in simple plate geometries, as shown in Figure 1. Schmachtenberg's experiments display a relative fiber-free flow front followed by a section of increased fiber content.

**Figure 1.** Relative change of fiber content distribution inside a compression molded sheet molding compound (SMC) plate from the melt front (1) to the initial charge position (17) with different flow length 's' [5].

Recent experiments by Kuhn et al. [6] show a similar behavior with glass mat reinforced thermoplastics (GMT), as shown in Figure 2. Kuhn et al. further showed the increasing fiber-matrix separation in complex rib regions.

**Figure 2.** Change of fiber content inside a plate region with glass mat reinforced thermoplastics (GMT) GF30 and long fiber GMT GF30.

The effect of fiber-matrix separation leads to significant decreases in mechanical properties and eventually to unexpected component failure. There are simulative models available by Morris and Boulay [7] to predict local changes of fiber content through shear migration in the normal direction of a polymer flow, which are unable to predict global changes in different component regions. Currently, these global effects can only be predicted with complex particle-level simulations, such as the Mechanistic Model, as displayed in Figure 3, where the fibers are simulated as beam elements inside the polymer flow, capable of displaying complex mechanical interactions with the polymer melt, other fibers and the mold walls [8–14].

**Figure 3.** Mechanistic model approach with fibers displayed as beam elements inside the polymer flow [6].

Peréz et al. applied the Mechanistic Model to compression molding with long fiber reinforced composites with fiber distribution results generally complying with experiments, as shown in Figure 4 [9].

**Figure 4.** Mechanistic model predictions of fiber volume fraction over distance from the part center 'X' during biaxial flow in compression molding with long fiber reinforced plastics [9].

While the Mechanistic Model is capable to predict the effect of fiber matrix-separation very accurately, it is currently inapplicable to most industrial applications due to the simulative complexity, with extensive pre-and post-processing and long calculation time even with high computational

capacity. Therefore, a simplified model for fiber-matrix separation is necessary to account for these effects during manufacturing, comparable to existing fiber orientation and fiber attrition models.

#### **2. General Fiber Retardation Model**

The general interactions implemented in the complex Mechanistic Model are applied to a two-phase flow model to generate a simplified fiber retardation model (FRM). Initially, a two-phase flow consisting of the polymer phase and the fiber phase is generated. In order to allow relative motion between the fiber and polymer melt phase, each phase respectively travels at its own velocity, as shown in Figure 5. For the application of the model for commercial software, it is assumed that the velocity of the melt phase *u* is provided.

**Figure 5.** Two-phase model for relative motion between fiber and polymer melt phase.

Regarding the fiber phase in Figure 6, it is obvious that any movements is governed by the surrounding polymer melt. The movement of the fiber phase is based on the forces onto the fiber phase, the hydrodynamic melt forces during polymer flow Fhydro, and the counter-acting forces of the fiber-network (or inter-fiber) forces FFN and the fiber interactions with the mold wall Fint. According to experimental results and Mechanistic Model simulations, the hydrodynamic forces lead to the fiber phase following the melt flow, while increasing fiber interactions lead to the slowing, or retardation, of the fiber phase.

**Figure 6.** Fiber phase force balance.

The ratio of the interacting forces to the hydrodynamic forces can be described with the retardation factor K. For minimal fiber interaction forces and high hydrodynamic forces K → 0 and no fiber matrix separation occurs and the fiber phase follows the hydrodynamic forces. With increasing fiber interaction or low hydrodynamic forces, K increases, leading to maximum fiber-matrix separation for values of K ≥ 1, which would lead to a full halt of the fiber phase movement while the polymer melt phase advances.

$$\mathbb{K} = \frac{F\_{FN} + F\_{INT}}{F\_{hydro}}$$

This implies that with negligible inertia of the fiber bed inside a stationary flow, the retardation factor K (for 0 < K < 1) can be used to calculate the relative velocity of the fiber bed *v* with a given velocity of the melt phase *u*.

$$v = (1 - \mathbb{K}) \times u$$

The detailed information on the force balance calculation is described in detail in the following section.

#### **3. Force Balance Calculation**

The three main forces during processing as explained earlier are the hydrodynamic forces, leading to the initial fiber phase movement and the counteracting fiber forces–fiber interaction and fiber network forces.

The hydrodynamic forces are calculated using Stokes Equation [15] regarding the pressure gradient *dp*/*dx* through a medium with porosity *P* based on the fluid viscosity η and the fluid velocity *u*. This further complies with Darcys law.

$$\frac{dp}{d\mathbf{x}} = -\frac{\eta}{P}\boldsymbol{w}$$

Bhakarev and Tucker [16] describe the porosity of a fiber bed *P* with the fiber volume content φ and the fiber diameter *d*.

$$P = 0.00025 \phi^{2.4} d^2$$

The fiber network forces have to be applied to the respective type of fiber network inside the composite. One solution for nonwoven felts, as generally used in glass mat reinforced thermoplastics (GMT), is the Cox model [17]. Cox describes the Youngs modulus of the fiber network with the single fiber modulus Ef, the total fiber length per unit area ρ*geom* and the cross section area Ω.

$$E(\varepsilon = \frac{1}{3} E\_f \Omega \rho\_{\text{geom}})$$

The fiber network forces are important during the initial stages of compression molding, where the fiber network experiences the first deformation and fibers are pulled apart. If the hydrodynamic forces are significantly lower than the network forces, e.g., through low viscosity levels, the polymer seeps through the fiber bed without creating any fiber movement. This sponge-like phenomena is called "bleeding out". Traditionally, the general type of fiber network is matched to the viscosity of the polymer material to avoid this effect, hence fiber networks with sheet molding compounds are generally weaker than the needled felts applied in GMTs.

The fiber interaction forces with the mold walls are the main drivers for fiber-matrix separation as observed in experiments and Mechanistic Model simulations. Londono et al. [18] applied the fiber interaction forces inside a gap with the general bending forces of a fiber. The force required for fiber bending is calculated with the fiber bending constant *C*, deflection of the fiber δ, the Young's modulus *E* and momentum of inertia *I* and the free bending distance *L*.

$$F\_{f\text{hor}} = C \frac{\delta EI}{L^3}.$$

During molding, the fiber bending only occurs when the individual fiber is in interaction with the mold walls. With a given fiber distribution function at a specific gap height *h* and a fiber length *l*, the amount of fibers in interaction can be calculated with a critical angle of interaction ψ*crit*.

$$
\psi\_{crit}(h) = \sin^{-1} \frac{h(t, \propto)}{l}
$$

The geometrical bending of the fiber and the resulting fiber force is calculated on the angle and length of the fiber in Figure 7. With given fiber length *l* at an orientation angle ψ, a straight fiber would technically extend through the mold wall. The bending deflection of the fiber δ is then necessary so the fiber geometrically fits inside the gap.

**Figure 7.** Fiber bending inside a gap during compression molding.

The overall fiber interaction force is then calculated by the summation of fiber interaction forces of all individual fibers.

#### **4. Compression Molding Example**

The fiber retardation model is applied in a simple compression molding simulation with a unidirectional flow front and the use of glass fibers, where a compression charge with initial dimensions *h*0, *x*0, constant width *B* and a constant glass fiber volume fraction φ is compressed with the constant closing speed *a*, as shown in Figure 8.

**Figure 8.** Compression molding setup.

During constant compression, the velocity of the melt front advances over time.

$$\mathbf{x}'(t) = \boldsymbol{\mu}(t) = \frac{a \times h0 \times \mathbf{x}0}{\left(h0 - at\right)^2}$$

The resulting gap height and flow velocity development over time during compression molding are presented in Figure 9.

**Figure 9.** Gap height and flow velocity during compression molding.

The fiber orientation in commercially available molding software is calculated at every time step and therefore considered provided. For our simple demonstration, three fiber orientation distribution profiles are chosen from literature, as displayed in Figure 10, which were taken from final compression molding experiments with different initial mold coverage of 33%, 50% and 67% [19,20]. In addition, a fully random orientation distribution was added. The three orientation distribution functions are considered constant during our simulation, which of course does not comply with the realistic fiber orientation development during molding. In this first trial, the simplification is applicable.

**Figure 10.** Fiber orientation distribution frequency profiles over fiber orientation angles in compression molded plates with 33%, 50% and 67% initial mold coverage.

The different fiber distribution functions are then applied to calculate the fiber interaction forces at every gap height *h* for all fibers at the specific orientation angle ψ. Fiber network forces are not considered during this simulation and will be focused on in a later publication. The hydrodynamic forces are calculated accordingly. All resulting forces are displayed in Figure 11.

**Figure 11.** Development of fiber forces with different orientation distributions (blue) and hydrodynamic forces (red) during compression molding.

The development of the retardation factor K, the ratio of fiber interaction forces to hydrodynamic forces during compression molding is shown in Figure 12. It is observed that K initially remains at zero until the first fiber interactions with the mold wall occur. K increases significantly with time until the increasing molding speed leads to significantly higher hydrodynamic forces. This implies that a fiber phase of constant porosity *P*, would first follow the polymer melt velocity, then follow the melt phase up to 16% relatively slower than the governing polymer melt and then catch up with the melt velocity again towards the end of molding.

**Figure 12.** Fiber retardation factor K during the compression molding of a random fiber network φ = 0.4.

This generally proves that the fiber retardation model is able to predict a slowing of the fiber phase caused by increasing fiber interactions during compression molding.

Additional evaluations regarding the influence of different properties have been conducted. Values which are only applied in either force calculation, hydrodynamic and fiber interaction, have a clear influence of either increasing or decreasing the forces accordingly. Other factors, which are implemented in both calculations are further evaluated in details. The influence of different fiber contents which have an influence on fiber interactions as well as hydrodynamic forces, are displayed in Figure 12. It is shown that in Figure 13, the fiber retardation factor K increases significantly with lower fiber contents and decreases with higher fiber contents. Simulations with a fiber content of 10 Vol.−% even display a K ≥ 1, which implies a complete halt of the fiber phase inside the melt flow.

**Figure 13.** Development of retardation factor K in compression molding simulations with different fiber contents φ at random orientation.

The results on the influence of fiber content shows that the hydrodynamic forces are the leading factor regarding fiber-matrix-separation during compression molding. Furthermore, it shows that higher fiber contents display less fiber-matrix-separation, which generally agrees to earlier results from compression molding experiments with GF30 and GF40 by Kuhn.

#### **5. Discussion and Outlook**

This paper proposes a novel force-balanced approach to predict fiber-matrix-separation behavior during the processing of long fiber reinforced composites. Earlier experimental findings show increasing interaction forces of the fiber phase during flow, counteracting the hydrodynamic forces of the melt phase and eventually leading to a slowing of the fiber phase inside the melt flow. The proposed model takes these fundamental forces acting on the fiber bed into relation, enabling a slowing of the fiber phase if interaction forces increase compared to the hydrodynamic forces of the melt flow.

While the current stage of the model employs a number of generalized model assumptions and simplified mechanics, it demonstrates a general proof of concept for the realistic prediction of the "fiber-phase retardation". The model results generally comply with experimental findings and literature, whereas increasing fiber interaction forces during processing lead to a force counteracting the hydrodynamic forces of the melt. This increase in counteracting force reduces the fiber bed velocity relative to the melt velocity, eventually leading to fiber accumulation inside complex component regions and, in drastic cases, to a complete stop of the fiber bed and a bleeding out of the resin. As mentioned above, the governing factors are the fiber and melt properties, which also complies with literature findings.

While the current state of the fiber retardation model has shown its general suitability to predict fiber-matrix-separation effects, it has not been fully evaluated with experimental studies regarding the choice of simulative boundary conditions and assumed parameters, which is planned for future publications. In addition, the first calculations were applied to a single fiber region with constant fiber content in order to show the slowing of the fiber phase relative to the melt flow. The next steps

regarding the fiber retardation model is to apply it to a model of multiple fiber regions, in which fiber transport and accumulation due to retardation is possible. This will make the application of the fiber retardation model more suitable for industrial application in phenomenological models in compression molding.

The authors are aware of the current simplicity of the model, with specifically chosen simple model assumptions and simplified boundary conditions. It is clear that extensive work both on the simulative and experimental side is required to ultimately create an accurate fiber content distribution model for commercial implementation. In this context, the model's current simplicity of employing a simple force-balance between fiber and melt phase is a grea<sup>t</sup> starting point for future improvements. Within the force balance equation, it is effortlessly possible to include more complex models, boundary conditions and formulations to utilize existing process information of for example, the surrounding flow field, local fiber orientation, fiber content and fiber length settings. Hence, it is possible to increase the predictability and efficiency of the force balance model significantly. Comparable to the evolution of fiber orientation models in the last decades, new improvements can be added to the force-balanced fiber retardation model step by step.

**Author Contributions:** Data curation, S.W.; Investigation, C.K.; Writing–original draft, C.K. and S.W.; Writing–review & editing, C.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
