Random Fiber Array Generation Considering Actual Noncircular Fibers with a Particle-Shape Library

: In this work, we generated a set of random representative volume elements (RVEs) of unidirectional composites considering actual noncircular cross-sections and positions of ﬁbers with the aid of a shape-library approach. The cross-section of the noncircular carbon ﬁber was extracted from the M55J / M18 composite using image processing and a signed-distance-based mesh trimming scheme, and they were stored in a particle-shape library. The obtained noncircular ﬁbers randomly chosen from the particle-shape library were applied to random ﬁber array generation algorithms to generate RVEs of various ﬁber volume fractions. To check the randomness of the proposed RVEs, we calculated spatial and physical metrics, and concluded that the proposed method is su ﬃ ciently random. Furthermore, to compare the e ﬀ ective elastic properties and the maximum von Mises stress in the matrix, it was applied to composite materials with di ﬀ erent relative ratios of elastic moduli of M55J / M18 and T300 / PR319. In the case of T300 / PR319 having a high RR T (relative ratio of the transverse elastic moduli), simulation results were deviated up to about 5% in the e ﬀ ective elastic properties and 13% in the maximum von Mises stress in the matrix according to the ﬁber shapes.


Introduction
Random fiber generators can provide arbitrary arrangements of fibers or fillers in a two-dimensional (2D) or three-dimensional (3D) space, which approximate the shape of the microstructure of a composite material. As a result, the representative volume elements (RVEs) simulating the real microstructure can be constructed; they have been commonly used in studying the mechanical behavior of composite materials and structures. RVEs can be an excellent micromechanical-based model that predicts equivalent physical properties, such as the strength, stiffness, coefficient of thermal expansion, and so forth of composites. However, for efficient modeling of these RVE models, there are some situations in which the shape of the fiber/filler is excessively simplified into a circle, sphere, cylinder, or disc. Sometimes, these approximations could distort the effect of stress concentration between fibers, resulting in less accurate prediction of the equivalent stiffness and strength. In addition, the actual fiber has a more complex shape due to various manufacturing defects such as waviness [1,2], misalignment [3,4], crater, and so forth [5][6][7]. In particular, in the case of unidirectional composites, unlike the cross-sections of glass fibers, which are almost circular, those of carbon fibers have been reported to have various shapes, such as circular, triangular, C-shaped, kidney-shaped, and so forth, so that the stress concentration between the fibers is complicated [8][9][10][11][12][13][14][15][16].
shape library. Furthermore, it can be seen that the actual fibers of the M55J/M18 composite are concave and convex.  From the final FE model, a set of the outer nodes of the intact fibers was stored in the particle-shape library of fibers as shown in Figure 2. In this work, we stored 482 fiber samples in the particle-shape library. Furthermore, it can be seen that the actual fibers of the M55J/M18 composite are concave and convex.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 21 shape library. Furthermore, it can be seen that the actual fibers of the M55J/M18 composite are concave and convex.

RVE Generation with Diverse V f
Many RVE samples were generated with various V f values of 60%, 55%, 45%, 35%, 25%, 15%, and 5% by randomly choosing the geometry of actual fibers from the particle-shape library created as described in the previous subsection. To arrange the fibers properly in a square domain, the RSE and RFR algorithms were employed with proper modification as described in the following subsection.

RVE Generation Using the RSE Algorithm
The procedure for applying the RSE algorithm [23] and its modification is introduced as follows (see Figure 3). An arbitrary rotation procedure of fibers is added to the beginning of the original RSE algorithm. A flowchart of the algorithm is shown in Figure 4.

1.
After a fiber shape is randomly chosen from the particle-shape library, it is rotated at arbitrary angles around its centroid to give a random orientation to the fiber. Here, n f denotes the current number of fibers used in RVE.

2.
For the first fiber (n f = 1), it is placed in any position in the square window of L × L (see Figure 3a). The first fiber is set as the reference fiber (F re f ).

3.
Then, step 1 is repeated to choose the new fiber to be placed from the particle-shape library.
The new fiber is placed to satisfy the random minimum distance d (l min ≤ d ≤ l max ) and random orientation angle θ (0 ≤ θ < 2π) between the new fiber and F re f . Additionally, it is necessary to ensure that the minimum distance between the new fiber and the existing fibers is at least l min . The values l min and l max can be determined according to the requested fiber volume fraction (V req f ). If the new fiber is located across the boundary of the window, the fiber is placed on the opposite side of the window for geometric periodicity (see Figure 3b). 4.
Step 3 is repeated until there is no more space around F re f for new fibers (see Figure 3c,d). To determine whether there is space for new fibers, the number of attempts to place new fibers around the F re f was set to 300. At the end of repeating step 3, a new F re f is chosen by the fiber placed next to the original F re f .

5.
Steps 3 and 4 are repeated with respect to the new F re f .

6.
This process is repeated until the current V f (= V cur f ) is close to the V req f or there is no room for a new fiber in the window.     Diverse fiber arrangements were made to have V f = 60%, 55%, 45%, 35%, 25%, 15%, and 5% in the square window of 10 × 10, and 100 samples were generated for each V f . To avoid clustering of fibers, which routinely occurs with the RSE algorithm, variables l min and l max for the minimum distance between the fibers were determined through extensive trial and error as summarized in Table 1.
Here, r f denotes the effective radius of the average area of the fibers, and it was 0.4549 for this work. Note that while an RVE including 30 fibers with V f = 50% has been proved to be adequate to represent the microstructure of unidirectional composites [26], the size of RVEs and the number of fibers were set to contain about 78 fibers for V f = 50% in this work.
After fiber placement, the parts of the fiber outside the window were trimmed and a 2D FE mesh was constructed using three-noded triangular elements assuming perfect bonding between fiber and matrix. At this time, the FE mesh was generated to ensure convergence through a mesh refinement study. Then, a 3D FE mesh consisting of six-noded prism elements was constructed by dragging them by 0.16 with two element layers in the longitudinal direction (1-direction), where 2-and 3-directions are the transverse direction.

RVE Generation Using RFR Algorithm
To apply the RFR algorithm with various V f , a master RVE with the highest V f is required. For this purpose, we used the RSE algorithm to generate the master RVE with V f = 65.57% containing 102 fibers. At this time, the parameters for the minimum distance between the fibers were set to l min = 0.07r f and l max = 0.08r f . The 3D FE model of the master RVE finally was created similarly as mentioned in Section 2.2.1. The total number of nodes and elements were 46,047 and 60,668, respectively. The procedure of applying the RFR algorithm to this master RVE is summarized as follows.

1.
First, the value of V req f is set.

2.
One of the fibers is randomly selected to be removed, and its fiber volume fraction is calculated. 3.
The elements of the fiber selected in step 2 are replaced by those of the matrix. 4.

Check the Randomness of the Centroids with Statistical Spatial Metrics
To check the randomness of the RVEs generated by the proposed method, we calculated the spatial metrics [27], such as the nearest neighbor orientation, Ripley's K function, and pair distribution function, using the centroids of fibers of the RVEs, and compared them with those of the completely spatial random (CSR). Furthermore, as seen in Figure 1, experimental data of fiber positions were extracted by placing a square window at the random position in the microscopic image of M55J/M18 to compare them with those of other generated RVEs.

Nearest Neighbor Orientation
Nearest neighbor orientation is a cumulative distribution function for the orientation angle of the line between each fiber and its nearest fiber [24]. Figure 5a shows the nearest neighbor orientation of experimental data, RFR, RSE, and CSR, according to V f = 65% and 45%. To compute the mean value (µ) and error bar of the nearest neighbor orientation according to the angle, 10 experimental samples with V f = 65% and 100 RVE samples with V f = 45% for each case of RFR and RSE were prepared (note that for V f = 65%, one master RVE with noncircular fibers was used). The results show that the nearest neighbor orientation of experimental data, RFR, and RSE is close to the CSR pattern. Furthermore, the RSE results are closer to experimental data and CSR than those of RFR. This finding is consistent for the case of circular fibers as reported in other works [23,24].

Comparison of the Results between Actual Noncircular Fibers and Circular Fibers
In this section, we compared the performance of the proposed RVEs of actual noncircular fibers with conventional RVEs of circular fibers. To generate RVEs consisting of circular fibers, the RFR and RSE algorithms were also used. The minimum distance between fibers used in the RSE algorithm was the same as that of noncircular fibers. Furthermore, the radius of circular fibers was set to f r .

Ripley's K Function
Ripley's K function [23], also termed the second-order intensity function, is defined as the ratio of points within an arbitrary radius distance r to the unit square area in Equation (1) as where A is the area of the window, N is the number of all points in the window, and d ij is the distance between points i and j. Additionally, I( ) means an indicator function with a value of 1 if the expression in parentheses is true, and a value of 0 otherwise. Here, w(i, j) is a weight function that represents the ratio of the circumference inside the window to the total circumference of the circle passing through point j, centered on point i. Ripley's K function in CSR is expressed in Equation (2).
If the K(r) corresponding to any of the points is higher than the K r (r) of the CSR pattern, it means that the points are clustered. On the other hand, being lower than the CSR pattern indicates that the points are somewhat regular [23]. Figure 5b shows Ripley's K function of experimental data, RFR, RSE, and CSR according to V f = 65% and 45%. Here, we divided r by the r f for normalization. It can be seen that the results for experimental data, RFR, and RSE are close to the CSR pattern.

Pair Distribution Function
The pair distribution function, termed radial distribution function, indicates the probability that points will exist in an annular region with an inner diameter r and an outer diameter r + dr; it can be written as Substituting Equation (2) into Equation (3), the pair distribution function of CSR is g(r) = 1. Therefore, as r increases, the pair distribution function for the fiber distribution approaches 1, indicating that the fiber distribution is randomly arranged [28]. Figure 5c shows the pair distribution function of experimental data, RFR, RSE, and CSR according to V f = 65% and 45%. The results of the experimental data, RFR, and RSE show that the larger r is close to g(r) = 1, the CSR pattern. In particular, in the case of V f = 65%, the results of RSE were almost same to the experimental results. It can be seen that the proposed RVEs are proper to simulate the cross-section of actual fibers.

Comparison of the Results between Actual Noncircular Fibers and Circular Fibers
In this section, we compared the performance of the proposed RVEs of actual noncircular fibers with conventional RVEs of circular fibers. To generate RVEs consisting of circular fibers, the RFR and RSE algorithms were also used. The minimum distance between fibers used in the RSE algorithm was the same as that of noncircular fibers. Furthermore, the radius of circular fibers was set to r f . Figure 6 shows the RVE samples of V f = 45% according to RFR and RSE algorithms with conventional circular fibers and actual noncircular fibers by a shape-library approach. Table 2 shows the µ, standard deviation (SD, σ), and relative standard deviation (RSD, σ/µ) of V f of the RVE samples corresponding to V req f , and it can be seen that V req f and its corresponding V f of RVE samples are not exactly identical. In particular, the RVE samples with noncircular fibers show some deviation from those with circular fibers because they have large RSD values in comparison to the RSD values of circular fibers. Given that the volume of each noncircular fiber is not identical, RVE samples with the same V f cannot be produced. The V req f will be approximately satisfied. Figure 6 shows the RVE samples of f V = 45% according to RFR and RSE algorithms with conventional circular fibers and actual noncircular fibers by a shape-library approach.     To calculate the effective elastics properties of RVEs, we employed a popular computational homogenization scheme [29] with a set of periodic boundary conditions (PBCs) showing high convergence rate [30] on RVEs. Its details are summarized in Appendix A. As a result, longitudinal modulus E 1 , transverse modulus E 2 and E 3 , out-of-plane shear modulus G 12 , in-plane shear modulus G 23 , and Poisson's ratio ν 12 and ν 23 were calculated by using Abaqus and MATLAB. Furthermore, for each RFR; RSE algorithm; and noncircular, circular fiber shape, 100 RVEs were used for each V f . As constituent materials of the composite, two composite materials: M55J/M18 [22] and T300/PR319 [31] were used. Fibers are transversely isotropic materials and matrices are isotropic materials as shown in Table 3. We chose these two materials due to the large difference of the relative ratios of elastic modulus RR between fiber and matrix defined as where RR L = {RR 1 , RR 2 } and RR T = {RR 3 , RR 4 , RR 5 }. Here, RR L indicates the relative ratio of the longitudinal elastic moduli, and RR T indicates the relative ratio of the transverse elastic moduli.  Figures 7 and 8 show the effective elastic properties (E 1 , E 2 , G 23 , G 12 ) according to V f from 5% to 60% of M55J/M18 and T300/PR319 by the computational homogenization, respectively. Figures 9  and 10 show the histogram and normal distribution in terms of µ and σ of the effective elastic properties for V f = 60% of M55J/M18 and T300/PR319, respectively. Tables 4 and 5 show the µ and relative error of the effective elastic properties for V f = 60%, 55%, and 45% of M55J/M18 and T300/PR319, respectively. From the results, it can be seen that both materials show little deviation in E 1 , but the deviation of RVE with noncircular fiber is greater than that of RVE with circular fibers as seen in Figures 9a and  10a. This is because the deviation in V f of RVE with noncircular fiber is larger than that of RVE with circular fiber as discussed in Section 3.2.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 21    Additionally, it should be noted that, in the case of M55J/M18, the deviations of the transverse moduli E 2 , G 23 are minimal regardless of the fiber shape and the algorithm of random RVEs. However, the other transverse modulus G 12 varies according to the fiber shape and the algorithm of random RVEs. In the case of T300/PR319, all transverse modulus E 2 , G 23 , and G 12 show large deviations according to the RFR and RSE algorithms, and as V f increases, the deviation increases. For example, for V f = 60% with the RFR algorithm, the deviation of G 12 according to fiber shape is 4.33%, which is larger than that of M55J/M18 (=2.58%). As mentioned earlier, it could be because RR T of the M55J/M18 composite is not significant compared to that of T300/PR319. Although the difference in the effective elastic properties according to the fiber shape of M55J/M18 is small, the actual noncircular fiber shapes should be considered because of the transverse tensile strength accompanying crack propagation. The reason is that, according to our previous work [18], it was reported that the result of the transverse crack propagation simulation showed a different crack pattern from the test result when the circular fibers were assumed instead of the actual fibers. In addition, as RR T of the composite is larger, the deviation of E 2 , G 23 , and G 12 values of RFR are larger than those of RSE. This is because the minimum distance between the fibers of the RFR is closer than that of the RSE, as reported in [24].
To check randomness using the effective elastic properties of RVEs with noncircular fibers in the transverse direction, the symmetry of the compliance matrix was checked as follows: In addition, the anisotropic ratio was confirmed as follows: Table 6 shows the results for Equations (5) and (6). It can be seen that all values are close to unity and are sufficiently random in the transverse direction.

Comparison of the Maximum Von Mises Stress in Matrix
In this section, the distribution of the von Mises stress (σ ν ) in the matrix is compared when 1% of the tensile strain is given with PBCs in the transverse direction (2-direction) by using Abaqus. At this time, the stress was calculated at the centroid of the element. Furthermore, as in Section 3.2.1, for each RFR; RSE algorithm; and noncircular, circular fiber shape, 100 RVEs were used for each V f . Figure 11 shows an example of a contour sample of von Mises stress in the matrix of T300/PR319 with V f = 60%. Figure 12 and Table 7 show the maximum von Mises stress value in terms of the µ and error bar in the matrix for diverse V f of M55J/M18 and T300/PR319, respectively. The results show that the RVE generated by RFR has a higher maximum von Mises stress value in the matrix than those generated by RSE. As mentioned in Section 3.2.1, it is because the minimum distance between the fibers of RFR is smaller than that of the RSE, so that it has larger maximum von Mises stress in the matrix. Furthermore, as the V f increases, µ of the maximum von Mises stress tends to increase. However, the variation of µ is sensitive to the minimum distance between the fibers as reported in other work [32]. In this work, we set the minimum distance of RFR as l min = 0.07r f and l max = 0.08r f and that of RSE as in Table 1. Thus, the trend of the maximum stress could be altered by adjusting the minimum distance. Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 21  In addition, as RR T increases (i.e., T300/PR319), the RVE composed of noncircular fibers tends to a larger maximum von Mises stress in the matrix and wider deviation than the RVEs composed of circular fibers. In the case of T300/PR319 with V f = 60%, the deviations of RFR and RSE according to fiber shape are 8.83%, 9.68%, respectively. Furthermore, it can be seen that the maximum deviation is 12.69% when V f = 45% with the RSE algorithm in T300/PR319. However, in the case of M55J/M18 with a lower RR T , the maximum deviation of von Mises stress is quite minimal, less than 2%. Thus, it is concluded that the effect of shape on the maximum stress is critical as the RR T is larger.

Conclusions
In this work, we generated RVEs considering actual noncircular fibers randomly selected from a particle-shape library, including the features of M55J. The RSE and RFR algorithms were employed to create RVEs various V f along with proper modification to consider the noncircular fibers.
To check the randomness of the proposed RVEs, we calculated spatial and physical metrics, concluded that the proposed method is sufficiently random, and reproduced experimental data in terms of the shapes and positions of actual fibers. In addition, to investigate the effect of the fiber shapes on the stiffness and strength of the composite, RVEs with circular fibers were also prepared. The material properties of M55J/M18 and T300/PR319 were adopted, and the effective elastic properties were compared. Based on the results, it is concluded that M55J/M18 having a low RR T showed minimal deviation according to the fiber shapes. However, in the case of T300/PR319, which has a high RR T , the maximum deviation of G 12 according to the fiber shapes was about 5%. Finally, to investigate the effect on the material fracture strength, the maximum von Mises stress in the matrix was compared with a tensile strain of 1% in the transverse direction according to the fiber shapes. The results showed that the maximum difference according to the fiber shapes of M55J/M18 was up to about 2%, but that of T300/PR319 was up to about 13%. Thus, from a practical point of view, when predicting the strength of composite using micromechanical models having a large RR T of the composites, it is necessary to consider the actual noncircular fiber shapes and positions in the micromechanics model.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Computational Homogenization Schemes
The stress and strain relationship in a unidirectional lamina is typically expressed as where σ i , ε j , and D ij are stress, strain, and stiffness matrices in the form of the vector notation, respectively, and x is the position vector at an arbitrary point in an RVE. Due to the transversely isotropic symmetry of the stiffness D ij , it is well known that there are only five independent elastic moduli: E 1 (= E L ), E 2 = E 3 (= E T ), G 12 = G 13 (= G L ), G 23 (= G T ), and ν 12 = ν 13 (= ν L ). In the case of random fiber arrays, due to an anisotropic ratio, one more parameter, ν 23 = ν 32 (= ν T ), is needed. To determine the elastic properties of unidirectional composites in a theoretical manner, computational homogenization in relation to the present study is briefly summarized.
To predict the equivalent elastic properties of a unidirectional lamina in the computational homogenization scheme, an RVE comprising fibers and a matrix with a predefined fiber-volume ratio is first modeled with a set of FE meshes. Proper periodic boundary conditions in Equation (A2) are then applied to satisfy Hill-Mandel or macro-homogeneity condition [33,34]: u r+ p − u r− p = ε pq ∆x r q , (p, q, r = 1, 2, 3) (A2) where the superscripts r+ and r− indicate the positive and negative sides of the RVE in the r-direction, respectively. That is, the displacements u p for each periodic pair are constrained by Equation (A2). In addition, ∆x r q indicates the difference between the coordinates x q for the periodic pair in the r-direction. Only in the case of r = q, ∆x r q corresponds to the length of the RVE with respect to the r-direction, i.e., ∆x 1 1 = X 1 , ∆x 2 2 = X 2 , ∆x 3 3 = X 3 , respectively; and the others are zero when the periodic boundary surfaces of the RVE are assumed to be perpendicular to the Cartesian coordinates axes. The prescribed strain imposed on the RVE is denoted by ε pq .
We consider six independent macroscopic periodic deformations of the RVE, which correspond to six strain components in the form of vector notation, respectively. First, one of the six strain components, corresponding to the prescribed ε pq , is non-zero, while the others are zero. For example, when ε 23 (= ε 32 ) is prescribed, the corresponding constraints are given as Then, the FE analysis was conducted imposing the above constraints at all the nodes on boundary surfaces of the RVE. Similarly, the FE analyses for the other five strain components were carried out. For each of the macroscopic periodic deformations, the average stresses can be obtained by < σ i >= 1 Ω Ω σ i (x)dΩ, (i = 1, 2, . . . , 6) where Ω is the RVE volume. The macroscopic strain-stress relationship in the form of the vector notation is given with the compliance matrix S ij , as follows: < ε i >= S ij < σ j >, (i, j = 1, 2, . . . , 6) (A5) where the macroscopic averaged strain < ε i > corresponds to the prescribed strain ε pq . The compliance matrix S ij can be thus computed from Equations (A5) and (A6). Finally, all elastic moduli are determined by comparing the computed compliance matrix with the compliance matrix of a unidirectional lamina as follows: