**2. Theory**

SPH was developed independently by Lucy [39] and Gingold & Monaghan [40] for astrophysical problems in 1977. Since then, multiple formulations were developed and were applied to various fields. The core idea of SPH is the conversion of a partial differential equation (PDE) for a continuum into a set of ordinary differential equations (ODE) for multiple Lagrangian particles. The ODEs are then integrated for each particle to determine their properties (such as mass density or velocity) and position for the next time step. Since particles carry the mass and are moved according to the velocity field, this method conserves mass exactly and treats pure advection correctly. Any arbitrary property *α*(*x*) in an n-dimensional domain Ω ⊂ R*n* at a position **x**, **x** ∈ Ω can be described using an integral of the form

$$\alpha(\mathbf{x}) = \int\_{\Omega} \alpha(\mathbf{x}') \delta(\mathbf{x} - \mathbf{x}') d\mathbf{x}' \tag{3}$$

employing the Dirac distribution *<sup>δ</sup>*(**x**). The two fundamental approximations of SPH are the description of the continuous integral over the domain using a sum over interpolation points, that can be interpreted as particles, and the replacement of the Dirac distribution with a smooth kernel *<sup>W</sup>*(**x** − **x** , *h*). Here, *h* describes the smoothing length that controls the kernel size and is usually chosen similar to the average initial distance between particles. A very common smooth kernel is the cubic spline kernel defined as

$$\mathcal{W} = \beta\_{\text{tr}} \begin{cases} (2-q)^3 - 4(1-q)^3 & 0 \le q < 1 \\ (2-q)^3 & 1 \le q < 2 \\ 0 & 2 \le q \end{cases} \tag{4}$$

with *q* = **<sup>x</sup>**−**<sup>x</sup>** *h* and a dimension-dependent normalization factor *β<sup>n</sup>*. Employing the kernel *<sup>W</sup>*(**<sup>x</sup>***i* − **<sup>x</sup>***j*, *h*) = *Wij*, Equation (3) can be approximated as

$$\mathfrak{a}(\mathbf{x}) \approx \mathfrak{h}(\mathbf{x}\_i) = \mathfrak{k}\_i = \sum\_{j \in \Omega} \mathfrak{a}\_j \mathcal{W}\_{ij} \frac{m\_j}{\rho\_j} \tag{5}$$

denoting a particle's property *α* at position *xj* as *αj* and by expressing the associated volume as ratio between mass *mj* and density *ρj* of the particle. The key is the usage of an analytically differentiable kernel such that

$$\frac{\partial \mathbb{K}\_i}{\partial \mathbf{x}} = \sum\_{j \in \Omega} \alpha\_j \frac{m\_j}{\rho\_j} \nabla\_i \mathcal{W}\_{ij} \tag{6}$$

can be computed analytically with the kernel gradient *∂W∂***x** (**x** − **<sup>x</sup>***j*, *h*)|*i* = <sup>∇</sup>*iWij*, hence turning the PDE to an ODE. The domain Ω = Ωm ∪ Ωs is the set of all particles in the model with subsets for the fluid particles Ωm and solid particles Ωs = Ωw ∪ Ωf consisting of wall particles Ωw and fiber particles Ωf. There is a wide range of variants and actual implementations for this concept and Monaghan gives a comprehensive review in one of his later publications [41].

The modelling of the fluid phase and its interaction with solid particles follows the work of Adami et al. [42,43]. The basic fiber model is adapted from Yang et al. [38] and extended to account for viscous surface traction as well as fiber interactions.

#### *2.1. Fluid Model*

In this work, a Transport Velocity Formulation [43] is used to model the suspending fluid, since it is fairly robust against stability issues such as the tensile instability [44]. Adami et al. [43] separate momentum velocity **v** and advection velocity ˜**v** leading to a Navier-Stokes equation of the form

$$\frac{d\mathbf{\hat{d}}(\rho\mathbf{v})}{dt} = -\nabla p + \eta\nabla^2 \mathbf{v} + \rho\mathbf{g} + \nabla \cdot (\rho \mathbf{v} \otimes (\bar{\mathbf{v}} - \mathbf{v})) \tag{7}$$

with density *ρ*, pressure *p*, viscosity *η* and body accelerations **g**. The last term is a momentum convection that compensates the deviation between advection velocity and momentum velocity. The difference is based on a virtual background pressure *p*b that effectively suppresses tensile instability, but does not influence the actual momentum. Hence, the term *σ*<sup>A</sup> = *ρ***v** ⊗ (**v**˜ − **v**) is called an artificial stress. The SPH-discretized version of (7) used in this work models the acceleration of particle *i* in the fluid domain by

$$\begin{split} \bar{\mathbf{a}}\_{i} &= \frac{1}{m\_{i}} \sum\_{j \in \Omega} \left( V\_{i}^{2} + V\_{j}^{2} \right) \left( -\frac{\rho\_{j}\eta\_{j} + \rho\_{i}\eta\_{j}}{\rho\_{i} + \rho\_{j}} \nabla\_{i} \mathbb{W}\_{ij} \right) \\ &+ \frac{1}{m\_{i}} \sum\_{j \in \Omega\_{o}} \left( V\_{i}^{2} + V\_{j}^{2} \right) \left( \frac{2\eta\_{i}\eta\_{j}}{\eta\_{i} + \eta\_{j}} \frac{\nabla\_{i} \mathbb{W}\_{ij} \cdot (\mathbf{x}\_{i} - \mathbf{x}\_{j})}{\|\mathbf{x}\_{i} - \mathbf{x}\_{j}\| + \varepsilon} (\mathbf{v}\_{i} - \boldsymbol{\Psi}\_{j}) \right) \\ &+ \frac{1}{m\_{i}} \sum\_{j \in \Omega\_{f}} \left( V\_{i}^{2} + V\_{j}^{2} \right) \left( \frac{2\eta\_{i}\eta\_{j}}{\eta\_{i} + \eta\_{j}} \frac{\nabla\_{i} \mathbb{W}\_{ij} \cdot (\mathbf{x}\_{i} - \mathbf{x}\_{j})}{\|\mathbf{x}\_{i} - \mathbf{x}\_{j}\| + \varepsilon} (\mathbf{v}\_{i} - \mathbf{v}\_{j}) + \frac{1}{2} (\mathbf{\sigma}\_{i}^{\Lambda} + \mathbf{\sigma}\_{j}^{\Lambda}) \nabla\_{i} \mathbb{W}\_{ij} \right), \quad i \in \Omega\_{\mathbf{m}} \tag{8} \end{split}$$

with the volume attributed to each particle *Vi* = *mi*/*ρi* and the velocity of an adjacent solid **v**ˆ*j*, which is explained in the next section. Indices are used to refer to the particle based density *ρi*, pressure *pi*, viscosity *ηi*, mass *mi*, position **x***i* and velocity **<sup>v</sup>***i*. The parameter *ε* is a small value to avoid singularities in the formulation. The last term is a momentum convection that compensates the deviation between advection velocity and momentum velocity.

Finally, the discrete conservation of mass is written as

$$\rho\_{\dot{i}} = m\_{\dot{i}} \sum\_{j \in \Omega} \mathcal{W}\_{\dot{i}j}, \quad \dot{i} \in \Omega\_{\text{m}} \tag{9}$$

and an equation of state is used to relate mass density to pressure. The equation of state has the form

$$p\_i = p\_0 \left( \left( \frac{\rho\_i}{\rho\_0} \right)^{\gamma} - 1 \right), \quad i \in \Omega\_{\text{m}} \tag{10}$$

where *ρ*0 describes the nominal density and *p*0 is a reference pressure chosen large enough to keep the density variation small. The latter is achieved by setting *p*0 = *<sup>ρ</sup>*0*c*<sup>2</sup> *γ* with the speed of sound *c* set to the tenfold of the maximum speed in the flow, thus limiting density variation to approximately 1%. Following Adami et al. [43], the parameter *γ* is set as *γ* = 1.

#### *2.2. Interaction between Fluid and Solid Particles*

The summation of the first term in Equation (8) includes fiber particles and solid wall particles. Adami et al. [42] determined the pressure of a solid particle from a force balance along the centerline of a solid-fluid particle pair as

$$p\_i = \frac{\sum\_{j \in \Omega\_{\text{m}}} \rho\_j \mathcal{W}\_{ij} + (\mathbf{g} - \mathbf{a}\_i) \cdot \sum\_{j \in \Omega\_{\text{m}}} \rho\_j (\mathbf{x}\_i - \mathbf{x}\_j) \mathcal{W}\_{ij}}{\sum\_{j \in \Omega\_{\text{m}}} \mathcal{W}\_{ij}}, \quad i \in \Omega\_{\text{s}} \tag{11}$$

with a prescribed acceleration of the solid **a***i*. The corresponding density may be computed by inverting (10). To ensure a no-slip condition, the velocity of solid particles is modified before insertion in (8) to

$$\mathfrak{v}\_{i} = 2\upsilon\_{i} - \frac{\sum\_{j \in \Omega\_{\mathfrak{m}} \upsilon\_{j} W\_{ij}}}{\sum\_{j \in \Omega\_{\mathfrak{m}} W\_{ij}}}, \quad i \in \Omega\_{\mathfrak{s}} \tag{12}$$

where the fluid velocity field is extrapolated and subtracted to ensure zero velocity at the interface between solid particles and fluid particles [42].

In this work, wall particles are represented by three particle layers to ensure full kernel support and move with a constant prescribed velocity

$$\mathbf{v}\_{i} = \text{const}, \quad \mathbf{a}\_{i} = 0, \quad i \in \Omega\_{\text{W}}.\tag{13}$$

Fibers are represented by a single chain of particles that experience hydrodynamic forces from the fluid, elastic forces from neighbors within the fiber and contact forces from other fibers. The first contribution to acceleration is a hydrodynamic interaction

$$\mathbf{a}\_{l}^{\text{hydro}} = \frac{1}{m\_{l}} \sum\_{j \in \Omega\_{f}} \left( V\_{l}^{2} + V\_{f}^{2} \right) \left( -\frac{\rho\_{f} p\_{j} + \rho\_{l} p\_{l}}{\rho\_{l} + \rho\_{f}} \nabla\_{l} \mathcal{W}\_{lj} + \frac{2 \eta\_{l} \eta\_{j}}{\eta\_{l} + \eta\_{j}} \frac{\nabla\_{l} \mathcal{W}\_{lj} \cdot (\mathbf{x}\_{l} - \mathbf{x}\_{j})}{\|\mathbf{x}\_{l} - \mathbf{x}\_{j}\| + \varepsilon} (\boldsymbol{\Psi}\_{l} - \mathbf{v}\_{f}) \right), \quad i \in \Omega\_{l} \tag{14}$$

that balances the momentum together with (8). Modelling the fiber as a single layered chain of SPH particles has also been applied by other authors [38,45], but it has some implications, which are discussed further in Section 2.4. One implication is that the particle spacing Δ*x* is directly related to the fiber diameter by

$$d = \frac{2\Delta x}{\sqrt{\pi}}\tag{15}$$

to ensure that fiber particles and fluid particles describe equal volumes in the beginning of a simulation.

#### *2.3. Basic Model for Flexible Fibers*

Besides the description of the fluid phase with SPH, a suitable model for the elastic fiber is needed. Thus, the straight-forward linear elastic bead chain formulation of Yang et al. [38] for tensile forces

$$\mathbf{F}\_{ij}^{t} = EA \left( \frac{||\ \mathbf{x}\_{ij} ||}{\mathbf{x}\_{ij}^{0}} - 1 \right) \frac{\mathbf{x}\_{ij}}{||\ \mathbf{x}\_{ij} ||} \tag{16}$$

and bending moment

$$\mathbf{M}\_{i}^{\mathbf{b}} = \frac{EI}{2} \frac{\mathbf{\theta}\_{i} - \mathbf{\theta}\_{i}^{0}}{\|\mathbf{x}\_{i(i-1)}\| + \|\mathbf{x}\_{i(i+1)}\|} \tag{17}$$

can be used as a basis for further model development. Here, *E* describes Young's modulus, while *A* and *I* are the fiber's cross sectional area and second moment of area respectively. The vector between two neighboring particles with indices *i* and *j* is **<sup>x</sup>***ij* and the enclosed angle is denoted as *θi*. Here, the vector notation of *θi* indicates that its direction resembles the axis of rotation. The undeformed reference configuration is denoted with a superscript (·)<sup>0</sup> in both, Equations (16) and (17). The bending moment can be converted to pairs of forces

$$\mathbf{F}\_{ij}^{\mathbf{b}} = \frac{1}{2} \frac{\mathbf{M}\_i^{\mathbf{b}} \times \mathbf{x}\_{ij}}{||\!\!\!/ \mathbf{x}\_{ij}||^2} \tag{18}$$

that act on the particle and its neighbors. Figure 1 illustrates these forces and it can been seen, that generally **F**b*ij* = **F**b*ji*. It is assumed that torsional torque of the fiber is of minor importance to the orientation evolution investigated in this work. Finally, the acceleration on a particle *i* due to elastic forces can be summarized as

$$\mathbf{a}\_{i}^{\text{elastic}} = \frac{1}{m\_{i}} \left( \mathbf{F}\_{i(i+1)}^{\text{t}} + \mathbf{F}\_{i(i-1)}^{\text{t}} \right. \tag{19}$$

$$+ \mathbf{F}\_{i(i-1)}^{\text{b}} + \mathbf{F}\_{i(i+1)}^{\text{b}} - \mathbf{F}\_{(i-1)i}^{\text{b}} - \mathbf{F}\_{(i+1)i}^{\text{b}} \right).$$

**Figure 1.** Force pairs representing bending moment on segments between fiber particles.

If the angle between two adjacent particle pairs exceeds a certain critical value *θ*c, the neighborhood relation between these particles may be revoked permanently to model fracture of the fiber. Such a criterion is motivated by brittle fiber fracture, as it is typical for glass fibers.

#### *2.4. Viscous Traction at Fiber Surface*

A fiber modeled as a particle chain cannot capture a variation of a property in thickness direction of the fiber, as it stores properties at the center line only. Hence, a fiber placed horizontally in a shear flow with periodic boundary conditions, as depicted in Figure 2, would not experience any moment caused by friction forces on its surface. In order to model a physical thickness of the fiber, this work uses an analytical derivation to apply appropriate moment from surface friction to the fiber segment. Figure 3 illustrates a cylindrical fiber segmen<sup>t</sup> of length Δ*L* and diameter *d* with the orientation of its centerline **p**. One exemplary surface normal **n** is shown with its parametrization angle *ψ*.

The fiber direction **p** and any arbitrary surface normal **n**<sup>0</sup> are perpendicular unit-vectors. Any other surface normal can be constructed from this arbitrary normal by rotating it around **p**. The surface normal can be parameterized using **n**<sup>0</sup> and *ψ* employing a rotational tensor **R** around axis **p** [46] . The parameterized normal becomes

$$\mathbf{n}(\psi) = \mathbf{R}(\psi)\mathbf{n}\_0. \tag{20}$$

Using this normal, the viscous traction on the surface can be expressed as

$$\mathbf{t}(\boldsymbol{\psi}) = \left(-p\mathbf{I} + \eta \nabla \mathbf{v}\right) \mathbf{n}(\boldsymbol{\psi}) = -p\mathbf{n}(\boldsymbol{\psi}) + \eta \nabla \mathbf{v} \mathbf{n}(\boldsymbol{\psi})\tag{21}$$

if Newtonian viscosity is assumed. This term represents the force acting on each infinitesimal area of the fiber surface. The resulting moment can be computed by integrating **t** with its corresponding leverage *d*2**<sup>n</sup>**(*ψ*) as

$$\mathbf{M}^{\rm v} = \Delta L \int\_0^{2\pi} \frac{d}{2} \mathbf{n}(\boldsymbol{\psi}) \times \mathbf{t}(\boldsymbol{\psi}) \frac{d}{2} \mathbf{d}\boldsymbol{\psi} \tag{22}$$

with a constant fiber diameter *d*. The term *d*2d*ψ* represents an infinitesimal circumferential line segmen<sup>t</sup> on the cylinder surface. As the cross product of a vector with itself vanishes, this can be simplified to

$$\mathbf{M}^{\rm v} = \Delta L \int\_0^{2\pi} \frac{d^2}{4} \mathbf{n}(\psi) \times \eta \nabla \mathbf{v} \mathbf{n}(\psi) \mathrm{d}\psi. \tag{23}$$

The diameter is finite and thus provides some leverage for the traction to generate a moment. The discrete evaluation of (23) is explained in Appendix A. The equation for angular momentum is then multiplied with the distance to its neighbors as a cross product leading to the acceleration of individual particles due to viscous friction

$$\mathbf{a}\_{i}^{\text{traction}} = \frac{1}{2} \frac{\mathbf{M}\_{i}^{\text{V}} \times \mathbf{x}\_{i\bar{j}}}{\bar{f}}, \quad j \in [i-1, i+1] \tag{24}$$

Here, *J* denotes the moment of inertia for a cylindrical body around its first principle axis of inertia and the factor 12 is chosen to represent the moment by two equal forces at both neighboring particles. This acceleration is then applied to the two neighboring particles and implies hereby a rotational acceleration of a fiber segment, consisting of three particles Δ*L* = 3Δ*x*. The central particle is used to evaluate the velocity gradient of Equation (21) as

$$
\nabla \mathbf{v}\_i = -\frac{1}{\rho\_i} \sum\_j m\_j (\mathbf{v}\_i - \mathbf{v}\_j) \nabla\_i \mathcal{W}\_{ij}. \tag{25}
$$

It is assumed that the velocity gradient is approximately constant within each segment. In theory, the velocity gradient could be determined at each point of the cylindrical surface from kernel functions, but this comes at much higher computational costs and the difference in the resulting moment is expected to be small.

**Figure 2.** A fiber is placed horizontally with *φ* = *π*/2 in a shear flow. The top and bottom walls have a prescribed velocity and a no-slip condition, the boundary conditions in *x*1 direction are periodical.

**Figure 3.** Cylindrical fiber segmen<sup>t</sup> with length Δ*L*, orientation **p** and an arbitrary surface normal **n**.

#### *2.5. Fiber Interactions*

If suspended objects come in close contact (10–50% of the radius [26,47]), lubrication forces oppose the relative velocity between fibers. Sundararajakumar and Koch [24] showed that lubrication forces alone do not necessarily prevent penetration at higher fiber volume fractions and added contact forces. The pressure gradient computed from SPH is generally not sufficient to counteract the accumulated forces on the fiber bead chain and does not necessarily prevent penetration of fibers.

It is too simple to apply contact forces directly bead to bead, because this would lead to entangled fibers that interlock at two bead gaps. For the determination of the contact properties at each potential contact pair (*i*, *j*) within the kernel radius, different cases need to be considered:

**Surface-Surface:** If both interacting particles are located at the center of the fiber (e.g., they have two neighbor particles each, compare Figure 4 at *t*0), the normal direction of contact pair (*i*, *j*) can be computed using the cross product

$$\mathbf{n}\_{ij} = \begin{bmatrix} \mathbf{p}\_i \times \mathbf{p}\_j \end{bmatrix} \tag{26}$$

of the involved fiber direction vectors **p***i* and **p***j*. The operator [[·]] = (·) · is used to conveniently denote the normalization of a vector. Solving the small linear system of equations

$$[\mathbf{p}\_{i\prime}, \mathbf{n}\_{i\prime\prime} - \mathbf{p}\_{j}][\mathcal{P}\_{i\prime}, D\_{i\prime}, \mathcal{P}\_{j\prime}]^{\top} = \mathbf{x}\_{i\prime} \tag{27}$$

with its adjugate matrix leads to the solution for the distance between the fibers *Dij* and the projections to source and destination vectors <sup>P</sup>*ij* and <sup>P</sup>*ji*, respectively.

**Figure 4.** Snapshots of two fibers in contact. At contact initiation *t*0, **p***i* and **p***j* denote unit vectors for the fiber directions and **<sup>x</sup>***ij* is the vector between particles of one active contact pair (*i*, *j*). The contact forces are indicated by arrows which scale with contact force magnitude and rotate in the subsequent time steps (*<sup>t</sup>*1, *t*2, *t*3).

**Surface-End and End-Surface:** If a particle of a fiber end interacts with a central particle of another fiber, the vector between these two particles **<sup>x</sup>***ij* can be used to obtain the normal direction by projection. It is assumed that **p***i* describes a unit vector in fiber direction at one fiber particle at position A. Let **<sup>x</sup>***ij* be the vector from another fibers' end particle B to the point A. The closest point to B on a line with direction **p***i* is denoted as C and can be used to define the normal direction as

$$\mathbf{n}\_{ij} = \|\overrightarrow{\mathcal{BC}}\| = \|\overrightarrow{\mathcal{AC}} - \overrightarrow{\mathcal{AC}}\|. \tag{28}$$

Using the definitions above and the fact that C is the projection of B to the line with direction **p***i*, Equation (28) can be rewritten as

$$\mathbf{n}\_{ij} = \left[\mathcal{P}\_{ij}\mathbf{p}\_i - (-\mathbf{x}\_{ij})\right] = \left[\mathbf{x}\_{ij} + \mathcal{P}\_{ij}\mathbf{p}\_i\right].\tag{29}$$

The projection on the destination fiber is given as <sup>P</sup>*ij* = −**p***i* · **<sup>x</sup>***ij* and the contact distance is computed as *Dij* = **<sup>x</sup>***ij* + <sup>P</sup>*ij***p***i* .

**End-End:** The simplest case is the interaction of two fiber ends. Here, the vector between those two particles can be simply determined by

$$\mathbf{n}\_{ij} = \begin{bmatrix} \mathbf{x}\_{ij} \end{bmatrix} \tag{30}$$

with the corresponding distance *Dij* = **<sup>x</sup>***ij* .

A penalty approach is proposed to prevent fiber penetration, if the distance between fibers falls below the fiber diameter.

The penalty force is formulated as a Hertzian contact force between two cylinders [48]

$$F\_{ij}^c = \begin{cases} \frac{4}{3} E^\* R^\* (d - D\_{ij})^{1.5} & D\_{ij} < d \\ 0 & \text{else} \end{cases} \tag{31}$$

with

$$E^\* = \frac{E}{2(1 - \nu^2)} \quad \text{and} \quad R^\* = \sqrt{\frac{d}{4}} \tag{32}$$

where *E* is Young's modulus and *ν* denotes the Poisson ratio. Strictly, the contact force varies slightly at fiber ends due to different contact areas. However, the exact pressure distribution in the contact area is not the focus of this work and for high fiber aspect ratios, the portion of fiber end contact becomes relatively small. Thus, the contact force at fiber ends may be rather interpreted as a penetration penalty.

Finally, the particle acceleration due to contact forces is computed

$$\mathbf{a}\_{i}^{\text{contact}} = \sum\_{j} \frac{w\_{ij} F\_{ij}^{\text{c}}}{m\_{i}} \mathbf{n}\_{ij} \tag{33}$$

 as

for each proximity point *j* with the contact normal **<sup>n</sup>***ij* and a weighting factor *wij*. The weighting factor is necessary to distribute the force at a contact point between particles associated to this contact point and is defined as

$$w\_{ij} = \begin{cases} \frac{d\_{\mathbb{P}} - \mathcal{P}\_{ij}}{d\_{\mathbb{P}}} & 0 < \mathcal{P}\_{ij} < d\_{\mathbb{P}}\\ \frac{d\_{\mathbb{n}} + \mathcal{P}\_{ij}}{d\_{\mathbb{n}}} & -d\_{\mathbb{n}} < \mathcal{P}\_{ij} \le 0\\ 1 & \text{two fiber ends} \\ 0 & \text{else} \end{cases} \tag{34}$$

with the distance to the previous fiber particle *d*p = **x***i* − **<sup>x</sup>***i*−<sup>1</sup> and next particle *d*n = **x***i* − **<sup>x</sup>***i*+<sup>1</sup> . Friction in tangential direction is neglected and fiber surfaces are assumed to be smooth. However, friction could be easily incorporated at this point, if the friction coefficient is available. The acceleration is computed for all particles that might possibly contact any other particle and this way, equal force magnitudes on destination and source fibers are ensured.

#### *2.6. Time Integration and Implementation*

The time integration is an extension of the kick-drift scheme used in the original transport velocity formulation [43]. The advection velocities are computed for a half step as

$$\mathbf{v}\_{i}^{n+\frac{1}{2}} = \mathbf{v}\_{i}^{n} + \frac{\Delta t}{2} \left[ \mathbf{\bar{a}}\_{i} + \mathbf{g} \right]\_{n-\frac{1}{2}} \tag{35} \\ \tag{35} \\ \tag{35}$$

$$\mathbf{v}\_{i}^{n+\frac{1}{2}} = \mathbf{v}\_{i}^{n+\frac{1}{2}} + \frac{\Delta t}{2} \left[ \frac{p\_{\rm b}}{m\_{i}} \sum\_{j \in \Omega\_{\rm m}} \left( V\_{i}^{2} + V\_{j}^{2} \right) \nabla\_{i} \mathcal{W}\_{ij} \right]\_{n-\frac{1}{2}} \qquad \qquad , \quad i \in \Omega\_{\rm m} \tag{36}$$

$$\mathbf{\tilde{v}}\_{i}^{n+\frac{1}{2}} = \mathbf{\tilde{v}}\_{i}^{n} + \frac{\Delta t}{2} \left[ \mathbf{a}\_{i}^{\text{hydro}} + \mathbf{a}\_{i}^{\text{elastic}} + \mathbf{a}\_{i}^{\text{transition}} + \mathbf{a}\_{i}^{\text{contact}} + \mathbf{g} \right]\_{n-\frac{1}{2}} \qquad \qquad , \quad i \in \Omega\_{\text{f}} \tag{37}$$

based on the previous accelerations at step *n* − 12 . Equation (36) utilizes the background pressure *p*b to move fluid particles such that no agglomerations or voids form. This is done by modifying the momentum velocity **v***n*<sup>+</sup> 12 *i* to the advection velocity **v**˜ *n*+ 12 *i* . The difference in momentum velocity and advection velocity is corrected by the artificial stress *σ*A*i* in the momentum balance. Fiber particles do not experience a background pressure, as they should not be used to fill voids etc. and thus, their advection velocity is equal to their physical velocity. Then, all particles are moved according to

$$\mathbf{x}\_{i}^{n+1} = \mathbf{x}\_{i}^{n} + \Delta t \mathbf{\bar{v}}\_{i}^{n+\frac{1}{2}} \tag{38} \\ \tag{39}$$

for the full step. Finally, the velocities are updated as

$$\mathbf{v}\_{i}^{n+1} = \mathbf{v}\_{i}^{n+\frac{1}{2}} + \frac{\Delta t}{2} \left[ \mathbf{\bar{a}}\_{i} + \mathbf{g} \right]\_{n+\frac{1}{2}} \tag{39}$$

$$\mathbf{v}\_{i}^{n+1} = \mathbf{\tilde{v}}\_{i}^{n+\frac{1}{2}} + \frac{\Delta t}{2} \left[ \mathbf{a}\_{i}^{\text{hydro}} + \mathbf{a}\_{i}^{\text{elastic}} + \mathbf{a}\_{i}^{\text{transition}} + \mathbf{a}\_{i}^{\text{contact}} + \mathbf{g} \right]\_{n+\frac{1}{2}} \qquad \qquad , \quad i \in \Omega\_{\text{f}}.\tag{40}$$

As this is an explicit time integration scheme, it is only conditionally stable. The maximum time step is computed by

$$
\Delta t = \min \left( 0.4 \frac{\Delta x}{1.1 \varepsilon\_0}, 0.125 \frac{\Delta x^2}{\eta}, 0.25 \sqrt{\frac{\Delta x}{g}}, 0.5 \Delta x \sqrt{\frac{\rho\_0}{E}}, 0.5 \Delta x^2 \sqrt{\frac{\rho\_0 A}{2EI}} \right) \tag{41}
$$

as the minimum of the CFL condition, a viscous condition, a body force condition, a tensile elastic condition and a elastic bending condition. The implementation was realized in PySPH [49] due to its flexibility in implementing the additional equations on top of the transport velocity scheme. The code is publicly available as fork of the original PySPH project (https://github.com/nilsmeyerkit/pysph/ tree/fibers).
