An algorithm for generating microstructures of fiber‐reinforced composites with long fibers

We describe a sequential addition and migration (SAM) algorithm for generating microstructures of fiber‐reinforced composites with a direct control of the magnitude of curvature of the fibers. The algorithm permits to generate microstructures with fibers that are significantly longer than the edge lengths of the underlying cell. Industrially processed short and long‐fiber composites naturally feature a high volume fraction, which needs to be reflected by state‐of‐the‐art microstructure generation tools. Nowadays, it is well understood that digital twins of the microstructure of composites are essential for reliable computational multiscale methods. The original SAM algorithm was shown to reliably generate microstructures for short and straight cylindrical fibers. Digital volume images reveal, however, that the fibers in such fiber‐reinforced composites may show significant curvature, in particular for long fibers. The work at hand introduces an extension of the original SAM approach to curved fibers. More precisely, curved fibers are considered as sequences of straight fibers which are joined at their respective ends and whose level of bending is controlled by the angle between adjacent fiber segments. We discuss how to efficiently implement the novel method and how to select the crucial numerical parameters. We compare the introduced methodology to the original SAM algorithm for short fibers and demonstrate the superiority of the novel strategy for long fibers.

ites naturally feature a high volume fraction, which needs to be reflected by state-of-the-art microstructure generation tools. Nowadays, it is well understood that digital twins of the microstructure of composites are essential for reliable computational multiscale methods. The original SAM algorithm was shown to reliably generate microstructures for short and straight cylindrical fibers. Digital volume images reveal, however, that the fibers in such fiber-reinforced composites may show significant curvature, in particular for long fibers. The work at hand introduces an extension of the original SAM approach to curved fibers.
More precisely, curved fibers are considered as sequences of straight fibers which are joined at their respective ends and whose level of bending is controlled by the angle between adjacent fiber segments. We discuss how to efficiently implement the novel method and how to select the crucial numerical parameters. We compare the introduced methodology to the original SAM algorithm for short fibers and demonstrate the superiority of the novel strategy for long fibers.

K E Y W O R D S
curved fibers, long-fiber reinforced thermoplastics, microstructure generation, representative volume element, sequential addition and migration Computational multiscale methods, on the other hand, rest upon three fundamental pillars. For a start, they require a geometry which resolves the heterogeneity on the lower scale. 4 Second, efficient computational tools [5][6][7] need to be available to compute the mechanical response of the considered microstructure provided expressive material models are given. Last but not least, for nonlinear material behavior it is necessary to provide means of upscaling the results in a form that permits computing on component scale in reasonable time. [8][9][10] The work at hand is focused on generating microstructures for short and long fiber reinforced composites. Although digital images are typically available for such composites, for instance via microcomputed tomography, 11,12 they are typically insufficient unescorted and should be complemented by synthetic digital twins. Indeed, only a limited number of such micrographs is available and their statistics is a consequence of the rather involved production process. In particular, it is typically impossible to create a desired state, for example, a desired filler fraction and fiber orientation, also due to the stochastic influences. Moreover, extracting geometrical data from digital images may be difficult, for instance as a result of the segmentation process. [13][14][15] As an example, it is difficult to estimate the fiber-length distribution from micrographs, see Hessman et al. 16 for a discussion. Also, the influence of the fiber orientation is pronounced. [17][18][19] Algorithms for generating discontinuous composites essentially come in two flavors: sequential insertion algorithms and collective rearrangement algorithms. The former handle each fiber individually, adding the fibers one by one to the volume. The most widely used algorithm for generating fiber-filled microstructures belongs to this class and is called random sequential adsorption/addition. 20,21 Applied to cylindrical fibers, it samples the fiber directions randomly, typically based on a prescribed fiber-orientation distribution, and samples the fiber positions sequentially conditioned on having no overlap with previously placed fibers. The fiber-volume fractions which can be achieved by RSA depend on the fiber orientation and the (average) fiber aspect ratio (length/diameter). 22,23 Aligned fibers can be packed quite densely, whereas the packing achievable for isotropic fibers decreases as one over the aspect ratio. 22 Still, the RSA method is quite popular, presumably due to its simplicity of realization, and triggered quite a number of extensions and modifications. [24][25][26] Another insertion algorithm is based on sequential deposition, 27,28 where fibers are dropped and curve as a result of "gravity." This type of method realizes larger volume fractions than RSA for straight cylinders, but is restricted to essentially planar fiber-orientation states.
For sequential insertion algorithms, a fiber is never moved or rotated once it is placed. This fact is responsible for the inability of this class of methods to reach high filler fractions. The second class of fiber-generation methods are collective rearrangement algorithms. For this type of methods, all fibers move and rotate simultaneously. In particular, this strategy appears to be closer to most manufacturing processes prevalent for fiber composites like injection or compression molding.
Collective rearrangement algorithms may be classified into two categories, depending on how they treat the non-interpenetration constraint. The first category considers arrangements of non-interpenetrating fibers and targets rearranging them in such a way that the volume fraction is increased, but the non-interpenetration condition continues to hold. For instance, Williams and Philipse 29 present the mechanical contraction method designed for spherocylinders, that is, cylinders with spherical half-caps attached. They prepare a nonoverlapping configuration at low volume fraction by RSA. Then, the cell is continuously shrunk but the condition of nonoverlap is maintained via a suitable algorithm which displaces and rotates the spherocylinders if necessary. To the same category of approaches belong methods inspired by molecular dynamics and which originate from the Lubachevsky-Stillinger (LS) algorithm 30 that was devised for spheres. These methods 31,32 start with a nonoverlapping configuration at dilute concentration, for instance prepared by RSA. To each particle, a velocity and an angular velocity is assigned. Then, the particles follow Newton's equation of motion unless they collide with another particle, which causes the colliding particles to exchange momentum. Moreover, the considered cell shrinks with time, so that the volume content increases, as well. Such methods are typically implemented in an event-driven fashion, and stopped when the time between consecutive collision events drops below a prescribed threshold. LS-type methods require an efficient method for resolving collisions between particles to be fast. The collision detection for ellipsoidal particles is rather costly, limiting the potential of the approach. 31,32 Surprisingly, these types of collective rearrangement particles are not much better than plain RSA. In fact, the analysis of Toll 23 reveals that, both for RSA and nonoverlapping collective rearrangement algorithms, the volume fraction to be reached depends on the number of contacts which preclude the movement of the fiber. For straight fibers, this number may be assumed independent of the fiber length. In turn, the average number of contacts per fiber increases with fiber length. Thus, except for the aligned state, the achievable volume fraction decreases with increasing fiber aspect-ratio.
On the other hand, it is known that long-fiber reinforced thermoplastics (LFT) materials can be processed to a much higher volume fraction than predicted by the theory. 33,34 The analysis of Toll 23 is true, so where is the catch? Actually, there are (at least) two ways to reach higher volume fraction. A real fiber will not remain perfectly straight during the manufacturing process, that is, it will bend as a result of the experienced deformation. When fibers are considered as flexible, the number of average contacts with other fibers will remain the same, but the number of contacts required to restrain the motion will be significantly larger. Thus, it appears imperative to consider flexible fibers.
Fliegener et al. 35 demonstrate that the mechanical contraction method does indeed work for LFT when elastic deformation of the fibers is permitted. These authors account for the peculiarities of LFT, that is, their predominantly planar fiber orientation, by contracting the volume element only in out-of-plane direction. Unfortunately, their approach requires a full-fledged finite element analysis involving hundreds or thousands of accurately resolved fibers with an associated time integration and contact resolution.
Another way to circumvent the restrictions of Toll's analysis 23 is to consider collective rearrangement algorithms which permit overlap during the algorithm's inner workings. Indeed, a state of nonoverlapping is only required at convergence. To understand the underlying mechanisms think of the mechanical contraction method applied to a configuration with two crossing fibers. Suppose it would be possible to pack the fibers more tightly by moving one fiber across the other one. Such a move is not possible for nonoverlapping configuration preserving algorithms, but perfectly feasible when overlap is permitted during the runtime.
Altendorf and Jeulin 36 present such a force-based overlap-removal algorithm for short-fiber composites. They consider flexible fibers described as a chain of spheres. They sample the initial positions, directions and curvatures of the fibers at the beginning, ignoring overlap constraints. Then, they apply a force-based overlap removal algorithm, which is similar in spirit to the overlap-removal strategy of Williams and Philipse, 29 but augmented by constraints on the distance between adjacent spheres and on the angles between consecutive fiber segments. Altendorf and Jeulin 36 demonstrate that their method is able to generate extremely high fiber-volume fractions for arbitrary fiber-orientation states provided the curvature constraint is relaxed. In a subsequent study, 37 the effective elastic and thermal properties for the developed model of short-fiber reinforced composites was studied.
Another representative for this class of algorithms was introduced by Schneider. 38 In contrast to Altendorf and Jeulin, 36 he considered straight fibers, but used a similar force-based overlap-removal strategy, augmented by forces ensuring the fidelity in the realized fiber-orientation state. The computational experiments demonstrated that it was possible to generate microstructures of short-fiber reinforced composites for industrial volume fractions (up to 30%) and industrial aspect ratios (up to 100) in reasonable time with fibers that are straight. The analysis in Altendorf et al. 37 revealed that the fiber curvature had little influence of the effective elastic and thermal properties for small and moderate fiber aspect-ratio, so fidelity in the microstructures generated by Schneider, 38 via the sequential addition and migration (SAM) algorithm, is appropriate.

Contributions
In this work, we address the problem of generating fiber-filled volume elements with long fibers. We want to satisfy two requirements. First, the size of the volume elements should be computationally reasonable. Second, we aim for a periodic description of the geometry. In fact, it has been known since a long time that nonperiodic geometries may adversely affect the accuracy of the computed apparent properties. 39,40 This fact is exemplified for straight, overlapping and infinitely long fibers by Dirrenberger et al. 41 who show that in such a situation, representative volume elements (RVEs) are expected to be truly gigantic as a result of the slow decorrelation of the ensemble over large distances. A similar behavior is expected for straight cylindrical fibers if the volume elements are not large enough relative to the fiber length. 42,43 For this purpose, we extend the SAM algorithm 38 to curved fibers. More precisely, we discretize each fiber as a chain of straight fiber segments. 44,45 The original SAM algorithm is augmented by connectivity constraints between consecutive line segments and restrictions on the allowed fiber curvature, encoded by a maximum allowed angle between adjacent segments. We discuss implementation details, and conduct expressive computational studies. After identifying the proper resolution and RVE size, we compare the effective elastic properties of straight and curved short fibers. Moreover, we demonstrate the capabilities of the algorithm to generate volume elements of long-fiber reinforced composites at an industrially relevant volume fraction. Specifically, we demonstrate the possibility to include fibers whose fiber length significantly exceeds the volume-element size.

Describing a microstructure with curved fibers
Suppose a rectangular cell is given which contains N curved cylindrical fibers, see Figure 1A. We assume that all fibers have the same diameter D and that we may describe the centerline of each fiber (i = 1, … , N) in terms of a twice continuously differentiable curve where L i > 0 stands for the length of the fiber and the curve c i is assumed to be parameterized by arc length, that is, the identity see Figure 1B. We assume periodic boundary conditions for the cell Q, and interpret the parameterization (2) to respect periodicity. In particular, fibers are allowed to wrap around the cell, see Figure 1A.
For each fiber c i , the curvature is defined as see Figure 1B. Assuming that the cylindrical fibers with diameter D described by Equation (2) are nonintersecting, the volume of the ith fiber computes as by (the generalized) Pappus' theorem, 46 that is, it coincides with the volume of a straight cylinder. In case all of the N fibers are nonoverlapping, the fiber-volume fraction inside the cell Q computes as Moreover, the volume-weighted fiber-orientation tensors of order two and four compute as and respectively. The formulas (5) and (6) generalize the classical expressions 47,48 for straight fibers and are typically identified on digital-volume images anyway. [13][14][15] In particular, the usual conditions hold, as a consequence of the parameterization by arc length (3).
For the article at hand, we target generating such fiber-reinforced microstructures for a given cell Q, fiber lengths L i , volume fraction and fourth-order fiber-orientation tensor A. Moreover, we impose a bound on the fiber curvature

Discretized problem description
We discretize a curved fiber as follows. 44,45 We consider the centerline to be described in terms of a polygonal chain. The volumetric fiber as a set then consists of the straight cylinders described by the straight lines in the polygonal chain and spheres at the joints whose diameter coincides with the fiber's diameter. More precisely, we suppose that the ith fiber (i = 1, … , N) consists of n i segments of length i = L i ∕n i , where L i refers to the length of the ith fiber. Each segment (a = 1, … , n i ), in turn, is described by its centroid m i a ∈ Q and its direction p i a ∈ S 2 , an element of the unit sphere. The setup is illustrated in Figure 2.
The condition of forming a polygonal chain may be algebraically encoded by the n i − 1 constraints where d Q denotes the periodic distance on the cell Q. We refer to Figure 3A for a schematic. The discretized fourth-order fiber-orientation tensor (6) computes as F I G U R E 2 Representing the ith fiber in terms of straight cylindrical segments with spherical joints (left) and underlying description in terms of a polygonal chain (right) We approximate the curvature constraint (7) by finite differences For convenience, we rewrite the latter condition in the form Put differently, the discretized version of the curvature constraint (7) reads expressed in terms of the angle i ∈ [0, ], implicitly determined by the condition Actually, it turns out to be convenient to prescribe the angle i in the condition (10) directly, see Figure 3B. Such a strategy is reasonable as long as the segment length i remains fixed. Moreover, the fibers should neither overlap themselves nor each other. To simplify the overlap detection, we attach spherical caps to the cylindrical segments, and consider them as spherocylinders, 29 see Figure 4. Suppose that segment a of fiber i and segment b of fiber j with either i < j or i = j and a < b − 1 are considered. Then, the nonoverlap condition Equation (11) reads min s i a ,s see Figure 4.

An SAM algorithm
Let us collect the centroids and the directions of the fiber segments into (large) vectors ⃗ m and ⃗ p. To satisfy the conditions, the nonoverlapping requirement (11), the connectivity constraint (8), the restriction on the curvature (10) and the specification of the fiber-orientation tensor we utilize an optimization framework, extending the groundwork. 38 More precisely, we introduce a function F which is non-negative and whose roots characterize an admissible fiber configuration. To construct the function F, we encode the constraints in terms of suitable functions. Equality constraints are handled via quadratic penalty terms. For the inequality constraints, we use the old trick in terms of the Macaulay bracket ⟨h⟩ + = max(0, h). Then, to encode the nonoverlap condition (11), we introduce the quantity for i < j and all a and b as well as a for the curvature constraint (10). The connectivity constraint (8) is encoded by the quantity , a = 1, … , n i − 1.
With this notation at hand, we define the function F as follows in terms of suitable positive weights w , w A , and w . By construction, the roots of the function F encode nonoverlapping configurations which satisfy the curvature constraint, the connectivity constraint and match the prescribed fourth-order fiber-orientation tensor. Some care has to be taken with the weights. Their relative magnitude governs the importance assigned to the individual constraints. Moreover, the weights do not all have the same dimension. The function F has dimensions (length ) 2 , as dictated by the first summand in Equation (14). The second summand in itself has the same dimensions, so that we can select dimension one for the weight w . We found it beneficial to choose w = 0.75. In this way, overlap removal is prioritized over satisfaction of the connectivity constraint between consecutive segments. The third and the fourth summand in Equation (14) have dimension one, so that the weights w A and w need to have dimensions (length ) 2 . In the groundwork, 38 choosing the weight was suggested, where denotes the segment length. As we will stick to a uniform segment length throughout this work, i ≡ , we may retain this formula, as well. By the same reasoning as before, we select w = 0.75 w A , prioritizing the fiber-orientation tensor over the angle constraint.
The optimization problem is nonlinear for several reasons. For a start, the inequality-encoding strategy (13) introduces a nonlinearity. Moreover, the variables live in a nonlinear space. Indeed, the centroids ⃗ m take their values in a Cartesian product of the cell Q, which is nonlinear due to the periodic boundary conditions. The directions ⃗ p live in a Cartesian product of the unit sphere S 2 , which is also a nonlinear space.
Yet, the objective function (14) is continuously differentiable (but not twice continuously differentiable), and we may employ a gradient-descent scheme to find critical points. Optimality of a critical point of the program (15) is easily checked. As the objective function (14) is a sum of non-negative terms, any root ( ⃗ m, ⃗ p) of the function F is a global minimum. This minimum is not unique unless in specially constructed situations. Indeed, usually one may rattle (via a rigid-body motion) at specific fibers and retain optimality.
Some care has to be taken with gradient descent on curved spaces, as naive (i.e., extrinsic) gradient descent may leave the space. Instead, we use an intrinsic approach via the Riemannian exponential mapping. 49(Ch3) For a general (geodesically complete) Riemannian manifold (M, g), the exponential mapping exp takes any point x ∈ M, together with a tangent vector v ∈ T x M and moves along the geodesic, that is, the locally shortest line, emanating from the point x in direction of the tangent vector v. In our case, the Riemannian manifold is composed of (products of) the flat torus Q and the (positively curved) unit sphere S 2 . Geodesics on the flat torus Q correspond to straight lines, accounting for Q-periodicity. On the unit sphere S 2 , geodesics correspond to great circles.
With this machinery at hand, the SAM algorithm computes iterates of the form in terms of a suitable step size s. We follow previous work 38 and set s = 0.501. The rationale behind this choice is that if only two fibers (or segments) overlap, a step size s = 0.5 will push both objects apart to resolve the overlap, see Figure 4. Choosing the step size s to be slightly larger provides a small comfort zone beyond removing the overlap. Choosing the step size significantly larger than s = 0.5 precludes reaching high volume fractions. In turn, selecting a smaller step size leads to a less aggressive algorithm which tends to be (significantly) slower. For the function (14) at hand, the components of the gradient compute as 38 where k ij ab denotes the vector realizing the minimum distance (11) between the line segments and s ij ab ∈ [−1, 1] refers to the minimizing parameter, see Figure 4. Similarly, v i a denotes the vector realizing the minimum distance (8) between successive line segments (and where we formally set v i n i ≡ 0). We also scale the orientation part by a prefactor of 0.3 to facilitate segment translation over segment rotation.
We are led to an extended SAM algorithm, 38 Add fiber-coherence terms v i a − v i a+1 to Equation (17) 6: Add fiber-angle penalization-term i a p i a+1 − i a−1 p i a−1 to Equation (17) 7: Add orientation-penalization term to Equation (17) 8: Project the direction term onto the tangent space via Id − p i a ⊗ p i a and rescale by 1∕ 9: Perform the update (16) by addition in space and along the great circle for the direction 10: end while 11: end while

Algorithmic choices
This section provides some essential details which make the implementation work. We also added some explanations for the reasons behind some choices.

Prescribed fiber orientation
For straight fibers which all have identical length it is known both theoretically and from computational experiments that the second-order fiber-orientation tensor (5) is insufficient for fixing the elastic behavior of the composite, see Müller and Böhlke. 50 Thus, it is imperative to incorporate the fourth-order fiber-orientation tensor (6), at least. Computational experiments 38,50 suggest that choosing the fourth-order fiber-orientation tensor is sufficient for this purpose, although the situation is less clear in the context of the effective inelastic response. In this context, we refer to work on the higher-order fiber-orientation tensor closures. 51,52 For straight fibers of variable length, the situation is less clear. It is known that the volume-weighted fiber-orientation tensor of fourth order is not sufficient to determine the effective elastic properties at high contrast, see Bertóti et al. 53(eq.(3.6)) Instead, it would be necessary to consider also higher moments in the length. Unfortunately, computing such moments is not a simple task for digital volume images of fiber-reinforced composites, and the state of the art computes volume-weighted fiber-orientation tensors. 13,[54][55][56] For curved fibers, we again follow this line of thought, and prescribe the fiber-orientation data which is available, that is, arises from digital image processing. 13,[54][55][56] In this way, we were led to considering Equations (5) and (6).
Working with fourth-order fiber-orientation tensors is possible, but tedious. For instance, it is not clear (to the author) which fourth-order tensors arise as fourth-order fiber-orientation tensors of an actual fiber-orientation distribution, see also Linn 57( §2) for further details. This problem may be circumvented if results of digital image processing are available. Indeed, those arise from fiber-orientation distributions by construction.
In other cases, it is convenient to prescribe the second-order fiber-orientation tensor (5) and to compute the fourth-order fiber-orientation tensor (6) via a closure approximation, 58 that is, a function which associates a reasonable fourth-order fiber-orientation tensor to any second-order fiber-orientation tensor. Of course, the data encoded by the fourth-order fiber-orientation tensor cannot be fully recovered by the closure approximation. 50,52 Still, this strategy is still used for convenience and practical relevance. A variety of closure approximations is available, see Breuer et al. 52 for a recent overview. We select the exact closure approximation of Montgomery-Smith et al., 59,60 as it arises from an exact solution of (a simplified model for) fiber-orientation dynamics.

Sampling
For generating the initial configuration ( ⃗ m, ⃗ p), we rely upon a two-phase strategy. In the first phase, the centroids ⃗ m and the directions ⃗ p are sampled randomly. We consider the fibers to be straight initially. The centroids of the fibers are sampled randomly from a uniform distribution on the cell Q, whereas the directions ⃗ p are drawn from the angular central distribution 61 which gives rise to the exact closure. We refer to Ospald et al. 62 for details. Once the fibers are sampled, they are decomposed into individual segments, retaining the original shape of the sampled fiber.
Due to the random sampling, the error between the targeted fiber-orientation tensor A and the realized fiber-orientation tensor A N may be large(r than desired), in particular for a small number of fibers. Indeed, the error of a classical Monte Carlo sampling decays as N −1∕2 for N fibers. 63( §7.6) Thus, to reach a relative error of 1%, on the order of 10,000 fibers are required. To compensate, we set up a second phase during sampling which reorients the fibers to match the prescribed fiber-orientation tensor A to the desired accuracy. For this purpose, we use an optimization approach based on the objective function which coincides with the function (14) in the SAM algorithm except for the overlap removal, which is responsible for the bulk of the effort of the SAM algorithm. In this way, the fiber orientation is corrected accounting for the connectivity and angle constraints (8) and (10). When a small number of fibers with many segments is sampled, it turns out to be convenient to add small random perturbations to the directions upon sampling. Indeed, if only a small number of fibers is sampled, the prescribed fiber-orientation tensor may not be realizable. As the directions of all segments of a fiber coincide, there is no possibility to compensate this shortcoming due to bending. This shortcoming is removed by the small random perturbations.

Minimum distance
For elastic and inelastic analyses, it is convenient to prescribe a minimum distance between the fibers. Indeed, high strains and stresses may form when two fibers get close, and a certain degree of mesh resolution is required to provide accurate results. Classically, 20% of the fibers' diameter is necessary and sufficient for the minimum distance in elastic analyses and voxel-based discretizations. 64 A minimum inter-fiber distance is easily realized within the framework at hand by enlarging the fiber diameter D for the overlap detection in an artificial way. In postprocessing, the fibers are furnished with the original fiber-orientation tensor.

Overlap detection
The bulk of the effort in the SAM algorithm is hidden in evaluating the inter-segment distances (11). Classically, for two given segments in Euclidean space, the problem reduces to a quadratic program in two variables with four linear constraints. General-purpose methods are inefficient for the case at hand, and specific solutions are available. We utilize a modification of the Vega-Lago algorithm 65 presented in Pournin et al. 31 Some care has to be taken for parallel fibers, as the connecting vector is not unique in this case and the corresponding parameters s ij ab and s ji ba influence the overlap-removal algorithm. Moreover, the condition of being parallel needs to be quantified for finite numerical precision.
On a periodic cell Q of the form (1), only a single inter-segment distance needs to be computed if the segment length is smaller than min(Q 1 , Q 2 , Q 3 )∕2 by correcting the vector connecting the centers as for periodic distances between points. 66 Longer segments are handled by dividing the segments into sufficiently short parts.

Periodic boundary conditions
Periodicity is automatically enforced by working with the periodic distance d Q , used, for instance, in Equation (8). On the cell Q, see Equation (1), the periodic distance between two points x, y ∈ Y = [0, It is computationally more convenient to use the equivalent representation with the function for ∈ [−L, L] and L > 0. For fiber segments whose length is smaller than the half of the minimum of the Q i (i = 1, 2, 3), it is therefore appropriate to replace the distance between the centers by the version modified with  component-wise. For longer fiber segments, a decomposition approach may be used, described in Mehta and Schneider. 67( §3.3)

Cell-linked lists and nested Verlet lists
For M segments, checking all inter-segment distances naively requires on the order of M 2 overlap checks between segments. If the segments are rather short, cell-linked lists 68 provide a strategy for reducing the effort to O(M). In this approach, the cell Q is decomposed into a regular grid of rectangular subcells which can contain a segment fully. Each segment is associated to the cell which contains its centroid. Then, overlap checks are only made between adjacent cells and their containing segments. Unfortunately, for "long" segments, this approach may be slow despite its optimal complexity. Indeed, the runtime is determined by the number of segments K which are contained in a subcell (on average). The overlap checks with the neighboring cells leads to an O(K 2 ) effort. If K is not small, this effort may be significant as well. Therefore, we utilize Verlet lists 69 on top of the cell-linked lists. More precisely, we use nested Verlet lists. 38

Setup
The algorithm presented in Section 2.3 was implemented in Python with Cython extensions and parallelized with OpenMP. Compared to the C++ implementation, 38 the current implementation is slower (by about a factor of two) when executed with straight fibers. The slow-down of the current implementation is a consequence of three factors. For a start, C++ implementations tend to admit a higher degree of optimization than Cython code. Second, a number of corrections and safeguards were added to the current implementation. Last but not least, parts of the speed were sacrificed for maintaining (future) extensibility of the code. The effective elastic properties were computed by FFT-based computational homogenization methods 71,72 which exploit the benefits of regular grids. More precisely, we use a discretization on a staggered grid 73 and solve with a conjugate gradient method 74,75 up to a relative error of 10 −5 . 76( §3.6) We always compute the full effective stiffness tensor using the standard six load cases, and subsequently approximate the effective elasticity tensor respecting the underlying symmetry class of the material. For the materials, we use the material parameters of E-glass fibers and of polybutylene terephtalate (PBT), obtained from mechanical tests, 70 see Table 1. Different matrix materials may be used. However, their Young's modulus is typically in the range of 1 GPa to slightly above 3 GPa.
The introduced SAM algorithm and the elastic computations were run on a desktop computer with a 6-core Intel i7 CPU and 32GB RAM.
As the reference for our studies on fiber-reinforced composites, we consider fibers with a length of 250 μm and a diameter of 10 μm, filling 18% per volume. We use a minimum distance of 2 μm, segments with a length of 50 μm and a maximum angle of = 15 • . We will rely upon the exact closure approximation 59,60 throughout, parameterizing the fiber-orientation distribution via the second-order fiber-orientation tensor.
For the SAM algorithm, we enforce the nonoverlapping condition (11) strictly, whereas the fiber-orientation constraint (12) is solved up to a relative tolerance of 10 −4 . The angle constraint (10) is enforced up to an absolute error of 10 −2 , whereas the connectivity constraint (8) is required not to exceed 10 −3 D, depending on the fiber diameter D.

Resolution study
We start our investigations by studying the resolution which is necessary to provide meaningful results for the effective elastic properties. For this purpose, we consider fibers with a length L = 250 μm and a diameter D = 10 μm, enforcing an inter-fiber distance of 2 μm. The volume element is chosen as a cube with edge length Q 1 = Q 2 = Q 3 = 256 μm. Thus, more than two fibers fit along any edge of the cell. Such a cell size will turn out to be (more than sufficient to be) representative for the results in Section 3.3. We consider a fiber orientation described by the second-order fiber-orientation tensor that was determined experimentally in Magino et al.. 70 At this point, we would like to reiterate that we restrict to the exact closure approximation 59,60 for our analysis. We consider three voxel lengths h ∈ {4 μm, 2 μm, 1 μm}.
A generated volume element with a fiber-volume content of = 18% is shown in Figure 5 with varying degrees of mesh width. Curved fibers are readily apparent, as a result of the individual segments with a length of 50 μm, whose maximum angle is 15 • . For the coarse resolution, h = 4 μm, only 2.5 voxels lie across a typical fiber. For the medium resolution, five voxels resolve a typical fiber diameter, whereas ten voxels spread across a diameter for the fine resolution, h = 1 μm.
Due to the orthotropy of the second-order fiber-orientation tensor (22) and the fact that the exact closure preserves orthotropy, we expect the effective elastic properties to be orthotropic as well. Thus, we compare the engineering constants  of the orthotropic approximation of the effective stiffness and monitor the orthotropic approximation error err orth in addition. 77 The fiber-orientation tensor (22) reveals that the fibers lie predominantly in the x-y-plane, with a clear preference of the x-direction over the y-direction. This fact is also reflected in the effective Young's moduli shown in Table 2. The Young's modulus in x-direction is about twice as high as those in the transverse directions. Similarly, the shear modulus in the x-y-plane, G 12 exceeds the other two shear moduli significantly. By the way, discussing the orthotropic engineering constants is meaningful, as the orthotropic approximation error err orth is rather low. For the coarsest resolution, it is slightly below two percent. Doubling the resolution halves the orthotropic approximation error. For the finest resolution, we obtain an error below half a percent.

TA B L E 2 Results for the resolution study
For the coarsest resolution, the difference to the finest resolution considered is noticeable. The largest error is observed for the Young's modulus in x-direction, E 1 , with 8.1%. For a resolution of h = 2 μm, the relative error decreases to 1.7%. The other engineering constants are affected less severely. For instance, the transverse Young's moduli are only affected in the third significant digit.
Thus, we select a resolution of five voxels per diameter. This result is consistent with previous studies on straight fibers. 38,64 This indicates that a moderate curvature of the fibers does not necessitate a finer resolution than for straight fibers.

Representativity study
RVEs form a cornerstone of multiscale mechanics of materials with a random microstructure. For a stationary and ergodic medium, effective properties emerge on volume elements of sufficient size. [78][79][80] To assess the representativity of a considered volume, two quantities need to be monitored. The effective elastic properties of volume elements form random variables. The standard deviation of these effective properties on cells of fixed size constitute the dispersion of the effective properties, 80 also called random error in the mathematical literature. 81 This quantity measures how much different realizations of the random material differ in terms of their apparent properties on this cell size. The dispersion depends on the cell size and decreases for increasing cell volume. [82][83][84][85] In addition to the dispersion there is another error source when computing the effective properties. Indeed, let us consider the expectation of the effective elastic properties on cells of fixed size. This expectation is, in general, different from the infinite volume limit, for instance as a result of artificial correlations introduced by the finiteness of the cell. The difference between this expectation and the effective properties gives rise to the bias 80 or systematic error 81 in the mathematical literature.
To quantify bias and dispersion, we consider the following setup. We examine ten realizations for different cell sizes and compute the mean and the standard deviation of the orthotropic engineering constants. By coincidence, a sample size of ten also gives rise to the 99% confidence interval for Student's t-distribution. Figure 6 shows examples of the generated cells. The results are listed in Table 3. We observe that the dispersion is rather small for all considered volume elements, and that it decreases with increasing volume. Moreover, for the shear moduli, the bias is extremely small, affecting only the third significant digit. The same holds for the out-of-plane Young's modulus E 3 . The most affected is the Young's modulus in principal fiber direction. The two cell sizes whose edge length is slightly below and above a single fiber length, (192 μm) 3 and (256 μm) 3 , show the highest relative error of about 1.4% and 1.1%, respectively. The order of magnitude of this error is on the same level as the orthotropic approximation error err orth and the resolution error studied in the previous section. In particular, we come to the surprising conclusion that the smallest considered cell is representative already for typical engineering accuracy. Considering even smaller volume elements turned out to be difficult as only 64 fibers are contained in a cell of size (192 μm) 3 . Let us emphasize that the smallness of the necessary RVEs is certainly the result of two aspects. For a start, we consider periodic boundary conditions for periodic microstructures. It is well-known in the engineering and in the mathematical community that periodized ensembles lead to reduced representativity errors. 39,40,80 The second ingredient rests upon the fact that we enforce a number of factors which are known to influence the effective properties, for instance the fiber-volume content and the fiber-orientation tensor A. In a way, such a strategy may be considered as a variance-reduction strategy for stochastic homogenization. 86

Straight versus curved fibers
In this section, we wish to compare the effective properties of straight and curved short fibers. For this purpose, we study varying fiber orientation. More precisely, we parameterize the second-order fiber-orientation tensor by a parameter t ∈ [0, 1]. Via this parameterization, all transversely isotropic fiber-orientation states which may be represented via the chosen closure approximation are covered. For t = 0, a planar isotropic state in the y-z-plane emerges. For t = 1∕3, we obtain an isotropic fiber orientation, whereas the case of unidirectional, that is, aligned, fibers is represented by t = 1. All the other parameters were chosen to coincide with those introduced in Section 3.1.
We choose a resolution of h = 2 μm and a cell size of (320 μm) 3 , which turned out to be confidently representative in the previous section.
We compare volume elements with straight fibers and with curved fibers. Figure 7 shows a realization in either case for varying fiber orientation, parameterized via Equation (23). Please notice that for the unidirectional state, the potentially curved fibers need to be straight, as well. Indeed, suppose the second-order fiber-orientation tensor A has the form A = p ⊗ p for a unit vector p. Let q be a unit vector which is orthogonal to p. Then, it must hold with the second-order fiber-orientation tensor compare Equation (9). As all the summands are non-negative, we conclude that is, each segment is orthogonal to the direction q. As this equation holds for all directions q orthogonal to the unit direction p we started out with, each segment needs to be parallel to the direction p, that is Thus, for a unidirectional fiber-orientation state, all fibers need to be straight and aligned.
Let us return to the study at hand. We generated ten realizations for both straight and curved fibers associated to the parameter values We chose t slightly larger than zero to avoid geometric difficulties with strictly planar orientations, see Schneider 38( §3.2.2) for a discussion. We also selected a value of the parameter t slightly below one to see whether there are differences between straight and curved fibers close to the unidirectional state.
The engineering constants 77 of the effective transversely isotropic material are shown in Figure 8A. The transversely isotropic approximation error is close to 1% for all considered cases, in agreement with the observations in the previous section made for the orthotropic approximation error. We generated ten samples for each configuration to assess the dispersion of the results. However, it turned out that the dispersion is rather small. The largest absolute dispersion is realized for curved fibers with t = 0.98 and the longitudinal Young's modulus E L , and amounts to 0.054 GPa. In terms of the corresponding absolute value of E L = 11.986 GPa, the relative dispersion is 0.45%. In particular, there is a high confidence in the computed results, and we will only discuss the mean values of the effective properties computed for the ten realizations. We observe in Figure 8A that the various Young's moduli and the shear moduli almost coincide for straight and curved fibers. The differences between all the considered moduli is about 1%. In particular, from an engineering point of view, this deviation may be neglected. Our observations confirm those of Altendorf et al. 37 who rely upon a slightly different microstructure-generation strategy for curved fibers.
The observation may also be interpreted differently and used as a(n a posteriori) justification for the various computational studies of short-fiber reinforced composites conducted with straight fibers.
Last but not least let us look at the runtime of the introduced algorithm for generating fiber-filled volumes, see Figure 8B, for identical parameter setups and the same (Python + Cython) implementation. We observe that generating straight fibers leads to a smaller standard deviation than generating curved fibers. Moreover, except for the borderline cases t = 0.02 and t = 0.98, the runtimes are virtually independent of t and are almost equal for straight and curved fibers. The average runtime for curved fibers is slightly higher. However, for these cases the runtime never exceeds 7 s, which we consider acceptable for applications. For the almost aligned case t = 0.98, generating curved fibers takes 11.6 s on average, which is about three times higher than for the straight case. For the almost planar state, t = 0.98, generating curved fibers takes 51.4 s, which is about six times higher than for straight fibers. Still, such runtimes are acceptable taking into account the runtime for the computational homogenization.

Long fibers
The previous section showed that it is not necessary to consider fiber bending when considering the effective elastic properties of typical short-fiber reinforced composites. The situation changes, however, when longer fibers need to be considered. To illustrate the basic problem with straight fibers, let us consider a volume fraction of = 13% and a fiber orientation which are typical for LFT, 35,88 and take a look at Figure 9. When generating straight fibers, the minimum edge length of the volume elements should exceed the fiber length. In principle, the SAM algorithm can handle fibers which are longer than the edge length of the volume element. However, whenever such a fiber is displaced, it will be displaced everywhere, leading to a strongly nonlocal interaction of the fiber with its environment whenever the fiber wraps around the cell multiple times. Another consequence of considering straight fibers is the emergence of bundles at high volume fractions or for long fibers, see Figure 9D. Although fiber bundles are commonly observed in digital volume images of LFTs, within the SAM framework for straight fibers the emergence of bundles is a consequence of the algorithm, and is not associated to any modeling. Moreover, the formation of individual layers with roughly aligned fibers is observed. In any case, when increasing the fiber length, the considered volume elements grow quickly, that is, cubically in the fiber aspect-ratio. In particular, extremely high aspect ratios cannot be handled in this framework. Moreover, LFTs are often used for components with fine features and thin parts. The edge length of large cubic volume elements may then exceed the cross section of the component, which is not realistic.
For curved fibers, the situation is entirely different. Due to the ability to bend, an overlap does not necessitate a global displacement or rotation. Rather, a local change may be sufficient to remove the overlap. In this way, volume elements may be considered which are significantly shorter than the considered fibers. To illustrate the matter, Figure 9A shows a volume element with fibers which have an aspect ratio of 150. This is significantly larger than the longest straight fibers considered, see Figure 9D.
In this section, we fix the volume-element size 600 μm × 600 μm × 256 μm, as shown in Figure 9A, and increase the fiber length.
Realizations of such volume elements with increasing aspect ratio are shown in Figure 10. Up to an aspect ratio of 150, the volume elements could be generated routinely at once, that is, inserting all fibers as a single portion into the volume element and running the overlap removal algorithm (16). For higher aspect ratio, this strategy turned out to be not (always) successful. Thus, we added the fibers in 0.5% steps, which worked also for higher aspect ratios.
Visually, the shown volume elements do not differ much, probably due to the identical size. To illustrate the differences, we highlighted a single fiber for each configuration. For an aspect ratio of r a = 30, the fiber length equals half of the cell edge-length. With increasing aspect ratio, the fibers wrap around the volume element several times. For the highest aspect ratio considered here, each fiber wraps around the volume element four times (on average, and only in plane). Please notice that, for the different aspect ratios, the number of segments is always the same. Only the number of fibers is different. For r a = 30, 508 fibers are contained in the element. For r a = 240, only 63 fibers are necessary.
As the fiber-orientation tensor (24) is orthotropic (the most general case which can be realized for second-order tensors), we expect the effective elastic properties to be orthotropic as well. In particular, we are interested in the orthotropic engineering constants. For all considered cases, both curved and straight, the orthotropic approximation error does not exceed 0.7%.
The effective directional Young's moduli are shown in Figure 11A. We observe the ordering  For the Young's modulus E 1 , the differences are larger. For an aspect ratio of 30, the Young's modulus is 7.1% lower than for r a = 240. These errors decrease to 2.6% and 1.2% for aspect ratios of 60 and 90, respectively. For the shear moduli, shown in Figure 11B, the situation is similar. The shear moduli involving the z-direction, G 13 and G 13 remain almost unaffected by the increase in the aspect ratio. Only the in-plane shear modulus G 12 is affected strongly up to an aspect ratio of 90, where the relative deviation is only 0.53%.
The results are based on ten realizations for each aspect ratio except for r a = 240, where only three realizations were considered. The dispersion is low, in general. The highest absolute standard deviation is observed for r a = 210 and E 1 with a value of 0.013 GPa, which corresponds to a relative value of 0.18%. To get an assessment of the bias, we generated volume elements with a size of 1200 μm × 1200 μm × 512 μm, that is, doubling the edge length for each coordinate direction, for the aspect ratios r a = 60, 120, 180, 240, based on three realizations. The corresponding (mean) apparent properties are included in Figure 11 by + marks. We observe negligible bias. Thus, we were able to present representative results for volume elements which are significantly smaller than the fiber length. In particular, this novel technology opens the door to treating composites with long fibers, like LFT or sheet molding compound (SMC) composites. 89 Moreover, we wish to compare the results obtained for curved fibers to the effective properties computed for microstructures with straight fibers, see Figure 9B-D. Due to the drastic increase in (necessary) volume-element size, only aspect ratios up to 90 could be computed in reasonable time. Of course, advanced multi-processing implementations 90,91 could offer mitigation. Still going to rather high aspect ratios requires excessive computing resources. For the straight-fiber geometries, the orthotropic approximation error was also small, around 0.7%. For the considered engineering constants, see Figure 11, we observe that the Young's moduli E 1 and E 2 as well as the shear modulus G 12 are slightly higher than for the curved fibers. The shear modulus G 13 is slightly smaller, and the moduli E 3 and G 23 are the same in engineering precision.
Inspecting the generated geometries in Figure 9 clarifies the differences. For straight fibers, the fibers prefer to form bundles and to arrange in layers. Both strategies tend to increase the stiffness in the direction of the bundles, as the load can be transferred more efficiently. By the way, for an aspect ratio of r a = 90, the dispersion increased suddenly for the Young's modulus E 1 with a mean of 7.253 GPa and a standard deviation of 106 MPa, which is about 1.46%. Thus, compared to the case of curved fibers, the dispersion is about one order of magnitude higher. In particular, it might be necessary to increase the volume-element size to obtain representative results. This observation adds further weight to considering curved fibers for modeling fiber composites with rather long fibers.
Last but not least, we would like to study how the generated microstructures compare to real ones in terms of the effective properties. For this purpose, we furnish the matrix material with the isotropic elastic properties of polypropylene with Young's modulus E = 1.8 GPa and Poisson's ratio = 0.35. 35,88 Computations of the effective elastic properties on the samples with fiber aspect-ratio r a = 240 yields the directional Young's moduli E 1 = 6.22 GPa, E 2 = 3.28 GPa, and E 3 = 2.71 GPa.
Compared to the experimental results, 35 E 1 = 6.65 GPa and E 2 = 3.07 GPa, we observe that the first Young's modulus is underestimated slightly whereas the second Young's modulus is overestimated slightly. The differences, however, are still within the experimental standard deviation. Moreover, there is also some uncertainty concerning the out-of-plane components of the fiber-orientation tensor (24).

CONCLUSION
This work was devoted to an extension of the SAM algorithm 38 to curved fibers described as chains of straight fiber segments. We showed that for short fibers, considering the fiber curvature is negligible, confirming previously reported results 37 for a different type of fiber-generation algorithm. We demonstrated the ability of the algorithm to generate periodic microstructures with fibers that were much larger than the edge length of the considered volume elements. Moreover, the generated volume elements show very little dispersion and bias. We speculate that this is a result of the flexibility of the fibers. Indeed, the considered model permits to resolve a potential overlap between fibers by a local move, whereas a global movement of the fiber is necessary for straight fibers. In particular, the correlation length (or integral range 92,93 ) is much smaller for our model of curved fibers than for straight fibers. As the integral range is closely tied to the RVE size, 80,94 this coincidence explains the comparatively small volume elements which emerged as representative.
Still, we observed a high number of iterations (on the order of 10 6 or 10 7 ) for high aspect ratios and/or volume fractions, or even the inability to generate suitable microstructures in reasonable time. It remains to be identified whether this is a restriction of the model or whether a more clever algorithmic realization may mitigate these shortcomings.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.