Multiscale Analysis of Mechanical Properties of 3D Orthogonal Woven Composites with Randomly Distributed Voids

Abstract

Voids are common defects in 3D woven composites because of the complicated manufacturing processes of the composites. In this study, a micro–meso multiscale analysis was conducted to evaluate the influence of voids on the mechanical properties of three-dimensional orthogonal woven composites. Statistical analysis was implemented to calculate the outputs of models under the different scales. A method is proposed to generate the reasonable mechanical properties of the microscale models considering randomly distributed voids and fiber filaments. The distributions of the generated properties agree well with the calculated results. These properties were utilized as inputs for the mesoscale models, in which void defects were also considered. The effects of these defects were calculated and investigated. The results indicate that tensile and shear strengths were more sensitive to the microscale voids, while the compressive strength was more influenced by mesoscale voids. The results of this study can provide a design basis for evaluating the quality of 3D woven composites with void defects.

1. Introduction

The past few decades have witnessed significant advances in the application of 3D woven composites owing to their superior mechanical properties, such as high specific stiffness/strength, high damage tolerance and high energy absorption capability [1,2,3]. The wide and extended use of these composites requires higher mechanical property standards. However, defects such as voids in the matrix are almost inevitable during the molding process because of the complex spatial structure of the woven composite. The presence of random void defects can significantly reduce the mechanical properties and cause the fluctuation of mechanical properties [4,5,6]. Moreover, the randomly distributed fiber filaments in microscale models also significantly influence the mechanical performance of composites [7]. Thus, this paper presents a micro–mesoscale analysis to evaluate the influence of randomly distributed voids on the mechanical behaviors of three-dimensional orthogonal woven composites (3DOWCs), considering randomly distributed fiber filaments in bundles.

Three-dimensional woven composites have complicated spatial structures, and their heterogeneity mainly reflects in micro and mesoscales. In the microscale model, the fiber bundles, also treated as unidirectional (UD) composites, of 3D woven composites consist of a matrix and thousands of randomly distributed fiber filaments [7,8,9]. In the mesoscale model, the weave architecture, orientations of different fiber bundles and the matrix are incorporated [10,11,12,13,14]. In both scales, the representative volume element (RVE) models are treated as periodic structures and utilized to calculate the mechanical properties by applying periodic boundary conditions (PBCs). Outputs from microscale models are used as inputs for mesoscale models. For example, Wang et al. [15] used a hexagonal-array microscale model to predict effective elastic properties of fiber bundles. The outcomes were used as inputs for predicting the moduli of a 3D orthogonal woven fiber-reinforced polymer matrix composite. Zhou et al. [16] conducted a progressive damage analysis for 2D plain woven composites. The authors utilized RVE models of hexagonally arranged fiber filaments in microscale analysis, calculated the elastic properties and strength properties, and then applied the results to mesoscale models. In recent years, an increasing number of researchers have used randomly distributed fiber models instead of regularly arranged models in microscale analysis. Representative volume elements with more fibers can capture the variation in strength, especially in transverse directions. Vajari [17,18] used randomly distributed models to study the influence of micro-voids in the matrix on biaxial strength properties and failure modes. Wang [19] studied the transverse tensile strengths of UD composites, considering thermal residual stress, using randomly distributed fiber models, and found variations due to the distributions. The arrangement of fibers is uncertain in advance. Randomly distributed fibers in micro-scale models will result in different elastic and strength properties. Thus, Monte-Carlo simulations are implemented in this paper to study the influence of randomness on the mechanical behaviors. To obtain a stable result, a process should be repeated numerously because of uncertainties [20,21,22]. In addition to the randomly distributed fibers, there are many other uncertainty factors, such as material properties, volume fractions, random cross sections and randomly distributed voids. Tao et al. [7,23] proposed a multiscale simulation to quantify the effects of numerous uncertainties on mechanical properties for 3DOWCs, including fiber and void distributions, modulus of the matrix, fiber volume contents and fiber bundle dimensions. Zhou et al. [24] coupled the multiscale method and perturbation to study the variability in the effective elastic properties of composites with random material properties. Shi et al. [14] used Weibull distribution material properties to conduct a multiscale damage analysis for 2.5D fiber-reinforced ceramic matrix composites. The results showed significant relationships existed between the constituent properties and the effective mechanical properties. Roham and Mohammad [25] conducted a progressive damage analysis for composite vessels, considering several uncertainties, including fiber volume fraction, winding angle and material strength properties. The importance of considering uncertainties is emphasized using statistical analysis. Guo et al. [26] proposed a random weft cross-section model based on experimental observations. The predicted results agreed well with the experimental results.

Void defect is also a common uncertainty generated from manufacturing processes, and it can mainly be classified into two types: microscale and mesoscale void defects. Recently, several researchers have focused on the effects of these defects on the mechanical performance of composites. For microscale defects, Dong [27] provided a practical method for predicting the effects of process-induced voids on the mechanical properties of carbon fiber/epoxy composites using the RVE model. Tensile and interlaminar shear strengths were calculated and compared with available experimental data. Jiang et al. [28] compared the axial stiffness and strength properties of a single-fiber bundle in UD composites with and without void defect based on a three-representative unit cell model. Carrera [29] used a 1D finite element model based on the Carrera unified formulation to evaluate the influence of voids on UD composites. Hyde et al. [30] conducted comprehensive research on the effects of micro-voids on the strength of UD composites. The effects of void shapes, void volume fractions and void orientations were considered. For mesoscale defects, Dong and Huo [31] studied the elastic properties of 3D braided composites with internal defects using a two-scale finite element model. Huang and Gong [32] used a similar method to predict the void effects on the effective elastic properties of 3DOWC. Gao et al. [13] recently predicted the mechanical properties of 3D braided composites with void defects, establishing a constitutive model and a finite element model of 3D braided composites. The results provided the design basis for evaluating the influence of void defects on mechanical behaviors.

An integral part of the multiscale finite element research is to build a microscale RVE that considers uncertainties. Monte-Carlo simulations are then conducted to obtain the statistical distributions of required mechanical properties. Some researchers used mean values as inputs for higher-scale analysis, while other researchers considered the distributions of these parameters and then generated distributed inputs for higher-scale analysis. However, most researchers ignored the relationships between these mechanical properties; the relationships are not stated in their papers. Thus, the main goal of this paper is to devise a multiscale approach to evaluate the influence of microscale and mesoscale void defects on the mechanical properties of 3DOWCs. The relationships between microscale properties are considered. The rest of this paper is organized as follows: Section 2 describes the multiscale models and boundary conditions. Section 3 introduces the constitutive models and failure criteria used in this paper. Section 4 provides the results and discussion on the multiscale analysis. Section 5 presents the concluding remarks and future work.

2. Multiscale Models and Methods

A 3D orthogonal woven fiber-reinforced polymer matrix composite was investigated in this study, and the material properties of the basic components are listed in Table 1. The 3DOWC had a one-by-one weave architecture. In the thickness direction, there were two layers in the warp direction and three layers in the weft direction. Figure 1 illustrates the framework of the multiscale stochastic model used in this study. Multiscale models with randomly distributed fiber filaments and voids were generated. Statistical analysis was conducted to obtain the correlations between the mechanical properties, and the results are described in Section 4. Then, based on the above correlations, new mechanical properties were generated as the input parameters for mesoscale models. Finally, statistical analysis was conducted to evaluate the effect of voids on the mechanical properties of 3DOWCs.

2.1. Microscale Model for UD Composites

In the microscale model, randomly distributed fibers were considered. Since the focus of this study is to evaluate the influence of void defects on mechanical properties, randomly distributed voids were also implemented in the matrix. Therefore, the material properties of fibers and matrix remained constant, and the randomly distributed fibers and voids were the variables. CATIA V5R21 (Vélizy-Villacoublay, France) was used to generate the randomly distributed geometry model of UD composites considering the periodicity, and then HyperMesh 14.0(Altair Engineering, Inc., Troy, MI, USA) was used for meshing. Corresponding nodes were on opposite surfaces to ensure the following application of PBCs. There are two ways to generate the void defects: One is considering the real geometry of voids and distributions and then establishing the finite element model; the other one is randomly selecting the matrix elements, putting them into one set and modifying the material properties of these elements. These two methods result in different void densities or numbers of voids, depending on the grid density. According to Vajari [17,18], increasing the number of matrix voids or void density in the model has no significant effect on the strengths or moduli but alters the crack path. In this study, the latter method was adopted, and its feasibility has been verified by many researchers [28,29,31,32]. To generate the random voids in the matrix, a Python (Python Software Foundation, Beaverton, OR, USA) script was used to modify the input file for ABAQUS 2020(Dassault Systemes Simulia Corp., Johnstone, RI, USA). Random elements were selected from the matrix set to the void set until the content of voids reached the set value. This way, 150 models were generated. Illustrations are shown in Figure 2. The fiber diameter was 7 $\mathsf{\mu }\mathrm{m}$, and the RVE model was 90 $\mathsf{\mu }\mathrm{m}$ in length and width and 2 $\mathsf{\mu }\mathrm{m}$ in thickness. Over one hundred fibers filaments are included in microscale model. Fibers, voids and matrix were meshed with eight-node linear hexahedral reduced integration elements (C3D8R). There were 84,519 elements and 113,888 nodes in this model. After the mesh sensitivity analysis, accurate mechanical properties could be obtained using the mesh size.

2.2. Mesoscale Model for 3DOWCs

The 3DOWC has a relatively complicated weave architecture, which may cause difficulty in the periodic meshing process. In this study, the textile geometry models of 3DOWC RVEs were generated using TexGen software (version 3.10.0, University of Nottingham, Nottingham, UK) [33]. Compared with the consistent meshing method that considers the real geometry of the fiber bundle section, the voxel meshing method has advantages in periodic meshing with a balance of efficiency and accuracy. Liu [34] studied the strength properties of 3DOWCs based on the voxel model, and the results agreed well with the experimental results. Owing to the efficiency and accuracy of this method, the voxel meshing method for the textile geometry model generated from TexGen was utilized in this study. Figure 3a,b present the mesoscale geometry model and voxel meshing for fiber bundles, respectively. The generation technique of voids in mesoscale model of 3DOWCs is the same with that in microscale model as is discussed in Section 2.1. The material orientations for each element of fiber bundles were automatically assigned by TexGen. The numbers of mesh seeds were 40, 80 and 32 in the warp, weft and thickness directions, respectively. Eight-node linear hexahedral reduced integration elements were assigned for fiber bundle and matrix components. A total of 102,400 elements and 109,593 nodes were included in this model.

2.3. Periodic Boundary Conditions

To satisfy the boundary displacement and boundary stress continuity of the RVE, PBCs should be employed to obtain accurate mechanical behavior. The PBCs should satisfy the following equations:

(1)$\{\begin{array}{l}{u}_i(a_1,x_2,x_3)-u_i(-a_1,x_2,x_3)=2a_1{\epsilon }_{i1}^0\begin{array}{c}-a_2\le x_2\le a_2\\ -a_3\le x_3\le a_3\end{array} \\ u_i(x_1,a_2,x_3)-u_i(x_1,-a_2,x_3)=2a_2{\epsilon }_{i2}^0\begin{array}{c}-a_1\le x_1\le a_1\\ -a_3\le x_3\le a_3\end{array}\\ u_i(x_1,x_2,a_3)-u_i(x_1,x_2,-a_3)=2a_3{\epsilon }_{i3}^0\begin{array}{c}-a_1\le x_1\le a_1\\ -a_2\le x_2\le a_2\end{array}\end{array}$

Here, $a_1$, $a_2$ and $a_3$ represent the RVE length in the first, second and third directions, respectively; ${\epsilon }_{i1}^0$ denotes the strain in the $i$ direction applied on the RVE surface in the first direction; and $u_i$ denotes the displacement in the $i$ direction.

The PBCs described in this paper were applied using nodes constraint equations by Python script in ABAQUS. To obtain the RVE mechanical properties, including moduli and strengths, separate PBCs should be considered. Figure 4 presents the program of applying the PBCs for the microscale model. In the first step, the effective moduli of the RVE are calculated in the post-processing stage using $\overline{\mathsf{\sigma }}={\displaystyle {\int }_{\mathsf{\Omega }}\frac{\mathsf{\sigma }dV}{V}}$. Here, $I_6$ is the unit matrix. Details for computing tensor $C$ and the related elastic moduli can refer to Barbero [35]. Then, PBCs for axial and biaxial strength analysis are applied with the Poisson effect. The normal and shear stresses of the RVE are computed from the resultant normal and tangential forces acting on the surfaces divided by the cross section at every sub-step in ABAQUS. Finally, the strength properties are extracted from $\sigma -\epsilon$ curves in the direction of interest. For the microscale RVE model, PBCs in all directions are applied, while for mesoscale analysis, PBCs in the thickness direction are ignored because of the full-thickness RVE model [36]. One corner for each RVE model is fixed in all directions to remove the rigid motion.

2.4. Distance Function for Generating Mechanical Properties of Microscale Model

The randomness of filament and void distributions will result in different mechanical properties for microscale models. In order to evaluate these parameters under certain void content, fitting functions are implemented to show the distributions of these mechanical properties. However, some of these parameters are positive correlated or negative correlated, see Section 4.1. Thus, it is inappropriate to directly use the fitting functions to generate the parameters as inputs for higher scale models.

To tackle this problem, a distance function is proposed here. After sufficient simulations, all elastic and strength properties are obtained. Thus, the scattered diagrams can be drawn between every two above parameters. As shown in Figure 5, we assume the fitting line between the two parameters is $y=a+bx$. Next, two points A and B which forms the error line should be found manually. It should be ensured that nearly all the scatter points are below the error line. Define the distance between the two lines at arbitrary point $x_1$ as $D\left(x_1\right)$. Then the relationship between $x$ and $y$ can be written as follows:

(2)$y=a+bx+rand(-D(x),D(x)),$

where x and y represent the mechanical parameters, a and b are the linear fitting parameters of the fitting line, and the random function $rand(c,d)$ returns a random value between $c$ and $d$. $D\left(x\right)$ is the distance function, which returns the linear distance value in y-direction between the linear fitting line and the error lines.

3. Failure Criteria and Constitutive Models

3.1. Failure Initiation Criteria

For the microscale model, the fiber component is regarded as a transversely isotropic material. In actual situations, no damage behavior is observed in transverse directions; thus, the maximum stress failure criteria are used to describe the failure behavior of fiber filaments because of their brittle damage behavior, as shown in Equation (3):

(3)$\begin{array}{l}{f}_{f1t}={\sigma }_{11}/X_{ft}\ge 1\\ f_{f1c}=\left|{\sigma }_{11}\right|/X_{fc}\ge 1^{\prime }\end{array}$

where $X_t$ and $X_c$ are failure stresses for fiber tensile and compressive failure modes, respectively, and ${\sigma }_{11}$ denotes the stress of the fiber component in the fiber longitudinal direction.

Matrix materials in fiber bundles (microscale) and between fiber bundles (mesoscale) were modeled as an isotropic material, and the modified Drucker–Prager yield model developed by Lubliner et al. [37] and Lee and Fenves [38] was applied to estimate the failure as follows:

(4)$f_m=\frac{1}{1-\alpha }\left(\sqrt{3J_2}+\alpha I_1+\beta ⟨{\sigma }_1⟩\right)-{\sigma }_{xc}\ge 0,$

where $I_1$ is the first invariant of the stress tensor; $J_2$ is the second invariant of the deviatoric stress tensor; ${\sigma }_1$ is the maximum principal stress; $⟨x⟩=\frac{x+\left|x\right|}{2}$; $\alpha$ is the pressure-sensitivity parameter and the value 0.13 is adopted in is paper; and $\beta$ is a function of tensile (${\sigma }_{xt}$) and compressive (${\sigma }_{xc}$) yield stress and is defined as follows:

(5)$\beta =\frac{{\sigma }_{xc}}{{\sigma }_{xt}}\left(1-\alpha \right)-\left(1+\alpha \right)$

The sudden stiffness degradation law is used to describe the damage behavior after failure initiation for the fiber filament and matrix components.

For the mesoscale model, the fiber bundles can be regarded as fiber-reinforced UD composites. Thus far, numerous failure criteria have been proposed to describe the failure initiation and damage evolution behavior for UD composites. Puck’s criterion, developed by Puck and Schürmann [39], ranks high among the several failure criteria in “World Wide Failure Exercise I and II” and was adopted in this study. Moreover, Gu and Chen [40] extended Puck’s inter-fiber failure (IFF) criterion, and the extended criterion was also considered in this study. Puck’s theory can be classified into two sets of equations: fiber failure (FF) and IFF.

The refined FF criteria were selected to identify the failure initiation for fiber damage under tensile and compressive loads as follows [41]:

(6)$f_{t(c)}=\left|\frac{{\sigma }_{11}}{X_{t(c)}}-(\frac{{\upsilon }_{12}}{X_{t(c)}}-{\upsilon }_{12f}\frac{E_1}{E_{1f}}{m}_{\sigma f})({\sigma }_{22}+{\sigma }_{33})\right|\ge 1$

where $X_t$ and $X_c$ are the longitudinal tensile and compressive strengths of UD composites, respectively; $E_1$ and $E_{1f}$ are the longitudinal moduli of UD composites and fiber, respectively; ${\upsilon }_{12}$ and ${\upsilon }_{12f}$ are the Poisson ratios of UD composites and fiber, respectively; and $m_{\sigma f}$ is a magnification factor for the matrix stress caused by the mismatch in the moduli of matrix and fibers; in this study, $m_{\sigma f}=1.1$ was adopted, as suggested in [42,43].

Puck’s IFF criterion is based on the Mohr–Coulomb theory, and failure will occur on the fracture plane where only shear stresses ${\tau }_{nt}$ and ${\tau }_{nl}$ and normal stress ${\sigma }_n$ exist, as shown in Figure 6. The IFF criterion is classified into two modes: matrix tension (${\sigma }_n(\theta )\ge 0$) and matrix compression (${\sigma }_n(\theta )<0$), which can be introduced as follows:

(7)$f_{IFF}=\{\begin{array}{ll}\sqrt{{\left(\frac{1}{R_{\perp }^{A+}}-\frac{p_{\perp \psi }^{(+)}}{R_{\perp \psi }^A}\right)}^2{\sigma }_n^2(\theta )+{\left(\frac{{\tau }_{nt}(\theta )}{R_{\perp \perp }^A}\right)}^2+{\left(\frac{{\tau }_{nl}(\theta )}{R_{\perp \parallel }}\right)}^2}+\frac{p_{\perp \psi }^{(+)}}{R_{\perp \psi }^A}{\sigma }_n(\theta )& \mathrm{for} {\sigma }_n(\theta )\ge 0\\ \sqrt{{\left(\frac{p_{\perp \psi }^{(-)}}{R_{\perp \psi }^A}{\sigma }_n(\theta )\right)}^2+{\left(\frac{{\tau }_{nt}(\theta )}{R_{\perp \perp }^A}\right)}^2+{\left(\frac{{\tau }_{nl}(\theta )}{R_{\perp \parallel }}\right)}^2}+\frac{p_{\perp \psi }^{(-)}}{R_{\perp \psi }^A}{\sigma }_n(\theta )& \mathrm{for} {\sigma }_n(\theta )<0\end{array}$

With:

(8)$\begin{array}{l}\frac{p_{\perp \psi }^{(+,-)}}{R_{\perp \psi }^A}=\frac{p_{\perp \perp }^{(+,-)}}{R_{\perp \perp }^A}{\cos }^2\phi +\frac{p_{\perp \parallel }^{(+,-)}}{R_{\perp \parallel }^A}{\sin }^2\phi ’{\cos }^2\phi =\frac{{\tau }_{nt}^2(\theta )}{{\tau }_{nt}^2(\theta )+{\tau }_{nl}^2(\theta )}\\ R_{\perp \perp }^A=\frac{Y_c}{2(1+p_{\perp \perp }^{-})},R_{\perp \parallel }^A=S_{21}\end{array}$

(9)$\begin{array}{l}{\sigma }_n(\theta )={\sigma }_2{\cos }^2\theta +{\sigma }_3{\sin }^2\theta +2{\tau }_{23}\sin \theta \cos \theta \\ {\tau }_{nt}(\theta )=({\sigma }_3-{\sigma }_2)\sin \theta \cos \theta +{\tau }_{23}({\cos }^2\theta -{\sin }^2\theta )\\ {\tau }_{nl}(\theta )={\tau }_{12}\cos \theta +{\tau }_{13}\sin \theta \end{array}$

The seven parameters $R_{\perp \perp }^A,R_{\perp \parallel }^A,R_{\perp }^{A+},p_{\perp \parallel }^{(+)},p_{\perp \parallel }^{(-)},p_{\perp \perp }^{(+)},p_{\perp \perp }^{(-)}$ need to be determined before Puck’s criterion is implemented. According to Gu and Chen [40], UD composites are divided into three categories, resulting in different formulas to obtain the above parameters. Details can be found in ref. [40].

In Equation (7), $f_{IFF}$ is a function of the potential fracture angle $\theta$, namely the fracture plane. The plane with the highest $f_{IFF}$ is the actual fracture plane. Thus, for an arbitrary stress state and uncertain material properties, the potential fracture angle should be determined before the failure onset is predicted. However, it is extremely difficult to obtain the analytical solution of fracture angle for such a 3D problem. Thus, the semi-analytical procedure is utilized. In this study, an extension and combination algorithm of the selective range golden section search (SRGSS) [44] and the semi-analytical algorithm (SAA) [45] was developed to accurately determine the fracture angle. The algorithm implemented in this paper is more accurate than the SRGSS and SAA. However, a detailed description of the algorithm is beyond the scope of this paper and is thus not presented here. Once the fracture angle is obtained, the global maximum value $f_{IFF}$ is calculated, and the matrix damage is initiated when $f_{IFF}$ reaches 1.

Once the fiber damage and matrix damage starts, damage evolution processes begin and the equivalent stresses and strains of the fiber and matrix are defined as follows [46]:

(10)$\begin{array}{l}{\sigma }_{eq}^f=\sqrt{{\sigma }_{11}^2+{\sigma }_{22}^2+{\sigma }_{33}^2}, {\sigma }_{eq}^m=\sqrt{{⟨{\sigma }_n⟩}^2+{\tau }_{nt}^2+{\tau }_{nl}^2}\\ {\epsilon }_{eq}^f=\sqrt{{\epsilon }_{11}^2+{\epsilon }_{22}^2+{\epsilon }_{33}^2}, {\epsilon }_{eq}^m=\sqrt{{⟨{\epsilon }_n⟩}^2+{\epsilon }_{nt}^2+{\epsilon }_{nl}^2}\end{array}$

With:

(11)$\begin{array}{l}{\epsilon }_n={\epsilon }_2{\cos }^2\theta +{\epsilon }_3{\sin }^2\theta +2{\epsilon }_{23}\sin \theta \cos \theta \\ {\epsilon }_{nt}=({\epsilon }_3-{\epsilon }_2)\sin \theta \cos \theta +{\epsilon }_{23}({\cos }^2\theta -{\sin }^2\theta )\\ {\epsilon }_{nl}={\epsilon }_{13}\sin \theta +{\epsilon }_{13}\cos \theta \end{array}$

With ${\epsilon }_{eq}^0$ and ${\epsilon }_{eq}^f$ as the equivalent strains of damage initiation and complete damage, respectively, the equivalent strain ${\epsilon }_{eq}^0$ for FF and IFF can be written as:

(12)${\epsilon }_{eq}^{f0}=\frac{{\epsilon }_{eq}^f}{f_{t(c)}},$

(13)${\epsilon }_{eq}^{m0}=\frac{{\epsilon }_{eq}^m}{f_{IFF}}.$

As for the equivalent strain of complete damage for IFF, the energy-based damage law under mixed-mode is built with the element characteristic length to alleviate mesh dependency during material softening:

(14)${\left(\frac{g_{n}L}{G_{mt(c)}}\right)}^{\xi }+{\left(\frac{g_{nt}L}{G_{23c}}\right)}^{\xi }+{\left(\frac{g_{nl}L}{G_{12c}}\right)}^{\xi }=1’$

where $g_n,g_{nt}$ and $g_{nl}$ are the strain energy density of the corresponding stress component, respectively; $G_{mt(c)}$ is the energy dissipated during damage for transverse tensile and compressive directions, respectively; $G_{23c}$ and $G_{12c}$ are in-plane and out-of-plane dissipated energies for the corresponding directions; $L$ is the characteristic length; $\xi$ is a material parameter, set as 2 [46] in this paper.

The strain energy density associated with each effective stress component due to complete damage is defined as:

(15)$g_j^f={\displaystyle {\int }_0^{{\epsilon }_j^f}{\sigma }_j}d{\epsilon }_j\approx \frac{1}{2}{\sigma }_j^0{\epsilon }_j^f=\frac{1}{2}{\sigma }_j^0\beta {\epsilon }_{eq}^{mf}, j=n,nt,nl,$

where ${\sigma }_j^0$ is the stress component on the fracture plane when the damage is initiated, and ${\beta }_j$ represents the mix-mode ratio and can be expressed as:

(16)${\beta }_n=\frac{⟨{\epsilon }_n⟩}{{\epsilon }_{eq}^m}, {\beta }_{nt}=\frac{{\epsilon }_{nt}}{{\epsilon }_{eq}^m}, {\beta }_{nl}=\frac{{\epsilon }_{nl}}{{\epsilon }_{eq}^m}.$

Thus, the equivalent failure strain ${\epsilon }_{eq}^{mf}$ can be formulated by substituting Equations (15) and (16) into Equation (14):

(17)${\epsilon }_{eq}^{mf}=\frac{2}{L}{\left[{\left(\frac{{\sigma }_n^0{\beta }_n}{G_{mt(c)}}\right)}^{\xi }+{\left(\frac{{\sigma }_{nt}^0{\beta }_{nt}}{G_{23c}}\right)}^{\xi }+{\left(\frac{{\sigma }_{nl}^0{\beta }_{nl}}{G_{12c}}\right)}^{\xi }\right]}^{-1/\xi }$

The IFF damage variable is given by the following expression based on the bilinear model:

(18)$d_m=\frac{{\epsilon }_{eq}^{mf}\left({\epsilon }_{eq}^m-{\epsilon }_{eq}^{m0}\right)}{{\epsilon }_{eq}^m\left({\epsilon }_{eq}^{mf}-{\epsilon }_{eq}^{m0}\right)}$

Similarly, the equivalent strain of complete damage for FF can be obtained through the same procedure.

In the damage evolution process, the generalization of the Duvaut–Lions regularization model [47] is adopted for the mesoscale model to improve the convergence of the numerical calculation and smoothen the stiffness degradation. Thus, the viscous damage variable is defined as follows:

(19)$\stackrel{·}{d_i^v}=\frac{1}{\eta }(d_i-d_i^v)$

3.2. Damage Constitutive Model

In this paper, the damage modes for fiber bundles are mainly classified into four categories: longitudinal tensile damage $D_{ft}$, longitudinal compressive damage $D_{fc}$, transverse tensile damage $D_{mt}$ and transverse compressive damage $D_{mc}$. When the material is damaged, its stiffness is degraded. If the damage occurs only in the $y$ direction (${\theta }_{fr}=0$), the material can still bear load in the $z$ direction. The degraded stiffness matrix can therefore be expressed as follows:

(20)$C_D=\left[\begin{array}{cccccc}{C}_{11}(1-D_f)& C_{12}(1-D_f)A_1{A}_4& C_{13}(1-D_f)A_2{A}_3& & & \\ & C_{22}{A}_1{A}_2& C_{23}(1-D_{mt})(1-D_{mc})& & 0& \\ & & C_{33}{A}_3{A}_4& & & \\ & & & C_{44}(1-D_f)A_1{A}_2& & \\ & SYM& & & C_{55}(1-D_f)A_3{A}_4& \\ & & & & & C_{66}(1-D_{mt})(1-D_{mc})\end{array}\right],$

where, $D_{\mathrm{f}}=1-(1-D_{ft})(1-D_{fc})$, $A_1=1-D_{mt}\cos {\theta }_{fr}$, $A_2=1-D_{mc}\cos {\theta }_{fr}$, $A_3=1-D_{mt}\sin {\theta }_{fr}$ and $A_4=1-D_{mc}\sin {\theta }_{fr}$.

Thus, the stresses are calculated using the following equation:

(21)$\mathsf{\sigma }=C_D:\mathsf{\epsilon }$

The Jacobian matrix can be formulated by differentiating Equation (21) as follows:

(22)$\frac{\partial \Delta \mathsf{\sigma }}{\partial \Delta \mathsf{\epsilon }}=C_D+\left({\displaystyle \sum _i\frac{\partial C_D}{\partial D_i}\frac{\partial D_i}{\partial {\epsilon }_{eq}}\frac{\partial {\epsilon }_{eq}}{\partial \epsilon }}\right):\mathsf{\epsilon }\left(i=ft,fc,mt,mc\right)$

4. Results and Discussion

4.1. Microscale Model

Figure 7 displays the typical stress–strain curves of UD composites with 3% void content under different uniaxial loadings. The modulus, strength and energy release rate properties were all obtained from the curves. In this study, 150 randomly distributed microscale models were established. For each non-zero void content (1%, 3% and 5%), three models with randomly distributed voids were established for each microscale model. Therefore, there were a total of 150 random models without voids and 1350 random models with voids. Figure 8 gives the maximum principal stress plots for one microscale model with different void contents under transverse tensile loads. The fiber component is hidden in the figure. The presence of voids does not affect the general stress distribution but will cause discontinuities of the distribution. Meanwhile, a higher stress concentration occurs between two adjacent fiber filaments in the load direction, at which the matrix elements will be damaged compared with the void-free case. Figure 9 presents the counts distribution histograms of $E_1$ and $E_2$ of the UD composite with various void contents, and Gaussian distributions are implemented to fit the probability density functions, which will be utilized for generating mechanical properties in the mesoscale models. It should be noted when the Gaussian distribution is implemented in mesoscale model, the generated values should be in the minimum-maximum range, obtained from finite element simulations, in order to avoid the appearance of too large or too small values. The statistical parameters of Gaussian distributions for all mechanical properties are listed in Appendix A (Table A1).

Even though all the mechanical parameters follow the Gaussian distributions, they cannot be directly utilized as inputs for the mesoscale analysis because there are correlations among these parameters. For instance, Figure 10 displays the relations between $\left(G_{12},E_2\right)$, $\left(G_{23},E_2\right)$,$\left({\sigma }_{yt},E_2\right)$,$\left({\sigma }_{yc},{\sigma }_{yt}\right)$,$\left(G_{IC}^{yt},{\sigma }_{yt}\right)$ and $\left(G_{IC}^{23},{\tau }_{23s}\right)$ when the void content is 3%. In Figure 10a, the blue scattered bubbles denote the distributions of transverse modulus and longitudinal shear modulus. The solid line denotes the linear fitting curve, and the dashed error line is used to measure the deviation, which is also reflected in the bubble color intensity. The further away from the fitting curve, the lighter the bubble color. Upward and downward trends can be observed; therefore, the mechanical parameters cannot be directly and separately obtained from Gaussian distributions in the case of unreasonable situations such as UD composites with high transverse tensile strength and low transverse compressive strength.

Owing to the above reasons, Gaussian distributions were utilized for parts of the mechanical properties $E_1,E_2,{\sigma }_{xt},{\sigma }_{xc},G_{IC}^{xt}$ and $G_{IC}^{xc}$ as basic inputs for mesoscale analysis. Other mechanical properties were calculated based on the above parameters using Equation (2).

To predict the mechanical properties for unknown microscale models, Python scripts were utilized to generate all the microscale properties using the abovementioned method. Figure 11a,b displays the distributions of calculated results and generated results. Figure 11c gives the box plots of the distributions of mechanical properties. From the figures, the distributions of the generated mechanical properties agree well with the numerical results. Hence, the method proposed in this paper is accurate enough for predicting the mechanical properties considering randomly distributed fibers and void defects.

In the next scale analysis, the UD composites will be regarded as transverse isotropic materials, and Puck’s criterion will be adopted to predict the failure onset. Thus, biaxial loading simulations for UD composites were conducted to verify the applicability of Puck’s criterion for predicting the strength properties under uncertain models and complex loading conditions. Figure 12 displays the failure envelopes in ${\sigma }_2-{\tau }_{12}$ stress space under different void contents. In this case, 10 random models were selected for each void content case, and 12 stress ratios were considered. The numerical results agree well with the analytical results, verifying the applicability of Puck’s criterion.

4.2. Mesoscale Model

The elastic properties and strength properties of 3DOWC mesoscale models in both warp and weft directions and in-plane shear properties were analyzed considering the randomly distributed voids in fiber bundles and between bundles. For the mesoscale models, the bundles were homogenized from the microscale model. Stochastic properties of bundles were generated based on the microscale models discussed above. Considering that the voids were randomly distributed in fiber bundles and that the fiber distribution is different at different places, the properties of fiber bundle elements were different from each other. In this study, we assumed that the distributions of voids and fibers differed among the bundle elements, which results in different mechanical properties. Thus, each bundle element had its own sections, and each section had unique mechanical properties generated from the constructed microscale models in the given void content case. Figure 13 shows the random material RVE model. Different color represents different materials or sections. A total of 50 groups of models were established for each load case to characterize the stochasticity of the mechanical properties. The above process was conducted using Python script by modifying the input files generated from ABAQUS.

Figure 14 depicts the typical stress–strain curves of 3DOWCs with no voids under different loading directions. From the image, the initial slopes of tensile and compressive stress–strain curves were consistent; therefore, the elastic moduli were obtained from the tensile stress–strain curves. Figure 15 presents the statistical results for elastic properties and strength properties in the warp, weft and in-plane shear directions. The void contents in fiber bundles were set as 0%, 1%, 3% and 5%, which are the same as those of the microscale models. Void contents of 0% and 2% were considered between bundles. The fiber volume content remained constant in all models. From the figures, the modulus and strength properties all exhibited a decreasing trend as the void content increased. Figure 15a–c illustrates the relationships between the modulus and void content. Even though the presence of voids reduced the moduli in the warp, weft and in-plane shear directions, the decreasing amplitude was unnoticeable. Moreover, the variance of the moduli was not significant. Taking the void-free case as an example, the minimum and maximum moduli in the warp direction were 33,679.79 MPa and 33,687.56 MPa, respectively. In the microscale model, the standard deviation of the longitudinal modulus was small, while that of the transverse modulus was relatively large. The void had nearly no influence on the longitudinal modulus, while the transverse modulus reduced by 9.35% when the void content was increased to 5%. In the mesoscale model, the 3DOWCs had longitudinal fiber tows in both the warp and weft directions, and the longitudinal fiber tows were the main load-bearing components. Different from the results for modulus, the models with randomly distributed element properties showed significant variation in strength. With the uneven materials distributed in fiber bundles, some elements had high tensile or compressive strengths, while others had lower ones. Low-strength properties resulted in the early failure of the elements. Evaluation using the Puck’s criterion revealed that the elements exhibited reduced stiffness, lost the load-bearing capacity and featured stress concentration, which may further propagate the damage process. However, the elements with low strength properties may not be located at the key positions or in the load directions; if so, the overall strengths of the 3DOWCs will be relatively high. Moreover, the variations in mechanical properties, including moduli and strengths, were generally remarkable when voids were randomly distributed in the matrix between bundles because of the inclusion of one more uncertainty source.

From Figure 15, the randomly distributed voids in fiber bundles and between bundles caused different degrees of reduction in modulus and strength. To further evaluate the influence of voids on these mechanical properties, the reductions were calculated using the mean values. The total content of fiber bundles in the mesoscale model was 55.56%; thus, the void content in the microscale model should be translated to the mesoscale model by multiplying the void content by 0.5556. Assume $V_1$ is the void content in fiber bundle, then $0.5556V_1$ denotes the void content in the fiber bundles relative to the whole 3DOWC model. $V_2$ means the void content between fiber bundles and it is relative to the whole composite. For the convenience of subsequent descriptions, solid, dashed, black, red and blue lines or symbols are used to describe the various cases. The solid lines and dashed lines denote the reduction and reduction efficiency of the corresponding mechanical parameter at different void content cases. Black symbols describe the $V_2=0$ cases, and the baseline is $\left(V_1,V_2\right)=(0,0)$. It denotes there is no void between fiber bundles and only voids inside the bundles are considered. Red symbols describe the $V_2=2\%$ cases, and the baseline is $\left(V_1,V_2\right)=(0,2\%)$. Blue symbols denote the $V_2=2\%$ cases, but the baseline is $\left(V_1,V_2\right)=(0,0)$. The difference between red and blue symbols is that they have different baselines. In Figure 16a–c, the square and circle solid lines almost coincide, which demonstrates that the voids between bundles did not influence the effect of voids in bundles on the elastic properties; that is, the modulus decreased to the same extent as the increase in the void content inside the fiber bundle for different $V_2$ cases. The decrease rate increased linearly with the increase in the void content inside the fiber bundles. In Figure 16, the blue solid line is almost parallel to the red and black solid line; this indicates that the reduction in modulus caused by the increase in the void content between fiber bundles was the same when the void content inside the fiber bundle was constant. The reduction efficiency is the corresponding mechanical parameter reduced per unit void content compared with the baseline. The void content here represents the total voids in 3DOWC including voids inside and between fiber bundles. Even though the pores between and inside the fiber bundles exhibited a linear reduction in the modulus, the efficiencies of the modulus reduction for both void contents were not the same. As illustrated in the figure, the black and red dashed lines do not differ much for a certain $V_2$, but the blue dashed line is significantly smaller than the black and red dashed lines, indicating that the effect of voids between fiber bundles on the 3DOWC modulus was smaller than the effect of pores inside fiber bundles with the same content on the 3DOWC modulus. Moreover, the blue dashed line tends to rise more slowly. The effect of the two types of voids on the strength properties was much more complicated than that on the modulus properties. For the tensile strength, ${\sigma }_{xt}$ increased with $V_1$ at a certain $V_2$, but the increase tended to decrease after 3% void content. As seen from the dashed line, the black line had the highest reduction efficiency, and the reduction in strength due to voids became less significant with the appearance of microscale void defects, indicating that the voids inside the fiber bundle had a greater effect on the 3DOWC tensile strength, and the internal defects were more likely to cause a reduction in the 3DOWC strength. However, for ${\sigma }_{xc}$, the blue curve has the highest reduction efficiency, which indicates that the voids between fiber bundles had a greater effect on the overall compressive strength under compression in the warp direction. Similar conclusions can be reached for ${\sigma }_{yc}$. According to [48,49], the tensile strength of 3D woven composites is strongly influenced by the fiber bundle properties, while the compressive properties are largely controlled by the matrix properties. In this study, the internal voids of fiber bundles caused the reduction of the properties of 3DOWC fiber bundle components, while the mesoscale voids caused the reduction of the properties of 3DOWC matrix components. These findings are consistent with the results in the literature [48,49] that intra-fiber bundle voids have a greater effect on tensile and shear properties, while the mesoscale defects have a greater effect on compressive strength properties. In the current study, for shear strength properties, the decrease in ${\tau }_{s12}$ with increasing $V_1$ was the same for the same $V_2$ and was approximately linear in the given range of void contents. At the same $V_2$, the ${\tau }_{s12}$ reduction efficiency tended to decrease with increasing $V_1$ and eventually tended to level off. At $V_2=2\%$, the overall ${\tau }_{s12}$ reduction efficiency was significantly lower than that at $V_2=0$, but the reduction efficiency featured a growing trend and eventually leveled off, which is inconsistent with the trend of the black line, indicating that coupling effects that affect the shear strength performance existed between the two kinds of voids.

In general, the void content inside the fiber bundle and between the fiber bundles both exhibited a roughly linear decreasing trend on the modulus, and the effect of the voids inside fiber bundles was higher than that of the voids between the fiber bundles. The two defects seemed to be independent of each other, with no significant coupling effect. Regarding strength properties, coupling relationships existed between the effects caused by the two defects, and the presence of voids between fiber bundles changed the pattern of the effect of voids inside fiber bundles on strength properties. Moreover, the internal voids had a greater effect on tensile and shear properties, but the compressive strength was more influenced by the mesoscale voids.

5. Conclusions

In this paper, a multiscale analysis is proposed to evaluate the influence of micro- and meso-void defects on the elastic properties and strength properties of 3DOWCs. Randomly distributed fiber filaments and voids were considered in the microscale model. The outputs from lower-scale models were utilized as inputs for higher-scale models. Non-uniform material properties were assigned to the fiber bundle elements of 3DOWCs. Statistical results are presented to provide an insight into the effect of voids on the mechanical performance of the composites.

In the microscale analysis, 150 randomly distributed fiber filaments models were established, and three models with randomly distributed voids were generated for each void-free model and each void content. Mean values and standard deviations were calculated to provide intuitive results of mechanical properties. The relationships between the two properties were obtained by plotting the bubble figures between any two parameters. The mean values and dispersity of the material properties could be considered to describe the relationships through linear fitting and the generation of random errors. The proposed method ensures the reasonableness of newly generated material properties. The generated parameters agree well with the original calculated results.

In the microscale analysis, a voxel finite element model for 3DOWC was established. Generated material properties using the previous method were assigned for the fiber bundle properties as uncertainty sources from microscale models. Another uncertainty source was the void defects between bundles. Extended Puck’s criterion was implemented to predict the failure initiation, and an energy-based damage evolution model was adopted for the progressive damage processes. The elastic properties and strength properties in the warp/weft tension/compression and in-plane shear directions were calculated. The results indicate the following: (1) the void defects will reduce the elastic properties, but the variations are not sensitive to the uncertainties; (2) the strength properties are sensitive to the uncertainties caused by the non-uniform distributed material properties of fiber bundles; (3) the elastic properties roughly linearly decreased with the void contents, including voids in and between fiber bundles; (4) coupling effects that affect the strength properties occurred between the two kinds of void defects; and (5) tensile and shear strengths were sensitive to the voids inside the bundles, while the compressive strength was sensitive to the voids between bundles.

The present research can provide the design basis for evaluating the influence of two kinds of void defects on 3DOWC mechanical properties. In future work, a similar procedure will be conducted for macroscale models, such as 3D woven structures, considering not only the present uncertainties but also the bundle section uncertainties.

Author Contributions

Conceptualization, Y.G., T.H. and S.Z.; software, Y.G. and Y.S.; methodology, Y.G., Y.S. and S.Z.; writing—original draft, Y.G.; investigation, resources, X.Z. and P.J.; supervision, X.Z.; writing—review and editing, T.H. and S.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China [NO. 51475372].

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Acknowledgments

The authors would like to thank all the companions in Office 415, School of Mechanics, Civil Engineering and Architecture, Northwestern Polytechnical University, for discussions and advices during the preparation of this paper.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A

Table A1

Mean value and standard deviation of all mechanical properties for micro-scale models.

	Void	Mean	Std		Void	Mean	Std		Void	Mean	Std
$E_1$ (MPa)	0	138,553.96	5.304	${\sigma }_{xt}$ (MPa)	0	2022.16	0.251	$G_{IC}^{xt}$ (N/mm)	0	15.16	5.8 × 10−4
	1%	138,513.33	5.4		1%	2005.18	2.664		1%	15.1	0.0236
	3%	138,433.2	5.583		3%	2002.2	0.158		3%	14.99	0.02594
	5%	138,353.58	5.736		5%	2001.07	0.243		5%	14.9	0.02473
$E_2$ (MPa)	0	7558.85	92.336	${\sigma }_{xc}$ (MPa)	0	1503.3	0.18	$G_{IC}^{xc}$ (N/mm)	0	8.256	2.81 × 10−4
	1%	7423.79	95.367		1%	1502.14	0.365		1%	8.245	0.00364
	3%	7149.63	96.21		3%	1500.29	1.807		3%	8.238	0.00876
	5%	6852.07	109.8		5%	1497.43	2.274		5%	8.235	0.01295
${\upsilon }_{12}$	0	0.318	4.038 × 10−4	${\sigma }_{yt}$ (MPa)	0	64.07	2.053	$G_{IC}^{yt}$ (N/mm)	0	0.426	0.0358
	1%	0.3172	4.007 × 10−4		1%	57.41	4.512		1%	0.394	0.0352
	3%	0.3156	4.019 × 10−4		3%	49.17	4.142		3%	0.346	0.0344
	5%	0.3139	4.454 × 10−4		5%	43.69	3.505		5%	0.315	0.032
${\upsilon }_{23}$	0	0.4755	0.00602	${\sigma }_{yc}$ (MPa)	0	147.2	5.299	$G_{IC}^{yc}$ (N/mm)	0	2.236	0.1241
	1%	0.4728	0.00619		1%	131.78	11.065		1%	1.809	0.1482
	3%	0.4664	0.00709		3%	112.3	10.002		3%	1.777	0.141
	5%	0.4597	0.00871		5%	99.53	8.314		5%	1.585	0.1409
$G_{12}$ (MPa)	0	3876.98	73.309	${\tau }_{12s}$ (MPa)	0	50.91	2.637	$G_{IC}^{12}$ (N/mm)	0	0.581	0.0548
	1%	3790.19	76.972		1%	48.92	2.695		1%	0.554	0.0549
	3%	3621.21	75.79		3%	45.54	2.705		3%	0.52	0.0493
	5%	3445.74	77.938		5%	41.9	2.389		5%	0.479	0.0475
$G_{23}$ (MPa)	0	2528.02	38.013	${\tau }_{23s}$ (MPa)	0	60.06	4.352	$G_{IC}^{23}$ (N/mm)	0	1.381	0.1438
	1%	2486.08	38.284		1%	55.17	4.032		1%	1.253	0.1214
	3%	2396.78	40.006		3%	50.35	3.117		3%	1.135	0.099
	5%	2306.22	47.622		5%	46.08	3.126		5%	1.026	0.0956

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Multiscale stochastic modeling framework.

Figure 2. Microscale finite element model of randomly distributed fibers with void contents: (a) 0%, (b) 1%, (c) 3% and (d) 5%.

Figure 3. Textile geometry (a) and voxel mesh (b) of the three-dimensional orthogonal woven composites.

Figure 4. Program of applying periodic boundary conditions for determining elastic properties and strength properties.

Figure 5. Diagram for determining the relationship between x and y.

Figure 6. Stress components acting on potential fracture plane. (Blue stresses are in global coordinate system and red stresses are in local coordinate system.)

Figure 7. Typical stress–strain curves of unidirectional composite (3% void content) under different uniaxial loadings.

Figure 8. Maximum principal stress plots of microscale matrix component under transverse tensile load: (a) 0%, (b) 1%, (c) 3% and (d) 5%.

Figure 9. Counts distributions of E1 and E2 under different void contents.

Figure 10. Relationships between the mechanical properties of microscale model with 3% void content.

Figure 11. Distributions of calculated and generated mechanical properties. (a,b) represent the examples of calculated and generated results; (c) shows the distributions of two mechanical properties.

Figure 12. Failure envelopes of microscale model in σ2−τ12 stress space under different void contents between calculated and theoretical results.

Figure 13. Representative volume element of mesoscale model with random material properties generated from microscale outputs.

Figure 14. Typical stress–strain curves of three-dimensional orthogonal woven composite with no voids under different uniaxial loadings.

Figure 15. Statistical elastic properties and strength properties of mesoscale model in warp, weft and in-plane shear directions. (a–c) represent the elastic moduli of 3DOWCs; (d–h) represent the strength properties of 3DOWCs.

Figure 16. Reduction and reduction efficiency of mechanical properties for mesoscale models. (a–c) represent the modulus reduction and reduction efficiency; (d–h) represent the strength reduction and reduction efficiency.

Table 1

Mechanical properties of fiber and matrix.

	${\mathit{E}}_{\mathbf{1}}/\mathbf{GPa}$	${\mathit{E}}_{\mathbf{2}}/\mathbf{GPa}$	${\mathit{\upsilon }}_{\mathbf{12}}$	${\mathit{G}}_{\mathbf{12}}/\mathbf{GPa}$	${\mathit{G}}_{\mathbf{23}}/\mathbf{GPa}$	${\mathit{\sigma }}_{\mathit{x}\mathit{t}}/\mathbf{MPa}$	${\mathit{\sigma }}_{\mathit{x}\mathit{c}}/\mathbf{MPa}$
Fiber	231	12.97	0.3	11.28	4.45	3350	2500
Matrix	3.5		0.35			103.4	241