Parameter Identification of Fiber Orientation Models Based on Direct Fiber Simulation with Smoothed Particle Hydrodynamics

Abstract

The behavior of fiber suspensions during flow is of fundamental importance to the process simulation of discontinuous fiber reinforced plastics. However, the direct simulation of flexible fibers and fluid poses a challenging two-way coupled fluid-structure interaction problem. Smoothed Particle Hydrodynamics (SPH) offers a natural way to treat such interactions. Hence, this work utilizes SPH and a bead chain model to compute a shear flow of fiber suspensions. The introduction of a novel viscous surface traction term is key to achieve full agreement with Jeffery’s equation. Careful modelling of contact interactions between fibers is introduced to model suspensions in the non-dilute regime. Finally, parameters of the Reduced-Strain Closure (RSC) orientation model are identified using ensemble averages of multiple SPH simulations implemented in PySPH and show good agreement with literature data.

1. Introduction

Molding of discontinuous reinforced polymer fiber suspensions, e.g., glass fibers in polymer melt, leads to fiber reorientation. Understanding and predicting the behavior of such fiber suspensions is crucial for achieving high quality products in common composite manufacturing processes such as injection molding and compression molding. Due to the high cost of production facilities and molds, it is desirable to simulate the flow in a computationally efficient and reliable way before running experiments. The fiber reorientation is also of high interest for subsequent anisotropic structural simulations in the framework of a continuous CAE chain [1,2].

Jeffery [3] was the first to analytically derive the motion of a single rigid, ellipsoidal shaped body in a viscous Newtonian flow without buoyancy or inertia. Jeffery’s work was extended by Folgar and Tucker [4] to account for fiber interactions by introducing a fiber interaction coefficient that models isotropic diffusivity of fiber orientation. Other phenomenological parameters were introduced to capture experimentally observed orientation delays in the SRF and RSC model [5] and to account for anisotropic diffusion in the Anisotropic Rotary Diffusion (ARD) model [6] or the improved Anisotropic Rotary Diffusion (iARD) model [7]. These phenomenological parameters account for interactions of multiple fibers in non-dilute suspensions and are fitted to experimental observations, but do not describe a two-phase suspension. Instead, they model the fiber orientation state with second and fourth order moments of the fiber orientation distribution function $\Psi \left(\mathbf{p}\right)$ introduced by Advani and Tucker [8] as fiber orientation tensors

$$\mathbf{A}={\int }_S\mathbf{p}\otimes \mathbf{p}\Psi \left(\mathbf{p}\right)\mathrm{d}p$$

(1)

and

$$\mathbb{A}={\int }_S\mathbf{p}\otimes \mathbf{p}\otimes \mathbf{p}\otimes \mathbf{p}\Psi \left(\mathbf{p}\right)\mathrm{d}p.$$

(2)

Here, $\mathbf{p}$ describes a fiber direction and $\mathrm{d}p$ is the surface element on a unit sphere $S:=\{\mathbf{p}\in {\mathbb{R}}^3:∥\mathbf{p}∥=1\}$. Typically, a closure approach for the fourth order fiber orientation tensor is employed to describe the evolution of the second order fiber orientation tensor. Whenever this work utilizes a closure approximation, the invariant-based optimal fitting (IBOF) closure is chosen [9]. Determining the parameters in these macroscopic models from experiments can be time consuming. Thus, a direct fiber simulation may be utilized to identify these parameters.

1.1. Point-Wise Interaction Methods

A common approach for the simulation of fiber suspensions is the treatment of fibers as slender bodies that interact with the fluid at discrete points. Exact solutions from Stokesian dynamics [10] or slender-body theory [11] are utilized to describe long-range hydrodynamic interaction between fluid and particles for creeping flows. Several authors proposed models for single flexible fibers suspended in a fluid. Hinch [12] started by modeling inextensible threads. Yamamoto et al. [13,14] developed a fiber model consisting of individual beads that experience Stokesian drag forces from the fluid. These bead chain models are computationally expensive and authors have combined multiple beads to rods leading to more efficient rod chains [15,16,17]. Alternatively, spheres [18,19] and spheroids [20] connected with hinges were suggested. Lindström and Uesaka [21] use a discrete field of point forces to ensure that the fluid is experiencing forces opposed to the fiber (two-way coupling). Several authors investigated rheological properties of multiple rigid fibers suspended in a fluid based on these models for single fibers. Yamane et al. [22] described the motion of multiple rigid rods with hydrodynamic drag force and torque on the individual rods. They added a lubrication force for rods that are in close proximity to each other in order to capture short range hydrodynamic effects. However, they did not account for long range hydrodynamic interaction between particles, which was addressed later by Fan et al. [23] using slender-body theory. Sundararajakumar and Koch [24] showed that lubrication forces alone do not prevent penetration of fibers and included contact forces. Suspensions of multiple flexible fibers were investigated with rod-chain models [25] and simplified bead chain models [26]. In addition to these rheological investigations, few direct simulations have been applied on component scale with flexible fibers or bundles to investigate effects such as fiber-matrix separation [27,28,29,30].

1.2. Resolved Methods

In contrast to models based on discrete point-wise interaction forces, the suspension flow might be also described directly by a two-phase flow, in which one phase represents fibers and the other one the suspending fluid. Solving this fully coupled two-phase system with mesh-based approaches such as Finite Elements or Finite Volumes raises the problem of large mesh deformations. A simulation using an immersed boundary method resolves the remeshing problem but requires a mesh resolution significantly smaller than the fiber diameter to track the interface and is therefore computationally expensive. An alternative approach for the simulation of two-phase fiber suspensions is the use of particle methods. Meshless particle methods can offer a natural way to represent fluid-structure interaction at the fiber surface and have been studied less than point-wise interaction methods. Bian and Ellero propose a splitting scheme for separate integration of long-range hydrodynamic interactions and short-range lubrication [31], which was applied in an SPH simulation of a concentrated spherical particle suspension [32]. The fine resolution of suspended particles leads to a high accuracy, but comes at high computational costs. Yashiro et al. [33,34] connected particles to rigid bodies for modeling the injection molding of dilute short-fiber-reinforced composites using a moving particle semi-implicit method. Recently, a very similar approach using SPH was reported by He et al. [35,36]. However, both showed only short time periods when comparing the simulation to Jeffery’s equation and did not investigate the interaction between fibers. Fiber suspensions were also investigated in the framework of dissipative particle dynamics with a focus on Boger fluids [37]. This investigation used a Langrange multiplier to constrain fiber extension and multiple parameters had to be calibrated to analytical solutions in order to achieve correct hydrodynamic forces on the fiber. Yang et al. [38] employed the SPH method in order to simulate a single flexible fiber in a viscous fluid with a focus on determining its drag coefficient. However, they limited their work to a single fiber without interactions.

The present work extends Yang’s model with a viscous surface traction term necessary to meet Jeffery’s equation precisely and allows for the simulation of non-dilute fiber suspensions through the introduction of fiber contact forces. First, rotation and bending of a single fiber is compared to Jeffery’s equation to validate the approach. Then, parameters of the RSC fiber orientation model are determined using SPH simulations.

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 $\alpha \left(x\right)$ in an n-dimensional domain $\Omega \subset {\mathbb{R}}^n$ at a position $\mathbf{x},{\mathbf{x}}^{\prime }\in \Omega$ can be described using an integral of the form

$$\alpha \left(\mathbf{x}\right)={\int }_{\Omega }\alpha \left({\mathbf{x}}^{\prime }\right)\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\mathrm{d}{\mathbf{x}}^{\prime }$$

(3)

employing the Dirac distribution $\delta \left(\mathbf{x}\right)$. 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 $W(\mathbf{x}-{\mathbf{x}}^{\prime },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

$$W={\beta }_n\left\{\begin{array}{cc}{(2-q)}^3-4{(1-q)}^3\hfill & 0\le q<1\hfill \\ {(2-q)}^3\hfill & 1\le q<2\hfill \\ 0\hfill & 2\le q\hfill \end{array}\right.$$

(4)

with $q=\frac{\parallel \mathbf{x}-{\mathbf{x}}^{\prime }\parallel }{h}$ and a dimension-dependent normalization factor ${\beta }_n$. Employing the kernel $W({\mathbf{x}}_i-{\mathbf{x}}_j,h)=W_{ij}$, Equation (3) can be approximated as

$$\alpha \left(\mathbf{x}\right)\approx \check{\alpha }\left({\mathbf{x}}_i\right)={\check{\alpha }}_i=\sum _{j\in \Omega }{\alpha }_j{W}_{ij}\frac{m_j}{{\rho }_j}$$

(5)

denoting a particle’s property $\alpha$ at position $x_j$ as ${\alpha }_j$ and by expressing the associated volume as ratio between mass $m_j$ and density ${\rho }_j$ of the particle. The key is the usage of an analytically differentiable kernel such that

$$\frac{\partial {\check{\alpha }}_i}{\partial \mathbf{x}}=\sum _{j\in \Omega }{\alpha }_j\frac{m_j}{{\rho }_j}{\nabla }_i{W}_{ij}$$

(6)

can be computed analytically with the kernel gradient $\frac{\partial W}{\partial \mathbf{x}}(\mathbf{x}-{\mathbf{x}}_j,h){|}_i={\nabla }_i{W}_{ij}$, hence turning the PDE to an ODE. The domain $\Omega ={\Omega }_{\mathrm{m}}\cup {\Omega }_{\mathrm{s}}$ is the set of all particles in the model with subsets for the fluid particles ${\Omega }_{\mathrm{m}}$ and solid particles ${\Omega }_{\mathrm{s}}={\Omega }_{\mathrm{w}}\cup {\Omega }_{\mathrm{f}}$ consisting of wall particles ${\Omega }_{\mathrm{w}}$ and fiber particles ${\Omega }_{\mathrm{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 $\mathbf{v}$ and advection velocity $\tilde{\mathbf{v}}$ leading to a Navier-Stokes equation of the form

$$\frac{\tilde{\mathrm{d}}\left(\rho \mathbf{v}\right)}{\mathrm{d}t}=-\nabla p+\eta {\nabla }^2\mathbf{v}+\rho \mathbf{g}+\nabla ·(\rho \mathbf{v}\otimes (\tilde{\mathbf{v}}-\mathbf{v}))$$

(7)

with density $\rho$, pressure p, viscosity $\eta$ and body accelerations $\mathbf{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_{\mathrm{b}}$ that effectively suppresses tensile instability, but does not influence the actual momentum. Hence, the term ${\sigma }^{\mathrm{A}}=\rho \mathbf{v}\otimes (\tilde{\mathbf{v}}-\mathbf{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{array}{cc}\hfill {\tilde{\mathbf{a}}}_i& =\frac{1}{m_i}\sum _{j\in \Omega }\left(V_i^2+V_j^2\right)\left(-\frac{{\rho }_j{p}_j+{\rho }_i{p}_i}{{\rho }_i+{\rho }_j}{\nabla }_i{W}_{ij}\right)\hfill \\ \hfill & +\frac{1}{m_i}\sum _{j\in {\Omega }_s}\left(V_i^2+V_j^2\right)\left(\frac{2{\eta }_i{\eta }_j}{{\eta }_i+{\eta }_j}\frac{{\nabla }_i{W}_{ij}·({\mathbf{x}}_i-{\mathbf{x}}_j)}{\parallel {\mathbf{x}}_i-{\mathbf{x}}_j\parallel +\epsilon }({\mathbf{v}}_i-{\widehat{\mathbf{v}}}_j)\right)\hfill \\ \hfill & +\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{W}_{ij}·({\mathbf{x}}_i-{\mathbf{x}}_j)}{\parallel {\mathbf{x}}_i-{\mathbf{x}}_j\parallel +\epsilon }({\mathbf{v}}_i-{\mathbf{v}}_j)+\frac{1}{2}({\mathit{\sigma }}_i^{\mathrm{A}}+{\mathit{\sigma }}_j^{\mathrm{A}}){\nabla }_i{W}_{ij}\right),i\in {\Omega }_{\mathrm{m}}\hfill \end{array}$$

(8)

with the volume attributed to each particle $V_i=m_i/{\rho }_i$ and the velocity of an adjacent solid ${\widehat{\mathbf{v}}}_j$, which is explained in the next section. Indices are used to refer to the particle based density ${\rho }_i$, pressure $p_i$, viscosity ${\eta }_i$, mass $m_i$, position ${\mathbf{x}}_i$ and velocity ${\mathbf{v}}_i$. The parameter $\epsilon$ 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 }_i=m_i\sum _{j\in \Omega }{W}_{ij},i\in {\Omega }_{\mathrm{m}}$$

(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),i\in {\Omega }_{\mathrm{m}},$$

(10)

where ${\rho }_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=\frac{{\rho }_0{c}^2}{\gamma }$ 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 $\gamma$ is set as $\gamma =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 }_{\mathrm{m}}}{\rho }_j{W}_{ij}+(\mathbf{g}-{\mathbf{a}}_i)·{\sum }_{j\in {\Omega }_{\mathrm{m}}}{\rho }_j({\mathbf{x}}_i-{\mathbf{x}}_j)W_{ij}}{{\sum }_{j\in {\Omega }_{\mathrm{m}}}{W}_{ij}},i\in {\Omega }_{\mathrm{s}}$$

(11)

with a prescribed acceleration of the solid ${\mathbf{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

$${\widehat{\mathbf{v}}}_i=2v_i-\frac{{\sum }_{j\in {\Omega }_{\mathrm{m}}{v}_j{W}_{ij}}}{{\sum }_{j\in {\Omega }_{\mathrm{m}}{W}_{ij}}},i\in {\Omega }_{\mathrm{s}},$$

(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=\mathrm{const},{\mathbf{a}}_i=0,i\in {\Omega }_{\mathrm{w}}.$$

(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}}_i^{\mathrm{hydro}}=\frac{1}{m_i}\sum _{j\in {\Omega }_f}\left(V_i^2+V_j^2\right)\left(-\frac{{\rho }_j{p}_j+{\rho }_i{p}_i}{{\rho }_i+{\rho }_j}{\nabla }_i{W}_{ij}+\frac{2{\eta }_i{\eta }_j}{{\eta }_i+{\eta }_j}\frac{{\nabla }_i{W}_{ij}·({\mathbf{x}}_i-{\mathbf{x}}_j)}{\parallel {\mathbf{x}}_i-{\mathbf{x}}_j\parallel +\epsilon }({\widehat{\mathbf{v}}}_i-{\mathbf{v}}_j)\right),i\in {\Omega }_{\mathrm{f}}$$

(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 $\Delta x$ is directly related to the fiber diameter by

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

(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}^{\mathrm{t}}=EA\left(\frac{\Vert {\mathbf{x}}_{ij}\Vert }{x_{ij}^0}-1\right)\frac{{\mathbf{x}}_{ij}}{\Vert {\mathbf{x}}_{ij}\Vert }$$

(16)

and bending moment

$${\mathbf{M}}_i^{\mathrm{b}}=\frac{EI}{2}\frac{{\mathit{\theta }}_i-{\mathit{\theta }}_i^0}{\Vert {\mathbf{x}}_{i(i-1)}\Vert +\Vert {\mathbf{x}}_{i(i+1)}\Vert }$$

(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 ${\mathbf{x}}_{ij}$ and the enclosed angle is denoted as ${\mathit{\theta }}_i$. Here, the vector notation of ${\mathit{\theta }}_i$ indicates that its direction resembles the axis of rotation. The undeformed reference configuration is denoted with a superscript ${(·)}^0$ in both, Equations (16) and (17). The bending moment can be converted to pairs of forces

$${\mathbf{F}}_{ij}^{\mathrm{b}}=\frac{1}{2}\frac{{\mathbf{M}}_i^{\mathrm{b}}\times {\mathbf{x}}_{ij}}{\Vert {\mathbf{x}}_{ij}{\Vert }^2}$$

(18)

that act on the particle and its neighbors. Figure 1 illustrates these forces and it can been seen, that generally ${\mathbf{F}}_{ij}^{\mathrm{b}}\ne {\mathbf{F}}_{ji}^{\mathrm{b}}$. 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

$$\begin{array}{cc}\hfill {\mathbf{a}}_i^{\mathrm{elastic}}& =\frac{1}{m_i}({\mathbf{F}}_{i(i+1)}^{\mathrm{t}}+{\mathbf{F}}_{i(i-1)}^{\mathrm{t}}\hfill \\ \hfill & +{\mathbf{F}}_{i(i-1)}^{\mathrm{b}}+{\mathbf{F}}_{i(i+1)}^{\mathrm{b}}-{\mathbf{F}}_{(i-1)i}^{\mathrm{b}}-{\mathbf{F}}_{(i+1)i}^{\mathrm{b}}).\hfill \end{array}$$

(19)

If the angle between two adjacent particle pairs exceeds a certain critical value ${\theta }_{\mathrm{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 segment of length $\Delta L$ and diameter d with the orientation of its centerline $\mathbf{p}$. One exemplary surface normal $\mathbf{n}$ is shown with its parametrization angle $\psi$.

The fiber direction $\mathbf{p}$ and any arbitrary surface normal ${\mathbf{n}}^0$ are perpendicular unit-vectors. Any other surface normal can be constructed from this arbitrary normal by rotating it around $\mathbf{p}$. The surface normal can be parameterized using ${\mathbf{n}}^0$ and $\psi$ employing a rotational tensor $\mathbf{R}$ around axis $\mathbf{p}$ [46]. The parameterized normal becomes

$$\mathbf{n}\left(\psi \right)=\mathbf{R}\left(\psi \right){\mathbf{n}}_0.$$

(20)

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

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

(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 $\mathbf{t}$ with its corresponding leverage $\frac{d}{2}\mathbf{n}\left(\psi \right)$ as

$${\mathbf{M}}^{\mathrm{v}}=\Delta L{\int }_0^{2\pi }\frac{d}{2}\mathbf{n}\left(\psi \right)\times \mathbf{t}\left(\psi \right)\frac{d}{2}\mathrm{d}\psi$$

(22)

with a constant fiber diameter d. The term $\frac{d}{2}\mathrm{d}\psi$ represents an infinitesimal circumferential line segment on the cylinder surface. As the cross product of a vector with itself vanishes, this can be simplified to

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

(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^{\mathrm{traction}}=\frac{1}{2}\frac{{\mathbf{M}}_i^{\mathrm{v}}\times {\mathbf{x}}_{ij}}{J},j\in [i-1,i+1]$$

(24)

Here, J denotes the moment of inertia for a cylindrical body around its first principle axis of inertia and the factor $\frac{1}{2}$ 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 $\Delta L=3\Delta 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{W}_{ij}.$$

(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.

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}=\left[[{\mathbf{p}}_i\times {\mathbf{p}}_j]\right]$$

(26)

of the involved fiber direction vectors ${\mathbf{p}}_i$ and ${\mathbf{p}}_j$. The operator $\left[[·]\right]=\frac{(·)}{\parallel ·\parallel }$ is used to conveniently denote the normalization of a vector. Solving the small linear system of equations

$$[{\mathbf{p}}_i,{\mathbf{n}}_{ij},-{\mathbf{p}}_j]{[{\mathcal{P}}_{ij},D_{ij},{\mathcal{P}}_{ji}]}^{\top }={\mathbf{x}}_{ij}$$

(27)

with its adjugate matrix leads to the solution for the distance between the fibers $D_{ij}$ and the projections to source and destination vectors ${\mathcal{P}}_{ij}$ and ${\mathcal{P}}_{ji}$, respectively.

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 ${\mathbf{x}}_{ij}$ can be used to obtain the normal direction by projection. It is assumed that ${\mathbf{p}}_i$ describes a unit vector in fiber direction at one fiber particle at position $\mathcal{A}$. Let ${\mathbf{x}}_{ij}$ be the vector from another fibers’ end particle $\mathcal{B}$ to the point $\mathcal{A}$. The closest point to $\mathcal{B}$ on a line with direction ${\mathbf{p}}_i$ is denoted as $\mathcal{C}$ and can be used to define the normal direction as

$${\mathbf{n}}_{ij}=\left[\left[\overrightarrow{\mathcal{B}\mathcal{C}}\right]\right]=\left[[\overrightarrow{\mathcal{A}\mathcal{C}}-\overrightarrow{\mathcal{A}B}]\right].$$

(28)

Using the definitions above and the fact that $\mathcal{C}$ is the projection of $\mathcal{B}$ to the line with direction ${\mathbf{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].$$

(29)

The projection on the destination fiber is given as ${\mathcal{P}}_{ij}=-{\mathbf{p}}_i·{\mathbf{x}}_{ij}$ and the contact distance is computed as $D_{ij}=\parallel {\mathbf{x}}_{ij}+{\mathcal{P}}_{ij}{\mathbf{p}}_i\parallel$.

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}=\left[\left[{\mathbf{x}}_{ij}\right]\right]$$

(30)

with the corresponding distance $D_{ij}=\parallel {\mathbf{x}}_{ij}\parallel$.

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}^{\mathrm{c}}=\left\{\begin{array}{cc}\frac{4}{3}{E}^{\ast }{R}^{\ast }{(d-D_{ij})}^{1.5}\hfill & D_{ij}<d\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.$$

(31)

with

$$E^{\ast }=\frac{E}{2(1-{\nu }^2)}\mathrm{and}{R}^{\ast }=\sqrt{\frac{d}{4}}$$

(32)

where E is Young’s modulus and $\nu$ 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 as

$${\mathbf{a}}_i^{\mathrm{contact}}=\sum _j\frac{w_{ij}{F}_{ij}^{\mathrm{c}}}{m_i}{\mathbf{n}}_{ij}$$

(33)

for each proximity point j with the contact normal ${\mathbf{n}}_{ij}$ and a weighting factor $w_{ij}$. 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}=\left\{\begin{array}{cc}\frac{d_{\mathrm{p}}-{\mathcal{P}}_{ij}}{d_{\mathrm{p}}}\hfill & 0<{\mathcal{P}}_{ij}<d_{\mathrm{p}}\hfill \\ \frac{d_{\mathrm{n}}+{\mathcal{P}}_{ij}}{d_{\mathrm{n}}}\hfill & -d_{\mathrm{n}}<{\mathcal{P}}_{ij}\le 0\hfill \\ 1\hfill & \mathrm{two}\mathrm{fiber}\mathrm{ends}\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.$$

(34)

with the distance to the previous fiber particle $d_{\mathrm{p}}=\parallel {\mathbf{x}}_i-{\mathbf{x}}_{i-1}\parallel$ and next particle $d_{\mathrm{n}}=\parallel {\mathbf{x}}_i-{\mathbf{x}}_{i+1}\parallel$. 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

$$\begin{array}{ccc}\hfill {\mathbf{v}}_i^{n+\frac{1}{2}}& ={\mathbf{v}}_i^n+\frac{\Delta t}{2}{\left[{\tilde{\mathbf{a}}}_i+\mathbf{g}\right]}_{n-\frac{1}{2}}\hfill & \hfill ,i\in {\Omega }_{\mathrm{m}}\end{array}$$

(35)

$$\begin{array}{ccc}\hfill {\tilde{\mathbf{v}}}_i^{n+\frac{1}{2}}& ={\mathbf{v}}_i^{n+\frac{1}{2}}+\frac{\Delta t}{2}{\left[\frac{p_{\mathrm{b}}}{m_i}\sum _{j\in {\Omega }_{\mathrm{m}}}\left(V_i^2+V_j^2\right){\nabla }_i{W}_{ij}\right]}_{n-\frac{1}{2}}\hfill & \hfill ,i\in {\Omega }_{\mathrm{m}}\end{array}$$

(36)

$$\begin{array}{ccc}\hfill {\tilde{\mathbf{v}}}_i^{n+\frac{1}{2}}& ={\tilde{\mathbf{v}}}_i^n+\frac{\Delta t}{2}{\left[{\mathbf{a}}_i^{\mathrm{hydro}}+{\mathbf{a}}_i^{\mathrm{elastic}}+{\mathbf{a}}_i^{\mathrm{traction}}+{\mathbf{a}}_i^{\mathrm{contact}}+\mathbf{g}\right]}_{n-\frac{1}{2}}\hfill & \hfill ,i\in {\Omega }_{\mathrm{f}}\end{array}$$

(37)

based on the previous accelerations at step $n-\frac{1}{2}$. Equation (36) utilizes the background pressure $p_{\mathrm{b}}$ to move fluid particles such that no agglomerations or voids form. This is done by modifying the momentum velocity ${\mathbf{v}}_i^{n+\frac{1}{2}}$ to the advection velocity ${\tilde{\mathbf{v}}}_i^{n+\frac{1}{2}}$. The difference in momentum velocity and advection velocity is corrected by the artificial stress ${\mathit{\sigma }}_i^{\mathrm{A}}$ 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

$$\begin{array}{ccc}\hfill {\mathbf{x}}_i^{n+1}& ={\mathbf{x}}_i^n+\Delta t{\tilde{\mathbf{v}}}_i^{n+\frac{1}{2}}\hfill & \hfill ,i\in \Omega \end{array}$$

(38)

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

$$\begin{array}{ccc}\hfill {\mathbf{v}}_i^{n+1}& ={\mathbf{v}}_i^{n+\frac{1}{2}}+\frac{\Delta t}{2}{\left[{\tilde{\mathbf{a}}}_i+\mathbf{g}\right]}_{n+\frac{1}{2}}\hfill & \hfill ,i\in {\Omega }_{\mathrm{m}}\end{array}$$

(39)

$$\begin{array}{ccc}\hfill {\tilde{\mathbf{v}}}_i^{n+1}& ={\tilde{\mathbf{v}}}_i^{n+\frac{1}{2}}+\frac{\Delta t}{2}{\left[{\mathbf{a}}_i^{\mathrm{hydro}}+{\mathbf{a}}_i^{\mathrm{elastic}}+{\mathbf{a}}_i^{\mathrm{traction}}+{\mathbf{a}}_i^{\mathrm{contact}}+\mathbf{g}\right]}_{n+\frac{1}{2}}\hfill & \hfill ,i\in {\Omega }_{\mathrm{f}}.\end{array}$$

(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.1c_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)$$

(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).

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.2L_{\mathrm{f}}$ and the dimension in flow direction is $L=2B$ with the length of a fiber $L_{\mathrm{f}}$. The Reynolds number

$$Re=\frac{{\rho }_0{L}_{\mathrm{f}}\frac{GB}{2}}{\eta }$$

(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 G{r}_{\mathrm{e}}^4}$$

(43)

for a given shear rate G and ellipsoidal aspect ratio $r_{\mathrm{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}=\mathit{\omega }\mathbf{p}+\xi \left(\mathbf{D}\mathbf{p}-\left(\mathbf{p}\otimes \mathbf{p}\otimes \mathbf{p}\right)\mathbf{D}\right)$$

(44)

with vorticity tensor $\mathit{\omega }$ and symmetric strain rate tensor $\mathbf{D}$. The shape factor $\xi$ is an alternative measure for the (equivalent) ellipsoidal aspect ratio $r_{\mathrm{e}}$ and is defined as

$$\xi \left(r_{\mathrm{e}}\right)=\frac{r_{\mathrm{e}}^2-1}{r_{\mathrm{e}}^2+1}.$$

(45)

The solution to Equation (44) is periodic with

$$T=\frac{2\pi }{G}\left(r_{\mathrm{e}}+\frac{1}{r_{\mathrm{e}}}\right)$$

(46)

being the time for a full rotation of a fiber. When solving Jeffery’s equation for a cylinder with geometric aspect ratio $r_{\mathrm{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_{\mathrm{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_{\mathrm{e}}\left(r_{\mathrm{p}}\right)=0.000035r_{\mathrm{p}}^3-0.00467r_{\mathrm{p}}^2+0.764r_{\mathrm{p}}+0.404$$

(47)

is used here.

Figure 5 depicts the orientation $\varphi$ of a single fiber in shear flow for fiber length $L_{\mathrm{f}}=5\Delta x$ and $L_{\mathrm{f}}=11\Delta 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.

Bending modes were shown in Yamamoto and Matsuka’s [13] numerical results and the corresponding SPH simulation in

Table 1. Bending modes of a fiber with length $L_{\mathrm{f}}=11\Delta x$ for varied dimensionless stiffness. One example with $S=10$ and critical fiber bending angle ${\theta }_{\mathrm{c}}=\frac{\pi }{4}$ is shown to demonstrate fiber fracture.

Strain	0	$\frac{1}{4}\mathit{TG}$	$\frac{1}{2}\mathit{TG}$	$\frac{3}{4}\mathit{TG}$	$\mathit{TG}$
$S=5,{\theta }_{\mathrm{c}}=\infty$
$S=10,{\theta }_{\mathrm{c}}=\infty$
$S=10,{\theta }_{\mathrm{c}}=\frac{\pi }{4}$
$S=20,{\theta }_{\mathrm{c}}=\infty$
$S=100,{\theta }_{\mathrm{c}}=\infty$

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.

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_{\mathrm{I}}$ and a phenomenological factor $\kappa$ 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}=\mathit{\omega }\mathbf{A}-\mathbf{A}\mathit{\omega }+\xi \left(\mathbf{D}\mathbf{A}+\mathbf{A}\mathbf{D}-2\left(\mathbb{A}+(1-\kappa )(\mathbb{L}-\mathbb{M}\mathbb{A})\right)\left[\mathbf{D}\right]\right)+2\kappa C_{\mathrm{I}}G\left(\mathbf{1}-3\mathbf{A}\right)$$

(48)

with $\mathbb{L}={\sum }_{i=1}^3{\lambda }_i{\mathbf{e}}_i\otimes {\mathbf{e}}_i\otimes {\mathbf{e}}_i\otimes {\mathbf{e}}_i$ and $\mathbb{M}={\sum }_{i=1}^3{\mathbf{e}}_i\otimes {\mathbf{e}}_i\otimes {\mathbf{e}}_i\otimes {\mathbf{e}}_i$ using the eigenvalues ${\lambda }_i$ and eigenvectors ${\mathbf{e}}_i$ of the second order fiber orientation tensor $\mathbf{A}$. Setting $\kappa =1$ would reduce this model to the Folgar-Tucker model and setting $C_{\mathrm{I}}=0$ reduces it to Jeffery’s Equation (44). The choice of feasible parameters $C_{\mathrm{I}}\in [0,0.1]$ and $\kappa \in [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\Delta 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 $x_1$-direction according to

$$x_2^{\prime }=\left(x_2\pm GLt\right)modL$$

(49)

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

$$v_2^{\prime }=\left(v_2\pm GL\right)$$

(50)

for the shift. Conventional periodic boundary conditions apply in $x_2$- and $x_3$-direction, as depicted in Figure 6. All fibers are initially oriented in $x_1$-direction and randomly positioned in the cube. This unidirectional arrangement 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 ${\langle ·\rangle }_{t_n}$ is defined as the mean across multiple realizations at the same time step $t_n$.

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

$$\underset{C_i,\kappa }{\mathrm{minimize}}\sum _n{∥\mathbf{A}\left(t_n\right)-{\langle \tilde{\mathbf{A}}\rangle }_{t_n}∥}^2$$

(51)

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

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

(52)

with each fiber’s orientation ${\mathbf{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 $\mathbf{A}$. However, the method is in no way limited to rigid fibers.

Figure 7 shows the non-trivial components of the ensemble average ${\langle \tilde{\mathbf{A}}\rangle }_{t_n}$ 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 $\mathbf{A}$ computed with optimal parameters according to the RSC model given in Equation (48). The simulation time is $10T$, 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 agreement 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 $x_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 $x_1{x}_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.

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 agreement with literature data and support the use of this approach for parameter identification.

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 agreement 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_{\mathrm{f}}=5\Delta 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.