The influence of microstructure randomness on prediction of fiber properties in composites

Abstract

Elastic properties were predicted for AS4, IM7, T300, and T650-35 graphite fibers. An inverse method was employed using lamina and epoxy properties taken from literature. Fiber properties predicted using finite element analysis based on a hexagonal microstructure, finite element analysis based on a random microstructure, and the Mori-Tanaka averaging scheme were compared. It was observed that the Mori-Tanaka averaging scheme and finite element analysis based on a hexagonal microstructure predicted nearly identical fiber properties. In contrast, randomness in the arrangement of fibers results in significantly different predictions for the longitudinal shear modulus of the fiber. The three microstructural models were also used to predict properties for isotropic E-glass 21xK43. It was observed that the three models predicted nearly the same properties for the glass fibers.

Keywords

Randomness fiber properties unit cell microscale finite element analysis inverse analysis representative volume element (RVE)

Introduction

Micromechanical analysis of fibrous polymer matrix composite (PMC) systems has been a field of great interest for some time. By examining the behavior of a composite at the microscale, where fibers and matrix are accounted for discretely, it becomes possible to directly model a variety of complex phenomena in composites. Some of these phenomena include predicting the onset and progression of matrix cracking and fiber-matrix debonding under transverse load,^1–3 predicting matrix plasticity under longitudinal shear,^3–5 and characterizing the stresses and strains that arise in the fibers and matrix due to thermal loading and moisture absorption.^6–9 Microscale models also provide a means for predicting how these various phenomena will interact. Homogenization techniques can then be utilized to bring information regarding the composite’s response from the microscale up to the lamina or tow scale and subsequently up to the structural scale. This hierarchical modeling approach makes it possible to avoid the costly process of empirically characterizing a composite system’s response for a wide variety of loading and environmental parameters through direct experimentation. In many cases, such an empirical approach may even be impossible due to the high dimensionality of the parameter space.

One of the fundamental challenges associated with microscale modeling of fiber-matrix composites is the determination of the constituents’ properties. While it is possible to directly measure some of these properties (e.g. the engineering constants for the neat matrix and the longitudinal Young’s modulus of the fibers), other properties, such as engineering constants associated with transverse and shear deformation of the fibers, might not be directly measurable due to the small transverse dimension of the fibers. These properties are typically obtained through the solution of an inverse problem. A micromechanics model is utilized to determine the constituent properties that yield homogenized properties which match values measured in experiments. This procedure has been conducted for a variety of material systems using both analytical and finite element-based models^3,10–12 These models either do not consider fiber arrangement (e.g. Mori-Tanaka averaging scheme) or they are based on regular arrays of fibers. The solution of this inverse problem using a finite element micromechanics model based on a random microstructure is examined in this paper.

The introduction of randomness to the fiber positions in the microstructure is desirable for a variety of reasons. It results in a more realistic microstructure. Randomness leads to fibers in very close proximity to one another, which gives rise to high stress concentrations. These concentrations will strongly influence extreme-based mechanisms such as failure and plasticity. The introduction of randomness can also be useful for the characterization of uncertainty in composite properties. The use of random microstructures to predict the behavior of composites is not new.^13–17 However, previous work with random microstructures has typically used fiber properties obtained from models that do not include randomness, creating an inconsistency. It is possible that the use of inconsistent fiber properties in a micromechanics model may introduce inaccuracies into the analysis that affect the homogenized response as well as the micro-scale stresses in a problematic way. The current work examines this possibility.

This paper provides a description of the microstructural models used to represent the microstructure of the fiber-matrix composite. This is followed by a description of the algorithm used to determine fiber properties and a discussion of the assumptions made concerning the material properties.

A description is given of each fiber-matrix material system examined. Two of these systems have the same fiber type but different epoxy matrices, which will provide insight into the consistency of the methods described. In addition to the graphite fiber material systems, a glass fiber material system is examined to explore the effects of the modifications in the approach required to predict properties for isotropic fibers. A discussion of noted inconsistencies in literature when reporting a material system’s fiber volume fraction, fiber’s longitudinal Young’ modulus, and lamina’s longitudinal Young’s modulus is also given, addressing some of the assumptions made when selecting properties.

Since multiple realizations are required for the random microstructure, the results begin with a study of how the RVE size and number of realizations affect the predicted properties in order to determine an appropriate random RVE. The predicted graphite fiber properties using various material systems, volume fractions, and micromechanics models are presented and compared to the values given in literature. After these findings, the method is applied to an isotropic material system.

Models and approach

This section describes the models used to represent the hexagonal and random RVEs, the Mori-Tanaka averaging scheme, the algorithm developed to determine fiber properties, and the assumptions associated with the material properties.

Microstructure models

Two types of microstructures were examined using finite element analysis (FEA) – hexagonal and random fiber arrays. Quasi-three-dimensional (3D) assumptions were utilized to reduce the analysis domain to two dimensions, reducing the number of degrees of freedom and the bandwidth of the stiffness matrix. These assumptions, given in Ref. [18], are exactly satisfied by the configuration currently being studied, so they do not introduce any inaccuracy as compared to traditional 3D analysis. In addition to FEA-based microstructural models, the Mori-Tanaka averaging scheme was utilized.

Hexagonal arrangement

A single unit-cell of a hexagonal array was generated for the fiber volume fraction examined. An example of this is shown in Figure 1. Periodic boundary conditions are applied to the unit cell by imposing the following relationship for the displacements at the unit cell boundary

u_{i} (x_{j} + d_{j}) = u_{i} (x_{j}) + á \frac{\partial u_{i}}{\partial x_{k}} ñ d_{k}

(1)

Repeated indices imply summation, and $d_{k}$ is a vector of periodicity that starts at some coordinate x_j in one unit cell and ends at the equivalent point in another unit cell. Angle brackets imply the volume-averaged value of a field variable. A more in-depth description of the imposition of periodic boundary conditions is given by Appendix 1 and Ref. [19].

Random arrangement

Figure 2(a) shows a micrograph of an actual lamina, illustrating that physical laminae do not exhibit a uniform distribution of fibers but rather have some random distribution. As a result, this work aims to model the lamina microstructure more accurately by introducing randomness into the arrangement of fibers, as Figure 2(b) shows.

This work made use of a procedure to create a periodic RVE that contains randomly positioned fibers.²⁰ Random positioning of fibers is accomplished by iteratively adjusting fiber positions to resolve spatial interference. This approach results in some fibers being in extremely close proximity to one another. The algorithm used to generate these geometries proceeds in the following manner:

The fiber volume fraction (V_f), the number of fibers, and the fiber radius are specified

RVE dimensions are calculated and fibers are initially placed at random coordinates distributed uniformly throughout the unit cell

Iterations are performed to resolve spatial interference between overlapping fibers while enforcing periodicity of the geometry

For the ith fiber, overlap is checked for each fiber with index j > i

If overlap is found, the ith and jth fibers are each moved away from one another by one-half the overlap distance plus the radius of the fiber divided by 500

When no overlap is found for any fiber, then iteration is stopped

The geometry is meshed using quadratic triangles, taking care to ensure that nodes on the RVE boundaries are positioned periodically

Periodicity in the mesh is obtained by periodic mesh seeding along the edges. Periodic boundary conditions are imposed in the same manner as for the hexagonal array. A more detailed description of the boundary conditions is shown in Appendix 1. The triangle meshing tool²¹ generates straight-sided quadratic triangular elements, so the mid-side nodes along the fiber/matrix interface are moved onto the geometric interface to obtain a more circular boundary around the fibers when meshing.

A minimum spacing between fibers can be enforced by creating a random microstructure using a larger volume fraction and fiber diameter than desired (see Figure 3(a)), then shrinking each fiber’s diameter (see Figure 3(b)).

For this work, the minimum spacing fraction is defined to be the minimum distance, s, between any two fibers divided by the radius of the fiber, r, as illustrated in Figure 4. The relationship between the minimum spacing fraction, ξ, the artificially large volume fraction, $V_{f}''$ , and the final volume fraction, V_f, is given by

ξ = 2 (\sqrt{\frac{V'_{f}}{V_{f}}} - 1)

(2)

The largest value of fiber spacing for a given volume fraction is limited by the largest possible value of $V_{f}''$ , which in theory can be about 90.7%, corresponding to hexagonal close packing arrangement for circles. In practice, the algorithm for generating the random RVE’s cannot achieve the theoretical limit due in part to the fact that a periodic hexagonal arrangement cannot be represented by a square RVE. As the size of the RVE increases, the square shape has less of an effect on the maximum volume fraction. The larger the RVE, the closer the algorithm can get to the theoretical volume fraction limit.

Mori-Tanaka averaging scheme

The Mori-Tanaka averaging scheme considers inhomogeneities placed into an infinite medium of the matrix material. It accounts for interactions between phases through averaging and the geometry of the inhomogeneities (cylinders in this case) through the Eshelby tensor. However, the Mori-Tanaka averaging scheme does not account for the size and locations of the inhomogeneities.²²

Inverse method for fiber property determination

The procedure for determining the fiber properties involves the solution of an inverse problem. A micromechanics model uses the fiber and matrix properties to determine the homogenized composite’s properties. In the case of the Mori-Tanaka averaging scheme with an isotropic matrix material, the relationship can be inverted and the fiber properties can be solved directly. However, it is not possible to invert this input-output relationship for finite element models. In these cases, the problem is solved by iteratively identifying fiber properties that cause the model to yield lamina properties that match the expected values. For convenience, the unique engineering constants of the fiber and lamina properties are expressed as vectors. The ordering of the properties within the vector is not important, but if they have a different number of independent engineering constants, the problem becomes a minimization problem rather than a root finding problem, which is addressed at the end of this section. For the case of transverse isotropy of the fiber and lamina properties, the constants are expressed as vectors in the following manner for the lamina, $Γ^{L}$ , and fiber, $Γ^{f}$ .

Γ^{L} = [E_{1}^{L}, E_{2}^{L}, G_{12}^{L}, G_{23}^{L}, ν_{12}^{L}] Γ^{f} = [E_{1}^{f}, E_{2}^{f}, G_{12}^{f}, G_{23}^{f}, ν_{12}^{f}]

(3)

The micromechanics model serves as a function which operates as follows

Γ^{L} = f (Γ^{f})

(4)

In the case of a hexagonal unit cell, this function represents a single finite element model. For a random RVE, this function represents the average of multiple realizations of the microstructure. The goal is to find a collection of fiber properties such that the residual of the lamina properties decreases to within a given tolerance.

ψ_{i} = Γ_{i}^{L} - {\bar{Γ}}_{i}^{L} ≅ 0

(5)

The overbar signifies the lamina properties from experiments. Because of the large differences in the magnitudes of the properties, the residuals are defined in terms of normalized properties, which are the properties divided by the initial guess of the fiber properties.

{\overset{\land}{ψ}}_{i} = \frac{ψ_{i}}{Γ_{i}^{f^{0}}} = \frac{Γ_{i}^{L} - {\bar{Γ}}_{i}^{L}}{Γ_{i}^{f^{0}}} (nosummationon i)

(6)

Similarly, the fiber properties are normalized by the initial fiber guess as follows

{\overset{\land}{Γ}}_{i}^{f} = \frac{Γ_{i}^{f}}{Γ_{i}^{f^{0}}} (nosummationon i)

(7)

The Newton-Raphson method is used to find the value of the fiber properties that result in the correct lamina properties. To avoid numerical problems, the Jacobian matrix is in terms of the normalized fiber and lamina properties.

J_{ij} = \frac{\partial {\overset{\land}{ψ}}_{i}}{\partial {\overset{\land}{Γ}}_{j}^{f}} ≅ \frac{Δ {\overset{\land}{ψ}}_{i}}{Δ {\overset{\land}{Γ}}_{j}^{f}}

(8)

The initial guess is obtained from the literature, similar material systems, or intuition, if no better source is available. It is only necessary that each engineering constant in the guess have a realistic order of magnitude, but accurate guesses can significantly help the time to convergence. For each Newton-Raphson iteration, a better approximation for the fiber properties is obtained by calculating the correction term, $Δ {\overset{\land}{Γ}}_{i}^{f}$ , as shown in equation (9).

(\frac{\partial {\overset{\land}{ψ}}_{i}}{\partial {\overset{\land}{Γ}}_{j}^{f}}) Δ {\overset{\land}{Γ}}_{i}^{f} = - {\overset{\land}{ψ}}_{j}

(9)

The correction term $Δ {\overset{\land}{Γ}}_{i}^{f}$ must be calculated for each iteration. When the number of fiber constants, $Γ^{f}$ , matches the number of lamina constants, $Γ^{L}$ , the calculation of the correction term in equation (9) requires the solution of a set of linear equations. However, when the number of fiber constants does not match the number of lamina constants, equation (9) requires a least squares solve for the correction term that minimizes the residual of the system of equations. Iteration on the fiber properties continue until the L₂ norm of the normalized residual decreases to a value of 1.0e-6 or less, at which point the fiber properties have been determined.

Several steps are taken to help expedite the solution process. The first is the use of the modified Newton’s method. The calculation of the gradient matrix is computationally expensive, requiring lamina properties to be determined for multiple different sets of fiber properties. One set corresponds to the current approximation of the fiber properties, while the other sets of properties perturb one of the fiber properties by 1% to obtain gradients. Due to this high cost, the gradient matrix is not recalculated for every iteration. Instead, it is reused as long as an acceptable rate of convergence is obtained (50% reduction in residual on each iteration). Also, it should be noted that some iterations lead to fiber properties that are not physically realistic. Therefore, the allowable fiber properties are bounded such that moduli remain positive and Poisson ratios remain between 0.01 and 1.0. If an increment results in a fiber property that lies outside these bounds, the property is set to the respective limit for that iteration. Figure 5 illustrates the described inverse method. The method to calculate the lamina properties depends on the microstructure used in the analysis. The methods for the microstructures considered in the work are shown in Figure 5. For the random RVEs, many realizations are used, so averaging is required. For the hexagonal RVE, only one solution is needed.

Each iteration involves the solution of one homogenization problem for each realization when the gradient matrix does not need to be updated or multiple homogenization problems for each realization when the gradient matrix does need to be updated. Consequently, it is common that more than a thousand FEA problems must be solved to obtain a single set of converged fiber properties. Due to this expense for the random case, the solution is expedited by initially calculating the effective properties for a small RVE size (two fibers for this work). This result is then used as the initial guess for larger RVE sizes, since the small RVE tends to give a fair prediction for the fiber properties and is much faster to run.

Investigations showed that this approach yields the same solution for a variety of initial fiber property guesses. Furthermore, consistency was tested in the following manner. Lamina properties were obtained for an arbitrary set of fiber properties. Then, the inverse method was applied using those lamina properties, and the process converged to the original fiber properties as expected.

Material properties

In this work, the properties of the unidirectional composites are assumed to be transversely isotropic. For the graphite fiber material systems, the fiber properties are assumed to be transversely isotropic, and for the glass fiber material system, the fiber properties are assumed to be isotropic.

For random microstructures with the RVE sizes currently being examined, a single realization is not expected to exhibit transverse isotropy or even orthotropy. In general, the anisotropic compliance matrix will follow the form as given in equation (10). For fully anisotropic materials, there are 21 unique constants that characterize the constitutive law.

S = [S_{11} S_{12} S_{13} S_{14} S_{15} S_{16} S_{12} S_{22} S_{23} S_{24} S_{25} S_{26} S_{13} S_{23} S_{33} S_{34} S_{35} S_{36} S_{14} S_{24} S_{34} S_{44} S_{45} S_{46} S_{15} S_{25} S_{35} S_{45} S_{55} S_{56} S_{16} S_{26} S_{36} S_{46} S_{56} S_{66}]

(10)

For orthotropic properties, the compliance matrix will follow the form as given in equation (11), for which there are nine unique engineering constants. To enforce the assumption or orthotropy, the non-orthotropic terms must be eliminated, and there are two methods to accomplish this. The effects of the two methods on the results will be compared.

S = [S_{11} S_{12} S_{13} 000 S_{12} S_{22} S_{23} 000 S_{13} S_{23} S_{33} 000000 S_{44} 000000 S_{55} 000000 S_{66}]

(11)

One method to eliminate non-orthotropic terms (the non-zero terms in equation (10) that are zero in equation (11)) assumes that for an applied normal strain, the shear stresses are assumed to be zero, and for an applied shear strain, all stresses are assumed to be zero except the respective shear stress. This corresponds to setting the non-orthotropic terms in the compliance matrix of equation (10) to zero.

The other method to eliminate non-orthotropic terms assumes that for an applied normal stress, the shear strains are assumed to be zero, and for an applied shear stress, all strains are assumed to be zero except the respective shear strain. This corresponds to eliminating the terms in the stiffness matrix rather than the compliance matrix.

For transversely isotropic properties, there are five unique engineering constants. Though the homogenized properties for a single realization are not expected to match the form for transversely isotropic materials, the unidirectional fiber matrix material properties are assumed to be transversely isotropic in this study. Therefore, in the random RVE case, the following average predicted lamina properties from all the realizations are averaged to enforce transverse isotropic: E₂ and E₃, $E_{12}$ and $E_{13}$ , and $v_{12}$ and $v_{13}$ . The percent differences between these properties that should be equal are used to measure how transversely isotropic the predicted lamina properties are and choose the RVE size and number of realizations.

Composite systems

This section describes the fiber types to be analyzed and the material systems used in the determination of fiber properties. It also discusses the inconsistencies in the properties reported by the literature for the material systems and fibers.

Systems considered

This study considers the graphite fibers of IM7, T650-35, T300, and AS4. Each of these are common graphite fibers used in fiber matrix composites for which measured lamina properties are given in the literature. The material systems used to predict the properties of the graphite fibers are IM7/8552, IM7/5250-4, T650-35/PMR-15, T300/BSL914C, and AS4/3501-6, shown in Table 1. In addition, a study of E-glass 21xK43/LY556 will be used to explore the fiber property prediction for isotropic fibers.

Table 1.

Lamina and matrix elastic properties taken from literature.

IM7/5250-4 Ref. [11]	IM7/8552 Ref. [23]	T650-35/PMR-15 Ref. [10]	T300/BSL914C Ref. [24]	AS4/3501-6 Ref. [24]	E-glass 21xK43/LY556 Ref. [24]
Rep. V_f	55.60%	60.00%	60.00%	62.00%
Calc. V_f	59.45%	59.09%	51.34%	59.29%	53.70%	65.40%
Lamina
E₁ (GPa)	165.5	165.0	126.7	138.0	126.0	53.48
E₂ (GPa)	10.34	11.38	8.255	11.00	11.00	17.70
G₁₂ (GPa)	5.792	5.120	4.200	5.500	6.600	5.830
G₂₃ (GPa)	3.315^a	3.826^a	2.810	3.929^a	3.929^a	6.321^a
ν₁₂	0.3100	0.3000	0.4100	0.2800	0.2800	0.2780
ν₂₃	0.5600	0.4870	0.4689^a	0.4000	0.4000	0.4000
Matrix
E	3.450	4.670	4.000	4.000	4.200	3.350
G	1.278^a	1.730^a	1.460	1.481	1.567	1.241
ν	0.3500	0.3500	0.3699^a	0.3500	0.3400	0.3500

Calculated using isotropic or transversely isotropic relation.

Inconsistencies in literature

For unidirectional composites, it is widely held that a rule of mixtures approximation using the longitudinal Young’s modulus of the lamina and matrix and the composite volume fraction provides a very good approximation for the fiber longitudinal Young’s modulus. Less than 0.1% difference was measured between the rule of mixtures and FEA for the material systems considered. Using a rule of mixtures approximation, the fiber volume fractions reported by some literature sources are inconsistent with the longitudinal Young’s modulus of the fibers used as reported by the manufacturers. This inconsistency leads to a dilemma between two possible approaches. One option is to take the reported fiber volume fraction, lamina properties, and matrix properties as consistent and allow the longitudinal Young’s modulus of the fiber to differ from the manufacturer’s reported value. The other is to take the lamina properties, matrix properties, and manufacture’s longitudinal Young’s modulus of the fiber as consistent and calculate the fiber volume fraction of the composite, allowing it to differ from the reported values in literature. Both fiber volume fractions will be considered in the paper allowing insight into the influence of these inconsistencies on the resulting fiber properties.

Results

The results are presented in four sections. The effect of assuming the non-orthotropic terms to be zero on the compliance matrix and strain field is explored. Next, the choice of an RVE size and number of realizations is made based on the variation of the averaged lamina properties, the transverse isotropy of the average compliance matrix, and the standard deviation of the predicted lamina properties. The predicted properties for the graphite fibers are then presented. Finally, the results for a glass fiber material system are presented.

Effect of non-orthotropic terms

For the case of the random RVE, both methods of eliminating the non-orthotropic terms of the averaged compliance matrices result in approximately the same properties with every term in the resulting compliance matrices varying by less than 0.9%. Furthermore, the non-orthotropic terms of the averaged compliance matrices have little influence on the strain states. This was shown by multiplying the full compliance matrix with a stress state to obtain the strain and by multiplying the orthotropic compliance matrix with the same stress state to obtain the approximate strain. The strains differed by less than 0.2% for the compliance matrices determined in this study. Since the method has little effect on the results, either method is valid for the RVE sizes considered. However, it is possible that very small RVE sizes (e.g. less than five fibers) may result in significant differences. For this work, the non-orthotropic terms will be eliminated from the averaged compliance matrix due to the convenience of implementation.

Determination of RVE size and realization count

The distribution of each predicted lamina property depends on the RVE size and number of realizations. Consequently, before solving for the predicted fiber properties, an appropriate size for the randomly arranged RVEs and the number of realizations were determined. The choice of an appropriate size and number of realizations was based on the average lamina properties, a measure of the transverse isotropy, and the standard deviation of the lamina properties. The T650-35/PMR-15 material system was used in the study to determine the RVE size and realization count.

Average lamina properties

To accurately characterize the fiber properties, the average of the predicted lamina properties from all the realizations must be determined within an acceptable tolerance. For various RVE sizes, the average lamina properties were determined as a function of the number of realizations. Figure 6 shows a typical graph of the average predicted lamina properties as a function of the number of realizations for a 30-fiber RVE. The 60 and 100 fiber RVEs demonstrated similar behavior and are not shown. The average lamina properties were calculated using up to about 900 realizations for the 30-fiber RVE, 600 for the 60-fiber RVE, and 450 for the 100-fiber RVE. Fewer realizations are used for larger RVE sizes since it is expected that as the RVE size is increased, a smaller number of realizations will be required to accurately characterize the average value. The value predicted by the maximum realization count for a given RVE size is known as the reference value for this work. Comparing the reference values for different RVE sizes shows that the average predicted properties vary by less than 0.5% for the RVE sizes considered when many realizations are used. This indicates that the average lamina properties are not sensitive to RVE size for RVEs with at least 30 fibers. However, this may not be the case for very small RVE sizes (such as RVEs containing two or three fibers).

The particular realizations depend on the seed used for the random number generator that determined the initial placement of the fibers. To check whether the predicted average properties were overly sensitive to the particular seeding, the process was repeated for the 30-fiber case with a different seed for the random number generator, and the results did not change significantly.

Figure 7 quantifies the variation of the average lamina properties for the 30-fiber RVE case shown in Figure 6. The averages shown in Figure 6 do not converge monotonically, which complicates the quantification of error in terms of the realization count. Instead of stating an error for a given realization count, it is better to identify the maximum error that occurs between N realizations and the maximum number of realizations. The error in this context was defined as the percent difference between the value predicted using a given number of realizations and the reference value. For an RVE size of 30 fibers, the average lamina properties do not vary more than 2% when at least 50 realizations are used. For an RVE size of 60 and 100 fibers, the variation was less than 2% for smaller realization counts (e.g. 30–50 realizations).

Transverse isotropy of lamina properties

The predicted lamina properties are assumed to be transversely isotropic. For a material to be considered transversely isotropic, the properties must not change when the compliance matrix is rotated about a single axis. In this case, the assumption means that transverse Young’s moduli (E₂ and E₃) should be equal, the longitudinal shear moduli ( $G_{12}$ and $G_{13}$ ) should be equal, and the longitudinal Poisson ratios ( $v_{12}$ and $v_{12}$ ) should be equal. In addition, $G_{23}$ should follow the relation given in equation (12).

G_{23} = \frac{E_{22}}{2 (1 + ν_{23})}

(12)

Figure 8 quantifies the percent difference between properties which should be equal if transverse isotropy holds. The difference is examined versus the number of realizations used to calculate the properties. The longitudinal shear moduli showed a much larger difference than the transverse shear moduli and longitudinal Poisson ratios. Interestingly, the percent difference between the transverse shear modulus and the value calculated using the isotropic relation (see equation (12)) did not converge towards zero as the number of realizations increased, but rather converged towards a non-zero value. However, the converged value decreased as the RVE size increased. This implies that transverse isotropy cannot be assured just by averaging many realizations but also requires a sufficiently large RVE size. The most likely reason for this is the square shape of the periodic RVE. Square arrays are fundamentally unable to yield transversely isotropic properties. The impact of the RVE’s shape can be reduced to the point of being negligible by increasing its size, but it cannot be eliminated entirely. It should be possible to obtain transversely isotropic properties with different RVE shapes (e.g. hexagons), but the development of such a model is beyond the scope of this work. For the examined RVE sizes, the deviation of $G_{23}$ from transverse isotropy is adequately small for the current investigation into fiber properties.

In order to obtain lamina properties that comply with the assumption of transverse isotropy within 1% (excluding $G_{23}$ ), about 200 realizations are required for 30- and 60-fiber RVEs, and about 50 realizations are required for 100-fiber RVEs.

Standard deviation of lamina properties

Accurately representing the distribution of lamina properties requires a sufficient number of realizations to determine both the average values of the properties as well as characterize the scatter in the properties. In the current study, the scatter is described using the standard deviation. The number of realizations to achieve an approximately constant standard deviation was shown to depend on the RVE size, with more realizations needed for the smaller RVE sizes. Similar to Figure 7, Figure 9 quantifies the maximum difference of coefficient of variation (the standard deviation normalized by the mean) that occurs between N realizations and the maximum number of realizations for the 30-fiber RVE. Comparing the data for the RVE sizes, it was shown that the standard deviation decreases as the RVE size increases.

Choice of RVE size and number of realizations

For this work, it was chosen that 200 realizations for an RVE size of 30 fibers, 100 realizations for an RVE size of 60 fibers, and 50 realizations for an RVE size of 100 fibers would be an appropriate RVE size and realization count. Since the degrees of freedom in the FEA mesh increases with the fiber count, the lower fiber counts take less time to solve. Each RVE size was timed for many realizations. The average time multiplied with the number of realizations for the respective size showed that all three combinations of RVE sizes and number of realizations require similar computation time. However, timings of FEA solutions are very dependent on the system and software used, so this behavior is only valid for this work. With the similar computation cost, the larger RVE size was chosen, since it exhibits more transversely isotropic behavior.

Predicted graphite fiber properties

Tables 2 –4 show the predicted fiber properties for the graphite fiber material systems using a hexagonal fiber arrangement, random fiber arrangement, and the Mori-Tanaka averaging scheme. As previously discussed, reported lamina properties generally present a dilemma in choosing whether to use the reported volume fraction or the volume fraction calculated from lamina properties and the manufacturer’s longitudinal fiber properties. As a result, the properties were calculated for both volume fractions when available. Some sources did not report a volume fraction, requiring its calculation using rule of mixtures. It was found that most properties are similar between the methods with the exception of the longitudinal shear modulus. The longitudinal shear modulus showed a higher sensitivity to the volume fraction for the Mori-Tanaka averaging scheme and hexagonal RVE than the random RVE, as seen in Tables 3 and 4.

Table 2.

Predicted fiber properties for IM7 and AS4 graphite fibers.

IM7 (5250-4)	IM7 (8552)			AS4
Calculated V_f	Calculated V_f			Calculated V_f	Reported V_f
Random	Mori-Tanaka	Hex	Random	Random	Random
E₁ (GPa)	276	276	276	276	231	207
E₂ (GPa)	26.0	26.5	24.9	22.4	29.9	22.6
G₁₂ (GPa)	20.7	19.6	19.3	12.0	32.2	20.6
G₂₃ (GPa)	7.55	9.07	8.33	7.53	12.13	8.96
ν₁₂	0.292	0.272	0.303	0.274	0.244	0.251

Table 3.

Predicted fiber properties for T300 graphite fibers.

T300
Calculated V_f			Reported V_f
Mori-Tanaka	Hex	Random	Mori-Tanaka	Hex	Random
E₁ (GPa)	230	230	230	227	227	227
E₂ (GPa)	33.5	29.4	24.3	31.8	28.2	23.4
G₁₂ (GPa)	99.9	91.5	16.0	71.4	66.7	15.2
G₂₃ (GPa)	15.14	12.12	9.90	14.13	11.50	9.60
ν₁₂	0.240	0.240	0.244	0.241	0.241	0.246

Table 4.

Predicted fiber properties for T650-35 graphite fibers.

T650-35
Calculated V_f			Reported V_f
Mori-Tanaka	Hex	Random	Mori-Tanaka	Hex	Random
E₁ (GPa)	243	243	243	225	225	225
E₂ (GPa)	18.0	17.2	16.0	15.4	14.9	14.2
G₁₂ (GPa)	49.7	48.8	13.7	21.1	20.9	11.0
G₂₃ (GPa)	7.04	6.59	6.07	5.84	5.55	5.28
ν₁₂	0.444	0.443	0.442	0.439	0.438	0.437

During this study, some material systems resulted in negative or infinite longitudinal shear moduli for the hexagonal unit cell and Mori-Tanaka averaging scheme. For these models, even rigid fibers will not yield a sufficient lamina shear modulus. This is due to the reduced interaction between fibers for the hexagonal array and Mori-Tanaka averaging scheme as compared to the random fiber model.

The IM7 properties obtained using the two different material systems showed significant differences, especially with the longitudinal shear modulus. The difference in resulting properties could be due to error in the material system properties reported by the literature or due to differences in the actual fiber properties between the composite systems.

Table 5 shows the values reported by literature for the fiber properties. These numbers are based on the reported volume fraction, if one was reported. It appears that the method used to predict the fiber properties for IM7 using the IM7/5250-3 material system had the same issue with the longitudinal shear modulus that FEA with a hexagonal RVE and the Mori-Tanaka averaging scheme encountered. With the exception of the longitudinal shear modulus, the values predicted using FEA with a random RVE were less than 6% different from the values in Ref. [11]. The properties for the AS4 and T300 fibers in literature were said to be assumed without an explanation of the method used to obtain them (fiber property determination was not the aim of the paper that reports these properties). The AS4 and T300 properties predicted using FEA with a random RVE differed from the reported values by 8–45%. The T650-35 properties from Ref. [10] were determined using the Mori-Tanaka averaging scheme, which matches the properties calculated in this work using the same method. The properties for T650-35 predicted using FEA with a random RVE differed from the reported values by less than 10% for all properties except the longitudinal shear modulus, which was 62% smaller for the random RVE. These comparisons to between models and results published in literature indicate that the predictions of the fiber’s Young’s moduli, Poisson’s ratios, and transverse shear modulus are somewhat insensitive to the assumed microstructure, but the predicted longitudinal shear modulus exhibits large differences between microstructures. This highlights the importance of consistently using the same microstructure to predict fiber properties and model the microstructural behavior of composite systems.

Table 5.

Fiber properties reported by the literature.

IM7 (5250-4) Ref. [11]	AS4 Ref. [24]	T300 Ref. [24]	T650-35 Ref. [10]
E₁ (GPa)	276	225	230	225
E₂ (GPa)	27.6	15	15	15.4
G₁₂ (GPa)	138	15	15	21.1
G₂₃ (GPa)	7.67	7	7	5.8
ν₁₂	0.30	0.2	0.2	0.44

Due to the large differences in the predicted longitudinal shear modulus for the microstructural models, the effect of spacing between fibers was investigated to understand the cause of the difference. For the material system of T650-35/PMR-15 with a volume fraction of 0.556, the minimum spacing fraction can theoretically vary from 0, which allows fibers to touch, to 0.554, which corresponds to hexagonal packing. For RVEs with 100 fibers, the algorithm used to generate the random RVEs cannot achieve volume fractions higher than about 83%. Consequently, the random RVE generation can only achieve a minimum spacing fraction between 0 and about 0.4 (See equation (2) with $V_{f} = 0.556 and V_{f}^{'} = 0.83)$ .

Figure 10 shows the variation of the predicted fiber properties as a function of the minimum spacing fraction. The longitudinal Young’s modulus, E₁, and longitudinal Poisson’s ratio, $v_{12}$ , were omitted since they varied by less than 0.3%. For all minimum spacing fractions considered, the transverse Young’s modulus, E₂, and transverse shear modulus, $G_{23}$ , varied by less than 5%. The longitudinal shear modulus, $G_{12}$ , was significantly affected by the spacing, increasing by 90%.

As the minimum spacing fraction is increased, the random microstructure becomes more uniform and forms fewer networks of fibers. As a result, the resulting lamina properties decrease with the spacing, causing the fiber properties needed to match the experimental lamina properties to increase with the spacing. This relationship helps to explain the significant difference between the properties predicted using a hexagonal microstructure and the random microstructure.

Predicted isotropic fiber properties

To validate the developed method, an epoxy matrix, glass fiber material system was analyzed. Literature reports the properties of the glass fibers to be isotropic with a Young's modulus of 80 GPa and Poisson's ratio of 0.2.²⁴ Only enforcing transverse isotropy for the fibers using the procedure previously described with a random microstructure does not result in isotropic predicted properties, as seen in Table 6. Similar to the problem found for some graphite fibers, the properties predicted using a hexagonal RVE and the Mori-Tanaka averaging scheme resulted in negative or infinite transverse Young’s and shear moduli. For these models, even rigid fibers will not yield sufficient lamina transverse moduli.

Table 6.

E-Glass 21×K43 predicted fiber properties using random RVEs.

E₁ (GPa)	80.0
E₂ (GPa)	63.0
G₁₂ (GPa)	15.4
G₂₃ (GPa)	25.8
ν₁₂	0.258

Enforcing isotropy of the fibers modifies the algorithm to have more outputs (five lamina properties) than inputs (two fiber properties). As a result, solving the correction term for the next iterative guess for fiber properties requires a least squares solution, making the method a modified Gauss-Newton algorithm. It is assumed that in general there is no guarantee that there exists a solution within tolerance when isotropy is enforced. The iterations are stopped when the residual fails to decrease. Effectively, this method finds the local minimum of the function space for the norm of the residual vector.

Again, several homogenization methods were used when forcing the glass fibers to be isotropic, including the Mori-Tanaka averaging scheme, FEA with the hex microstructure, and FEA with the random RVE microstructure. Table 7 shows the resulting properties. The Mori-Tanaka averaging scheme predicted the closest values to the literature-reported values of the Young’s modulus and Poisson’s ratio, and the hexagonal RVE was the closest to the reported value for the shear modulus. However, no two methods varied by no more than 7%. The predicted Poisson’s ratio for all three methods was quite a bit larger than the value reported by the literature. Compared to the properties predicted for transversely isotropic fiber material systems, the properties predicted for the isotropic fiber material system showed much less sensitivity to the model used.

Table 7.

E Glass 21×K43 fiber predicted properties forcing isotropy.

Mori-Tanaka	Hex	Random	Literature Ref. [24]
E (GPa)	79.7	81.5	79.0	80.0
G (GPa)	32.1	32.6	31.4	33.3
ν	0.242	0.251	0.260	0.2000

Conclusions

It was shown that the extent to which a set of random microstructures satisfies transverse isotropy depends on both the RVE size and the number of realizations. It was also shown that the average values of predicted lamina properties do not vary for the RVE sizes studied as long as a sufficient number of realizations are used, but the standard deviation of the properties (a measure of the distribution of properties within the set of realizations) decreases for larger RVE sizes. Based on these findings, an appropriate RVE size and number of realizations was identified to fully characterize both the average lamina properties and the standard deviation of lamina properties.

Comparing the predicted fiber properties using three methods, it was noted that the microstructure has a significant impact on the predicted longitudinal shear modulus of the fiber. By examining the relationship between the spacing of fibers and the predicted properties, the longitudinal shear modulus of the fiber was shown to increase approximately linearly with the minimum spacing fraction. This is because an increase in the minimum spacing fraction causes the random microstructure to be more uniform and form fewer networks, causing the lamina to be more compliant for higher spacing fractions and requiring a stiffer fiber to match experimental results. The properties predicted for the IM7 and T650-35 fibers showed close similarity to those presented in the literature with the exception of the longitudinal shear modulus, which was much lower for the random microstructure. The properties for AS4 and T300 found in the literature differed significantly from the properties predicted in this work; however, it was noted that the properties in the literature had inconsistencies.

For the glass material system, if isotropy of the fibers was not enforced, the resulting fiber properties were not isotropic. As a result, isotropy was enforced, and the fiber properties were determined using three different microstructural models. The predicted isotropic moduli matched the values reported in the literature quite well, although the Poisson’s ratio was higher for the methods used in this work than the value reported. If the fibers are known to be close to isotropic and isotropic properties are desired, the algorithm must be modified to enforce the isotropy of the fibers.

Footnotes

Acknowledgement

The authors also acknowledge the Texas A&M Supercomputing Facility () for providing computing resources useful in conducting the research reported in this paper.

Funding

This work was supported by the Multidisciplinary University Research Initiative grant FA9550-09-1-0686 from the Air Force Office of Scientific Research to Texas A&M University with Dr. David Stargel as the program manager.

Conflict of Interest

None declared.

References

Mollenhauer DH, Iarve EV, Hoos K, et al. Simulation of mode I fracture at the micro-level in polymer matrix composite laminate plies. In: 53nd AIAA/ASME/ASCE/AHS/ASC structures, Structural Dynamics, and Materials Conference, 2012.

Trias

Costa

Fiedler

. A two-scale method for matrix cracking probability in fibre-reinforced composites based on a statistical representative volume element. Compos Sci Technol 2006; 66: 1766–1777.

King

Blackketter

Walrath

. Micromechanics prediction of the shear strength of carbon fiber/epoxy matrix composites: the influence of the matrix and interface strengths. J Compos Mater 1992; 26: 558–573.

Adams

Crane

. Finite element micromechanical analysis of a unidirectional composite including longitudinal shear loading. Comput Struct 1984; 18: 1153–1165.

Goyal

Whitcomb

. Effect of fiber properties on plastic behavior of 2×2 biaxial braided composites. Compos Sci Technol 2008; 68: 969–977.

Fiedler

Gagela

Hobbiebrunken

. Modelling of the transverse strength of fibre reinforced epoxy composite at low and high temperature. Compos Interfaces 2005; 12: 379–394.

Fiedler

Hojo

Ochiai

. The influence of thermal residual stresses on the transverse strength of CFRP using FEM. Compos Part A Appl Sci Manufact 2002; 33: 1323–1326.

Jin

. Effect of fiber arrangement on residual thermal stress distributions in a unidirectional composite. J Compos Mater 2006; 41: 591–611.

Yang

Yan

. Effects of inter-fiber spacing and thermal residual stress on transverse failure of fiber-reinforced polymer–matrix composites. Computat Mater Sci 2013; 68: 255–262.

10.

Rupnowski

Gentz

Sutter

. An evaluation of the elastic properties and thermal expansion coefficients of medium and high modulus graphite fibers. Compos Part A Appl Sci Manufact 2005; 36: 327–338.

11.

Pagano

Schoeppner

Kim

. Steady-state cracking and edge effects in thermo-mechanical transverse cracking of cross-ply laminates. Compos Sci Technol 1998; 58: 1811–1825.

12.

Pagano

Tandon

. Thermo-elastic model for multidirectional coated-fiber composites: traction formulation. Compos Sci Technol 1990; 38: 247–269.

13.

González

LLorca

. Mechanical behavior of unidirectional fiber-reinforced polymers under transverse compression: microscopic mechanisms and modeling. Compos Sci Technol 2007; 67: 2795–2806.

14.

Brockenbrough

Suresh

Wienecke

. Deformation of metal-matrix composites with continuous fibers: geometrical effects of fiber distribution and shape. Acta Metall Mater 1991; 39: 735–752.

15.

Jin

K-K

. Interfacial strain distribution of a unidirectional composite with randomly distributed fibers under transverse loading. J Compos Mater 2005; 40: 759–778.

16.

Garnich MR, Fertig RS, Anderson EM, et al. Micromechanics of fatigue damage in unidirectional polymer composites. In: 53nd AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics, and materials conference, April 2012, pp.1–9.

17.

Anderson

. An automated finite element program for micromechanics modeling of random-wavy fiber composites, M.S. thesis, Dept. Mech. Eng., Univ. of Wyoming, Laramie, WY, 2010, 2010.

18.

Pipes

Pagano

. Interlaminar stresses in composite laminates under uniform axial extension. J Compos Mater 1970; 4: 538–548.

19.

Whitcomb

Chapman

Tang

. Derivation of boundary conditions for micromechanics analyses of plain and satin weave composites. J Compos Mater 2000; 34: 724–747.

20.

McLendon WR and Whitcomb JD. Micro-scale analysis for the prediction of strength under biaxial thermomechanical load. In: Proceedings of the American Society for Composites, 2012.

21.

Shewchuk

. Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator. Appl Comput Geometry 1996; 1148: 203–222.

22.

Cherkaoui

. Fundamentals of micromechanics of solids, 1st ed. Hoboken, New Jersey: Wiley, 2006, pp. 400–400.

23.

Huang

C-Y

Trask

Bond

. Characterization and analysis of carbon fibre-reinforced polymer composite laminates with embedded circular vasculature. J Royal Soc Interface 2010; 7: 1229–1241.

24.

Soden

Hinton

Kaddour

. Lamina properties, lay-up configurations and loading conditions for a range of fibre-reinforced composite laminates. Compos Sci Technol 1998; 58: 1011–1022.