Next Article in Journal
Progress of Bio-Calcium Carbonate Waste Eggshell and Seashell Fillers in Polymer Composites: A Review
Next Article in Special Issue
Parameter Identification of Fiber Orientation Models Based on Direct Fiber Simulation with Smoothed Particle Hydrodynamics
Previous Article in Journal
Remanufacturing of Woven Carbon Fibre Fabric Production Waste into High Performance Aligned Discontinuous Fibre Composites
Previous Article in Special Issue
Experimental Validation of a Direct Fiber Model for Orientation Prediction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

Fiber Orientation Predictions—A Review of Existing Models

1
Corporate Sector Research and Advance Engineering, Robert Bosch GmbH, 71272 Renningen, Germany
2
Polymer Engineering Center, University of Wisconsin-Madison, Madison, WI 53706, USA
*
Author to whom correspondence should be addressed.
J. Compos. Sci. 2020, 4(2), 69; https://doi.org/10.3390/jcs4020069
Submission received: 11 May 2020 / Revised: 29 May 2020 / Accepted: 2 June 2020 / Published: 8 June 2020
(This article belongs to the Special Issue Discontinuous Fiber Composites, Volume II)

Abstract

:
Fiber reinforced polymers are key materials across different industries. 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 requires therefore 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 modelling approaches. In this context, the article describes different macroscopic fiber orientation models such as the Folgar-Tucker, nematic, reduced strain closure (RSC), retarding principal rate (RPR), anisotropic rotary diffusion (ARD), principal anisotropic rotary diffusion (pARD), and Moldflow rotary diffusion (MRD) model. We discuss briefly about closure approximations, which are a common mathematical element of those macroscopic fiber orientation models. In the last section, we introduce some micro-scale numerical methods for simulating the fiber orientation phenomenon, such as the discrete element method (DEM), the smoothed particle hydrodynamics (SPH) method and the moving particle semi-implicit (MPS) method.

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 3 : = { p R 3 : | | p | | = 1 } along the fiber axis. All p S 3 can be defined by two angles ( ϕ , θ ) (Figure 1)
p 1 = sin θ cos ϕ
p 2 = sin θ sin ϕ
p 3 = cos θ
To describe the orientation of many fibers statistical methods are useful. The probability density function (PDF) ψ ( ) [1] is defined, such that ψ ( p , t ) d 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.
B ( ψ ) = [ 0 , 1 ]
ψ ( p , t ) d p = 1
ψ ( p ) = ψ ( p )
D ψ D t = s · ( ψ p ˙ ) .
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 i j = p i p j ψ ( p ) d p
a i j k l = p i p j p k p l ψ ( p ) d p .
a i n = p i p n ψ ( p ) d p .
Since the PDF is even (Equation (6)) all odd-ordered orientation tensors are zero. To simplify notation we introduce A = a i j , A = a i j k l . 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.
A i j = A j i
tr A = 1 .
The fourth order tensor is symmetric with respect to any pair of indices.
A i j k l = A j i k l = A i j k l = A k j i l = A l j k i = A i k j l = A i l k j
and all information of the second order orientation tensor can be retrieved from the fourth order tensor
A i j = A i j k k .
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 0 0 . 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:
p ˙ = W · p + ξ ( D · p D : p p p ) ,
where W = 1 2 ( ( u ) T u ) is the vorticity tensor, D = 1 2 ( ( u ) T + u ) the rate of strain tensor, ξ = r e 2 1 r e 2 + 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
p ˙ = W · p + ξ ( D · p D : p p p ) + C I γ ˙ · s ψ ,
where C I 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 C I is empirically determined and describes the rate of interaction. Setting C I = 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
D ψ D t = s · ( ψ ( W · p + ξ ( D · p D : p p p ) + C I γ ˙ · s ψ ) = s · ( ψ ( W · p + ξ ( D · p D : p p p ) ) ) + C I γ ˙ · s 2 ψ .
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
D A D t = A ˙ = A ˙ h + A ˙ d
A ˙ h = ( W · A A · W ) + ξ ( D · A + A · D 2 A : D )
A ˙ d = 2 C I γ ˙ ( 3 A ) . .
The influence of the phenomenological parameter C I is displayed in Figure 2. With decreasing diffusion (small C I ) 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.
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 ]
A ˙ = α ( A ˙ h + A ˙ d ) .
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 = i = 1 3 λ i e i e i with orthonormal eigenvectors. The modified growth rate has the following form
λ ˙ i RSC = κ λ ˙ i
e ˙ i RSC = e ˙ i ,
with the constant κ [ 0 , 1 ] . In the differential equation form, this equals
A ˙ = A ˙ RSC + κ A ˙ d
A ˙ RSC = W · A A · W + ξ [ D · A + A · D 2 ( A + ( 1 κ ) ( L M : A ) ) : D ]
L = i = 1 3 λ i e i e i e i e i
M = i = 1 3 e i e i e i e i .
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
λ ˙ i RPR = κ λ ˙ i + ( 1 κ ) β ( λ ˙ i 2 + λ ˙ j λ ˙ k )
e ˙ i RPR = e ˙ i ,
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).
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
A ˙ = A ˙ h + A ˙ ARD
A ˙ ARD = γ ˙ [ 2 C 2 tr ( C ) A 5 ( C · A + A · C ) + 10 A : C ] .
The different modeling approaches for the anisotropic rotary diffusion are listed below
C = C ( D , A ) = b 1 I + b 2 A + b 3 A 2 + b 4 γ ˙ D + b 5 γ ˙ 2 D 2 .
ARD model [5] with constants b i , i = 1 , , 5
C = C I 4 C M D 2 γ ˙
iARD model [4] with constants C I , C M
C = C I R A D 1 0 0 0 D 2 0 0 0 D 3 R A
pARD model [6] with constants C I , D 1 , D 2 , D 3 .
The eigenmatrix R A is defined by
A = R A λ 1 0 0 0 λ 2 0 0 0 λ 3 R A .
Tseng et al. [6] chose
D 1 = 1 , D 2 = c , D 3 = 1 c
to reduce the needed amount of parameters. This implies that the rotary diffusion factor in the first principal fiber orientation direction is defined by C I . The second principal fiber orientation direction is scaled with c · C I and in the third principal fiber orientation direction the smallest rotary diffusion factor with ( 1 c ) C I is applied. The influence od D 2 is displayed in Figure 4.
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
A ˙ = A ˙ h + A ˙ mARD
A ˙ mARD = 2 γ ˙ ( C tr ( C ) A ) .
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 D i a function of λ i or using Equation (32) with b 4 = b 5 = 0 and b i , i = 1 , 2 , 3 as scalar functions of tr A 2 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 C I on the volume fraction and aspect ratio of the fillers. Different experiments in different regimes even showed contradictory results [14,15,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 on physical descriptors. The model is defined by
A ˙ = A ˙ h + A ˙ nem
A ˙ nem = ( C I ( I 3 A ) + U 0 ( A · A A : A ) ) γ ˙ ,
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
A ˙ = A ˙ R S C + A ˙ ( p ) A R D R S C
A ˙ ( p ) A R D R S C = γ ˙ [ 2 [ C ( 1 κ ) M : C ] 2 κ ( tr C ) A 5 ( C · A + A · C ) + 10 [ A ( 1 κ ) ( L M : A ) ] : C ] .

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]
A ̲ = A 11 A 12 A 13 A 14 A 15 A 16 A 21 A 22 A 23 A 24 A 25 A 26 A 31 A 32 A 33 A 34 A 35 A 36 A 41 A 42 A 43 A 44 A 45 A 46 A 51 A 52 A 53 A 54 A 55 A 56 A 61 A 62 A 63 A 64 A 65 A 66 ,
where A α β = A i j k l for α , β { 1 , 2 , 3 , 4 , 5 , 6 } the indices α and β represent an index pair i j or k l in the following way:
α , β : 1 11 2 22 3 33 4 23 5 31 6 12 .
The first simple closure approaches have been introduced by Advani and Tucker [8]. The linear approach [8] is a summation of all products of a i j and δ i j . After applying symmetry and normalization condition the following linear approximation occurs
A i j k l l i n = 1 35 ( δ i j δ k l + δ i k δ j l + δ i l δ j k ) + 1 7 ( A i j δ k l + A i k δ j l + A i l δ j k + A k l δ i j + A j l δ i k + A j k δ i l ) .
The quadratic closure [8] omits all linear terms
A i j k l q u a d = A i j A k l .
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]
A i j k l h y b = ( 1 f ) A i j k l l i n + f A i j k l q u a d
f = 3 2 A i j A j i 1 2 .
Another approach to determine the function for the hybrid closure was determined by Advani and Tucker [9]
f = 1 27 · det A .
A more advanced approach is the natural closure (NAT). It is based on the relationship A = f ( A ) when no diffusion occurs ( C I = 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
A = S p p 4 π ( B : p p ) 3 2 d p
B = B · ( W + ξ D ) ( W + ξ D ) · B , B = I at t = 0 .
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
D A D t = C : D B D t , D B D t = D : D A D t .
hold true. The ordinary differential equations (ODEs) are solved simultaneously, which can be computed very efficiently. The special form used for this approach is
D A D t = C : F ( B ) + G ( A ) , D B D t = F ( B ) D : G ( A ) .
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
A ^ = λ 1 0 0 0 λ 2 0 0 0 λ 3 = R A T A R A ,
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
A ^ ̲ = A ̲ ^ 11 A ̲ ^ 12 A ̲ ^ 13 0 0 0 A ̲ ^ 12 A ̲ ^ 22 A ̲ ^ 23 0 0 0 A ̲ ^ 13 A ̲ ^ 23 A ̲ ^ 33 0 0 0 0 0 0 A ̲ ^ 44 0 0 0 0 0 0 A ̲ ^ 55 0 0 0 0 0 0 A ̲ ^ 66 .
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
A ̲ ^ 11 = f 1 ( λ 1 , λ 2 )
A ̲ ^ 22 = f 2 ( λ 1 , λ 2 )
A ̲ ^ 33 = f 3 ( λ 1 , λ 2 ) .
The remaining entries are defined by symmetry and normalization condition
A ̲ ^ 12 = A ̲ ^ 66
A ̲ ^ 23 = A ̲ ^ 44
A ̲ ^ 13 = A ̲ ^ 55
A ̲ ^ 55 + A ̲ ^ 66 = λ 1 A ̲ ^ 11
A ̲ ^ 44 + A ̲ ^ 66 = λ 2 A ̲ ^ 22
A ̲ ^ 44 + A ̲ ^ 55 = λ 3 A ̲ ^ 33 .
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
A ̲ ^ i i = f i ( λ 1 , λ 2 ) = C i 1 + C i 2 λ 1 + C i 3 λ 1 2 + C i 4 λ 2 + C i 5 λ 2 2 + C I 6 λ 1 λ 2
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 agreement with the distribution function and the NAT and better performance than the hybrid closure. For small C I the ORF closure shows oscillating behavior.
Chung and Kwon [24] improved the ORF closure to overcome the oscillation for small C I 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
A ̲ ^ i i = f i ( λ 1 , λ 2 ) = C i 1 + C i 2 λ 1 + C i 3 λ 1 2 + C i 4 λ 2 + C i 5 λ 2 2 + C I 6 λ 1 λ 2 + C i 7 λ 1 2 λ 2 + C i 8 λ 1 λ 2 2 + C i 9 λ 1 3 + C i 10 λ 2 3
for for i = 1 , 2 , 3 . The closures show stable behavior for a wide range of C I but tend to over predict the orientation.
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
A ̲ ^ 11 = f 1 ( λ 1 , λ 2 ) = λ 1 2
A ̲ ^ 22 = f 2 ( λ 1 , λ 2 ) = λ 2 2
A ̲ ^ 33 = f 3 ( λ 1 , λ 2 ) = ( 1 λ 1 λ 2 ) 2 .
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
A ̲ ^ i i = f i ( λ 1 , λ 2 ) = C i 1 + C i 2 λ 1 + C i 3 λ 1 2 + C i 4 λ 2 + C i 5 λ 2 2 + C I 6 λ 1 λ 2 + C i 7 λ 1 λ 2 ( 1 λ 1 λ 2 ) .
The data points U i , B i j , 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
A ̲ ^ i i = f i ( λ 1 , λ 2 ) = C i 1 + C i 2 λ 1 + C i 3 λ 1 2 + C i 4 λ 2 + C i 5 λ 2 2 + C I 6 λ 1 λ 2 .
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
  • fiber fiber interaction forces
  • elastic and bending forces (intra fiber forces)
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,34,35].
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 i j ( t ) = n = 1 N p n , j ( t ) p n , i ( t ) N ,
where N denotes the number of fibers, t the actual time, and p n ( 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].

6. Summary

The fiber orientation phenomenon can be simulated either on a macro- or micro-scale by using current numerical techniques and state-of-the-art modelling approaches. Macroscopic fiber orientation models are computationally efficient and are the preferred solution for estimating fiber orientation in large industrial applications. The limitation of those models is their phenomenological nature and the dependency on a set of fitting parameters. The challenge is typically the correct choice of parameters for a specific material. On the other hand, microscopic fiber models have a larger physical basis, but since fibers are explicitly discretized they are not suitable for large-scale simulations, because they are computational expensive. A combination of the models from the two scales can relieve the individual shortcomings and provides an interesting numerical solution for the virtual engineering design of fiber-reinforced polymer parts. One possibility is the prediction of optimal macroscopic parameters by a microscopic simulation, or the derivation of a macroscopic data-driven model using microscopic simulation results as input data.

Author Contributions

Conceptualization, S.K.K.; Writing—Original Draft Preparation, S.K.K.; Writing—Review and Editing, T.O., A.K. and C.C.; Visualization, S.K.K.; Supervision, T.O., A.K. and C.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
PDFprobability density function
ODEordinary differential equation
FTFolgar-Tucker
nemnematic
SRFstrain reduction factor
RSCreduced strain closure
RPRretarding principle rate
ARDanisotropic rotary diffusion
iARDimproved anisotropic rotary diffusion
pARDprincipal anisotropic rotary diffusion
MRDMoldflow rotary diffusion
NATnatural closure
FECfast exact closure
IBOFinvariant-based optimal fitted
EBOFeigenvalue-based optimal fitted
ORForthotropic fitted
OWEorthotropic fitted closure approximation for wide interaction coefficients
OWE3orthotropic fitted closure approximation for wide interaction
coefficients with third order polynomial approximation
SPHsmoothed particle hydrodynamic
DEMdiscrete element method
pFEAparticle finite element analysis
DNSdirect numerical simulation
MPSmoving particle semi-implicit
EBGelement bending group

References

  1. Folgar, F.; Tucker, C.L. Orientation Behavior of Fibers in Concentrated Suspensions. J. Reinf. Plast. Compos. 1984, 3, 98–119. [Google Scholar] [CrossRef]
  2. Latz, A.; Strautins, U.; Niedziela, D. Comparative numerical study of two concentrated fiber suspension models. J. Non-Newton. Fluid Mech. 2010, 165, 764–781. [Google Scholar] [CrossRef]
  3. Wang, J.; O’Gara, J.F.; Tucker, C.L. An objective model for slow orientation kinetics in concentrated fiber suspensions: Theory and rheological evidence. J. Rheol. 2008, 52, 1179–1200. [Google Scholar] [CrossRef]
  4. Tseng, H.C.; Chang, R.Y.; Hsu, C.H. Phenomenological improvements to predictive models of fiber orientation in concentrated suspensions. J. Rheol. 2013, 57, 1597–1631. [Google Scholar] [CrossRef]
  5. Phelps, J.H.; Tucker, C.L., III. An anisotropic rotary diffusion model for fiber orientation in short- and long-fiber thermoplastics. J. Non-Newton. Fluid Mech. 2009, 156, 165–176. [Google Scholar] [CrossRef]
  6. Tseng, H.C.; Chang, R.Y.; Hsu, C.H. The use of principal spatial tensor to predict anisotropic fiber orientation in concentrated fiber suspensions. J. Rheol. 2017, 62, 313–320. [Google Scholar] [CrossRef]
  7. Bakharev, A.; Yu, H.; Ray, S.; Speight, R.; Wang, J. Using New Anisotropic Rotational Diffusion Model to Improve Prediction of Short Fibers in Thermoplastic Injection Molding; ANTEC: Orlando, FL, USA, 2018. [Google Scholar]
  8. Advani, S.G.; Tucker, C.L. The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites. J. Rheol. 1987, 31, 751–784. [Google Scholar] [CrossRef]
  9. Advani, S.G.; Tucker, C.L. Closure approximations for three-dimensional structure tensors. J. Rheol. 1990, 34, 367–386. [Google Scholar] [CrossRef]
  10. Jeffery, G.B. The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 1922, 102, 161–179. [Google Scholar]
  11. Junk, M.; Illner, R. A New Derivation of Jeffery’s Equation. J. Math. Fluid Mech. 2007, 9, 455–488. [Google Scholar] [CrossRef]
  12. Huynh, H.M. Improved Fiber Orientation Prediction For Injection-Molded Composites. Master’s Thesis, University of Illinois Urbana-Champaign, Champaign County, IL, USA, 2001. [Google Scholar]
  13. Favaloro, A.J.; Tucker, C.L. Analysis of anisotropic rotary diffusion models for fiber orientation. Compos. Part A Appl. Sci. Manuf. 2019, 126, 105605. [Google Scholar] [CrossRef]
  14. Phan-Thien, N.; Fan, X.J.; Tanner, R.; Zheng, R. Folgar–Tucker constant for a fibre suspension in a Newtonian fluid. J. Non-Newton. Fluid Mech. 2002, 103, 251–260. [Google Scholar] [CrossRef]
  15. Bay, R.S. Fiber Orientation in Injection-Molded Composites: A Comparison of Theory and Experiment. Ph.D. Thesis, Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, Champaign County, IL, USA, 1991. [Google Scholar]
  16. Petrich, M.P.; Koch, D.L.; Cohen, C. An experimental determination of the stress–microstructure relationship in semi-concentrated fiber suspensions. J. Non-Newton. Fluid Mech. 2000, 95, 101–133. [Google Scholar] [CrossRef]
  17. Kuzmin, D. Planar and orthotropic closures for orientation tensors in fiber suspension flow models. SIAM J. Appl. Math. 2018, 78, 3040–3059. [Google Scholar] [CrossRef] [Green Version]
  18. Verley, V.; Dupret, F. Numerical prediction of the fiber orientation in complex injection molded parts. Trans. Eng. Sci. 1994. [Google Scholar] [CrossRef]
  19. Montgomery-Smith, S.; He, W.; Jack, D.; Smith, D. Exact tensor closures for the three-dimensional Jeffery’s equation. J. Fluid Mech. 2011, 680, 321–335. [Google Scholar] [CrossRef] [Green Version]
  20. Dinh, S.M.; Armstrong, R.C. A Rheological Equation of State for Semiconcentrated Fiber Suspensions. J. Rheol. 1984, 28, 207–227. [Google Scholar] [CrossRef]
  21. Montgomery-Smith, S.; Jack, D.; Smith, D.E. The Fast Exact Closure for Jeffery’s equation with diffusion. J. Non-Newton. Fluid Mech. 2011, 166, 343–353. [Google Scholar] [CrossRef]
  22. Chung, D.H.; Kwon, T.H. Invariant-based optimal fitting closure approximation for the numerical prediction of flow-induced fiber orientation. J. Rheol. 2002, 46, 169–194. [Google Scholar] [CrossRef] [Green Version]
  23. Cintra, J.S.; Tucker, C.L. Orthotropic closure approximations for flow-induced fiber orientation. J. Rheol. 1995, 39, 1095–1122. [Google Scholar] [CrossRef]
  24. Chung, D.H.; Kwon, T.H. Improved model of orthotropic closure approximation for flow induced fiber orientation. Polym. Compos. 2001, 22, 636–649. [Google Scholar] [CrossRef]
  25. Yashiro, S.; Okabe, T.; Matsushima, K. A Numerical Approach for Injection Molding of Short-Fiber-Reinforced Plastics Using a Particle Method. Adv. Compos. Mater. 2011, 20, 503–517. [Google Scholar] [CrossRef]
  26. Yashiro, S.; Sasaki, H.; Sakaida, Y. Particle simulation for predicting fiber motion in injection molding of short-fiber-reinforced composites. Compos. Part A Appl. Sci. Manuf. 2012, 43, 1754–1764. [Google Scholar] [CrossRef]
  27. He, L.; Lu, G.; Chen, D.; Li, W.; Lu, C. Three-dimensional smoothed particle hydrodynamics simulation for injection molding flow of short fiber-reinforced polymer composites. Model. Simul. Mater. Sci. Eng. 2017, 25, 055007. [Google Scholar] [CrossRef]
  28. Bertevas, E.; Férec, J.; Khoo, B.C.; Ausias, G.; Phan-Thien, N. Smoothed particle hydrodynamics (SPH) modeling of fiber orientation in a 3D printing process. Phys. Fluids 2018, 30, 103103. [Google Scholar] [CrossRef]
  29. Yamagata, N.; Ichimiya, M. Numerical Approach of Viscous Flow Containing Short Fiber by SPH Method. In Computational and Experimental Simulations in Engineering; Okada, H., Atluri, S.N., Eds.; Springer International Publishing: Cham, Switzerland, 2020; Volume 75, pp. 301–307. [Google Scholar] [CrossRef]
  30. Wu, K.; Wan, L.; Zhang, H.; Yang, D. Numerical simulation of the injection molding process of short fiber composites by an integrated particle approach. Int. J. Adv. Manuf. Technol. 2018, 97, 3479–3491. [Google Scholar] [CrossRef]
  31. Yang, X.; Liu, M.; Peng, S. Smoothed particle hydrodynamics and element bending group modeling of flexible fibers interacting with viscous fluids. Phys. Rev. E 2014, 90. [Google Scholar] [CrossRef] [Green Version]
  32. Yang, X.; Liu, M.B. Bending modes and transition criteria for a flexible fiber in viscous flows. J. Hydrodyn. 2016, 28, 1043–1048. [Google Scholar] [CrossRef]
  33. Challabotla, N.R.; Zhao, L.; Andersson, H.I. On fiber behavior in turbulent vertical channel flow. Chem. Eng. Sci. 2016, 153, 75–86. [Google Scholar] [CrossRef]
  34. Dotto, D.; Marchioli, C. Orientation, distribution, and deformation of inertial flexible fibers in turbulent channel flow. Acta Mech. 2019, 230, 597–621. [Google Scholar] [CrossRef]
  35. Njobuenwu, D.O.; Fairweather, M. Simulation of inertial fibre orientation in turbulent flow. Phys. Fluids 2016, 28, 063307. [Google Scholar] [CrossRef] [Green Version]
  36. Yamamoto, S.; Matsuoka, T. A method for dynamic simulation of rigid and flexible fibers in a flow field. J. Chem. Phys. 1993, 98, 644–650. [Google Scholar] [CrossRef]
  37. Yamamoto, S.; Matsuoka, T. Viscosity of dilute suspensions of rodlike particles: A numerical simulation method. J. Chem. Phys. 1994, 100, 3317–3324. [Google Scholar] [CrossRef]
  38. Yamane, Y.; Kaneda, Y.; Dio, M. Numerical simulation of semi-dilute suspensions of rodlike particles in shear flow. J. Non-Newton. Fluid Mech. 1994, 54, 405–421. [Google Scholar] [CrossRef]
  39. Yamane, Y.; Kaneda, Y.; Doi, M. The Effect of Interaction of Rodlike Particles in Semi-Dilute Suspensions under Shear Flow. J. Phys. Soc. Jpn 1995, 64, 3265–3274. [Google Scholar] [CrossRef]
  40. Yamamoto, S.; Matsuoka, T. Dynamic simulation of fiber suspensions in shear flow. J. Chem. Phys. 1995, 102, 2254–2260. [Google Scholar] [CrossRef]
  41. Thomasset, J.; Grmela, M.; Carreau, P.J. Microstructure and rheology of polymer melts reinforced by long glass fibres: Direct simulations. J. Non-Newton. Fluid Mech. 1997, 73, 195–203. [Google Scholar] [CrossRef]
  42. Sundararajakumar, R.R.; Koch, D.L. Structure and properties of sheared fiber suspensions with mechanical contacts. J. Non-Newton. Fluid Mech. 1997, 73, 205–239. [Google Scholar] [CrossRef]
  43. Skjetne, P.; Ross, R.F.; Klingenberg, D.J. Simulation of single fiber dynamics. J. Chem. Phys. 1997, 107, 2108–2121. [Google Scholar] [CrossRef]
  44. Ross, R.F.; Klingenberg, D.J. Dynamic simulation of flexible fibers composed of linked rigid bodies. J. Chem. Phys. 1997, 106, 2949–2960. [Google Scholar] [CrossRef] [Green Version]
  45. Fan, X.; Phan-Thien, N.; Zheng, R. A direct simulation of fibre suspensions. J. Non-Newton. Fluid Mech. 1998, 74, 113–135. [Google Scholar] [CrossRef]
  46. Harlen, O.G.; Sundararajakumar, R.R.; Koch, D.L. Numerical simulations of a sphere settling through a suspension of neutrally buoyant fibres. J. Fluid Mech. 1999, 388, 355–388. [Google Scholar] [CrossRef]
  47. Joung, C.G.; Phan-Thien, N.; Fan, X.J. Direct simulation of flexible fibers. J. Non-Newton. Fluid Mech. 2001, 99, 1–36. [Google Scholar] [CrossRef]
  48. Joung, C.G.; Phan-Thien, N.; Fan, X.J. Viscosity of curved fibers in suspension. J. Non-Newton. Fluid Mech. 2002, 102, 1–17. [Google Scholar] [CrossRef]
  49. Joung, C.G. Direct Simulation Studies of Suspended Particles and Fibre-Filled Suspensions. Ph.D. Thesis, School of Aerospace, Mechanical and Mechatronic Engineering The University of Sydney, Sydney, Australia, 2003. [Google Scholar]
  50. Switzer, L.H.; Klingenberg, D.J. Rheology of sheared flexible fiber suspensions via fiber-level simulations. J. Rheol. 2003, 47, 759–778. [Google Scholar] [CrossRef] [Green Version]
  51. Kromkamp, J.; Van Den Ende, D.T.M.; Kandhai, D.; Van Der Sman, R.G.M.; Boom, R.M. Shear-induced self-diffusion and microstructure in non-Brownian suspensions at non-zero Reynolds numbers. J. Fluid Mech. 2005, 529, 253–278. [Google Scholar] [CrossRef] [Green Version]
  52. Ausias, G.; Fan, X.; Tanner, R. Direct simulation for concentrated fibre suspensions in transient and steady state shear flows. J. Non-Newton. Fluid Mech. 2006, 135, 46–57. [Google Scholar] [CrossRef]
  53. Wang, G.; Yu, W.; Zhou, C. Optimization of the rod chain model to simulate the motions of a long flexible fiber in simple shear flows. Eur. J. Mech.-B/Fluids 2006, 25, 337–347. [Google Scholar] [CrossRef]
  54. Lindström, S.B.; Uesaka, T. Simulation of the motion of flexible fibers in viscous fluid flow. Phys. Fluids 2007, 19, 113307. [Google Scholar] [CrossRef]
  55. Lindström, S.B.; Uesaka, T. Simulation of semidilute suspensions of non-Brownian fibers in shear flow. J. Chem. Phys. 2008, 128, 024901. [Google Scholar] [CrossRef] [PubMed]
  56. Lindström, S.B.; Uesaka, T. A numerical investigation of the rheology of sheared fiber suspensions. Phys. Fluids 2009, 21, 083301. [Google Scholar] [CrossRef]
  57. Yamanoi, M.; Maia, J.M. Analysis of rheological properties of fibre suspensions in a Newtonian fluid by direct fibre simulation. Part1: Rigid fibre suspensions. J. Non-Newton. Fluid Mech. 2010, 165, 1055–1063. [Google Scholar] [CrossRef]
  58. Yamanoi, M.; Maia, J.; Kwak, T.s. Analysis of rheological properties of fibre suspensions in a Newtonian fluid by direct fibre simulation. Part 2: Flexible fibre suspensions. J. Non-Newton. Fluid Mech. 2010, 165, 1064–1071. [Google Scholar] [CrossRef]
  59. Yamanoi, M.; Maia, J.M. Analysis of rheological properties of fiber suspensions in a Newtonian fluid by direct fiber simulation. Part 3: Behavior in uniaxial extensional flows. J. Non-Newton. Fluid Mech. 2010, 165, 1682–1687. [Google Scholar] [CrossRef]
  60. Yamanoi, M.; Maia, J.M. Stokesian dynamics simulation of the role of hydrodynamic interactions on the behavior of a single particle suspending in a Newtonian fluid. Part 1. 1D flexible and rigid fibers. J. Non-Newton. Fluid Mech. 2011, 166, 457–468. [Google Scholar] [CrossRef]
  61. Andrić, J.; Fredriksson, S.T.; Lindström, S.B.; Sasic, S.; Nilsson, H. A study of a flexible fiber model and its behavior in DNS of turbulent channel flow. Acta Mech. 2013, 224, 2359–2374. [Google Scholar] [CrossRef] [Green Version]
  62. Andrić, J.; Lindström, S.; Sasic, S.; Nilsson, H. Rheological properties of dilute suspensions of rigid and flexible fibers. J. Non-Newton. Fluid Mech. 2014, 212, 36–46. [Google Scholar] [CrossRef]
  63. Do-Quang, M.; Amberg, G.; Brethouwer, G.; Johansson, A.V. Simulation of finite-size fibers in turbulent channel flows. Phys. Rev. E 2014, 89. [Google Scholar] [CrossRef]
  64. Mezher, R.; Abisset-Chavanne, E.; Férec, J.; Ausias, G.; Chinesta, F. Direct simulation of concentrated fiber suspensions subjected to bending effects. Model. Simul. Mater. Sci. Eng. 2015, 23, 055007. [Google Scholar] [CrossRef]
  65. Mezher, R.; Perez, M.; Scheuer, A.; Abisset-Chavanne, E.; Chinesta, F.; Keunings, R. Analysis of the Folgar & Tucker model for concentrated fibre suspensions in unconfined and confined shear flows via direct numerical simulation. Compos. Part A Appl. Sci. Manuf. 2016, 91, 388–397. [Google Scholar] [CrossRef]
  66. Wang, J.; Cook, P.; Bakharev, A.; Costa, F.; Astbury, D. Prediction of fiber orientation in injection-molded parts using three-dimensional simulations. AIP Conf. Proc. 2016, 1713, 040007. [Google Scholar] [CrossRef] [Green Version]
  67. Perez, M.; Scheuer, A.; Abisset-Chavanne, E.; Chinesta, F.; Keunings, R. A multi-scale description of orientation in simple shear flows of confined rod suspensions. J. Non-Newton. Fluid Mech. 2016, 233, 61–74. [Google Scholar] [CrossRef] [Green Version]
  68. Sasayama, T.; Inagaki, M. Simplified bead-chain model for direct fiber simulation in viscous flow. J. Non-Newton. Fluid Mech. 2017, 250, 52–58. [Google Scholar] [CrossRef]
  69. Kuhn, C.; Walter, I.; Taeger, O.; Osswald, T. Experimental and Numerical Analysis of Fiber Matrix Separation during Compression Molding of Long Fiber Reinforced Thermoplastics. J. Compos. Sci. 2017, 1, 2. [Google Scholar] [CrossRef] [Green Version]
  70. Kuhn, C. Analysis and Prediction of Fiber Matrix Separation during Compression Molding of Fiber Reinforced Plastics. Ph.D. Thesis, Friedrich-Alexander Universität Erlangen-Nürnberg, Erlangen, Germany, 2018. [Google Scholar]
  71. Meirson, G.; Hrymak, A.N. Two dimensional long-flexible fiber orientation simulation in squeeze flow. Polym. Compos. 2018, 39, 4656–4665. [Google Scholar] [CrossRef]
  72. Sasayama, T.; Inagaki, M. Efficient bead-chain model for predicting fiber motion during molding of fiber-reinforced thermoplastics. J. Non-Newton. Fluid Mech. 2019, 264, 135–143. [Google Scholar] [CrossRef]
  73. Sasayama, T.; Inagaki, M.; Sato, N. Direct simulation of glass fiber breakage in simple shear flow considering fiber-fiber interaction. Compos. Part A Appl. Sci. Manuf. 2019, 124, 105514. [Google Scholar] [CrossRef]
  74. Laurencin, T.; Laure, P.; Orgéas, L.; Dumont, P.; Silva, L.; Rolland du Roscoat, S. Fibre kinematics in dilute non-Newtonian fibre suspensions during confined and lubricated squeeze flow: Direct numerical simulation and analytical modelling. J. Non-Newton. Fluid Mech. 2019, 273, 104187. [Google Scholar] [CrossRef]
  75. Pérez, C. The Use of a Direct Particle Simulation to Predict Fiber Motion in Polymer Processing. Ph.D. Thesis, University of Wisconsin-Madison, Madison, WI, USA, 2017. [Google Scholar]
  76. Yun, M.; Argerich Martin, C.; Giormini, P.; Chinesta, F.; Advani, S. Learning the Macroscopic Flow Model of Short Fiber Suspensions from Fine-Scale Simulated Data. Entropy 2019, 22, 30. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Orientation of a single rigid fiber p .
Figure 1. Orientation of a single rigid fiber p .
Jcs 04 00069 g001
Figure 2. Influence of the phenomenological parameter C I on fiber orientation evolution.
Figure 2. Influence of the phenomenological parameter C I on fiber orientation evolution.
Jcs 04 00069 g002
Figure 3. Influence of the phenomenological parameter κ on fiber orientation evolution. The diffusion constant is set to C I = 0.01 .
Figure 3. Influence of the phenomenological parameter κ on fiber orientation evolution. The diffusion constant is set to C I = 0.01 .
Jcs 04 00069 g003
Figure 4. Influence of the phenomenological parameter D 2 on fiber orientation evolution. The diffusion is et to C I = 0.01 .
Figure 4. Influence of the phenomenological parameter D 2 on fiber orientation evolution. The diffusion is et to C I = 0.01 .
Jcs 04 00069 g004
Figure 5. Reference triangle for orthotropic closure.
Figure 5. Reference triangle for orthotropic closure.
Jcs 04 00069 g005
Table 1. Literature overview of element based simulations for fiber reinforced polymers in chronological order.
Table 1. Literature overview of element based simulations for fiber reinforced polymers in chronological order.
Discretization of FibersFlow FieldsFluid-Fiber InteractionFiber-Fiber InteractionFlexibilityRegarded Quantities
Yamamoto and Matsuoka 1993 [36]chain of beadsshearone-way coupled-flexiblesingle fiber movement
Yamamoto and Matsuoka 1994 [37]chain of beadsshearone-way coupled-flexibleviscosity of dilute solutions
Yamane et al., 1994 [38]rodsshearone-way coupledlubrication-semi dilute suspensions, orientation evolution, diffusion constant, shear viscosity
Yamane et al., 1995 [39]rodsshearone-way coupledlubrication-semi dilute suspensions, bounded and unbounded system
Yamamoto and Matsuoka 1995 [40]chain of beadsshearone-way coupledlubricationflexibleconcentrated suspension, viscosity, stresses
Thomasset et al., 1997 [41]rigid rodsvaries flow fieldsone-way coupledlubrication, mechanical and hydrodynamical contact, no friction-2D, effects of fiber motion and orientation
Sundararajakumar and Koch 1997 [42]rodsshearone-way coupledlubrication, mechanical and hydrodynamical contact, no friction-dilute: hydrodynamical contact most important, approaching higher concentration fiber contact
Skjetne et al., 1997 [43]prolate spheroidsshearone-way coupled-flexible, rigidsingle fiber movement
Ross and Klingenberg 1997 [44]prolate spheroidsshearone-way coupledrepulsive interactionsflexible and rigidsingle fiber movement, viscosity
Fan et al., 1998 [45]rodsshearone-way coupledlubrication, no friction, long range hydrodynamic interactions by slender body theory-orientation, viscosity, stresses, all regimes
Harlen et al., 1999 [46]rodsno imposed flow-mechanical contact, friction, long range hydrodynamic interactions by slender body theory-sphere settling through suspension of neutrally buoyant fibers, fiber contact has significant influence
Phan-Thien et al., [14]rodsshearone-way coupledlubrication, no friction, long range hydrodynamic interactions by slender body theory-FT constant, dilute and semi dilute
Joung et al., 2001 [47]chain of spherical beadsshear and extensional flowsone-way coupledlubrication, preventing from overlapping, long range hydrodynamic interactionsrigid, flexibleviscosity, orientation
Joung et al., 2002 [48]chain of spherical beadsshear and complex flowsone-way coupledlubrication, preventing from overlapping, long range hydrodynamic interactionsrigid, curvedviscosity for curved fibers
Joung et al., 2003 [49]chain of spherical beadsshear and complex flowsone-way coupledlubrication, preventing from overlapping, long range hydrodynamic interactionsflexibleJeffrey orbits for rigid and flexible fibers, relationship between fiber stiffness, and bulk viscosity, arbitrary particle shapes, dilute regime
Switzer and Klingenberg 2003 [50]chain of rodsshearone-way coupledmechanical interaction, frictionflexibleeffects of shape, friction, aspect ratio and stiffness, yield stress, rheology in flocculated systems
Kromkamp et al., 2005 [51]rodsshearcoupled, Lattice Bolzmann for fluid forces, particles as boundary surfaceslubrication correction, mechanical interaction, no friction-2D, effects of shear rate on flow behavior and micro structure, shear-induced self diffusion
Ausias et al., 2006 [52]rigid prolate spheroidsshearone-way coupledlubrication and interaction in normal direction, no friction, no long range interactions-orientation, viscosity, stresses, up to Φ = 11.5
Wang et al., 2006 [53]rod-chainshearone-way coupled flexible, rigidoptimal rod length for high accuracy and efficient calculation
Lindström and Uesaka 2007 [54]rod-chain modelshearcoarse two-way couplinglubrication and interaction in normal direction, frictionflexibleJeffrey orbits, curvature, regimes of motions for flexible fibers
Lindström and Uesaka 2008 [55]rod-chain modelshearcoarse two-way couplinglubrication and interaction in normal direction, frictionflexibleorientation, viscosity, dilute and semidilute regime,
Lindström and Uesaka 2009 [56]rod-chain modelshearcoarse two-way couplinglubrication and interaction in normal direction, frictionflexiblerheological properties
Yamanoi and Maia 2010 [57]chain of beadsshearone-way couplinglubrication, mechanical contact, long range hydrodynamic interactions-rheological properties, orientation
Yamanoi et al., 2010 [58]chain of beadsshearone-way couplinglubrication, mechanical contact, long range hydrodynamic interactionsflexiblenylon fiber, rheological properties, orientation, effect of flexibility
Yamanoi and Maia 2010 [59]chain of beadsuniaxel elongation flowone-way couplinglubrication, mechanical contact, long range hydrodynamic interactionsflexiblerheological properties, orientation, orientation tensor independent of aspect ratio, volume fraction
Yamanoi and Maia 2011 [60]chain of beadssheartwo way couplingsingle fiberrigid and flexiblehydrodynamic interaction in single fiber movement in shear
Andrić et al., 2013 [61]rod-chain modelturbulent flowtwo way coupling, DNS for fluid motionsingle fiberrigid and flexiblefiber-flow interaction for a single fiber
Andrić et al., 2014 [62]rod-chain modelsheartwo way coupling, DNS for fluid motion-rigid and flexibledilute solution, no interaction, rheological properties, orbit drifts
Do-Quang et al., 2014 [63]rod-chain modelturbulent flowtwo way coupling, entropy lattice Boltzmann for fluid, external boundary force methodlubrication and mechanical contactrigidcellulose fibers in water, accumulation effects
Mezher et al., 2015 [64]prolate spheroidsshearone-way couplinglubrication and interaction in normal direction, no friction, no long range hydrodynamicsflexibleconcentrated, orientation, normalized stresses, interactions, elastic energy
Mezher et al., 2016 [65]prolate spheroidsshearone-way couplinglubrication and interaction in normal direction, no friction, no long range hydrodynamicsflexibleconcentrated Φ = 7 18.2 , orientation, diffusion constants, confinement effects
Wang et al., 2016 [66]rod-chain modelshearone-way coupled-flexiblenew rod chain model, optimal rod length
Perez et al., 2016 [67]rodshearone-way couplingonly wall interaction-dilute, confinement effects
Sasayama and Inagaki 2017 [68]simplified bead-chain modelshearone-waymechanical, lubricationflexiblesimplified bead-chain model for hydrodynamic calculations
Kuhn et al., 2017 [69]rod-chain modelcomplex flow fieldsone-way couplingmechanical, friction, no long-range hydrodynamicflexiblefiber matrix separation in compression molding, LFRT
Kuhn et al., 2018 [70]rod-chain modelcomplex flow fieldsone-way couplingmechanical, friction, no long-range hydrodynamicflexiblerib filling
Meirson and Hrymak 2018 [71]rod-chain modelsqueezeone-way-flexible2D, fiber orientation and deformation
Wu et al., 2018 [30]bounded spherescomplex flow2-way coupled, SPH for fluid motionlinear contact-2D, fiber orientation, accumulation during injection molding
Sasayama and Inagaki 2019 [72]efficient bead-chain modelshearone-waymechanical, lubrication, frictionflexibleefficient bead-chain model for hydrodynamic calculations
Sasayama et al., 2019 [73]efficient bead-chain modeshearone-waymechanical, lubrication, frictionflexible, breakagefiber breakage
Laurentcin et al., 2019 [74]sphero- cylindersqueeze flow (lubricated, non-lubricated)one-way-rigidnon-Newtonian fluid, dilute regime, comparison between numerical analytical and experimental results

Share and Cite

MDPI and ACS Style

Kugler, S.K.; Kech, A.; Cruz, C.; Osswald, T. Fiber Orientation Predictions—A Review of Existing Models. J. Compos. Sci. 2020, 4, 69. https://doi.org/10.3390/jcs4020069

AMA Style

Kugler SK, Kech A, Cruz C, Osswald T. Fiber Orientation Predictions—A Review of Existing Models. Journal of Composites Science. 2020; 4(2):69. https://doi.org/10.3390/jcs4020069

Chicago/Turabian Style

Kugler, Susanne Katrin, Armin Kech, Camilo Cruz, and Tim Osswald. 2020. "Fiber Orientation Predictions—A Review of Existing Models" Journal of Composites Science 4, no. 2: 69. https://doi.org/10.3390/jcs4020069

APA Style

Kugler, S. K., Kech, A., Cruz, C., & Osswald, T. (2020). Fiber Orientation Predictions—A Review of Existing Models. Journal of Composites Science, 4(2), 69. https://doi.org/10.3390/jcs4020069

Article Metrics

Back to TopTop