Progressive damage for 3D printed composites with variable fiber volume content and ply location

Abstract

At present, 3D printed continuous carbon fiber reinforced composites (CCFRCs) are widely used in aviation and automotive fields. Most of the existing research on 3D printed CCFRCs focuses on fiber volume content and rarely explores the specific fiber ply location. However, the latter is one of the most important factors affecting the strength and stiffness properties of materials. In this paper, the constitutive model and progressive damage model of 3D printed CCFRCs are established, concentrating on the fundamental aspects of fiber volume content and fiber ply location. The mechanical properties of 3D printed CCFRCs with different fiber volume content and fiber ply location were studied by tensile tests. The mapping relationship between mechanical properties and fiber volume content and fiber ply location was obtained based on experimental data. In addition, the failure morphology and microstructure of the sample at the fracture section were analyzed, and the failure process under different fiber angles was studied to reveal the failure mechanism. The progressive damage model was used to predict and analyze the performance of 3D printed CCFRCs with different fiber volume content and fiber ply location. The results show that the feasibility and effectiveness of the present model are verified through the comparison between experiment and simulation.

Introduction

The momentum of 3D printing technology in the composite materials industry is increasing. ¹ Due to the advantages of excellent mechanical properties and strong designability, 3D printed composite materials are widely used in industries such as aerospace, defense, sports, and automobiles.^2–4 Continuous wire fabrication is one of the 3D printing technologies,^5,6 which is used through the desktop 3D printer MarkTwo by controlling the extrusion trajectory of the continuous carbon fibers in the part as well as the layup sequence. ⁷ It provides a new idea for the preparation of CCFRCs whose mechanical properties change with changing fiber volume content and fiber ply location. ⁸

At present, there is extensive research on the mechanical analysis of 3D printed CCFRCs. At the fiber level, some researches.^9–12 studied the tensile and bending properties of composite materials by changing the fiber volume content and filling method. Chaco N et al.^13,14 believe that the fiber ply location has a certain impact on the tensile mechanical properties of composite materials. Zhengming Huang et al. ¹⁵ systematically studied the constitutive equation, deformation, failure and strength of continuous fiber, and short fiber and particle reinforced composites.

In terms of printing technology, Kentaro Sugiyama et al. ¹⁶ conducted optimization research on continuous fiber trajectories to obtain composite materials with variable fiber volume fraction and stiffness. Seyed HR Sanei et al. ¹⁷ studied the effect of stress concentration during 3D printing on the tensile properties of carbon fiber nylon composites. Arda Özen et al.^18–21 prepared samples with different printing process parameters through using FDM technology, studied the effects of layer thickness and fiber overlap rate on the mechanical properties of the specimens, and conducted tensile analysis of the printed samples prepared under different standards. Mahmoud Moradi et al. ²² systematically investigated the effect of the filling mode (IPs) on the specific mechanical response of parts manufactured by melt deposition modeling (FDM), revealing the correlation between the mechanical properties of the material and the print filling mode. Nathalie Ramos et al. ²³ analyzed the effect of different filling characteristics on the heat transfer of the material, and accurately predicted the stress.

In terms of the structural characteristics of 3D printed CCFRCs, some researches.^24–26 evaluated the elastic properties of fiber reinforced composite materials and predicted the elastic properties using the average stiffness (VAS) method, and a new method is proposed to study the failure of engineering structures through comparative analysis of different 3D printing topologies. GOH et al. ²⁷ compared the damage and failure mechanisms of carbon fiber reinforced and glass fiber reinforced composite materials. Khalid Saeed et al. ²⁸ used the average stiffness (VAS) method to evaluate the elastic properties of fiber reinforced 3D printing composites with different fiber volume content.

The above studies have investigated the mechanical properties, stress concentration, damage and failure mechanisms, and interlayer adhesion of 3D printed CCFRCs from different perspectives. However, the arrangement of fiber ply location has not been integrated into the material constitutive model to study the effects of different fiber volume contents and fiber layer positions on the mechanical properties of 3D printed CCFRCs. In addition, the constitutive model of 3D printed composite materials with variable fiber volume content and location has not been established. ²⁹

Therefore, in this paper, a constitutive model of 3D printed CCFRCs with different fiber volume contents and different fiber ply location was established through extensive performance tests on the stiffness and strength. The tensile damage process of the specimens under different fiber orientations was investigated, and the interface morphology and failure mechanism were analyzed by observing the fracture interface and microscopic morphology of the tensile specimens. A progressive damage model with two-dimensional Hashin damage model was established to predict the damage evolution of materials, which provided support for the design and mechanical analysis of controllable 3D printing composites.

Experimental procedure

Materials and equipment

The MarkTwo 3D printer was used in this study, as shown in Figure 1. The working process of the extrusion head is shown in Figure 2. ³⁰ There are two print heads for short fiber reinforced nylon onyx (matrix material) and fiber reinforcements (carbon fiber, Kevlar, and glass fiber). In this paper, continuous carbon fiber is chosen as the reinforcing material with a density of 1.4 g/cm³ and a diameter of 0.35 mm. The matrix material is a mixture of nylon and microcarbon fiber with a diameter of 1.75 mm.^31–33 Mechanical loads and shocks can be applied as its high Young’s modulus and tensile strength. In addition, the onyx can withstand high temperature environments (300°C). Onyx combines the toughness and durability of nylon with the dimensional stability and high strength of composite materials and is considered one of the most powerful 3D printing plastics.^34,35

Figure 1.

MarkTwo 3D desktop printer.

Figure 2.

3D printer working principle. ³³

Internal layer structure of two materials is shown in Figure 3. Onyx and carbon fiber layups were printed at hot end temperatures of 252°C and 278°C, respectively, and each layer of the entire specimen was printed with a layer thickness of 0.125 mm.

Figure 3.

Internal layer structure.

In order to explore the mechanical properties of 3D printed CCFRCs with different fiber volume contents and different fiber ply location, the specimens were designed as shown in Table 1. Each mechanical property parameter was tested with four different fiber volume contents and four different fiber ply location, and there were at least three specimens in each category to exclude the influence of incidental errors. The centroid distance D refers to the distance from the centroid of the carbon fiber layer to the axis of symmetry in the thickness direction of the specimen layer, as shown in Figure 4.

Table 1.

Specimen parameters.

Sample ID	Number of carbon fiber layers	Fiber volume content (%)	Central axis (D)	Layer thickness (mm)	Total number of layers
1	2	8.3	0	0.125	24
2	2	8.3	1	0.125	24
3	2	8.3	2	0.125	24
4	2	8.3	3	0.125	24
5	4	16.7	0	0.125	24
6	4	16.7	1	0.125	24
7	4	16.7	2	0.125	24
8	4	16.7	3	0.125	24
9	6	25.0	0	0.125	24
10	6	25.0	1	0.125	24
11	6	25.0	2	0.125	24
12	6	25.0	3	0.125	24
13	8	33.3	0	0.125	24
14	8	33.3	1	0.125	24
15	8	33.3	2	0.125	24
16	8	33.3	3	0.125	24

Figure 4.

Fiber laying position.

Sample preparation and characterization

The dimensions of the specimens with the measured mechanical parameters were listed in Table 2. The specimen size of unidirectional fiber reinforced composite material is shown in Figure 5. Before the specimens were printed, CatiaV5R20 was used for 3D modeling, and the STL file was imported into the web slicing system Eiger® to set up the position of the carbon fiber ply location and the number of layers.

Table 2.

Specimen size.

Load type	Load direction	Test performance	3D printed sample sizes (mm)	Reference standard
Tensile test	Longitudinal	Modulus	230*12.5*3	GB/T 3354-2014
		Strength
	Transverse	Modulus	170*25*3
		Strength
Shear test	In-plane shear	Modulus	230*25*3	GB/T 3355-2014
		Strength

Figure 5.

Tensile specimen size.

According to the recommendation of GB/T 3354-2014, 12.5*50*1.5 tabs were bonded at both ends of the specimen using DP460 epoxy structural adhesive, as shown in Figure 6(a). The test specimen was obtained with electric blast drying oven for 1.5 h, as shown in Figure 6(b) and (c). There are three printing directions of carbon fibers in the specimen, which are 0°, 90° and ±45°, as shown in Figure 7. Longitudinal tensile strength and tensile modulus of elasticity were obtained with defining 0° fiber direction, transverse tensile strength and tensile modulus of elasticity were obtained with defining 90° fiber direction, and in-plane shear strength and shear modulus were obtained with defining ±45° fiber direction in the tensile test.

Figure 6.

Specimen preparation. (a) Reinforced tablets, (b) specimens, and (c) electric blast dryer.

Figure 7.

Specimen fiber direction.

Mechanical test setup

For the mechanical analysis of 3D printed CCFRCs, a significant number of uniaxial tensile tests are needed to characterize the mechanical properties of the material such as strength and elastic modulus. The tests are performed using a 100 kN electronic universal tensile machine of Labsans Material Testing (Labsans), as shown in Figure 8. In addition, the tensile test standard of unidirectional fiber reinforced composite material was used in the experiment, the gauge length is 50 mm, and the loading rate is 1 mm/min.

Figure 8.

Electronic universal testing machine.

In addition, the model of the electronic universal testing machine is LD26.105, and the strength and stiffness parameters of the specimen are tested according to the GB/T 3354-2014 standard. The calculation formula is as follows:

$\text{Tensile}\hspace{0.17em}\text{strength}\hspace{0.5em}{\sigma }_t=\frac{P_{\mathit{max}}}{wh}$

(1)

${\sigma }_t$ is tensile strength, measured in MPa. $P_{\mathit{max}}$ is the maximum load on the sample before failure, expressed in Newton (N). $w$ is the specimen width, expressed in millimeters. $h$ is the thickness of the specimen, in millimeters.

$\text{Tensile}\hspace{0.17em}\text{modulus}\hspace{1em}{\mathrm{E}}_t=\frac{\mathrm{\Delta }Pl}{wh\mathrm{\Delta }l}$

(2)

$E_t$ is the tensile modulus, measured in MPa. $l$ is the gauge length of the extensometer, in millimeters. $\mathrm{\Delta }P$ is the incremental load, expressed in Newton (N). $w$ is the specimen width, expressed in millimeters. $h$ is the thickness of the specimen, in millimeters. $\mathrm{\Delta }l$ is the increment of deformation over the length of the extensometer scale, expressed in millimeters.

Experimental results and discussion

Experimental results

Effect of fiber volume content and ply location on strength properties

In order to establish the mechanical analysis model of 3D printed CCFRCs, it is necessary to characterize the tensile strength and modulus characteristics of composite materials with different fiber volume content and fiber ply location.

The mapping relationship between tensile strength, fiber volume content, and fiber ply location was obtained through the test data, as shown in Figures 9–11. Longitudinal tensile strength has an advantage in the mechanical properties of CFRCS and is often used as the main load-bearing performance in design. Under the experimental conditions in Table 2, the average longitudinal tensile strength of 3D printed CCFRCs increases from 54.81 MPa to 264.37 MPa when the fiber volume content increases from 8.3% to 33.3% and the fiber centroid distance D increases from 0 to 3 (reference Figure 4), basically showing a linear increase, as shown in Figure 9(a). The centroid distance D has little effect on the longitudinal tensile strength when the fiber volume content increases from 8.3% to 25%. However, the maximum difference of longitudinal tensile strength is about 15 MPa with the change of D, as shown in Figure 9(b). Data of longitudinal tensile strength difference analysis for D is shown in Table 3.

Figure 9.

Longitudinal tensile strength properties of 3D printed CCFRCs. (a) The three-dimensional relationship diagram between longitudinal tensile strength, fiber volume content, and fiber ply location. (b) The relationship between longitudinal tensile strength and fiber volume content at different centroid distances.

Figure 10.

Transverse tensile strength properties of 3D printed CCFRCs. (a) The three-dimensional relationship between transverse tensile strength, fiber volume content, and fiber ply location. (b) The relationship between transverse tensile strength and fiber volume content at different centroid distances.

Figure 11.

Shear strength properties of 3D printed CCFRCs. (a) The three-dimensional relationship between shear strength, fiber volume content, and fiber ply location. (b) The relationship between shear strength and fiber volume content at different centroid distances.

Table 3.

Numerical analysis of longitudinal tensile strength by fiber volume fraction and centroid distance.

Fiber volume content (vol%)	8.3	16.7	25.0	33.3
Standard deviation	2.67	1.61	3.55	10.61
Data difference for D	5.93	4.00	8.05	29.16
True difference	3.26	2.39	4.50	18.55

In contrast, Figure 10(a) and (b) show that the transverse tensile strength of the material first increases and then decreases with the increase of fiber volume content and centroid distance. The transverse tensile strength reaches the maximum when the fiber volume content reaches 16.7%. Moreover, the transverse tensile strength is better when the centroid distance D is 0. Data of transverse tensile strength difference analysis for D is as shown in Table 4.

Table 4.

Numerical analysis of transverse tensile strength by fiber volume fraction and centroid distance.

Fiber volume content (vol%)	8.3	16.7	25.0	33.3
Standard deviation	0.34	0.68	0.69	0.16
Data difference for D	0.89	1.84	1.72	0.42
True difference	0.55	1.16	1.03	0.25

In Figure 11(a) and (b), the shear strength of the material increases with the fiber volume content when the centroid distance D is a constant. However, the shear strength decreases gradually with the increase of the centroid distance when the fiber volume content is 22%–33.3%, as shown in Table 5.

Table 5.

Numerical analysis of shear tensile strength by fiber volume fraction and centroid distance.

Fiber volume content (vol%)	8.3	16.7	25.0	33.3
Standard deviation	1.85	1.57	2.73	3.48
Data difference for D	4.55	4.28	6.78	9.42
True difference	2.71	2.72	4.05	5.93

Effect of fiber volume content and ply location on stiffness properties

Elastic modulus is an important performance parameter in engineering materials that measures the ability of objects to resist elastic deformation. Figures 12–14 show the relationship between the elastic modulus, the fiber volume content, and fiber ply location under different loads. The results show that the average longitudinal tensile modulus increases from 1766.47 MPa to 4454.01 MPa when the fiber volume content is from 8.3% to 33.3%, and the centroid distance D is from 0 to 3, as shown in Figure 12(a). The longitudinal tensile modulus is negatively correlated with the centroid distance when the fiber volume content is 8.3%–20% and is positively correlated with the centroid distance when the fiber volume content is 20%–33.3%, as shown in Figure 12(b). The centroid distance has a great influence on the longitudinal tensile modulus, which ranges from 130.38 MPa to 230.96 MPa, as shown in Table 6.

Figure 12.

Longitudinal tensile modulus properties of 3D printed CCFRCs. (a) The three-dimensional relationship diagram between longitudinal tensile modulus, fiber volume content, and fiber ply location. (b) The relationship between longitudinal tensile modulus and fiber volume content at different centroid distances.

Figure 13.

Transverse modulus properties of 3D printed CCFRCs. (a) The three-dimensional relationship between transverse modulus, fiber volume content, and fiber ply location. (b) The relationship between transverse tensile modulus and fiber volume content at different centroid distances.

Figure 14.

Shear modulus properties of 3D printed CCFRCs. (a) The three-dimensional relationship between shear modulus, fiber volume content, and fiber ply location. (b) The relationship between shear modulus and fiber volume content at different centroid distances.

Table 6.

Numerical analysis of longitudinal tensile modulus by fiber volume fraction and centroid distance.

Fiber volume content (vol%)	8.3	16.7	25.0	33.3
Standard deviation	84.73	93.19	129.75	63.70
Data difference for D	215.11	247.82	360.70	173.63
True difference	130.38	154.63	230.96	109.93

Figure 13(a) and (b) show that the transverse tensile modulus first increases and then decreases with the increase of fiber volume content and centroid distance. The transverse tensile modulus peaks when the fiber volume content is 23%. The transverse tensile modulus is reduced by 16.57% when the fiber volume content is 23% and the centroid distance D is from 0 to 2. The centroid distance has a great influence on the transverse tensile modulus, as shown in Table 7.

Table 7.

Numerical analysis of transverse tensile modulus by fiber volume fraction and centroid distance.

Fiber volume content (vol%)	8.3	16.7	25.0	33.3
Standard deviation	70.19	41.95	61.79	16.11
Data difference for D	178.69	101.66	171.00	40.32
True difference	108.50	59.71	109.21	24.21

In the shear test, Figure 14(a) and (b) show that the shear modulus increases significantly when the fiber volume content is between 6.7% and 25%, while the shear modulus decreases slightly when the fiber volume content is between 25% and 33.3%. Shear modulus trend is uncertain when D is from 0 to 3. However, the degree of influence of centroid distance on shear modulus decreases with the increase of fiber volume content, as shown in Table 8.

Table 8.

Numerical analysis of shear tensile modulus by fiber volume fraction and centroid distance.

Fiber volume content (vol%)	8.3	16.7	25.0	33.3
Standard deviation	46.74	27.01	29.32	22.87
Data difference for D	115.39	75.14	73.69	62.22
True difference	68.64	48.13	44.37	39.35

Failure morphology of 3D printed CCFRCs

Figure 15 shows the damage modes of 3D printed CCFRCs specimens from macro to micro. Figure 15(a) shows the fracture morphology of specimen with the fibers printed in the 0° direction under the tensile loading. The surface flakes off with the pulling out of the lamination between the fiber and the matrix. This is due to the process characteristics of 3D printed Fused Filament Fabrication (FFF), which leads the possibility of defects when stacking layers and the inability to obtain a tight interface between neighboring layers.

Figure 15.

Fracture failure diagram of the specimen. (a) Tensile failure in the 0° direction, (b) tensile failure in the 90° direction, and (c) shear failure in the ±45° direction.

In contrast with Figure 15(b), the specimen with fiber printed in 90° direction is neatly disconnected from along the boundary between the matrix (onyx) and the carbon fiber under tensile load. There are no separated bare fiber bundles in the cross section, which also indicates that the impregnation between the fiber and the matrix is excellent, and the failure mode is brittle failure. 3D printed CCFRCs under longitudinal tensile loading undergo fiber damage with fiber printed in ±45° direction and fiber peels from the matrix with the transverse cracking of the matrix, indicating that the fibers are the main bearers of material failure, as shown in Figure 15(c).

Discussion

Constitutive model of 3D printed CCFRCs

The fitting function between the mechanical property parameters, fiber volume content, and fiber position was obtained based on the experimental data, as shown in Table 9 where $V_f$ is the carbon fiber volume content and D is the fiber centroid distance. $E_1$ is the elastic modulus in the fiber direction, $E_2$ is the elastic modulus perpendicular to the fiber, and $G_{12}$ is the Shear modulus in the plane. Since it is difficult to obtain stable transverse strain in the experiment, Poisson’s ratio of ${\nu }_{12}$ with different fiber volume content is obtained by mixing law. ³⁶ The formula is shown in (3).

${\nu }_{12}={\nu }_{12f}+{\nu }_{12m}·\left(1-{\nu }_f\right)$

(3)

where ${\nu }_{12f}$ is selected as 0.3 and ${\nu }_{12m}$ is selected as 0.38.^37,38

Table 9.

Material properties of 3D printed continuous fiber reinforced composites.

Material properties	Value	Units	Determination coefficient
Modulus (fiber direction), $E_1$	$E_1\left(V_f,D\right)=143.3V_f+54.6D-0.55{V_f}^2-5.3V_{f}D+19.24D^2+592.5$	MPa	0.98
Modulus (transverse direction), $E_2$	$E_2\left(V_f,D\right)=25.24V_f+128.7D+1.3{V_f}^2$ $-15.76V_{f}D+5.1D^2+0.049{V_f}^3+0.3{V_f}^{2}D+0.33V_f{D}^2+429.9$	MPa	0.92
Shear modulus, $G_{12}$	$G_{12}\left(V_f,D\right)=-31.95V_f-139.4D+3.4{V_f}^2+$ $6.12V_{f}D+36.21D^2-0.065{V_f}^3-0.057{V_f}^{2}D-1.25V_f{D}^2+1057$	MPa	0.95
Poisson’s ratio, ${\nu }_{12}$	${\nu }_{12}\left(V_f\right)=0.3V_f+0.38\left(1-V_f\right)$		Mixing rule 39
Tensile strength (fiber direction), $X_T$	$X_T\left(V_f,D\right)=3.69V_f+0.6D+0.11{V_f}^2‐0.19V_{f}D+0.9D^2+16.42$	MPa	0.95
Tensile strength (transverse direction), $Y_T$	$Y_T\left(V_f,D\right)=-0.94V_f-0.0078D+0.08{V_f}^2$ $-0.11V_{f}D+0.37D^2-0.002{V_f}^3+0.003{V_f}^{2}D-0.007V_f{D}^2+27.3$	MPa	0.98
Shear strength, $S_{12}$	$S_{12}\left(V_f,D\right)=1.72V_f-1.64D-0.019{V_f}^2-0.105V_{f}D+0.64D^2+14.08$	MPa	0.94

The composite material studied in this paper is an anisotropic material, and the research is limited to the linear elastic stage. Therefore, the stress-strain relationship of the material is shown as formula (4). Since the thickness of 3D printed CCFRCs specimen in the Z-axis direction is extremely thin (the thickness size is much smaller than the length and width size), the stress in the z-direction changes little during the tensile testing. Assuming that z = 0, then the constitutive equation of composite under plane stress condition is equation (5).

$\left(\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\tau }_{33}\\ {\tau }_{31}\\ {\tau }_{12}\end{array}\right)=\left(\begin{array}{llllll}{Q}_{11}& Q_{12}& Q_{13}& 0& 0& 0\\ Q_{21}& Q_{22}& Q_{23}& 0& 0& 0\\ Q_{31}& Q_{32}& Q_{33}& 0& 0& 0\\ 0& 0& 0& Q_{44}& 0& 0\\ 0& 0& 0& 0& Q_{55}& 0\\ 0& 0& 0& 0& 0& Q_{66}\end{array}\right)\left(\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ {\gamma }_{23}\\ {\gamma }_{31}\\ {\gamma }_{12}\end{array}\right)$

(4)

When the thickness z is 0, consider ${\sigma }_3={\mathrm{\tau }}_{31}={\mathrm{\tau }}_{23}=0$

$\left(\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\tau }_{12}\end{array}\right)=\left(\begin{array}{lll}{Q}_{11}& Q_{12}& 0\\ Q_{21}& Q_{22}& 0\\ 0& 0& Q_{66}\end{array}\right)\left(\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\gamma }_{12}\end{array}\right)$

(5)

where

$Q_{11}=\frac{E_1}{1‐{\mathrm{\nu }}_{12}{\nu }_{21}}$

(6)

$Q_{12}=\frac{{\nu }_{12}{E}_1}{1‐{\mathrm{\nu }}_{12}{\nu }_{21}}$

(7)

$Q_{21}=\frac{{\nu }_{12}{E}_1}{1‐{\mathrm{\nu }}_{12}{\nu }_{21}}$

(8)

$Q_{22}=\frac{E_2}{1‐{\mathrm{\nu }}_{12}{\nu }_{21}}$

(9)

$Q_{66}={\mathrm{G}}_{12}$

(10)

$\frac{{\nu }_{21}}{E_2}=\frac{{\nu }_{12}}{E_1}$

(11)

In formula (5), ${\sigma }_1\hspace{0.17em}\text{and}\hspace{0.17em}{\mathrm{\epsilon }}_1$ represent longitudinal stress and longitudinal strain; ${\sigma }_2\hspace{0.17em}\text{and}\hspace{0.17em}{\mathrm{\epsilon }}_2$ represent transverse stress and transverse strain; and ${\tau }_{12}\hspace{0.17em}\text{and}\hspace{0.17em}{\mathrm{\gamma }}_{12}$ represent in-plane shear stress and in-plane shear strain, respectively. The constitutive equation of the 3D printed continuous fiber reinforced composites with different fiber volume content and fiber layup positions can be obtained through substituting the stiffness properties of the material $\left(E_1,E_2,{\nu }_{12},{\nu }_{21},\text{and}\hspace{0.17em}{G}_{12}\right)$ into equations (6)–(11), as shown in formula (12).

$\left(\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\tau }_{12}\end{array}\right)=\left(\begin{array}{ccc}\frac{E_1\left(V_f,D\right)}{1-{\nu }_{12}\left(V_f\right){\nu }_{21}\left(V_f\right)}& \frac{{\nu }_{12}\left(V_f\right)E_1\left(V_f,D\right)}{1-{\nu }_{12}\left(V_f\right){\nu }_{21}\left(V_f\right)}& 0\\ \frac{{\nu }_{12}\left(V_f,D\right)E_1\left(V_f,D\right)}{1-{\nu }_{12}\left(V_f\right){\nu }_{21}\left(V_f\right)}& \frac{E_2\left(V_f,D\right)}{1-{\nu }_{12}\left(V_f\right){\nu }_{21}\left(V_f\right)}& 0\\ 0& 0& G_{12}\left(V_f,D\right)\end{array}\right)\left(\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\gamma }_{12}\end{array}\right)$

(12)

Demonstration

In order to verify the constitutive equation, a specimen is designed and its printing parameters are shown in Table 10. The parameter properties of materials in Abaqus simulation are shown in Table 11. The validity of the constitutive model in this paper is verified by tensile test and finite element analysis.

Table 10.

Printing parameters.

	Value	Unit
Specimen size	240*12.5*3	mm
Carbon fiber orientation	0°	/
Number of carbon fiber layers	3	/
Fiber volume fraction	12.5	%
Centroid distance D	2.5	/
Layer thickness	0.125	mm
Total number of layers	24	/

Table 11.

Material properties.

	Value	Unit	Description	Reference
E1	2388.2	MPa	Homogenized lamina modulus of elasticity in direction 1 (L)	Table 9
E2	857.49	MPa	Homogenized lamina modulus of elasticity in direction 2 (T)	Table 9
V12	0.37	/	Poisson ratio in 12 plane	Table 9
G12	1011.5	MPa	Homogenized lamina shear modulus in 12 plane	Table 9

Figure 16(a) shows the simulated stress cloud map, where the maximum stress is 58.83 MPa, which is also the breaking limit stress. Figure 16(b) shows the displacement-load curve of the test specimen. The simulation results are basically consistent with the experimental results. In the simulation, the sample is basically linear before failure, that is, the tensile modulus remains unchanged because of the elastic phase in the materials. In the test, it was observed that when the load was increased to 1762.71 N, a small amount of matrix and interface of the sample began to crack. The specimen suddenly fractured when the load was added to 1886.15 N. In comparison, the failure load of the finite element simulation is 1847.18 N, which is 2.1% lower than the test results. The reason is that there are some corners that are not fully covered by the fiber layer in 3D printed CCFRCs. In the finite element analysis, the carbon fiber layer is treated as a completely covered uniform plane.

Figure 16.

Simulation model and experimental analysis: (a) stress nephogram and (b) test and simulation load-displacement curves.

The relative error of the final failure load of the test and simulation is less than 10%, which is in line with the expectation. The validity and correctness of the damage model are further verified.

Progressive damage analysis

Progressive damage model

Progressive damage analysis method is a finite element method based on the assumption that the damaged material structure can continue to bear the material properties after degradation. The progressive damage analysis method consists of three parts: stress analysis solution, material failure criterion, and degradation of damaged material properties. In this paper, Abaqus built-in 2D Hashin damage degradation criterion is used to study the progressive damage of 3D printed CCFRCs.

The 2D Hashin criterion can judge four types of failure, namely, fiber tensile failure, fiber compression failure, matrix tensile failure, and matrix compression failure. The expression is as follows:

$F_f^t={\left(\frac{{\sigma }_1}{X_T}\right)}^2+{\left(\frac{{\sigma }_6}{S_L}\right)}^2=1\hspace{1em}\text{Fiber}\hspace{0.17em}\text{tension}\hspace{0.17em}\left({\sigma }_1\ge 0\right)$

(13)

$F_f^c={\left(\frac{{\sigma }_1}{X_C}\right)}^2=1\hspace{1em}\text{Fiber}\hspace{0.17em}\text{compression}\left({\sigma }_1<0\right)$

(14)

$F_m^t={\left(\frac{{\sigma }_2}{Y_T}\right)}^2+{\left(\frac{{\sigma }_6}{S_L}\right)}^2=1\hspace{0.5em}\text{Matrix}\hspace{0.17em}\text{tension}\left({\sigma }_2\ge 0\right)$

(15)

$F_m^c={\left(\frac{{\sigma }_2}{2S_T}\right)}^2+\left[{\left(\frac{Y_C}{2S_T}\right)}^2-1\right]\left(\frac{{\sigma }_2}{Y_C}\right)+{\left(\frac{{\sigma }_6}{S_L}\right)}^2=1\hspace{1em}\text{Matrix}\hspace{0.17em}\text{compression}\left({\sigma }_2<0\right)$

(16)

where $X_T,X_C,\text{and}\hspace{0.17em}{Y}_T,Y_C$ are longitudinal tensile strength, compressive strength and transverse tensile strength, compressive strength, respectively; ${\sigma }_1,{\sigma }_2,\text{and}\hspace{0.17em}{\sigma }_6$ are the stress in the main axis direction of a certain layer in a certain cell; and $S_L\hspace{0.17em}\text{and}\hspace{0.17em}{S}_T$ are the longitudinal shear strength and transverse shear strength.

If any of the conditions in equations (13)–(16) are met, the material is judged to have suffered initial damage. After the initial damage occurs, it is necessary to update the stiffness matrix to simulate the degradation behavior of composite materials. In this paper, the stiffness coefficient degradation method based on fracture toughness is adopted. The damage constitutive relation of composites can be expressed as:

$\sigma =C_d\hspace{0.17em}\epsilon$

(17)

where $\epsilon$ is strain and $C_d$ is the damage elastic matrix, whose expression is:

$C_d=\frac{1}{D}\left[\begin{array}{lll}\left(1-d_f\right)E_1& \left(1-d_f\right)\left(1-d_m\right){\nu }_{21}{E}_1& 0\\ \left(1-d_f\right)\left(1-d_m\right){\nu }_{12}{E}_2& \left(1-d_m\right)E_2& 0\\ 0& 0& \left(1-d_s\right)GD\end{array}\right]$

(18)

$d_f=\{\begin{array}{ll}{d}_{ft}& {\sigma }_1\ge 0\\ d_{fc}& {\sigma }_1<0\end{array}{d}_m=\{\begin{array}{ll}{d}_{mt}& {\sigma }_2\ge 0\\ d_{mc}& {\sigma }_2<0\end{array}$

(19)

$d_s=1-\left(1-d_{ft}\right)\left(1-d_{fc}\right)\left(1-d_{mt}\right)\left(1-d_{mc}\right)$

(20)

$D=1-\left(1-d_f\right)\left(1-d_m\right){\nu }_{12}{\nu }_{21}$

(21)

where $d_f$ is the current state of fiber damage; $d_m$ is the current state of matrix damage; $d_s$ is the current state of shear damage; and $d_{ft},d_{fc},d_{mt},{and\hspace{0.17em}d}_{mc}$ are the damage variables corresponding to fiber stretching, fiber compression, matrix stretching, and matrix compression, respectively. $E_1$ is elastic modulus in the direction of the fiber; $E_2$ is elastic modulus in the vertical direction of the fiber; $G$ is shear modulus; and ${\nu }_{12}$ and ${\nu }_{21}$ are Poisson’s ratios.

The degradation relationship shown in Figure 17(a) is adopted for the evolution of each damage variable of the four failure modes, and the damage variable is calculated according to the following formula:

${\delta }_{eq}^{ft}=L_c\sqrt{{〈{\epsilon }_{11}〉}^2+{\epsilon }_{12}^2}\hspace{0.17em}\text{Fiber}\hspace{0.17em}\text{tension}\hspace{0.17em}\left({\sigma }_1\ge 0\right)\hspace{0.17em}{\sigma }_{eq}^{ft}=\frac{〈{\sigma }_{11}〉〈{\epsilon }_{11}〉+{\epsilon }_{12}}{\frac{{\delta }_{eq}^{ft}}{L_c}}$

(22)

${\delta }_{eq}^{fc}=L_c〈-{\epsilon }_{11}〉\text{Fiber}\hspace{0.17em}\text{compression}\hspace{0.17em}\left({\sigma }_1<0\right)\hspace{0.17em}{\sigma }_{eq}^{fc}=\frac{〈-{\sigma }_{11}〉〈-{\epsilon }_{11}〉}{\frac{{\delta }_{eq}^{fc}}{L_c}}$

(23)

${\delta }_{eq}^{mt}=L_c\sqrt{{〈{\epsilon }_{22}〉}^2+\alpha {\epsilon }_{12}^2}\hspace{0.17em}\text{Matrix}\hspace{0.17em}\text{tension}\hspace{0.17em}\left({\sigma }_2\ge 0\right)\hspace{0.17em}{\sigma }_{eq}^{mt}=\frac{〈{\sigma }_{22}〉〈{\epsilon }_{22}〉+{\tau }_{12}{\epsilon }_{12}}{\frac{{\delta }_{eq}^{mt}}{L_c}}$

(24)

${\delta }_{eq}^{mc}=L_c\sqrt{{〈{\epsilon }_{22}〉}^2+\alpha {\epsilon }_{12}^2}\text{Matrix}\hspace{0.17em}\text{compression}\hspace{0.17em}\left({\sigma }_2<0\right)\hspace{0.17em}{\sigma }_{eq}^{mc}=\frac{〈-{\sigma }_{22}〉〈-{\epsilon }_{22}〉+{\tau }_{12}{\epsilon }_{12}}{\frac{{\delta }_{eq}^{mc}}{L_c}}$

(25)

here $L_c$ is the characteristic length, in this paper is the square root of the unit area and ${\delta }_{eq}^{ft}$ and ${\delta }_{eq}^{fc}$ are equivalent displacement of fiber tension and compression. ${\delta }_{eq}^{mt}$ and ${\delta }_{eq}^{mc}$ are the equivalent displacement of tension and compressive matrix. ${\sigma }_{eq}^{ft}$ and ${\sigma }_{eq}^{fc}$ are the equivalent stress of fiber tension and compression. ${\sigma }_{eq}^{mt}$ and ${\sigma }_{eq}^{mc}$ are the equivalent stress of tension and compressive matrix. $〈〉$ is a symbolic operator. For example, $〈\alpha 〉=\frac{\left(\alpha +\left|\alpha \right|\right)}{2}$ .

$d_i=\frac{{\delta }_{eq}^f\left({\delta }_{eq}-{\delta }_{eq}^0\right)}{\hspace{0.17em}{\delta }_{eq}\left(\hspace{0.17em}{\delta }_{eq}^f-{\delta }_{eq}^0\right)}$

(26)

In the formula, $i$ is the four damage variables of ft, fc, mt, and mc. ${\delta }_{eq}$ and ${\delta }_{eq}$ are equivalent displacement and equivalent stress. ${\delta }_{eq}^f$ is the failure displacement corresponding to complete failure. ${\delta }_{eq}^0$ is the equivalent displacement when some failure criterion reaches exactly 1. The relationship between damage variable and equivalent displacement is shown in Figure 17(b).

Figure 17.

Damage degeneration relation. (a) Relationship between equivalent displacement and equivalent stress. (b) The relationship between damage variable and equivalent displacement.

In addition, the specimen used in progressive damage has a small hole with a diameter of 1 in the middle, as shown in Figure 18. The printing parameters are shown in Table 12. In the process of establishing the finite element model, the idea of simplification is adopted, and the actual 24 layers are simplified to 8 layers. The top and bottom layers are the base material onyx, and the remaining 6 layers are carbon fiber reinforced layers (0°/45°/90°/–45°/90°/0°), as shown in Figure 19.

Figure 18.

Specimen size of progressive damage.

Table 12.

Printing parameters.

	Carbon fiber	Onyx
Fill mode	Isotropic	Isotropic
Fill density (%)	100	100
Fiber orientation	0°/45°/90°/–45°/90°/0°	±45°
Number of layers	6	18
Layer thickness (mm)	0.125	0.125
Centroid distance D	0	/
Fiber volume fraction (%)	25	/

Figure 19.

Structural diagram. (a) Specimen structure diagram and (b) finite element structure diagram.

The parametric properties of carbon fiber and matrix material layering in the finite element model are shown in Table 13 and 14. The mesh of the specimen X-Y plan is shown in Figure 20. The origin of coordinates is selected at the bottom left end of the laminate. The X-axis is the direction of external load, the Y-axis is the direction of orthogonal X-axis in the laminate surface, and the Z-axis is the direction of thickness of superposition of laminate layering. Due to the stress concentration around the circular hole and the local edge of the hole is the initial place of damage, the local mesh around the circular hole is encrypted. Reinforcement branch constraints are applied to the left end face, and displacement loads are applied to the right end, as shown in Figure 20.

Table 13.

Onyx lamina properties.

	Value	Unit	Reference
Elastic modulus	950	MPa	21
Poison’s ratio (υ)	0.38	/	21
Tensile strength	37.83	MPa	21
Fiber density	1.2	(g/cm3)	12

Table 14.

Carbon fiber lamina properties.

	Value	Unit	Description	Reference
E1	11486.1	MPa	Homogenized lamina modulus of elasticity in direction 1 (L)	Table 9
E2	2557.23	MPa	Homogenized lamina modulus of elasticity in direction 2 (T)	Table 9
V12	0.38	/	Poisson ratio in 12 plane	Table 9
G12	3186.6	MPa	Homogenized lamina shear modulus in 12 plane	Table 9
XT	159.14	MPa	Tensile strength in the fiber direction	Table 9
XC	108.46	MPa	Compressive strength in the fiber direction	a
YT	48.72	MPa	Tensile strength perpendicular to the fiber direction	Table 9
YC	79.56	MPa	Compressive strength perpendicular to the fiber direction	a
S12	64.94	MPa	In-plane shear strength	Table 9

^aValues acquired experimentally.

Figure 20.

Mesh image.

Result analysis

Figures 21 and 22 show the progressive damage distribution of fiber and matrix elements during tensile failure of the whole structure. This indicates that the unit has failed when the output variables HSNFTCRT and HSNMTCRT of the damaged unit are 1. If the output variable is less than 1, the failure state has not been reached. It can be seen from the fiber damage distribution diagram in Figure 21 that the maximum number of fiber failures occurs in the 0° layup when the ultimate load is reached. The 0° fiber tensile damage extends to the boundary of the composite laminates, the fibers are completely broken, and the whole structure has been penetrated. The fiber failure units in 45° and −45° layings are similar because the fiber laying Angle is rotated by 180°. Moreover, the fiber failure elements extend to the boundary of the laminate, indicating that the layer is completely broken. There are very few fiber failure units in the 90 ply, indicating that the 90° fiber does not bear a load in the X-direction stretch.

Figure 21.

Distribution of fiber failures in each layer during overall damage.

Figure 22.

Distribution of matrix failures in each layer during overall damage.

As can be seen from the damage distribution diagram of the matrix in Figure 22, a large area of the matrix fails in the ±45° paving layer, and the damage units have penetrated the whole structure, so the matrix will fail in large numbers when carrying loads. The matrix failure in 0° layering is relatively small, because the fiber is the main load-bearing material in 0° layering. As long as the fiber is not damaged in a large area, the matrix in the layering will not break. In the 90° layup, the matrix damage occurs before the fiber damage, so the failure matrix element that permeates the whole laminate appears.

Figure 23 shows the load-displacement curve of finite element calculation and test. In the test, the ultimate tensile load when the specimen was damaged was 2355.2 N, but after the ultimate load, the specimen suddenly broke and the load decreased instantaneously. The progressive damage method was used in the finite element calculation, and the damage degradation began when the load was 2246.29 N. However, the structure can still bear a little load after the failure, and the load-displacement curve changes greatly after the complete failure. The error of the finite element calculation and test results is 4.6%, and the error is less than 10%, which is in line with the expectation. Figure 24 shows the fracture diagram of the test piece. The main failure position is on both sides of the hole, and the failure mode is tensile failure. In addition, in the test, normal layered damage appeared at other positions of the hole, and the damage was not in a single form, and various damages affected each other. This is also the next step to improve the research content.

Figure 23.

Load-displacement curves for testing and simulation.

Figure 24.

Fracture diagram of test piece.

The comparison of the whole results shows that the analysis method is accurate and the analysis results are reliable. The rationality and accuracy of the progressive damage model established in this paper are verified.

In this paper, the tensile strength and modulus parameters of a large number of different fiber characteristics are obtained by testing 16 fiber arrangements, as shown in Table 15. As can be seen from the table, there are 3 independent variables involved, which are 4 fiber volume fractions, 4 fiber layup sequences and 4 fiber ply location. There are 4 dependent variables involved, namely, transverse and longitudinal tensile strength and modulus, which can obtain more comprehensive mechanical parameters than other literatures.^39–41

Table 15.

Comparative study of relevant literatures.

Mechanical property	Independent variable			Dependent variable
	Fiber content (%)	Ply location	Ply angle	Tensile strength	Tensile modulus	Flexural strength	Flexural modulus
This study	8.3, 16.7, 25, and 33.3	4	0°/90/±45°	Longitudinal and transverse	Longitudinal and transverse	/	/
Bartosz Nowinka1 42	15, 25,37, and 40	/	0°/45°	Longitudinal	Longitudinal	/	/
M.A. Caminero 43	0, 27.2, and 73.2	2	0°	Shear	/	/	/
Zhongsen Zhang 25	18.9, 19.5, 19.7, and 28.0	/	0°	Longitudinal	Longitudinal	√	√

The concept of centroid distance was first introduced in this paper. Compared with the specific ply location mentioned in other literature, this study covers a wider range of laying methods and has broader applicability. For the first time, this study employs an experimental data fitting method to derive a constitutive equation for the mechanical properties of 3D printed fiber reinforced composite materials in relation to fiber volume content and ply location. Compared with other literature that obtains the mechanical properties of fiber reinforced composite materials through constitutive equations, this study’s equation takes into account not only the influence of fiber volume content on mechanical properties but also establishes a bivariate mechanical parameter (fiber volume content and ply location) constitutive equation. This achieves a dual-variable influence factor prediction for the mechanical properties of 3D printed fiber reinforced composite materials regarding fiber volume content and ply location, making the predictions widely applicable,^25,28,42,43 as shown in Table 16.

Table 16.

Comparative study of relevant literatures.

	Factors affecting mechanical properties
	Fiber volume content	Fiber ply location	Fiber type
Present study	√	√	/
Garrett W. Melenka 28	√	/	/
E. Polyzos 42	√	/	/
Andrew N. Dickson 43	√	/	√
Vairavan Manikandan 25	/	√	/

In addition, the constitutive equations including fiber volume fraction (8.3%–33.3%) and fiber ply location (0–3) are established, which provide rich data support for 3D printing composites with variable stiffness and strength. In the study of progressive damage of 3D printed CCFRCs, the damage evolution process of specimens with fiber volume fraction of 25% and centroid distance of 0 was analyzed. Although only one fiber configuration was selected for damage analysis, it is highly representative and reveals the damage change law of the fiber and matrix units of each layer of the 3D printed composite.

Conclusions

In this paper, the mechanical properties of 3D printed CCFRCs with different fiber volume content and fiber ply location were experimentally studied. The specific findings and conclusions are as follows:

(1) In the longitudinal tensile test, the strength and stiffness of the material gradually increase with the increase of fiber volume content, but the change of centroid distance D has little effect on its properties. In the transverse tensile test, the strength of the material decreases with the increase of fiber volume content and centroid distance D.

(2) In the transverse tensile test, the strength of the material decreases with the increase of fiber volume content and centroid distance D. The stiffness increases first and then decreases with the increase of fiber volume content, and peaks when the fiber volume content is 23%. In addition, the stiffness is negatively correlated with the centroid distance D when the fiber volume content is constant.

(3) In the shear test, the strength decreases significantly with the increase of the centroid distance, and the strength reaches the maximum when the centroid distance is 0.

(4) Establishing the constitutive relationship between different fiber volume contents, ply locations, and the mechanical properties of composite materials allows for more precise design to meet specific performance requirements, enabling the optimization of material design. Researchers and engineers can predict the mechanical properties of new materials under specific conditions, reducing the number of experimental iterations and saving time and costs. For varying application demands, such as aerospace, automotive, and sports equipment, it’s possible to customize the performance of composite materials by adjusting the fiber content and ply location strategy to suit different loads and environmental conditions.