**1. Introduction**

Fiber reinforced polymers are key materials across different industries. For example short fiber reinforced thermoplastics are wildly used in the automotive industry to reduce weight. The manufacturing processes of those materials have typically strong impact on their final microstructure, which at the same time controls the mechanical performance of the part. A reliable virtual engineering design of fiber-reinforced polymers therefore requires considering the simulation of the process-induced microstructure.

One relevant microstructure descriptor in fiber-reinforced polymers is the fiber orientation. This work focuses on the modeling of the fiber orientation phenomenon and presents a historical review of the different modeling approaches. In this context, the article describes different modeling approaches such as the addition of a scalar diffusion by Folgar and Tucker [1], the nematic potential approach [2], the modeling of a retarding rate in the reduced strain closure (RSC) [3] and the retarding principle rate (RPR) model [4], and lastly the anisotropic rotary diffusion approach (ARD) by Phelps et al. [5]. Additionally, reduced parameters models like the improved anisotropic rotatory diffusion (iARD) [4], principal anisotropic rotary diffusion (pARD) [6], and Moldflow rotary diffusion (MRD) [7] will be introduced. The mentioned models are provided in commercial injection molding software such as Autodesk Moldflow®, Moldex 3D®, Sigmasoft® and Cadmould®. For example,

Autodesk Moldflow® provides the Folgar-Tucker, RSC and ARD-RSC model and Moldex 3D® the Folgar-Tucker, ARD, and iARD-RPR model.

Furthermore, we briefly discuss closure approximations, which are a common mathematical element of those macroscopic fiber orientation models. Simple closure approximation like the linear, quadratic, and hybrid closure [8,9] will be introduced. Then, we focus on exact closure approximations and fitted closure approximations. In the field of fitted closures we distinguish between orthotropic and invariant based closures.

Afterwards, we introduce micro-scale numerical methods for simulating the fiber orientation phenomenon, such as the discrete element method (DEM), smoothed particle hydrodynamics (SPH) and moving particle semi- implicit MPS. The focus will be on DEM based method and existing approaches will be looked at under the points of fiber discretization, imposed flow fields, fluid-fiber interaction, and fiber-fiber interaction.

The last section focuses on combining the advantages of both scales, for example through parameter fitting or machine learning approaches.

#### **2. Fiber Orientation**

The orientation of a single fiber can be characterized by a unit vector **p** ∈ S<sup>3</sup> := {**p** ∈ R<sup>3</sup> : ||**p**|| = 1} along the fiber axis. All **p** ∈ S<sup>3</sup> can be defined by two angles (*φ*, *θ*) (Figure 1)

$$p\_1 \quad = \quad \sin \theta \cos \phi \tag{1}$$

$$p\_2 = -\sin\theta\sin\phi\tag{2}$$

$$p\_3 = -\cos\theta \tag{3}$$

**Figure 1.** Orientation of a single rigid fiber **p**.

To describe the orientation of many fibers statistical methods are useful. The probability density function (PDF) *ψ*() [1] is defined, such that *ψ*(**<sup>p</sup>**, *<sup>t</sup>*)<sup>d</sup>**p** is the probability that a fiber is directed between **p** and **p** + d**p** at time *t*. The PDF has the following distinct properties.

$$\mathcal{B}(\psi) \quad = \quad [0,1] \tag{4}$$

$$\oint \psi(\mathbf{p},t)d\mathbf{p} \quad = \quad 1\tag{5}$$

$$
\psi(\mathbf{p}) \quad = \quad \psi(-\mathbf{p}) \tag{6}
$$

$$\frac{D\psi}{Dt} = -\nabla\_s \cdot (\psi \dot{\mathbf{p}}).\tag{7}$$

Equation (6) is valid under the assumption that the fibers are cylindrical and have no preferred end. Equation (7) is called the continuity equation. The operator ∇*s* represent the gradient operator on the surface of the unit sphere.

Since the PDF is defined on the unit sphere, computation is expensive and numerically difficult. For that reason, moments of the PDF are commonly used for computational efficiency [8]. The orientation tensors are defined by

$$a\_{ij} = \int p\_i p\_j \psi(\mathbf{p}) d\mathbf{p} \tag{8}$$

$$a\_{i\bar{j}kl} = \int\_{\cdot} p\_i p\_j p\_k p\_l \psi(\mathbf{p}) d\mathbf{p}.\tag{9}$$

$$a\_{i\dots n} = \int p\_i \cdot \cdots \cdot p\_n \psi(\mathbf{p}) d\mathbf{p}.\tag{10}$$

Since the PDF is even (Equation (6)) all odd-ordered orientation tensors are zero. To simplify notation we introduce **A** = *aij*, A = *aijkl*. The orientation tensors have important properties [8], which are described here for the second and fourth moment only. The second order tensor is symmetric and has a unit trace.

$$\mathbf{A}\_{\ddot{\imath}\ddot{\jmath}} = \mathbf{A}\_{\ddot{\jmath}\ddot{}}\tag{11}$$

$$\text{tr}\mathbf{A} \; \; \; \;=\; \; 1. \tag{12}$$

The fourth order tensor is symmetric with respect to any pair of indices.

$$\mathbf{A}\_{\rm ijkl} = \mathbf{A}\_{\rm jikl} = \mathbf{A}\_{\rm ijkl} = \mathbf{A}\_{\rm kjil} = \mathbf{A}\_{\rm ljki} = \mathbf{A}\_{\rm ikjl} = \mathbf{A}\_{\rm ilkj} \tag{13}$$

and all information of the second order orientation tensor can be retrieved from the fourth order tensor

$$\mathbf{A}\_{ij} = \mathbf{A}\_{ijkk}.\tag{14}$$

The use of orientation tensors simplifies the computation because no discretization of the unit sphere is necessary, but it is impossible to distinguish between certain orientations. For example, 

the orientation tensor for a bipolar and planar random orientation is identical with **A** = ⎛⎜⎝ 0.5 0 0 0 0.5 0 0 00 ⎞⎟⎠.

An advantage is the objectivity, that means the equation is independent of the coordinate system.

#### **3. Macroscopic Fiber Orientation Models**

Macroscopic fiber orientation models are used to predict fiber orientation in parts. The models are integrated in commercial software such as Autodesk Moldflow®, Moldex 3D®, Sigmasoft® and Cadmould®. Additionally open source software such as OpenFoam can be used to implement the models.

#### *3.1. Macroscopic Fiber Orientation Models in the Dilute Regime*

The first description of the motion of a single fiber was developed by Jeffrey [10]. This description is based on the following assumptions: The fluid is Newtonian and has no turbulences (rot **u** = 0). The particle is an ellipsoid and perfectly rigid, such that no bending or breaking occurs. The velocity field of the fluid is not influenced by the particle and there exists a perfect contact between the particle and the fluid. Under these assumptions the model is accurate up to order two. This has been proven in a more general way by Junk and Illner [11] under the same assumptions. Jeffrey's model has the form:

$$
\dot{\mathbf{p}} = \mathbf{W} \cdot \mathbf{p} + \tilde{\boldsymbol{\zeta}} (\mathbf{D} \cdot \mathbf{p} - \mathbf{D} : \mathbf{ppp}),
\tag{15}
$$

where **W** = 1 2 (( <sup>∇</sup>**u**)*<sup>T</sup>* − ∇**u**) is the vorticity tensor, **D** = 1 2 (( <sup>∇</sup>**u**)*<sup>T</sup>* + ∇**u**) the rate of strain tensor, *ξ* = *r*2 *e* −1 *r*2 *e* +1 the particle shape function and ∇ the nabla operator. Applying it to injection molding simulation, it has to be highlighted that the polymer melt is non-Newtonian and fibers are not perfectly rigid. In fact they can break and bend. Furthermore for highly filled polymers, the fluid is influenced by the fiber [2]. Coupling between fluid and fiber orientation will not be considered in this review.

#### *3.2. Macroscopic Fiber Orientation Models in the Concentrated Regime*

In the concentrated regime fiber interaction is dominant for fiber orientation. All macroscopic modeling approaches published up to today, accounting for interaction, are phenomenological.

Folgar and Tucker [1] added a diffusion term to account for fiber interaction to predict the orientation in semi-dilute and concentrated solutions. The model is valid under the following additional assumptions: The fibers are rigid cylinders, uniform in length and diameter. Moreover, they are sufficient large such that Brownian motion is negligible. The matrix is incompressible and sufficient viscous, such that particle inertia and buoyancy is negligible. The center of mass of the fibers are randomly distributed and no external forces or torques act on the suspension.

The fiber motion of a single fiber can then be described by

$$\dot{\mathbf{p}} = \mathbf{W} \cdot \mathbf{p} + \zeta (\mathbf{D} \cdot \mathbf{p} - \mathbf{D} : \mathbf{ppp}) + C\_l \dot{\gamma} \cdot \nabla\_s \psi,\tag{16}$$

where *CI* describes the interaction coefficient and *γ*˙ = √2**D** : **D** the scalar magnitude of the rate of strain tensor. The equation is, in contrast to Jeffrey's equation (15), not reversible. The interaction coefficient *CI* is empirically determined and describes the rate of interaction. Setting *CI* = 0 retrieves Jeffrey's equation (15). Interactions of fibers cause random orientation. The Fokker-Planck equation (17) [1] expresses the rate of change for the PDF, using the Folgar-Tucker equation for a single fiber (16) and the continuity equation (7). Fibers are modeled as independent, identically distributed, random variables with zero mean. Each interaction causes an orientation change in both fibers

$$\begin{split} \frac{D\boldsymbol{\psi}}{Dt} &= \quad -\nabla\_{\boldsymbol{s}} \cdot \left( \boldsymbol{\psi} (\mathbf{W} \cdot \mathbf{p} + \boldsymbol{\xi} (\mathbf{D} \cdot \mathbf{p} - \mathbf{D} : \mathbf{pp} \mathbf{p}) + \mathbf{C}\_{I} \boldsymbol{\gamma} \cdot \nabla\_{\boldsymbol{s}} \boldsymbol{\psi} \right) \\ &= \quad -\nabla\_{\boldsymbol{s}} \cdot \left( \boldsymbol{\psi} (\mathbf{W} \cdot \mathbf{p} + \boldsymbol{\xi} (\mathbf{D} \cdot \mathbf{p} - \mathbf{D} : \mathbf{pp} \mathbf{p})) \right) + \mathbf{C}\_{I} \boldsymbol{\gamma} \cdot \nabla\_{\boldsymbol{s}}^{2} \boldsymbol{\psi} .\end{split} \tag{17}$$

The concept introduced by Folgar and Tucker offers advantages in the prediction of fiber orientation but is difficult to solve numerically and can only be solved with high computational effort. Based on the Fokker-Planck equation (17), Advani and Tucker [8] developed an equation for the rate of change of the second order orientation tensor

$$\frac{D\mathbf{A}}{Dt} = \dot{\mathbf{A}}\quad =\quad \dot{\mathbf{A}}^h + \dot{\mathbf{A}}^d\tag{18}$$

$$\mathbf{A}^{\mathrm{h}} = \begin{pmatrix} \mathbf{W} \cdot \mathbf{A} - \mathbf{A} \cdot \mathbf{W} \end{pmatrix} + \prescript{\mathrm{\pi}}{\zeta} (\mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2\mathbf{A} : \mathbf{D}) \tag{19}$$

$$\mathbf{A}^d \quad = \ 2C\_I \dot{\gamma} (\mathbf{I} - 3\mathbf{A}) \dots \tag{20}$$

The influence of the phenomenological parameter *CI* is displayed in Figure 2. With decreasing diffusion (small *CI*) fibers are more aligned.

Evolution equations of the second order orientation tensor **A** contain the fourth order fiber orientation tensor A. It is possible to derive an equation for the fourth order orientation tensor, but it will contain the sixth order orientation tensor. This leads to an infinite series of evolution equations. Therefore, it is common to truncate the series at the level of the second order orientation tensor. Since the fourth order orientation tensor is then not explicitly computed, a so-called closure approximation is necessary for calculating the fourth order fiber orientation tensor. In Section 3.3 different closure approximations will be introduced.

˙

**Figure 2.** Influence of the phenomenological parameter *CI* on fiber orientation evolution.

Three model enhancements have been made since the simulated fiber orientation shows deviation to fiber orientation determined experimentally.

To slow down the evolution speed, Huynh [12] introduced the strain reduction factor (SRF) 1*α* by multiplying equation (18) with *α* ∈ [0, 1]

$$
\mathbf{A} = \boldsymbol{\kappa} (\mathbf{A}^{\mathbf{h}} + \mathbf{A}^{\mathbf{d}}).\tag{21}
$$

The SRF model violates the material objectivity, so Wang et al. [3] developed the reduced strain closure (RSC), based on the eigenvalue *λ* and eigenvector **e** decomposition, that is, **A** = ∑<sup>3</sup>*i*=<sup>1</sup> *λi***e***i***e***i* with orthonormal eigenvectors. The modified growth rate has the following form

$$
\dot{\lambda}\_i^{\mathsf{RSC}} = \kappa \dot{\lambda}\_i \tag{22}
$$

$$
\dot{\mathbf{e}}\_{i}^{\mathsf{RSC}} = \quad \dot{\mathbf{e}}\_{i\prime} \tag{23}
$$

with the constant *κ* ∈ [0, 1]. In the differential equation form, this equals

$$\mathbf{A} \quad = \mathbf{A}^{R \mathbf{S} \mathbf{C}} + \kappa \mathbf{A}^{d} \tag{24}$$

$$\mathbf{A}^{\mathsf{RSC}} = \underset{\cdot}{\mathbf{W} \cdot \mathbf{A}} - \mathbf{A} \cdot \mathbf{W} + \mathfrak{J} \left[ \mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2(\mathbf{A} + (1 - \kappa)(\mathbb{L} - \mathbf{M} : \mathbf{A})) : \mathbf{D} \right] \tag{25}$$

$$\begin{array}{rcl} \mathbf{L} &=& \sum\_{i=1}^{3} \lambda\_{i} \mathbf{e}\_{i} \mathbf{e}\_{i} \mathbf{e}\_{i} \mathbf{e}\_{i} \\\\ \end{array} \tag{26}$$

$$\|M\| = \sum\_{i=1}^{3} \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i \mathbf{e}\_i. \tag{27}$$

Figure 3 displays the influence of the phenomenological constant *κ* with a smaller retarding rate (smaller *κ*) a slower orientation evolution is reached. The steady state is unchanged by *κ*. Tseng et al. [4] introduced at the same time the retarding principal rate (RPR), modifying the growth rate of the eigenvalue by

$$
\dot{\lambda}\_i^{\text{RPR}} = \kappa \dot{\lambda}\_i + (1 - \kappa) \beta (\dot{\lambda}\_i^2 + \dot{\lambda}\_j \dot{\lambda}\_k) \tag{28}
$$

$$
\dot{\mathbf{e}}\_{i}^{\mathsf{RPR}} = \quad \dot{\mathbf{e}}\_{i\prime} \tag{29}
$$

with a fine-tuning parameter *β*. For *β* = 0 this equals the RSC model equation (22). The differential equation form (similar to (24) with PRR) can be developed based on equations (28) and (29).

**Figure 3.** Influence of the phenomenological parameter *κ* on fiber orientation evolution. The diffusion constant is set to *CI* = 0.01.

The second model approach, added by Phelps and Tucker [5], introduced an anisotropic rotary diffusion (ARD) term by replacing *C*I with a rotary diffusion tensor **C**. This allows spatially non-uniform rotary diffusion, which makes the rotary diffusion effect a function of the orientation state

$$
\dot{\mathbf{A}}^{\prime} = \dot{\mathbf{A}}^{\prime} + \dot{\mathbf{A}}^{\prime \text{ARD}} \tag{30}
$$

$$\mathbf{A}^{\text{ARD}} \quad = \quad \dot{\gamma} [2\mathbf{C} - 2\text{tr}(\mathbf{C})\mathbf{A} - 5(\mathbf{C} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{C}) + 10\mathbf{A} : \mathbf{C}].\tag{31}$$

The different modeling approaches for the anisotropic rotary diffusion are listed below

$$\mathbf{C} = \mathbf{C}(\mathbf{D}, \mathbf{A}) = b\_1 \mathbf{I} + b\_2 \mathbf{A} + b\_3 \mathbf{A}^2 + \frac{b\_4}{\dot{\gamma}} \mathbf{D} + \frac{b\_5}{\dot{\gamma}^2} \mathbf{D}^2. \tag{32}$$

ARD model [5] with constants *bi*, *i* = 1, ..., 5

$$\mathbf{C} = \mathbf{C}\_{\text{l}} \left( \mathbf{I} - 4 \mathbf{C}\_{\text{M}} \frac{\mathbf{D}^2}{\dot{\gamma}} \right) \tag{33}$$

iARD model [4] with constants *C*I, *CM*

$$\mathbf{C} = \mathbf{C}\_{\text{I}} \mathbf{R}\_{\text{A}} \begin{pmatrix} D\_1 & 0 & 0 \\ 0 & D\_2 & 0 \\ 0 & 0 & D\_3 \end{pmatrix} \mathbf{R}\_{\text{A}}^{\top} \tag{34}$$

pARD model [6] with constants *C*I, *D*1, *D*2, *D*3.

> The eigenmatrix **R**A is defined by

$$\mathbf{A} = \mathbf{R}\_{\mathsf{A}} \begin{pmatrix} \lambda\_1 & 0 & 0 \\ 0 & \lambda\_2 & 0 \\ 0 & 0 & \lambda\_3 \end{pmatrix} \mathbf{R}\_{\mathsf{A}}^{\top}. \tag{35}$$

Tseng et al. [6] chose

$$D\_1 = 1,\\ D\_2 = c,\\ D\_3 = 1 - c \tag{36}$$

to reduce the needed amount of parameters. This implies that the rotary diffusion factor in the first principal fiber orientation direction is defined by *CI*. The second principal fiber orientation direction is scaled with *<sup>c</sup>*·*CI* and in the third principal fiber orientation direction the smallest rotary diffusion factor with (1 − *c*)*CI* is applied. The influence od *D*2 is displayed in Figure 4.

**Figure 4.** Influence of the phenomenological parameter *D*2 on fiber orientation evolution. The diffusion is et to *CI* = 0.01.

A slightly different approach to add anisotropic rotary diffusion is defined in the MRD model [7]. The rotary diffusion is modeled according to equation (34), but not the full ARD equation (31) is used. The reduced form is given by

$$
\dot{\mathbf{A}}^{\quad} = \dot{\mathbf{A}}^{\text{h}} + \dot{\mathbf{A}}^{\text{mVARD}} \tag{37}
$$

$$\mathbf{A}^{\mathsf{mARD}} \quad = \ \ \mathbf{\hat{}}\dot{\gamma}(\mathbf{C} - \text{tr}(\mathbf{C})\mathbf{A}).\tag{38}$$

Recently, Favaloro and Tucker [13] published a framework to compare anisotropic rotary diffusion approaches and compared the mentioned approaches in shear and elongation flows. They added a suggestion for a more general but stable model by using equation (34) but making *Di* a function of *λi* or using equation (32) with *b*4 = *b*5 = 0 and *bi*, *i* = 1, 2, 3 as scalar functions of tr**A**<sup>2</sup> and tr**A**3. Latz et al. [2] added a nematic potential to the Folgar-Tucker model to account for excluded volume effects. So far, there is no clear dependency of *CI* on the volume fraction and aspect ratio of the fillers. Different experiments in different regimes even showed contradictory results [14–16]. Latz et al. [2] stated that the excluded volume mechanism may be a possible explanation for the observed effects. The integration of a second constant *U*0 could decouple this effects and perhaps give a clear dependence onphysicaldescriptors.Themodelisdefinedby

$$\mathbf{A}^{\hspace{1cm}} = \mathbf{A}^{\hspace{1cm}} + \mathbf{A}^{\text{nem}} \tag{39}$$

$$\dot{\mathbf{A}}^{\text{new}} = \left( \mathbb{C}\_{l} (\mathbf{I} - 3\mathbf{A}) + lL\_{0} (\mathbf{A} \cdot \mathbf{A} - \mathbf{A} : \mathbb{A}) \right) \dot{\gamma},\tag{40}$$

with one additional constant *U*0.

In injection molding simulation combination of the methods are used. Combining Equations (24) and (31) leads to the following (p)ARD-RSC model

$$\mathbf{A} \quad = \quad \mathbf{A}^{RSC} + \mathbf{A}^{(p)ARD - RSC} \tag{41}$$

$$\begin{split} \mathbf{A}^{(p)ARD-RSC} &= \begin{aligned} ^{\prime}\dot{\mathbf{c}} [\mathbf{C} - (1 - \kappa)\mathbf{M} : \mathbf{C}] - 2\mathbf{x}(\mathbf{trC})\mathbf{A} - \mathbf{5}(\mathbf{C} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{C}) \\ + 10[\mathbf{A} - (1 - \kappa)(\mathbb{L} - \mathbf{M} : \mathbb{A})] : \mathbf{C}]. \end{aligned} \tag{42}$$

#### *3.3. Closure Approximations*

Since the fourth order orientation tensor is used in all models, closure approximations are needed. Any fourth order tensors with the symmetric properties stated in equation (13) can be represented by a 6 × 6 matrix with at most 36 independent entries [17]

$$
\Delta = \begin{pmatrix}
\mathbb{A}\_{11} & \mathbb{A}\_{12} & \mathbb{A}\_{13} & \mathbb{A}\_{14} & \mathbb{A}\_{15} & \mathbb{A}\_{16} \\
\mathbb{A}\_{21} & \mathbb{A}\_{22} & \mathbb{A}\_{23} & \mathbb{A}\_{24} & \mathbb{A}\_{25} & \mathbb{A}\_{26} \\
\mathbb{A}\_{31} & \mathbb{A}\_{32} & \mathbb{A}\_{33} & \mathbb{A}\_{34} & \mathbb{A}\_{35} & \mathbb{A}\_{36} \\
\mathbb{A}\_{41} & \mathbb{A}\_{42} & \mathbb{A}\_{43} & \mathbb{A}\_{44} & \mathbb{A}\_{45} & \mathbb{A}\_{46} \\
\mathbb{A}\_{51} & \mathbb{A}\_{52} & \mathbb{A}\_{53} & \mathbb{A}\_{54} & \mathbb{A}\_{55} & \mathbb{A}\_{56} \\
\mathbb{A}\_{61} & \mathbb{A}\_{62} & \mathbb{A}\_{63} & \mathbb{A}\_{64} & \mathbb{A}\_{65} & \mathbb{A}\_{66}
\end{pmatrix}, \tag{43}
$$

where <sup>A</sup>*αβ* = <sup>A</sup>*ijkl* for *α*, *β* ∈ {1, 2, 3, 4, 5, 6} the indices *α* and *β* represent an index pair *ij* or *kl* in the following way:

$$a, \beta: \begin{cases} 1 & \to & 11 \\ 2 & \to & 22 \\ 3 & \to & 33 \\ 4 & \to & 23 \\ 5 & \to & 31 \\ 6 & \to & 12 \end{cases} \tag{44}$$

The first simple closure approaches have been introduced by Advani and Tucker [8]. The linear approach [8] is a summation of all products of *aij* and *<sup>δ</sup>ij*. After applying symmetry and normalization condition the following linear approximation occurs

$$\begin{split} \mathsf{A}\_{ijkl}^{\mathrm{lin}} &= -\frac{1}{35} (\delta\_{ij}\delta\_{kl} + \delta\_{ik}\delta\_{jl} + \delta\_{il}\delta\_{jk}) \\ &+ \frac{1}{7} (\mathbf{A}\_{ij}\delta\_{kl} + \mathbf{A}\_{ik}\delta\_{jl} + \mathbf{A}\_{il}\delta\_{jk} + \mathbf{A}\_{kl}\delta\_{ij} + \mathbf{A}\_{jl}\delta\_{ik} + \mathbf{A}\_{jk}\delta\_{il}). \end{split} \tag{45}$$

The quadratic closure [8] omits all linear terms

$$\mathbf{A}\_{ijkl}^{quad} = \mathbf{A}\_{ij}\mathbf{A}\_{kl}.\tag{46}$$

The quadratic closure does not preserve the symmetry of the fourth order tensor (equation (8)), but applying it to equation (19), it does preserve the symmetry of the second order tensor.

It is also possible to combine both approaches. This leads to the hybrid closure [8]

$$\mathbb{A}\_{ijkl}^{l \text{j} \text{j}} \; \;= \; \; (1 - f)\mathbb{A}\_{ijkl}^{l \text{in}} + f\mathbb{A}\_{ijkl}^{q \text{mid}} \tag{47}$$

$$f\_{\parallel} = \frac{3}{2} \mathbf{A}\_{ij} \mathbf{A}\_{\parallel i} - \frac{1}{2}. \tag{48}$$

Another approach to determine the function for the hybrid closure was determined by Advani and Tucker [9]

$$f = 1 - 2\mathcal{T} \cdot \det \mathbf{A}.\tag{49}$$

A more advanced approach is the natural closure (NAT). It is based on the relationship A = *f*(**A**) when no diffusion occurs (*CI* = 0).

Verley and Dupret [18] stated that there exists an exact closure in the case that the orientation is at one time isotropic and used a numerical approximation in the 3-D case to obtain manageable computational cost.

Montgomery-Smith et al. [19] determined the exact formulation in Cartesian coordinates on the sphere using the Carlson form of elliptic integrals. They used the analytic solution of the Jeffrey equation presented by Dinh and Armstrong [20]. The approach is based on the assumption of isotropic orientation at *t* = 0. Then the second order orientation tensor can be expressed by

$$\mathbf{A} = \int\_{S} \frac{\mathbf{pp}}{4\pi (\mathbf{B} : \mathbf{pp})^{\frac{3}{2}}} \mathrm{d}\mathbf{p} \tag{50}$$

$$\mathbf{B}\_{\perp} = -\mathbf{B} \cdot \left( \dot{\mathbf{W}} + \ddot{\xi} \mathbf{D} \right) - \left( -\mathbf{W} + \xi \mathbf{D} \right) \cdot \mathbf{B}\_{\prime} \ \mathbf{B} = \mathbf{I} \text{ at } t = 0. \tag{51}$$

The method presented uses high computational effort, so Montgomery-Smith et al. [19] introduced the fast exact closure (FEC). The FEC introduced a computationally efficient way to compute the closure. Instead of computing **B** by inverting the integral, equations (18) and (51) are solved simultaneously. If the initial data is not isotropic, **B** has to be computed for the initial condition.

Later Montgomery-Smith et al. [21] also include the anisotropic rotary diffusion of Phelps and Tucker [5]. In this case (**C** = 0) the closure is not exact. The key idea is to introduce a matrix B and define two conversion tensors C and D such that

$$\frac{D\mathbf{A}}{Dt} = -\mathbb{C} : \frac{D\mathbf{B}}{Dt}, \frac{D\mathbf{B}}{Dt} = -\mathbb{D} : \frac{D\mathbf{A}}{Dt}.\tag{52}$$

hold true. The ordinary differential equations (ODEs) are solved simultaneously, which can be computed very efficiently. The special form used for this approach is

$$\frac{D\mathbf{A}}{Dt} = -\mathbb{C} : F(\mathbf{B}) + G(\mathbf{A}),\\\frac{D\mathbf{B}}{Dt} = F(\mathbf{B}) - \mathbb{D} : G(\mathbf{A}).\tag{53}$$

They showed that it can also by applied to the reduced strain closure in equation (24) and proved that the solution stays physical for example, that **A** stays positive definite with tr **A** = 1.

A large family of closure are fitted closures. Chung and Kwon [22] stated two ways to develop fitted closure. The method depends on the coordinate system, global or eigenspace, and the representation of the approximation, invariants or eigenvalues, of the orientation tensor. The orthotropic closures use the eigenspace coordinate system and eigenvalues as representatives. They are the most used family of closures. Cintra and Tucker [23] stated that an objective closure has to be orthotropic. Orthotropic means that the principal axes must match those of the second order tensor and each principal fourth order component is a function of just two principal values of the second order tensor.

The second order orientation tensor can be transformed in the eigenraum representation

$$
\hat{\mathbf{A}} = \begin{pmatrix} \lambda\_1 & 0 & 0 \\ 0 & \lambda\_2 & 0 \\ 0 & 0 & \lambda\_3 \end{pmatrix} = \mathbf{R}\_{\mathsf{A}}^T \mathbf{A} \mathbf{R}\_{\mathsf{A}\prime} \tag{54}
$$

with non-negative eigenvalues *λ*1 ≥ *λ*2 ≥ *λ*3 ≥ 0 and *λ*1 + *λ*2 + *λ*3 = 1. Figure 5 shows the reference triangle for orthotropic closures.

Kuzmin [17] showed that regarding all symmetries of the fourth order orientation tensor the orthotropic representation has the form

$$
\hat{\underline{\dot{\Delta}}} = \begin{pmatrix}
\underline{\hat{\underline{\dot{\Delta}}}}\_{11} & \underline{\hat{\underline{\dot{\Delta}}}}\_{112} & \underline{\hat{\underline{\dot{\Delta}}}}\_{13} & 0 & 0 & 0 \\
\underline{\hat{\underline{\dot{\Delta}}}}\_{12} & \underline{\hat{\underline{\dot{\Delta}}}}\_{22} & \underline{\hat{\underline{\dot{\Delta}}}}\_{23} & 0 & 0 & 0 \\
\underline{\hat{\underline{\dot{\Delta}}}}\_{13} & \underline{\hat{\underline{\dot{\Delta}}}}\_{23} & \underline{\hat{\underline{\dot{\Delta}}}}\_{33} & 0 & 0 & 0 \\
0 & 0 & 0 & \underline{\hat{\underline{\dot{\Delta}}}}\_{44} & 0 & 0 \\
0 & 0 & 0 & 0 & \underline{\hat{\underline{\dot{\Delta}}}}\_{55} & 0 \\
0 & 0 & 0 & 0 & 0 & \underline{\hat{\underline{\dot{\Delta}}}}\_{66}
\end{pmatrix}.\tag{55}
$$

It has been proven by Cintra and Tucker [23] and Kuzmin [17] that the fourth order orientation tensor can be represented by three independent components. This is due to the orthotropic properties, the symmetries of the second order tensor, the symmetry with respect to any pair of indices for the fourth order tensor and the normalization condition. The three components can be expressed as functions of the two largest eigenvalues of the second order tensor

$$
\hat{\underline{A}}\_{11} = \begin{array}{c} f\_1(\lambda\_1, \lambda\_2) \\ \end{array} \tag{56}
$$

$$\underline{\mathbb{A}}2^{\underline{\mathbf{a}}} = \,\_1f\_2(\lambda\_1, \lambda\_2) \tag{57}$$

$$
\hat{\underline{A}}\_{33} \quad = \quad f\_3(\lambda\_1, \lambda\_2). \tag{58}
$$

The remaining entries are defined by symmetry and normalization condition

ˆ

$$
\begin{array}{rcl}
\underline{\mathring{A}}\_{12} &=& \underline{\mathring{A}}\_{66} \\
\underline{\mathring{A}}\_{2} &=& \underline{\mathring{A}}\_{66}
\end{array}
\tag{59}
$$

$$\begin{array}{rcl} \underline{\mathbb{A}}\_{23} &=& \underline{\mathbb{A}}\_{44} \\ \underline{\mathbb{A}}\_{13} &=& \underline{\mathbb{A}}\_{55} \end{array} \tag{60}$$

$$
\hat{\mathbb{A}}\_{\mathsf{55}} + \hat{\mathbb{A}}\_{\mathsf{66}} = \begin{array}{c} \mathsf{42} \\ \end{array}\_{\mathsf{66}} = \begin{array}{c} \mathsf{42} \\ \end{array}\_{\mathsf{64}} \tag{62}
$$

$$
\underline{\mathbb{A}}\_{44} + \underline{\mathbb{A}}\_{66} = \lambda\_2 - \underline{\mathbb{A}}\_{22} \tag{63}
$$

$$
\underline{\hat{A}}\_{44} + \underline{\hat{A}}\_{55} = \begin{array}{c} \lambda\_3 - \underline{\hat{A}}\_{33}. \end{array} \tag{64}
$$

Cintra and Tucker [23] used the eigenspace coordinate system and eigenvalues as representatives to fit the closure (EBOF). Complete second order polynomials were used to approximate the components (ORF). The orthotropic fitted (ORF) closure and the NAT closure are based on identical assumptions. In conclusion they can be transferred to each other. In fact they are mathematically equivalent, but the ORF closure is numerically more stable for repeated eigenvalues [19] The three remaining components are fitted by

=

$$\hat{\Delta}\_{\text{iii}} = f\_i(\lambda\_1, \lambda\_2) = \mathbb{C}\_i^1 + \mathbb{C}\_i^2 \lambda\_1 + \mathbb{C}\_i^3 \lambda\_1^2 + \mathbb{C}\_i^4 \lambda\_2 + \mathbb{C}\_i^5 \lambda\_2^2 + \mathbb{C}\_I^6 \lambda\_1 \lambda\_2 \tag{65}$$

for for *i* = 1, 2, 3.

The components are fitted with different flow types (simple shear, two shearing/stretching flows, uniaxial elongation, biaxial elongation) using a least squares routine. The closure approximation shows good agreemen<sup>t</sup> with the distribution function and the NAT and better performance than the hybrid closure. For small *CI* the ORF closure shows oscillating behavior.

Chung and Kwon [24] improved the ORF closure to overcome the oscillation for small *CI* values and introduced the orthotropic fitted closure approximation for wide interaction coefficients (OWE) and orthotropic fitted closure approximation for wide interaction coefficients with third order polynomial approximations (OWE3) closure. The OWE closure is fitted with two additional flow fields (shear/planar elongation, balanced shear/biaxial elongation) to cover the orientation triangle *K* ˆ (Figure 5) more closely. Consequently the only difference between the ORF and OWE closure are the values of the fitted parameters. The OWE3 closure uses the same flows for the parameter fit but approximates the coefficients values by a third order polynomial expression

$$\begin{aligned} \hat{\Delta}\_{ii} &= \; f\_i(\lambda\_1, \lambda\_2) = \mathcal{C}\_i^1 + \mathcal{C}\_i^2 \lambda\_1 + \mathcal{C}\_i^3 \lambda\_1^2 + \mathcal{C}\_i^4 \lambda\_2 + \mathcal{C}\_i^5 \lambda\_2^2 + \mathcal{C}\_i^6 \lambda\_1 \lambda\_2 \\ &+ \mathcal{C}\_i^7 \lambda\_1^2 \lambda\_2 + \mathcal{C}\_i^8 \lambda\_1 \lambda\_2^2 + \mathcal{C}\_i^9 \lambda\_1^3 + \mathcal{C}\_i^{10} \lambda\_2^3 \end{aligned} \tag{66}$$

for for *i* = 1, 2, 3. The closures show stable behavior for a wide range of *CI* but tend to over predict the orientation.

**Figure 5.** Reference triangle for orthotropic closure.

Kuzmin [17] introduced a mathematical concept to develop orthotropic closures. A concept for planar orientations was developed and extended to the 3D case. In this work only the 3D case is explained, for the planar case refer to Reference [17] The linear and smooth closures were stated in the orthotropic state [17]. An orthotropic version of the quadratic closure was developed

$$
\underline{\hat{A}}{}\_{11} \quad = \quad f\_1(\lambda\_1 \lambda\_2) = \lambda\_1^2 \tag{67}
$$

$$
\hat{\underline{A}}{22}\_2 = \begin{array}{c} f\_2(\lambda\_1, \lambda\_2) = \lambda\_2^2 \end{array} \tag{68}
$$

$$\underline{\mathbb{A}}\_{33} \quad = \; \, f\_{\mathbb{S}}(\lambda\_1 \lambda\_2) = (1 - \lambda\_1 - \lambda\_2)^2. \tag{69}$$

Since there does not exist an analytical form of the standard NAT, natural closures based on extended quadratic and piecewise linear interpolation have been developed.

The extended quadratic fit is fitted on cubic polynomials of the form

$$\begin{aligned} \underline{\hat{A}}\_{ii} &= \; f\_i(\lambda\_1, \lambda\_2) = \mathcal{C}\_i^1 + \mathcal{C}\_i^2 \lambda\_1 + \mathcal{C}\_i^3 \lambda\_1^2 + \mathcal{C}\_i^4 \lambda\_2 + \mathcal{C}\_i^5 \lambda\_2^2 + \mathcal{C}\_I^6 \lambda\_1 \lambda\_2 \\ &+ \mathcal{C}\_i^7 \lambda\_1 \lambda\_2 (1 - \lambda\_1 - \lambda\_2). \end{aligned} \tag{70}$$

The data points *Ui*, *Bij*, *i* = 1, 2, 3 *j* = *i* + 1 for *i* = 1, 2, *j* = 1 for *i* = 3 and *T* are used, using the planar natural closure and the triaxial orientation state at T. For extrapolation the extended quadratic finite element method is used.

For the exact midpoint fit quadratic interpolation polynomials are used

$$\underline{\mathbb{A}}\_{il} = f\_i(\lambda\_1, \lambda\_2) = \mathbb{C}\_i^1 + \mathbb{C}\_i^2 \lambda\_1 + \mathbb{C}\_i^3 \lambda\_1^2 + \mathbb{C}\_i^4 \lambda\_2 + \mathbb{C}\_i^5 \lambda\_2^2 + \mathbb{C}\_I^6 \lambda\_1 \lambda\_2. \tag{71}$$

The values at points *U*, *B*, *T* and the midpoints *M*1, *M*2, *M*3 are validated using the exact closure.

In contrast to all eigenspace based closures, Reference [22] introduced the invariants-based optimal fitted (IBOF) closure in the global coordinate system with invariants as representatives. This closure is computational more efficient than the orthotropic closures.

#### **4. Microscopic Fiber Simulation**

Simulation methods on the microscopic level can be used to simulate fiber movement for discretized fibers. In contrast to the phenomenological macroscopic fiber orientation models, modeling approaches on the microscopic scale enable the approximation of the physical behavior more accurately. However, they are computationally very expensive and are, in most cases, not suitable for the simulation of actual parts. Mostly, they are used on reference elements to investigate physical effects and quantities. Different modeling approaches can be used on the micro level. They can coarsely be divided in particle based methods, where the polymer matrix and the fibers are treated as particles for example the smoothed particle hydrodynamics (SPH) and moving particle semi-implicit (MPS) method, or element based methods, which treat fibers as particles and the matrix is treated as a continuous media.

One big advantage of particle based method, is the computationally cheap coupling between fluid motion and fibers in two ways. Yashiro et al. [25,26] developed a method based on MPS to predict fiber movement in injection molding. They simulated complex flow fields and investigated orientation in a T-shaped bifurcation.

The SPH method is also applied to polymer composites. He et al. [27] simulated the 3D injection molding process for short-fiber reinforced polymers; Bertevas et al. [28] simulated the 3D printing process of fiber-reinforced polymer and Yamagata and Ichimiya [29] used the method to simulate the solidification process for injection molded short-fiber reinforced parts.

The SPH method can also be combined with the discrete element method (DEM) [30], or the element bending group (EBG) method [31,32]. The fluid is then model by the SPH method and the fibers by the respective method.

Element based methods are mostly developed on the principle of the DEM. In contrast to particle based methods, which are per se two-way coupled, DEM based simulations are often solved one-way coupled. The backcoupling from the fiber motion on the fluid is computationally expensive, but can be integrated in a coarse or refined way. The resulting fluid equation can be solved by multiple approaches such as the direct numerical simulation (DNS), the Latice-Bolzmann or the particle finite element analysis (pFEA). In the DEM, fibers are considered as particles and the movement of each particle is calculated by the solution of the force and torque balance acting on each particle. Fibers are either discretized as chain of spheres, prolate spheroids or chain of rods. The forces acting on a fiber are


Hydrodynamic forces are exerted from the fluid on the fiber. The fluid motion can be considered undisturbed by the fiber motion or disturbed. In the second case a backcoupling is necessary. The interaction forces can be divided in two cases: long-range hydrodynamic interaction and short range interaction. The short range interactions can then be divided in three regimes: short range lubrication forces, transition and mechanical contact. Fibers can be modeled flexible by using chains of beads or rods connected by joints. Breaking can also be incorporated at the defined joints. The modeling of the forces, the discretization of the fibers and the fluid motion varies between the different approaches and the following Table 1 gives an overview of published approaches, without claim for completeness. Only inertia free models are considered. Examples for models incorporating inertia are References [33–35].





calculation

*J. Compos. Sci.* **2020**, *4*, 69


**Table 1.** *Cont.*



The listed approaches show, that depending on the evaluated quantities, the matrix and fiber material and the volume fraction, different modeling approaches show more promising results. In case of rigid fibers, long range hydrodynamic interactions show the highest influence in semidilute solution, whereas mechanical interaction gets dominant in the concentrated regime and long-range hydrodynamic interaction can be neglected [42,46,57]. The backcoupling has neglectable influence on rigid single fiber movement [60,75]. Once flexibility of the fiber is higher the backcoupling cannot be neglected [60]. Due to the high computational effort it has not be used in the concentrated regimen. If there is a significant influence in the concentrated regime can not be answered.

Based on the movement of single fibers, orientation evolution curves can be calculated. With the discrete position of each fiber the second order orientation tensor can be calculated in each time step by

$$A\_{ij}(t) = \sum\_{n=1}^{N} \frac{p\_{n,j}(t)p\_{n,i}(t)}{N} \, \tag{72}$$

where *N* denotes the number of fibers, *t* the actual time, and *pn*(*t*) the position of fiber *n* at time *t*.

#### **5. Using Microscopic Models for an Enhanced Prediction on the Macroscopic Scale**

A combination of two scales can enhance the prediction on the macroscopic scale, while being computationally cheap. A straight forward approach is to use the microscopic simulation for parameter definition of existing macroscopic model. This has been done by many authors, for example Reference [65,75]. Microscopic fiber orientation evolution results can also lead to new macroscopic models [67]. A different approach is to create orientation data with a microscopic simulation and use a machine learning based approach on the macroscopic level [76].
