Microstructure generation of long fiber reinforced hybrid composites using the fused sequential addition and migration method

Abstract

We present a fused sequential addition and migration (fSAM) algorithm to generate synthetic microstructures for long fiber reinforced hybrid composites. To incorporate the bending of long fibers, we model the fibers as polygonal chains and use an optimization framework to search for a non-overlapping fiber configuration with the desired properties. As microstructures of hybrid composites consist of materials with varying characteristics, e.g., different fiber orientation distributions or fiber volume fractions, the unit cell needs to be divided in subcells. Then, the objective function is required to account for the characteristics of the subcells individually. To ensure that the location of a fiber is restricted to its respective subcell, we apply box constraints during the optimization procedure. Depending on the materials and the manufacturing processes, the interfaces within hybrid composites show either strict separation or interlocking of the different materials. Hence, we enable selecting the box-constraint type to model the interface according to the desired composite. We provide a detailed discussion of the extensions of the original fSAM algorithm to generate hybrid composites, considering two alternative procedures to handle the box constraints during the optimization procedure. As the selection of the box-constraint type influences the synthetic microstructures, we investigate the shape of the interfacial areas and the computed effective stiffnesses. Last but not least, to validate the capability of the fSAM algorithm to generate representative microstructures, we compare the resulting mechanical behavior with experimental data.

Keywords

Long fiber reinforced hybrid composites curved fibers microstructure generation representative volume elements fused sequential addition and migration

Introduction

State of the art

Driven by the need for materials which feature high specific stiffnesses and permit to be mass-produced, fiber-reinforced composites play an essential role for lightweight technologies.^1–3 Of central importance are fiber-reinforced hybrid composites, which consist of different fiber-reinforced material systems, as they enable tailoring the material to specific requirements. In the field of hybrid composites, the continuous-discontinuous fiber-reinforced polymers (CoDicoFRP) represent an innovative material class which combines the high load-bearing capacity of the continuous (Co) reinforcement with the design freedom of the discontinuous (Dico) reinforcement.⁴ To use CoDicoFRP for lightweight applications, tools to predict the mechanical properties are necessary which account for the heterogeneity, randomness and anisotropy of the material system. Characterizing the effective properties by experiments requires a significant expense in time and resources, especially to cover the various fiber orientation states and fiber volume fractions. As an efficient alternative, computational multiscale methods^5–7 based on the theory of homogenization^8,9 may be used. However, such methods require a geometrical description of the underlying microstructure to be given, for a start.

For fiber-reinforced composites, microstructural images may be obtained from micro-CT scanning.^10–13 However, using these real digital images as a starting point for computational multiscale methods suffers from several disadvantages. First, the expense to obtain a single image is rather high. Secondly, each image represents a specific combination of microstructural descriptors, e.g., the fiber orientation state and the fiber volume fraction. For desired descriptors, a corresponding microstructure needs to be manufactured, which may be difficult to realize with desired accuracy. Last but not least, the extracted geometries are not periodic, which leads to boundary artifacts and decreased representativity of the unit cells.^14–16 Due to these disadvantages, microstructure generation tools are commonly used to complement the real digital images with synthetic microstructures for various microstructural descriptors and periodic geometries.

Microstructure generation tools for particle composites may be classified in two categories. On the one hand, sequential insertion algorithms place the particles consecutively, keeping the location of the particles fixed after insertion. The most widespread algorithm in this class is the random sequential addition (RSA) algorithm,^17,18 placing the particles sequentially within the unit cell under the condition that a newly inserted particle does not overlap with previously placed particles. The packing is considered successful if the cell is packed to the desired fiber volume fraction. For fiber-reinforced composites with cylindrical fibers, the direction, the length and the midpoint of each fiber are sampled from pre-defined distribution functions. Based on the RSA algorithm, various extensions are available.^19–24 However, in the context of long fiber reinforced composites, sequential insertion algorithms fail to achieve industrial fiber volume fractions for general orientation states. Besides the RSA algorithm, another representative for sequential insertion algorithms is the sequential deposition algorithm,^25,26 where the fibers descend into the unit cell driven by Newton’s law with gravitational force. In contrast to the RSA algorithm, fiber bending is typically considered. As a result, larger fiber volume fractions may be realized by such algorithms. However, the sequential deposition algorithm is restricted to rather planar fiber orientation states.

On the other hand, collective rearrangement algorithms position all particles in the cell first and subsequently move the fillers simultaneously.^27–29 The mechanical contraction method (MCM)³⁰ generates isotropic microstructures with spherocylinders, i.e., cylinders with half-caps at their ends, as fillers. Therefore, first a microstructure with small fiber volume fraction is prepared by the RSA. Subsequently, the cell dimensions are reduced without changing the size of the fillers. The resulting collisions are removed with an optimization procedure. Based on the MCM, the sequential addition and migration (SAM) algorithm³¹ for short fibers with uniform length is capable of realizing more general fiber orientation states. First, a fiber configuration is sampled according to the desired characteristics, i.e., the fiber volume fraction and the fiber orientation distribution, where the fibers may overlap. Then, an optimization framework is used to find roots of a non-negative objective function which encodes the non-overlap condition and may account for additional criteria, e.g., the fiber orientation state. Further extensions of the SAM algorithm account for fiber length distributions,³² typical for discontinuous fiber-reinforced composites,^33–36 and coupling between the fiber orientation and the fiber length data.³⁷

Whereas representing short fibers by straight cylinders is a realistic modeling assumption, the bending of fibers plays a key role for long fiber reinforced composites. Additionally, tremendous unit cell sizes are necessary to realize long fiber reinforced composites with straight fibers as the unit cell edge-lengths are required to be at least as large as the mean fiber length. To account for the fiber bending and decrease the unit cell size, the SAM algorithm for long fibers³⁸ models a fiber as a polygonal chain with spherocylinders as segments. During the optimization problem, the segments are moved separately and the connectivity of adjacent segments is ensured at convergence by an additional penalty term in the objective function. However, this procedure leads to convergence problems for complex microstructures with, e.g., larger fiber aspect ratios at higher fiber volume fractions. To overcome this limitation, the fused sequential addition and migration (fSAM) algorithm³⁹ represents an alternative approach to the SAM algorithm for long fibers. The key aspect of the fSAM algorithm is an iterative optimization scheme which accounts for a fused fiber movement, i.e., the polygonal chain is connected for each iterate of the algorithm. Thus, the optimization space is described by the curved, i.e., non-flat, manifold of fused polygonal chains. As the standard gradient descent approach is only admissible for flat optimization spaces, the fSAM algorithm is based on an adapted gradient descent approach moving along the geodesics,^40,41 i.e., the shortest path between two points on a manifold. Thus, the connectivity constraint is fulfilled in every iterative step and an additional penalty term responsible for keeping the chains connected within the objective function is not necessary. Due to the physically meaningful formulation of the fiber movement, the fSAM algorithm proves to generate microstructures for industrial long fiber reinforced composites in a computationally efficient manner.

Contributions

The microstructure of hybrid fiber-reinforced composites consists of subcells with different characteristics. However, current microstructure generation tools prescribe the desired quantities, like the fiber volume fraction or the fiber orientation distribution, for the entire microstructure without distinguishing between specific areas. Hence, these tools are only applicable to composites with location independent descriptors, e.g., to pure discontinuous fiber-reinforced composites with a single fiber volume fraction, fiber length and orientation distribution, but not to materials with location-dependent descriptors. In this manuscript, we address the challenge to generate microstructures for hybrid fiber-reinforced composites with continuous and discontinuous long fiber reinforcement, where we target unidirectional fiber orientation states for the continuous fiber-reinforcement. For the microstructure generation tool, we define two additional requirements. First, we aim to minimize the necessary unit cell size for representativity. Therefore, we target microstructures with periodic geometry and accurately fulfilled descriptive components, i.e., the fiber volume fraction or the fiber orientation tensor^42,43 of desired order. Modeling the fibers with straight geometries leads to unit cell edge-lengths which need to be at least as large as the mean fiber length. As this may be restrictive for long fiber reinforced composites, we account for the fiber bending of long fibers, which permits to decrease the unit cell size compared to straight fibers. As a second task, we need to model the areas between adjacent materials. Depending on the materials and the manufacturing processes, these interfaces may feature a strict separation of the areas or an interlocking where fibers reach into the adjacent layers.

For this purpose, we extend the fused sequential addition and migration algorithm³⁹ to generate hybrid long fiber reinforced composites. The classical fSAM algorithm³⁹ models the fibers as polygonal chains and positions the fibers within a periodic unit cell without further restrictions on the fibers’ location. To distinguish between different subcells of a hybrid composite, we assign each fiber to a specific subcell and restrict the location of the fiber by enforcing box constraints. In this manuscript, we introduce two types of box constraints. On the one hand, hard box constraints require that a fiber is completely located within its subcell, which leads to a strict separation between adjacent subcells. On the other hand, soft box constraints enforce that only the midpoints of the segments are within the respective subcell, which enables an interlocking between adjacent subcells. For soft box constraints, the degree of interlocking decreases with decreasing segment lengths. To ensure that interlocking is still possible for small segment lengths, we additionally consider an adapted soft constraint type where only specific midpoints of the segments are restricted.

The original fSAM algorithm³⁹ solves an optimization problem to find an admissible fiber configuration where the fibers are in a non-overlapping configuration and further criteria are fulfilled. To ensure that the iterates of the gradient descent are admissible points on the curved, i.e., non-flat, configuration space of polygonal chains, an adapted iteration rule is used. The key aspect of this iteration rule is a movement along the geodesics, i.e., the shortest path between two points on the configuration space. The equations governing the geodesics may be formulated as mechanical problems with holonomic constraints and no external forces in the form of a d’Alembert type constrained mechanical system. Thus, a movement along the geodesic is obtained by integrating the d’Alembert type constrained mechanical system.

During the optimization procedure of the fSAM algorithm, the box constraints may be considered in two ways. On the one hand, only the solution configuration at convergence fulfills the box constraints, which may be realized by adding a penalty term to the objective function. Hence, we call this procedure the penalty approach. On the other hand, the box constraints may be fulfilled by using an iteration rule which computes admissible iterates on the configuration space of constrained polygonal chains. For the latter case, the box constraints are enforced via an intrinsic approach, and a penalty term is not necessary. Due to the box constraints, the configuration space of a constrained polygonal chain is a manifold with boundary. However, the original iteration rule of the fSAM algorithm³⁹ is only applicable to smooth manifolds including equality constraints in their description. To obtain an iteration rule to move along the geodesics respecting the box constraints, we extend the integration scheme of the fSAM algorithm for manifolds with boundary. Therefore, we formulate the box constraints as inequality constraints for the numerical integration of the d’Alembert type constrained mechanical system, using the Karush-Kuhn-Tucker (KKT) conditions. Then, we solve the discretized system of equations with a non-smooth Newton scheme^44–46 in every time step. For the penalty and the intrinsic approach, we discuss the integration into the framework of the fSAM algorithm. Additionally, we provide information on further extensions of the fSAM algorithm to generate microstructures of hybrid composites, e.g., the adaption of the objective function to account for subcells with varying desired characteristics.

We start the computational investigations by studying the performance of the microstructure generator handling the box constraints with either the penalty or the intrinsic approach. Subsequently, we investigate the necessary resolution and RVE size for a two-layered CoDicoFRP. Furthermore, we compare the interface between adjacent subcells as well as the effective properties of the entire microstructure when using the hard or soft box-constraint type. Last but not least, we validate the computed effective elastic properties of a CoDicoFRP with discontinuous core and continuous shell layers by comparing the stiffness with experimentally measured data.

Notation

This manuscript uses a direct tensor notation or matrix-vector notation with orthonormal bases ${e_{1}, \dots, e_{n}} (n \in N)$ . Scalars are indicated by non-bold letters, e.g. b, vectors in matrix-vector notation use non-cursive bold lowercase letters, e.g. b, and vectors in direct tensor notation use cursive bold lowercase letters, e.g. b . Matrices are represented by bold non-cursive uppercase letters, e.g. B, second-order tensors by bold cursive uppercase letters, e.g. B , and fourth order tensors are denoted by $B$ . For example, with respect to a three-dimensional vector space, a second-order tensor in diagonalized form is indicated by, e.g., $B \hat{=} d i a g (b_{1}, b_{2}, b_{3})$ . The following mathematical operations are used for the direct tensor notation as well as for the matrix-vector notation. We denote a transposed second-order tensor by, e.g., B ^T, a linear mapping of a first-order tensor to a second-order tensor by, e.g., a = Cb and the linear mapping with full contraction involving higher-order tensors by, e.g., $a = B [C_{⟨ 3 ⟩}] \hat{=} B_{ijkl} C_{jkl}$ . The The scalar product is given with the notation, e.g., $A \cdot B = t r (A B^{T})$ and the Frobenius norm with, e.g., $‖ B ‖ = {(B \cdot B)}^{1 / 2}$ . For the dyadic product, we use the formulation, e.g., a ⊗ b and for the ℓ times repeated dyadic product of a vector, e.g., b ^⊗ℓ = b ⊗ b ⋯ ⊗ b . The unit sphere in $R^{3}$ is denoted by S². For the elementary multiplication of two vectors, we use the notation, e.g.,

a ⊙ b \hat{=} d i a g (a) b .

(1)

Optimization on the configuration space of curved and constrained fibers

Parametrization of curved and constrained fibers

We consider a curved fiber of length L and diameter D, which we model as a polygonal chain with n congruent spherocylinders, i.e., cylinders with hemispherical ends.^38,39 Hence, each segment has a segment length of ℓ = L/n. For the parametrization of such a polygonal chain, following Lauff et al.,³⁹ we normalize the starting point x⁰ of a fiber

{\hat{x}}^{0} = {\hat{Q}}^{- 1} x^{0} with the matrix \hat{Q} = d i a g (Q_{1}, Q_{2}, Q_{3}),

(2)

comprising the cell dimensions Q_i (i = 1, 2, 3). Then, a fiber may be parametrized via the coordinate vector

q = {[\begin{matrix} {({\hat{x}}^{0})}^{T} & {(p^{1})}^{T} & \dots & {(p^{n})}^{T} \end{matrix}]}^{T}

(3)

with dimension

n_{q} = 3 (n + 1),

(4)

where the unit vectors p^a ∈ S² (a = 1, …, n) denote the directions of the individual segments. For an illustration of the core idea of the used discretization, see Figure 1. The considered fiber is modeled with n = 5 segments and parametrized via its starting point and the five directions of its segments. With the coordinate vector q, the midpoints of the segments may be expressed explicitly in the form

x^{a} (q) = \hat{Q} {\hat{x}}^{0} + \sum_{b = 1}^{a - 1} ℓ p^{b} + \frac{ℓ}{2} p^{a}, a = 1, \dots n .

(5)

As the directions p^a are assumed to be unit vectors, the configuration space is described by the set

R = \{q \in R^{n_{q}} |Φ_{1} (q) = 0\}

(6)

with the n-dimensional vector-valued constraint function

Φ_{1} (q) = {[\begin{matrix} Φ_{1}^{1} (q) & \dots & Φ_{1}^{n} (q) \end{matrix}]}^{T} with Φ_{1}^{a} (q) = \frac{1}{2} (‖ p^{a} ‖^{2} - 1), a = 1, \dots, n .

(7)

The configuration space in equation (6) describes fibers residing in the Euclidean space

R^{3}

without further restrictions on their location. Also for periodic unit cells, the description in equation (6) is applicable with little extra effort, see Lauff et al.³⁹ An example for a two-dimensional periodic unit cell with fibers as inclusions is shown in Figure 2. However, when considering fibers with restrictions on their location, i.e., due to the layered structure of a continuous-discontinuous fiber reinforced polymer (CoDicoFRP), the configuration space in equation (6) needs to be adapted according to the underlying box constraints.

Figure 1.

Discretization of a curved fiber as a polygonal chain.

Figure 2.

Fibers in a two-dimensional periodic unit cell.

In this manuscript, we aim to consider hard and soft box constraints. In case of hard box constraints, an entire fiber is restricted to its respective cell

C = [B_{1}, T_{1}] \times [B_{2}, T_{2}] \times [B_{3}, T_{3}], 0 \leq B_{i} \leq T_{i} \leq Q_{i}, i = 1,2,3,

(8)

e.g., to the bottom or the upper layer of a two-layered CoDicoFRP. An admissible configuration of a fiber within its cell is shown in Figure 3(a). Hard box constraints are enforced by restricting the starting point of a fiber and the endpoints of the individual segments to lie within the unit cell. Such a condition may be formulated via the vectorial constraint function

Φ_{H} (q) = {[\begin{matrix} {Φ_{H}^{0} (q)}^{T} & {Φ_{H}^{1} (q)}^{T} & \dots & {Φ_{H}^{n} (q)}^{T} \end{matrix}]}^{T}

(9)

with

\begin{aligned} Φ_{H}^{0} (q) = [\begin{matrix} - x^{0} \\ + x^{0} \end{matrix}] + [\begin{matrix} + B_{1} \\ + B_{2} \\ + B_{3} \\ - T_{1} \\ - T_{2} \\ - T_{3} \end{matrix}], Φ_{H}^{a} (q) = [\begin{matrix} - (x^{a} (q) + ℓ p^{a} / 2) + [\begin{matrix} B_{1} \\ B_{2} \\ B_{3} \end{matrix}] \\ + (x^{a} (q) + ℓ p^{a} / 2) - [\begin{matrix} T_{1} \\ T_{2} \\ T_{3} \end{matrix}] \end{matrix}], a = 1, \dots, n . \end{aligned}

(10)

Figure 3.

Illustrative unit cells with hard box constraints (a), soft box constraints with control parameter n_SBC = 0 (b) and soft box constraints with control parameter n_SBC = 2 (c), where the black dots represent the constrained fiber points.

Thus, hard constraints are fulfilled if each component of the vectorial constraint function (9) is smaller or equal than zero.

In case of soft constraints, only the segments’ midpoints of a fiber are restricted to their respective unit cell, see Figure 3(b), and, e.g., an interlocking between two adjacent layers may be realized. Then, the vectorial constraint function is defined by

Φ_{S} (q) = {[\begin{matrix} {Φ_{S}^{1} (q)}^{T} & \dots & {Φ_{S}^{n} (q)}^{T} \end{matrix}]}^{T} with Φ_{S}^{a} (q) = [\begin{matrix} - x^{a} (q) + [\begin{matrix} B_{1} \\ B_{2} \\ B_{3} \end{matrix}] \\ + x^{a} (q) - [\begin{matrix} T_{1} \\ T_{2} \\ T_{3} \end{matrix}] \end{matrix}], a = 1, \dots, n .

(11)

If a fiber is modeled with short segments, most parts of the fiber will not penetrate the cell boundaries. Hence, the difference between hard and soft box constraints decreases for smaller segment lengths. To enable an interlocking also for small segment lengths, we adapt the soft box constraints such that not all midpoints are constrained to their cell, see Figure 3(c). Therefore, we introduce the parameter n_SBC ≥ 0 to control the amount of constrained midpoints. Depending on this parameter, the vectorial constraint function in equation (11) accounts for the midpoints x^a with

\{\begin{matrix} a = 1, \dots, n, & if n_{SBC} = 0 \\ a = n_{SBC} + 1, \dots, (1 + 2 k) n_{SBC} + 1 \leq n, & if 0 < n_{SBC} < n / 2 \\ a = ⌊ n / 2 ⌋, & else \end{matrix}\},

(12)

where ⌊⋅⌋ denotes the floor function. Thus, e.g., for a fiber with n = 30 segments, in case of n_SBC = 3 the midpoints x^a (a = 4, 10, 16, 22, 28) are restricted and in case of n_SBC ≥ 15 only the midpoint x¹⁵ is restricted.

We permit to mix the box-constraint types coordinate-wise, e.g., by using no box constraints for the e ₁-direction, soft box constraints for the e ₂-direction and hard constraints for the e ₃-direction. Hence, the constraint functions (9) and (11) may not be active in all directions and the inactive parts need to be deleted from the formulations. To shorten the expressions, we merge the box constraints to

Φ_{BC} (q) = [\begin{matrix} Φ_{H} (q) \\ Φ_{S} (q) \end{matrix}] \leq 0

(13)

with the dimension n_BC, where the inequality is to be understood component-wise. Then, the configuration space of a curved and constrained fiber is given by the manifold with boundary

R^{BC} = \{q \in R^{n_{q}} |Φ_{1} (q) = 0, Φ_{BC} (q) \leq 0\} .

(14)

Adapted gradient descent approach on the configuration space of curved and constrained fibers

We consider an objective function

F : R^{BC} \to R, q \mapsto f (q),

(15)

optimized in the configuration space

q \in R^{BC}

(14). To minimize the objective function, we aim to use a gradient descent approach. However, the standard, i.e., Euclidean, gradient descent method may compute iterates that violate the constraints in equation (14). To ensure an optimization procedure on the optimization space, the gradient descent approach needs to compute iterates along the geodesics, i.e., the locally shortest paths between two points on a manifold.⁴¹ The equations governing these geodesics follow from minimizing the kinetic energy functional on the manifold, which may be understood as a constrained mechanical system without external forces and driven solely by the kinetic energy.^47,48 For the manifold of a double spherical pendulum, Betsch⁴⁷ derived the equations governing the geodesics in the form of a d’Alembert type constrained mechanical system. Based on this formulation, Lauff et al.³⁹ extended the equations for a polygonal chain parametrized via

q \in R

, i.e., without box constraints, which read

\begin{aligned} {P (q)}^{T} G \ddot{q} & = 0, \\ Φ_{1} (q) & = 0 \end{aligned}

(16)

with the null space matrix P(q) and the constant symmetric positive-definite metric matrix G. For a detailed derivation of equation (16) and the explicit formulas for the matrices P and G, we refer to Lauff et al.³⁹ Proceeding from equation (16), the geodesics on the configuration space

R^{BC}

may be obtained by enforcing that at each point along the path fulfills the box constraints (13).

With the equations governing the geodesics at hand, we may compute the exponential mapping

{exp}_{q} (v) = q_{v} (1)

(17)

for prescribed points

q \in R^{BC}

on the manifold and admissible tangent vectors

v \in T_{q} R^{BC}

. Here, q_v(1) refers to the endpoint of the geodesic

q_{v} : [0,1] \to R^{BC}

with the starting point q_v(0) = q and the initial tangent vector

{\dot{q}}_{v} (0) = v

. Using the exponential mapping to ensure admissible iterates, we adapt the intrinsic gradient descent approach³⁹ with the iterative procedure

q^{k + 1} = \exp_{q^{k}} (- τ_{k} \nabla_{q} F (q^{k})), k \in N,

(18)

where τ_k stands for a selected positive stepsize.

Exponential mapping on the configuration space of curved and constrained fibers

In this section, we introduce a numerical integration scheme for computing the exponential mapping of the manifold with boundary $R^{BC}$ (17). Following previous publications which consider the exponential mapping of manifolds describing a double spherical pendulum⁴⁷ and a polygonal chain,³⁹ we use the energy conserving time integration scheme presented by Gonzales.⁴⁹

We consider the time interval $I = [0,1]$ , which we decompose into n_t subintervals of equal time steps Δt = 1 / n_t. For each subinterval $I_{j} = [t_{j}, t_{j + 1}]$ starting at time t_j, the coordinate vector $q (t) \in R^{BC}$ and its time derivative $\dot{q} (t) \in T_{q} R^{BC}$ are approximated by the linear approximation formulas

\begin{aligned} q (t) \approx q_{j} + \frac{t - t_{j}}{Δ t} (q_{j + 1} - q_{j}), \\ \dot{q} (t) \approx v_{j} + \frac{t - t_{j}}{Δ t} (v_{j + 1} - v_{j}) . \end{aligned}

(19)

With this time discretization scheme, the discretized form of the differential-algebraic equation (16) takes the form

\begin{aligned} {P (q_{j + \frac{1}{2}})}^{T} \hat{H} (q_{j}, q_{j + 1}) = 0, \\ Φ_{1} (q_{j + 1}) = 0 \end{aligned}

(20)

with the vector-valued function

\hat{H} (q_{j}, q_{j + 1}) = \frac{2}{Δ t} \hat{G} [q_{j + 1} - q_{j}] - 2 \hat{G} v_{j}

(21)

under the additional inequality constraint accounting for the box constraints (13). The matrix P(q_j+1/2) denotes the null space matrix evaluated for the coordinate vector

q_{j + \frac{1}{2}} = \frac{q_{j} + q_{j + 1}}{2}

(22)

and

\hat{G}

represents the normalized metric matrix

\begin{aligned} \hat{G} = [\begin{matrix} {\hat{Q}}^{- 2} [n {\hat{Q}}^{2}] & {\hat{Q}}^{- 2} [f_{1} ℓ \hat{Q}] & {\hat{Q}}^{- 2} [f_{2} ℓ \hat{Q}] & \dots & {\hat{Q}}^{- 2} [f_{n} ℓ \hat{Q}] \\ {\hat{Q}}^{- 2} [f_{d, 1} 1_{3 \times 3}] & {\hat{Q}}^{- 2} [f_{2} ℓ^{2} 1_{3 \times 3}] & \dots & {\hat{Q}}^{- 2} [f_{n} ℓ^{2} 1_{3 \times 3}] \\ ⋱ & ⋮ \\ ⋱ & {\hat{Q}}^{- 2} [f_{n} ℓ^{2} 1_{3 \times 3}] \\ sym & {\hat{Q}}^{- 2} [f_{d, n} 1_{3 \times 3}] \end{matrix}] \\ with f_{a} = n - a + \frac{1}{2} and f_{d, a} = ε (ℓ, D) + ℓ^{2} (n - a + \frac{1}{4}), \end{aligned}

(23)

where

ε (ℓ, D)

denotes the volume-specific moment of inertia of a spherocylinder.⁵⁰ For more details on the derivation of the discretized differential-algebraic equation (20), we refer to Betsch⁴⁷ and Lauff et al.³⁹ The additional inequality constraint accounting for the box constraints (13) leads to the Karush-Kuhn-Tucker (KKT) conditions

\begin{aligned} {P (q_{j + \frac{1}{2}})}^{T} [\hat{H} (q_{j}, q_{j + 1}) + {\frac{\partial {\hat{Φ}}_{BC} (q_{j + 1})}{\partial q}}^{T} λ_{j}] = 0, \\ Φ_{1} (q_{j + 1}) = 0, \\ λ_{j} \geq 0, {\hat{Φ}}_{BC} (q_{j + 1}) \leq 0, λ_{j} ⊙ {\hat{Φ}}_{BC} (q_{j + 1}) = 0, \end{aligned}

(24)

where we use the normalized vectorial constraint function

{\hat{Φ}}_{BC} (q) = [\begin{matrix} {\hat{Q}}^{- 1} & 0_{3 \times 3} & \dots & 0_{3 \times 3} \\ 0_{3 \times 3} & {\hat{Q}}^{- 1} & \dots & 0_{3 \times 3} \\ ⋮ & ⋮ & ⋱ & 0_{3 \times 3} \\ 0_{3 \times 3} & 0_{3 \times 3} & \dots & {\hat{Q}}^{- 1} \end{matrix}] Φ_{BC} (q) .

(25)

In equation (24), $λ \in R^{n_{BC}}$ denotes the additional vector of unknown Lagrangian multipliers and the notation () ⊙ () represents the element-wise multiplication of two vectors (1). For any fixed parameter vector $μ \in {(R^{> 0})}^{n_{BC}}$ , the conditions

λ_{j} \geq 0, {\hat{Φ}}_{BC} (q_{j + 1}) \leq 0, λ_{j} ⊙ {\hat{Φ}}_{BC} (q_{j + 1}) = 0

(26)

hold if and only if the equation

{⟨λ_{j} + μ ⊙ {\hat{Φ}}_{BC} (q_{j + 1})⟩}_{+} = λ_{j}

(27)

is satisfied, where ⟨h⟩₊ = max(0, h) denotes the Macaulay bracket, applied component-wise to a vector. Hence, the system of equations and inequalities (24) may be equivalently rewritten as a system of equations

\begin{aligned} {P (q_{j + \frac{1}{2}})}^{T} [\hat{H} (q_{j}, q_{j + 1}) + {\frac{\partial {\hat{Φ}}_{BC} (q_{j + 1})}{\partial q}}^{T} λ_{j}] = 0, \\ Φ_{1} (q_{j + 1}) = 0, \\ {⟨λ_{j} + μ ⊙ {\hat{Φ}}_{BC} (q_{j + 1})⟩}_{+} - λ_{j} = 0 . \end{aligned}

(28)

For each subinterval $I_{j} = [t_{j}, t_{j + 1}]$ , the coordinate vector $q_{j} \in R^{BC}$ and the velocity vector $v_{j} \in R^{n_{q}}$ are given quantities at time t = t_j. The coordinate vector q_j+1 is obtained by solving the system of equation (28) with a non-smooth Newton scheme.^44–46 As initial guess for the coordinate vector q_j+1, we choose

q_{j + 1} = q_{j} + Δ t v_{j} .

(29)

At convergence of the non-smooth Newton scheme, the velocity vector for the next time step is obtained via

v_{j + 1} = \frac{2}{Δ t} (q_{j + 1} - q_{j}) - v_{j} .

(30)

For the parameter vector μ , we use $μ = 10 \cdot 1_{n_{BC}}$ . As convergence criteria of the Newton scheme, we choose that the norm of the residual vector is below 10⁻⁸ and that the inequality constraints are strictly satisfied. Following Lauff et al.,³⁹ we initially consider the complete time interval I with a single time step Δt = 1. If the non-smooth Newton scheme does not converge in ten iterations, we repeat the time integration with doubled subintervals. For Newton’s method applied to find roots of a sufficiently smooth function, backtracking⁵¹ is a popular strategy to globalize the convergence properties of Newton’s method, which is only characterized by local convergence in its traditional undamped form. For non-smooth problems, however, backtracking is not guaranteed to help upgrading the local convergence behavior of Newton’s method to global convergence. Therefore, backtracking is only used in the case that the inequality constraints are not violated in a time step at all. If backtracking is applied, we use the backtracking parameter β = 0.9, taken from literature.⁵² Thus, in each backtracking step, the stepsize is moderately reduced by 10%. For the fixed backtracking parameter β = 0.9, we have investigated the influence of the second backtracking parameter $α \in (0,0.5)$ on the runtime to solve the system of equation (28). The study has shown that the effect of the parameter α is negligible and we select α = 1/3. To still ensure global convergence of the non-smooth Newton’s method, we rely on a homotopy method, i.e., the previously described adaptive reduction of the time-step size, which permits us to remain in the basin of attraction of the sought root.

Please note that considering the unknown Lagrangian multipliers in equation (28) increases the dimension of the system equations significantly. Hence, it is important for the performance to consider only those components which are active in the current time step. If a fiber does not touch the boundary of its containing cell, the inequality constraints are considered inactive and should be deleted from the formulation to decrease the dimension and thus the runtime for solving the equations.

Microstructure generation of long fiber reinforced hybrid composites

Description of long fiber reinforced hybrid composites

We consider periodic microstructures comprising a rectangular cell $Q = [0, Q_{1}] \times [0, Q_{2}] \times [0, Q_{3}]$ , which is divided into m non-overlapping subcells with rectangular cell dimensions

C^{c} = [B_{1}^{c}, T_{1}^{c}] \times [B_{2}^{c}, T_{2}^{c}] \times [B_{3}^{c}, T_{3}^{c}], c = 1, \dots, m .

(31)

Hence, the volume of a subcell computes as

V_{c} = (T_{1}^{c} - B_{1}^{c}) (T_{2}^{c} - B_{2}^{c}) (T_{3}^{c} - B_{3}^{c}) .

(32)

In each cell, N_c long, curved fibers of constant diameter D_c without fiber-overlap are given. Thus, the total number of fibers inside the entire unit cell Q equals the sum

N = \sum_{c = 1}^{m} N_{c} .

(33)

The ordering of the fibers is chosen such that the first N₁ fibers belong to the first cell, the next N₂ fibers to the second cell and so on. The smallest fiber-index corresponding to the c^th-cell is denoted by i_c.

The fibers are modeled as polygonal chains with spherocylindrical segments and are parametrized via the coordinate vector q, see equation (3). The fiber lengths follow a prescribed fiber length distribution, e.g, the uniform, Gamma or Weibull distribution. In this manuscript, we consider cells filled with discontinuous or continuous fibers, where we assume a unidirectional fiber arrangement for the continuous fibers. To model the discontinuous fibers, we select the number of segments of the i^th-fiber such that the segment lengths ℓ_i is below a defined maximum segment length ${\bar{ℓ}}_{c}$

n_{i} = ⌈ \frac{L_{i}}{{\bar{ℓ}}_{c}} ⌉, ℓ_{i} = \frac{L_{i}}{n_{i}}, i = 1, \dots, N,

(34)

where the functional ⌈⋅⌉ computes the smallest integer larger or equal to a given real number. For the continuous and unidirectional fibers, we use a single segment, i.e., we model the fiber as a straight cylindrical geometry.

The fiber volume fraction within a subcell computes as

ϕ_{c} = \frac{π D_{c}^{2}}{4 V_{c}} L_{c}^{t o t a l} with the total length L_{c}^{t o t a l} = \sum_{i = i_{c}}^{i_{c} + N_{c}} L_{i}

(35)

according to the Pappus’ theorem.⁵³ The fiber orientation state is typically represented by the fiber orientation tensors,^42,43 where the realized volume-weighted fiber orientation tensors of second

A_{c}^{r}

and fourth order

A_{c}^{r}

of each cell compute as

A_{c}^{r} = \frac{1}{L_{c}^{t o t a l}} \sum_{i = i_{c}}^{i_{c} + N_{c}} \sum_{a = 1}^{n_{i}} ℓ_{i} {(p_{i}^{a})}^{\otimes 2}, A_{c}^{r} = \frac{1}{L_{c}^{t o t a l}} \sum_{i = i_{c}}^{i_{c} + N_{c}} \sum_{a = 1}^{n_{i}} ℓ_{i} {(p_{i}^{a})}^{\otimes 4},

(36)

where the vector

p_{i}^{a}

denotes the direction of the a^th-segment of the i^th-fiber. For fiber reinforced composites, the fiber orientation tensor of second order may be identified via micro-CT imaging.^54,55,33 To obtain an adequate approximation of the fiber orientation tensor of fourth order, closure approximation^56–62 may be used. For cells with unidirectional fiber arrangement, e.g., oriented in e ₁-direction, the fiber orientation tensor of second order is given by

A_{c}^{r} \hat{=} d i a g (1,0,0)

. In Figure 4, an example for a hybrid composite with two layers in e ₃-direction is shown. The bottom layer is reinforced with discontinuous fibers and the upper layer is reinforced with continuous fibers in a unidirectional fiber arrangement in e ₁-direction.

Figure 4.

Hybrid composite with two layers in e ₃-direction.

Extension of the fused sequential addition and migration (fSAM) algorithm

In this section, we provide information on the extensions to the original fused sequential addition and migration (fSAM) algorithm³⁹ to generate periodic hybrid microstructures composed in multiple cells with different prescribed characteristics, e.g., varying fiber length and orientation distributions or fiber volume fractions.

Comment on the fiber orientation term within the objective function

The fSAM algorithm aims to minimize the objective function

F (f_{1}, \dots, f_{N}) \hat{=} \frac{1}{2} \sum_{i, j = 1}^{N} \sum_{a = 1}^{n_{i}} \sum_{b = 1}^{n_{j}} {(δ_{i j}^{a b})}^{2} + \frac{w_{A}}{8} N ‖ A - A^{r} ‖^{2} + \frac{w_{ρ}}{2} \sum_{i = 1}^{N} \sum_{a = 1}^{n_{i} - 1} {({⟨\cos \bar{α} - p_{i}^{a + 1} \cdot p_{i}^{a}⟩}_{+})}^{2},

(37)

where the first term accounts for the non-overlap condition between the fiber segments, the second term considers the fiber orientation tensor and the third term restricts the angle between adjacent segments to the maximum angle

\bar{α}

. For more details on the terms and the weights, we refer to Lauff et al.³⁹ The coordinate vector f_i denotes an alternative parametrization of a fiber

f_{i} = {[\begin{matrix} {(x_{i}^{1})}^{T} & \dots & {(x_{i}^{n_{i}})}^{T} & {(p_{i}^{1})}^{T} & \dots & {(p_{i}^{n_{i}})}^{T} \end{matrix}]}^{T},

(38)

which is obtained from the coordinate vector q_i by computing the midpoints of the segments

x_{i}^{a}

, see equation (5).

Typically, we prescribe a tolerance for the relative error of the fiber orientation tensor $ϵ_{A}$ as a termination criterion for the optimization problem. However, due to the formulation in equation (37), the fiber orientation will be further corrected during the gradient descent scheme even if the orientation error is below the tolerance, but the other criteria are not fulfilled. As the fiber orientation term leads to non-zero gradients of the objective function for all fiber directions, it is necessary to numerically integrate the d’Alembert constrained mechanical system for all fibers. Thus, the fiber orientation term may result in the computation of the exponential mapping for fibers which already satisfy the maximum angle condition and do not intersect with other fibers, although the fiber orientation tensor is already satisfying its convergence criterion. As a result, the runtime increases significantly, especially at the end of the optimization procedure, as the fiber orientation condition converges first most of the time. To decrease the runtime, we formulate the objective function

F_{{⟨ \cdot ⟩}_{+}} (f_{1}, \dots, f_{N}) \hat{=} \frac{1}{2} \sum_{i, j = 1}^{N} \sum_{a = 1}^{n_{i}} \sum_{b = 1}^{n_{j}} {(δ_{i j}^{a b})}^{2} + \frac{w_{A}}{8} N {⟨‖ A - A^{r} ‖^{2} - ϵ_{A}⟩}_{+} + \frac{w_{ρ}}{2} \sum_{i = 1}^{N} \sum_{a = 1}^{n_{i} - 1} {({⟨\cos \bar{α} - p_{i}^{a + 1} \cdot p_{i}^{a}⟩}_{+})}^{2}

(39)

such that the term for the fiber orientation tensor is only be considered if the orientation error exceeds the prescribed tolerance

ϵ_{A}

Extension of the objective function for multiple subcells with different characteristics

The objective function in equation (39) only accounts for a single fiber orientation distribution and a uniform maximum angle within the entire cell. As we aim to consider m subcells with different characteristics, we need to distinguish between these areas. Hence, the objective function in equation (37) turns to

\begin{aligned} F_{hybrid} (f_{1}, \dots, f_{N}) \hat{=} \frac{1}{2} \sum_{i, j = 1}^{N} \sum_{a = 1}^{n_{i}} \sum_{b = 1}^{n_{j}} {(δ_{i j}^{a b})}^{2} + \frac{w_{A}}{8} \sum_{c = 1}^{m} N_{c} {⟨‖ A_{c} - A_{c}^{r} ‖^{2} - ϵ_{A}⟩}_{+} \\ + \frac{w_{ρ}}{2} \sum_{c = 1}^{m} \sum_{i = i_{c}}^{i_{c} + N_{c}} \sum_{a = 1}^{n_{i} - 1} {({⟨\cos {\bar{α}}_{c} - p_{i}^{a + 1} \cdot p_{i}^{a}⟩}_{+})}^{2}, \end{aligned}

(40)

where

A_{c}

denotes the desired fiber orientation tensor of fourth order of the c^th-subcell and

{\bar{α}}_{c}

the maximum angle of the c^th-subcell. The term for the non-overlap condition does not distinguish between different subcells, as the fibers may reach into adjacent subcells, depending on their box constraints. Thus, fibers of different subcells may intersect during the optimization procedure which is not admissible for convergence.

Implementation of soft and hard box constraints

Let us discuss the implementation of the soft and hard box constraints as part of the fSAM algorithm. Therefore, we consider two different procedures. On the one hand, the box constraints may be fulfilled in every iterative step, which we call the intrinsic approach. Then, the gradient steps are restricted on the adapted configuration space of the fibers $R^{BC}$ , see equation (14). Hence, moving along the geodesics of this manifold requires the numerical integration of the d’Alembert type constrained mechanical system in equation (16) under the condition that the box constraints (13) are fulfilled. In our investigations, it turned out that the term for the fiber orientation tensor converges slowly if fibers are close to their boundaries. Thus, we increase the weight of this term for all fibers which violate the inequality conditions during the previous numerical integration of the d’Alembert constrained mechanical system by a factor of 2.5.

Alternatively, the box constraints may only be enforced via the convergence criterion of the optimization problem, but not in every iterative step. Hence, we consider the d’Alembert type constrained mechanical system in equation (16) without the box constraints during the exponential mapping. To enforce the box constraints, we extend the objective function by a term accounting for the soft and hard box constraints

F_{penalty,BC} (f_{1}, \dots, f_{N}) \hat{=} F_{hybrid} (f_{1}, \dots, f_{N}) + \frac{w_{BC}}{2} \sum_{i = 1}^{N} {‖Φ_{BC} (q_{i})‖}^{2},

(41)

where the coordinate vector q_i is computed from the coordinate vector f_i. For this change of fiber parametrization, the normalized starting point computes as

{\hat{x}}_{i}^{0} (f_{i}) = {\hat{Q}}^{- 1} (x_{i}^{1} - \frac{ℓ_{i}}{2} p_{i}^{1}) .

(42)

The objective function

F_{penalty,BC} (f_{1}, \dots, f_{N})

has dimension (length)². Hence, the weight of the box constraints term w_BC needs to be dimensionless and we select w_BC = 5. Due to the additional penalty term within the objective function, we call this procedure the penalty approach.

Sampling and pre-optimization step

Before the main optimization problem is solved, a sampling step is conducted. Therefore, we first sample fiber lengths according to the prescribed fiber length distribution of each cell until the desired fiber volume fraction is reached. Then, we sample the midpoints and the directions of the fibers, assuming straight fibers first. Therefore, we follow Schneider,³⁸ by sampling the midpoints from a uniform distribution within the subcells and by sampling the directions from the ACG closure,^63,64 accounting for the prescribed fiber orientation tensors of second order. Subsequently, the segment number of each fiber is determined, see equation (34), and the fibers decompose into segments without changing their location in the cell.

Due to the sampling procedure, the realized fiber orientation tensors may differ from the desired ones. Additionally, the fibers may violate their box constraints, as only the fiber’s midpoints are surely located within their respective cell. To stabilize the convergence behavior of the algorithm, we aim to ensure that the initial fiber orientation tensor of the main optimization problem is within the prescribed tolerance. Furthermore, moving along a geodesic requires the starting configuration to be admissible, i.e., to lie on the manifold of curved and constrained fibers. To obtain a starting configuration for the main optimization according to these requirements, we consider a pre-optimization step subsequent to the sampling step. Therefore, we solve an optimization problem with the objective function

\begin{aligned} F_{pre} (f_{1}, \dots, f_{N}) \hat{=} \frac{w_{A}}{8} \sum_{c = 1}^{m} N_{c} {⟨‖ A_{c} - A_{c}^{r} ‖^{2} - ϵ_{A}⟩}_{+} \frac{w_{ρ}}{2} \sum_{c = 1}^{m} \sum_{i = i_{c}}^{i_{c} + N_{c}} \sum_{a = 1}^{n_{i} - 1} {({⟨\cos {\bar{α}}_{c} - p_{i}^{a + 1} \cdot p_{i}^{a}⟩}_{+})}^{2} \\ + \frac{w_{BC}}{2} \sum_{i = 1}^{N} {||Φ_{BC} (q_{i})||}^{2}, \end{aligned}

(43)

which coincides with the objective function of the penalty approach, see equation (41), except for the term of the non-overlap condition. For the weights, we follow the main optimization problem, except for the weight of term for the box constraints w_BC = 0.5. Solving the optimization problem in equation (43) results in a fiber arrangement with accurate fiber orientation realization where the fibers are in admissible configurations on the spaces

R_{i}^{BC}

, but may be overlapping.

Handling of the continuous and unidirectional fibers

Continuous and unidirectional fibers are modeled with a single straight cylindrical segment, which is arranged in one of the three coordinate directions e _i. The segment has the length of the parallel edge of its cell $(T_{i}^{c} - B_{i}^{c})$ . Due to this special resulting configuration, the ends of the segments are in touch when periodically repeating the subcell in fiber direction. Hence, the fiber is continuous.

The fSAM algorithm²⁴ prescribes a minimum distance to avoid stress peaks between close fibers. Additionally, a self-intersection scheme³² is used to detect long fibers intersecting themselves. Thus, the algorithm would detect self-intersections for the continuous fibers, although they shall touch at their end. To overcome this limitation, the self-intersection scheme is turned off for continuous fibers.

For subcells reinforced with continuous fibers, we prescribe a unidirectional fiber orientation state. Hence, we know that all segment directions need to be unidirectional for convergence. Thus, we do not change the fiber directions within the main optimization problem. To find a fiber arrangement with no fiber overlap, only the midpoints of the fibers are changed.

Computational investigations

Setup

We use an implementation of the fSAM algorithm written in Python with Cython extensions, where the collision checks between the segments and the numerical integration of the exponential mapping are parallelized with OpenMP. The runtimes were measured on a desktop computer with a 8-core Intel i7 CPU and 64 GB RAM.

As our reference setup, we work with a two-layered CoDicoFRP comprising a discontinuously fiber-reinforced bottom layer and a continuously fiber-reinforced upper layer. For both subcells, we consider a polypropylene (PP) matrix reinforced by E-glass fibers. Table 1 records the corresponding elastic moduli.^65,66 The fiber diameter is assumed to be D = 10 μm for all fibers. For the discontinuous phase, the fibers are modeled with a uniform fiber length of L = 1500 μm, i.e., a fiber aspect ratio of 150, the maximum segment length is set to

\bar{ℓ} = 50 μ m

and the maximum angle is chosen as

\bar{α} = 60 °

. Between all fibers, we ensure a minimum distance of 20% of the fiber diameter, i.e., 2.0 μm.

Table 1.

Material properties for the PP matrix and the E-glass fibers.^65,66

E-glass fibers	PP matrix
E = 72.0 GPa	E = 1.25 GPa
ν = 0.22	ν = 0.35

The fiber orientation tensor of fourth order is obtained by a closure approximation from a prescribed fiber orientation tensor of second order. As we use the exact closure^63,64 and fiber orientation tensors of second order are orthotropic, in general, the fiber orientation tensors of fourth order are orthotropic as well.⁶⁷ For the convergence of the optimization procedure of the fSAM algorithm, it is necessary that no fiber overlap is detected. Additionally, the relative error in the fiber orientation tensor of fourth order needs to be below 10⁻⁴ and the absolute error of the angle constraint is enforced to be below 10⁻². For the penalty approach, we enforce that the constrained fiber points exceed their boundaries 10⁻³D at most.

To compute the effective elastic properties, we use an FFT-based computational homogenization software^68,69 with a discretization on a staggered grid⁷⁰ and a conjugate gradient solver.^71,72 As termination criterion for the solver, we choose a relative tolerance of 10⁻⁵. The effective elasticity tensor is computed using six independent load cases. Based on the effective elastic stiffness, we approximate the effective orthotropic engineering constants^73,74 due to the orthotropic material symmetry of the microstructures.

On the algorithmic handling of the box constraints

In this section, we compare the intrinsic and the penalty approach enforcing the box constraints. Therefore, we generate two-layered microstructures with a cubic unit cell and cell dimensions Q_i = 500 μm. The bottom layer with discontinuous fiber-reinforcement has cell dimension $C^{1} = [0,500] \times [0,500] \times [0,400] μ m^{3}$ and we consider the fiber volume fractions

ϕ_{1} \in {20 %, 21 %, 22 %, 23 %, 24 %, 25 %} .

(44)

We prescribe the fiber orientation tensor of second order

A_{1} \hat{=} d i a g (0.77, 0.17, 0.06)

, experimentally measured for long fiber reinforced thermoplastics.⁷⁵ The upper layer with continuous and unidirectional fiber-reinforcement has cell dimension

C^{2} = [0,500] \times [0,500] \times [400,500] μ m^{3}

and a fiber volume fraction of ϕ₂ = 35%. The principal direction of the continuous fibers is the e ₁-direction. For both layers, we consider no box constraints for the e ₁-direction and e ₂-direction. Hence, the fibers may wrap around the unit cell due to the periodic boundary conditions. For the e ₃-direction, we apply hard constraints at height 0 μm and 500 μm, as well as soft box constraints at height 400 μm. A generated microstructure is shown in Figure 4.

For each fiber volume fraction and procedure, we consider ten generated microstructures and report on the measured runtimes in Figure 5. We observe an increase in runtime for packings with higher fiber volume fractions. For the lowest fiber volume fractions ϕ₁ ∈ {20%, 21%}, both procedures have similar runtimes, e.g., the penalty procedure generates microstructures with ϕ₁ = 21% in about 43s on average. However, for higher fiber volume fractions the runtime of the intrinsic approach is always below the runtime of the penalty approach. Additionally, the penalty approach has outliers with extremely high runtimes. For the fiber volume fractions 22%, 23%, 24% and 25%, the intrinsic approach is 1.23, 1.33, 2.26 and 2.21 times faster than the penalty approach, when including the outliers to this consideration.

Figure 5.

Runtimes for the microstructure generation using an approach with intrinsic box constraints control or via a penalty term.

Let us focus on the packing with the highest fiber volume fraction of ϕ₁ = 25% to have a detailed look on the runtimes. For this case, the minimum measured runtime is 86s for the intrinsic approach and 114s for the penalty approach. The maximal measured runtime is 221s for the intrinsic approach and 849s for the penalty approach. Hence, we observe that the minimum runtime observed for both approaches is rather close. However, this is not the case for the maximal runtime as the penalty approach suffers from extreme outliers. For both approaches, the variability in runtime mainly depends on the starting configuration of the fibers, which results from both the sampling and the pre-optimization step. If the starting configuration is challenging with respect to the main optimization problem, e.g., due to a high number of collisions or high angles between the segments, the runtime will increase. Whereas the intrinsic approach is capable of dealing with challenging starting configurations, the penalty approach suffers from convergence problems, especially with increasing complexity, e.g., higher fiber volume fractions. However, still the penalty approach leads to reasonable runtimes for the microstructure generation. Due to the better performance of the intrinsic approach, we will use this procedure for the following investigations.

Resolution study

To predict the effective properties with FFT-based computational homogenization methods, a mesh size needs to be chosen which is sufficiently fine to enable an accurate prediction of the effective properties.^72,76 For the resolution study, we consider the setup from the previous subsection with the fiber volume fraction ϕ₁ = 25% and material parameters from Table 1. The cubic cell with cell dimensions Q_i = 500 μm turns out to be representative, see Table 3. To select an adequate mesh size, we compute the effective properties for the four voxel edge-lengths h = 4.00 μm, 2.00 μm, 1.25 μm and 1.00 μm. Compared to the fiber diameter D = 10 μm, a fiber is resolved with 2.5, 5, 8, 10 voxels for decreasing mesh size. As the effort of the homogenization increases with the total voxel number, the finest mesh size with 500³, i.e., 125 ⋅ 10⁶, voxels lead to a significantly higher runtime than the coarsest mesh size with 125³, i.e., about 2 ⋅ 10⁶, voxels. In Figure 6, the same microstructure resolved with the different voxel-edge lengths h = 4.00 μm and h = 1.25 μm is shown.

Figure 6.

Two-layered microstructures resolved with the voxel edge-lengths h = 4.00 μm (a) and h = 1.25 μm (b).

We report on the approximated orthotropic engineering constants for all four voxel-edge lengths in Table 2. To assess the quality of the orthotropic approximation, also the orthotropic approximation error is listed. For the two coarsest mesh sizes, the orthotropic approximation error is 2.76% and 1.07%, dropping below 1% for the two finer mesh sizes. Hence, it turns out that the approximation of the orthotropic engineering constants is adequate. Due to the fiber-reinforcement in e ₁-direction of the upper layer and the preferred fiber arrangement in e ₁-direction of the bottom layer, the Young’s modulus E₁ exceeds the remaining Young’s moduli by multiples. To select a sufficiently fine mesh size, we compare the elastic moduli of the three coarser mesh sizes to the finest mesh size. For the shear modulus G₁₃, the highest relative difference is measured with 16.50% for the coarsest mesh size, reducing to 4.85% and 0.97% for the two finer mesh sizes. As the computed orthotropic approximation error and the resolution induced error is below 1% for the voxel-edge length h = 1.25 μm, we select this voxel-edge length as our mesh size.

Table 2.

Approximated orthotropic engineering constants for the voxel edge-lengths h = 4.00 μm, h = 2.00 μm, h = 1.25 μm and h = 1.00 μm.

h	E ₁	E ₂	E ₃	G ₂₃	G ₁₃	G ₁₂	${e r r}_{o r t h}$
μm	GPa	GPa	GPa	GPa	GPa	GPa	%
4.00	12.53	3.13	2.76	1.03	1.20	1.62	2.76
2.00	13.32	3.02	2.57	0.94	1.08	1.49	1.07
1.25	13.53	3.00	2.54	0.91	1.04	1.47	0.56
1.00	13.60	3.00	2.53	0.90	1.03	1.46	0.44

RVE study

For computational homogenization of random materials, the selected unit cell size is of key importance to ensure representativity. Due to practical reasons the unit cell needs to be finite, and the apparent effective properties of random materials include a certain degree of randomness as well. However, to obtain the effective properties as deterministic descriptors, the unit cell size needs to be large enough to cover the entire material statistics and ensure that the boundary conditions do not impact the results.^77,78 Such a unit cell is called representative volume element (RVE). As the computational effort for the homogenization depends on the unit cell size, the smallest possible RVE size should be selected.

To assess the representativity of a unit cell, two errors need to be monitored.^14,79 On the one hand, the random error, or dispersion, measures the variability of the realized apparent properties for a fixed unit cell size, i.e., may be monitored by the standard deviation of multiple realizations. On the other hand, the systematic error, or bias, quantifies the difference between the mean apparent properties and the effective properties. Typically, the effective properties are not known. However, to quantify the systematic error, the mean apparent properties of increasing unit cell sizes may be monitored.

For the RVE study, we work with the two-layered microstructure already used in the previous two sections. To study the representativity of different unit cell sizes, we consider cubic cells with cell dimensions Q_i = 300 μm, 500 μm and 700 μm, where the boundary between the two layers is always at height 0.8 Q₃, i.e., at 240 μm, 300 μm and 560 μm, respectively, for increasing cell dimensions. To ensure a sufficiently fine resolution, we use a voxel edge-length of h = 1.25 μm, according to the resolution study in Table 2. Thus, 240³, i.e., about 14 ⋅ 10⁶, voxels resolve the smallest unit cell and 560³, i.e., 176 ⋅ 10⁶, resolve the largest unit cell. For each unit cell size, a generated microstructure is shown in Figure 7. To investigate the representativity errors, we consider ten microstructures for each unit cell size. The mean and the standard deviation of the approximated engineering constants are listed in Table 3. We observe that the orthotropic approximation is suitable for all three unit cell sizes as the error is always below 1%. Also, the standard deviation is low for all unit cell sizes. For the systematic error, we report the highest relative error for the smallest unit cell size Q_i = 300 μm and the shear modulus G₁₃ with 2.86%, dropping to 0.95% for the unit cell size Q_i = 500 μm. Hence, even the smallest unit cell size has rather small representativity errors. To ensure that both the systematic and the random error are well below 1%, we select the unit cell size with cell dimensions Q_i = 500 μm for the further investigations.

Figure 7.

Generated microstructures for three different cubic cell-sizes Q_i.

Table 3.

Approximated orthotropic engineering constants for the cubic cell-sizes Q_i = 300 μm, Q_i = 500 μm and Q_i = 700 μm.

Q _i	E ₁	E ₂	E ₃	G ₂₃	G ₁₃	G ₁₂	${e r r}_{o r t h}$
μm	GPa	GPa	GPa	GPa	GPa	GPa	%
300	13.48 ± 0.04	3.02 ± 0.04	2.53 ± 0.01	0.90 ± 0.00	1.02 ± 0.01	1.47 ± 0.02	0.69 ± 0.20
500	13.54 ± 0.05	3.02 ± 0.03	2.54 ± 0.00	0.91 ± 0.00	1.04 ± 0.01	1.48 ± 0.01	0.60 ± 0.08
700	13.55 ± 0.04	3.02 ± 0.02	2.56 ± 0.01	0.92 ± 0.00	1.05 ± 0.00	1.48 ± 0.01	0.57 ± 0.04

We remark that the fSAM algorithm generates periodic microstructures, regardless which box constraints are applied to the fibers. Already in previous publications^38,24 it is emphasized that periodic microstructures with accurately realized characteristics, i.e., the fiber volume fraction or the fiber orientation tensor of fourth order, combined with periodic boundary conditions for the displacement field during the homogenization are advantageous for decreasing the RVE size.

On the difference between soft and hard box constraints

In this section, we study the influence of the selected box constraints on the generated microstructures. Therefore, we work with the two-layered CoDicoFRP already used in the previous sections, see Figure 4, and change the box-constraint type between the bottom and the top layer. For the study, we compare four different box-constraint types. First, we use the hard box constraints, where the complete fiber is enforced to be within its cell. Then, we consider the soft box-constraint type with three different parameter selections of n_SBC, which controls the amount of restricted midpoints. In Table 4, the three considered parameter selections of n_SBC and the corresponding restricted midpoints x^a are listed.

Table 4.

Parameter selections of n_SBC for the soft box-constraint type and the corresponding restricted midpoints.

n _SBC	0	1	2
Restricted midpoints x^a	a = 1, 2, 3, …	a = 2, 4, 6, …	a = 3, 7, 11, …

For each box-constraint type, one of the ten generated microstructures is shown in Figure 8. We observe that for the hard box constraint in Figure 8(a) the boundary between the two layers is extremely straight and all fibers are completely in their cell, as enforced. For the soft box constraint with n_SBC = 0, the microstructure in Figure 8(b) looks quite similar. Due to the small maximum segment length of $\bar{ℓ} = 50 μ m$ , the fibers are mainly located in their respective layer and interlocking is suppressed. During the compression molding of CoDicoFRP, the material flow leads to ply migration between the two phases.⁸⁰ To model such a ply migration, we also consider higher parameter selections of n_SBC, restricting less segment midpoints to their cell, see Table 4. Figure 8(c) shows a microstructure generated with n_SBC = 1. We observe that the discontinuous fibers reach into the adjacent layer. Hence, this parameter selection may be suitable to realize the interlocking of the CoDicoFRP layers. When further increasing the parameter n_SBC to 2, the discontinuous fibers wrap around the continuous fiber-reinforcement, see Figure 8(d). Whereas this is not representative for CoDicoFRP, such a strong connection may appear for other use cases, i.e., the layer boundaries of discontinuous fiber-reinforced materials with a skin-core–skin structure.^81,33,82

Figure 8.

Two-layered continuous-discontinuous fiber-reinforced composite with different box constraints at the layer interface.

Besides a visual comparison of the generated microstructures, we also study the influence of the box constraints on the effective elastic properties. Therefore, we consider ten microstructures for each box constraint and report on the approximated orthotropic engineering constants in Table 5. It turns out that the effective stiffness increases slightly with higher interlocking, where the differences are rather small and below 1% for all cases. However, the influence of the box constraints may be higher, e.g., when considering the damage behavior of the inter-laminar region.^83,84

Table 5.

Approximated orthotropic engineering constants for different box constraints.

	E ₁	E ₂	E ₃	G ₂₃	G ₁₃	G ₁₂
	GPa	GPa	GPa	GPa	GPa	GPa
Hard box constraint	13.53 ± 0.03	2.99 ± 0.02	2.54 ± 0.00	0.91 ± 0.00	1.05 ± 0.00	1.47 ± 0.01
Soft box constraint	13.54 ± 0.05	3.02 ± 0.03	2.54 ± 0.00	0.91 ± 0.00	1.04 ± 0.01	1.48 ± 0.01
Soft box constraint	13.59 ± 0.05	3.05 ± 0.04	2.54 ± 0.01	0.91 ± 0.00	1.04 ± 0.00	1.50 ± 0.01
(ℓ_SBC = 150 μm)
Soft box constraint	13.58 ± 0.03	3.05 ± 0.03	2.56 ± 0.01	0.91 ± 0.00	1.05 ± 0.00	1.50 ± 0.02
(ℓ_SBC = 250 μm)

Comparison to experimental data

Hemp-fiber reinforced DicoFRP

We start the validation section by comparing the effective properties of microstructures with either straight or curved fibers to experimental data. Therefore, we consider a hemp-fiber reinforced DicoFRP with polypropylene (PP) as the matrix. The elastic parameters of the materials are listed in Table 6. According to Bourgogne,⁸⁵ the fibers have a diameter of D = 30 μm, come with a fiber length of L = 250 μm and are processed with a fiber volume content of 17.3%. Hence, we consider short fibers with a fiber aspect ratio of r_a ≈ 8.3. Following the suggestion by Bourgogne,⁸⁵ we use an isotropic fiber orientation state for the considered composite. For the microstructure generation, we enforce a minimum distance between the fibers of 20%, i.e., 6.0 μm, and use a maximum segment length of 40 μm in case of curved fibers. Straight fibers are realized with a single segment per fiber. We select a cubic cell with cell dimensions Q_i = 750 μm and resolve the microstructures with a voxel edge-length of 3.75 μm, i.e., according to the resolution study, eight voxels per diameter. In Figure 9, microstructures with straight and curved fibers are shown. Based on ten realizations, the mean runtime to generate a microstructure with straight fibers is 2s and to generate a microstructure with curved fibers is 19s.

Table 6.

Elastic parameters for the PP matrix and the hemp fibers.⁸⁵

	E	ν
	GPa	-
PP matrix	1.18	0.36
Hemp fibers	35.00	0.25

Figure 9.

Microstructures with straight (a) and curved (b) fibers.

The experiments⁸⁵ revealed a Young’s modulus E₁ of about 2.06 GPa for the hemp-fiber reinforced composite. To compare the mechanical behavior of the generated microstructures with the experimental data, we also report on the mean Young’s modulus E₁ obtained from ten realizations. For the microstructures with straight fibers, we obtained a mean Young’s modulus of E₁ = 2.03 GPa and for the curved fibers we obtained a mean Young’s modulus of E₁ = 1.98 GPa. Hence, we observe that the influence of the curvature is rather small with a relative decrease of 2.5% compared to straight fibers. Additionally, the results show a good agreement with the experimental data for both cases.

Carbon-fiber reinforced CoDicoFRP

Besides the hemp-fiber reinforced DicoFRP, we consider a CoDicoFRP with two continuous fiber-reinforced shell layers and a core layer with discontinuous fiber-reinforcement. For the materials, we use a polyamide 6 (PA6) matrix (Technyl^®star, DOMO Chemicals GmbH) reinforced with carbon fibers (PX35, ZOLTEK). The considered elastic parameters for both materials are listed in Table 7. We model the matrix with isotropic material parameters, taking the Young’s modulus E from the data sheet of a commercially available polyamide⁸⁶ and the Poisson’s ratio ν from the literature.⁸⁷ For the carbon fibers, we assume a transversely isotropic material symmetry and use the elastic modulus E₁ from the data sheet of the PX35 carbon fibers.⁸⁸ For the remaining parameters, we use the material data of HTA carbon fibers (Toho Tenax)⁸⁹ with similar elastic modulus E₁. Fiedler et al.⁸⁹ report on the orthotropic engineering constants of the HTA carbon fibers. Based on the corresponding orthotropic stiffness, the transversely isotropic engineering constants in Table 7 are obtained by a projection on the transverse isotropic stiffness space.^73,74 The diameter of the carbon fibers is D = 7.2 μm.⁸⁸

Table 7.

Elastic parameters for the PA6 matrix^86,87 and the carbon fibers.^88,89

	E	ν
	GPa	-
PA6 matrix	3.7	0.4

	E ₁	E ₂	G ₁₂	ν ₂₃	ν ₁₂
	GPa	GPa	GPa	-	-
Carbon fibers	242	28	32.4	0.23	0.23

The CoDicoFRP parts were manufactured with the long fiber reinforced thermoplastic direct (LFT-D) compression molding process^90,91 at the Fraunhofer ICT in Pfinztal, Germany, co-molding the LFT plastificate with the continuous and unidirectional fiber-reinforced tapes. The tapes include a fiber mass fraction of 60%,⁹² i.e., a fiber volume fraction of 48%, and are fiber-reinforced in the flow direction, which is the e ₁-direction. For the Dico material, fiber content and fiber length measurements were conducted at the Fiber Institute Bremen (FIBRE), Germany, see⁹³ for more details. The fiber volume content measurement revealed a fiber volume fraction of 23.5% on average. With respect to the fiber lengths, a number-weighted mean fiber length of L = 1550 μm and a median value of L = 540 μm was measured. Additionally, the fiber orientation tensors of the Dico material were analyzed and tensile tests for the manufactured CoDicoFRP parts conducted at the IAM, KIT, Germany. For more details on the procedures for the fiber orientation determination and the tensile tests, we refer to Scheuring et al.^75,93 Compression molding leads to parts which may be separated in two areas, depending on the initial charge of the LFT plastificate, the charge and the flow area. For the fiber orientation determination and the tensile tests, at least five samples were taken from both areas. In this investigation, we consider the results of the flow area. The mean fiber orientation tensor of second order is

A_{Dico} \hat{=} [\begin{matrix} 0.85 & 0.07 & 0.00 \\ 0.07 & 0.14 & 0.00 \\ 0.00 & 0.00 & 0.01 \end{matrix}] .

(45)

The experimental results for the planar Young’s modulus are shown in Figure 10, where the angle φ represents the deviation from the flow direction. The measuring points

φ \in {0 °, 22.5 °, 45 °, 90 °, 135 °, 157.5 °}

(46)

were considered. Due to the continuous fiber-reinforcement in flow direction and the highest eigenvalue a₁ of the fiber orientation tensor (45), the Young’s modulus in flow direction at φ = 0° is the highest modulus with 43.14 GPa. The smallest planar Young’s modulus is at φ = 90° with 7.95 GPa.

Figure 10.

Comparison of the Young’s modulus in the e ₁- e ₂-plane of the CoDicoFRP.

We aim to generate RVEs for the considered CoDicoFRP. Therefore, we select the subcell dimensions

\begin{aligned} C^{1} = [0,600] \times [0,600] \times [0,25] μ m, \\ C^{2} = [0,600] \times [0,600] \times [25,275] μ m, \\ C^{3} = [0,600] \times [0,600] \times [275,300] μ m, \end{aligned}

(47)

such that the height of the cell (and of the subcells) is ten times smaller than the manufactured parts. According to the manufactured part, we pack the shell layers with a fiber volume fraction of 48% and prescribe a unidirectional fiber arrangement in e ₁-direction. For the core layer, we prescribe the measured fiber volume fraction of 23.5% and the second-order fiber orientation tensor in equation (45). To stay between the measured mean and median fiber length, we use a uniform fiber length of 720 μm, i.e., a fiber aspect ratio of 100. With respect to the algorithmic choices, we use the soft box-constraint type with n_SBC = 1 to model the interface between the subcells. Due to the decreased diameter of D = 7.2 μm, we consider a smaller minimum distance of 1.44 μm. In Figure 11, a synthetic microstructure for a CoDicoFRP is shown, including 1073 fibers with a total number of 11167 segments. For the ten generated microstructures, the mean runtime for the microstructure generation is 156min.

Figure 11.

Synthetic microstructure of a CoDicoFRP using experimentally measured characteristics.

Let us compare the effective stiffness of the synthetic microstructures with the manufactured parts. For the resolution, we use h = 0.8 μm, i.e., the diameter is resolved by 9 voxels and the entire cell consists of 750² × 375 ≈ 211 ⋅ 10⁶ voxels. Based on the computed stiffness of the ten generated microstructures, we compute the direction dependent planar Young’s modulus⁹⁴ and report on the mean values in Figure 10. We observe that the highest Young’s modulus is not computed at the angle φ = 0° but at the angle φ ≈ 1.5°. This shift results from the compression molding process of the LFT plastificates and is known in the literature.^95,96,75,93 Notice that the generated microstructures are still orthotropic as we prescribe a fiber orientation tensor of second order and use an exact closure^63,64 to obtain the fiber orientation of fourth order.⁶⁷ However, the considered bases { e ₁, e ₂, e ₃} does not represent the orthotropic axes as the second-order fiber orientation tensor of the Dico phase (45) is not given in its principle axis system, i.e., its matrix-representation is not diagonalized. The orthotropic axis system of the generated microstructures may be obtained with a normal right rotation of the bases { e ₁, e ₂, e ₃} by about 1.5° around the e ₃-axes. Compared to the experimental data, it turns out that the shape of the planar Young’s modulus shows good agreement but the quantitative values are slightly underestimated. For the Young’s modulus at φ = 0°, the computed value is 40.41 GPa, which is a relative difference to the experiments of 6.3%, and at φ = 90°, the computed value is 7.28 GPa, i.e., a relative underestimation of 8.4%.

Summary and conclusion

This manuscript deals with an extension of the fused sequential addition and migration (fSAM)³⁹ algorithm to generate hybrid composites with long fiber reinforcement. The fSAM algorithm models the fibers as polygonal chains and uses an optimization procedure on the configuration space of a polygonal chain to generate microstructures. To distinguish between different cells of hybrid composites, an adapted configuration space of a polygonal chain is used to account for box constraints, which restrict the location of a fiber to its respective subcell. Depending on the considered material system, the subcells may be strictly separated or interlocked by fibers reaching into adjacent subcells. To enable different interfaces, the fSAM algorithm works with multiple box-constraint types. In case of the hard box constraints, the entire fiber is restricted to its cell, whereas for the soft box constraints only the segment midpoints of a polygonal chain are enforced to be within their cell. Hence, for the latter box constraint, an interlocking between adjacent subcells is enabled. However, when considering small segment lengths, the degree of interlocking is rather small as the fibers are mainly located within their subcell. To ensure interlocking also for decreasing segment lengths, an adapted soft box-constraint type may be selected where only selected segment midpoints are restricted.

Due to the connectivity constraint of a polygonal chain, the optimization problem of the original fSAM algorithm is defined on a curved, i.e., non-flat, configuration space. To compute admissible iterates, an adapted gradient descent approach is used, moving along the geodesics, i.e., the shortest path between two points on the configuration space. To account for the box constraints, the extended fSAM algorithm provides two alternative approaches. On the one hand, the penalty approach uses the original iteration rule which leads to iterates violating the box constraints. To ensure that the solution configuration at convergence is admissible, an additional penalty term accounting for the box constraints is added to the objective function. On the other hand, the intrinsic approach is based on an iteration rule which accounts directly for the box constraints. As the original iteration rule of the fSAM algorithm is only capable of handling equality constraints, a novel strategy is included in this manuscript to move along the geodesics of manifolds whose descriptions include inequality constraints as well, e.g., the box constraints. Hence, each iterate of the intrinsic approach is admissible on the configuration space of a curved and constrained fiber, i.e., a penalty term for the box constraints is not necessary. A detailed derivation of the extended iteration rule is given within the work at hand. For both strategies, the manuscript includes a discussion on the implementation within the framework of the fSAM algorithm. Additionally, information on further adaptions, which are necessary to generate hybrid composites, is provided, concerning the sampling, the pre-optimization step and the adaption of the objective function.

The computational investigations lead to the following findings:

• Due to the faster microstructure generation, the intrinsic approach is advantageous compared to the penalty approach to handle the box constraints. Especially for complex microstructures with high fiber volume fractions, the penalty approach suffers from convergence problems in case of challenging starting configurations, e.g., with many fiber collisions.

• The fSAM algorithm for hybrid composites generates fiber microstructures with high confidence in the desired properties, e.g., the fiber volume fraction or the fiber orientation tensor. Thus, even small unit cell sizes are representative and may be used as RVEs. As the runtime for the microstructure generation and the computational homogenization increases for larger unit cell sizes, the small RVE sizes ensure a computationally efficient prediction of the effective mechanical behavior of the material.

• The shape of the interfaces between adjacent subcells strongly depends on the selected box-constraint type. Whereas the hard box-constraint type leads to a strict separation, the soft box-constraint type enables interlocking, which increases with less restricted segment midpoints. Although the visual comparison results in significant differences, the effective elastic stiffness is in a close range for all implemented box-constraint types. However, the box-constraint type may have a higher influence, e.g., when investigating the damage behavior of the inter-laminar region.

• The microstructures generated by the fSAM algorithm are capable of representing the mechanical behavior of long fiber reinforced hybrid composites. The computed effective stiffness of a synthetic CoDicoFRP with shell-core-shell structure shows a good agreement with experimental measured data. The computed Young’s moduli E₁ and E₂ are 40.41 GPa and 7.28 GPa, slightly underestimating the experimental results with a relative difference of 6.3% and 8.4%, respectively.

To conclude, the studies demonstrate the capability of the extended fSAM algorithm to generate representative volume elements for long fiber reinforced hybrid composites in an efficient manner. Due to the possibility of choosing the box-constraint type, the interface may be modeled according to the considered material system. Hence, either strictly separated areas or interlocked areas can be realized by the microstructure model. With respect to the computational effort of the subsequent homogenization, the fSAM algorithm is a resource-saving method due to the small RVE sizes. By comparing the mechanical behavior of synthetic and manufactured CoDicoFRPs, the novel methodology permits to predict the measured elastic properties with high fidelity.

Besides CoDicoFRPs with long, curved fibers, as shown in Figure 12(a), the fSAM algorithm is also applicable to other material systems. Let us give two examples. For CoDicoFRPs with short fibers in the discontinuous fiber-reinforced phase, bending is typically negligible. Thus, the fSAM algorithm generates the corresponding microstructures by modeling the fibers with one cylindrical segment, see Figure 12(b). Another field of application are purely discontinuous fiber-reinforced materials with a skin-core–skin structure,^81,33,82 where the descriptors vary in each layer, see Figure 12(c).

Figure 12.

Three examples for hybrid or multilayer fiber-reinforced composites generated with the fSAM algorithm.

As a future development, one could account for the generation of interphases between the fibers and the matrix material,⁹⁷ especially when investigating interphase damage.⁹⁸ Additionally, strategies are necessary to generate microstructures with higher fiber volume fractions and fiber aspect ratios as, e.g., highly aligned fiber-reinforced microstructures are manufactured with fiber volume fractions above 50% and fiber aspect ratios above 500.^99,100 Throughout this manuscript, only the linear elastic material behavior of the composites was investigated. In perspective, it also might be of interest to study the damage behavior^52,83,84 of the generated microstructures and to take the nonlinear thermoviscoelastic behavior of the matrix material¹⁰¹ into consideration.

As a key feature, the fSAM algorithm models the fibers as flexible structures including curvature. Different studies^{102,38,103,85} report that microstructures with varying degrees of curvature show differences in their mechanical properties. Hence, it is important to study the influence of the fiber curvature on the computed effective stiffness. The fSAM algorithm permits only to decide whether fibers can be curved or not. Prescribing the degree of curvature is not (yet) possible. Additionally, for long fibers, the fiber volume content that can be achieved with straight fibers is extremely limited. Hence, for high fiber aspect ratios and fiber volume fractions, the influence of the curvature cannot be studied at all. To overcome these limitations, it could be part of future research to extend the fSAM algorithm by curvature control such that the influence of the fiber curvature is assessable and the degree of curvature may be prescribed as additional desired characteristic.

Footnotes

Acknowledgements

We thank the anonymous reviewers for their insightful comments and valuable suggestions which led to improvements of the manuscript. The research documented in this manuscript was funded by the German Research Foundation (DFG), project number 255730231, within the International Research Training Group “Integrated engineering of continuous-discontinuous long fiber reinforced polymer structures” (GRK 2078). The support by the German Research Foundation (DFG) is gratefully acknowledged. Support from the European Research Council within the Horizon Europe program - project 101040238 - is gratefully acknowledged by Matti Schneider.

Authors contribution

Celine Lauff: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – review & editing, Visualization, Project administration. Matti Schneider: Conceptualization, Methodology, Writing – review & editing, Supervision, Project administration, Funding acquisition. Thomas Böhlke: Resources, Writing – review & editing, Supervision, Project administration, Funding acquisition. All authors read and approved the manuscript.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Deutsche Forschungsgemeinschaft (255730231) and HORIZON EUROPE European Research Council (101040238).

ORCID iDs

Celine Lauff

Thomas Böhlke

References

Görthofer

Meyer

Pallicity

, et al. Virtual process chain of sheet molding compound: development, validation and perspectives. Compos B Eng 2019; 169: 133–147.

Schüle

Eyerer

. Polymer engineering. Berlin: Springer, 2020, 2.

Meyer

Gajek

Görthofer

, et al. A probabilistic virtual process chain to quantify process-induced uncertainties in Sheet Molding Compounds. Compos B Eng 2023; 249: 110380.

Böhlke

Henning

Hrymak

, et al. Continuous—discontinuous fiber-reinforced polymers — an integrated engineering approach. Munich: Carl Hanser, 2019.

Matouš

Geers

MGD

Kouznetsova

, et al. A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. J Comput Phys 2017; 330: 192–220.

Fish

Wagner

Keten

. Mesoscopic and multiscale modelling in materials. Nat Mater 2021; 20: 774–786.

Elmasry

Azoti

El-Safty

, et al. A comparative review of multiscale models for effective properties of nano- and micro-composites. Comput Mech 2022; 132: 101022.

Kozlov

. Averaging of differential operators with almost periodic rapidly oscillating coefficients. Math USSR Sbornik 1978; 107 (149)(2 (10)): 199–217.

Papanicolaou

Varadhan

SRS

. Boundary value problems with rapidly oscillating random coefficients. Random fields, Vol. I, II (Esztergom, 1979), Colloq. Math. Soc. János Bolyai. Amsterdam-New York: North-Holland, 1981, Volume 27, pp. 835–873.

10.

de Paiva

Bisiaux

Lynch

, et al. High resolution x-ray tomography in an electron microprobe. Rev Sci Instrum 1996; 67(6): 2251–2256.

11.

Shen

Nutt

Hull

. Direct observation and measurement of fiber architecture in short fiber-polymer composite foam through micro-CT imaging. Compos Sci Technol 2004; 64(13–14): 2113–2120.

12.

Buck

Brylka

Müller

, et al. Two-scale structural mechanical modeling of long fiber reinforced thermoplastics. Compos Sci Technol 2017; 117: 159–167.

13.

Kuhn

CIW

Taeger

, et al. Experimental and numerical analysis of fiber matrix separation during compression molding of long fiber reinforced thermoplastics. J Compos Sci 2017; 1(1): 2.

14.

Kanit

Forest

Galliet

, et al. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solid Struct 2003; 40(13–14): 3647–3679.

15.

Sab

Nedjar

. Periodization of random media and representative volume element size for linear composites. Compt Rendus Mec 2005; 333(2): 187–195.

16.

Schneider

Josien

Otto

. Representative volume elements for matrix-inclusion composites - a computational study on the effects of an improper treatment of particles intersecting the boundary and the benefits of periodizing the ensemble. J Mech Phys Solid 2022; 158: 104652.

17.

Widom

. Random sequential addition of hard spheres to a volume. J Chem Phys 1966; 44(10): 3888–3894.

18.

Feder

. Random sequential adsorption. J Theor Biol 1980; 87(2): 237–254.

19.

Tian

Zhou

, et al. Representative volume element for composites reinforced by spatially randomly distributed discontinuous fibers and its applications. Compos Struct 2013; 131(7): 366–373.

20.

Chen

Zhou

, et al. Study of the effectiveness of the RVEs for random short fiber reinforced elastomer composites. Fibers Polym 2019; 20(7): 1467–1479.

21.

Bahmani

Willett

, et al. Three-dimensional micromechanical assessment of biomimetic composites with non-uniformly dispersed inclusions. Compos Struct 2019; 212: 484–499.

22.

Tian

Chao

, et al. An advanced method for efficiently generating composite RVEs with specified particle orientation. Compos Sci Technol 2021; 205: 108647.

23.

Bahmani

Nooraie

Willet

, et al. A sequential mobile packing algorithm for micromechanical assessment of heterogeneous materials. Compos Sci Technol 2023; 237(2): 110008.

24.

Lauff

Schneider

Montesano

, et al. An orientation corrected shaking method for the microstructure generation of short fiber-reinforced composites with almost planar fiber orientation. Compos Struct 2023; 323: 117352.

25.

Coelho

Thovert

Adler

. Geometrical and transport properties of random packings of spheres and aspherical particles. Phys Rev 1997; 55(2): 1959–1978.

26.

Pan

Iorga

Pelegri

. Numerical generation of a random chopped fiber composite RVE and its elastic properties. Compos Sci Technol 2006; 68: 56–66.

27.

Lubachevsky

Stillinger

. Geometric properties of random disk packings. J Stat Phys 1990; 60(5–6): 561–583.

28.

Torquato

Jiao

. Dense packings of polyhedra: platonic and archimedean solids. Phys Rev 2009; 80(4): 041104.

29.

Torquato

Jiao

. Robust algorithm to generate a diverse class of dense disordered and ordered sphere packings via linear programming. Phys Rev 2009; 82(6): 061302.

30.

Williams

Philipse

. Random packings of spheres and spherocylinders simulated by mechanical contraction. Phys Rev 2003; 67(5): 1–9.

31.

Schneider

. The Sequential Addition and Migration method to generate representative volume elements for the homogenization of short fiber reinforced plastics. Comput Mech 2017; 59: 247–263.

32.

Mehta

Schneider

. A sequential addition and migration method for generating microstructures of short fibers with prescribed length distribution. Comput Mech 2022; 70(4): 829–851.

33.

Hessman

Riedel

Welschinger

, et al. Microstructural analysis of short glass fiber reinforced thermoplastics based on x-ray micro-computed tomography. Compos Sci Technol 2019; 183: 107752.

34.

Fliegener

Luke

Gumbsch

. 3D microstructure modeling of long fiber reinforced thermoplastics. Compos Sci Technol 2014; 104: 136–145.

35.

Nguyen

Bapanapalli

Holbery

, et al. Fiber length and orientation in long-fiber injection-molded thermoplastics — Part I: modeling of microstructure and elastic properties. J Compos Mater 2008; 42(10): 1003–1029.

36.

Senior

Osswald

. Measuring fiber length in the core and shell regions of injection molded long fiber-reinforced thermoplastic plaques. J Compos Sci 2020; 4(3): 104.

37.

Mehta

Schneider

. A maximum-entropy length-orientation closure for short-fiber reinforced composites. Comput Mech 2024; Online First: 1–26.

38.

Schneider

. An algorithm for generating microstructures of fiber-reinforced composites with long fibers. Int J Numer Methods Eng 2022; 123(24): 6197–6219.

39.

Lauff

Schneider

Montesano

, et al. Generating microstructures of long fiber reinforced composites by the fused sequential addition and migration method. Int J Numer Methods Eng 2024; 125(22): e7573.

40.

do Carmo

. Riemannian geometry. Basel: Birkhäuser, 1992.

41.

Petersen

. Riemannian geometry. In: Graduate texts in mathematics. Heidelberg: Springer, 2016, vol 171.

42.

Kanatani

. Distribution of directional data and fabric tensors. Int J Eng Sci 1984; 22(2): 149–164.

43.

Advani

Tucker III

. The use of tensors to describe and predict fiber orientation in short fiber composites. J Rheol. 1987; 31: 751–784.

44.

Sun

. A nonsmooth version of Newton’s method. Math Program 1993; 58(1–3): 353–367.

45.

Robinson

. Newton’s method for a class of nonsmooth functions. Set-Valued Anal 1994; 2(1): 291–305.

46.

Ito

Kunisch

. Lagrange multiplier approach to variational problems and applications. Philadelphia: SIAM Publishing, 2008.

47.

Betsch

. The discrete null space method for the energy consistent integration of constrained mechanical systems: Part I: holonomic constraints. Comput Methods Appl Mech Eng 2005; 194(50–52): 5159–5190.

48.

Angleses

Lee

. The modelling of holonomic mechanical systems using a natural orthogonal complement. Trans Can Soc Mech Eng 2005; 13(4): 81–89.

49.

Gonzales

. Mechanical systems subject to holonomic constraints: differential-algebraic formulations and conservative integration. Physica D 1999; 132(132): 165–174.

50.

García

. Random packings via mechanical contraction. Utrecht: Universiteit Utrecht, 2015.

51.

Boyd

Vandenberghe

. Convex optimization. Cambridge: Cambridge University Press, 2004.

52.

Meyer

Gajek

Görthofer

, et al. A convex anisotropic damage model based on the compliance tensor. Int J Damage Mech 2022; 31(1): 43–86.

53.

Goodman

. Generalizations of the theorems of pappus. Am Math Mon 1969; 76(4): 355–366.

54.

Krause

Hausherr

Burgeth

, et al. Determination of the fibre orientation in composites using the structure tensor and local X-ray transform. J Mater Sci 2010; 45(4): 888–896.

55.

Wirjadi

Schladitz

Easwaran

, et al. Estimating fibre direction distributions of reinforced composites from tomographic images. Image Anal Stereol 2016; 35(3): 167–179.

56.

Folgar

Tucker III

. Orientation behaviour of fibers in concentrated suspensions. J Reinforc Plast Compos. 1984; 3: 98–119.

57.

Cintra

Tucker III

. Orthotropic closure approximations for flow-induced fiber orientation. J Rheol. 1995; 39(6): 1095–1122.

58.

Dray

Gilormini

Régnier

. Comparison of several closure approximations for evaluating the thermoelastic properties of an injection molded short-fiber composite. Compos Sci Technol 2007; 67(7–8): 1601–1610.

59.

Breuer

Stommel

Korte

. Analysis and evaluation of fiber orientation reconstruction methods. J Compos Sci 2019; 3(3): 67.

60.

Kugler

Kech

Cruz

, et al. Fiber orientation predictions - a review of existing models. J Compos Sci 2020; 4(2): 69.

61.

Karl

Gatti

Frohnapfel

, et al. Asymptotic fiber orientation states of the quadratically closed Folgar-Tucker equation and a subsequent closure improvement. J Rheol 2021; 65(5): 999–1022.

62.

Karl

Schneider

Böhlke

. On fully symmetric implicit closure approximations for fiber orientation tensors. J Non-Newtonian Fluid Mech 2023; 318: 105049.

63.

Montgomery-Smith

Jack

, et al. Exact tensor closures for the three-dimensional Jeffery’s equation. J Fluid Mech 2011; 680: 321–335.

64.

Montgomery-Smith

Jack

Smith

. The fast exact closure for Jeffery’s equation with diffusion. J Non-Newtonian Fluid Mech 2011; 166: 343–353.

65.

DOW^® Chemical Company . C711–70RNA polypropylene resin. Product information, 2003.

66.

PPG fiber glass TufRov^® 4575. Data Sheet, 2013.

67.

Bauer

Böhlke

. On the dependence of orientation averaging mean field homogenization on planar fourth-order fiber orientation tensors. Mech Mater 2022; 170: 104307.

68.

Moulinec

Suquet

. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes Rendus de l’Académie des Sciences Série II 1994; 318(11): 1417–1423.

69.

Moulinec

Suquet

. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput Methods Appl Mech Eng 1998; 157: 69–94.

70.

Schneider

Ospald

Kabel

. Computational homogenization of elasticity in a staggered grid. Int J Numer Methods Eng 2016; 105(9): 693–720.

71.

Zeman

Vondřejc

Novak

, et al. Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients. J Comput Phys 2010; 229(21): 8065–8071.

72.

Brisard

Dormieux

. Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites. Comput Methods Appl Mech Eng 2012; 217–220: 197–212.

73.

Cowin

. The relationship between the elasticitiy tensor and the fabric tensor. Mech Mater 1985; 4(2): 137–147.

74.

Görthofer

Schneider

Ospald

, et al. Computational homogenization of sheet molding compound composites based on high fidelity representative volume elements. Comput Mater Sci 2020; 174: 109456.

75.

Blarr

Sabiston

Krauß

, et al. Implementation and comparison of algebraic and machine learning based tensor interpolation methods applied to fiber orientation tensor fields obtained from CT images. Comput Mater Sci 2023; 228: 112286.

76.

Schneider

. Convergence of FFT-based homogenization for strongly heterogeneous media. Math Methods Appl Sci 2015; 38(13): 2761–2778.

77.

Hill

. Elastic properties of reinforced solids: some theoretical principles. J Mech Phys Solid 1963; 11(5): 357–372.

78.

Drugan

Willis

. A micromechanics-based nonlocal constitutive equations and estimates of representative volume element size for elastic composites. J Mech Phys Solid 1996; 44: 497–524.

79.

Gloria

Otto

. An optimal variance estimate in stochastic homogenization of discrete elliptic equations. Ann Probab 2011; 39(3): 779–856.

80.

Corbridge

Harper

De Focatiis

DSA

, et al. Compression moulding of composites with hybrid fibre architectures. Composites Part A 2017; 95: 87–99.

81.

Bay

Tucker III

. Fiber orientation in simple injection moldings. Part I: theory and numerical methods. Polym Compos. 1992; 13(4): 317–331.

82.

Breuer

Spickenheuer

Stommel

. Statistical analysis of mechanical stressing in short fiber reinforced composites by means of statistical and representative volume elements. Fibers 2021; 9(5): 32.

83.

Ernesti

Schneider

. Computing the effective crack energy of heterogeneous and anisotropic microstructures via anisotropic minimal surfaces. Comput Mech 2021; 69(3): 45–57.

84.

Ernesti

Schneider

. Accounting for weak interfaces in computing the effective crack energy of heterogeneous materials using the composite voxel technique. Arch Appl Mech 2017; 93(10): 3983–4008.

85.

Bourgogne

QCP

. Influence of fibre curvature and orientation on the mechanical properties of injected irregular short hemp fibres reinforced polypropylene. J Thermoplast Compos Mater 2024; 37(9): 2963–2986.

86.

Altech^®. Technical datasheet ALTECH PA6 A 1000/677, 2020.

87.

Tschoegl

Knauss

Emri

. Poisson’s ratio in linear viscoelasticity-a critical review. Mech Time-Dependent Mater 2002; 6: 3–51.

88.

ZOLTEK Corporation . Technical datasheet ZOLTEK^TM PX35 vinyl ester compatible tow, 2019.

89.

Fiedler

Holst

Hobbiebrunken

, et al. Modelling of the initial failure of CFRP structures by partial discretisation: a micro / macro-mechanical approach of first ply failure. Compos Adv Mater 2004; 13: 219–226.

90.

Krause

Henning

Tröster

, et al. LFT-D — a process Technology for large scale production of fiber reinforced thermoplastic components. J Thermoplast Compos Mater 2003; 289–302(4): 117352.

91.

Christ

Scheuring

Schelleis

, et al. Characterization and simulation of the interface between a continuous and discontinuous carbon fiber reinforced thermoplastic by using the climbing drum peel test considering humidity. Polymers 2024; 16(7): 976.

92.

DOMO . Technical datasheet TECHNYL^® LITE C130 C60, 2023.

93.

Scheuring

Christ

Blarr

, et al. Experimental and homogenized orientation-dependent properties of hybrid long fiber-reinforced thermoplastics. Int J Mech Sci 2024; 280: 109470.

94.

Böhlke

Brüggemann

. Graphical representation of the generalized Hooke’s law. Tech Mech 2001; 21(2): 145–158.

95.

Bondy

Pinter

Altenhof

. Experimental characterization and modelling of the elastic properties of direct compounded compression molded carbon fibre/polyamide 6 long fibre thermoplastic. Mater Des 2017; 122: 184–196.

96.

Schelleis

Scheuring

Liebig

, et al. Approaching polycarbonate as an LFT-D material: processing and mechanical properties. Polymers 2023; 15(9): 2041.

97.

Riaño

Joliff

. An Abaqus^TM plug-in for the geometry generation of representative volume elements with randomly distributed fibers and interphases. Compos Struct 2019; 209: 644–651.

98.

Bourgogne

QCP

Bouchart

Chevrier

, et al. Numerical investigation of the fiber/matrix inter-phase damage of a PPS composite considering temperature and cooling liquid ageing. SN Appl Sci 2021; 3: 1–15.

99.

Cender

Simacek

Gillespie

JJW

, et al. Microstructural evolution of highly aligned discontinuous fiber composites during longitudinal extension in forming. Compos Sci Technol 2024; 254: 110649.

100.

Yarlagadda

Deitzel

Heider

, et al. Tailorable universal feedstock for forming (TUFF): overview and performance. In: Proceedings of the International SAMPE Technical Conference. Charlotte, NC, USA, 2019.

101.

Kehrer

Keursten

Hirschberg

, et al. Dynamic mechanical analysis of PA 6 under hydrothermal influences and viscoelastic material modeling. J Thermoplast Compos Mater 2023; 36(11): 4630–4664.

102.

Zhang

Yang

Sheng

. Performance-driven 3D printing of continuous curved carbon fibre reinforced polymer composites: a preliminary numerical study. Composites Part B 2018; 151: 164–256.

103.

Dai

Wei

Tao

, et al. Mechanical model of stiffness coefficients prediction of curved fiber reinforced composites considering fiber distribution and aggregation. Compos Struct 2023; 321: 117277.