**3. Results**

#### *3.1. Rotation and Bending in a Simple Shear Flow*

This section presents simulation results for a 3D shear flow (cf. Figure 2) with periodic boundaries in *x*1 and *x*3. The lateral dimensions of the modelled fluid domain are *B* = 1.2*L*f and the dimension in flow direction is *L* = 2*B* with the length of a fiber *L*f. The Reynolds number

$$Re = \frac{\rho\_0 L\_f \frac{GB}{2}}{\eta} \tag{42}$$

is set to *Re* = 0.5 to be small enough for a quasi-creeping flow, but also as large as possible to increase the time increment. A dimensionless measure for the fiber stiffness is

$$S = \frac{E\pi}{4\eta Gr\_e^4} \tag{43}$$

for a given shear rate *G* and ellipsoidal aspect ratio *r*e [50]. Bending of fibers starts with an increased aspect ratio, higher shear rates and reduced bending stiffness [51]. The dimensionless stiffness *S* summarises these effects in one parameter.

First, a stiff fiber with *S* = 100 is considered. Such a fiber spins in a rigid manner without significant bending and can be compared to the solution of Jeffery's equation

$$\frac{D\mathbf{p}}{Dt} = \omega \mathbf{p} + \zeta \left( \mathbf{D} \mathbf{p} - \left( \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \right) \mathbf{D} \right) \tag{44}$$

with vorticity tensor *ω* and symmetric strain rate tensor **D**. The shape factor *ξ* is an alternative measure for the (equivalent) ellipsoidal aspect ratio *r*e and is defined as

$$f\_{\mathbf{e}}(r\_{\mathbf{e}}) = \frac{r\_{\mathbf{e}}^2 - 1}{r\_{\mathbf{e}}^2 + 1}. \tag{45}$$

The solution to Equation (44) is periodic with

$$T = \frac{2\pi}{G} \left( r\_{\text{e}} + \frac{1}{r\_{\text{e}}} \right) \tag{46}$$

being the time for a full rotation of a fiber. When solving Jeffery's equation for a cylinder with geometric aspect ratio *r*p instead of an ellipsoid, Jeffery's equation can still be applied, but an equivalent aspect ratio has to be used. Such an equivalent aspect ratio *r*e was derived by Cox [52] based on slender-body theory and Zhang et al. [53] utilizing Finite Element Analysis. The latter is applicable for small aspect ratios and therefore Zhang's cubic fit

$$r\_{\rm e}(r\_{\rm p}) = 0.000035r\_{\rm p}^3 - 0.00467r\_{\rm p}^2 + 0.764r\_{\rm p} + 0.404\tag{47}$$

is used here.

Figure 5 depicts the orientation *φ* of a single fiber in shear flow for fiber length *L*f = 5Δ*x* and *L*f = 11 Δ*x* (e.g., the fiber is represented by 5 or 11 particles in a chain). The rotation period obtained with the present SPH approach is sligthly faster than the solution to Jeffery's equation. One possible reason for the difference is the finite Reynolds number and finite simulation domain with periodic boundaries, whereas Jeffery used an infinite domain with strictly no effect of inertia. Another reason might be the coarse resolution with only one layer of particles for the fiber. Due to the averaging nature of SPH, the exact flow field close to the suspended particles cannot be resolved exactly. The presented approach solves the entire fluid field, but the low resolution and smoothing makes it less accurate than e.g., Stokesian Dynamics simulations. However, the simplicity and computational efficiency make it attractive for engineering applications, such as the parameter fitting presented later.

**Figure 5.** Orientation angle *φ* for fiber in a shear flow. The fiber length *L*f is expressed in multiples of the particle spacing Δ*x*. The dashed gray line represents the solution to Jeffery's equation with Zhang's [53] fit. The solid black line represents the solution obtained with this SPH implementation. If the viscous surface traction term (Equation (24)) is neglected, the fiber stops rotating after one half rotation, as shown with the black dotted line.

Bending modes were shown in Yamamoto and Matsuka's [13] numerical results and the corresponding SPH simulation in Table 1 agree well with their observations. Differences can be attributed to the fact that Yamamoto and Matsuka used an inextensible fiber, while the fibers simulated in this work experience stretching, because the Young's modulus is taken into account for tensile stiffness as well.


**Table 1.** Bending modes of a fiber with length *L*f = 11Δ*x* for varied dimensionless stiffness. One example with *S* = 10 and critical fiber bending angle *θ*c = *π*4 is shown to demonstrate fiber fracture.

This section illustrated that the implementation based on SPH can reproduce the rotation periods quantitatively. Furthermore, numerically obtained bending modes can be described qualitatively. In addition, a fully coupled solution for the fluid field is computed. It can be noted that the fluid field obtained for these single fiber setups does not significantly differ from the ideal field. This is

reasonable, since a single flexible fiber does not offer much resistance to a highly viscous flow. However, a significant difference can be expected as soon as multiple fibers interact in a concentrated suspension. In that scenario, the presence of suspended fibers and their interactions are expected to raise the macroscopically observed effective viscosity and fiber interactions affect the fiber orientation evolution.

#### *3.2. Parameter Identification for the Orientation Evolution in a Non-Dilute Short Fiber Suspensions*

A fully resolved computation of all fibers is often not feasible for full components made from composite material. It can be sufficient to give a reasonable description of the fiber orientation function in terms of a fiber orientation tensor, if these are accurate and scale-separation applies. A common two-parameter model is the RSC model [5] with fiber interaction coefficient *C*I and a phenomenological factor *κ* that models a delay to compensate an over-prediction in the change of orientation observed in the classical Folgar-Tucker model. In tensor notation, the RSC model reads

$$\frac{D\mathbf{A}}{Dt} = \mathbf{\omega}\mathbf{A} - \mathbf{A}\boldsymbol{\omega} + \boldsymbol{\xi}\left(\mathbf{D}\mathbf{A} + \mathbf{A}\mathbf{D} - 2\left(\mathbf{A} + (1-\kappa)(\mathbb{L}-\mathbb{M}\mathbf{A})\right)\mathbf{D}\right) + 2\kappa\mathbf{C}\_{\mathrm{I}}\mathbf{G}\left(\mathbf{1} - 3\mathbf{A}\right) \tag{48}$$

with L = ∑<sup>3</sup>*i*=<sup>1</sup> *λi***e***i* ⊗ **e***i* ⊗ **e***i* ⊗ **e***i* and M = ∑<sup>3</sup>*i*=<sup>1</sup> **e***i* ⊗ **e***i* ⊗ **e***i* ⊗ **e***i* using the eigenvalues *λi* and eigenvectors **e***i* of the second order fiber orientation tensor **A**. Setting *κ* = 1 would reduce this model to the Folgar-Tucker model and setting *C*I = 0 reduces it to Jeffery's Equation (44). The choice of feasible parameters *C*I ∈ [0, 0.1] and *κ* ∈ [0, 1.0] remains. Hence, this section computes the orientation evolution in terms of the second order fiber orientation tensor for a small set of fibers in a 3D shear flow and compares the solution obtained with SPH to macroscopic models.

The investigated domain is a cube with edge length *L* = 15Δ*x* and Lees-Edwards boundary conditions [54,55] are employed to induce a shear rate *G* on the periodic fluid domain. Essentially, these boundary conditions shift dummy particles and particles leaving the domain in *<sup>x</sup>*1-direction according to

$$\mathbf{x}\_2' = (\mathbf{x}\_2 \pm GLt) \bmod L \tag{49}$$

with the sign depending on the direction of the shift and the modulo operator mod. Additionally, the velocity is adjusted

$$v\_2' = (v\_2 \pm GL) \tag{50}$$

for the shift. Conventional periodic boundary conditions apply in *x*2- and *<sup>x</sup>*3-direction, as depicted in Figure 6. All fibers are initially oriented in *<sup>x</sup>*1-direction and randomly positioned in the cube. This unidirectional arrangemen<sup>t</sup> is chosen because it can be easily achieved without a micro-structure generator. Instead of analyzing larger representative volumes [26], this work performs multiple realizations of the random process to create a statistical representative behavior. The advantage of this approach is that scatter and standard deviation between random realizations on the micro scale can be observed. The ensemble average of a property ·*tn* is defined as the mean across multiple realizations at the same time step *tn*.

To obtain optimal parameters at different volume fractions, a least squares fit

$$\underset{\mathcal{C}\_{i},\mathbf{x}}{\text{minimize}} \quad \sum\_{n} \left\| \mathbf{A}(t\_{n}) - \langle \tilde{\mathbf{A}} \rangle\_{t\_{n}} \right\|^{2} \tag{51}$$

is applied to minimize the squared difference of **A** obtained by Equation (48) and the ensemble average of the discrete second order fiber orientation tensor **A ˜** of the SPH analysis. This tensor can be computed as

$$\mathbf{\tilde{A}} = \frac{1}{N} \sum\_{i=1}^{N} \mathbf{p}\_i \otimes \mathbf{p}\_i \tag{52}$$

with each fiber's orientation **p***i*. A flexible fiber has different tangential orientations and the tensor's definition becomes ambiguous then. Consequently, the following examples use fibers with a high stiffness *S* = 100 to ensure enough rigidity for an unambiguous interpretation of **A**. However, the method is in no way limited to rigid fibers.

**Figure 6.** Setup for the 3D shear with fibers in a cube of edge length *L* = 15Δ*x*. Lees-Edwards boundary conditions [54] are employed to induce a shear rate *G*. Therefore dummy particles (light gray) and particles leaving the domain in *<sup>x</sup>*1-direction are shifted periodically in *x*2 during each domain update. Conventional periodic boundary conditions apply to all other sides of the cube. The initial state is generated from fibers aligned in *<sup>x</sup>*1-direction at unique random positions in the entire volume. (Some fibers appear longer in this figure due to other overlapping fibers behind it.)

Figure 7 shows the non-trivial components of the ensemble average **A˜** *tn* computed from simulations with five different initial random realizations as a solid green line. The standard deviation at each time step is depicted as light filled area in the background. The gray solid line represents orientation tensor components **A** computed with optimal parameters according to the RSC model given in Equation (48). The simulation time is 10*T*, which is the time of 10 full rotations of a single rigid fiber. The corresponding strain is approximately 150.

The simple parameter fitting approach proposed in (51) works well and generally shows a good agreemen<sup>t</sup> of macroscopic models with the SPH micro-model. All results show a decreasing orientation amplitude and a trend towards a stationary state with a significant non-zero component in *<sup>x</sup>*3-direction. This trend arises from a combination of fiber contacts and long range pertubations of the flow field that push fibers out of their original trajectories. Figure 8 shows, how fibers are oriented after 30 strains in the case of 10 % fiber volume fraction. Several fibers have left the sheared *<sup>x</sup>*1*<sup>x</sup>*2-plane due to interactions at this point. Eventually, fiber interactions and shear-induced reorientation balance each other and lead to a stationary orientation state. The stationary orientation state is reached faster for increased fiber volume fractions.

**Figure 7.** Ensemble average of fiber orientation tensor components for fibers with aspect ratio *r*p = 4.43 (*L*f = 5) in a 3D shear flow and its comparison to the reduced strain model with optimal parameters. Each simulation result was obtained from five independently sampled initial configurations. The standard deviation is indicated by a light green filled area for each volume fraction.

0.011405, *κ* = 0.74

(**c**) *φ* = 30%, 202 fibers, *C*I =

The deviations between different initial configurations decrease for increasing fiber volume fraction, as the sample size increases. The deviation between individual realizations at 1% volume fraction is large and it might not be appropriate to describe such a system with a macroscopic fiber orientation model. This highlights that a macroscopic description requires a sufficient scale separation and a sufficient number of fibers to provide reliable results. The proposed SPH simulation can be used to quickly evaluate different configurations and may be used as a tool to not only determine parameters, but also quantify deviations from macromodels, if the underlying conditions of such models are in question for a specific application.

The obtained parameters of the interaction coefficient are compared to Folgar and Tucker's original work [4] and the values obtained by Phan-Thien et al. [47] in Figure 9. The results from SPH simulations show a good agreemen<sup>t</sup> with literature data and support the use of this approach for parameter identification.

**Figure 8.** Snapshot of fibers at 10% fiber volume fraction after 30 strains. The colors are introduced to distinguish fiber particles and represent the particle ID.

**Figure 9.** Comparison of *C*I values to the original work of Folgar and Tucker [4] and a fit based on simulation results by Phan-Tien et al. [47]. The values obtained in the presented work are reasonably close to these literature results.
