A Force-Balanced Fiber Retardation Model to Predict Fiber-Matrix-Separation during Polymer Processing

: The use of discontinuous ﬁber reinforced composites in injection and compression molding faces a number of challenges regarding process-induced changes in microstructure, which have a signiﬁcant inﬂuence on the mechanical properties of the ﬁnal component. The changes in ﬁnal microstructure are caused by complex ﬁber movements, such as ﬁber orientation, attrition and accumulation during ﬂow. While there are existing phenomenological prediction models for both ﬁber orientation and attrition, the prediction of ﬁber accumulation due to ﬁber-matrix separation is currently only possible with a complex mechanistic particle simulation, which is not applicable in industrial simulations. A simpliﬁed phenomenological model, the ﬁber retardation model (FRM), for the prediction of ﬁber-matrix separation in commercially available software tools is presented in this paper. The model applies a force balance onto an interacting two phase ﬂow of polymer melt and ﬁber phase and applies a retardation factor K to calculate the slowing and accumulation of the ﬁber phase. The general model is successfully applied to a simple compression molding simulation.


Introduction
Discontinuous fiber reinforced polymer composites in compression molding have been applied successful in a variety of industries due to their cost-efficient processing of complex parts and their great 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][2][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][2][3][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.  (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   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][9][10][11][12][13][14].   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][9][10][11][12][13][14]. 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][9][10][11][12][13][14].  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   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 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.

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. account for these effects during manufacturing, comparable to existing fiber orientation and fiber attrition models.

General Fiber Retardation Model
The general interactions implemented in the complex Mechanistic Model are applied to a twophase 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. 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. The ratio of the interacting forces to the hydrodynamic forces can be described with the retardation factor Κ. For minimal fiber interaction forces and high hydrodynamic forces Κ→0 and no fiber matrix separation occurs and the fiber phase follows the hydrodynamic forces. With increasing fiber interaction or low hydrodynamic forces, Κ increases, leading to maximum fiber-matrix 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 F hydro , and the counter-acting forces of the fiber-network (or inter-fiber) forces F FN and the fiber interactions with the mold wall F int . 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. account for these effects during manufacturing, comparable to existing fiber orientation and fiber attrition models.

General Fiber Retardation Model
The general interactions implemented in the complex Mechanistic Model are applied to a twophase 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. 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. The ratio of the interacting forces to the hydrodynamic forces can be described with the retardation factor Κ. For minimal fiber interaction forces and high hydrodynamic forces Κ→0 and no fiber matrix separation occurs and the fiber phase follows the hydrodynamic forces. With increasing fiber interaction or low hydrodynamic forces, Κ increases, leading to maximum fiber-matrix 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.
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.
The detailed information on the force balance calculation is described in detail in the following section.

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. dp dx = − η P u 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φ 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 E f , the total fiber length per unit area ρ geom and the cross section area Ω.
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.
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 .
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. 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. The overall fiber interaction force is then calculated by the summation of fiber interaction forces of all individual fibers.

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 h0, x0, constant width B and a constant glass fiber volume fraction ϕ is compressed with the constant closing speed a, as shown in Figure 8. The resulting gap height and flow velocity development over time during compression molding are presented in Figure 9. The overall fiber interaction force is then calculated by the summation of fiber interaction forces of all individual fibers.

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. 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. The overall fiber interaction force is then calculated by the summation of fiber interaction forces of all individual fibers.

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 h0, x0, constant width B and a constant glass fiber volume fraction ϕ is compressed with the constant closing speed a, as shown in Figure 8. The resulting gap height and flow velocity development over time during compression molding are presented in Figure 9. During constant compression, the velocity of the melt front advances over time.
The resulting gap height and flow velocity development over time during compression molding are presented in Figure 9. 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. 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. 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. 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. 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. 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 Κ increases significantly with 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. 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. 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 Κ increases significantly with 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. 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.

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 great 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.