**4. Conclusions**

Fiber suspensions are treated as flexible bead chains of SPH particles surrounded by other particles representing the fluid domain. The bead chains are connected by elastic tension- and bending forces and interact with the fluid particles in a two-way coupled manner. A novel viscous surface traction term is introduced to compensate the missing fiber thickness that is introduced by the line representation of a fiber. In addition, contact forces are introduced to model fiber interactions in a non-dilute suspension.

The fiber orientation evolution of a single stiff fiber shows good agreemen<sup>t</sup> with the rotation periods based on Jeffery's equation thanks to the introduction of a new surface traction term. Bending modes of single fibers are consistent with results reported in literature.

A periodic domain with Lees-Edwards boundary conditions and suspended fibers is subjected to shear. The investigation of suspensions with different volume fractions of fibers can be used to directly compute the fiber orientation tensor. If several of these computations are evaluated in a statistical sense, the ensemble average can be used to fit optimal parameters of fiber orientation models. For relatively stiff fibers of length *L*f = 5Δ*x*, good parameters of the RSC fiber orientation model are found based on the SPH simulation.

In future, the authors plan to extend the investigation to arbitrary initial configurations, flexible and breaking fibers as well as other flow types beyond shear flows. Further, the assessment of standard deviations may enable modeling of uncertainties related to the fiber orientation models.

**Author Contributions:** N.M.: writing–original draft preparation, methodology, software, visualization; O.S. and M.H.: conceptualization, writing–review and editing; A.N.H.: supervision, conceptualization, writing–review and editing; F.H.: supervision, funding acquisition; L.K.: supervision, methodology, writing–review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research documented in this manuscript has been funded by the German Research Foundation (DFG) within the International Research Training Group "Integrated engineering of continuous-discontinuous long fiber-reinforced polymer structures" (GRK 2078). The support by the German Research Foundation (DFG) is gratefully acknowledged.

**Acknowledgments:** We acknowledge support by the KIT-Publication Fund of the Karlsruhe Institute of Technology.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Evaluation of the Surface Traction Integral**

An arbitrary normal to the fiber direction **p***i* named **n**0 must be generated first. This can be achieved by forming the cross product of **p***i* with an arbitrary other vector that is not parallel to **p***i*. Chosing the arbitrary vector [1, 0, 0] , the parametrized normal in (20) becomes

$$\mathbf{n}\_i(\boldsymbol{\psi}) = \begin{pmatrix} -\sin(\boldsymbol{\psi})\sqrt{p\_i^2 + p\_i^2} \\ \frac{p\_{i1}\sin(\boldsymbol{\psi})p\_{i2} + \cos(\boldsymbol{\psi})p\_{i3}}{\sqrt{p\_i^2 + p\_{i3}^2}} \\ -\frac{-p\_{i1}\sin(\boldsymbol{\psi})p\_{i3} + \cos(\boldsymbol{\psi})p\_{i2}}{\sqrt{p\_i^2 + p\_{i3}^2}} \end{pmatrix}. \tag{A1}$$

The integration in (23) leads to

**M**<sup>v</sup> *i* = 3Δ*x d*2 4 *ηπ pi* 2 2 + *pi* 2 3 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ (*pi*1 *pi* 2 2 *pi*3 + *pi*1 *pi* 3 3) *∂<sup>v</sup>*2 *∂<sup>x</sup>*1 + (−*pi* 2 1 *pi*2 *pi*3 + *pi*2 *pi*3) *∂<sup>v</sup>*2 *∂<sup>x</sup>*2 + (−*pi* 2 1 *pi* 2 3 − *pi* 2 2) *∂<sup>v</sup>*2 *∂<sup>x</sup>*3 + (−*pi*1 *pi* 3 2 − *pi*1 *pi*2 *pi* 2 3) *∂<sup>v</sup>*3 *∂<sup>x</sup>*1 + (*pi* 2 1 *pi* 2 2 + *pi* 2 3) *∂<sup>v</sup>*3 *∂<sup>x</sup>*2 + (*pi* 2 1 *pi*2 *pi*3 − *pi*2 *pi*3) *∂<sup>v</sup>*3 *∂<sup>x</sup>*3 (−*pi*1 *pi* 2 2 *pi*3 − *pi*1 *pi* 3 3) *∂<sup>v</sup>*1 *∂<sup>x</sup>*1 + (*pi* 2 1 *pi*2 *pi*3 − *pi*2 *pi*3) *∂<sup>v</sup>*1 *∂<sup>x</sup>*2 + (*pi* 2 1 *pi* 2 3 + *pi* 2 2) *∂<sup>v</sup>*1 *∂<sup>x</sup>*3 + (−*pi* 4 2 − 2*pi* 2 2 *pi* 2 3 − *pi* 4 3) *∂<sup>v</sup>*3 *∂<sup>x</sup>*1 + (*pi*1 *pi* 3 2 + *pi*1 *pi*2 *pi* 2 3) *∂<sup>v</sup>*3 *∂<sup>x</sup>*2 + (*pi*1 *pi* 2 2 *pi*3 + *pi*1 *pi* 3 3) *∂<sup>v</sup>*3 *∂<sup>x</sup>*3 (*pi*1 *pi* 3 2 + *pi*1 *pi*2 *pi* 2 3) *∂<sup>v</sup>*1 *∂<sup>x</sup>*1 + (−*pi* 2 1 *pi* 2 2 − *pi* 2 3) *∂<sup>v</sup>*1 *∂<sup>x</sup>*2 + (−*pi* 2 1 *pi*2 *pi*3 + *pi*2 *pi*3) *∂<sup>v</sup>*1 *∂<sup>x</sup>*3 + (*pi* 4 2 + 2*pi* 2 2 *pi* 2 3 + *pi* 4 3) *∂<sup>v</sup>*2 *∂<sup>x</sup>*1 + (−*pi*1 *pi* 3 2 − *pi*1 *pi*2 *pi* 2 3) *∂<sup>v</sup>*2 *∂<sup>x</sup>*2 + (−*pi*1 *pi* 2 2 *pi*3 − *pi*1 *pi* 3 3) *∂<sup>v</sup>*2 *∂<sup>x</sup>*3 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (A2)

with the velocity gradient from Equation (25). For the special case that **p***i* is equivalent to [1, 0, 0] , the moment is computed in the same way with a different arbitrary initial direction, e.g., [0, 1, 0] .
