Numerical Assessment of Effective Elastic Properties of Needled Carbon/Carbon Composites Based on a Multiscale Method

Abstract

Needled carbon/carbon composites contain complex microstructures such as irregular pores, anisotropic pyrolytic carbon, and interphases between fibers and pyrolytic carbon matrices. Additionally, these composites have hierarchical structures including weftless plies, short-cut fiber plies, and needled regions. To predict the effective elastic properties of needled carbon/carbon composites, this paper proposes a novel sequential multiscale method. At the microscale, representative volume element (RVE) models are established based on the microstructures of the weftless ply, short-cut fiber ply, and needled region, respectively. In the microscale RVE model, a modified Voronoi tessellation method is developed to characterize anisotropic pyrolytic carbon matrices. At the macroscale, an RVE model containing hierarchical structures is developed to predict the effective elastic properties of needled carbon/carbon composites. For the data interaction between scales, the homogenization results of microscale models are used as inputs for the macroscale model. By comparing these against the experimental results, the proposed multiscale model is validated. Furthermore, the effect of porosity on the effective elastic properties of needled carbon/carbon composites is investigated based on the multiscale model. The results show that the effective elastic properties of needled carbon/carbon composites decrease with the increase in porosity, but the extent of decrease is different in different directions.

1. Introduction

As a lightweight advanced material, carbon/carbon (C/C) composites have high strength-to-density and stiffness-to-density ratios, low thermal expansion, and excellent wear resistance [1,2,3]. During the fabrication of C/C composites, the chemical vapor infiltration method is usually used due to its low cost and the fact that it does little damage to the mechanical properties of carbon fibers [4,5]. In the aspect of the carbon fiber preform, the 2D carbon cloth, 3D braided structure, and needled structure are the three main types of preparing C/C composites. Among these preforms, 2D carbon cloth has low interlaminar shear strength and out-of-plane tensile strength and thus is applied in structures that sustain only an in-plane stress field [6]. The 3D braided preform can be used to produce complex near-net-shape composites, but it costs a lot [7]. Therefore, the needled preform becomes a moderate choice for preparing high-performance C/C composites. So far, needled C/C composites have been widely used as some structural components in the field of aerospace and aviation, such as rocket nozzles, throat liners, and aircraft brake disks [8,9,10]. Elastic properties are important criteria for evaluating the vibration of structural components and are also the foundation for analyzing the damage behaviors of needled C/C composites [7,11]. Thus, it is crucial to establish the relationship between the microstructures and elastic properties of needled C/C composites.

Currently, experimental, analytical, and numerical methods are commonly used to solve some issues about the relationship between microstructures and mechanical properties of composites. Among them, the experimental method is the most reliable and convincing, but it needs special equipment, a long period, and cannot be cost-effective [12]. Thus, the analytical and numerical methods become alternative and sought-after in many research fields. The analytical method mainly includes Eshelby’s theory [13,14,15], boundary theory [16,17], and empirical formula [18]. Although these theories have been developed for many years, they cannot continue to be used to obtain local stress/strain distributions in composites and are not applied in complex composites. Compared to the analytical method, numerical modeling is accepted to obtain more accurate results for the needled C/C composites with complex microstructures.

Efforts have been made in the numerical calculation of the effective elastic properties of needled C/C composites. Xu [6] and Tan [19] established a hierarchical multiscale model to predict effective elastic properties based on the microstructural characterization of needled C/C composites. Meng [10] considered the effects of needled fibers to develop a fiber-based representative volume element (RVE) model for estimating the elastic constants of needled C/C composites. Xie [20] constructed four typical RVE models to predict the effective elastic properties of needling and un-needling regions. Meanwhile, a periodic RVE model including needling and un-needling regions was built. However, these models are established based on the assumption of the isotropic pyrolytic carbon and neglecting the effect of the interphase between the fiber and matrix.

In this paper, a multiscale model is proposed with the consideration of microstructures of needled C/C composites. At the microscale, RVE models of the weftless ply, short-cut fiber ply, and needled fiber region are constructed by considering interphases and the morphology of the pyrolytic carbon, respectively, where the Voronoi tessellation method is modified to characterize the pyrolytic carbon morphology. At the macroscale, a macroscale RVE model is established to predict the effective elastic properties of needled C/C composites. The data interaction between the two scales is achieved by the sequential multiscale calculation method. Compared to the experimental results, the proposed multiscale model is validated. Based on this model, the effect of porosity on the effective elastic properties of needled C/C composites is further explored.

2. Experimental Procedure and Microstructural Analysis

The carbon fiber preforms provided by Jiangsu Tianniao High Technology Co., Ltd. (Yixing, China) have a needling density of 27 punch/cm² and a fiber volume fraction of 33% in this paper. The preform is alternately stacked with 0° weftless plies, short-cut fiber plies, and 90° weftless plies to a certain thickness, and then fixed by the needling technique in the stacked direction, as shown in Figure 1. During the needling process, short-cut carbon fibers of about 50~100 mm in length are introduced into the stacked direction by the needle plate to achieve the stitching of the preform.

In each ply, 12K PAN-based T700 carbon fibers are used, so the weftless, short-cut fiber ply and needled fiber region have the same fiber diameter. The scanning electron microscope (SEM) is used to obtain plenty of microstructural images of weftless plies, of which a representative image is given in Figure 2a. According to these SEM images, the ImageJ software is employed to obtain the statistical distribution of the fiber diameter, as shown in Figure 2b. The results show that the fiber diameter obeys the Gauss distribution, and the average fiber diameter is about 7 μm. Thus, this manuscript utilizes the average value to calculate the effective elastic properties of needled C/C composites.

In this study, needled carbon/carbon composites are prepared by the isothermal chemical vapor infiltration (ICVI) method [21,22]. When the gaseous precursor passes through the interior of the heated porous needled preform in the ICVI process, it will undergo the pyrolysis reaction to deposit pyrolytic carbon on the surface of the carbon fibers to realize the densification [4,23]. Figure 3 gives the schematic diagram of the ICVI device, where the red color indicates the heating elements, and the thermal couple is used to measure the temperature in the chamber. Before the carbon deposition, the fiber preform with the size of 100 mm × 100 mm × 10 mm is fixed using a porous graphite mold and placed in the center area of the chamber. A gas-tightness check is required to ensure an oxygen-free environment. During the densification, methane is used as the precursor, argon serves as the protective gas, the infiltration temperature is 800 °C~1200 °C, the infiltration pressure is 8~12kPa, and the infiltration time is 80~120 h.

Figure 4 gives some details about the microstructures of needled carbon/carbon composites. It is found that needled fibers pass through the weftless and short-cut fiber plies in Figure 4a, where t_s is the thickness of the short-cut fiber ply, t_u is the thickness of the weftless ply, and D_F is the diameter of the needled region. Based on the abundance of microstructural images, the average thickness of weftless and short-cut fiber plies, and the diameter of the needled region can be determined by measurement. In Figure 4b, the short-cut fiber ply lies between the 0° and 90° weftless ply, and many pores exist in the weftless and short-cut fiber ply. Figure 4c presents the morphology of the pyrolytic carbon under the polarized light microscope. It is found that the pyrolytic carbon grows outward around the fiber. By measurement, the extinction angle of the pyrolytic carbon is about 21°, which proves that the pyrolytic carbon is high-textured. According to Refs. [24,25], the carbon plane of the high-textured pyrolytic carbon is parallel to the fiber axis, which makes the pyrolytic carbon anisotropic. Figure 4d–f are grayscale histograms corresponding to Figure 4a–c, respectively. Based on these histograms, pore locations can be determined accurately.

The Archimedes drainage method is easy to operate for measuring the porosity, and it is convenient, safe, and environmentally friendly. However, it cannot be used to determine the pore volume distribution and is applied to samples attacked by water. Because the pore shape and distribution have little effect on the elastic properties of composites [26,27,28], the Archimedes drainage method is adopted to obtain the total porosity of needled C/C composites. According to the standard of ASTM-C20, the measured porosity was 8.9% by the Equation (1).

$$V_p={\displaystyle \frac{m_1-m_0}{m_1-m_2}}$$

(1)

where m₀ and m₁ denote the mass of the dried and water-infiltrated sample, respectively; m₂ is the mass of the sample in the water.

The tensile moduli of needled carbon/carbon composites in both directions (perpendicular to and along the needled fiber) are tested according to ASTM C1275-95 [23], where the specimen size is 75 mm × 6.5 mm × 3.5 mm and the gauge length is 15 mm, as shown in Figure 5. The CMT5340-30kN universal testing machine is used, and the loading speed is 0.5 mm/min during the static tensile process. To ensure the accuracy of the test results, the number of specimens in each direction is five.

3. Multiscale Modeling Scheme

Needled carbon/carbon composites consist of carbon fibers, pyrolytic carbon, pores, and interphases between the fibers and matrices, and there are different microstructural characterizations among the different plies. In the weftless ply, the carbon fibers are unidirectionally distributed, so the composite in this ply exhibited very high modulus and strength in the fiber axis direction. In the short-cut fiber ply, the carbon fibers have a random distribution in space, which makes the composite generally show isotropic material properties. In the needled region, the carbon fibers come from the short-cut fiber ply and are aligned along the needle direction. In this paper, a sequential multiscale method is proposed to accurately calculate the effective elastic properties of the needled C/C composite. This method is a bottom-up calculation method to achieve the prediction of effective properties of composites by bringing homogenization results of microscale models into the macroscale model, as shown in Figure 6.

To obtain the effective properties of models at different scales, numerical homogenization methods are used, which are suitable for dealing with composites with complex microstructures. Given that the composite consists of N-phase components and the volume occupied by the nth component is V_r, the average micro-field variable (strain, stress, etc.) for each phase is expressed as follows:

$${〈\mathit{f}〉}_r={\displaystyle \frac{1}{V_r}}{\displaystyle \underset{V_r}{\int }{\mathit{f}}_{r}dV},$$

(2)

where f denotes the micro-field variable and $〈\cdot 〉$ denotes the effective property. Since the total volume of the composite satisfies the equation $V={\displaystyle \sum _{r=0}^{N-1}{V}_i}$, the effective property of the composite is as follows [29]:

$$〈\mathit{f}〉={\displaystyle \frac{1}{V}}{\displaystyle \underset{V}{\int }\mathit{f}dV}={\displaystyle \sum _{r=0}^{N-1}\frac{1}{V}{\displaystyle \underset{V_r}{\int }{\mathit{f}}_{r}dV}}={\displaystyle \frac{V_r}{V}}{\displaystyle \sum _{r=0}^{N-1}\frac{1}{V_r}{\displaystyle \underset{V_r}{\int }{\mathit{f}}_{r}dV}}={\displaystyle \sum _{r=0}^{N-1}\frac{V_r}{V}{〈\mathit{f}〉}_r}={\displaystyle \sum _{r=0}^{N-1}{c}_r{〈\mathit{f}〉}_r},$$

(3)

where c_r denotes the volume fraction of the rth component that satisfies the equation ${\displaystyle \sum _{r=0}^{N-1}{c}_r}=1$.

Based on the numerical homogenization method, the average stress and average strain of the RVE model are expressed as follows in the post-processing of the software Abaqus2020.

$$〈{\sigma }_{ij}〉={\displaystyle \frac{1}{V}}{\displaystyle \sum _{k=1}^{N_{ele}}{V}_e\left({\displaystyle \sum _{s=1}^{N_{\mathrm{int}}}{\sigma }_{ij}(y_I)\cdot J(y_I)\cdot W(y_I)}\right)}={\displaystyle \frac{1}{V}}{\displaystyle \sum _{s=1}^{N_{\mathrm{int}}}{\sigma }_{ij}(y_I)\cdot IVOL(y_I)},$$

(4)

$$〈{\epsilon }_{ij}〉={\displaystyle \frac{1}{V}}{\displaystyle \sum _{k=1}^{N_{ele}}{V}_e\left({\displaystyle \sum _{s=1}^{N_{\mathrm{int}}}{\epsilon }_{ij}(y_I)\cdot J(y_I)\cdot W(y_I)}\right)}={\displaystyle \frac{1}{V}}{\displaystyle \sum _{s=1}^{N_{\mathrm{int}}}{\epsilon }_{ij}(y_I)\cdot IVOL(y_I)},$$

(5)

where V is the RVE model volume, N_ele is the number of elements, N_int is the number of integration points, V_e is the elemental volume, J(y_I) and W(y_I) are the Jacobian matrix, and weight at the integration point y_I in the element e, respectively, and IVOL(y_I) denotes the volume of the integration point y_I, satisfying the equation IVOL(y_I) = V_e·J(y_I)·W(y_I). According to Hooke’s law, the effective stiffness of the composite is as follows:

$$〈C_{ijkl}〉={\displaystyle \frac{〈{\sigma }_{ij}〉}{〈{\epsilon }_{kl}〉}}$$

(6)

3.1. Ply Models Considering the Interphase and Pyrolytic Carbon Morphology

To accurately characterize the microstructures of needled C/C composites at the microscale, the interphase between the fiber and matrix and pyrolytic carbon morphology should be considered. According to our previous work [30], the thickness and modulus of the interphase are 0.163 μm and 5.84 GPa, respectively. The Poisson’s ratio of the interphase is assumed to be 0.35, which is the same as the Poisson’s ratio of the matrix. For the pyrolytic carbon, it has a feature of growing outwards around the fiber. In addition, pores are distributed in pyrolytic carbon matrices. To numerically characterize the pyrolytic carbon matrix, the treatment of transferring the pore–matrix system into an equivalent matrix is performed. Combining the Voigt and Reuss methods, the effective property of the equivalent matrix is calculated by Equation (7) [31].

$${〈C〉}_{em}=0.5\times \left({〈C〉}_V+{〈C〉}_R\right)=0.5\times \left(\left(V_0{C}_0+V_1{C}_1\right)+{\left(V_0{C}_0^{-1}+V_1{C}_1^{-1}\right)}^{-1}\right)$$

(7)

where the subscripts em, V, R, 0, and 1 represent the equivalent matrix, Voigt method, Reuss method, pyrolytic carbon, and pore.

3.1.1. Weftless Ply

The fibers in the weftless ply show a unidirectional distribution, which belongs to a kind of unidirectional fiber-reinforced composite system. The main methods currently used to create RVE models for unidirectional fiber-reinforced composites include the random sequence adsorption algorithms [32,33], Monte Carlo methods [34,35], and nearest neighbor algorithms [36]. Among these methods, the nearest neighbor algorithm can maximize the use of space and ensure the randomness of fiber distribution to estimate the RVE model with a fiber volume fraction above 50%. Generally, the actual fiber volume fraction in unidirectional composites is lower than the maximum fiber volume fraction obtained by the nearest neighbor algorithm. Therefore, a modified random fiber removal algorithm is used in this study to establish the RVE model for the weftless ply. This algorithm consists of two main steps: (1) constructing the RVE model with the maximum fiber volume fraction based on the nearest neighbor algorithm, and (2) removing some fibers arbitrarily from the generated fibers to reach the target fiber volume fraction. The detailed completion process of the random fiber removal algorithm can be seen in Ref. [30].

For the morphology of the pyrolytic carbon matrix, a Voronoi tessellation method is used [31]. In this method, a fiber is considered as a seed in a Voronoi cell. Let l₀ represent a reference fiber, l₁ is any fiber except the reference fiber. If a matrix element p lies in a Voronoi cell whose seed is the reference fiber, the distance between this element and the reference fiber needs to meet Equation (8).

$$d(p,l_0)<d(p,l_1).$$

(8)

If two endpoints of a fiber are defined as P₀(x₀, y₀, z₀) and P₁(x₁, y₁, z₁), respectively, and the center point of the matrix element is M(x_m, y_m, z_m), the distance d can be calculated by the Equations (9) and (10):

$$\lambda ={\displaystyle \frac{\overrightarrow{P_{0}M}·\overrightarrow{P_0{P}_1}}{{\left|\overrightarrow{P_0{P}_1}\right|}^2}}$$

(9)

$$d=\sqrt{{(x_m-\lambda x_1+(\lambda -1)x_0)}^2+{(y_m-\lambda y_1+(\lambda -1)y_0)}^2+{(z_m-\lambda z_1+(\lambda -1)z_0)}^2}$$

(10)

Based on Equations (8)–(10), the equivalent matrix is partitioned into some Voronoi cells by selecting the elements. Figure 7 gives the RVE model of the weftless ply and the morphology of the equivalent matrix after treatment, where the model size is L₁ × L₂ × L₃ = 56 μm × 56 μm × 42 μm [37].

3.1.2. Short-Cut Fiber Ply

The short-cut fiber ply is composed of carbon fiber webs, interphases, pyrolytic carbon matrices, and pores, in which the pores and pyrolytic carbon are homogenized into an equivalent matrix. The carbon fibers are obtained by cutting long fibers into short fibers with a certain length. The short fibers in the short-cut fiber ply have a random distribution, so the short-cut fiber ply can be regarded as an isotropic composite. In this section, the random sequential adsorption algorithm is used to establish the RVE model of the short-cut fiber ply [38]. According to Ref. [31], the model length is larger than two times the fiber length for random fiber-reinforced composites. Thus, the model size of the short-cut fiber ply is selected as L × L × L = 84 μm × 84 μm × 84 μm, as shown in Figure 8a. In addition, a Voronoi tessellation method is applied to characterize the morphology of the pyrolytic carbon. The process of the Voronoi tessellation method has been given in Section 3.1.1. Figure 8b gives the illustration of the partitioned equivalent matrix. It is found that each fiber only has a Voronoi cell, and each element in this cell is the nearest to this fiber. It is worth noting that the material orientation of each Voronoi cell is the same as that of the carbon fiber in this cell, in order to characterize the anisotropic pyrolytic carbon.

3.1.3. Needled Fiber Region

During the needling process, some fibers of the short-cut fiber ply are brought into the stacking direction by the needle plate, so that the needled fibers exhibit unidirectional distribution. In addition, the length of the short-cut fibers is generally between 3 and 25 mm, which is larger than the needle thickness. Therefore, the needled region can be regarded as a unidirectional composite. In this region, the fibers are tightly arranged and exhibit a square or square-hexagonal distribution. The previous results [39] showed that the RVE model constructed based on the square-hexagonal arrangement is closer to the experimental results. Thus, the RVE model with hexagonally arranged fibers is established, as shown in Figure 9. To avoid the poor convergence caused by the contact or too small distance between fibers, the minimum distance between fibers is defined as 0.07d (d denotes the fiber diameter) in this section. According to the calculation, the RVE model size of the needled fiber region is L₁ × L₂ × L₃ = 13.72 μm × 7.82 μm × 42 μm.

3.2. Macroscale Model

Needle density is the number of needle punches per square centimeter and can be calculated by the following equation:

$$D_n={\displaystyle \frac{D_p\times f}{v}},$$

(11)

where D_n denotes the needling density, D_p is the implantation density of the needle on the needle plate (punch/cm), f is the needling frequency (punch/min), and v denotes the fabric speed (cm/min) of the machine. Considering the complexity and randomness of the needling process, this paper assumes that the needling density and the implantation density are equal, and the other parameters are defaulted to 1. Proximity to the actual microstructure is achieved by randomly generating the location of the needled fiber region, but it is necessary to ensure that there is no overlap between the needle punches. Figure 10 represents the schematic diagram of the macroscale model of needled C/C composites, where h_d is the needling depth, t_u is the thickness of the weftless ply, t_s is the thickness of the short-cut fiber ply, H is the thickness of the macroscale model, and L is the length of the macroscale model. The values of L and H are determined in Section 5.1 by the parameter investigation.

4. Material Property, Ply Parameters, and Boundary Conditions

The material, component information, and boundary conditions need to be determined before calculating the effective elastic properties of needled C/C composites. In addition, the ply parameters containing the fiber volume fraction and porosity of each part, the thicknesses of the weftless and short-cut fiber plies, and the diameter of the needled region also need to be defined.

4.1. Material Property

Needled C/C composites consist of carbon fibers, pyrolytic carbon, pores, and interphases between fibers and matrices. According to Hill’s law, the inner continuity of the RVE model needs to be ensured when using the average stress/strain method [40]. Therefore, pores are considered soft materials with very low stiffness to disregard the effects of pores on the effective elastic properties of needled C/C composites [41].

Table 1. Elastic constants of each component [6,42,43,44].

	Elastic Modulus		Shear Modulus		Poisson’s Ratio
	EL (GPa)	ET (GPa)	GLT (GPa)	GTT (GPa)	νLT	νTT
T700	237	15.8	23	15	0.01	0.37
Pyrolytic carbon	27	13	9.18	4.81	0.47	0.35
Pore	2.5 × 10−5	2.5 × 10−5	1.0 × 10−5	1.0 × 10−5	0.25	0.25
Interphase	5.84	5.84	2.16	2.16	0.35	0.35

lists the material properties of each component.

4.2. Ply Parameters

The weftless ply has a surface density of 240.5 g/m², the short-cut fiber ply has a surface density of 120 g/m², and the density of the needled preform is 0.44 g/cm³. By derivation, the fiber volume fraction in the weftless ply is as follows:

$$V_{\mathit{f}w}={\displaystyle \frac{{\rho }_{sur\mathit{f}ace\_w}}{{\rho }_0}}{V}_{\mathit{f}0},$$

(12)

where V_f0 represents the fiber volume fraction in the needled C/C composites, ρ_{surface_w} is the surface density of the weftless ply, and ρ₀ is the density of the needled preform. The needled fibers come from the short-cut fiber ply. Before needling, the proportion of fibers in the short-cut fiber ply to the volume of the short-cut fiber ply is

$$V_{\mathit{f}s0}={\displaystyle \frac{{\rho }_{sur\mathit{f}ace\_s}}{{\rho }_0}}{V}_{\mathit{f}0},$$

(13)

where ρ_{surface_s} is the density of the short-cut fiber ply. For the needled preform, the needle plate starts to work when one short-cut fiber ply and one weftless ply are stacked. Thus, if the effect of needled fibers is considered, the actual fiber content in the short-cut fiber ply is expressed as

$$V_{\mathit{f}s}={\displaystyle \frac{V_{\mathit{f}s0}\times t_s\times L_0\times L_0-(t_s+t_u)\times \pi \times \frac{D_{\mathit{f}}^2}{4}\times V_{\mathit{f}n}}{(t_s\times L_0\times L_0-t_s\times \pi \times \frac{D_{\mathit{f}}^2}{4})}}$$

(14)

where V_fn is the fiber volume fraction in the needled fiber region, the value of which is 75% in this study [6]; D_F is the diameter of the needled region; and L₀ is the side length of the region containing a single needled fiber region, one short-cut fiber ply, and one weftless ply in the direction perpendicular to the needled fiber. Equation (15) is used to obtain the value of L₀, and N_D is the needling density.

$$L_0=\sqrt{{\displaystyle \frac{100}{N_D}}}.$$

(15)

The fiber volume fraction of needled C/C composites is 25%. According to Equations (12)–(15), the fiber volume fraction in the weftless ply is 41.4% and the fiber volume fraction in the short-cut fiber ply is 11.2%. By the optical images at different locations, the average thickness of the weftless ply is 0.28 mm, the average thickness of the short-cut fiber ply is 0.34 mm, and the average diameter of the needled region is 0.6 mm. In addition, the porosity in the short-cut fiber ply is 10.3%, the porosity in the weftless ply is 4.1%, and the porosity in the needled fiber region is 3.1% based on SEM images using the image processing technique [45], as shown in

Table 2. Geometrical parameters of each part of needled C/C composites.

	Weftless Ply	Short-Cut Fiber Ply	Needled Fiber Region
Fiber volume fraction	41.4%	11.2%	75%
Porosity	4.1%	10.3%	3.1%
Thickness or diameter	0.28 mm	0.34 mm	ϕ0.6 mm

.

4.3. Boundary Conditions

Generally, it is assumed that the composites consist of representative cells arranged periodically [46]. Thus, periodic boundary conditions are used to constrain the freedom of the two opposite nodes on the boundary of the RVE model at different scales to ensure displacement and force continuity. When the displacement in the periodic boundary conditions is constrained, the force constraint condition is automatically satisfied [47]. For two opposite nodes in the RVE model, the displacement boundary condition can be expressed as

$$u_i^{j+}-u_i^{j-}={\overline{\epsilon }}_{ik}(x_k^{j+}-x_k^{j-}),$$

(16)

where the subscripts i and k denote the direction and are defined as 1, 2, and 3. The superscript j is the jth node pair on the opposite surfaces. u and $\overline{\epsilon }$ are the displacement and average strain. To realize the displacement boundary condition in Equation (16), the linear multipoint constraint equation [38] is used in this study, as shown in Equation (17).

$$a_1{u}_i^N+a_2{u}_j^M+\dots +a_n{u}_l^Z=0,$$

(17)

where $u_i^N$ is the nodal displacement of node N at the degree of freedom and a_n is the coefficient that constrains the relative motion of the nodes. To further implement the periodic constraints in the Abaqus software, reference points need to be introduced for defining the external loads, as shown in Equation (18).

$$u_i^{j+}-u_i^{j-}=u_i^{R{P}_j}$$

(18)

where the superscript RP denotes the reference point. In this paper, three reference points (RP₁, RP₂, and RP₃) are introduced, where RP₁ is used to impose tensile deformations in three directions, and RP₂ and RP₃ are used to impose three shear deformations, respectively. If a degree of freedom is constrained, this degree of freedom is deleted to avoid over-constraint.

5. Results and Discussion

Needled C/C composites have hierarchical structures containing weftless plies, short-cut fiber plies, and needled fiber regions. To validate the proposed multiscale model, the model’s effectiveness at different scales needs to be investigated. In addition, the relationship between porosity and the elastic properties of needled C/C composites is quantitatively characterized. It is worth noting that the following slash symbol indicates that the value of the variable is the same on both sides of the symbol.

5.1. Model Validation

(i)

Microscale models;

Before validation, mesh sensitivity needs to be determined. The mesh size has a significant effect on computational convergence and accuracy. For microscale models, the previous study [48] showed that the mesh size L_e = L/45 (L is the maximum side length of RVE) can provide a very close result and acceptable computation cost. Thus, the microscale models are meshed using this mesh size in this paper. The computation time of each model is about 20 min.

The weftless ply and the short-cut fiber ply are regarded as unidirectional C/C composites and short-cut C/C composites, respectively. Thus, experimental results of the unidirectional C/C composites and short-cut C/C composites are used to verify the proposed methods. According to the experimental results from Refs. [30,31], the transverse elastic modulus of unidirectional C/C composites is 7.37 GPa; and the elastic modulus of short-cut C/C composites is 8.14 GPa. Based on the proposed methods in this paper, the predicted transverse effective elastic modulus of unidirectional C/C composites is 7.58 GPa, and the predicted effective elastic modulus of short-cut C/C composites is 8.80 GPa. The results show that the maximum relative error for these two materials is 8.11%, which is in good agreement with the experimental results.

To validate the RVE model of the needled region, the Halpin–Tsai method, which has now been widely used to evaluate the elastic properties of composites with unidirectional fiber distribution [49,50], is used. The results of the needled model and Halpin–Tsai method are given in

Table 3. Comparison between the RVE model and Halpin–Tsai results.

Properties	Model (GPa)	Halpin–Tsai (GPa)	Err
E11	162.32	174.96	7.22%
E22/E33	11.46	11.43	0.26%
G12/G13	6.79	7.19	5.56%
G23	5.38	4.15	29.64%

, where the property in the 2–3 plane is isotropic. It is found that the predicted results of the proposed model agree with the Halpi–Tsai results, which indicates that the proposed model can predict the elastic properties of the needled fiber region.

(ii)

Macroscale model

Before validating the macroscale model, the mesh size and model size need to be determined by the parametric study. For the mesh sensitivity analysis, the macroscale model size is kept as L × L × H = 5 mm × 5 mm × 1.24 mm, while the mesh size is set to 0.1 mm, 0.15 mm, 0.2 mm, 0.25 mm, and 0.3 mm, respectively. The effects of the mesh size on the effective elastic properties of the macroscale model are given in Figure 11. The results show that when the mesh size exceeds 0.25 mm, the elastic properties tend to be stable. The relative error in Figure 11d is the average relative deviation between the five realizations of RVE models and the average value, which is expressed as follows:

$$Error={\displaystyle \frac{\max \left({\varphi }_i-\frac{1}{n}{\displaystyle \sum _{i=1}^n{\varphi }_i}\right)}{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\varphi }_i}}},$$

(19)

where ϕ denotes the predicted elastic properties of each realization, and n denotes the number of realization of models. According to the results in Figure 11d, it is found that the mesh size does not have a significant effect on the predicted results of the macroscale model, in which the maximum relative error is 1.89%. Therefore, considering the computational efficiency, the model mesh size is chosen as 0.25 mm.

When exploring the effect of model length (L) on the effective elastic properties of the macroscale model, the model thickness (H) is 1.24 mm, the mesh size is 0.25 mm, and L is set to 1.9 mm, 5 mm, 10 mm, 15 mm, and 20 mm, respectively. The results in Figure 12a–c show that the effective elastic properties of the macroscale model are kept unchanged when L exceeds 10 mm. Figure 12d demonstrates the relative error of the macroscale model with different L. The results show that the effects of the model length on the relative error can almost be disregarded. Therefore, the macroscale model length is selected as 10 mm in this paper.

While investigating the effect of the model thickness (H) on the effective elastic properties of the macroscale model, the length of the model is 10 mm, the mesh size is 0.25 mm, and the thickness varies from 1.24 mm to 6.20 mm. The variation of the effective elastic modulus of the macroscale model with the thickness is given in Figure 13a. The results show that when the model thickness is more than 3.72 mm, the effective elastic modulus almost remains unchanged. In Figure 13b, the effective shear modulus of the macroscale model increases with the increase in the model thickness but remains stable when the model thickness exceeds 3.72 mm. Figure 13c shows the variation in the Poisson’s ratio with the model thickness. It is found that the effective Poisson’s ratio remains unchanged when the model thickness exceeds 3.72 mm. In addition, the results in Figure 13d show that the model thickness does not have a significant effect on the relative error. Thus, the model thickness is 3.72 mm in this paper.

According to the above parametric investigation, the macroscale model size of the needled C/C composite is 10 mm × 10 mm × 3.72 mm, and the mesh size is selected as 0.25 mm. The computation time is about 1 min. Based on the proposed multiscale model of the needled C/C composite, the in-plane (the direction perpendicular to the needled fiber) modulus is 26.05 GPa, and the out-of-plane (the direction along the needle fiber) modulus is 17.38 GPa. By the in-plane and out-of-plane tensile tests, the average in-plane tensile modulus of the needled C/C composites is 27.62 ± 8.32 GPa and the average out-of-plane tensile modulus is 19.04 ± 5.24 GPa. According to Equation (20), the predicted relative error for the in-plane tensile modulus is 5.68% and that of the out-of-plane tensile modulus is 8.72% The results show that the proposed multiscale model for needled C/C composites is reliable in predicting the effective elastic properties.

$$Error={\displaystyle \frac{\left|R_{exp}-R_{\mathit{f}E}\right|}{R_{exp}}}$$

(20)

where R_exp and R_FE are the experimental and predicted results, respectively.

5.2. Effects of Porosity

Needled C/C composites contain plenty of irregular pores. These pores are not conducive to the improvement of the elastic properties of needled C/C composites. To quantify the effect of porosity on the elastic properties of needled C/C composites, this section explores the relationship between the porosity and the effective elastic properties of each model. In Figure 14, the effective elastic and shear moduli of each model decrease with the increase in porosity, showing a linear relationship with the porosity. When the porosity in the needled region increases from 0.017 to 0.059, the E₁₁ of the needled region decreases by 2.56%, E₂₂/E₃₃ decreases by 7.05%, G₁₂/G₁₃ decreases by 8.36%, and G₂₃ decreases by 10.65% in Figure 14(a1,a2). When the porosity of the weftless ply changes from 0.023 to 0.078, the E₁₁ decreases by 1.45%, E₂₂/E₃₃ decreases by 6.74%, G₁₂/G₁₃ decreases by 6.76%, and G₂₃ decreases by 7.72% in Figure 14(b1,b2). When the porosity of the short-cut fiber ply varies from 0.058 to 0.197, the elastic modulus of the short-cut fiber ply decreases by 15.3% and the shear modulus decreases by 15.29% in Figure 14(c1,c2). When the porosity of the composite varies from 0.05 to 0.17, E₁₁/E₂₂ of the composite decreases by 5.1%, E₃₃ decreases by 9.28%, G₁₂ decreases by 10.26%, and G₁₃/G₂₃ decreases by 11.85% in Figure 14(d1,d2). From the above analysis, it is shown that the porosity has little effect on the elastic properties along the fiber direction in the needled region and the weftless ply is small, while it has a significant effect on the remaining elastic properties. The pores cause an equal decrease in the elastic modulus and shear modulus of the short-cut fiber ply, which is because the fibers and pores show a random distribution in the needled C/C composites. In the macroscale model, the porosity has less effect on the E₁₁ and E₂₂, which is because these two directions are the fiber directions of the 0° and 90° weftless plies, respectively.

6. Conclusions

In this paper, a multiscale model is proposed to predict the effective elastic properties of needled C/C composites. At the microscale, the microstructures of the weftless ply, the short-cut fiber ply, and the needled region are analyzed to establish the RVE model of each ply. Meanwhile, the matrix region is divided by the modified Voronoi tessellation method to realize the numerical characterization of the pyrolytic carbon morphology and its anisotropic properties. At the macroscale, the homogeneous properties of the microscale model are introduced into the macroscale model by the sequential multiscale calculation method to predict the effective elastic properties of needled C/C composites. Considering that the needling process is complex and random, the macroscale model size is determined by parametric investigation. In addition, by exploring the effect of porosity on the effective elastic properties of different scale models, it is found that the effective elastic properties of each model decrease with the increase in porosity, but the extent of decrease in each direction is not consistent, while the porosity causes the same extent of decrease in the elastic and shear moduli of the short-cut fiber ply.