Influence of flow–fiber coupling during mold-filling on the stress field in short-fiber reinforced composites

The anisotropic elastic properties of injection molded composites are fundamentally coupled to the flow of the fiber suspension during mold-filling. Regarding the modeling of mold-filling processes, both a decoupled and a flow–fiber coupled approach are possible. In the latter, the fiber-induced viscous anisotropy is considered in the computation of the flow field. This in turn influences the evolution of the fiber orientation compared to the decoupled case. This study investigates how flow–fiber coupling in mold-filling simulation affects the stress field in the solid composite under load based on the final elastic properties after fluid–solid transition. Furthermore, the effects of Newtonian and non-Newtonian polymer matrix behavior are investigated and compared. The entire process is modeled micromechanically unified based on mean-field homogenization, both for the fiber suspension and for the solid composite. Different numerical stabilization methods of the mold-filling simulation are discussed in detail. Short glass fibers with a typical aspect ratio of 20 and a volume fraction of 20% are considered, embedded in polypropylene matrix material. The results show that the flow–fiber coupling has a large effect on the fiber orientation tensor in the range of over ± 30% with respect to the decoupled simulation. As a consequence, the flow–fiber coupling affects the stress field in the solid composite under load in the range of over ± 10%. In addition, the predictions based on a non-Newtonian modeling of the matrix fluid differ significantly from the Newtonian setup and thus the necessity to consider the shear-thinning behavior is justified in a quantifiable manner.


Motivation and state of the art
Short-fiber reinforced composites are commonly used in lightweight design. The special case of adding fibers into a polymer matrix allows the mass-production of lightweight components with complex shapes while retaining the desired mechanical stiffness and strength. During injection moldingthe manufacturing process through which such components are produced-the fibers flow within the polymer matrix and thereby change their orientation based on the experienced flow conditions. At the same time, a local change of the fiber orientation affects the overall rheological properties of the fluid, which become inhomogeneous anisotropic, and in turn the flow characteristics. This mutual flow-fiber interaction is subject of recent numerical research efforts addressed in the following, focussing on the effects of flow-fiber coupling on the fluid side. In the context of this work, flow-fiber coupling expresses that the evolution of fiber orientation is influenced by the effective anisotropic viscosity of the fiber suspension. Three main points can be identified in this context: The effect of the flow-fiber coupling on the flow field, the influenced fiber orientation evolution due to the changed flow conditions and finally the change of the mechanical properties of the composite, which depend on the local fiber orientation. How the flow-fiber coupling during mold-filling affects the stress state in the solid composite under load after fluid-solid transition is a question that received less attention in literature and is the objective of the present work. Since the correct prediction of the stress field in the manufactured composite is an essential part in engineering applications, the need for a flow-fiber coupled simulation has to be studied. In this context, the question also arises which influence the two-phase simulation of the actual mold filling has both on the fiber orientation and on the local anisotropic properties.
The fiber-induced anisotropic viscosity within the flowfiber coupled approach causes changes of the velocity field during mold-filling. Latz et al. [1] investigated the flow-fiber coupling from a fluid mechanics perspective in a channel flow, a flow around a cylinder and in a contraction flow. Different coupling intensities were considered with the result, that the coupling effects are significant in the channel and in the contraction flow. It is shown that (e.g. for the channel flow) the stronger the coupling, the flatter the velocity profile. This flattening was already reported in the early studies of Altan et al. [2] and Tang and Altan [3]. The importance of accounting for flow-fiber coupling in mold-filling simulations was demonstrated in a low Reynolds number flow through a tapered channel by Krochak et al. [4]. Both isotropic and aligned orientation states were considered at the inlet, showing that the flow changes towards a plug flow because of flow-fiber coupling. Flow-fiber coupled moldfilling simulation based on probability density function was carried out by Mezi et al. [5] in 2D geometries and by Férec et al. [6] in 3D axisymmetric geometries. Both studies conclude that the influence of flow-fiber coupling on the velocity field is significant and the velocity profile flattens, as mentioned before. In addition, the difference between Newtonian and non-Newtonian matrix fluid is addressed by Mezi et al. [5]. Mezi et al. [7] observed a flattening of the velocity profile in the context of flow-fiber coupling for die swell flows. The results show that the swell ratio of the flow is affected by flow fiber coupling. By using a scalar rheologial model, Li and Luyé [8] studied the flow-fiber coupling effects based on mold-filling simulations of a rectangular plate. In addition, three different orientation evolution models were used in their study. They observed, that the coupling effects strongly depend on the chosen orientation evolution model. In this context, the classical Folgar-Tucker model [9,10] showed less coupling effects than the reduced-strain closure model [11]. The aforementioned flattening of the velocity profile also occured. In addition, the study of Li and Luyé [12] shows that results based on flow-fiber coupling improve parameter optimization of fiber orientation models in the core region of the flow. Tseng and Favaloro [13] introduced the informed isotropic viscosity model in order to improve convergence compared to the classic models [14,15]. It is shown for simple and complex geometries, that this model is capable of reproducing experimental findings and that the coupling effects on the flow field are not negligible. The informed isotropic viscosity model was used by Huang and Lai [16] to show that flow-fiber coupling influences the melt front behavior and that the coupling effects are locally different. Lee et al. [17] also used the informed isotropic viscosity model and studied the compression molding process of composite material. Sommer et al. [18] investigated the flow-fiber coupling effects in squeeze flow. The results show that only a flowfiber coupled analysis is capable of predicting the physically correct anisotropic behavior of the suspension. Based on two-phase flow, Witteman et al. [19] implemented the flowfiber coupling approach with a non-Newtonian, temperatureand curing-dependent viscosity model for fiber-laden thermosets. The coupled simulation leads to a better agreement of the pressure field with experimental measurements than the decoupled approach.
The flow-fiber coupling effects on the velocity field addressed above have a direct influence on the fiber orientation evolution. Chung and Kwon [20,21] studied the flow-fiber coupling effects on both flow field and fiber orientation. Their results show that flow-fiber coupling clearly influences the fiber orientation prediction. Latz et al. [1] found that the stronger the coupling, the larger the region of isotropic orientation in the channel center region. Li and Luyé [8] observed, that the fiber orientation evolution is clearly influenced by the changed flow kinematics for the studied mold-filling process of a rectangular plate. Wang and Smith [22,23] investigated the flow-fiber coupling regarding an additive manufacturing process. The fiber orientation evolution referred to the orthotropically closed Folgar-Tucker model and the coupling was realized based on the fourthorder fiber orientation tensor. The results indicate that the fiber-fiber interaction coefficient strongly affects the influence of flow-fiber coupling. Significant differences in the fiber orientation tensor were observed due to flow-fiber coupling. In contrast, Mezi et al. [5,7] and Férec et al. [6] did not observe significant differences regarding the fiber orientation state between the decoupled and the coupled approach in their simulation setups.
Karl et al. [24] studied the effect of flow-fiber coupling on both velocity and fiber orientation. A simple channel flow and a flow over a backward-facing step was considered with a detailed analysis of the aforementioned velocity flattening. Consistent with the studies already addressed, large differences in fiber orientation became evident. Based on polar plots of the direction-dependent Young's modulus [25] computed with the Mori-Tanaka model [26,27], the flow-fiber coupling effects on the local elastic anisotropy are shown. In this context, the results of Wang and Smith [22,23] demonstrate significant differences due to flow-fiber coupling in the engineering constants describing the effective elastic anisotropy. Wang and Smith [28] present a finite element procedure regarding flow-fiber coupling with the result that the differences caused by flow-fiber coupling can be sig-nificant. Differences up to 25-50% in the velocity field and differences up to 8% in the tensile modulus were observed, depending on the chosen parameters and geometries. Wang [29] improved the algorithm of Wang and Smith [28] by optimizing the scalar viscosity in order to consider higher fiber aspect ratios with less numerical difficulties.

Originality
This study investigates the differences between the decoupled and the coupled approach of fiber suspensions with focus on the stress field of the solid composite based on the final elastic properties depending on the local fiber orientation state. In addition, the influence of flow-fiber coupling on one side and the Newtonian and non-Newtonian modeling of the matrix fluid on the other side is investigated. A detailed analysis of the flow conditions is omitted, since the coupling effects were already studied in this context in the literature cited in Sect. 1.1. The content of this study can be summarized as follows: • Application of a unified mean-field homogenization procedure for both viscous fiber suspensions and solid fiber reinforced composites. • Detailed analysis of the two-phase mold-filling simulation effects on the local fiber orientation with respect to the alternative steady state approach. • Quantification of the flow-fiber coupling effects on the stress field of the solid composite representing the local anisotropic elastic properties considering both Newtonian and non-Newtonian matrix behavior during moldfilling. • Detailed documentation of the numerical procedure focusing on stabilization methods of mold-filling simulation. • In order to establish a uniform open-source simulation tool, the use of OpenFOAM ® also for solid composite simulations besides Abaqus ® is addressed.

Notation
Throughout this manuscript, a symbolic tensor notation is used. Scalar quantities are denoted by a, b, . . . and a,

Governing equations
For both viscous fiber suspensions and solid fiber reinforced composites, the governing equations, material models and assumptions are addressed in the following. It should be noted that all governing equations must be interpreted in the context of the finite volume method. Regarding the numerical procedure, the reader is referred to Sect. 3.

Volume of fluid approach for viscous fiber suspensions
Throughout this study, the fiber orientation is described based on the second-order and fourth-order fiber orientation tensor N and N of the first kind [10,30] N( In Eq. (1) the probability density function f depends on the direction n defined on the unit sphere S, on the time t and the spatial position x. The surface element dS on S is given by dS = sin(θ )dθ dϕ/4π , with θ ∈ [0, π] and ϕ ∈ [0, 2π). During mold-filling the fiber orientation evolution is modeled based on the Folgar-Tucker equation [9,10] with the effective strain rate tensorD, the effective spin tensorW , the fiber interaction coefficient C I and the second-order identity tensor I. The scalar shear rate is defined byγ = D ·D/2 and the shape parameter refers to ξ = (α 2 − 1)/(α 2 + 1), with α denoting the fiber aspect ratio. The scalar shear rate is introduced consistently to previous studies [24,31]; other definitions may be applied along with changing C I accordingly to ensure the same fiberfiber interaction strength. Since Eq. (2) depends on N, an adequate closure method is required which can be selected from numerous approaches given in the literature, e.g., [32][33][34]. In this study the invariant-based optimal fitting closure (IBOF) [34], whose implementation details are addressed in Appendix A, is adopted. It should be noted that neutrally buoyant, isotropically interacting rigid short fibers are considered embedded in an incompressible matrix fluid without body forces [9,35,36]. Regarding further discussion, parameterizations, and admissible parameter ranges of second-and fourth-order orientation tensors of the first kind, the reader is referred to the recent work of Bauer and Böhlke [37].
In order to realize the flow-fiber coupling, the effective viscosity tensorV of the fiber suspension is formulated in terms of fiber orientation tensors N and N. Throughout this study, the effective Cauchy stress tensorσ refers to a generalized incompressible Newtonian fluid. Furthermore, the Mori-Tanaka model [26,[38][39][40] is used readinḡ The effective pressure is denoted byp, V M stands for the viscosity tensor of the matrix fluid, c F and c M = 1 − c F refer to the fiber (F) and matrix (M) volume fraction. The polarization tensor P 0 contains information about the fiber geometry and the matrix material and can be computed analytically for spheroidal fibers. In this work the expressions given in Ponte Castañeda and Willis [41] are adapted for viscosity and simplified for incompressible matrix behavior following Karl and Böhlke [40]. The inversion P −1 0 is defined on the symmetric deviatoric subset (symDev) for the considered case of incompressible suspensions with P −1 0 P 0 = P 2 representing the identity on symDev for second-order tensors. The operator · OA denotes the orientation averaging procedure (OA) [10] leading to the following expression simplified for the special case of incompressibilitȳ with the coefficients b 1 , . . . , b 5 depending on the orientation averaged tensor P −1 0 [10]. As discussed by Karl et al. [24], the formulation of the fiber-induced viscous stress based on Eqs. (3) and (4) in general is not traceless for anisotropic fiber orientation, even for tracelessD in the incompressible case. In the framework of the decoupled approach (c F = 0), onlȳ V = V M is used and, therefore, the flow field is not influenced by the evolution of the fiber orientation distribution.
The mold-filling process is based on the volume of fluid method (VOF) covering the two-phase flow of fiber suspension and air. Details about this established method are not given here, but reference is made to the literature [42,43]. Both phases are considered immiscible, isothermal, incompressible and without phase transition and body forces. The location of both phases is described by the scalar phase parameter ψ with the value ψ = 1 indicating the fiber-laden melt and ψ = 0 indicating the air phase. The following system of equations (5) is solved for ψ, the effective velocity fieldv andp with the additional force vector f being present only in the fiber-laden melt or suspension (susp) for the flow-fiber coupled approach Details about the implementation are given in Sect. 3. The surface tension of the fiber suspension is denoted by σ and the curvature κ is defined as follows [44] κ = −div grad(ψ) grad(ψ) .
Furthermore, the effective mass densityρ and the viscosity μ are defined by the following adapted mixture rules representing phase-dependent fields It should be noted that the fibers are considered withinρ but not within μ and that ρ susp = c M ρ M + c F ρ F corresponds to Voigt's average leading to the exact suspension's density. In Eq. (3) the matrix viscosity V M is considered both Newtonian and non-Newtonian in the context of a generalized Newtonian fluid. Within the present study, the shear-thinning Cross-model [45] for the polymer matrix is implemented as follows to be compared against the Newtonian behavior [40,46] The Newtonian viscosity is denoted by μ 0 and m, n represent model parameters. It should be noted that the air phase is always modeled as a Newtonian fluid. In Eq. (9) the assumption of rigid fibers D F = 0 is used in order to expressγ M with respect toD and c F usingD = (1 − c F ) D M . Since the flow-fiber coupling affects the constitutive behavior, γ M is used instead ofγ introduced in the context of the Folgar-Tucker equation (2). It should be noted that the fiber orientation evolution is considered only dependent on the macroscopic kinematic fieldsD,W andγ . As a consequence, the coupling effects on the fiber orientation evolution are only present based on changes in the flow field compared to the decoupled approach. This procedure ensures a comparable fiber orientation evolution in order to study the coupling effects purely due to viscous anisotropy.

Solid fiber reinforced composites
After fluid-solid transition, effective linear elastic behavior σ =C[ε] is assumed for the solid fiber reinforced composite, withε denoting the effective strain andC representing the effective stiffness tensor. It is assumed that the geometry and fiber orientation do not change during the fluid-solid transition. For clarity, the basics regarding mean-field modeling are omitted here and the reader is referred to further literature, e.g., Nemat-Nasser and Hori [47] or Kanaun and Levin [48]. Based on the Mori-Tanaka model [26], the effective stiffness tensor readsC = C M +R with the abbreviationR defined as follows [27,40,49] It should be noted that only homogeneous and isotropic phases are considered. In Eq. (10), the abbreviation δC refers to C F − C M with the stiffness of the fiber material C F and of the matrix material C M . The polarization tensor P 0 can be computed analytically based on the expressions given in Ponte Castañeda and Willis [41]. Analogously, · OA refers to the orientation averaging procedure given in Advani and Tucker [10]. The stationary balance of linear momentum without body force densities 0 = div(σ ) is used as the governing equation for the effective displacement fieldū The bulk modulus is represented by K and the shear modulus by G. Based on the right minor symmetry ofR, the mappingR[grad(ū)] inherently uses only the symmetric part of grad(ū) referring toε. The tensorR contains the information about the local fiber orientation state based on orientation averaging. After solving Eq. (11) forū, the effective strain is computed directly by the kinematic relationε = sym grad(ū) and the effective stress follows from σ =C[ε].

Scale separation and requirements for homogenization
Fiber orientation tensors inherently represent the averaged orientation state, e.g., on every numerical grid cell in the context of coarse graining. Therefore, labeling N and N as effective fields by( ·) is omitted throughout this manuscript. All other fields are labeled as effective fields to emphasize that they are averages regarding the scale on which the orientation tensor description is valid. In this context, every grid cell used for numerical computation is seen as a statistically representative volume across which the described mean-field homogenization is performed. The spatial distribution of the fibers is assumed to be statistically homogeneous and the ergodic hypothesis is applied. In addition, the considered microstructure does not contain cracks or voids, the Hill-Mandel condition [50,51] is assumed to be valid and phase-wise constant properties are used. Further information can be found in the literature, e.g., in Nemat-Nasser and Hori [47] or in Torquato [52]. Regarding the uniform meanfield homogenization of viscous fiber suspensions and solid short-fiber reinforced composites, reference is made to the recent work of Karl and Böhlke [40].

Mold-filling simulation using OpenFOAM ®
For both fluid and solid simulations the 2D rib geometry shown in Fig. 1 is used. This generic geometry is inspired by the bifurcation flow of Haagh and van de Vosse [53] and represents a classic thin-walled component. The length h is 2.5 · 10 −3 m. The geometry is chosen such that the two ribs experience different flow states towards the end of the mold-filling process, which is discussed in detail in Sect. 4.1. It is shown that the two-phase simulation considering the actual mold-filling is necessary in order to predict the fiber orientation state correctly.
The implemented VOF solver is based on the inter-IsoFoam solver of OpenFOAM ® (version 2106) for the isothermal mold-filling of two immiscible and incompressible fluids with Newtonian or non-Newtonian viscosity [54]. The solver in use differs from the interFoam solver by using the isoAdvector scheme improving the advection towards a sharp interface [55][56][57][58]. Within the fluid solver, volume forces due to gravity are neglected in the balance of linear momentum analogous to the solid solver. The implemented methods of stabilizing the mold-filling simulations are addressed in Sect. 3.4. Regarding the temporal discretization, the second-order implicit backward scheme is used with the time step t determined by the maximum Courant number max(Co) = 0.02. This value corresponds to the order of magnitude of the value of 0.05 proposed by Larsen et al. [58] for the related interFoam solver. To ensure sufficient accuracy, t is limited by the maximum value of 10 −3 s. All simulations are stopped at the filling level of 99.999%. In addition, only second-order spatial discretization schemes are used. For solving the systems of equations the multigrid solver GAMG is chosen with GaussSeidel as the smoothing method. In addition,p andv are computed based on the PISO algorithm with 10 correction steps and no momentum predictor step. As a result of a grid study 80 orthogonal and uniformly distributed cells per h are set to resolve the fields with sufficient accuracy.
Throughout this study, polypropylene (PP) is considered as the thermoplastic polymer matrix with the melt density ρ M = 728 kg/m 3 (≈ 0.8 solid density [59] of 910 kg/m 3 [60]). The Newtonian kinematic viscosity is set to ν M = 0.779 m 2 /s which corresponds to ν 0 in the shearthinning Cross model given in Eq. (9) with the parameter values n = 0.749 and m = 17.366 · 10 −3 [61]. Please note that the kinematic viscosity ν and the density ρ are used as input in OpenFOAM ® instead of the dynamic viscosity μ = νρ. The surface tension is set to σ = 0.028 N/m [62]. For the air phase the default parameters ν air = 1.48 · 10 −5 m 2 /s and ρ air = 1 kg/m 3 are used. In this study glass fibers are considered with the density ρ F = 2540 kg/m 3 [60]. By defining h,v b , ν M , ρ M , σ as the characteristic quantities, the Reynolds number (Re) and the capillary number (Ca) describing the problem have the following values Based on the value of Re, the mold-filling process approximately corresponds to creeping flow conditions and Ca shows that viscous forces are larger than surface tension forces. The aspect ratio of the short glass fibers equals α = 20 with C I = 0.01 as the fiber-fiber interaction parameter affecting the orientation tensor evolution given in Eq. (2). The fibers are considered with a volume fraction c F = 0.2 which corresponds to 41.1 wt − % based on the solid densities. Analogous to Karl et al. [24], the thermodynamic consistency is checked during mold-filling based on the isothermal dissipation inequalityD ·V[D] ≥ 0. The initial and boundary conditions for every field are summarized in Table 1 with the boundaries wall, inlet and outlet sketched in Fig. 1. The zeroGradient condition refers to a vanishing gradient normal to the respective boundary. The zeroGradient condition for ψ corresponds to a contact angle of 90 • at the wall. The boundary condition fixedFlux-Pressure adjusts the pressure gradient in order to get a flux on the boundary which is consistent with the velocity boundary condition [54]. At the inlet a parabolic velocity profile is set with the bulk velocityv b = 0.1 m/s [53]. By choosing the standard interIsoFoam solver, pinning of the contact line is observed. However, experiments in comparable settings show that form-filling flows are achieved with the selected process data described above [63,64]. The observed contact line pinning is most likely related to the large capillary number given in Eq. (12). Under large Ca number conditions it is known that a large deformation of the free surface area can occur before the occurrence of a moving contact line. As a result, the advancing contact angle, which is related to Ca through the Cox-Voinov law [65,66], basically decouples from the static contact angle [67]. This physical process cannot be captured by the employed VOF formulation in which only the static contact angle is prescribed. In order to enable the corresponding simulations at high Ca with a moving contact line in the standard interIsoFoam solver a novel boundary condition, referred to as localAirVent, is implemented. The localAirVent boundary condition-which is inspired by a freely escaping air phase in previous mold flow simulations [68][69][70]-enables an outflow of air in an area around the contact point. The boundary condition airVentBC allows the air phase to leave the domain across the outlet until the boundary cells are almost filled with polymer. Then, the condition switches linearly from zeroGradient to the standard no-slip condition and prevents the suspension from flowing out. Both localAirVent and airVentBC are based on the following equation with the implemented blending function g(ψ) and the outward pointing unit normal vector n at each boundary The mixedPressureOut boundary condition used by, e.g., Ospald [70] refers to a constant reference pressure if the boundary cells are completely filled with air. If the boundary cells are completely filled with suspension, then the zero- Gradient condition is applied. In between, a mixed boundary condition is applied until the cells are completely filled with suspension following the condition with the reference pressure chosen as p ref = 0. All other boundary conditions listed in Table 1 are self-explaining and OpenFOAM ® standard. The 2D flow is achieved by assigning the empty boundary condition to the front and back faces.

Solid composite simulation using Abaqus ®
Regarding the solid composite simulations, the Young's modulus of 1.6 GPa is used for PP and 73 GPa for the glass fibers [60]. The Poisson's ratio of PP equals 0.4 and 0.22 for the glass fibers [60]. With these parameter values, the equations (10) and (11) are uniquely determined. The inlet is clamped with suppressed displacement and rotation. The displacementū 1 = 1.5 · 10 −4 m is applied at the two outlets of the ribs, which results in a bending of the whole structure. All details about the boundaries are illustrated in Fig. 1.
It should be noted that the fluid-solid transition during and after the mold-filling process is assumed to occur instantly without modification of the fiber orientation. The solid composite simulations throughout this study are mainly carried out with the commercial finite element software Abaqus ® CAE Standard/Explicit (version 2020). The anisotropic stiffness corresponds to the algorithmic tangent within a user-material (UMAT) subroutine implemented in Fortran 90. The fiber orientation data corresponds to a predefined field as an additional input of the UMAT. Based on the 2D load case, the 4-node bi-linear plane strain element CPE4 is used for meshing the geometry shown in Fig. 1. As a result of a grid study 80 orthogonal and uniformly distributed cells per h are used to resolve the stress field sufficiently. Moreover, the mapping of the orientation data from the fluid to the solid mesh is done based on the nearest-neighbor procedure used by Gajek et al. [71] extended by the trace-preserving correction scheme of Kuzmin [72]. It should be noted that the orientation data is mapped from the cell centers of the openFOAM ® mesh to the nodes of the Abaqus ® mesh. Since the UMAT is used at the element integration points, the nodal orientation values are interpolated by Abaqus ® . It is observed that elements with quadratic trial functions lead to negative eigenvalues of the system matrix. Therefore, using elements with linear trial functions is suggested.

Solid composite simulation using OpenFOAM ®
As a basis of the present study the solids4foam toolbox [73][74][75] for OpenFOAM ® is used. For further information regarding the finite volume method in solid mechanics the reader is referred to Cardiff and Demirdzić [76]  The material parameters are given in the section above. At the inlet, the fixedDisplacement boundary condition is used realizing the clamping. For the two ribs, the solidDirection-Mixed boundary condition is used withū 1 = 1.5 · 10 −4 m as the horizontal displacement and free vertical displacementū 2 . For all other boundaries the solidTraction condition is set as zeroGradient. The plane strain state is achieved by assigning the empty boundary condition to the front and back faces of the geometry. For solving the system of equations the multigrid solver GAMG is chosen with diagonal-based incomplete Cholesky preconditioner as the smoothing method. All chosen spatial discretization schemes are of second order and 80 orthogonal and uniformly distributed cells per h are used based on a grid study. Please note that in this study the same mesh is found to be sufficiently accurate for the mold-filling simulations, the finite element and for the finite volume simulations. In order to correct possible interpolation errors, which may occur during the transition between the fluid and the solid mesh, the tracepreserving correction scheme of Kuzmin [72] is implemented analogously to the finite element procedure described above. In contrast to the finite element implementation, methods for stabilization must be actively implemented in this context, which is described in Sect. 3.5.

Fluid solver stabilization methods
In the following subsections methods of stabilizing the twophase fluid solver for computing the flow-fiber coupled mold-filling process are addressed in detail. The discussion of the methods is adapted to the use of OpenFOAM ® and is intended to support engineering practice implementing flowfiber coupled computations. All of the shown methods are implemented to enable parallel computations.

Implicit-explicit splitting of the coupling term
In order to improve the stability, as many terms of f given in Eq. (6) as possible should be discretized implicitly regarding the primary fieldv. It turns out to be numerically unfavorable that gradients can only be discretized explicitly in OpenFOAM ® and that the coupling term is dominated by velocity gradients [70]. The coupling term f given in Eq. (6) can be rewritten as follows Only one part of the second term on the right-hand side of Eq. (15) can be discretized implicitly by using the product rule and the incompressibility condition leading to div By using Eq. (16) and Eq. (15), the implemented splitting of f reads The implicit-explicit splitting of all other terms in the system of governing equations (5) for computing ψ,p andv is not changed and corresponds to the implementation of the basic solver interIsoFoam described in Sect. 3.1.

Treatment of the viscous stress
As already discussed by Karl et al. [24] the fiber-related part of the viscous stress is not traceless in the case of anisotropic fiber orientation, even for tracelessD. As a consequence, additional contributions to the spherical stress − p I are present, which can be spatially strongly inhomogeneous with negative effects on the pressure solver convergence. Since the pressure acts as a Lagrangian multiplier for satisfying the incompressibility constraint, one can use, without loss of generality, the following expression for the fiber-related viscous stress Since the mass balance div(v) = 0 is fulfilled in terms of fluxes over the cell faces, the cell-centeredD is not necessarily traceless due to interpolation errors. Simulation results show that this error might be large close to the interface with the error decreasing to tolerable values as the distance from the interface increases. In order to ensure, that only traceless D is used within the viscous stress tensor, the following final expression is implemented [39]

Spatial restriction of the coupling term
In the following, the spatial restriction of f given in Eq. (6) is described in detail. As discussed in the literature, e.g., [53,[77][78][79], the three-phase contact point (TPCP) shown in Fig. 2 is critical in view of stress singularities. Furthermore, the large material contrast of the air and melt phases combined with the coupling in the area of the phase front has an unfavorable effect. In order to obtain convergence, the coupling is switched on successively with increasing distance from the phase front. Therefore, anisotropic viscosity effects directly at the phase front cannot be taken into account with this procedure. But this is not critical since the analyzed orientation tensor further downstream is not significantly influenced by this gradual introduction of the flow-fiber coupling behind the phase front.
In order to circumvent interface-induced coupling instabilities, the following procedure is implemented with all areas shown in Fig. 2. The area V susp occupied by the fiber suspension is divided into three sub-areas. In the area V decoupled close to the interface, no flow-fiber coupling is present avoiding large explicit coupling terms. Next, in the transition zone V trans the coupling strength is increased linearly with increasing distance from the interface. The full flow-fiber coupling is present in V coupled . As an input, d min is set as the minimum distance between the phase front (interface) and the beginning of flow-fiber coupling. In addition, d trans is set as an input representing the width of the transition area. In every time step, the positions of the interface grid cells x interface are detected by grad(ψ) > 5000 m −1 , which turns out to It should be noted that Eq. (21) refers to both the implicit and the explicit part of the coupling term. Within this study the distances are set to d min = d trans = h/10.

Under-relaxation of the explicit coupling term
Under-relaxation is a common method in computational fluid dynamics improving numerical convergence [80]. Throughout this study the explicit part of the coupling term is under-relaxed (ur) as follows with the relaxation factor 0 < r < 1, the coupling term of the previous time step f n expl and f n+1 expl referring to the explicit coupling term depending on the updated fiber orientation tensor and flow kinematics. As a result of under-relaxation, the explicit coupling term is not used entirely updated in the linear-momentum balance. As described in the literature [80], the relaxation parameter r depends on the flow problem and, therefore, cannot be chosen a priori. Within this study it is observed that strong under-relaxation is necessary for strong flow-fiber coupling (c F > 0.1, α > 10), even for kinematically simple channel flows. Throughout this study, r = 0.05 is fixed for the simulation setup. There are no negative effects on the convergence due to this small value of r since both the fiber orientation and the flow field quickly reach stationary values. The simulations results show that the difference f n+1 expl − f n expl is only present around the transition region (see Sect. 3.4.3), in which the coupling is not considered fully present anyway.

Stabilizing the Folgar-Tucker equation
As discussed by Paschkewitz et al. [81], an explicit solution procedure of Eq. (2) may lead to instabilities, which are manifested by a loss of positive semi-definiteness and by invalid values of the orientation components. Therefore, as many terms as possible are discretized implicitly regarding N, as shown in Appendix A in order to make the manuscript self-contained. Throughout this study, no numerical issues in view of solving the Folgar-Tucker equation occur by applying the procedure addressed in Appendix A. Nevertheless, further possibilities for stabilizing and correcting the Folgar-Tucker equation are summarized below, which are useful for the reader in case of numerical issues.
Since the Folgar-Tucker equation represents a purely convective equation, an artificial diffusion term R art = β N can be used to stabilize the solution procedure with β representing the artificial diffusion. Mezi et al. [5] used this stabilization method within the Fokker-Planck equation for the probability density function f . Alternatively, artificial diffusion can be restricted to regions where N loses its positive semi-definiteness [81,82] The parameter β can be modeled depending on the local grid spacing and a proper value depending on the flow problem [81,82]. It should be noted that unsuitable values of β lead to numerical issues and, therefore, a parameter study is suggested [82]. The trace condition tr(N) = 1 is not violated by adding artificial diffusion terms to the Folgar-Tucker equation. This can be easily shown by expanding the evolution equation for s:=tr(N) given in Linn [83] by the artificial diffusion term as follows Please note that the derivation of Eq. (24) requires N[I] = N to be fulfilled by the closure in use. Since always physical initial conditions for the orientation state are used with s = 1, the right-hand side of Eq. (24) vanishes. As a consequence, s = 1 is preserved. In addition to stabilizing the Folgar-Tucker equation, nonphysical orientation tensor results can be corrected. The correction schemes refer to eigenvalues or diagonal components of N violating the physical interval [0, 1]. Verweyst and Tucker [84] developed a procedure to correct the eigenvalues of N after solving the Folgar-Tucker equation and before the orientation state is used for flow-fiber coupling. For a tracepreserving correction based on setting non-physical negative eigenvalues to zero, the reader is referred to Kuzmin [72] and the references therein for further correction schemes.

Remarks on the mesh and time discretization
Within this study, it is observed that homogeneous and orthogonal mesh cells are the most suitable for the two-phase mold-filling simulations in view of computation time. This is due to the fact that non-orthogonal cells require a smaller time step to obtain the same Courant number for similar mesh fineness. Furthermore, with respect to the coupling term, coarse grids lead to stable solution behavior for spatially strongly inhomogeneous fields. This is due to the fact that gradients are smoothed in this way and numerically unfavorable local peaks in the coupling term are thus prevented. Despite the stability, coarse meshes are not recommended for improving stability due to spatial inaccuracy.
As already described in Sect. 3.1 the second-order implicit backward scheme is used for time discretization. Although this is an unbounded scheme [54,85], no instabilities in the context of time discretization are observed in this study. In case of stability issues, the bounded implicit Euler scheme [54] may be applied to stabilize the time integration. For progressive waves the numerical damping of the Euler scheme is addressed by Larsen et al. [58]. It should be noted that this first-order accurate method may only be applied if sufficient temporal resolution is ensured or if the integration leads to a steady state solution.

Solid solver stabilization methods
In order to make the manuscript self-contained, the standard implementation of Eq. (11) in solids4foam is explained below, which is extended by the Mori-Tanaka model for the simulation of anisotropic fiber-reinforced composites. In its most compact form, the linear momentum balance reads In order to improve convergence, Laplacian terms are added to Eq. (25) as follows with both implicit and explicit treatment [74,75] As described by Shayegh [75], the first two terms do not influence the converged solution. The prefactor k = 4G F /3 + K F is used [74,86,87] with the elastic properties of the fiber phase. During the study it is observed that a smaller chosen k, for example calculated with the elastic constants of the matrix, may have a negative effect on the convergence. In order to damp non-physical oscillations, the Rhie-Chow stabilization method [74,76,[87][88][89] already implemented in solids4foam is applied. As a consequence, the linear momentum balance reads with the explicitly treated Rhie-Chow stabilization term [75,76] with a representing a scaling factor. A parameter study shows that the solids4foam default value a = 0.1 is sufficient to damp the occurring numerical oscillations to an acceptable magnitude. It should be noted that both ū and div grad(ū) represent Laplacians, but div grad(ū) is used to express a different discretization compared to ū. In addition, the Laplacian term G M ū being part of div(C[grad(ū)]) as given in Eq. (11) can be treated implicitly to further stabilize Eq. (27). Based on a comparison of a purely explicit treatment, no advantage is observed in terms of stability and computational time. Therefore, the explicit discretization of div(C[grad(ū)]) shown in Eq. (27) is implemented.

Coupling effects on the fiber orientation
In Fig. 3 the fiber orientation field N is shown in the area of interest defined in Fig. 1. For simplicity, only the results regarding the non-Newtonian case are shown with the decoupled approach in the left column and the coupled approach in the right column, respectively. The plots in Fig. 3 serve to clarify the orientation state as a basis for the solid composite simulations of Sect. 4.2, which are based on the fiber-induced anisotropic stiffness. This is important because the coupling effects will be analyzed in form of difference fields. By comparing both columns in Fig. 3 against each other one can see that the coupling effects (see Sect. 1.1) are reproduced, e.g., the distinct isotropic orientation in the horizontal channel center region caused by the flattening of the velocity profile. In the considered 2D flow, mainly N 11 and N 22 are influenced by the kinematics, whereas N 33 evolves due to trace preservation. The component N 12 provides information about the eigensystem of N with N 12 = 0 referring to an eigensystem coinciding with the spatial axes. Regarding the 2D flow consideration, N 13 = N 23 = 0 holds. The area of interest shown in Fig. 3 can be divided into two main parts, the lower horizontal area and the vertical ribs. Both parts can be simplified as a simple channel flow, in which the fibers align in the main shear direction. For the horizontal part, this direction is e 1 represented by large values of N 11 . For the vertical ribs, this direction is e 2 leading to large values of N 22 .
By comparing the final fiber orientation states below and inside both ribs at the filling level of 99.999%, the influence of the different filling behavior is clearly visible: The left rib is completely filled first, while the right rib and the horizontal channel are not completely filled at this particular time step. Once the left rib is filled, no fiber orientation evolution is present inside and the kinematics in the area below change into a shear-dominated one leading to an alignment in e 1direction. This horizontal shear is not present in the area below the right rib since the horizontal channel is filled before the right rib. Therefore, N 22 is the favored fiber orientation below the right rib caused by vertical shear and a deflection of the flow. In the following, the filling behavior of the two ribs is investigated further in a comparative manner.
In the top-right plot of Fig. 3 two points are marked for which the temporal fiber orientation evolution is investigated. For the point in the left rib, the diagonal components of N are plotted over time in Fig. 4a and for the right rib in Fig. 4b, respectively. The analysis is restricted to the coupled non-Newtonian approach. Starting from the isotropic initial state marked with 1 the fiber orientation begins to evolve once the phase front (represented by ψ) passed the two considered points. In time interval 2 shortly behind the phase front, a transient orientation state is shown for both points. Subsequently, the orientation evolves towards an asymptotic state at 3 , which describes the mold-filling process of both ribs and the horizontal channel without any part of the cavity being completely filled. It should be noted that the asymptotic state only depends on the given kinematics, the aspect ratio of the fibers and the fiber-fiber interaction strength [31]. Since the kinematics are different at the two points during time interval 3 , also the orientation state differs. Note that the orientation state observed during time interval 3 corresponds to the state that would occur if the flow through the cavity at the two points were stationary, or in other words, if no two-phase simulation of the actual mold-filling process is carried out. At approximately t = 0.58 s the left rib is completely filled with suspension and, as already described above, the kinematics change into a shear-dominated one in the area around the considered point. This also affects the orientation evolution in the right rib at this time step since the flow is accelerating in the horizontal channel. Already a short time later at approximately t = 0.62 s the horizon-tal channel is also completely filled leading to a complete deflection of the flow into the right rib. Until the cavity is completely filled at t = 0.69 s, this change in kinematics leads to a new orientation state at 4 defining the anisotropy of the manufactured part after fluid-solid transition. Since level 4 is clearly different from 3 , the necessity of simulating the actual mold-filling becomes apparent. A simulation based purely on the steady-state flow through the geometry cannot represent this change in kinematics due to the filling process and predicts the fiber orientation differently. Therefore, the filling of individual areas of the cavity must necessarily be taken into account in order to correctly predict the evolution of the fiber orientation and the associated anisotropy.
In order to further quantify the mentioned differences between 3 and 4 in view of fiber orientation, the local anisotropic viscosity is investigated in the coupled approach. For this purpose, the following two measures e V and e V,aniso are defined, which quantify the difference between 3 (left rib fills quasi-stationary at t = 0.5 s) and 4 (shear-dominated flow below left rib at t = 0.69 s) The total viscosity difference is represented by e V , whereas e V,aniso addresses the change of anisotropy, both as a scalar measure for simplicity. The anisotropic part ofV is defined as follows with P 1 representing the identity on spherical second-order tensors. In Fig. 5 the viscosity measures defined in Eq. (28) are shown for the left rib. The results indicate that the differ- Non-dimensional stress difference σ * and ranges of values with respect to the different mold-filling approaches ence in orientation between 3 and 4 significantly affects the local effective viscosity and also significantly changes the degree of anisotropy. In the following, the analysis is restricted to the comparisons shown in Fig. 6 with the fiber orientation differences defined by Karl et al. [24]. In Fig. 6, vertical differences address the effects of flow-fiber coupling for both viscosity models, whereas horizontal differences account for the effect of the viscosity models for both decoupled and flowfiber coupled approach. For all four cases, the range of values is given with respect to the area of interest.
Based on Fig. 15 in Appendix B.1, where additional fiber orientation results are provided, the local orientation differences with respect to different viscosity models in the coupled approach range from − 40 to 41%, regarding the possible interval [0, 1] of N 11 , N 22 and N 33 . In addition, Fig. 14 in Appendix B.1 shows the effect of the viscosity model for the decoupled approach with local orientation differences − 27 to 30%. Therefore, one can state that for both the decoupled and the coupled approach, the viscosity model has a great influence on the local fiber orientation state. Furthermore, the relevance of the viscosity model changes with the flowfiber coupling, as the range of values increases compared to the decoupled case. Please note that the latter statement only holds for flow cases in which the shear rate is sufficiently large to cause significant shear-thinning in the matrix fluid. This is given under the selected process conditions and is addressed in Appendix B.1 in the context of Fig. 16.
For the sake of clarity, only the difference fields regarding the pure effects of flow-fiber coupling for each viscosity model are provided in the present section. In Fig. 7 the local orientation difference between the decoupled and the coupled approach is shown for the Newtonian case. The same difference is given in Fig. 8 for the non-Newtonian matrix fluid. Areas highlighted in red indicate an alignment which is larger in the coupled case and blue areas indicate a stronger alignment in the decoupled case. By comparing Figs. 7 and 8 one can see that the difference fields show similar characteristics for both considered matrix viscosity models. Regarding the possible interval [0, 1] of N 11 , N 22 and N 33 , the local effect of flow-fiber coupling ranges from − 36 to 31% for the Newtonian and from − 37 to 38% for the non-Newtonian case and, therefore, can be evaluated as significant. It can be seen that the flow-fiber coupling has approximately the same effect for both viscosity models. In the non-Newtonian case, the differences are slightly more distinct. However, this does not mean that the fiber orientation in the coupled case is almost identical for both viscosity models. In the coupled case, for example, the choice of viscosity model has a larger influence on the local fiber orientation than the flow-fiber coupling. This can be seen in Fig. 6 in view of the given ranges of values. Therefore, not only the coupling is crucial for the local evolution of the fiber orientation, but also the consideration of the shear-thinning behavior.

Coupling effects on the stress field
The fiber orientation states investigated above are now used within finite element simulations for solid composites. Independent of the fluid model during mold-filling, the solid composite is modeled linear elastic with its anisotropy completely defined by N and N. Thus, this study docks directly to the results of Karl et al. [24], which analyze the coupling effects on the local elastic anisotropy, but not on the stress field in the solid composite under load. Analogous to the previous section, Fig. 9 addresses the four investigated cases and the corresponding comparisons. In the following, the non-dimensional stress difference field σ * is considered with the respective definitions given in Fig. 9. Note that these definitions must be applied for each stress component i j = 11, 12, . . . individually. Again, horizontal differences address the influence of the viscosity model during mold filling for both decoupled and coupled approach, whereas vertical differences consider the effect of coupling for both viscosity models.
As a supplement, Fig. 17 in Appendix B.2 shows the stress field both for the decoupled and the coupled approach regarding the non-Newtonian case during mold filling. The corresponding stress difference field is given in Fig. 11 with the range of values from − 10.59 to 14.92%. Analogous to the previous section, only the non-dimensional stress difference fields regarding the pure effects of flow-fiber coupling for each viscosity model are considered in the section at hand. For the Newtonian viscosity Fig. 10 addresses these coupling effects with the range − 12 to 12.9%. All components not shown are zero due to the plane strain state and N 13 = N 23 = 0. Both ranges indicate, that the coupling effect on the fiber orientation evolution during mold-filling has a significant influence on the local stress of the solid composite under loading after fluid-solid transition. Please note that normalization is done with each maximum stress component, which emphasizes the importance of the coupling in the context of the specified ranges. If stress peaks occur locally due to the geometry, even small differences in fiber orientation caused by flow-fiber coupling during mold-filling can lead to large differences in the stress field. This is the case in the area of the left rib, in combination with large differences in fiber orientation. There, kinematically complex flow conditions are additionally present, which lead to large differences in the fiber orientation, since viscous anisotropy effects are fully present. Even a small load without stress peaks can lead to significant differences in local stress in such areas, as is the case, e.g., in the area under the right rib. In summary, also against the background of the stress distribution, both the necessity of coupling and the consideration of shear-thinning matrix behavior are confirmed.
The difference fields regarding the influence of the viscosity model can be found in Appendix B.2. For the decoupled case the range of the local stress difference is − 6.04 to 4.61% as shown in Fig. 18 and for the coupled case − 4.74 to 6.51% as shown in Fig. 19. As a consequence, one can state that different viscosity models during mold-filling affect the local stress significantly both for the decoupled and the coupled approach. Compared to the previous section, the relevance of the viscosity model does not change remarkably with the flow-fiber coupling, as the range of values is comparable to the decoupled case. However, the difference between the stress fields based on the different viscosity models is not negligible. More details can be found in Appendix B.2.

Abaqus ® vs. OpenFOAM ® for solid composite simulation
In this section the stress predictions based on the extended finite volume solver OpenFOAM ® (FVM in the following) are presented. Furthermore, a comparison is made with the results based on finite element simulations with Abaqus ® (FEM in the following). The analysis is limited to the coupled approach for non-Newtonian matrix viscosity since both flow-fiber coupling and shear-thinning behavior significantly affect the orientation and the stress field. In Fig. 12    large and strongly inhomogeneous gradients are present. To investigate the differences between FVM and FEM in more detail, Fig. 13 shows the components of the effective stress tensor along two paths. The two paths are shown in the upper left plot in Fig. 12 representing a cut through the horizontal channel at x 2 = h/2 and through the left rib at x 1 = 5h, respectively. In the considered area of interest, the results of FVM and FEM show an overall good agreement.
For the problem given in this study, the results show that compared to the FEM, the extended FVM solids4foam toolbox is suitable for calculating the stress field in the solid composite. This makes it possible to carry out the simulation chain starting with the mold-filling up to the loading of the solid composite in the open-source code OpenFOAM ® . Further advantages are the uniform environment, which enables standardized workflows, and the high flexibility, i.e. the possibility to expand the contained basic solvers. A disadvantage is that, depending on the problem and the code extension, suitable stabilization methods must be found and any parameters must generally be selected depending on the given problem. This goes hand in hand with a larger implementation effort compared to the implementation of an Abaqus ® user-material subroutine.
As any further comparison between FVM and FEM is not the focus of the study, the reader is referred to the review article of Cardiff and Demirdžić [76] and the references therein. The work of Demirdzić [90], taken as an example, compares FVM with FEM and evaluates FVM as outperforming FEM and as an equal alternative, depending on the problem. However, the disadvantage is listed that solid mechanics based on FVM is not widely used.

Summary and conclusion
In this study the effects of flow-fiber coupling during moldfilling are investigated in view of the stress field in the manufactured solid composite after fluid-solid transition. The two-phase mold-filling simulation of the fiber suspension is carried out with the finite volume solver OpenFOAM ® based on an extended VOF solver interIsoFoam [54]. Various methods for numerical stabilization are implemented. The fiber orientation taken from mold-filling simulations is used for the solid composite defining its local anisotropic elastic behavior. The fiber orientation itself during mold-filling is simulated both with Newtonian and non-Newtonian matrix behavior and also with and without flow-fiber coupling. Both the viscous fiber suspension and the solid fiber reinforced composite are modeled micromechanically unified with the Mori-Tanaka mean-field model [26,27,[38][39][40]49]. The simulation of the solid composite is carried out with the finite element solver Abaqus ® using an implemented subroutine for the anisotropic elastic behavior. In addition, the OpenFOAM ® solids4foam toolbox [73][74][75] is extended for simulating the solid composite and the differences to the results from Abaqus ® are discussed briefly.
The present study is based on the following main assumptions. The fiber geometry is assumed to be constant and the fibers are treated as rigid bodies during mold-filling. Based on the Folgar-Tucker equation with isotropic fiber-fiber interaction [9,10], the fiber orientation tensor of the first kind and of second order [30] is assumed to be a sufficient measure of anisotropy. The closure problem is treated by using the IBOF closure [34] both for the Folgar-Tucker equation and for the orientation average procedure. The mold-filling process is modeled incompressible, isothermal and without fluidsolid transition. To account for two different polymer matrix behaviors during mold-filling, the polymer matrix is modeled as a generalized Newtonian fluid with rigid fibers, whereas the solid composite and the fibers within are modeled linear elastic. Both the fiber suspension and the solid composite are homogenized on the scale of a grid cell, which is seen as a representative volume element. Voids and cracks are excluded and the Hill-Mandel condition [50,51] is assumed to be fulfilled on the grid cells. The present study is based on a double-rib geometry as a generic injection molding geometry.
Regarding the implemented methods for numerical stabilization the implicit-explicit splitting is concluded as important. If it is possible in a solver to treat the coupling term completely implicitly in the linear momentum balance, the greatest possible stability will be obtained. But this is not possible in OpenFOAM ® , which is why further methods have to be implemented to enable suspension simulations with large fiber volume fractions and aspect ratios. The underrelaxation of the explicit part of the coupling term is observed to be a powerful method for flows in which low relaxation factors can be used without negatively influencing the solution behavior. If only relatively large relaxation factors are allowed to be used, a combination with the method of spatial restriction of the coupling term is recommended. Furthermore, the described treatment of the viscous stress can be used to improve the convergence of the pressure solver. It can be concluded that stabilizing the code strongly depends on the given problem and that a coordination of the methods is necessary. Within this study, the combination of all described stabilization methods leads to a robust numerical behavior, provided that the adaptation of the methods is problem-specific. Since this is difficult a priori, extensive preliminary studies are recommended.
The simulation results can be summarized as follows and the associated conclusions can be drawn: • For both the decoupled and the coupled case, the choice of the viscosity model has a non-negligible influence on the fiber orientation (decoupled: − 27 to 30%, coupled: − 40 to 41%) and on the stress field in the solid composite (decoupled: − 6 to 4.6%, coupled: − 4.7 to 6.5%). Since real polymer melts show a shear-thinning behavior, the necessity of this modeling can be quantifiably justified. • The differences in fiber orientation between the viscosity models are larger in the coupled case than in the decoupled case. In this context, differences in the stress field in the decoupled and coupled case are comparable. Thus, it can be concluded that the relevance of the shear-thinning behavior increases with coupling. • For both viscosity models considered, it is shown that the coupling has a significant influence on the fiber orientation (Newtonian: − 36 to 31%, non-Newtonian: − 37 to 38%) and on the stress field (Newtonian: − 12 to 12.9%, non-Newtonian: − 10.6 to 14.9%). Based on the fact that the anisotropic viscosity model describes the real fiber suspension more physically than the neglect of fiber-induced viscous anisotropy, the implementation of flow-fiber coupling is recommended. • The influence of flow-fiber coupling is quite comparable for both viscosity models. Only the comparison between the two viscosity models reveals the necessity of the shear-thinning behavior. • The two-phase simulation approach leads to distinct differences in the local fiber orientation compared to the steady-state flow through the geometry. Therefore, changes in kinematics due to the actual filling process have to be predicted by the simulation setup in order to estimate the local anisotropic properties. • A brief comparison between FVM and FEM shows that a uniform simulation environment in OpenFOAM ® is suitable for both mold-filling and for solid composite simulations. The stress field predicted by the two numerical methods shows an overall good agreement in the considered area of interest.
Finally, it should be noted that the assumptions described above are aimed at considering only the effects of flow-fiber coupling on the fiber orientation and thus on the stress field in the solid. Influences that are more complex than the viscous behavior of the matrix are therefore not taken into account, such as the effects of temperature, fiber-induced anisotropic heat conduction and fluid-solid transition. In this context, the recent publication of Dietemann and Bierwisch [91] indicates that accurate rheological modeling has a larger impact on fiber orientation (and thus on the resulting stress field) than complex modeling of orientation dynamics. In addition, the anisotropic behavior is estimated based on mean-field homogenization which has advantages in terms of computational effort. In particular the mold-filling simulation requires a computational effective estimation of the anisotropic viscosity on each grid cell. On the other hand, critical stress peaks within the highly heterogeneous microstructure are not resolved by mean-field methods. As a consequence, the conclusions listed above have to be interpreted against this background. Moreover, it is noted that injection molding processes are often modeled as flows through narrow gaps, as discussed by, e.g., Tucker [92]. In the present work, this socalled Hele-Shaw simplification (see, e.g., Spurk and Aksel [93]) is deliberately omitted because the geometry under consideration is complex. Another point is that the fiber orientation dynamics can be influenced by small wall distance in such a way that the movement of the fibers is restricted. In the literature this is referred to as confined flow, for which the modeling of fiber orientation evolution was adapted in the studies of, e.g., Perez et al. [94] and Scheuer et al. [95].
Acknowledgements The partial financial support by the Friedrich und Elisabeth Boysen-Stiftung (Grant BOY-163) is gratefully acknowledged. The authors acknowledge support by the state of Baden-Württemberg through bwHPC.
Funding Open Access funding enabled and organized by Projekt DEAL.

Data availability
The data of this study are available from the corresponding author upon reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Implicit-explicit splitting of the Folgar-Tucker equation
In this supplementary section, the implicit-explicit splitting of the Folgar-Tucker equation is addressed. The implicitexplicit splitting of the IBOF closure terms [34] is given based on the work of Heinen [96] with respect to the symmetrypreserving representation of Karl et al. [24]. It should be noted that the incompressibility condition tr(D) = 0 is used within the closure terms. The Folgar-Tucker equation is solved for cells with ψ > 0.99 avoiding anisotropic effects in cells, which could neither be clearly assigned to the air phase nor to the suspension. In addition, the Folgar-Tucker equa-tion is solved in every time step before solving the system of governing equations (5) with the coefficients β 1 , . . . , β 6 depending on the invariants of N as given in Chung and Kwon [34] and the abbreviation S = N 2 or equivalently S i j = N ik N k j . Both terms of N[D] impl were linearized (N n ·D n )N n+1 and (S n ·D n )N n+1 as suggested by Heinen [96], with n representing the time step index.

B.1 Fiber orientation
In addition to Sect. 4.1, Figs. 14 and 15 address the local fiber orientation difference with respect to viscosity modeling. Figure 14 shows the difference between the non-Newtonian and Newtonian matrix fluid for the decoupled case, whereas Fig. 15 address the same difference for the coupled case. The fields in both figures look similar. This is because for the decoupled case, the non-Newtonian fluid has a flatter velocity profile than the Newtonian fluid. This phenomenon is exactly the same in the coupled case, in which the flowfiber coupling flattens the velocity profile (see, e.g., [5,24]).
In the coupled case, however, the non-Newtonian velocity profile flattens more than in the Newtonian case, which is why the characteristics of the difference field correspond to the decoupled approach.  The differences regarding the two matrix fluid behaviors are caused by different shear rates. Especially in the channel center region, flattening causes less shear leading to less fiber alignment in the main shear direction. This phenomenon can be seen for N 11 in the horizontal channel, where the Newtonian fluid results in larger values for N 11 close to the center line compared to the non-Newtonian case (blue color). The same can be seen for N 22 in inside the right rib. Furthermore, the difference fields vanish close to the channel walls, since the asymptotic orientation state being dominant there only depends on the kinematics, fiber aspect ratio α and the fiberfiber interaction strength C I and not on the fluid model of the matrix. For further information regarding asymptotic states of the Folgar-Tucker equation, the reader is referred to Karl et al. [31].
The differences between the viscosity models addressed above depend on the shear-thinning behavior of the matrix fluid in the non-Newtonian case. In order quantify this, Fig. 16a shows the non-dimensional viscosity μ(γ M )/μ 0 based on the Cross-model given in Eq. (9). In addition, the field of the effective matrix shear rateγ M is shown in Fig. 16b with respect to the coupled non-Newtonian simulation during mold-filling (t = 0.5 s). It can be seen that the shear rate is locally so large that significant shear-thinning is present.

B.2 Stress field
In addition to the results of Sect. 4.2 the local effective stress field is given in Fig. 17. The left column refers to the stress analysis for the decoupled simulation and the right column refers to the coupled approach, both for non-Newtonian matrix behavior during mold filling. Since the elastic limit of the material in use (PA-GF40) is 60.2-87.9 MPa [97], the deformation is purely elastic.
The non-dimensional stress difference field regarding both viscosity models for the decoupled approach is given in Fig. 18. For the coupled approach, this difference field is shown in Fig. 19. Due to the prescribed boundary conditions (see Sect. 3.1), effects of the viscosity model mainly can be limited to σ * 11 inside the horizontal channel and to σ * 22 inside the ribs. These differences are directly linked to the previous section and caused by the flattening of the velocity profile. This flatting affects the fiber alignment by a decreased shear rate in the center region of the channels which influences the local stiffness. These local stiffness differences are represented by the local stress differences regarding the prescribed displacement boundary conditions.