Multiscale model construction and modulus prediction of needle-punched carbon/carbon composites

Abstract

In this paper, the modulus prediction of NP-C/C is completed by establishing a multiscale finite element model, and the influence weight analysis of the modulus is carried out by using the analytical homogenization method. Firstly, three Micro-RVE, including unidirectional composites, short fiber reinforced composites and porous matrix, are established to characterize the layers of NP-C/C. Then, the finite element model is used for homogenization and passed upward into Meso-RVE, and the modulus prediction is completed. Finally, the analytical homogenization method is used combined with the MATLAB^® script to analyze the influence weight of the modulus by adjusting the parameters. The results show that the average errors of stiffness coefficient prediction values obtained based on the analytical homogenization method are only 6.0% relative to the experimental values, which are lower than 25.3% of the finite element model. The modulus prediction method given in this study does not need to establish a micromechanical model with complex structure, and can be directly completed by computer script when the prediction accuracy is met, which is simple and feasible.

Keywords

Needle-punched carbon/carbon composite representative volume element modulus prediction analytical homogenization method multiscale finite element model

Introduction

Carbon/Carbon (C/C) composite are considered to be one of the most promising high temperature materials¹ due to their high specific strength and specific modulus, excellent ablation resistance and thermal insulation properties.² Needle-punched carbon fiber preforms are made by laying short-chopped fiber felt between layers of unidirectional fiber bundles and then densely needling;³ the purpose of the needle-punched is to hook the short fibers in the fiber felt into the thickness direction,⁴ thereby increasing the strength of the preform thickness direction. Chemical Vapor Infiltration (CVI) is used for the densification of preforms, which refers to the placement of needle-punched carbon fiber preforms in hydrocarbon gases and vapor deposition until the inter-fiber pores are filled;⁵ Liquid-phase impregnation refers to the preform impregnated in the resin and heated to generate a pyrolytic carbon matrix^;6 The two processes are combined to form Needle-Punched Carbon/Carbon (NP-C/C) composite with some strength in the thickness direction. The needling process improves the interlaminar strength of 2D composites, and 2.5D NP-C/C is simpler to process, cheaper to produce, and easier to batch produce than 3D woven composites.⁷

Needle-punched Carbon/Carbon composites require multiscale analysis because of their properties that combine fiber arrangement, needling process, matrix, interface and pore distribution. In terms of experimental observations: Roy et al.⁸ observed the effects of needling density and needling depth on microscale fiber orientation, and analyzed the effect of needling process on tensile strength. Drach et al.⁹ determined the shape distribution and porosity of pores in the pyrolysis carbon matrix according to the specimens made by CVI, and gave the influence of pores on elastic properties. Bradley et al.¹⁰ provided shear failure experimental measurements in three ply forms, with shear failure stresses ranging from 30–45 MPa in-plane and 8–22 MPa out-of-plane, and shear modulus ranging from 6–16 GPa in-plane and 1–3 GPa out-of-plane. Based on the results of X-ray scan, Piat et al.¹¹ simplified the layered form and material properties, established a micromechanical model to describe the elastic behavior, and determined the elastic modulus of each layer and the stiffness coefficient ofNP-C/C. Cao et al.¹² studied the failure mode of the composite, and observed that the interface disbonding under shear load was one of the main failure modes. Microscopic toughening was also analyzed,¹³ and after using carbon nanotubes to toughen the interface, the shear strength measured by nanoindentation characterization experiments increased by 80%; Macro-Sample have also been tested,¹⁴ and the failure mode of sandwich composite blades depends on core stiffness and surface texture.

In terms of Finite Element (FE) model: Yu et al.¹⁵ based on the discrete grayscale image of Micro-CT, established a three-dimensional model that can accurately reflect the needling effects and microstructure, which can better reveal the failure mechanism under tensile and compressive loads; however, shape reconstruction is only for specific specimens, and the generalizability of prediction method is poor. Xie et al.¹⁶ and¹⁷ selected and established four typical RVEs for predicting the effective performance of single needling regions and overlapping needling region by observing the microstructure of NP-C/C, taking into account local geometric features such as fiber fracture and deflection. Due to the increased complexity of the finite element model due to the high needling density and the presence of overlapping needled fibers, it is not possible to further analyze the influence of geometric parameters on macroscopic properties. Xu et al.¹⁸ and Zhang et al.¹⁹ used RVE to establish a micromechanical layered model to express the elastic properties of NP-C/C. In order to accurately calculate the material properties, the model needs to express the combined action results of fiber, matrix and porosity at the same time, but the calculation process is more complicated. Based on Micro-CT scan results, Meng et al.²⁰ established Meso-RVE and predicted the elastic modulus of NP-C/C by hierarchical modeling that took into account material changes caused by needle-punch, including fiber dislocation, fiber distribution, and matrix effects.

Through the above literature review, it is found that the modulus prediction of NP-C/C can be accurately completed by establishing multiscale models combined with finite element analysis. However, due to the complexity of the microstructure, it is necessary to repeatedly check the finite element model to improve the prediction accuracy. In addition, analytical homogenization methods²¹ are used to predict modulus with the help of a series of numerical formula, including the Mori-Tanaka method,²² the Halpin-Tsai formula,²³ the Chamis empirical formula,²⁴ and the classical laminate theory.²⁵ The analytical homogenization method used in this paper does not require complex meshing and impose the periodic boundary conditions, which can reduce the complexity of model construction, and can provides an effective method for predicting the modulus and improves the applicability of the method. At the same time, while using this method for predicting modulus, we further evaluate the influence of each parameter on the modulus in combination with the MATLAB^® script, which can be further used for model structure optimization.

Model construction of needle-punched carbon/carbon composites

Construction of representative volume elements

Referring to the Micro-CT scan of Meng et al.,²⁰ it can be seen that NP-C/C has obvious layered structure and periodic reproducibility, which can be divided into three sub-regions: unidirectional fiber bundle (or called non-woven cloth), short-chopped fiber felt and needling fiber. Unidirectional fiber bundles are regarded as Unidirectional (UD) composite, short-chopped fiber felt is regarded as Short Fiber Reinforced Composite (SFRP), the properties of needling fibers are simplified to be same as unidirectional fiber bundles.

The multiscale finite element model established in this paper includes three cubic Micro-RVE and one laminate Meso-RVE. When the edge length of the micromechanical model is in the range of 100–500 μm, the influence of porosity on the matrix can be effectively analyzed;⁹ Based on this, a Micro-RVE with a side length of 100 μm is established to describe the porous matrix. The preliminary analysis of the RVE size by Melro et al.²⁶ proposed that the results of the finite element model tend to converge when the RVE side length is greater than 10 times the fiber radius. Since the fiber diameter is considered to be a constant 7 μm, a Micro-RVE with a side length of 35 μm is established to describe UD Composite, and a Micro-RVE with a side length of 100 μm is established to describe SFRP. When the thickness of NP-C/C is less than 10 mm, the needle puncture characteristics of each layer are basically uniform and similar;²⁷ from this, Meso-RVE containing a needling fiber is established, the side length of which is defined as 9.08 mm,¹⁸ which was sufficient to capture the characteristics produced by the needle and could reduce the computational cost. Using the sequential multiscale method, the mechanical properties of Mico-RVE are transferred upward to Meso-RVE and to Macro-Sample, as shown in Figure 1.

Figure 1.

Multiscale model construction and analytical homogenization methods of NP-C/C (Reprinted from [20]. with copyright permission from Elsevier).

This section establishes a multiscale model for finite element and analytical homogenization; The volume fraction and material properties of each phase of the NP-C/C composite are taken from Ref.²⁰, and the lower modulus of the pores is given according to the Ref.¹⁸, and the material properties of each phase of the material are shown in Table 1.

Table 1.

The Material properties of each phase of the NP-C/C composite.

Material	Parameter	Value
Fiber²⁰	Longitudinal modulus E_f11 (GPa)	200
	Transverse modulus E_f22,E_f33 (GPa)	30
	Longitudinal shear modulus G_f12,G_f13(GPa)	9.0
	Transverse shear modulus G₃₁ (GPa)	11.881
	Longitudinal Poisson’s ratio v_f12, v_f13	0.20
	Transverse Poisson’s ratio v_f23	0.27
	Volume fraction of fiber bundles V_cf	0.2531
	Volume fraction of fiber felt V_sf	0.1097
	Volume fraction of needling fiber V_nf	0.0108
	Fiber diameter d_f (μm)	7
	Fiber length l_f (μm)	800
Matrix²⁰	Elastic modulus E_m (GPa)	15.9
	Shear modulus G_m (GPa)	5.7
	Poisson’s ratio v_m	0.395
Pore¹⁸	Elastic modulus E_p (GPa)	2.5 × 10⁻⁵
	Shear modulus G_p (GPa)	1 × 10⁻⁵
	Poisson’s ratio v_p	0.25
	Volume fraction of fiber bundles V_cp	0.09
	Volume fraction of fiber felt V_sp	0.20
	Volume fraction of needling fiber V_np	0.00

Micro-RVE homogenization of porous matrix

The porous matrix can be regarded as ellipsoid inclusion problem after simplification,⁹ where the pores are arbitrarily oriented in spatial position and angle, so the equivalent matrix can be considered as an isotropic material. This section also calculates the modulus of equivalent matrix referring to the Mori-Tanaka method²² shown by formula (1).

{\begin{matrix} \begin{array}{c} G_{r} = \frac{3 (2 G_{m} + K_{m})}{10 G_{m} (4 G_{m} + 3 K_{m})}, & K_{r} = \frac{1}{3 (4 G_{m} + 3 K_{m})} \end{array} \\ \begin{array}{c} G_{e m} = G_{m} + \frac{G_{p} - G_{m}}{1 + 4 V_{m} G_{r} (G_{p} - G_{m})} V_{p}, & K_{e m} = K_{m} + \frac{K_{p} - K_{m}}{1 + 9 V_{m} K_{r} (K_{p} - K_{m})} V_{p} \end{array} \end{matrix}

(1)

where the subscripts m, p, and em represent the matrix, pore, and equivalent matrix, respectively; G and K are the shear modulus and the bulk modulus; V denotes the volume fraction.

Due to the different pore volume fractions in the fiber felt and fiber bundle, the finite element model and analytical homogenization method (abbreviated as analytical method) are used to obtain the modulus prediction value of the equivalent matrix. It can be seen from Table 2 that when there is a small amount of pore inclusions (V_p = 9%), the relative error of the modulus prediction value under the two methods is 1.4%, which is well fitted. The errors for medium inclusions (V_p = 20%) are 10.2% and 9.6%, both within the acceptable range. The results show that the Mori-Tanaka method can be used for modulus prediction, and the relative error obtained under this method is positively correlated with the volume fraction of the inclusion.

Table 2.

Predicted moduli of equivalent matrix by porous matrix.

Equivalent matrix source	Pore volume fraction, %	Elastic moduls	FE model	Analytical method	Relative error, %
Fiber bundles	9	E(GPa)	14.04	13.84	1.4
Fiber bundles	9	G(GPa)	5.03	4.96	1.4
Fiber felt	20	E(GPa)	12.37	11.22	10.2
Fiber felt	20	G(GPa)	4.43	4.04	9.6

Micro-RVE homogenization of unidirectional fiber bundles

The homogenization of the fiber bundle lies in its equivalence with the transverse isotropic material. According to the stiffness matrix form, the periodic strain is first applied to the Micro-RVE to calculate its stress situation. Then, the stress-strain equation is solved, the equivalent stiffness coefficient is obtained, and the stiffness matrix calculation of the homogenized layer is completed. Finally, the stress and strain differential elements of a single fiber or matrix can be derived by formula (2):

\begin{array}{l} {\bar{σ}}_{i j} = E_{i j} {\bar{ε}}_{i j} \\ {\bar{σ}}_{i j} = \frac{1}{V} \int_{Ω} σ_{i j} (x, y, z) d v = \frac{1}{V} \sum_{m = 1}^{N_{m}} σ_{i j}^{m} \\ {\bar{ε}}_{i j} = \frac{1}{V} \int_{Ω} ε_{i j} (x, y, z) d v = \frac{1}{V} \sum_{m = 1}^{N_{m}} ε_{i j}^{m} \end{array}

(2)

where E_ij is the stiffness matrix;

{\bar{σ}}_{i j}

is the average stress tensor;

{\bar{ε}}_{i j}

is the average strain tensor; m represents the fiber or matrix number, N_m represents the total number of fibers,

σ_{i j}^{m}

and

ε_{i j}^{m}

represents the fiber stress and strain numbered m, V represents the total volume of the homogenization layer.

Continue to homogenize the fiber bundles using the Chamis empirical formula,²⁴ shown as formula (3). The properties of UD composites are affected by fiber shape, arrangement and load form, that is, there is a Hashin-Shtrikman boundary²¹ in the elastic modulus; the principle of minimum potential energy gives the upper limit of Voigt and the principle of minimum residual energy gives the lower limit of Reuss. The accuracy of predicting the longitudinal modulus (E₁₁ and G₁₂) of UD composites by using the Chamis empirical formula is high, but the prediction accuracy of the transverse modulus (E₃₃ and G₂₃) needs to be further judged according to Micro-RVE.

{\begin{matrix} \begin{array}{c} E_{11} = V_{f} E_{f 11} + (1 - V_{f}) E_{m}, & E_{22} {= E}_{33} = \frac{E_{m}}{1 - \sqrt{V_{f}} (1 - E_{m} / E_{f 22})} \end{array} \\ \begin{array}{c} G_{12} {= G}_{13} = \frac{G_{m}}{1 - \sqrt{V_{f}} (1 - G_{m} / G_{f 12})}, & G_{23} = \frac{G_{m}}{1 - \sqrt{V_{f}} (1 - G_{m} / G_{f 23})} \end{array} \\ \begin{array}{c} v_{12} = v_{13} = V_{f} v_{f 12} + (1 - V_{f}) v_{m 12}, & v_{23} = E_{22} / (2 G_{23}) - 1 \end{array} \end{matrix}

(3)

where the subscripts f and m represent the fiber phase and the matrix phase, respectively; E, G, and v represent the Young’s modulus (GPa), shear modulus (GPa), and Poisson’s ratio, respectively; V_f is the fiber volume fraction.

Referring to random UD composites,²⁸ regular UD composites is established for comparison. The finite element models under different fiber arrangement and load form are established, and the internal stress distribution under lateral load of UD composites can be obtained, as shown in Figure 2.

Figure 2.

Stress field diagram of UD composites: (a) Regular;(b)Random; (c)-(d) Transverse tension; (e)-(f) Transverse shear.

Figure 2(c) and (e) present that the stress distribution of regular UD composites is uniform and periodically replicable. Figure 2(d) and (f) show that random UD composites have stress concentrations between fibers, with maximum stress values increasing by 14% and 20% under transverse tension and shear loads. The modulus deviation error of the two UD composites is calculated at the same time (see Table 3), and the results show that the maximum deviation error is generated by the transverse shear modulus G₂₃ (1.6%), and the rest of the deviation error is within 1.4%. Therefore, the UD composites analyzed in this section can ignore the modulus deviation error caused by the fiber arrangement and load form. Next, the modulus of the fiber bundle is obtained using the Chamis empirical formula and compared with the literature value,²⁰ shown as Table 3. It shows that the average relative error of the predicted value of Young’s modulus is 8.6%, shear modulus is 3.0%. It can be seen that the analysis homogenization method can also be used for approximate predictive analysis of modulus under the condition of satisfying the prediction accuracy.

Table 3.

Prediced moduli of fiber bundles with literature values.

Elastic modulus	Finite elements model			Literature value²⁰	Analytical method	Relative error
Elastic modulus	Random	Regular	Error	Literature value²⁰	Analytical method	Relative error
E₁₁(GPa)	59.524	60.373	1.4%	58.589	61.065	4%
E₂₂(GPa)	18.910	18.708	−1.1%	20.975	19.119	−9%
E₃₃(GPa)	18.910	18.708	−1.1%	22.035	19.119	−13%
G₁₂(GPa)	5.795	5.809	0.2%	6.267	6.450	3%
G₂₃(GPa)	6.033	6.134	1.6%	6.757	7.068	5%
G₃₁(GPa)	5.795	5.809	0.2%	6.401	6.450	1%
v ₁₂	0.349	0.347	−0.6%	0.3466	0.346	0%
v ₂₃	0.515	0.521	1.2%	0.4858	0.453	−9%
v ₃₁	0.349	0.347	−0.6%	0.3199	0.346	8%

Micro-RVE homogenization of short-chopped fiber felt

In this section, Halpin-Tsai formula²³ is used to complete the modulus prediction of fiber felt. When giving the fiber length, slenderness ratio and volume fraction, the modulus prediction of unidirectional short fiber composite can be derived using the following:

\begin{array}{c} \frac{E_{1}}{E_{m}} = \frac{1 + (2 l_{f} / d_{f}) η_{1} V_{f}}{1 - η_{1} V_{f}}, & η_{1} = \frac{E_{f} / E_{m} - 1}{E_{f} / E_{m} + 2 (l_{f} / d_{f})} \\ \frac{E_{2}}{E_{m}} = \frac{1 + 2 η_{2} V_{f}}{1 - η_{2} V_{f}}, & η_{2} = \frac{E_{f} / E_{m} - 1}{E_{f} / E_{m} + 2} \end{array}}

(4)

where E₁ and E₂ are the transverse and longitudinal modulus, l_f and d_f are the length and diameter of the unidirectional short fibers.

For fiber felt composed of random short fibers (Considered SFRP), Young’s modulus E_s and shear modulus G_s are calculated as formula (5).

E_{s} = \frac{3}{8} E_{1} + \frac{5}{8} E_{2}, G_{s} = \frac{1}{8} E_{1} + \frac{1}{4} E_{2}

(5)

This section analyzes the accuracy of modulus predictions based on the literature results from Meng et al.²⁰ (based on finite element model) with the current Halpin-Tai formula predictions, see Table 4. Fiber felt in literature used the sequential random adsorption method²⁹ generation, which allows the given directional angle range of short fibers, and the model can be regarded as a transverse isotropic material after homogenization; However, the random growth method³⁰ is used in this paper, and the spatial position and orientation angle of the short fibers are completely random, so the transverse shear modulus (G₂₃ and G₃₁) calculated by the literature values has a certain error with the analytical method .

Table 4.

Prediced moduli of short-chopped fiber felt with literature values.

Elastic modulus	Literature value²⁰	Analytical method	Relative error
E₁₁(GPa)	20.971	20.468	−2%
E₂₂(GPa)	20.971	20.468	−2%
E₃₃(GPa)	19.864	20.468	3%
G₁₂(GPa)	7.533	7.430	−1%
G₂₃(GPa)	6.022	7.430	23%
G₃₁(GPa)	6.022	7.430	23%
v ₁₂	0.3920	0.3774	−4%
v ₂₃	0.3624	0.3774	4%
v ₃₁	0.3625	0.3774	4%

The fiber felt is equivalent to the isotropic material to obtain the predicted elastic modulus. Except for the large relative error of the in-plane shear modulus (caused by the different forms of short fiber generation), the relative error of the predicted value of the modulus of fiber felt is within 4%. Since fiber felt is not used as a reinforced phase, the distribution mode of short fiber has little influence on the overall modulus of NP-C/C. The result of Table 4 also shows that, the proposed analytical homogenization method can be used to predict the modulus of fiber felt.

Meso-RVE modulus calculation

Modulus calculation based on finite element model

Establishing Meso-RVE with symmetry of structure and load reduces model size and analysis time. However, Meso-RVE is a layered structure, the material properties are not continuous, and the problems of stress discontinuity and strain incoordination are prone to occur during loading, which will lead to deviation when solving the stiffness matrix. To ensure the coordination of deformation under load, we define and apply the periodic boundary conditions.¹⁹

{\begin{array}{c} \vec{u} (0, x_{2}, x_{3}) - \vec{u} (L_{1}, x_{2}, x_{3}) = \vec{U_{1}} \\ \vec{u} (x_{1}, 0, x_{3}) - \vec{u} (x_{1}, L_{2}, x_{3}) = \vec{U_{2}} \\ \vec{u} (x_{1}, x_{2}, 0) - \vec{u} (x_{1}, x_{2}, L_{3}) = \vec{U_{3}} \end{array}

(6)

where L_i(i = 1,2,3) is the length of the RVE, U_i is the relative displacement of the two nodes on the opposite surface.

First, periodic strain is applied to the finite element model, equating it with a homogeneous elastomer. Meso-RVE can be regarded as an orthotropic material after homogenization, thus we can derive the stiffness matrix by formula (7).

[\begin{array}{c} {\bar{σ}}_{x} \\ {\bar{σ}}_{y} \\ {\bar{σ}}_{z} \\ {\bar{τ}}_{x y} \\ {\bar{τ}}_{x z} \\ {\bar{τ}}_{y z} \end{array}] = [\begin{array}{c} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{13} & C_{23} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{array}] [\begin{array}{c} {\bar{ε}}_{x} \\ {\bar{ε}}_{y} \\ {\bar{ε}}_{z} \\ {\bar{ε}}_{x y} \\ {\bar{ε}}_{x z} \\ {\bar{ε}}_{y z} \end{array}]

(7)

Secondly, the equivalent stiffness coefficient in the loading direction is calculated according to formula (8), and the loading is carried out multiple times according to the number of independent stiffness coefficients until the stiffness matrix is solved. Loads several times according to the strain number to solve all the stiffness coefficients.

E_{V}^{1} = \frac{1}{2} {(\bar{σ^{1}})}^{T} \bar{ε^{1}} V = \frac{1}{2} {(\bar{ε^{1}})}^{T} [C] \bar{ε^{1}} V = \frac{1}{2} C_{11} V \Rightarrow C_{11} = \frac{2 E_{V}^{1}}{V}

(8)

Then, multiple sets of periodic strains are applied to the finite element model according to Table 5, and the deformation is amplified according to the deformation scaling factor of 10, and the resulting deformation and stress field are analyzed.

Table 5.

Strain case and stiffness coefficient calculations.

Load case	Periodic strain	Stiffness coefficient
1	ε¹ = [0.1 0 0 0 0 0]^T	$C_{11} = 2 E_{V}^{1} / V$
2	ε² = [0 0.1 0 0 0 0]^T	$C_{22} = 2 E_{V}^{2} / V$
3	ε³ = [0 0.1 0 0 0 0]^T	$C_{33} = 2 E_{V}^{3} / V$
4	ε⁴ = [0 0 0 0.1 0 0]^T	$C_{44} = 2 E_{V}^{4} / V$
5	ε⁵ = [0 0 0 0 0.1 0]^T	$C_{55} = 2 E_{V}^{5} / V$
6	ε⁶ = [0 0 0 0 0 0.1]^T	$C_{66} = 2 E_{V}^{6} / V$
7	ε⁷ = [0.1 0.1 0 0 0 0]^T	$C_{12} = C_{21} = (E_{V}^{7} - E_{V}^{1} - E_{V}^{2}) / V$
8	ε⁸ = [0.1 0 0.1 0 0 0]^T	$C_{13} = C_{31} = (E_{V}^{8} - E_{V}^{1} - E_{V}^{3}) / V$
9	ε⁹ = [0 0.1 0.1 0 0 0]^T	$C_{23} = C_{32} = (E_{V}^{9} - E_{V}^{1} - E_{V}^{3}) / V$

According to Figure 3, it can be observed that the model has obvious stress stratification under the in-plane tensile strain(ε_x and ε_y), and the fiber bundle in the same direction as the load is used as the reinforced phase. Under shear strain, the stress of each layer is relatively continuous, and there is no obvious reinforced phase.

Figure 3.

Deformation field and stress field diagram of FE model.

All stiffness coefficients can be calculated according to multiple sets of periodic strains, and the inverse is obtained to obtain the flexibility matrix [S], seen as formula (9). The modulus of NP-C/C is solved based on the finite element model, and the analysis results are shown in Table 6 below.

[S] = {[C]}^{- 1} = [\begin{array}{c} 1 / E_{11} & - v_{21} / E_{22} & - v_{31} / E_{33} & 0 & 0 & 0 \\ - v_{12} / E_{11} & 1 / E_{22} & - v_{32} / E_{33} & 0 & 0 & 0 \\ - v_{13} / E_{11} & - v_{23} / E_{22} & 1 / E_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 / G_{23} & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 / G_{31} & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 / G_{12} \end{array}]

(9)

Table 6.

Comparison of measured and predicted values of composite specimens.

Number	Stiffness factor (N/mm)	C ₁₁	C ₂₂	C ₃₃	Average error
Specimen 1	Experimental value¹¹	38.94	40.58	13.25
	Finite element model	32.33	32.33	18.44
	Relative error	−17%	−20%	39%	25.3%
	Analytical method	41.65	41.65	14.25
	Relative error	7%	3%	8%	6.0%
Specimen 2	Experimental value¹¹	51.12	8.94	15.21
	Finite element model	45.32	13.13	18.12
	Relative error	−11%	47%	19%	25.6%
	Analytical method	55.32	7.06	14.79
	Relative error	8%	−10%	−3%	7.0%

Finally, we can draw some important conclusions about the internal stress distribution of the FE model. Taking the x-direction tensile strain (ε_x = 0.1) as an example, there is an obvious stress concentration at needling fiber, and the xy plane stress field diagram shows a symmetrical distribution at the center, shown in Figure 4(a). The 0° fiber bundle with the same load is used as the reinforced phase, and the 90° fiber bundle carries most of the remaining stresses, and the average stress value σ₁₁ in the fiber direction is three times the average stress value σ₂₂ in the matrix direction (σ₁₁/σ₂₂ = 3.17), as shown in Figure 4(b).

Figure 4.

Stress field diagram of FE model: (a) In plane direction; (b) Out of face direction.

Figure 5.

Stress field diagram of interface: (a) Numbering form; (b) Interface-1; (c) Interface-2; (d) Interface-3.

The interface is alternately distributed between the fiber felt and fiber bundle, with three layers and numbered sequentially. The interface is mainly used to describe damage at the contact between two-phase materials, and stiffness values need to be given to ensure deformation coordination and prevent material embedding. This section focuses on predicting the effective modulus of elasticity of NP-C/C, without considering interlaminar damage, so only the interfacial stiffness is assigned to analyze the stress field. According to the research of Vaughan et al.,³¹ the displacement continuity can be maintained when the interface stiffness is greater than 10⁵ N/mm². Continue to analyze the interfacial stress of each layer using the x-direction tensile strain as an example (ε_x = 0.1).

Figure 5 indicates that the maximum stress values of the three interfaces are: Interface 1 is 5.3 MPa, Interface 2 is 1.4 MPa, while Interface 3 is 1.5 MPa; since the 0° fiber bundle is used as the reinforced phase, interface 1 carries the highest stress value. The interfacial stress field diagram is still symmetrically distributed along the center, and the stress concentration area is basically consistent with the direction of the continuous fibers in the fiber bundle. The stress concentration at the interface occurs at the contact between the needling fiber and the fiber bundle, which is consistent with the work of Nie et al.,³² who found that the needle-punched damages caused the failure of the composite.

Modulus calculation based on analytical homogenization method

This section continues to complete the calculation of the modulus based on the analytical homogenization method. Based on the above analysis, three Micro-RVE are established to characterize unidirectional fiber bundles, short-chopped fiber felts, and needling fibers in NP-C/C. Using Meso-RVE as the smallest hierarchical unit of NP-C/C, the layered form shown in Figure 6 is obtained. At the same time, according to the direction and distribution of fiber phase, matrix phase and pore phase, the respective homogenization process is completed.

Figure 6.

Meso-RVE model: (a) Layered form; (b) Geometric parameters.

where l_n is the side length of Meso-RVE, t_s is the thickness of the short-chopped fiber felt, t_u is the thickness of the unidirectional fiber bundles, and ϕ_n is the diameter of the needling fiber.

Since the thickness of each layer of NP-C/C is constant, Meso-RVE without needling fibers can be considered laminates. According to the classical laminate theory,²⁵ the calculation formula for the total stiffness matrix [C]_u of fiber bundles under different laying sequences is:

{[C]}_{u} = \sum_{i}^{n} R_{z} {(θ_{i})}^{T} \times {[C]}_{u}^{i} \times R_{z} (θ_{i})

(10)

where the subscript i represents the i-th layer of fiber bundles, θ_i represents layup angle, R_z(θ_i) represents the rotation matrix (rotated by θ_i along the z-axis), and n is the total number of layups.

The stiffness matrix of fiber bundle and fiber felt is linearly superimposed according to the layer thickness ratio, and the laminates stiffness matrix [C]_L (where [C]_L = [S]⁻¹_L) is calculated. Because the needling fibers follow the z direction, the laminate flexibility matrix is used to stack linearly with the needling fiber flexibility matrix based on the area ratio. The Meso-RVE flexibility matrix [S]_R is calculated, shown as in the Formula (11). The modulus of NP-C/C can be calculated by the classical laminate theory, and the analysis results are shown in Table 6.

{\begin{array}{c} {[C]}_{L} = \frac{t_{u}}{t_{u} + t_{s}} {[C]}_{u} + \frac{t_{s}}{t_{u} + t_{s}} {[C]}_{s} \\ {[S]}_{R} = (1 - \frac{π ϕ_{n}^{2} / 4}{l_{n}^{2}}) {[S]}_{L} + \frac{π ϕ_{n}^{2} / 4}{l_{n}^{2}} {[S]}_{n} \end{array}

(11)

where the subscripts n, u, s, L, and R represent needling fiber, unidirectional fiber bundles, short-chopped fiber felt, Laminate and Meso-RVE, respectively.

Modulus comparison and analysis of needle-punched carbon/carbon composites

Modulus comparison based on analytical method and finite element model

As mentioned earlier, a multiscale finite element model is established and the modulus of NP-C/C is predicted using the analytical homogenization method. According to the distribution mode of fiber, matrix and pores, Micro-RVE and Meso-RVE were selected to establish a multiscale model; and through a series of numerical formula, we complete the analytical homogenization of the micromechanical model.

In this section, using the finite element model described above to compare with the experimental results of the Ref.¹¹, the prediction ability of the analytical homogenization method can be assessed and verified. In order to increase the test data, the uniaxial tensile test values of two specimens were selected for analysis; Specimen 1 adopts traditional orthogonal laying, that is, the fiber bundle is laid along 0° and 90°; All fiber bundles are in the same direction in Specimen 2, and the rest of the material parameters are the same. Table 7 shows the setting of the NP-C/C geometric parameters, and the corresponding stiffness coefficient predictions and measurements are shown in Table 6.

Table 7.

Geometric parameters of needle-punched carbon/carbon composites.

	Unidirectional fiber bundles	Short-chopped fiber felt	Needling fiber	Meso-RVE
Geometric parameters¹⁸	t_u = 0.94 mm	t_s = 0.42 mm	ϕ_n = 1.45 mm	l_n = 9.07 mm

Table 7 indicates that the average errors of the predicted values of stiffness coefficients of the two specimens obtained based on the analytical method are 6.0% and 7.0%, which are lower than the average relative errors of 25.3% and 25.6% of finite element model. The results show that the predicted values of the three stiffness coefficients of the two specimens under the analytical method are in good agreement with the experimental values.

The weight analysis of the influence of each parameter on the modulus

In the previous section, the prediction of NP-C/C modulus was realized based on the analysis homogenization method. In this section, the influence weight of each parameter on the modulus was analyzed by adjusting the structural parameters with the help of MATLAB script. Generally speaking, the parameters that can be adjusted in the manufacturing process are: thickness ratio of fiber felt and fiber bundles (t_s/t_u), fiber volume content (including V_cf and V_sf), and needling density ρ_n. In this section, we focus on analyzing the Young’s modulus (E₁₁ and E₃₃) and shear modulus (G₁₂ and G₂₃) of NP-C/C, by characterizing the influence weight of each parameter on the modulus (the total influence is 100%), to provide a reference for the optimization of mechanical properties.

Relying on the Monte Carlo simulation³³ shown in Figure 7, the geometric parameters are used as input variables and entered into the established analytical homogenization method to obtain the output response (modulus of NP-C/C). According to the Ref.³⁴, the initial input variable is mutated in a proportion of 10% and re-input to obtain a new output response value; then normalizing the input variables to [-1,1], the model is fitted by the ordinary least squares estimation to obtain the regression coefficient a_ij, which can be derived using the following estimation formula (12).

Y_{j} = a_{0}^{j} + \sum_{i = 0}^{n_{v}} a_{i}^{j} X_{i}, j = 1,2,3, \dots n_{r}

(12)

where a₀^j is the initial elastic constant, X_i is the input variable; Y_j is the output response; a_i^j is the linear regression coefficient; n_v, n_r represent the number of input variables and output responses, respectively.

N_{k} = a_{k}^{j} / \sum_{i = 0}^{n_{v}} | a_{i}^{j} |, k = 1,2,3, \dots n_{r}

(13)

Figure 7.

Schematic diagram of Monte Carlo simulation.

Finally, according to formula (13), each linear regression coefficient a_i^j is converted into the contribution rate N_k of the input variable to the output response, which characterizes the sensitivity of each parameter to the structural response.

Take the traditional [0/w/90/w]_n layered specimen 1 as an example, and its structural parameters are shown in Table 7. The initial structural parameters (including: t_s/t_u, V_cf, V_sf, and ρ_n) were each subjected to 10% mutation using MATLAB, the contribution rate of each parameter to the NP-C/C modulus is shown in Figure 8.

Figure 8.

Analysis of the contribution rate of each parameter to the modulus.

The results show that, increasing the needling density ρ_n favors the out-of-plane properties (E₃₃ and G₂₃) of NP-C/C, but reduces the in-plane properties (E₁₁ and G₁₂). Fiber volume fractions (including V_cf and V_sf) has the greatest influence weight on the modulus, but the impact is different: V_cf mainly affects the in-plane properties, and V_sf mainly affects the out-of-plane properties. Increasing the fiber volume fraction can significantly improve the modulus in all directions. At the same time, it can be found that the increase of the layer thickness ratio (t_s/t_u) will reduce the in-plane modulus E₁₁, but the influence weight on the remaining modulus is less than 10%.

Conclusion

In order to better predict the modulus of Needle-punched Carbon/Carbon composites, this paper complete modulus prediction based on the multiscale finite element model, and carefully analyze the influence of each parameter on the modulus by employing the analytical homogenization method and the MATLAB script. The research shows that, the analysis homogenization method used in this paper can directly program all numerical equations into a single computer script without modeling, and we can realize geometric parameter adjustment through Monte Carlo simulation, which simplifies the complexity of model construction and even enables structural optimization. The relevant conclusions of the study are as follows:

(1) The modulus prediction value of NP-C/C based on the analytical homogenization method is in good agreement with the experimental value; the average errors of the predicted values of stiffness coefficient and the experimental values are 6.0%, which are lower than 25.3% of the finite element model.

(2) Increasing the needling density ρn is beneficial to improve the (positively correlated) out-of-plane properties of NP-C/C (E₃₃ and G₂₃), but will reduce the (negatively correlated) in-plane properties of NP-C/C (E₁₁ and G₁₂).

(3) Fiber volume fraction (including V_cf and V_sf) has the greatest influence weight on NP-C/C modulus, but V_cf mainly affects the in-plane properties, and V_sf mainly affects the out-of-plane properties. Increasing the fiber volume fraction can significantly improve the modulus of NP-C/C in all directions.

Footnotes

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) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Yang Liu

References

Thakre

Soni

. Modeling of burr size in drilling of aluminum silicon carbide composites using response surface methodology. Engineering Science and Technology, an International Journal 2016; 193: 1199–1205.

Baghery Borooj

Mousavi Shoushtari

Haji

, et al. Optimization of plasma treatment variables for the improvement of carbon fibres/epoxy composite performance by response surface methodology. Composites Science and Technology 2016; 128: 215–221.

Han

Zhou

Silberschmidt

. Mesoscale damage analysis of needle-punched carbon/carbon composite considering randomness of inherent defects. Composites Science and Technology 2019; 183: 107821.

Yuan

Zhao

. Mechanical properties of a novel Tri-directional carbon-flax-aramid fiber reinforced composite. Compos Sci Technol 2021; 213(2021): 108923.

Han

Shi

Han

, et al. Preparation and deposition mechanism of pyrolytic carbon by CVI using 3D Ni/wood-carbon catalyst. Carbon Lett 2022; 321: 143–151.

Wei

Guo

Cui

, et al. Microstructure evolution and performance improvement of silicon carbide ceramics via impregnation method. Mater 2022; 155: 1717.

Shadhin

Jayaraman

Rahman

. Effect of non‐woven flax mat manufacturing parameters and consolidation pressure on properties of composites manufactured using vacuum‐assisted resin transfer molding. Polym Compo 2021; 4211: 5743–5754.

Roy

Ishtiaque

. Influence of punching parameters on fibre orientation and related physical and mechanical properties of needle punched nonwoven. Fibers Polym 2019; 201: 191–198.

Drach

Tsukrov

Gross

, et al. Numerical modeling of carbon/carbon composites with nanotextured matrix and 3D pores of irregular shapes. Int J Solids Struct 2011; 4818: 2447–2457.

10.

Bradley

Bowen

McEnaney

, et al. Shear properties of a carbon/carbon composite with non-woven felt and continuous fibre reinforcement layers. Carbon 2007; 4511: 2178–2187.

11.

Piat

, et al. Modeling of effective elastic properties of carbon/carbon laminates Proceedings of seventeenth international conference on composite materials. Edinburgh, UK, 2009.

12.

Cao

Malakooti

Kulkarni

, et al. Nanoindentation measurement of core–skin interphase viscoelastic properties in a sandwich glass composite. Mech Time Depend Mater 2021; 253: 353–363.

13.

Wang

Andrade

MJd

, et al. The interfacial shear strength of carbon nanotube sheet modified carbon fiber composites. Challenges in Mechanics of Time Dependent Materials 2021; 2. Cham: Springer, 2021, pp. 25–32.

14.

Cao

Malakooti

Kulkarni

, et al. The effect of resin uptake on the flexural properties of compression molded sandwich composites. Wind Energy 2022; 251: 71–93.

15.

Zhou

Zhang

. A micro-image based reconstructed finite element model of needle-punched C/C composite. Compos Sci Technol 2017; 153: 48–61.

16.

Xie

Fang

Chen

, et al. Numerical and experimental studies on scattered mechanical properties for 3D needled C/C-SiC composites. Compos Struct 2018; 192: 545–554.

17.

Xie

Liang

Fang

, et al. Effect of needling parameters on the effective properties of 3D needled C/C-SiC composites. Compos Sci Technol 2015; 117: 69–77.

18.

Zhang

, et al. Hierarchically modeling the elastic properties of 2D needled carbon/carbon composites. Compos Struct 2015; 133: 148–156.

19.

Zhang

Wan

, et al. Impact compressive behavior and failure modes of four-step three-dimensional braided composites-based meso-structure model. Int J Damage Mech 2015; 246: 805–827.

20.

Meng

Song

, et al. Predicting the effective properties of 3D needled carbon/carbon composites by a hierarchical scheme with a fiber-based representative unit cell. Compos Struct 2017; 172: 198–209.

21.

Kamin´ ski

Pawlak

. Various approaches in probabilistic homogenization of the CFRP composites. Compos Struct 2015; 133: 425–437.

22.

Mori

Tanaka

. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica 1973; 215: 571–574.

23.

Affdl

JCH

Kardos

. The Halpin‐Tsai equations: a review. Polym Eng Sci 1976; 165: 344–352.

24.

Chamis

. Mechanics of composite materials: past, present and future. 21st annual meeting of the society for engineering science, 1984.No. E-3936.

25.

Öchsner

. Foundations of classical laminate theory. Cham: Springer, 2021.

26.

Melro

Camanho

Pinho

. Influence of geometrical parameters on the elastic response of unidirectional composite materials. Compos Struct 2012; 94(11): 3223–3231.

27.

Fang

Wang

, et al. An improved analytical method for calculating stiffness of 3D needled composites with different needle-punched processes. Compos Struct 2020; 237: 111938.

28.

El Moumen

Kanit

Imad

. Numerical evaluation of the representative volume element for random composites. European Journal of Mechanics - A/Solids 2021; 86: 104181.

29.

Pan

Iorga

Pelegri

. Numerical generation of a random chopped fiber composite RVE and its elastic properties. Compos Sci Technol 2008; 68: 2792–2798.

30.

Wang

Sha

, et al. Effect of randomly distributed fibers on residual stress of SiC/GH4738 composites. Composites Part C: Open Access 2021; 4(2021): 100101.

31.

Vaughan

McCarthy

. Micromechanical modelling of the transverse damage behaviour in fibre reinforced composites. Compos Sci Technol 2011; 71(3): 388–396.

32.

Nie

Zhang

, et al. Microstructure and tensile behavior of multiply needled C/SiC composite fabricated by chemical vapor infiltration. J Mater Processing Technol 2009; 2091: 572–576.

33.

Liu

. Optimal Monte Carlo method in estimating areas. Results in Applied Mathematics 2021; 12: 100205.

34.

Klopp

. A tutorial on testing the equality of standardized regression coefficients in structural equation models using Wald tests with lavaan. PsyArXiv 2019; 10.