Prediction of the failure behavior of pseudo-ductile composites using a multi-scale finite element model

Abstract

The hybridization technique has recently been used to produce a new generation of composites called pseudo-ductile composites, which have shown higher failure strain compared to conventional composites, minimizing the risks of the occurrence of a catastrophic failure. The pseudo-ductility behavior in these composites is obtained by hybridization of fibers with high and low failure strains. In this study, a multi-scale finite element (FE) model incorporating micro and macro-scales is proposed to predict the failure behavior of pseudo-ductile composites. A micro-scale representative volume element (RVE), consisting of randomly distributed fibers, was generated using a Python code. Periodic boundary conditions (PBCs) were applied to the RVE generated with a periodic geometry. To account for fiber failure and ply fragmentation, the tensile strength of fibers was distributed based on the Weibull distribution function and a user-defined UMAT subroutine was developed. Tensile loading was then applied to the RVE to simulate the composite’s mechanical behavior. For validation, an RVE was developed based on experimental data from recent research on thinply and conventional thickness composites. Numerical results were compared to experimental data, demonstrating acceptable agreement. In the final step, following a sequential multi-scale modeling approach, a macro-scale model was constructed based on the outputs of the micro-scale model subjected to tensile and shear loads. The results were compared with experimental data, revealing good agreement. The proposed model allows for the optimization of pseudo-ductile composite structures to achieve a desired set of mechanical properties without the need for conducting extensive experimental material tests.

Keywords

Pseudo-ductile composite finite element model micro-mechanical model multiscale model

Introduction

In recent decades, unidirectional composites have been applied in different areas such as aerospace, civil engineering, and automotive industries.¹ However, these materials are usually brittle and may experience unexpected catastrophic failure which restricts their application.² Hybridization has been a solution to optimize the mechanical properties of composites while reducing the production costs.^3–5 This technique has recently been used to produce pseudo-ductile composites which have shown higher failure strain compared to conventional composites. Thus, a sudden collapse after overloading is not observed in these materials which minimizes the risks of the occurrence of a catastrophic failure.^6–8 The pseudo-ductility behavior in composites is obtained by hybridization of fibers with high and low failure strains, such as carbon and glass fibers, respectively.^9–12 Due to the difference in strain energy between plies, the middle ply fails first and then the created crack propagates through other plies until the whole structure fails. This gradual failure results in a nonlinear behavior similar to that in ductile materials. To observe the pseudo-ductility behavior, the thin-ply method can be employed in which a very thin high strain ply (i.e. the component with high failure strain) is used for construction of the composite.¹³ In the thin-ply method, the strength and failure strain of the pseudo-ductile composites can be improved by optimizing the sequence of fiber orientation of plies.¹⁴ Fabrication of such thin ply composites is significantly more time consuming and expensive compared to composites with the conventional thickness.¹⁵ A few research studies have been recently performed to produce pseudo-ductile composites with conventional thickness using carbon fibers with different elastic modulus and strain failure values.^7,15 For instance, Sapozhnikov et al.¹⁵ used pitch-based carbon fiber as the low failure strain, high elastic modulus fiber, and PAN-based carbon fiber, known for its high failure strain, as the high failure strain material. They manufactured various composite samples with different sequence layups, but a number of them exhibited a pseudo-ductile behavior.¹⁵

The fabrication methods for pseudo-ductile composites are typically more expensive than those used for conventional composites, such as hand lay-up or vacuum infusion plastic (VIP), due to the requirement for specialized equipment. Therefore, it would be very expensive to experimentally optimize the structure of pseudo-ductile composites for a specific application using trial and error techniques. Thus, developing a numerical or analytical model that can predict the mechanical properties of these composites based on the properties of their components would be beneficial to obtain an optimized design before fabrication. The analytical model proposed in Jalalvan et al. study,¹⁶ considers different failure possibilities of a pseudo-ductile composite and calculates the values of stress and strain in each scenario. Macro or micro-scale numerical models can also be developed to predict the mechanical behavior of pseudo-ductile composites. In a macro scale model, the exact geometry of each ply is considered and a cohesive-zone layer can be modeled between the plies to simulate the delamination damage.¹⁷ However, to accurately analyze the failure mechanisms of different composite components, a microscale model would be beneficial.

Despite the critical need for accurate models to predict the mechanical behavior of pseudo-ductile composites, current research is limited, and no numerical models exist for predicting their failure response under tensile loading. Pseudo-ductility, primarily influenced by the damage progression in different composite plies, has only been observed in composites with specific carbon fiber types and layup sequences.¹⁵ Determining the exact fiber and matrix characteristics that result in a pseudo-ductile response through trial and error experiments will be costly and inefficient. Therefore, the development of an experimentally validated computational framework that considers the micro-scale features of these composites will not only help in predicting accurate micro-mechanical failure, but will also aid in the optimization of the design of such composites.

This study aims to introduce a novel multi-scale FE model to predict damage progression in pseudo-ductile composites under tensile loads in the fiber direction, considering various material and geometric parameters. To this end, representative volume elements (RVEs) of three-ply all-carbon composite with varying fiber volume fractions in the low-strain ply were created at the micro-scale. Periodic boundary conditions were applied to the RVEs with periodic geometry and micro-mechanical simulations were performed on them. The accuracy of the micro-scale model was assessed by comparing the predicted tensile responses of both thin-ply and conventional thickness pseudo-ductile composites with those obtained from the experiments and an analytical model.¹⁶ In the next step, the elastic constants and strength parameters obtained from the micro-scale models were used in a macro-scale model to simulate the homogenized tensile response using the Hashin damage criteria.¹⁸ The tensile responses obtained from the initial micro-scale model and the subsequent multi-scale model were finally compared, and the advantages of each approach were discussed. The robust numerical tool developed in this study can be used for optimizing the mechanical behavior of pseudo-ductile composites, contributing to major advancements in the field of composite material engineering.

Pseudo-ductile composite

Pseudo-ductile composites include two different fibers with low and high failure strain as depicted in Figure 1. The high strain fibers used in the top and bottom layers (labeled as high-strain plies) have higher tensile strength and lower elastic modulus compared to the low strain fibers used in the middle ply of the composite (labeled as the low-strain ply). Gradual failure of the low-strain ply is the main reason for observing pseudo-ductile behavior in the composite. Two types of pseudo-ductile composite material systems exist, that is thin-ply composites and conventional thickness composites. To fabricate thin-ply composites, high pressure is applied during the manufacturing process, necessitating the use of an autoclave.⁸ However, for manufacturing conventional thickness pseudo ductile composites which have been recently introduced, employing a hot press would be sufficient.¹⁵ In this study, our focus is mainly on the latter type of composites due to their easier manufacturing process and also the mechanical disadvantages of the former type such as low impact resistance¹⁹ and notch tensile strength.²⁰ However, to present the capabilities of the proposed computational model, the model’s accuracy is assessed for both types of pseudo-ductile composites. Each type of pseudo-ductile composite employs specific low-strain and high-strain fibers and epoxy matrices tailored to achieve optimal pseudo-ductility, the details of which are described in the following sections of the paper.

Finite element model

Generation of the representative volume element (RVE)

To perform the micro-mechanical analysis of the three-ply hybrid composite shown in Figure 1, an RVE of the composite is constructed in which fibers in the micro-scale are distributed randomly through each ply with a specific fiber volume fraction. The fibers can be distributed using two different techniques. The first technique is the hardcore method, where fibers are distributed based on a random function. However, with this technique, generating an RVE with a high fiber volume fraction, such as 60%, is challenging.²¹ The second technique involves using the nearest neighbor algorithm. Following the initial random placement of the first fiber, subsequent fiber positions are determined with consideration for specific distance spacing and orientation relative to the first fiber.²² This technique allows for the generation of an RVE with a high fiber volume fraction. In this study, a Python code was developed to create the RVE geometry based on both fiber distribution techniques.

The input parameters to this code are the RVE dimensions including the thicknesses of high and low strain materials (2t_H and 2t_L) and width and depth of the RVE, fiber volume fraction of each ply (i.e. represented by V_f−L and V_f−H for low and high-strain materials, respectively), fiber radius (R_f), and a tolerance value. Fibers are distributed based on a random function, and in each iteration, the intersection between fibers is checked by ensuring that the minimum distance between the surface of each newly generated fiber and the previous fibers is greater than the tolerance value defined by the user. Values of the parameters inputted to the Python code are shown in Table 1.

Table 1.

Geometrical parameters of the model.

$R_{f} (m m)$	$V_{f - L}$	$V_{f - H}$	${2 t}_{H} = {2 t}_{L} (m m)$	width = depth (mm)
0.0065	0.1, 0.2, 0.3, 0.4	0.45	0.15	15-40 × $R_{f}$

In this study, the effect of V_f−L on the pseudo-ductile behavior of the composite under tensile load was investigated by considering four different fiber volume fractions in the middle ply namely, 10%, 20%, 30%, and 40% (Figure 2).

Material properties

In this study, the 736LT epoxy was considered for the matrix part of the composite and T800 and K63720 carbon fibers were considered as high and low failure strain materials, respectively. The mechanical properties of the carbon fibers and the epoxy used in the modelled composites are given in Table 2. According to Sapozhnikov et al. work,¹⁵ using certain volume fractions of these carbon fibers in a hybrid composite can result in a pseudo-ductile tensile response. In this study, the effect of volume fraction of carbon fibers on the mechanical response of composite is investigated to estimate the amount of carbon fibers required for obtaining a pseudo-ductile tensile response.

Table 2.

The fibers and matrix used for the simulated pseudo-ductile composite.^7,15

Property	T800	K63720	736LT epoxy
Elastic modulus (GPa)	294	620	3.15
Ultimate tensile strength (MPa)	5490	2650	101
Poisson’s ratio	0.22	0.25	0.35
Weibull modulus	4.8	5.61	-

On the fibers surface, some flaws exist which are made during the production process. Distribution of flaws on the surface of fibers can be different leading to a difference between the tensile strength of individual fibers. Thus, a suitable strength distribution function should be considered for the tensile strength of fibers. It has recently been shown that Weibull statistical function is suitable for distribution of the tensile strength of brittle fibers such as carbon and glass fibers.^23–25 The Weibull modulus is a fundamental parameter within the Weibull function that governs the range of strength distribution. A low Weibull modulus indicates a significant degree of variability in fiber strengths, thereby resulting in a wide spectrum of failure behaviors. Conversely, a high Weibull modulus signifies a more uniform distribution of fiber strengths, leading to more predictable patterns of failure. In accordance with Weibull statistics, the strength of each fiber can be determined as follows:

σ = σ_{0} {[\ln (\frac{1}{1 - η})]}^{\frac{1}{m}}

(1)

where η is a random number in [0,1], σ₀ and m are average ultimate tensile strength of fibers and Weibull modulus, respectively, and their values for the simulated RVEs are given in Table 2.

To model the failure behavior of each component of the composite, a user-defined material subroutine (UMAT) was developed in Fortran that can be linked with the ABAQUS software package²⁶ used for the FE analysis. In the UMAT code, the incremental stress and strain values are read from the ABAQUS FE solver and once the stress in each element exceeds the tensile strength value, the element is removed. The first step in the procedure involves updating the stress tensor, $σ_{i j}$ , in each element based on the current strain increment, ${∆ ε}_{i j}$ . The stress update follows the constitutive relations for the material, that is $σ_{i j} = C_{i j k l} {∆ ε}_{k l}$ , where $C_{i j k l}$ is the stiffness tensor of the material, and ${∆ ε}_{k l}$ represents the strain increment tensor.

To determine if an element should be deleted, the stress tensor components are compared with a predefined tensile strength limit, $σ_{t}$ . If any component of the stress tensor exceeds $σ_{t}$ , the element is considered to have failed.

A local damage variable, $d_{f}$ , is also introduced into the constitutive model which evolves based on the stress state of the element. Initially, $d_{f}$ is set to 0, indicating no damage. The damage variable is updated using the following rule:

d_{f} = {\begin{cases} 1 i f \max (σ_{i j}) \geq σ_{t} \\ 0 o t h e r w i s e \end{cases}

(2)

If the maximum stress exceeds $σ_{t}$ , $d_{f}$ is set to 1, signifying that the element has failed and should be deleted. This damage variable in incorporated into the stress tensor as below:

σ_{i j} = (1 - d_{f}) σ_{i j}

(3)

This micro-scale stress adjustment ensures that once the element in each fiber reaches the tensile strength limit that is d_f = 1, its contribution to the overall material stiffness and load-bearing capacity is nullified.

FE implementation

A composite structure consists of several RVEs that are in contact with each other, which means that the left, right, front, and back sides of an RVE are surrounded by other RVEs (see Figure 3). Thus, the displacement of nodes on each side of the RVE is controlled by that of the corresponding nodes on the parallel side. This can be simulated by applying periodic boundary conditions (PBCs) to the nodes on different sides of the RVE. It is important to note that the thickness of the RVE is the same as the thickness of the composite at the macro scale and thus PBCs are not assigned to the top and bottom sides of the RVE. To further clarify the concept behind a periodic RVE geometry, the two-dimensional surfaces of three similar RVEs stitched together are shown in Figure 4. It can be seen that the fibers on each edge of the middle RVE are complemented by those on the opposite edge and also all the fibers on the edges of the middle RVE form a complete fiber when RVEs are stitched together.

To apply the PBCs to the RVE, a displacement boundary condition, which specifies the displacement differences between each pair of nodes on parallel boundaries will be used as shown in equation (4):

U_{i}^{j +} - U_{i}^{j -} = {\bar{ϵ}}_{i k} [X_{k}^{j +} - X_{k}^{j -}] = {\bar{ϵ}}_{i k} Δ X_{k}^{j} (i, j, k = x, y, z)

(4)

where U_i is displacement in ith direction,

X_{k}^{j +}

is the coordinates of the node and ϵ is the macro strain tensor. j+ and j− represent the parallel boundary surfaces on opposite sides of an RVE of a composite. When an RVE is subjected to a specific strain of

{\bar{ϵ}}_{i k}

, the expression

X_{k}^{j +} - X_{k}^{j -}

remains constant for the parallel boundary surfaces, the righthand side of equation (4) will be constant. Such nodal displacement constraint is implemented in FE analysis to apply PBCs to the RVE.

Additionally, a displacement-controlled tensile load was applied to the RVE in the fiber direction. The RVE was then discretized into eight-node linear brick and reduced integration ABAQUS elements (C3D8R). In Section ’Mesh-sensitivity analysis’ the dependency of the results on the mesh size was evaluated through a mesh-sensitivity analysis and the element size length of 0.003 mm was selected to ensure a negligible difference between the results.

Convergence study

The tensile responses of several RVEs were simulated using the proposed numerical model to investigate the effect of RVE size, mesh size and the fiber volume fraction of the middle ply on the predicted results. results of these sensitivity analyses are discussed in the following sections.

Effect of the RVE size: To investigate the effect of RVE size on the results of the proposed numerical model, four RVEs with different width and depth sizes ranging from 15R × 5R to shows the cross section view of these RVEs which contain 40% fiber volume 40 × R were constructed.^27,28 Thickness of all RVEs is equal to 0.45 mm. Fraction in the middle ply. The failure strength and strain results obtained from the RVEs with different sizes (Figures 5 and 6) show that the RVE side length of 0.2 mm is large enough to have acceptable failure strength and strain. In the next step, this RVE size was considered to assess the mesh-sensitivity of the model.

Mesh-sensitivity analysis: To determine the proper element size for meshing the RVE, four RVEs with element sizes of 0.0045 mm to 0.003 mm were simulated by the proposed numerical model (Figure 7). While the proposed model is not quite mesh sensitive, the tensile failure strength and strain predicted by RVEs meshed with different element sizes (Figure 8) indicate that elements with side length of 0.003 mm can be considered for mesh generation of the composite models.

Results and discussion

As mentioned, two types of carbon fiber (i.e. T800 and K63720) with the 736LT epoxy resin were utilized to simulate a three-ply pseudo-ductile composite. The subsequent section presents the results of the finite element (FE) model, which will be compared with the results of the analytical model introduced in Jalalvan et al. work.¹⁶ Then, the accuracy of the model is assessed by simulating the tensile response of two distinct pseudo-ductile composites with thin-ply¹⁴ and conventional thickness.¹⁵ In the end, using a sequential multi-scale approach, the results of micro-mechanical simulations of RVEs subjected to tensile and shear loads were inputted to a macro-scale model to predict the tensile response of the conventional thickness pseudo-ductile composite based on the Hashin damage model.¹⁸

Failure process of the composite

In order to observe pseudo-ductile behavior in a composite structure, fiber fragmentation must first occur in the low-strain ply which is primarily due to the varying failure strain in each ply. When the low-strain ply is subjected to tensile load, due to its higher elastic modulus compared to the top and bottom high strain plies, it experiences greater stress concentration. As a result of the low failure strain of fibers in the low-strain ply, these fibers begin to fail individually, leading to a nonlinear pseudo-ductile stress-strain curve in the composite.

In this section, the stress distribution and progression of damage within an RVE of a pseudo-ductile composite with middle-ply volume fraction of 40% and side length of 0.2 mm are discussed. Figure 9 shows the longitudinal stress distribution within this RVE at different stages of the displacement-controlled tensile loading. At tensile strain (ϵ₂₂) of ϵ₂₂ = 0.17% (Figure 9(a)) due to the elastic modulus difference between the middle and top/bottom plies of the pseudo-ductile composite, stress concentration occurs in the fibers of this ply. Figure 9(b) then shows that at ϵ₂₂ = 0.58%, some fibers have failed sooner than the others, which is because of the tensile strength difference among fibers of the middle ply. Once the low-strain ply is completely failed, the stress shifts to the high-strain plies. At this stage of loading, where ϵ₂₂ = 1.33% (Figure 9(c)), fibers of the high-strain plies indicate a significant reinforcement role against external loads. Finally, the high-strain plies begin to fail when the stress in fibers reaches the tensile strength (Figure 9(d)).

Effect of fiber volume fraction of the middle ply: As explained in the previous section, the pseudo-ductile behavior is observed due to the occurrence of fiber fragmentation in the low-strain ply leading to a nonlinear stress-strain curve in the composite. Once the middle ply is completely failed, the stress is transferred to other plies, resulting in a linear increase in stress values.

In this section, to determine the range of the middle-ply volume fraction (V_f−L) that results in a pseudo-ductile stress-strain response, a number of RVEs with side length of 0.2 mm and different V_f−L were simulated using the proposed numerical model. The stress-strain curve obtained for the composite with V_f−L of 10% (Figure 10(a)) indicates that this composite does not show a pseudo-ductile behavior as all fibers of the middle ply fail almost at the same time. The stress-strain curves shown in Figure 10 indicate that increasing V_f−L from 10% to 30% led to the middle ply gradual failure and an increase in the first yield stress and initial elastic modulus of the composite (Figure 10(c) and (d)). As a result, the composites with V_f−L above 30% exhibited a pseudo-ductile tensile behavior as can be observed in their stress-strain curves. Since the composite with V_f−L = 40% exhibited a clear pseudo-ductile stress-strain curve, this composite was used for most numerical simulations of this paper.

Comparison between numerical and analytical results

In this section, the results obtained from the proposed numerical model are compared with those predicted by the analytical model explained in Appendix A.

The analytical model proposed four scenarios for failure of the composite, and the third and fourth scenarios can be recognized as pseudo-ductility behavior of the composite. This study will focus on modelling the failure behavior shown in the third scenario. The tensile stress-strain curves were obtained from the analytical model (Appendix A) for two three-ply composites with middle ply fiber volume fractions of 30% and 40%. The equations were proposed by Chamis²⁹ were employed to obtain the elastic modulus (E_c) and tensile strength (σ) of each ply, which calculates the properties of each ply based on the volume fraction of each component (Equations (5) and (6)).

E_{c} = E_{f} V_{f} + E_{m} (1 - V_{f})

(5)

σ_{c} = σ_{f} V_{f}

(6)

where, c, f, and m subscripts represent the composite ply, fiber, and matrix, respectively. The values of each parameter obtained using Chamis’s model are shown in Table 3.

Table 3.

Values of mechanical parameters for low-strain (L) and high-strain (H) plies used Appendix A and equation (7)

Parameter	30%	40%
E_L (MPa)	194205	257890
E_H (MPa)	134032	134032
S_L(MPa)	865.7	1120.60
S_H(MPa)	2526.1	2526.1
Kt ³⁰	2	2
G_II(N/mm)	0.267	0.405

To observe a pseudo-ductile behavior in a hybrid composite, it is essential to assess the following criterion as outlined in8.

G_{I I} = \frac{ϵ_{L b}^{2} E_{L} t_{L} (2 E_{H} t_{H} + E_{L} t_{L})}{8 E_{H} t_{H}} < G_{I I C}

(7)

In equation (7), G_II and G_IIC represent the energy release rates during the failure of the low-strain ply and the fracture toughness of mode II interface, respectively. ϵ_Lb denotes the failure strain of the low-strain ply. The energy release rate for each model is provided in Table 3.

According to Czél et al. study,⁸ the G_II value for a carbon/carbon composite should be less than 0.5 N/mm to observe a pseudo-ductile behavior. This criterion is met for both models with $V_{f - L}$ = 30% and $V_{f - L}$ = 40%.

Figure 11 compares the tensile stress-strain curves obtained from the numerical and analytical models and shows a good agreement between them. It is important to note that the analytical model considers a simple linear stress-strain response for each segment, whereas the numerical model can predict a nonlinear stress-strain behavior as it considers fiber fragmentation.

Validation of the model against experimental results

To validate the proposed model, two distinct types of pseudo-ductile composites were modeled, one manufactured using the conventional thickness methods¹⁵ and another using the thin-ply method.⁸ While the primary focus of this study is on the former type of composite, the latter is also modelled to assess the model’s capability in predicting the tensile response of other types of pseudo-ductile composites and also compare the failure process in these two types of composites. The conventional thickness composite labeled as T3DT3, consists of three plies with T800 carbon fibers located at the top and three at the bottom, serving as the high-strain plies, and also a single middle ply with K63720 carbon fibers as the low-strain ply. The matrix for fabricating this composite was ER450 epoxy.¹⁵ The thin-ply composite labeled as M55 in Czél et al. paper,⁸ is composed of two plies with T1000 carbon fibers at the top and two at the bottom, serving as the high-strain plies, and one middle-ply with M55 carbon fibers as the low-strain ply. The matrix used for fabricating these composite samples is EPHTg-402 epoxy.⁸ The properties of the carbon fibers and matrices used in these two types of composites are provided in Table 4.

Table 4.

Properties of the fibers and epoxy used for fabrications of the conventional thickness¹⁵ and thin-ply⁸ pseudo-ductile composites.

Composite type	Fiber	Elastic modulus (GPa)	Strain failure (%)	Strength (MPa)	Diameter (µm)
Thin ply
	T1000⁸	294	2.2	6370	5
	M55⁸	540	0.8	4020	5
	ER450 epoxy	3.56	3.09	110	-
Conventional thickness
	T800¹⁵	294	2.0	5490	5
	K63720¹⁵	640	0.4	2650	11
	EPHTg-402 epoxy	3.44	2.57	62	-

Since the middle ply gradual failure and the inter-ply delamination are the primary reasons for the nonlinear behavior exhibited by the conventional thickness and thin-ply pseudo-ductile composites, in this section, in order to save the computational costs while maintaining a high level of accuracy, only the middle layer is modelled with detailed fiber characteristics (micro-scale model), and the top and bottom layers are modeled as a macro-scale model. The mechanical properties of the top and bottom high-strain plies are obtained based on a homogenization approach using the Chamis’s formulas (Equations (5) and (6)).²⁹ The fiber volume fraction in the middle ply of the conventional thickness and thin-ply composites are 55% and 51.3%, respectively. The nearest neighbor algorithm²² employed to generate fibers within the matrix (Section FE implementation) ensured that each successive fiber was positioned in close proximity to the initial one and thus such high fiber volume fractions (above 45%) are feasible. The Weibull modulus (m) for all fibers was assumed to be 6.6.

In the failure process of a thin-ply pseudo-ductile composite, in addition to the low-strain ply gradual failure, the inter-ply delamination significantly contributes to the nonlinear pseudo-ductile response. Therefore, a cohesive zone was incorporated between the low-strain ply and high-strain plies in the micro-scale. The parameters utilized for modeling the cohesive zone, such as the cohesive strength, energy dissipation, and cohesive stiffness, are provided in Table 5.³¹ Figures 12 and 13 show the RVE of the conventional thickness¹⁵ and thin-ply⁸ pseudo ductile composites, respectively.

Table 5.

Properties of cohesive zone between high and low-strain plies.³¹

$G_{I c} (N / m m)$	$G_{I I c} = G_{I I I c} (N / m m)$	σ_n⁰ = τ_tn⁰ = τ_sn⁰ (MPa)	K (MPa/mm)
0.199	0.5	67	10⁵

After generating these RVE, the load and PBCs mentioned in Section FE implementation were applied to the meshed RVEs. The tensile behavior of these RVEs was then analyzed using the ABAQUS FE software linked with our UMAT Fortran subroutine.

Figure 14 shows the tensile stress-strain curves obtained from the proposed model and experimental results in Sapozhnikov et al.¹⁵ and Czél G et al.⁸ studies. Figure 14 shows that the proposed numerical model accurately predicts the tensile response of pseudo-ductile composites, irrespective of their manufacturing methods, that is thin-ply and conventional thickness. Figure 14 also shows that by incorporating the cohesive zone to the model, it can effectively capture delamination between low-strain and high-strain plies in the thin-ply composite, resulting in improved alignment between numerical and experimental results.

Figure 15 shows the von Mises stress distribution in both types of pseudo-ductile composites at different strain values during the displacement-controlled tensile loading. In Figure 15(c) the delamination and fiber failure are both the reason for gradual failure in thin-ply pseudo-ductile composite. However, fiber failure of the V_f−L is the main reason for gradual failure of conventional thickness composite (15-f).

When the RVE of the thin-ply pseudo-ductile composite is subjected to tensile loading, fibers within the middle ply experience elevated stress concentrations, leading to their initial failure. As the fibers in the low-strain ply fail subsequently, stress transmission occurs among the various plies, resulting in heightened stress concentration at the ply interfaces. This heightened stress level augments the interlaminar energy between plies, ultimately leading to delamination. The delamination between low and high-strain plies of thin-ply pseudo-ductile composite is shown in Figure 16. In Figure 16, ‘COPEN’ in the clearance between two contacted surfaces in mm.

Effect of the middle-ply fibers’ tensile strength variation: It was shown in the previous sections that fragmentation of the low-strain material is the primary cause for observing pseudo-ductile behavior in pseudo-ductile composite with conventional thickness. In addition, variation of the tensile strength in the fibers of the middle ply has a critical role in the low-strain ply fragmentation. In the proposed model, the Weibull modulus is capable of altering the strength distribution among fibers of the middle ply, resulting in a change in the pseudo-ductile behavior of the composite. Therefore, in this section, five different Weibull moduli were considered for the low-strain material in the T3DT3 composite (Figure 17), and the tensile response predicted by the proposed model are discussed. Variations of the fiber strength for RVEs modeled using different Weibull moduli are presented in Figure 17 and the stress-strain curves predicted by the proposed model and experiments are shown in Figure 18.

For RVEs with a high scatter in the tensile strength values of the middle-ply, that is a Weibull modulus of 1.6 in our model (Figure 17(a)), each middle-ply fiber fails at very different strain values resulting in a linear stress-strain response and not a pseudo-ductile behavior (Figure 18(a)). However, as this scatter is reduced meaning that there is less difference between the tensile strength of middle-ply fibers, that is Weibull moduli of 3.6 and 6.6 in our model (Figure 17(b) and (c)), the fibers fail at closer strain values. This results in an ideal pseudo-ductile behavior in the tensile stress-strain response of the composite (Figure 18(b) and (c)). In the case that the strength values of the middle-ply fibers are very close to each other, that is Weibull moduli of 10.6 and 15.6 in our model (Figure 17(d) and (e)), a large number of fibers will fail at almost the same strain value, leading to a sudden drop in stress values, as observed in Figure 18(d) and (e). Although in these cases (Figure 17(d) and (e)), the pseudo-ductile behavior is still observed, the sudden drop in the horizontal load-sharing plateau of the tensile stress-strain curve might not be ideal for certain structural and mechanical engineering load-bearing applications.

The analysis performed in this section indicates that in order to observe an ideal pseudo-ductile behavior, variation of the fibers’ tensile strength in the middle ply should be within a certain range. Because this variation is mainly controlled by the type of fibers used in the composite and the manufacturing process, it is important that it is optimized prior to the fabrication. The proposed micro-mechanical model can therefore be considered as a valuable tool as it is capable of controlling the fibers’ tensile strength variation by selecting the appropriate Weibull modulus.

Incorporation of micro-mechanical results into a macro-scale model: In this section, elastic constants of a T3DT3 composite (with fiber volume fraction of 55% in each layer) were calculated using a multi-scale modelling technique. In this approach, micro-scale RVEs were first created for each of the low-strain (‘D’) and high strain (‘T’) plies with periodic geometries (Figure 19). The elastic constants used for modelling the carbon fibers are shown in Table 6 and those used for matrix are EPHTg-402 epoxy. The RVEs were discretized using ABAQUS C3D8R elements (element size = 0.0025 mm for D ply and element size = 0.001 mm for T ply). PBCs were employed to the boundaries of each RVE, and tensile and shear loads were applied on different planes and along different directions to determine the elastic constants of the composite.

Table 6.

Transverse isotropic properties of T800 and K63720 carbon fibers.

Fiber	E₁₁(GPa)	E₂₂(GPa)	E₃₃(GPa)	ν12	ν13	ν23	G₁₂(GPa)	G₁₃(GPa)	G₂₃(GPa)
T800	294	30	30	0.22	0.22	0.05	27.3	27.3	3
K63720	640	40	40	0.22	0.22	0.05	37.3	37.3	4

The procedure for obtaining the failure properties of each ply from the micro-scale RVEs is similar to that outlined in Section Material properties. Upon applying tensile and shear loading conditions to the RVEs, the longitudinal and transverse strengths of the composite are calculated for use in the macro-scale model. Figure 20 shows the deformed geometry of micro-scale RVEs under tensile and shear loads.

After obtaining the elastic and failure properties of ‘D’ (low strain) and ‘T’ (high strain) plies from the micro-scale models (RVE Dl, RVE Dt, RVE Ds, RVE Tl, RVE Tt, and RVE Ts shown in Figure 20), the fracture energy of each RVE was determined by calculating the area under the predicted stress-strain curve until failure. Table 7 shows the values of elastic constants and fracture energies along different directions for ‘D’ and ‘T’ plies.

Table 7.

Outputs of the meso models.

Ply	E₁₁(MPa)	E₂₂(MPa)	ν12	G₁₂(MPa)	G₁₃(MPa)	G₂₃(MPa)
T	161356	8919	0.14	4363	4428	2010
D	320972	9339	0.14	4383	4449	2012
	X_T (MPa)	X_C (MPa)	Y_T (MPa)	Y_C (MPa)	S_L (MPa)	S_T (MPa)
T	2311	2311	55.10	55.10	49.91	36.43
D	1370	1370	52.04	52.04	44.61	34.74
	G_LT (N/mm)	G_LC (N/mm)	G_TT (N/mm)	G_TC (N/mm)
D	1.8	1.8	0.51	0.51

The ’D’ ply significantly contributes to the pseudo-ductile behavior of the composite. To model the damage evolution of the ’D’ ply, the fracture energy of this ply at different directions was utilized as input data for the macro model.

For the macro-scale model, a rectangular shell geometry measuring 150 mm × 15 mm was generated in ABAQUS FE software by considering the appropriate layer sequences (i.e. T3/D/T3). The Hashin damage criterion¹⁸ was employed to simulate the damage of each ply. The Hashin damage initiation criterion considers four potential failure modes, each associated with a specific index: (i) fiber breakage in tension ( $F_{f}^{t}$ ), (ii) fiber buckling in compression ( $F_{f}^{c}$ ), (iii) matrix cracking in tension ( $F_{m}^{t}$ ), and (iv) matrix crushing in compression ( $F_{m}^{c}$ ). The corresponding equations for these failure modes are as follows:

F_{f}^{t} = {(\frac{σ_{11}}{X_{T}})}^{2} + α {(\frac{σ_{12}}{S_{L}})}^{2} \leq 1.0 a n d σ_{11} \geq 0

(8)

F_{f}^{c} = {(\frac{σ_{11}}{X_{c}})}^{2} a n d σ_{11} \leq 0

(9)

F_{m}^{t} = {(\frac{σ_{22}}{Y_{T}})}^{2} + {(\frac{σ_{12}}{S_{L}})}^{2} \leq 1.0 a n d σ_{22} \geq 0

(10)

F_{m}^{c} = {(\frac{σ_{22}}{{2 S}_{T}})}^{2} + [{(\frac{Y_{C}}{{2 S}_{T}})}^{2} - 1] \frac{σ_{22}}{Y_{C}} + {(\frac{σ_{12}}{S_{L}})}^{2} \leq 1.0 a n d σ_{22} < 0

(11)

In the given equations, X_T, Y_T, X_C, and Y_C denote the tensile and compressive strength along the longitudinal and transverse directions, respectively. The values of G_LT, G_LC, G_TT, and G_TC from Table 7 represent the tensile and compressive fracture energies along the longitudinal and transverse directions.

The tensile response of the macro-scale model is illustrated in the stress-strain curve in Figure 21. In Figure 21, the progressive damage of the D ply at different points is depicted on the stress-strain diagram. As previously mentioned, the gradual failure of the low-strain ply (D ply) is the key factor contributing to the pseudo-ductile behavior of the conventional thickness composites, a phenomenon clearly evident in Figure 21. In the initial segment of the curve, all fibers within the D ply remain undamaged (point 1). As strain values increase, fiber damage initiates (point 2). Further strain elevation leads to the sequential failure of additional fibers, gradually encompassing the entire model in the fiber damage contour (points 3 and 4). Following the complete failure of fibers in the D layer, the T plies start to carry the tensile load (point 5).

Furthermore, this phenomenon was depicted in the micro-scale model as well. For a comparative analysis of the gradual failure of the D ply in both the micro and macro-scale models, the damage progression of this ply under different strain values is shown in Figure 22. In this figure, HSNMTCRT and HSNFTCRT represent the damage values of the matrix and fiber, respectively. These parameters range between 0 and 1, where a value of 0 signifies no damage and one indicates complete damage in the component.

At low tensile strain values, both fibers and matrix experienced tensile stress with no damage occurring in the fibers (Figure 22(a)). As the tensile strain increases, gradual fiber failure begins in the ’D’ ply, resulting in a gradual decrease in the tensile stress values and the observation of a plateau in the tensile stress-strain diagram (Figure 22(b) to (d)). This gradual failure persists until all fibers in the ’D’ ply fail, while the matrix remains undamaged (Figure 22(e)). Upon complete fiber failure in this ply, stress transfers to other plies while the matrix in the ’D’ ply remains intact (Figure 22(f)).

Comparison of the stress-strain curves predicted from the micro-scale and the macro-scale model generated with micro-mechanical model inputs shows that the error values in the nonlinear segment are more pronounced in the macro model, as depicted in Figure 23. Since the micro-scale model accurately considers the components of each layer and also assumes different tensile strengths for each fiber through the Weibull distribution, it can accurately simulate the damage progression in each ply. In contrast, the macro-scale model relies on equations and fracture energy for calculating the damage of each component. Additionally, inherent errors in the FE method contribute to modeling inaccuracies, and these errors may be exacerbated when analyzing at two scales. This indicates a micro-level error during the calculation of elastic constants or damage parameters, which is then compounded with macro-level errors, resulting in a significant overall discrepancy. Using a macro-scale model developed by a sequential multi-scale approach is advantageous because it allows us to test the composite under different loading conditions, such as bending or impact, with lower computational costs compared to the micro-scale model.

Conclusion

This paper presents a numerical multi-scale finite element model (micro and macroscales) for predicting the damage behavior of different types of pseudo-ductile composites, with a focus on those with conventional thickness. Different tensile strengths, based on the Weibull distribution function, were assigned to fibers in each ply. The proposed model was validated by comparing the numerical results with those obtained from recent experiments conducted on thin-ply and conventional thickness pseudo-ductile composites. The numerical results were also compared with those predicted by an analytical model¹⁶ described in the Appendix. Different sensitivity analyses were performed and the effect of the volume fraction and also the tensile strength variation of the fibers in the low-strain ply (D-ply), as the most effective parameters, were investigated. Based on the outcomes of this research, the following conclusions can be drawn:

• Pseudo-ductile behavior can be observed in composites with 30% and 40% fiber volume fractions of the middle ply, while composites with low fiber volume fractions (10% and 20%) showed linear behavior leading to failure.

• Increasing the fiber volume fraction of the middle ply leads to a higher initial elastic modulus and pseudo-ductile behavior.

• The comparison between numerical and analytical results shows a good agreement in predicting the tensile failure of the composite. For composites with low fiber volume fractions (10% and 20%), the division is increased after the middle ply fracture due to low fiber amounts.

• The middle-ply fibers’ tensile strength variation plays a critical load in the tensile pseudo-ductile behavior of composites: The composites with a large scatter of the fibers’ tensile strength in the middle ply (low Weibull modulus in our model) do not exhibit a pseudo-ductile behavior. Reducing this scatter (increasing the Weibull modulus in our model) results in a nonlinear and load-sharing horizontal plateau in the tensile stress-strain diagram of the composite.

In the end, using a sequential multi-sale modeling approach the results of the micro-mechanical simulations were inputted to a macro-scale model to predict the pseudo-ductile behavior of the conventional thickness composite. Such a multi-scale FE model offers acceptable accuracy while minimizing the computational costs. This approach will be particularly of interest if dynamic loads such as impact are to be considered where the three-dimensional state of stress needs to be accurately modeled.

The outcome of this study provides valuable insights into the design and engineering of pseudo-ductile composites. The proposed model enables us to design an optimized composite structure tailored to the final application prior to fabrication, thereby significantly reducing production costs by minimizing the requirement for extensive experimental tests. Future work could explore the mechanical response of these composites under other loading scenarios such as impact and fatigue loads which are highly experienced by them in different engineering applications.

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

Behzad Abdellahi

Appendix

References

Netravali

Chabba

. Composites get greener. Mater Today 2003; 4(6): 22–29.

Richardson

Wisheart

. Review of low-velocity impact properties of composite materials. Compos Appl Sci Manuf 1996; 27(12): 1123–1131.

Atiqah

Maleque

Jawaid

, et al. Development of kenaf-glass reinforced unsaturated polyester hybrid composite for structural applications. Compos B Eng 2014; 56: 68–73.

Fukuda

Chou

. Stress concentrations in a hybrid composite sheet, 1983.

Manders

Bader

. The strength of hybrid glass/carbon fibre composites. J Mater Sci 1981; 16(8): 2246–2256.

Sapozhnikov

Kheruvimov

Bezmel’nitsyn

. Design of pseudoductile unidirectional hybrid composites based on the model of a dry bundle. Mech Compos Mater 2019; 55(2): 231–240.

Sapozhnikov

Lomov

Swolfs

, et al. Deformation and failure of pseudo-ductile quasiisotropic all-carbon hybrid frps with an open hole under tension. Compos B Eng 2022; 237: 109870.

Czél

Jalalvand

Wisnom

, et al. Design and characterisation of high performance,’ pseudo-ductile all-carbon/epoxy unidirectional hybrid composites. Compos B Eng 2017; 111: 348–356.

Wisnom

Mechanisms to create high performance pseudo-ductile composites. IOP Conference Series: Materials Science and Engineering 2016. Vol. 139, No. 1. IOP Publishing.

10.

Ribeiro

Sena-Cruz

Branco

, et al. Hybrid effect and pseudo-ductile behaviour of unidirectional interlayer hybrid frp composites for civil engineering applications. Construct Build Mater 2018; 171: 871–890.

11.

Suwarta

Fotouhi

Czel

, et al. Fatigue behaviour of pseudo-ductile unidirectional thin-’ ply carbon/epoxy-glass/epoxy hybrid composites. Compos Struct 2019; 224: 110996.

12.

Leshkov

Sapozhnikov

. Modeling the nonlinear deformation and damage of carbonaramid fabric composites in tension. Mech Compos Mater 2020; 56(5): 591–600.

13.

Czél

Wisnom

. Demonstration of pseudo-ductility in high performance glass/epoxy’ composites by hybridisation with thin-ply carbon prepreg. Compos Appl Sci Manuf 2013; 52: 23–30.

14.

Czél

Jalalvand

Wisnom

. ‘Design and characterisation of advanced pseudoductile unidirectional thin-ply carbon/epoxy–glass/epoxy hybrid composites. Compos Struct 2016; 143: 362–370.

15.

Sapozhnikov

Swolfs

Lomov

. Pseudo-ductile unidirectional high modulus/high strength carbon fibre hybrids using conventional ply thickness prepregs. Compos B Eng 2020; 198: 108213.

16.

Jalalvand

Czél

Wisnom

. Damage analysis of pseudo-ductile thin-ply ud hybrid’ composites–a new analytical method. Compos Appl Sci Manuf 2015; 69: 83–93.

17.

Jalalvand

Czél

Wisnom

. Numerical modelling of the damage modes in ud thin’ carbon/glass hybrid laminates. Compos Sci Technol 2014; 94: 39–47.

18.

Hashin

. Fatigue failure criteria for unidirectional fiber composites. J Appl Mech 1980; 47(4): 329–334.

19.

Sasikumar

Trias

Costa

, et al. Mitigating the weak impact response of thin-ply based thin laminates through an unsymmetrical laminate design incorporating intermediate grade plies. Compos Struct 2019; 220: 93–104.

20.

Furtado

Arteiro

Catalanotti

, et al. Selective ply-level hybridisation for improved notched response of composite laminates. Compos Struct 2016; 145: 1–14.

21.

Ohser

Mucklich

. Statistical analysis of microstructures in materials science.¨ Pract Metallogr 2001; 38(9): 538–539.

22.

Vaughan

McCarthy

. A combined experimental–numerical approach for generating statistically equivalent fibre distributions for high strength laminated composite materials. Compos Sci Technol 2010; 70(2): 291–297.

23.

Bunsell

Gorbatikh

Morton

, et al. Benchmarking of strength models for unidirectional composites under longitudinal tension. Compos Appl Sci Manuf 2018; 111: 138–150.

24.

Mesquita

Bucknell

Leray

, et al. Single carbon and glass fibre properties characterised using large data sets obtained through automated single fibre tensile testing. Compos Appl Sci Manuf 2021; 145: 106389.

25.

Thomason

. On the application of weibull analysis to experimentally determined single fibre strength distributions. Compos Sci Technol 2013; 77: 74–80.

26.

Dassault Systemes Simulia Corp . `ABAQUS/Standard user’s manual, 2016.

27.

Tavares

Melro

Bessa

, et al. Mechanics of hybrid polymer composites: analytical and computational study. Comput Mech 2016; 57(3): 405–421.

28.

Swolfs

Gorbatikh

Romanov

, et al. Stress concentrations in an impregnated fibre bundle with random fibre packing. Compos Sci Technol 2013; 74: 113–120.

29.

Chamis

. Simplified composite micromechanics equations for strength, fracture toughness, impact resistance and environmental effects. In: Lewis research center Cleveland, NASA TM83696, 1984.

30.

Watanabe

Tanaka

Higuchi

, et al. A study of stress concentrations around fiber breaks in unidirectional cf/epoxy composites using double-fiber fragmentation tests. Adv Compos Mater 2018; 27(6): 575–587.

31.

Ollic

Sheibanian

Fazlali

, et al. Modelling the damage evolution in unidirectional all-carbon hybrid laminates. Composites Science and Technology 2024; 1247: 110420.