The effect of model dimensionality on compression strength of fiber reinforced composites

Abstract

Strength of fiber reinforced polymers (FRPs) under compression loads is typically limited by a shear localization failure mode called microbuckling which is highly sensitive to fiber misalignment. In addition to the magnitude of fiber misalignment, the dimensionality of fiber misalignment also plays a prominent role in the prediction of compression strength. Therefore, a comparison between 1D, 2D, and 3D fiber misalignment is carried out in a finite element setting with a homogenized representation of fiber and matrix materials. In real FRP structures, fiber misalignment is spread in a correlated random manner throughout the material volume resulting in a distribution of compression strength. Spectral representation method is used for developing the volumetric representation of fiber misalignment in numerical models, thus preserving the spatial correlations of fiber misalignment. As an input to the spectral representation method, two different functional forms of spectral density of the fiber misalignment are considered. The results of model series based on functional forms of spectral density are also compared against a reference model series based on experimental measurements of the fiber misalignment. Finally, a simple relation is proposed for prediction of compression strength at different percentiles of distribution of failure strengths with regard to scaling of mean square spectral density.

Keywords

Fiber reinforced polymers fiber misalignment strength probabilistic methods microbuckling

Introduction

Continuous fiber reinforced polymers (FRPs) are increasingly becoming a material of choice for advanced structural components because of their exceptionally high strength and stiffness to weight ratios. The fields of applications of FRPs include highly loaded components such as fins, rudders and fuselage panels of modern aircraft, blades of wind turbines, and structural components of modern rockets among others.^1,2 The compression strength of unidirectional FRPs can be as low as 60% of their tensile strength and shows uncertainty.^3,4 Therefore, prediction of compression strength becomes highly relevant for reliable design practices.

Compression strength of FRPs is dictated by a phenomenon referred to as microbuckling (MB), also known as kinking, in the literature. Microbuckling is divided into an extensional/symmetrical/tensile mode and a shear mode.^5–7 However, shear mode MB is of interest as it occurs in most of the industrial scale, high volume fraction FRPs. Other failure modes such as longitudinal splitting may compete with MB in low volume fraction FRPs, however, MB is the dominant failure mode in high volume fraction FRPs.^5,7–9 The MB failure is controlled by the matrix shear nonlinearity and the manufacturing induced fiber misalignment. The MB failure is initiated when geometric softening at a location of a critical fiber misalignment outpaces the material hardening inside the incipient shear band.^10–13 The MB failure mode is characterized by rotation of fibers in a localized band of shear deformation which leads to a sudden catastrophic reduction in the far field compression load. Further application of load matures the band to its final state before ultimate fracture of fibers at the ends of the band.

Different analytical and numerical models presented in the literature indicate the complexity and the difficulty associated with the accurate prediction of compression failure. Initial approaches to modeling compression failure mainly considered a single valued fiber misalignment angle ψ and presented analytical solutions.^6,14–16 The analytical model by Budiansky, which can be found in Ref.¹⁰ as equation (3), is widely used for prediction of compression strength. It considers compression strength σ^c to be dependent on the initial fiber misalignment ψ and the matrix shear nonlinearity τ₁₂(γ) as a function of the shear angle γ, see equation (1).

σ^{c} = \max_{γ} (\frac{τ_{12} (γ)}{ψ + γ})

(1)

Although analytical models are valuable to highlight the importance of matrix nonlinearity and the fiber misalignment, such as one presented above which takes a single value of the fiber misalignment as an input, they are limited in prediction of the compression strength of real structures.

Improving on the assumptions of the analytical modeling, different numerical models were introduced in literature over the decades. Typical numerical approaches consider a sinusoidal undulation of the fiber misalignment. The common assumption of an infinite band, in which imperfection is considered uniform in transversal directions and only varies in longitudinal direction, is commonly applied in 1D numerical models.^17,18 2D models considering variation of the fiber misalignment in longitudinal as well as in transversal direction with idealized forms of the misalignment can also be found in the literature.^19–22 In the same spirit, 3D models^23–29 have been used to explore different features associated with MB failure such as mechanisms leading to failure, band angle, and band width.

Aforementioned models have mostly considered idealized forms of the fiber misalignment. The fiber misalignment in real structures is correlated random in nature.^30–34 The random nature of the fiber misalignment observed in measurements requires modeling approaches where such variations of the fiber misalignment can be accounted for in the numerical models. Some micromechanical approaches incorporate misalignment variations, such as the investigation by Varandas et al.³⁵ with the help of a stochastic process and an optimization procedure, however, the difficult model generation and high computational costs deem such an approach relatively unsuitable for simulating a large number of model realizations. To overcome these problems, homogenized approaches offer a better alternative. In these approaches, fiber and matrix are modeled using a smeared material representation. In the homogenized approaches, a single material model accounting for the anisotropy represents the fibers and the matrix in a smeared sense. Material nonlinearities, such as matrix shear behavior which is essential for modeling microbuckling, are included in such models along with geometrical nonlinearities through different approaches such as plasticity or damage modeling. The first such approach in 1D was introduced by Slaughter et al.³⁶ Controlling the correlation of the fiber misalignment distribution using the spectral representation method, distributions of the compression strength were predicted. The spectral density of a given distribution of the fiber misalignment such as an experimentally measured distribution can be used to generate virtual distributions with the spectral representation method.^37–39 The spectral representation method offers a direct bridge between the measured spatial distributions of fiber misalignment and their virtually generated counterparts. Using discrete spectral densities as an input to the spectral representation method, the topologies of the fiber misalignment in the form of spatial distributions of misalignment angles or amplitudes are generated. These spatial distributions are in turn used to define local material orientations in the homogenized approaches. It is ensured this way, that not only the statistical information of the fiber misalignment but also the spatial distribution is represented in a model. Spectral representation method models all the wavelengths in the range of the discretization length (minimum wavelength) and the model size (maximum wavelength). Extending Slaughter’s approach to 2D, Liu et al. investigated the effects of different modeling parameters on the prediction of compression strength.⁴⁰ Sutcliffe⁴¹ presented an alternate simplified approach to Liu et al. where correlation of the fiber misalignment distribution was controlled in an apparent artificial manner. Another approach was presented by Allix et al.⁴² in 3D albeit with a random spread of the misalignment with no correlation.

Realizing the potential benefits of easier modeling and lower computational costs, a homogenized modeling approach using the finite element (FE) method is used in the current contribution. Previous homogenized investigations in the literature with random representation of the fiber misalignment are limited to 1D,³⁶ 2D,^40,41 and 3D but with a 2D fiber misalignment distribution.⁴² The fiber misalignment in FRP laminates is in reality three dimensional. To highlight the effect of dimensionality of the fiber misalignment on the prediction of the compression strength, an investigation using simple sinusoidal form of the misalignment is performed using a homogenized modeling approach. Realizing the importance of 3D modeling for realistic representation of FRP laminates, 3D distributions of the fiber misalignment calibrated from 3D spectral densities are used in the current investigation. For this purpose, the misalignment distributions are generated using an extended form of the spectral representation method. Square and exponential forms of the spectral densities are considered and their effects on the distributions of the predicted strength are explored. Furthermore, the effects of scaling the mean square spectral density, which is the area under the spectral density curve and is a material quantity, are investigated. A simple relation for prediction of compression strength of models with scaled misalignment is also provided.

Model description

For the prediction of the compression strength, accurate modeling of the fiber misalignment as well as an accurate description of the material model are necessary. The modeling methodology follows the approach presented in an earlier contribution of the authors, see Ref.⁴³ Therefore, only a brief description of the material model and the geometric model will be summarized in this section. Virtual distributions of the fiber misalignment are generated using the spectral representation method. The use of different functional forms of the spectral density for generation of misalignment distributions will be presented in detail.

Material model

As concluded from the state of the art review of modeling for compression failure of FRPs, homogenized approaches offer ease of modeling over larger sizes with lower computational effort. The anisotropy is usually incorporated by the orientation of the preferred material axis. In a numerical model, the fiber misalignment is included via orientation of the vector representing the preferred material axis at the integration point of the finite elements. For the high volume fraction FRPs having a tight hexagonal packing of the fibers in the matrix materials, the material behavior is usually assumed to be transversely isotropic.⁴ Hence, in current investigation the homogenized fiber-matrix material having a 55% fiber volume fraction was modeled using a transversely isotropic material model. The anisotropy is represented using a so-called structural tensor which is defined based on the vector of the initial fiber direction. The material nonlinearities were accounted for using plasticity. To predict realistic plastic deformations, a nonassociative flow rule was assumed. The material model has been implemented in the commercial software Abaqus, and is an extension of the small strain transversely isotropic material model for unidirectional FRPs.⁴⁴ The process for determining the parameters of the material model is given in Ref.⁴⁵ The geometrical nonlinearities were accounted for using a co-rotational formulation.⁴⁵

Wisnom⁴⁶ highlighted a limitation of using a homogenized material representation in the commercial software Abaqus as it ignores the effect of shear strains on the material axes rotation. By default, Abaqus ignores the effect of shear strains on the material axes rotation and only considers the material axes rotation to be identical with the rigid body rotations. Since the effect of shear strains on the material axes rotation during loading is crucial in case of compression, the default implementation in Abaqus gives inaccurate results. To overcome this problem, the material model from Ref.⁴⁵ was modified to account for the rotation of the material axes due to shear strains. The shear strains from the each iteration of the solution were considered as Euler angle rotations to update the current material axis orientation. A verification against Budiasnky’s analytical solution using a single element model over a range of initial misalignment angles with and without the modification was presented in Ref.⁴³ It was shown that without this correction, the numerical model overpredicted compression strength.

A question that arises in this context is, how many angles would be needed to accurately model the three-dimensional misalignment of a fiber in volumetric space from its mean direction in unidirectional FRPs. Following the Euler angles notation in combination with the realization that fibers’ rotation around their own axis does not affect MB failure, it is sufficient to have two independent angles for modeling the fiber misalignment in FRPs. These two angles correspond to the in-plane misalignment angle θ and the out-of-plane misalignment angle ϕ. This line of argumentation is coherent with the transversely isotropic behavior observed in unidirectional FRP materials. Therefore, these two angles θ and ϕ are sufficient for a complete volumetric representation of the total misalignment ψ from the nominal fiber direction i.e. 1-direction. A schematic illustration of the misaligned preferred direction vector a′ (in black) with the initial in-plane misalignment θ and the out-of-plane misalignment ϕ contributing to the total misalignment ψ is depicted in Figure 1. Here the nominal, i.e. perfectly aligned fiber direction, corresponds to the 1-direction.

Figure 1.

Schematic of a fiber (in black color) misaligned from the nominal fiber direction 1 to direction 1′.

Material model used in all the investigations of this contribution was calibrated from the experimental data gathered in Ref.⁴⁷ There, a prepreg UD300 carbon fibers with M79 epoxy resin Hexcel system was considered. The material characterization was carried out following the three standard tests in Ref.⁴⁷ These tests were performed following the established standards i.e. axial tension test following DIN EN 2561 standard for the material characterization in fiber direction, EN ISO 14129, ASTM D3518, DIN 65466 standard for the in-plane shear characterization, and transverse tesntion test following DIN EN 2597 standard for the measurement of transverse stiffness modulus. The material properties obtained from Ref.⁴⁷ are summarized in Table 1. There, E₁₁, E₂₂ = E₃₃, G₁₂ = G₁₃, and ν₁₂ = ν₁₃ represent axial modulus, transverse modulus, in-plane shear modulus, and major Poisson’s ratio, respectively. The value of major transverse Poisson’s ratio ν₂₃=0.5 was taken from the material datasheet provided by the suppliers.

Table 1.

Material properties of carbon fiber UD300-M79 epoxy matrix with a fiber volume fraction v_f of 55%.

E11[GPa]	E22 = E33[GPa]	G12 = G13[GPa]	ν12 = ν13[−]	ν23[−]
118.43	7.03	3.88	0.299	0.5

Misalignment topologies

Reference topologies driven from measurement of the misalignment

A detailed description of the procedure for acquiring the misalignment data may be found in Ref.⁴⁷, and only a brief synopsis is presented here. Ideally an ensemble averaged value of the spectral density has to be calculated by measuring misalignment experimentally in a statistical manner for it to be a material characteristic quantity. However, only two specimens were scanned for characterization of the spectral densities because of the steep costs per CT-scan at the required resolution.⁴⁷ The specimens were scanned by an X-ray microscope with a resolution of 1.902 µm. A voxel resolution of 1.902 µm was selected to ensure a fine enough resolution at which the algorithms for the measurement of fiber angles can differentiate between matrix material and fibers. The selected resolution limited the scanned cylinder diameter to approximately 2 mm. The volume of the composite in which the algorithm was able to detect fiber and matrix materials separately for both the in-plane θ and the out-of-plane ϕ angles was approximately 6.657 × 1.1412 × 0.951 mm. A schematic of a cuboid volume selected from a cylindrical scanned volume is shown in Figure 2. The viewing planes in 1–2 and 1–3 planes, i.e., the in-plane and the out-of-plane, are marked in green and red colors, respectively. Using the image processing technique “High Resolution Misalignment Analysis” (HRMA) proposed by Wilhelmsson and Asp,⁴⁸ the voxel data obtained from the CT-scans was used to determine both the θ and the ϕ angles. Image processing rendered 70 × 12 × 10 matrices for θ_k,l,m and ϕ_k,l,m over a volume of 6.657 × 1.1412 × 0.951 mm from each specimen. Subsequently, a discrete Fourier transformation rendered matrices of 70 × 12 × 10 for the corresponding spectral densities $S_{k, l, m}^{θ}$ and $S_{k, l, m}^{ϕ}$ . The spectral densities calculated from angle measurements of both samples were averaged. The plots of experimentally calculated spectral densities of each specimen showed rather spiky data but averages thereof were comparitevly smooth, cf. Figure 2 in Ref.⁴³ It indicates that further CT-scans of specimens would lead to even smoother data in averaged spectral densities. However, it is to be noted that the spiky data of the spectral densities does not affect the accuracy of the evaluation. Use of the average characterized spectral densities in generation of the fiber misalignment topologies, which resulted in smooth contours of the fiber misalignment angles of the reference model shown in Figure 3, supports this argument. The mean square values Ξ^θ and Ξ^ϕ of spectral densities are a measure of the spread of distributions of the misalignment angles θ_k,l,m and ϕ_k,l,m in frequency domain, and the average values from measurements of both samples are given in Table 2. Average values of standard deviations s^θ and s^ϕ of misalignment angles are also provided whereas the mean of each distribution is zero, per definition. The relation between s^θ and Ξ^θ is given in equation (2), and the same holds for s^ϕ and Ξ^ϕ.

Ξ^{θ} = \frac{L_{1} L_{2} L_{3}}{{(2 π)}^{3}} {(s^{θ})}^{2}

(2)

where L₁, L₂, and L₃ are lengths in directions 1, 2, and 3, respectively.

Figure 2.

Sketch of a cuboid volume selected from a scanned cylindrical volume. The cylindrical scanned volume was obtained over a cuboid shape [0]₄ unidirectional FRP laminate by performing a CT scan.

Figure 3.

Measured and virtually generated 3D contour plots of the in-plane θ_k,l,m and the out-of-plane ϕ_k,l,m misalignment angles from the reference model series. The reference model series was based on measured misalignment distributions using CT scans.

Table 2.

Average mean square spectral densities and standard deviations of fiber misalignment angles obtained from CT measurements.⁴³

	In-plane θ_k,l,m	Out-of-plane ϕ_k,l,m
Ξ^θ [ $°^{3}$ ] and Ξ^ϕ[ $°^{3}$ ]	0.0139	0.0071
s^θ[°] and s^ϕ[°]	0.6920	0.4920

In an earlier numerical contribution by the authors, which can be found in Ref.⁴³, the spectral representation method was extended to 3D. Using the extended spectral representation method with input of discrete matrices of spectral densities $S_{k, l, m}^{θ}$ and $S_{k, l, m}^{ϕ}$ calculated from the measurements of the fiber misalignment θ_k,l,m and ϕ_k,l,m, 250 realizations of virtual distributions of θ_k,l,m and ϕ_k,l,m were generated. These sets of virtual distributions of θ_k,l,m and ϕ_k,l,m were given as an input to the homogenized FE models, with dimensions of 6.657 × 1.1412 × 0.951 mm, for developing a reference model series. A particular realization each of θ_k,l,m and ϕ_k,l,m for the reference model series is shown in Figure 3. It is to be noted that the distribution of ϕ_k,l,m shown in Figure 3(d) is more uniformly spread in contrast to the distribution of θ_k,l,m shown in Figure 3(b). These differences were also measured in the experimentally scanned specimens in Ref.,⁴⁷ and can be transferred to the virtually generated counterparts with the help of spectral representation method.

Idealized topologies

In order to analyze the effect of the fiber misalignment dimensionality on the strength prediction, a comparative analysis using simple 2D and 3D models containing only single sinusoidal undulation are considered here. Four models, two in 2D with 4 × 2 mm in dimension and two in 3D with 4 × 2 × 2 mm in dimension, were generated in the preliminary series for the purpose of comparing the effect of the underlying undulation dimensionality on the compression strength prediction. The first 2D model contains a half period of a sine wave undulation, in which the undulation is variable in only 1-direction. This model corresponds to the infinite band assumption and may be regarded as a 1D-misalignment topology, i.e. of dimensionality 1. The resulting contour plot of the undulation is shown in Figure 4(a). The second 2D model contains an undulation with two dimensional sine waves, i.e. the dimensionality of the misalignment is now increased to 2. The resulting contour plot of the undulation is shown in Figure 4(b).

Figure 4.

Contour plots of 1D and 2D sine waves for 2D modeling, and 2D and 3D sine waves for 3D modeling of the in-plane misalignment angle θ. The maximum misalignment in each case is 2.5° and the minimum misalignment in each case is −2.5°.

The first 3D model also contains a 2D undulation i.e. it contains a two dimensional sine wave, see Figure 4(c). Hence, the dimensionality of the misalignment is, again, 2. The second 3D models contains a fully three dimensional sinusoidal undulation i.e. a 3D misalignment. The resulting contour plot of the undulation is shown in Figure 4(d). The misalignment of both 3D models with a single sinusoidal undulation are given as the in-plane misalignment and zero out-of-plane misalignment in the FE models for a direct comparison with the 2D models. The undulations in the models of the preliminary series correspond to fundamental frequencies in each direction, i.e. the wavelength is equal to the model length in that particular direction. The maximum and minimum values of the angle in the undulation in all the models of the preliminary series are 2.5° and −2.5°, respectively. The corresponding statistics of the fiber misalignment angles and the spectral densities of each of the four models are given in Table 3.

Table 3.

Statistics of fiber misalignment angles and corresponding spectral parameters. Even though s^θ values for 2D model with 2D sine wave and 3D model with 2D sine wave are identical, the corresponding Ξ^θ values are different. This is because of different volumes in 2D and 3D.

Model	Ξ^θ [ $°^{2}$ ] or [ $°^{3}$ ]	s^θ [°]	max \|θ\| [°]
2D model with 1D sine wave	0.6355	1.7720	2.5
2D model with 2D sine wave	0.3178	1.2532	2.5
3D model with 2D sine wave	0.1012	1.2532	2.5
3D model with 3D sine wave	0.0674	1.0221	2.5

Generated equivalent topologies

Previous authors considered the idea of a functional form of the spectral density for generating virtual distributions of fiber misalignment. Herein, the previous assumptions applied to 1D by Slaughter et al.³⁶ and 2D by Liu et al.⁴⁰ are extended to 3D. Moreover, further aspects of such an assumption for the spectral density are explored.

Spectral density functions

Two forms of spectral densities were considered for generating both θ_k,l,m and ϕ_k,l,m distributions namely, the square function S_square given in equation (3) and the exponential function S_exp given in equation (4).

S_{s q u a r e} = S (ω_{1}, ω_{2}, ω_{3}) = {\begin{cases} S_{o} if | ω_{1} | \leq ω_{c 1} and | ω_{2} | \leq ω_{c 2} and | ω_{3} | \leq ω_{c 3} \\ 0 if | ω_{1} | > ω_{c 1} or | ω_{2} | > ω_{c 2} or | ω_{3} | > ω_{c 3} \end{cases}

(3)

S_{\exp} = S (ω_{1}, ω_{2}, ω_{3}) = S_{o} \exp - ({(\frac{ω_{1}}{ω_{c 1}})}^{2} + {(\frac{ω_{2}}{ω_{c 2}})}^{2} + {(\frac{ω_{3}}{ω_{c 3}})}^{2})

(4)

where ω_c1, ω_c2, and ω_c2 are the cutoff frequencies in direction 1, 2, and 3, respectively. The mean square values Ξ^θ and Ξ^ϕ of the spectral densities were used for comparing S_square and S_exp based models with each other, and the subsequent comparison with the reference model driven from the CT measurements of the fiber misalignment. A section of the spectral density functions in plane ω₁/ω₂ with ω₃=0 and anisotropic cutoff values of the frequencies, i.e. where ω_c1 was not equal to ω_c2, are illustrated in Figure 5. The difference between the square and the exponential forms of the spectral density can be observed in Figure 5(a) and (b). The square form gives an equal weightage to all the underlying frequencies-ω (wavelengths - λ) in the virtually generated distributions of the fiber misalignment derived from this functional form. On the other hand, the exponential form of spectral density gives a higher weightage to smaller frequencies (larger wavelengths) in virtually generated distributions of the misalignment which is closer to the experimentally observed spectral densities, see.⁴³

Figure 5.

Plots of assumed functional forms of spectral densities. The figures show spectral densities in plane ω₁/ω₂ with anisotropic cutoff values of the frequencies i.e. ω_c1 not equal to ω_c2. The cutoff frequencies ω_c1 and ω_c2 as well as the initial value of spectral density S_o are also marked in the figures. In these figures ω₃ = 0, and is disregarded to allow 3D plots.

Resultant topologies

In order to understand the effects of the parameters as well as the shape of the assumed spectral density on the prediction of the compression strength, different virtual distribution of the in-plane misalignment θ_k,l,m and the out-of-plane misalignment ϕ_k,l,m were generated. The dimensions of the models were kept the same as the reference model i.e. 6.657 × 1.1412 × 0.951 mm with equal discretization of 0.0951 mm in each direction.

The virtual distributions of the misalignment for θ_k,l,m and ϕ_k,l,m were generated by scaling the functional forms of spectral densities in such a way that the mean square spectral densities Ξ^θ and Ξ^ϕ for each misalignment correspond to the experimentally measured value given in Table 2. The functional forms of spectral densities were scaled by the parameter S_o, called the initial spectral density. The first model series considered no cutoff frequency which implies that the model contained all possible frequencies from fundamental frequency (largest wavelength) to the Nyquist frequency (smallest wavelength) in each direction. Figure 6(a) and (c) show 3D contour plots of θ_k,l,m for a realization each using square and exponential forms of the spectral density without a cutoff frequency, respectively. Similarly, contour plots ofϕ_k,l,m for a realization each are shown in Figure 7(a) and (c) using square and exponential forms of spectral density, respectively. It can be seen that without using a reasonable cutoff in the discrete sampling of the spectral density, the resulting virtual distributions of the misalignment show an unreasonable noise in the pattern which is unlike the experimentally measured distributions and subsequent reference model distributions shown in Figure 3. The noise in the pattern results in much smaller regions of misalignment compared to the measured and virtually generated 3D reference topologies shown in Figure 3. This is due to the effect of very small wavelengths contributing to the pattern of the virtually generated distributions.

Figure 6.

Virtually generated 3D contour plots of the in-plane misalignment angle θ_k,l,m using functional forms of the spectral density as an input. The spectral representation method was used for the generation of virtual distributions of misalignment angles. Although contour plots based on S_square and S_exp without a cutoff appear quite similar, the underlying wave behavior is different. These differences lead to differences in the prediction of compression strength distributions.

Figure 7.

Virtually generated 3D contour plots of the in-plane misalignment angle ϕ_k,l,m using functional forms of the spectral density as an input. The spectral representation method was used for the generation of virtual distributions of the misalignment angles.

To improve on the distributions generated using functional forms of spectral density without a cutoff, a second model series with certain cutoffs was considered. For this purpose, cutoff frequencies of ω_c1=6.28 (corresponding to ${λ_{1}}_{\min}$ = 1 mm) and ω_c2=ω_c3=12.56 (corresponding to ${λ_{2}}_{\min}$ = ${λ_{3}}_{\min}$ = 0.5 mm) were used for the in-plane misalignment θ_k,l,m. For the out-of-plane misalignment ϕ_k,l,m, cutoff frequencies of ω_c1 = ω_c2 = ω_c3 = 12.56 (corresponding to ${λ_{2}}_{\min}$ = ${λ_{3}}_{\min}$ = 0.5 mm) were used. These values of the cutoff frequencies are approximately equal to the values observed in the experimentally measured distributions of the fiber misalignment, see Figure 3. A realization each from the resulting virtually generated distributions of θ_k,l,m using square and exponential functional forms of spectral density are plotted in Figure 6(b) and Figure 6(d), respectively. Similarly, Figure 7(b) and (d) show a realization each of ϕ_k,l,m using square and exponential spectral density, respectively. It can be observed that by using appropriate cutoff frequencies when discretizing the functional forms of the spectral density, the generated distributions resembles the experimental distributions shown in Figure 3 quite closely. Additionally, the anisotropic behavior in θ_k,l,m of the resulting pattern in the virtually generated distributions can now be modeled. It shows that the cutoff frequencies are not necessarily equal in all directions as assumed by the previous authors in literature.

Model geometry

Different series of models were considered in the current investigation. A first preliminary series is used for investigating the effect of misalignment dimensionality on the compression strength. It comprises 2D plane stress models with 4 × 2 mm in dimension and 3D models in the form of cuboids of 4 × 2 × 2 mm in dimension, each containing simple sinusoidal undulation. These dimensions are chosen for sake of simplicity as the purpose of first series is comparison of 2D vs 3D modeling. The second series of models are intended for performing probabilistic analyses under axial compression load. The probabilistic series comprise reference models whose topology is based on experimental measurements of misalignment, see reference topologies section. It also contains models which are based on functional forms of spectral densities, see generated topologies section. All models of the probabilistic series are in 3D in the form of cuboids, representing a region of a unidirectional FRP [0]₄-laminate. Figure 8 illustrates the schematic of the 3D FE model geometry from the probabilistic series. The geometry of models from the preliminary series is similar. Nominal fiber direction was parallel to the 1-axis. The 2 and the 3-directions were the in-plane and the out-of-plane transverse directions, respectively. The FE model dimensions and the discretization were kept same as the volumetric dimensions and the discretization of the virtually generated topologies presented in misalignment topologies section. The nominal dimensions of the FE models for the probabilistic series were 6.657 × 1.1412 × 0.951 mm in 1-2-3 directions, respectively.

Figure 8.

Schematic of the FE model geometry used for the probabilistic series. (a) depicts a [0]₄ laminate with large in-plane dimensions compared to the thickness in 3-direction. The laminate was stacked in direction 3. The laminate consists of representative volumes which are periodic in directions 1 and 2. One such volume is marked in dark blue color. (b) illustrate the geometry of the representative volume used in numerical modeling. The contour of an in-plane distribution of the fiber misalignment angles is plotted on the back planes of the cuboid shaped volume. The direction 1 is the nominal fiber direction, and 2 and 3 are the in-plane and the out-of-plane transverse directions, respectively. The FE mesh is shown with the grid on the back planes of the cuboid with blue lines.

3D models were meshed using cubic three dimensional fully integrated 20 node C3D20 elements in the commercial software Abaqus 2016, whereas 2D models were meshed using 8 node CPS8 plane stress elements. The resulting mesh for the preliminary series of models comprised 40 × 20 square elements of 0.1 mm edge length and 40 × 20 × 20 cubic elements with 0.1 mm edge length for 2D and 3D models, respectively. The resulting mesh for the probabilistic series comprised 70 × 12× 10 cubic elements with 0.0951 mm edge length, each. The distributions of the fiber misalignment generated using the spectral representation method were mapped onto the integration points of each element of the corresponding mesh as the material orientations. This was achieved by defining a local coordinate system for each finite element. The local coordinate systems were defined in such a way that they were rotated in-plane and out-of-plane from the global coordinate system by the corresponding values of misalignment angles at those spatial locations. A script written in Python Scripting Interface of Abaqus was used for this purpose.

Periodic boundary conditions were applied in directions 1 and 2 for the preliminary series of 2D and 3D models, whereas plane stress condition was applied on 3D models of the preliminary series in direction 3 for a direct comparison with 2D models. This was achieved by applying multi point constraints (MPC) on surfaces in direction 3 of the 3D models of the preliminary series using the method employed by Ref.⁴⁹ As the model schematic in Figure 8 shows, the FE model the probabilistic series represents a section of a large plate extending in directions 1 and 2. The model thickness is same as the laminate, hence, there is no periodicity in direction 3. Therefore, the 3D FE models of the probabilistic series were representative of a cutout region from a unidirectional FRP laminate based panel with same thickness but larger in-plane dimensions, hence, periodic boundary conditions were used in directions 1 and 2 only. Both the surfaces perpendicular to direction 3 for the probabilistic series had traction free boundary conditions. To suppress rigid body motions in 3D models of both series corner node C5 was constrained in all directions and corner node, C8 and C6 were fixed in both transverse directions and in-plane transverse direction, respectively. Constraints and loads on 2D models were applied in a similar fashion. The axial compression load was applied in the form of displacements u1 at the corner node C1. Node C1 was coupled to the whole of the right surface/edge, therefore, the load applied on node C1 got transferred to whole of the right surface/edge. Since the surfaces/edges in direction 1 were periodic, the coupling of node C1 with the right surface/edge ensured the periodic application of load on the right and left surfaces/edges. Additionally, node C1 was constrained in transverse directions 2 and 3.

It is well established in the literature that MB is the failure mode under axial compression load for unidirectional FRPs. The global response of the model until the MB-failure is approximately linear. Locally at the location of the critical fiber misalignment, there is significant matrxi hardening before MB failure. When the geometric softening of the material due to compression load outpaces the materix hardening at a local location, the material fails in an unstable manner with a snap-back response. The failure is defined as the point after which the material suddenly sees a drastic reduction in the global load carrying capacity. Thus, the failure point of the model was defined as the peak load. Since the peak load is of interest, the analyses were terminated when reaction force RF₁ in the direction 1 at the control node C1 dropped by 5% from the peak load with the help of URDFIL user subroutine of Abaqus⁵⁰ in order to speed up. Moreover, the solution procedure was performed in two load steps for higher efficiency. The 1^st step up to approximately 40% of the expected mean failure strain was solved with larger sub-steps as the response was linear in this range. Both the implicit dynamic and the arc-length based Riks′ solver gave the same prediction up to peak load. The major advantage of the arc-length method is the prediction of the snap-back response, which is out of scope of the current work. The arc-length method tends to suffer from convergence problems whereas the implicit dynamic solver is more reliable and shows reliable convergence past the peak load. With the implicit dynamic method, there was the added advantage of easier post processing because the reaction force RF₁ and displacement u1 from multiple steps could be extracted simultaneously, which was not possible in the Riks′ arc-length method in Abaqus. Therefore, the implicit dynamic solver was deemed a better choice for probabilistic simulations. Hence, the geometrically non-linear dynamic implicit solver with quasi static load application was employed. By using the quasi static application of load in Abaqus, it was ensured that the static analysis does not get any effect of the inertial forces. The quasi static condition was guaranteed by the observation of negligible values of the kinetic energy over the complete solution until peak load.

Engineering definitions of (small) strains ɛ and stresses σ are used in this work. Since the displacement load u1 was applied at the reference node C1, far field strain was defined as the displacement of control node C1 divided by the initial length in the nominal fiber direction L₁ i.e. ɛ₁₁=u₁/L₁. The length L₁= 4 mm for models of first group whereas the length L₁= 6.657 mm for the model series used in the second group. To calculate percentage strains, they were multiplied by the factor 100. Similarly, far field stresses were calculated from reaction forces RF₁ at the node C1 divided by the initial cross sectional area perpendicular to the nominal fiber direction A₁ i.e. σ₁₁=|RF₁|/A₁. Compression stresses were defined positive in this work. The area A₁=2 × 1 and 4 × 2 mm² for 2D and 3D models of the preliminary series, respectively. The area A₁=1.1412 × 0.951 mm² for the models of probabilistic series. For the sake of brevity, the peak load for the axial compression load case was given a special symbol i.e., σ^c: = max(σ₁₁).

Results and discussions

Effect of the misalignment dimensionality on the compression strength prediction

Four models were considered in the preliminary series in order to analyze the effect of the fiber misalignment dimensionality on the strength prediction, see Figure 4. The resulting far field stress-strain behavior of all four models is plotted in Figure 9. Strength prediction using the Budiansky’s analytical model given in equation (1) depends on a single misalignment angle ψ. Moreover, equation (1) is consistent with the so-called infinite band assumption. The angle ψ is frequently identified with the largest misalignment in a given domain, i.e. ψ = 2.5° for the topologies considered in Figure 4. However, it will be shown in the following that this assumption is a simplification and it held only in case of the infinite band topology. Accepting this simplification of ψ for now and solving the maximization problem of equation (1) for ψ = 2.5° along with the consideration of the experimentally characterized nonlinear shear behavior of the material, Budiansky’s model gives a prediction of 721 MPa for compression strength. It is plotted as a dotted line in Figure 9. The result of the 2D model with a 1D misalignment (Figure 4(a)) coincides with the prediction of the analytical model, as it emulates the misalignment assumptions of the analytical model. Counter-intuitively, the 2D model with same maximum and minimum values of the misalignment angle, i.e. 2D model with a 2D undulation shown in Figure 4(b) results in an increase in strength with a value to 1051 MPa. Similarly, the 3D model with an in-plane 2D undulation and in-plane boundary conditions with a value of 1063 MPa shows an equivalent prediction to the 2D model with a 2D undulation. The 3D model with an in-plane 3D undulation and in-plane boundary conditions shows the strongest response among the four models with a strength value of 1326 MPa. The increase in strengths of the models with increasing misalignment dimensionality, albeit with the same absolute value of the maximum angle, can be understood by considering the statistics of the misalignment. The standard deviations s^θ of the misalignment angle from the models containing higher dimensional undulation are smaller, see Table 3. This is due to the fact that a region with different superposing waves contributing to the misalignment results in a dilution of misalignment in space. In the spatial sense, it means that only a few fibers are misaligned with a bigger portion of relatively straight fibers supporting the misaligned region.

Figure 9.

Comparison of results from 2D and 3D models of preliminary series, each containing a sinusoidal undulation with 1D, 2D and 3D misalignment characteristics.

This analysis shows that the use of an absolute maximum value of the misalignment angle alone for the prediction of compression strength has significant drawbacks. The statistics of the misalignment such as the standard deviation of the misalignment angle also need to be taken into account. Moreover, these results highlight the limitations of modeling the misalignment in FRPs in 2D, and point towards the need for realistic modeling. In real composite structures, the misalignment has 3D characteristics and a single angle based representation is inadequate. The 3D misalignment of fibers from a FRP based structure can be modeled using the in-plane θ and the out-of-plane ϕ angles as shown schematically in Figure 1. The measurements of the misalignment angles θ_k,l,m and ϕ_k,l,m shown in Figure 3(a) and (c) depict that not only the statistics of the misalignment but also topology needs to be taken into account. The spectral representation method was used to generate such topologies for probabilistic analyses, the results of which will be discussed in the next sections.

Effect of cutoff frequencies on the compression strength prediction

In order to investigate the effect of the cutoff frequencies ω_c1, ω_c2, and ω_c3 on the resulting numerical prediction of the compression strength of models containing realistic representation of the fiber misalignment, 250 realizations of each model from the probabilistic series were simulated. The resulting probabilities of survival of each model based on the functional forms of the spectral densities are plotted in Figure 10. Additionally, the probability of survival of the reference model discussed in reference topologies section, which was based on the spectral densities calculated from the experimental measurements of the fiber misalignment, is also plotted in Figure 10. At a sufficiently low load level, all the realizations survive. As the applied far field load σ₁₁ increases, the proportion of the failed realizations increases. At a sufficiently high value of σ₁₁, all realizations fail giving a P(σ^c) of 0. The mean and standard deviation values of all the models are provided in Table 4.

Figure 10.

Empirical probability of survival P(σ^c) over the compression strength σ^c for 3D models with cutoff frequencies and 3D models without cutoff frequencies. The survival probability using the reference model is also plotted for comparison.

Table 4.

Mean and standard deviation of the compression strength σ^c from different models of the probabilistic series. In addition, the shape m and scale R_o parameters for a 2-parameter Weibull distribution are also given as they will be needed for comparison to the experimental results in later sections.

Model	Mean[MPa]	s[MPa]	shape m[-]	Scale R_o[MPa]
Reference model	1424.54	61.96	25.77	1452.96
S_square based model no cutoff	1158.76	62.27	21.35	1187.39
S_square based model: with cutoff	1402.43	51.53	30.73	1426.51
S_exp based model: no cutoff	1164.54	62.30	20.83	1193.70
S_exp based model: with cutoff	1394.39	62.86	24.78	1423.43

A glance at the Figure 10 already shows that even though the mean square spectral densities Ξ^θ and Ξ^ϕ, and thus the standard deviation of the misalignment angle, of the models with and without cutoff are equal, the resulting strength predictions vary widely. The models without any cutoff frequencies have a sizable proportion of larger frequencies sampled from the functional forms of the spectral densities, which result in a sizable proportion of small wavelengths in each direction. It results in a large number of weak regions, which consequently result in a decrease of approximately 16–18% in the mean compression strength compared to the reference model based on CT scanned driven spectral densities. The models with the appropriate cutoff frequencies on the other hand show a probability of survival of compression strength very close to the reference model. This results highlights that in addition to the mean square spectral density of the misalignment, the frequency range of the misalignment also needs to be preserved in virtually generated models.

Looking at the range of the compression strength results using the corresponding standard deviation of each model’s strength distribution, some further observations can be made. All the models show standard deviations of predicted compression strength close to the reference model, except the model based on S_square with cutoff. Seemingly, when there is no cutoff, the results of the model series based on S_square and S_exp are almost identical. When appropriate cutoffs are applied, the mean strengths of both model series are still very close but values of the standard deviations and the Weibull parameter m vary significantly. The standard deviation and the Weibull parameter m values of the model series based on S_exp are almost identical to those of the reference model. But the standard deviation is lower and the parameter m value is higher of the model series based on S_square, showing a lower spread of the compression strength prediction. The reason behind this peculiar behavior can be understood in the discretized spectral densities. The experimentally measured spectral densities show larger values at lower frequencies. The model series based on S_exp emulates the observed allocation of higher intensity to lower frequencies. On the other hand, the model based on S_square discretizes the spectral density with equal weightage for all frequencies. This apparently translates into lower spread of the resulting distribution of the numerically predicted compression strength. It is concluded based on this observation that in addition to the statistical parameters of the distributions, the spectral parameters including the shape of spectral density curve need to be preserved for accurate prediction of the distribution of the compression strength.

Numerical predictions against experimental results using a scale law

For comparison of the model against experimental results, values of strength under axial compression load are taken from an earlier contribution of the authors, see Ref. ⁴⁷ The dimensions of the gauge section of the experimental specimens were 5 × 5 mm with a thickness of 1.15 mm. The experiments resulted in a mean strength of 1140 MPa and a standard deviation of 188.2 MPa. The volume of the CT-scans for measuring the fiber misalignment was limited due to the fine voxel resolution of 1.902 µm. The gauge section dimensions of the numerical models were hence limited by the available scan data for measuring the fiber misalignment to 6.657 × 1.1412 × 0.951 mm, see reference topologies section. Since the volume of the gauge section of the experimental specimens V_exp = 28.75 mm³ and the volume of the numerical models V_model = 7.2247 mm³ were significantly different, a direct comparison between numerical and experimental results cannot be made. Different experimental^51,52 as well as numerical^41,53 investigations have shown that the compression strength due to MB failure decreases with increasing sizes and vice versa in homogeneously loaded composites under axial compression load. This decrease in compression strength is attributed to the increasing number of the fiber misalignment induced defects in larger sizes based on the Weibull statistic driven weakest link model.⁵⁴ Similar change in strength with changing sizes has been observed in brittle materials such as ceramics.⁵⁵ To obtain a mathematical model for the size effect predictions, a power law proposed by Weibull is combined with the weakest link model resulting in following relation:

\frac{σ_{model}^{c}}{σ_{\exp}^{c}} = {(\frac{V_{model}}{V_{\exp}})}^{- 1 / m}

(5)

where

σ_{m o d e l}^{c}

and

σ_{\exp}^{c}

are strengths from model predictions and experimental results, and m is the shape parameter of the Weibull distribution. Therefore, numerical predictions for the reference model and the models based on assumed forms of spectral densities with cutoffs were compared against the experimental results using the Weibull scaling law given in equation (5) for the mean strengths, i.e. at the cumulative probability of failure of 0.5.

The values of the parameter m for the reference model and the models based on assumed forms of spectral densities with cutoffs are given in Table 4. The value of m is 7.11 for the experimental results given in Ref.⁴⁷ Figure 11 shows the aforementioned two parameter form of the scale law for all three model series as well as experiments. The inclination of the scale law is based on the respective m parameter. The scaling law appears as a straight line because of the double logarithmic scale in the figure. To aid the comparison, another variant of the scaling law for 0.5 cumulative probability of failure is included in the plot with a dotted line. For this variant, the line starts at the mean strength of the experimental results and the parameter m was adjusted so that it passes through the mean strength of the reference model. The corresponding value for m is 7.75 which is very close to the value 7.11 of the distribution fitted to the experimental results. As expected based on the Weibull scaling law, the results in Figure 11 indicate that there is a significant increase in the predicted mean strength with smaller model sizes compared to the experiments. Given the mean experimental strength, the value m = 7.75 for the Weibull based scale law can predict the numerical mean strength of the reference model quite accurately. The m values from different series of the numerical models on the other hand are quite high due to less scatter in corresponding probability distributions. The high m values from numerical models can be partially attributed to periodic boundary conditions limiting edge driven failure. Previous numerical investigations in the literature⁴¹ observed similar high values of the m parameter. The complete distributions of the numerically predicted strengths and the experimental strengths are also plotted using error bars.

Figure 11.

Comparison of σ^c from the reference model and the models based on assumed forms of spectral densities with cutoffs against the experimental results using the Weibull scaling law given in equation (5).

Effect of scaling mean square spectral densities on the prediction of the compression strength

In all previous investigations the standard deviation of misalignment for both the in-plane and the out-of-plane misalignment were held equal to the experimentally measured values given in Table 2. Consequently, the mean square spectral densities Ξ^θ and Ξ^ϕ were also same for different models. In this section, however, the effect of scaling mean square spectral densities on the prediction of the compression strength for same model dimensions (6.657 × 1.1412×0.951 mm) are considered by scaling Ξ^θ and Ξ^ϕ. For this purpose, three distinct models with each of them using both square and exponential forms of spectral densities were considered. For each of the two functional forms of the spectral densities, the values of Ξ^θ and Ξ^ϕ were scaled by factors of 0.25, 4, and 9, resulting in six distinct models. The relation between Ξ and s is given in equation (2) i.e., the resulting scaling of s^θ and s^ϕ was by factors of 0.5, 2, and 3. 250 realizations of each scaled model were simulated. The resulting probabilities of survival P(σ^c) of scaled model series are plotted in Figure 12 along with the unscaled models based on the average experimentally determined values of Ξ^θ and Ξ^ϕ. The mean and standard deviation values of all these model series are provided in Table 5.

Figure 12.

Empirical probability of survival P(σ^c) over the compression strength σ^c for models with cutoff frequencies and different scaled values of the mean square spectral densities Ξ^θ and Ξ^ϕ.

Table 5.

Mean and standard deviations of the compression strength σ^c for models based on functional forms of spectral densities with cutoff frequencies. Mean squares of the spectral densities were scaled to observe effects of resulting mean and standard deviation of numerically predicted distributions of the compression strength.

Scaling of mean square spectral density	σ^c of S_exp based models		σ^c of S_square based models
Scaling of mean square spectral density	Mean[MPa]	s[MPa]	Mean[MPa]	s[MPa]
0.25 ×Ξ^θ, 0.25 ×Ξ^ϕ	1964.93	76.65	1977.90	63.64
Ξ^θ, Ξ^ϕ	1394.39	62.86	1402.43	51.53
4 ×Ξ^θ, 4 ×Ξ^ϕ	967.21	42.93	967.83	35.84
9 ×Ξ^θ, 9 ×Ξ^ϕ	782.85	31.23	778.59	26.67

Solid lines in Figure 12 depict results of models based on S_square forms of spectral density whereas dotted lines illustrate results of models based on S_exp for both θ_k,l,m and ϕ_k,l,m. When the values of Ξ^θ and Ξ^ϕ were increased, the resulting distributions of σ^c were shifted towards lower values and vice versa. This is because increasing values of Ξ^θ and Ξ^ϕ means higher values of misalignment angles θ_k,l,m and ϕ_k,l,m, resulting in lower values of the compression strength σ^c.

Some further observations can be made when analyzing the range of the results using the corresponding means and standard deviations of distribution of compression strength from each model. The mean and standard deviation statistics are given in Table 5. Not only does the mean value of each model decreases with an increase of Ξ^θ and Ξ^ϕ, but the standard deviations of compression strengths also decrease for both S_exp and S_square models with increasing mean square misalignment. The larger the misalignment is, the more severe certain regions of the misalignment become, resulting in a smaller range of the distributions of σ^c. The analytical solution of Budiansky’s formula given in equation (1) indicates a hyperbola prediction of σ^c against misalignment. The slope of this hyperbola at larger angles is smaller compared to the steep slope at very small angles, which implies that the strength variations for models with larger s^θ and s^ϕ is smaller as predicted here. Interestingly, the standard deviation of S_square based model becomes comparable with that of the reference model when Ξ^θ and Ξ^ϕ are scaled by a factor of 0.25. It again highlights that the exponential form of the spectral density is more suitable in light of the measured spectral density for the current material.

Using the results of Figure 12, a relation for different percentiles of distribution of compression strength ^pσ^c(α) as a function of scaling factor α of Ξ^θ and Ξ^ϕ was defined. There, left super script p indicates the percentile of survival. The relation is applicable to all the percentiles of the distribution of the compression strength and is written as:

^{p} σ^{c} (α) =^{p} σ_{r e f}^{c} α^{a}

(6)

where a is defined by 3/4^th of the product of standard deviations of the misalignment of a reference model, i.e. a = 3/4 × s^θ × s^ϕ, see Table 2 for values of s^θ and s^ϕ. The symbol

^{p} σ_{r e f}^{c}

represents the p^th percentile of the distribution of the compression strength of a reference model. In current case,

^{p} σ_{r e f}^{c}

is associated with the percentile of the distribution of the compression strength of the unscaled S_exp and S_square based models. The S_exp and S_square based models can be scaled using the p values of the corresponding reference S_exp and S_square models, respectively. The relation given by equation (6) is plotted in a double logarithmic plot for both S_exp and S_square based models in Figure 13 for 1^st, 50^th (median), and 99^th percentiles of the distribution of σ^c. The values of 1^st, 50^th (median), and 99^th percentiles of the distribution of σ^c from each modeled case are given in Table 6. These values are also plotted as scatter points on top of the curves of equation (6) in Figure 13. The defined empirical relation fits quite well with the numerical results of the Ξ^θ and Ξ^ϕ scaled models for both S_exp and S_square based models. It shows that the strength decreases with a power law with an increasing scaling of Ξ^θ and Ξ^ϕ.

Figure 13.

1^st, 50^th (median), and 99^th percentiles of the distribution of the compression strength σ^c as a function of Ξ based on results of models with cutoff frequencies and different scaled values of the mean square spectral densities Ξ^θ and Ξ^ϕ. The corresponding values from each cases are plotted as scatter points on each curve.

Table 6.

1^st, 50^th (median), and 99^th percentiles of the distribution of the compression strength σ^c for different models.

Model parameters	Percentiles of distribution of σ^c
Model parameters	1^st[MPa]	50^th[MPa]	99^th[MPa]
S_exp, 0.25 ×Ξ^θ, 0.25 ×Ξ^ϕ	1767.12	1973.73	2103.31
S_exp, Ξ^θ, Ξ^ϕ	1223.35	1400.44	1520.04
S_exp, 4 ×Ξ^θ, 4 ×Ξ^ϕ	863.29	966.46	1057.73
S_exp, 9 ×Ξ^θ, 9 ×Ξ^ϕ	707.44	785.52	846.28
S_square, 0.25 ×Ξ^θ, 0.25 ×Ξ^ϕ	1821.97	1987.84	2103.28
S_square, Ξ^θ, Ξ^ϕ	1269.96	1407.20	1505.05
S_square, 4 ×Ξ^θ, 4 ×Ξ^ϕ	890.82	970.99	1044.81
S_square, 9 ×Ξ^θ, 9 ×Ξ^ϕ	720.86	779.93	.832.57

Conclusion

For best utilization of FRP based structures, the prediction of failure uncertainties and their derivative analyses is a fundamental requirement. Microbuckling (MB) is the dominant failure mode under axial compression load for unidirectional FRPs. The underlying uncertainties in MB failure lead to conservative design practices. This motivated the use of a numerical methodology with homogenized representation of fiber and matrix materials for quantifying the uncertainties in MB failure. Different aspects of the problem were investigated in this regard. Effect of the misalignment dimensionality on the compression strength prediction was investigated using simpler sinusoidal misalignment distributions in the homogenized approach. It was shown that the maximum misalignment angle alone is not the deciding factor but rather the whole distribution of the misalignment modeled in appropriate dimension plays a major role in prediction of compression strength.

The assumption of functional forms of the spectral densities for generation of virtual distributions of the fiber misalignment was investigated. Results of the models based on the different functional forms of the spectral densities were compared against results of a reference model informed by the experimental measurements of the fiber misalignment. The results imply that functional forms of the spectral densities can be assumed as an approximation. However, it was shown that proper cutoff frequencies need to be used for accurate prediction of compression strength. The effect of scaling mean square spectral densities on the prediction of the compression strength were also analyzed and an empirical relation for scaled spectral densities is provided. It was shown that the exponential form of spectral density gave better prediction for whole distribution of the compression strength for the material investigated in the current work.

It is expected that the presented numerical methodology and corresponding results can help in the further exploration of uncertainty quantification in FRPs. Eventually, it may enable more reliable design practices. Additionally, the proposed methodology offers an easy path to the probabilistic investigations of new fiber reinforced materials under compression dominated loads.

Footnotes

Acknowledgements

Funding by the Deutsche Forschungsgemeinschaft (DFG - German Research Foundation) under the project International Research and Training Group (IRTG1627) is gratefully acknowledged.

Declaration of conflicting interests

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

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Deutsche Forschungsgemeinschaft (IRTG1627).

ORCID iDs

Nabeel Safdar

Benedikt Daum

References

Soutis

. Fibre reinforced composites in aircraft construction. Prog Aerosp Sci 2005; 41: 143–151. DOI: 10.1016/j.paerosci.2005.02.004

Boeing . 787 Dreamliner by Design: Advanced Composites Use. Boeing.

Budiansky

Fleck

. Compressive kinking of fiber composites: a topical review. Appl Mech Rev 1994; 47(6): S246–S270.

Kaddour

Hinton

. Input data for test cases used in benchmarking triaxial failure theories of composites. J Comp Mat 2012; 46(19–20): 2295–2312. DOI: 10.1177/0021998312449886

Daum

Feld

Allix

, et al A review of computational modeling approaches to compressive failure in laminates. Comp Sci Tech 2019; 181. DOI: 10.1016/j.compscitech.2019.05.020

Rosen

. Mechanics of Composite Strengthening. NTRS. 1965, https://ntrs.nasa.gov/search.jsp?R=19660035520

Changfang

Zhitan

Changxing

, et al. Research on compression properties of unidirectional carbon fiber reinforced epoxy resin composite (UCFREP). J Comp Mat 2021; 55(11): 1447–1458, DOI: 10.1177/0021998320972176

Wang

Potter

Wisnom

, et al. Failure Mechanisms under Compression Loading in Composites with Designed Out-Of-Plane Fibre Waviness. Plastics, Rubber and Comp 2013; 42: 231–238.

Prabhakar

Waas

. Interaction between kinking and splitting in the compressive failure of unidirectional fiber reinforced laminated composites. Comp Struc 2013; 98: 85–92, DOI: 10.1016/j.compstruct.2012.11.005

10.

Budiansky

Fleck

. Compressive failure of fibre composites. J Mech Phy Sol 1993; 41(1): 183–211. DOI: 10.1016/0022-5096(93)90068-Q

11.