Low-Velocity Impact Response of Hybrid Fiber Reinforced Composite Thin-Walled Structures

1  Introduction

Hybrid fiber reinforced plastic (HFRP) composites have emerged as pivotal structural materials in aerospace, automotive, and marine engineering, owing to their tailorable mechanical properties, exceptional weight-reduction potential, and enhanced cost-effectiveness compared to single-fiber systems [1–3]. Among various hybridization strategies, intra-layer hybrids combining carbon fibers (CFs) and glass fibers (GFs) exhibit distinctive synergistic effects. Specifically, CFs deliver high specific strength and stiffness [4–6], while GFs contribute superior impact resistance, corrosion tolerance, and economic viability [7,8]. This intra-layer hybrid configuration effectively mitigates the inherent brittleness of CFs while overcoming the low-modulus limitation of GFs, offering a promising solution for load-bearing structures under complex service conditions.

Nevertheless, the mechanical behavior of intra-layer carbon/glass hybrid fiber composites remains insufficiently characterized, primarily owing to the complex multiscale interactions between the two fiber types and the polymer matrix. The stress distribution, fiber-matrix interfacial behavior, and failure mechanisms (e.g., fiber fracture, matrix cracking, and delamination) of intra-layer hybrids differ significantly from those of inter-layer configurations, as the intimate mixing of fibers introduces additional microstructural heterogeneity [9,10]. As a common occurrence, low-velocity impact (LVI) constitutes a prevalent threat during manufacturing, transportation, and in-service operations, often inducing barely visible internal damage that severely impairs residual strength and long-term durability [11–13]. For thin-walled structures like hat-shaped beams, LVI-induced damage can critically undermine safety, underscoring the necessity for systematic investigation.

Regarding LVI performance, existing research has predominantly focused on inter-layer hybrids or single-fiber composites. Maier and Mandoc [14] observed that incorporating glass fiber layers improves the impact energy absorption of carbon-fiber composites, while Ahamed et al. [15] demonstrated that impact damage area depends on the thickness of glass fiber layers. However, the LVI behavior of intra-layer carbon/glass hybrid composite hat-shaped beams has scarcely been addressed. The complex geometry of such beams induces highly non-uniform stress distributions under impact, and intra-layer fiber heterogeneity further complicates damage initiation and propagation. Moreover, the multitude of design variables and high costs associated with experimental trial-and-error have impeded the widespread industrial adoption of these materials, representing a significant research gap.

Advances in computational materials science have established numerical simulation as a powerful tool to complement experimental studies. While macromechanical models have been employed to simulate inter-layer hybrids [16–18] with reasonable computational efficiency, they typically neglect the influence of key microstructural parameters on structural performance. Consequently, such models cannot adequately characterize the multi-scale features of HFRP structures, limiting the full exploitation of their mechanical and economic potential. To overcome these limitations, multi-scale modeling approaches have been developed, which integrate heterogeneous microstructural information into macroscopic constitutive models via localization and homogenization techniques. Although various multiscale frameworks have been proposed [19–21], most are limited to predicting elastic properties and fail to capture progressive damage and failure, thereby decoupling macroscopic material failure from its microscopic origins. For energy-absorbing component design. In this context, Ha et al. [22–24] presented a multi-scale framework derived from the micromechanics of failure (MMF) theory. By incorporating stress amplification factors (SAFs), this method establishes a direct linkage between micro- and macro-scale stresses, enabling the determination of macroscopic failure criteria and damage parameters from micromechanical analysis [25–27]. This framework provides a promising foundation for the multiscale design of HFRP structures.

On the whole, HFRP composites possess great potential regarding comprehensive performance and economic efficiency, but their energy-absorption mechanisms under LVI remain unclear, and the effects of critical parameters such as hybrid ratio, hybridization scheme, and weaving pattern require systematic exploration. Furthermore, an efficient, physics-based design methodology for HFRP structures is still unavailable. This work intends to address these gaps by developing a multi-scale modeling and integrated design framework for intra-layer carbon/glass hybrid composites from micro-meso-macro scales. In detail, a comprehensive experimental program involving tensile and LVI tests was conducted on hybrid composite panels and hat-shaped beams with different configurations to quantify the effects of hybridization schemes and ratios, with tensile data serving primarily to validate multiscale model inputs. Subsequently, a new multi-scale model grounded in micromechanical failure theory was implemented via a VUMAT user subroutine in Abaqus. The validated model was employed for comparative numerical analysis. Finally, an orthogonal discrete optimization is conducted with the carbon fiber bundle fraction at the mesoscopic and the ply angle at the macroscopic as the design variables.

2  Experimental Study

In this part, various laminates and hatshaped beams were manufactured and tested under tension and impact loading. According to the measured results, the mechanical properties and influencing factors of HFRP specimens were discussed.

2.1 Preparation of Test Specimens

This section elaborates on the manufacturing procedures of CFRP, GFRP (Glass fiber reinforced plastic), and HFRP composite laminates and beam specimens. First, carbon fiber and glass fiber bundles were woven into various plain-woven dry fabrics using a semi-automatic small-sample weaving machine, as shown in Fig. 1a. These fabrics included pure carbon, pure glass, and intra-layer hybrid dry fabrics with different hybrid ratios. Specifically, the carbon fiber content is 75% for C0/CG90, 50% for C0/G90 and CG, and 25% for CG0/G90. The center-to-center spacing between two adjacent fiber bundles is uniformly 1 mm for all hybrid configurations, detailed geometry is presented in Fig. 1b. Thereafter, the laminates were manufactured via the thermoforming technique. In detail, various types of fiber fabrics were cut into suitable sizes and then impregnated with epoxy resin. The resulting preforms were placed into a metal mold for thermal curing, conducted at 120°C under a pressure of 0.5 MPa; after molding. According to the composite testing standard ASTM D3039, the laminate is cut to dimensions of 250 mm × 25 mm × thickness. Table 1 presents the geometric dimensions of the composite laminates, with three replicate tests performed for each specimen configuration. As specified in Table 1, the 6-layer CFRP specimen is designated as “P-C6”, and the corresponding 6-layer GFRP specimen is denoted as “P-G6”. The letter “P” represents a composite plate. Besides, the intra-layer hybrid specimens were named “P-[C0/CG90]6”, “P-[C0/G90]6”, “P-[CG0/G90]6” and “P-CG6”, respectively. Where “0” denotes fibers oriented in the 0° direction and “90” denotes fibers oriented in the 90° direction. The typical composite plates are shown in Fig. 2.

Figure 1: Fabrication and geometric dimensions of the fiber fabric: (a) weaving pattern; (b) geometric dimensions.

Figure 2: Specimens for tensile tests.

In order to evaluate the impact resistance of HFRP structures, a set of single hat-shaped beams was manufactured through the thermoforming technique. Intra-layer hybrid HFRP beams were compared with single fiber composite materials to reveal their superiorities in energy absorption characteristics. The hot compression molding process of the beams is schematically illustrated in Fig. 3. After molding, the preforms were trimmed into specimens for impact testing, all beams had a length of 100 mm and a width of 70 mm. To prevent bottom-initiated collapse of the specimens during impact, a 45° chamfer was machined on their top edges. Values for thickness and mass are presented in Table 2, where the letter “B” denotes the composite beams. Each specimen was tested in duplicate to ensure result reliability.

Figure 3: Schematic diagram of hat-shaped beam fabrication.

2.2 Tensile Tests

To examine the tensile performance of HFRP laminates, a set of tensile tests was carried out, demonstrating their superior mechanical properties compared to single fiber composite laminates. The influence of the hybridization scheme on the tensile behavior is discussed in this section.

2.2.1 Experimental Set-Up

Tensile tests were performed using a MTS universal testing system in accordance with the ASTM D3039 standard. The tests were performed at a loading rate of 2 mm/min. For each mixing method, tensile tests were carried out on specimens oriented at 0° and 90°, with three replicate tests conducted for each group. To prevent damage to the composite specimens caused by the grips during testing, reinforcement tabs were bonded to both ends of the specimens. Additionally, the specimens were sprayed with a uniform black-and-white speckle pattern. During the tensile tests, a camera was used to measure the full-field deformation response of the specimens. By combining with the VIC-2D software, the true strain values of the specimens during the tests could be obtained. The force values collected by the universal testing machine correspond to engineering stress, and the true stress was calculated using the following equations:

σt=σneεt(1)

εt=ln⁡(1+εn)(2)

where σt denotes the true stress; σn represents the engineering stress measured by the testing machine; and εt and εn correspond to the true strain and engineering strain, which are obtained from the VIC-2D software.

2.2.2 Experimental Results and Discussion

Fig. 4 presents the stress-strain curves of all laminates. The elastic modulus and strength values of all specimens were calculated from these stress-strain curves, with the results summarized in Table 3. As observed from Fig. 4 and Table 3, a comparative analysis reveals slight variations among repeated tests. These variations may stem from incomplete resin flow during the resin brushing process; in addition, they could be caused by insufficient clamping force of the fixtures during testing or slippage of the reinforcement sheets relative to the specimens. Nevertheless, the errors in the test results remain within the acceptable range.

Figure 4: Stress-strain curves of tensile specimens: (a) The curves under tension along the 0° direction; (b) The curves under tension along the 0° direction.

Aligning carbon fibers along the tensile direction results in relatively high elastic modulus and strength for the P-[C0/G90]6 and P-[C0/CG90]6 configurations under 0° tension, highlighting the significance of hybrid designs with carbon fibers in the load-bearing direction. Conversely, when glass fibers are aligned along the tensile direction, the P-[C0/G90]6 and P-[CG0/G90]6 configurations exhibit inferior mechanical performance under 90° tension, with modulus and strength approaching those of GFRP. Notably, increasing carbon fiber content enhances mechanical properties in the load-bearing direction, while higher glass fiber content slightly improves tensile strength in the non-load-bearing direction. This behavior is attributed to the fiber hybridization effect, which elevates inter-bundle stress and induces stress concentration. Additionally, the superior tensile properties of carbon fibers and greater plasticity of glass fibers may further contribute to the observed mechanical response.

2.3 Impact Tests

For the purpose of exploring the impact resistance characteristics and failure mechanism of hat-shaped beams, a series of impact tests were conducted. This section will discuss the effect of different intra-layer hybrid weaving patterns on their impact behavior.

2.3.1 Key Performance Parameters

To quantitatively characterize the energy absorption performance of the tested impact specimens, multiple evaluation indicators were adopted in this study. As a core requirement for energy-absorbing components in engineering applications, a well-designed structure should effectively protect occupants during collisions by dissipating more kinetic energy [3,4]. Therefore, Energy Absorption (EA) is widely recognized as a fundamental indicator to assess the energy absorption capability of components, which is defined as:

EA=∫0hF(x)dx(3)

where F(x) corresponds to the crushing force at displacement x; and h stands for the overall crushing displacement.

In addition to EA, Peak Force (Fmax) was introduced as a key force-related indicator to evaluate the impact force characteristics. It corresponds to the maximum reaction force induced by the structure throughout the collision process.

Additionally, the specific energy absorption (SEA) is commonly employed to evaluate the energy-absorbing efficiency of components and is defined as follows:

SEA=EAm(4)

where m denotes the mass of the specimen.

Finally, to comprehensively assess the practical application value of the hat-shaped beam, the cost-efficiency indices, namely EA/Cost (energy absorption per unit cost), was also calculated and compared in this study.

2.3.2 Experimental Set-Up

Impact tests were performed using a drop-weight apparatus (INSTRON 9450), as shown in Fig. 5. For the tests, the impact system was configured to match the target impact energy of 200 J. The impact hammer was set with a mass of 20.1 kg, and the impact velocity was adjusted to 4.46 m/s. The hat-shaped structure was placed on the rigid flat plate of the testing machine’s workbench, where its bottom edge was left unconstrained to minimize the influence of boundary conditions on test results. A high-speed camera was used to record the deformation process of the specimens, and the force sensor equipped on the drop weight testing machine captured the contact reaction force during the impact.

Figure 5: Experimental setup of impact test.

2.3.3 Experimental Results and Discussion

The load-displacement curves from the axial low-speed impact tests of hat-shaped beams are presented in Fig. 6. Repeated tests exhibit high consistency, strongly confirming test reproducibility. Three distinct stages are identified via curve analysis: in the initial stage, load rises sharply with increasing displacement, reflecting the beam’s elastic deformation in the early impact phase and rapid enhancement of bearing capacity. Among all specimens, B-G6 exhibits the highest Fmax; B-[CG0/G90]6, B-[C0/G90]6, and B-CG6 follow with relatively close values; B-C6 and B-[C0/CG90]6 record the lowest. Overall, Fmax increases with decreasing carbon fiber content, rising from glass fibers’ higher ductility, which enhances structural resistance to initial deformation. Upon reaching this peak, each curve enters a descending stage, corresponding to the initiation and accumulation of plastic deformation, internal material damage, and even local failure, with a consequent decrease in bearing capacity. In the subsequent continuous displacement stage, each curve exhibits sustained fluctuating behavior, attributed to complex structural deformations and continuous material damage evolution under impact.

Figure 6: Force-displacement curves of impact specimens: (a) The curves of B-C6, B-[C0/CG90]6 and B-[C0/G90]6; (b) The curves of B-CG6, B-[C0/CG90]6 and B-G6.

Fig. 7 illustrates the damage deformation characteristics and performance correlations of the six hat-shaped beam specimens over the displacement range of 0–25 mm. Fig. 8 presents the specimens after the impact test. Specifically, subsequent to elastic deformation, the B-C6 specimen is dominated by fiber fragmentation, thereby introducing a risk of sudden failure. In contrast, both the B-[C0/CG90]6 and B-[C0/G90]6 specimens exhibit mild wrinkling during damage progression, which enables them to retain high energy absorption capacity while suppressing brittle failure. For the B-CG6 and B-[CG0/G90]6 specimens, the relatively low carbon fiber content in the loading direction results in the development of significant wrinkling at approximately 10 mm displacement. Notably, the B-G6 specimen undergoes premature buckling, ultimately leading to limited energy absorption efficiency.

Figure 7: Impact process of typical hat-shaped beams.

Figure 8: Typical final deformation mode of hat-shaped beams.

The energy absorption performance of all hat-shaped beam specimens is summarized in Table 4. The B-C6 specimen achieves the highest EA and SEA, while the B-G6 specimen exhibits the lowest. For intra-layer hybrid specimens, their EA and SEA values fall between those of B-C6 and B-G6 counterparts, validating that hybrid weaving balances the energy absorption characteristics of distinct fiber types. As illustrated in Table 4, the cost-efficiency index highlights practical application potential. B-C6 specimen has an EA/Cost of 33.80 J/CNY, whereas the B-[C0/G90]6 hybrid specimen reaches 42.82 J/CNY (a 27% enhancement over the B-C6 specimen). This improvement stems from the lower cost of glass fibers, which reduces overall material cost while the hybrid structure retains high energy absorption. In contrast, despite a relatively high EA/Cost, the B-G6 specimen’s remarkably low total EA renders it unsuitable for high-energy-dissipation scenarios. In conclusion, intra-layer hybrid weaving configurations significantly regulate the impact performance of hat-shaped beams, with the B-[C0/G90]6 specimen standing out as the optimal design due to its balanced load-bearing capacity, energy absorption efficiency, and cost-effectiveness.

3  Multi-Scale Modeling Strategy

To realize the multiscale collaborative design of intra-layer hybrid composite beams, a multiscale model is developed herein by incorporating the MMF criterion with newly proposed damage evolution laws. This section presents the detailed modeling methods.

3.1 Calculation of Microscopic Stress

HFRP composite materials exhibit prominent multiscale characteristics, with robust correlations in mechanical properties across disparate scales. The multiscale model developed in this study encompasses microscopic, mesoscopic, and macroscopic hierarchies, where cross-scale coupling between adjacent levels is accomplished through the stress amplification factors (SAFs) matrices Mmeso and Mmicro, as delineated in Eqs. (5) and (6) [22–25]:

σmeso=Mmesoσ¯macro(5)

σmicro=Mmicroσ¯meso(6)

where σmicro, σ¯meso and σ¯macro respectively denote microscopic, mesoscopic and macroscopic stresses, with their unified stress tensor expressions presented in Eq. (7):

{σmicro=[σ11σ22σ33σ12σ23σ13]Tσ¯meso=[σ¯11σ¯22σ¯33σ¯12σ¯23σ¯13]Tσ¯macro=[σ¯11σ¯22σ¯33σ¯12σ¯23σ¯13]T(7)

Through the integration of Eqs. (5) and (6), the σmicro can be obtained by σ¯macro and SAFs matrix Mmicro−macro, which can be defined as:

Mmicro−macro=Mmeso/Mmicro(8)

Acting as a pivotal bridging component between the macroscopic and microscopic scales, the specific values of the Mmicro−macro matrix can be quantified utilizing micro- and mesoscopic representative volume elements (RVEs) models.

3.2 Establishment of RVEs Models

To establish more accurate RVE models, the volume fractions of composite specimens were quantified via scanning electron microscopy (SEM) observations, as presented in Fig. 9. Image processing was employed to calculate the microscopic fiber volume fractions of CFs and GFs, yielding values of 65.18% and 70.04%, respectively. For modeling simplification, the fiber bundles were assumed to have an elliptical cross-sectional geometry and a sinusoidal path profile. The mesoscopic volume fraction of the intra-layer hybrid composite was determined to be 61% through image analysis.

Figure 9: Micrograph of fibers and fiber yarns: (a) SEM images of fibers; (b) Cross-sectional micrograph of fiber yarns.

For the microscopic RVE model, a cuboidal configuration was employed, with a length and width of 1 mm, respectively, and a height of 3 mm. The average diameters of CFs and GFs were determined as 0.42 and 0.44 mm, correspondingly. Furthermore, the mesoscopic RVE model is predominantly composed of warp yarns, fill yarns, and matrix. With respect to the mesoscopic RVE model, the length and width were both specified as 2 mm, while the thickness was set to 0.15 mm. The yarn width and spacing were defined as 0.96 and 0.04 mm, respectively. Tetrahedral meshes tend to undergo significant distortion; therefore, the fibers and matrix were discretized into hexahedral elements (C3D8R). Mesh sensitivity analyses were performed separately for the microscopic and mesoscopic RVE models to ensure the accuracy and convergence of the stress SAFs predictions and damage evolution simulations. The microscopic and mesoscopic RVE models are illustrated in Fig. 10.

Figure 10: Microscopic and mesoscopic unit cell models: (a) Microscopic RVE model; (b) Mesoscopic RVE model.

According to Eqs. (5) and (6), the SAFs shall be extracted from each node of the RVE models at both the microscopic and mesoscopic perspectives. In the present study, the k-means clustering algorithm was employed to classify the SAFs of all elements in the RVE model into different categories, and the average value of the stress amplification factor matrices in each category is taken as the representative of the corresponding elements and used as the reference point [28]. The workflow of the proposed method is depicted in Fig. 11.

Figure 11: K-means clustering procedure for reference point determination.

When conducting micromechanics-driven analytical investigations, to ensure stress continuity and deformation compatibility of the microscale RVE model under applied loading, periodic boundary conditions must be imposed on it [29–32]. For each node pair situated on the opposing boundary surfaces of the model, the corresponding displacement parameters can be accurately characterized via Eq. (9).

uij+(x,y,z)−uij−(x,y,z)=cij(i,j=1,2,3)(9)

where uij+ and uij− denote the i-directional displacements corresponding to the boundary faces with the j-direction normal vector, the superscripts “+” and “−” are used to distinguish the node pairs on the opposing boundary faces. Meanwhile, the constant cij represents the displacement differential between the two paired nodes.

Based on Eq. (5), the Mmicro−macro can be determined as follows: [22–25].

Mmeso/Mmicro=ψψ¯−1(10)

{ψ¯=[Ω¯1Ω¯2Ω¯3Ω¯4Ω¯5Ω¯6]6×6ψ=[Ω1Ω2Ω3Ω4Ω5Ω6]6×6(11)

{Ω¯i=[σ¯11σ¯22σ¯33σ¯12σ¯23σ¯13]T,Ωi=[σ11σ22σ33σ12σ23σ13]T,i=11,22,33,12,23,13(12)

where ψ and ψ¯ correspond to the microscale stress and macroscale stress under six independent loading cases, respectively; Ω¯i denotes the unit stress load exerted on the RVE model, while Ωi signifies the microscopic stress at the reference point. Given that ψ¯ is an unit matrix, Eq. (8) can be reformulated as Eq. (13)

Mmicro−macro=ψ=[Ω1Ω2Ω3Ω4Ω5Ω6]6×6(13)

Accordingly, Mmicro−macro represents the microscale stress intensity ψ6×6.

The SAFs and microscopic stress values associated with different constituents vary with damage states. Adopting SAFs determined from a single damage state as the element failure criterion is inconsistent with the actual mechanical response. Therefore, it is necessary to calculate these factors for all damage states of the individual material constituents. Given the wide variety of damage scenarios for different material components, direct computation for each damage case would incur an excessively high computational cost. Therefore, a linear fitting method is employed in this study to evaluate these parameters. The material properties of the fiber and matrix are pre-degraded, and the corresponding SAFs are subsequently obtained. By employing matrix damage as the independent variable and establishing a linear fitting equation, the expressions for the fiber and matrix under different damage states can be derived, which allows the microscopic stress amplification to be updated in real time with damage evolution.

3.3 Failure Criteria Based on Micromechanics

In the investigation of microscopic damage behavior in composite materials, prior studies have established the MMF criterion [22,23]: by elucidating the damage mechanisms of composites, this criterion enables the prediction of microscopic damage states in fiber yarns and resin phases, thereby providing a theoretical framework for quantitative microscopic damage analysis. In the present section, the MMF failure theory is applied to composite systems, facilitating the effective characterization of their microscopic damage behavior.

3.3.1 Criterion for Fiber Failure

Extensive studies have demonstrated that the axial strength of fibers is substantially superior to their radial strength, thus enabling fibers to bear the predominant portion of longitudinal loads. Therefore, the maximum stress failure criterion was adopted to characterize the microscopic failure behavior of fibers, as formulated in Eqs. (14) and (15) [23,24].

σ1Xf0,T=1(14)

−σ1Xf0,C=1(15)

In the above expressions, superscripts “T” and “C” denote tensile and compressive loading states acting on the fibers, respectively; σ1 represents the longitudinal microscopic stress component. Xf0,T and Xf0,C correspond to the intrinsic tensile and compressive strengths of the fibers, respectively. The onset of fiber damage is deemed to occur when the dimensionless ratios σ1Xf0,T and −σ1Xf0,C reach or exceed unity, which indicates that the microscopic stress has attained the critical strength threshold of the fiber phase.

3.3.2 Criterion for Matrix Failure

It is well-established that the failure behavior of the matrix phase is governed by hydrostatic stress-dependent criteria. Therefore, the modified Von Mises criterion was utilized to characterize the damage of the matrix material, as mathematically defined in Eqs. (16) and (17) [32,33].

(σvmσvmcr)2+I1I1cr=1,σvmcr=Ym0,TYm0,C,I1cr=Ym0,TYm0,CYm0,C−Ym0,T(16)

{I1=σ1+σ2+σ3σvm=12[(σ11−σ22)2+(σ22−σ33)2+(σ11−σ33)2]+3(σ122+σ232+σ132)(17)

where, σ11, σ22, and σ33 represent the microscopic normal stress components corresponding to the 11, 22 and 33 tensile directions in the constitutive matrix, while σ12, σ23 and σ13 denote the microscopic shear stress components associated with the 12, 23 and 13 shear directions, respectively; Ym0,T and Ym0,C correspond to the tensile and compressive strengths of the undamaged matrix, respectively; σvm and σvmcr represent the von Mises equivalent stress of the matrix and its critical value, respectively; and I1 and I1cr denote the first stress invariant of the matrix and its critical value, respectively. Matrix damage initiation is deemed to occur when the dimensionless term (σvmσvmcr)2+I1I1cr reaches or exceeds unity, signifying that the matrix stress has attained its critical failure threshold.

3.4 Damage Evolution Criterion

In the present section, the damage evolution criteria of fiber and matrix phases and the computational method for macroscopic damage parameters are thoroughly described, whereby the update of the macroscopic stiffness matrices are realized.

3.4.1 Fiber Damage Evolution Criterion

Upon the fiber stress satisfying the threshold of the relevant failure criterion, fiber stiffness exhibits a drastic reduction. Eq. (18) [26] formulates the damage evolution law for fibers, and the constitutive behavior of micro-damaged fibers is mathematically described by Eq. (19) [34].

{DfT=0.9,σf>TfDfT(C)=0,−Cf<σf<TfDfC=0.9,σf<−Cf(18)

σf=(1−DfT(C))Cfεf,Ef1d,T(C)=(1−DfT(C))Ef1,Ef2d,T(C)=(1−DfT(C))Ef2(19)

where σf and εf correspond to the microscopic stress and strain of the fibers, respectively; Cf denotes the initial stiffness matrix of fibers; Ef1d,T(C) and Ef2d,T(C) correspond to the longitudinal and transverse elastic moduli of post-damage fibers, whereas Ef1 and Ef2 represent the equivalent elastic moduli of the pristine fiber counterparts. The fiber damage variable is denoted by DfT(C).

3.4.2 Matrix Damage Evolution Criterion

Upon satisfying the tensile or compressive failure criterion, damage and failure initiate in the matrix. In contrast to fiber damage, the matrix exhibits a progressive pattern of damage evolution. The damage evolution process of the matrix is shown in Eqs. (20)–(22) [33–35].

σeq=(β−1)I1+(β−1)2I12+4βσvm22β,β=Ym0,CYm0,T(20)

κm=max(σeq)(21)

DmT(C)=1−exp⁡[γ(1−κmYm0,T)](22)

where κm denotes the maximum equivalent stress (σeq) sustained by the matrix following the initiation of damage. γ represents the damage shape parameter, assigned a value of 1.5 for the purposes of this study. DmT(C) refers to the matrix damage factor, a stress-state-dependent parameter bounded within the range of 0 to 1.

The constitutive relation governing damaged matrix is mathematically formulated as Eq. (23) [30,31].

σm=(1−DmT(C))Cmεm,Emd,T(C)=(1−DmT(C))Em(23)

where σm and εm represent the microscopic stress and strain of the matrix, respectively; Cm denotes the initial stiffness matrix of matrix; Em and Emd,T(C) correspond to the intrinsic initial elastic modulus and post-damage elastic modulus of the matrix phase, respectively.

3.4.3 Macroscopic Stiffness Matrix Update

In the present study, macroscopic damage variables are employed to account for the microscale effects on the macroscale behavior. Considering that the research focuses on plain-woven hybrid fiber composites, the definition and computation of the macroscopic damage variables are given by Eqs. (24) and (25) [33,36]:

{d1T(C)=1−(Ef1Cd,T(C)Vf1+Emd,T(C)Vm)/E11d2T(C)=1−(Ef1Gd,T(C)Vf2+Emd,T(C)Vm)/E22d3T(C)=1−Emd,T(C)/(1−Vf(1−Emd,T(C)/Ef2d,T(C)))/E33(24)

{d1=1−(1−d1T)(1−d1C)d2=1−(1−d2T)(1−d2C)d3=1−(1−d3T)(1−d3C)(25)

where E11, E22 and E33 correspond to the elastic moduli of the composite in the 11, 22, and 33 directions, respectively; Vf1, Vf2 and Vm denote the volume proportions of CFs, GFs, and matrix phase, with Vf representing the aggregate volume proportion of carbon and glass fibers. d1, d2 and d3 stand for the macroscopic damage parameters along the 11, 22, and 33 directions of the composite, respectively. Following the onset of material damage, these macroscopic damage parameters are computed from the microscopic damage factors, thereby facilitating the update of the macroscopic damaged stiffness matrix via Eq. (26) [37]:

Cd=[(1−d1)C110(1−d1)(1−d2)C120(1−d1)(1−d3)C130000(1−d1)(1−d2)C120(1−d2)C220(1−d2)(1−d3)C230000(1−d1)(1−d3)C130(1−d2)(1−d3)C230(1−d3)C330000000(1−d1)(1−d2)C440000000(1−d1)(1−d3)C550000000(1−d2)(1−d3)C660](26)

where Cij0 refers to the undamaged macroscopic stiffness coefficient.

Fig. 12 depicts the architecture of the developed multi-scale analytical approach. Initially, the volume fractions of CFs and GFs are quantified based on their spatial distribution, with the concurrent establishment of microscopic RVE models. Subsequently, the mesoscopic RVE models for intra-layer hybrid composites are constructed by incorporating the geometric characteristics of fiber yarns, including their path configurations and cross-sectional profiles. To establish inter-scale correlations, PBCs and six sets of macroscopic unit stress loads are imposed on the microscopic and mesoscopic RVE models, respectively. Reference points corresponding to fibers, matrix phase, and fiber yarns are identified via k-means clustering analysis, enabling the calculation of the Mmicro and Mmeso. Macroscopic stresses are transformed to microscopic stresses via SAFs. The microscopic damage states are further mapped to macroscopic damage variables. The multi-scale model achieves quantitative correlation between the microstructural morphology, damage characteristics, and macroscopic mechanical response.

Figure 12: Computation framework of the proposed multi-scale model.

4  Numerical Results and Analysis

4.1 Validation of RVE Models

The elastic properties of CFs, GFs, and matrix utilized for developing microscopic and mesoscopic RVE models were derived from low-load tensile test data and material parameters furnished by raw material manufacturers, as collated in Table 5. Table 6 presents the periodic boundary conditions and six sets of macroscopic unit stress excitations. Representative stress nephograms of microscopic RVE models corresponding to CFs and GFs are presented in Fig. 13, respectively. Furthermore, the characteristic stress distribution profiles of mesoscopic RVE models are illustrated in Fig. 14, correspondingly. The findings demonstrate that under varied loading scenarios, considerable discrepancies exist in the stress magnitudes and spatial distribution traits of the composites, with the stress fields displaying an approximately symmetric configuration. With the computational expressions provided in Table 7 adopted for analysis, the elastic parameters of the fiber bundles were acquired and tabulated in Table 8.

Figure 13: Microscopic RVE models under unit stress excitations.

Figure 14: Mesoscopic RVE models under unit stress excitations.

To verify the accuracy of the established RVE models, a comparison was conducted between the elastic modulus obtained from tensile tests and that predicted by the multi-scale model. As shown in Table 9, the relative error of the simulation results with respect to the experimental results is within 10%. Accordingly, the proposed RVE models show reliable predictive performance for the macroscopic elastic characteristics of hybrid composites, and thus provide adequate accuracy for the computation of the SAFs matrix Mmicro−macro.

4.2 Validation of the Multi-Scale Model

Fig. 15 presents the finite element model for impact simulation of the hat-shaped beam. The key modeling parameters and configurations are elaborated as follows: the hat-shaped beam was discretized using C3D8R hexahedral solid elements with a mesh size of 1.5 mm × 0.8 mm × 0.15 mm, and all layers were interconnected via cohesive contact behavior [38]. Model parameters adopted for cohesive contact model are presented in Table 10, with analogous materials adopted for the simulation. A 45° chamfer was incorporated at the top section of the beam to faithfully replicate the geometric characteristics of the physical prototype. Within the impact test assembly, the lower steel plate was subjected to fully constrained boundary conditions to simulate a rigid support configuration, while the upper impactor steel plate was assigned a mass of 20.1 kg. Based on the target impact energy of 200 J, an initial impact velocity of 4.46 m/s was derived via the kinetic energy principle and imposed vertically downward on the upper plate to achieve the predefined energy-driven impact loading. In addition to the material properties presented in Table 5, additional parameters required for the development of the proposed multi-scale model are summarized in Table 11.

Figure 15: Finite element model of hat-shaped beam under low-velocity impact.

Fig. 16 illustrates a comparison between the experimental and numerical results regarding the deformation patterns, force responses, and displacement behaviors of the B-[C0/G90]6. Fig. 16a presents a comparative analysis of the force-displacement responses obtained from experimental measurements and numerical simulations, with subtle discrepancies identified between the two datasets. These deviations may be partially attributed to inherent inconsistencies between the numerical model and the physical prototype—encompassing disparities in geometric characteristics, mismatches between parameters and the actual material characteristics, and other relevant factors. Furthermore, manufacturing imperfections in the test specimens constitute an additional contributing factor to such variations. Despite these minor deviations, the overall evolutionary trends remain in good agreement. Additionally, Fig. 16b presents the impact deformation evolution of hat-shaped beams, where the red elements denote the failed regions. Complex deformation modes are well reproduced by the multi-scale model.

Figure 16: Numerical and experimental result comparison.

As shown in Table 12, a comprehensive comparison was conducted between the experimental and simulated performance indicators of the hat-shaped beam, including Fmax, EA, SEA, and EA/Cost. The results demonstrate a strong correlation between the experimental measurements and numerical predictions. The relative errors of all indicators are well within 10%. This confirms that the established multi-scale model has high accuracy and reliability, which can effectively support the subsequent investigation on the mechanical response of hybrid composite beams.

4.3 Discrete Optimization of Woven Hybrid Fiber Composites

As revealed by the preceding experimental and multi-scale simulation results, the intra-layer hybridization parameters (warp carbon fiber content, weft carbon fiber content, and ply angle) exert a significant coupled influence on the low-velocity impact energy absorption and material cost of woven hybrid fiber composites. To achieve an optimal balance between energy absorption efficiency and economic feasibility, the hybridization parameters were further determined in this section to obtain the optimal design scheme. In practical engineering problems, design variables are generally discrete values. Thus, an orthogonal discrete optimization method [39,40] was adopted for the optimization, and the detailed optimization steps are shown in Fig. 17.

Figure 17: Orthogonal discrete optimization flow chart.

Based on practical engineering application requirements, the energy absorption per unit cost index (EA/Cost) was selected as the core optimization objective. EA was derived from the impact simulation results of woven fabrics (validated in Section 4.2), where the simulation process fully incorporated the inter-fiber interactions and yarn arrangement characteristics inherent to the woven structure. The material cost was calculated based on the mass fractions of carbon fiber (0.5 CNY/g), glass fiber (0.1 CNY/g), and resin (0.05 CNY/g), and revised by accounting for the material utilization rate of the weaving process. The design variables included the proportion of carbon fiber bundles in the warp and weft yarns (0%, 25%, 50%, 75%, 100%) and the ply angle of the fiber fabrics (0°, 30°, 45°, 60°, 90°). To ensure the optimized composites meet practical service requirements, explicit constraints were imposed on energy absorption, tensile strength, and material cost, respectively, based on the experimental data in Table 4 and engineering practice experience with woven structural components. This optimization problem can be mathematically characterized by Eq. (27).

{R=Max{EA(Vf1,Vf2,A)/C(Vf1,Vf2,A)}s.t.{EA ≥120Jσ≥400 MpaCost ≤4.6 CNYVf∈(0%,25%,50%,75%,100%)A∈(0∘,30∘,45∘,60∘,90∘)(27)

where R denoted the objective function; Vf1,Vf2 were the proportions of carbon fiber bundles in the warp and weft yarns; EA, σ and Cost presented the energy absorption, the tensile strength, and the material cost of the hybrid specimen, respectively. To solve the problem that there were sample points unsatisfying the constraints, the concept of penalty function was introduced to transform the optimization problem with constraints into an optimization problem without constraints, as shown in Eq. (28).

{Rnew=Max{EA(Vf1,Vf2,A)/C(Vf1,Vf2,A)+φ1(Vf1,Vf2,A)+φ2(Vf1,Vf2,A)+φ3(Vf1,Vf2,A)}Vf∈(0%,25%,50%,75%,100%)A∈(0∘,30∘,45∘,60∘,90∘)φ1(Vf1,Vf2,A)=s1*Max[0,EA(Vf1,Vf2,A)−120]φ2(Vf1,Vf2,A)=s2*Max[0,σ(Vf1,Vf2,A)−400]φ3(Vf1,Vf2,A)=s3*Max[0,C(Vf1,Vf2,A)−6](28)

where φ1,φ2,φ3 was the penalty function and s1,s2,s3 was the penalty factor. When the calculated value of the sample point exceeded the constraints, there will be a certain degree of revision on the objective function value Rnew according to the penalty factor.

An orthogonal experimental array was developed with a 3-factor, 3-level orthogonal experimental design. To avert the emergence of local optima in the optimization process, the median value of each design variable was set as the middle level in the initial iteration, with its adjacent values defined as Level 1 and Level 3, respectively. The level values of the design variables for the first iteration are given in Table 13. For each experimental group, the multi-scale simulation model was used to predict EA, σ and Cost, and then Rnew was calculated.

To further verify the statistical significance of the investigated factors, analysis of variance (ANOVA) was carried out based on the orthogonal experimental results, as presented in Table 14. In the statistical evaluation, a factor with a p-value less than 0.05 was regarded as statistically significant, whereas those with higher p-values were considered insignificant [41]. The ANOVA results demonstrate that the regression model is highly significant (p < 0.001), and the ply angle A exhibits a statistically significant effect (p < 0.001). These statistical findings verify that the identified factor ranking and optimization trends are statistically reliable rather than accidental.

The iterative process of Rnew is illustrated in Fig. 18. The optimal design was achieved at the third iteration, and since there was no further improvement of more than 5% compared to the third iteration, the iterative process was deemed to have converged. The optimal design at each iteration step is listed in Table 15. Table 16 presents the design variables and Rnew values for both the initial and optimal designs, from which it can be observed that EA/Cost increased from the initial value of 42.82 to 55.13, representing an improvement of 28.74%.

Figure 18: Discrete optimization iteration curve.

5  Conclusions

This paper presents a series of tensile and impact tests on carbon/glass HFRP plates and hat-shaped beams, and discusses the effects of different hybridization configurations on their mechanical properties. Based on the experimental results, a multi-scale damage evolution model grounded in the MMF criterion is established and implemented via Abaqus/VUMAT. Using the validated multi-scale model, orthogonal discrete optimization is conducted with the mesoscopic carbon fiber bundle fraction and the macroscopic ply angle as design variables. The influences of different design variables on the energy absorption efficiency per unit cost (EA/Cost) are further explored. Within these constraints, the main conclusions of this study can be summarized as follows:

(1)   The tensile properties of intra-layer hybrid composites are significantly dependent on fiber orientation and hybridization schemes. When carbon fibers are arranged along the loading direction, the longitudinal elastic modulus and tensile strength are superior to those of pure CFRP, demonstrating obvious hybrid synergy. In contrast, when glass fibers dominate the loading direction, the material properties are closer to pure GFRP.

(2)   For single hat-shaped beams, the impact performance of intra-layer hybrid composites lies between pure CFRP and GFRP. Among them, the B-[C0/G90]6 configuration shows the optimal comprehensive performance. It maintains high EA (157.25 J) and SEA (3.21 J/g), and its cost-efficiency index (EA/Cost = 42.82 J/CNY) is 17% higher than that of pure CFRP, attributing to the ductility of glass fibers that suppresses brittle failure and reduces material costs.

(3)   The established multi-scale model, integrating the micromechanics failure criterion and stress amplification factors, successfully bridges the micro-meso-macro scales. This model accurately predicts the impact force-displacement response and performance indicators of hat-shaped beams. Compared with traditional macro-mechanical models, it considers the influence of microscopic parameters on macroscopic failure behavior, realizing the quantitative mapping of microstructural characteristics to macroscopic performance.

(4)   The orthogonal discrete optimization results indicate that the ply angle has the most significant impact on the cost-efficiency index (EA/Cost), followed by the warp carbon fiber content, while the weft carbon fiber content has the least influence. This highlights that adjusting the ply angle is the most effective way to balance energy absorption and cost. Compared with the initial design, the EA/Cost of the optimized structure increases from 42.82 to 55.13 J/CNY, with an improvement of 28.74%.

(5)   The multi-scale model and optimal hybrid architecture obtained in this study can be generalized to other thin-walled energy-absorbing structures with similar loading and deformation mechanisms. The model can also be extended to higher-energy impact cases by adjusting impact parameters and boundary conditions. However, the present model adopts idealized yarn geometry and ignores strain rate effects, which need to be further improved for complex structures and extreme impact scenarios.

Acknowledgement: The authors appreciate the support from the research group and relevant institution during the completion of this work.

Funding Statement: This research is funded by Shaanxi Province Natural Science Foundation (2024JC-YBMS-267, 2025JC-YBMS-447), Shaanxi Province Postdoctoral Research Project Funding (2023BSHYDZZ124), China Postdoctoral Science Foundation (2023M740358).

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Chaoshuai Duan, Guohua Zhu; data collection: Yin Wang; analysis and interpretation of results: Chaoshuai Duan, Xiaotian Zhang; draft manuscript preparation: Jiale Wang, Zhen Wang. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author, Guohua Zhu, upon reasonable request.

Ethics Approval: Not applicable. This study did not involve any human participants or animal subjects, and no ethical approval was required.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Qin W, Fernand M, Yang ZJ, Georges I, Zakari MD. Effect of fiber volume ratio on flexural behavior of RC beam with novel HFRP rebar and steel rebar. Structures. 2022;43:1882–900. doi:10.1016/j.istruc.2022.06.023. [Google Scholar] [CrossRef]

2. Zhu G, Shao Z, Cao X-A, Feng J, Wang Z, Zhao X. Material and structure integrated design for 3D-printed honeycomb beams. Eng Struct. 2026;361:122938. doi:10.1016/j.engstruct.2026.122938. [Google Scholar] [CrossRef]

3. Wang Z, Zhang W, Luo Q, Zheng G, Li Q, Sun G. A novel failure criterion based upon forming limit curve for thermoplastic composites. Compos Part B Eng. 2020;202(9):108320. doi:10.1016/j.compositesb.2020.108320. [Google Scholar] [CrossRef]

4. Ishtiaq S, Saleem MQ, Naveed R, Harris M, Khan SA. Glass–Carbon–Kevlar fiber reinforced hybrid polymer composite (HPCpart (A) mechanical and thermal characterization for high GSM laminates. Compos Part C Open Access. 2024;14(2):100468. doi:10.1016/j.jcomc.2024.100468. [Google Scholar] [CrossRef]

5. Uddin MS. The failure analysis of carbon fiber-reinforced epoxy composites against impact loading with numerical and experimental investigations. J Polym Mater. 2025;42(4):1051–73. doi:10.32604/jpm.2025.070688; [Google Scholar] [CrossRef]

6. Fu Z, Du Y, Guo A, Zhou F. Effects of impregnation resin elastic modulus and slack length on the tensile properties of carbon fiber bundles. Constr Build Mater. 2025;503(14376):144574. doi:10.1016/j.conbuildmat.2025.144574. [Google Scholar] [CrossRef]

7. Milosan I, Bedo T, Gabor C, Pop MA. Mechanical characteristics of glass-fiber-reinforced polyester composite materials. Materials. 2025;18(15):3595. doi:10.3390/ma18153595. [Google Scholar] [CrossRef]

8. Chelot D, Ravichandran S, Upadhyaya P. A sustainable solution with improved chemical resilience using repurposed glass fibers for sewage rehabilitation pipes. Recycling. 2024;9(2):28. doi:10.3390/recycling9020028. [Google Scholar] [CrossRef]

9. Mou H, Chen Y, Xie J, Liu H. Crashworthiness analysis of single CFRP, hybrid CFRP/Al and fiber metal laminate C-channels: quasi-static crushing tests and numerical simulations. Thin Walled Struct. 2025;213(16):113226. doi:10.1016/j.tws.2025.113226. [Google Scholar] [CrossRef]

10. Ma B, Cao X, Feng Y, Song Y, Yang F, Li Y, et al. A comparative study on the low velocity impact behavior of UD, woven, and hybrid UD/woven FRP composite laminates. Compos Part B Eng. 2024;271(1):111133. doi:10.1016/j.compositesb.2023.111133. [Google Scholar] [CrossRef]

11. Lu R, Zhang Q, Zhang H, Meloni M, Feng J, Cai J. Quasi-static crushing response of a novel triaxial isotropy mechanical metamaterial with dual-platform property. Thin Walled Struct. 2025;206(Part A):112630. doi:10.1016/j.tws.2024.112630. [Google Scholar] [CrossRef]

12. Wang M, Shi L, Pan Z, Wu Z. Mechanical response and damage mechanism of carbon/Kevlar hybrid braiding composite laminates under quasi-static indentation and low-velocity impact loading. Thin Walled Struct. 2025;213(6):113291. doi:10.1016/j.tws.2025.113291. [Google Scholar] [CrossRef]

13. Kazemi ME, Shanmugam L, Dadashi A, Shakouri M, Lu D, Du Z, et al. Investigating the roles of fiber, resin, and stacking sequence on the low-velocity impact response of novel hybrid thermoplastic composites. Compos Part B Eng. 2021;207(8):108554. doi:10.1016/j.compositesb.2020.108554. [Google Scholar] [CrossRef]

14. Maier R, Mandoc AC. Investigation on layer hybridization of glass/carbon fibre woven reinforced composites subjected to low-speed impact. J Compos Sci. 2023;7(2):83. doi:10.3390/jcs7020083. [Google Scholar] [CrossRef]

15. Ahamed MAM, Dhakal HN, Zhang Z, Barouni A, Pillai JR. Low-velocity impact damage characteristics of flax/glass epoxy hybrid laminates on the influence of different temperatures: experimental and numerical analysis. Compos Struct. 2025;359:113337. doi:10.1016/j.compstruct.2024.118704. [Google Scholar] [CrossRef]

16. Bussadori BP, Schuffenhauer K, Scattina A. Modelling of CFRP crushing structures in explicit crash analysis. Compos Part B Eng. 2014;60(2):725–35. doi:10.1016/j.compositesb.2014.01.020. [Google Scholar] [CrossRef]

17. Czél G, Wisnom MR. Demonstration of pseudo-ductility in high performance glass/epoxy composites by hybridisation with thin-ply carbon prepreg. Compos Part A Appl Sci Manuf. 2013;52(12):23–30. doi:10.1016/j.compositesa.2013.04.006. [Google Scholar] [CrossRef]

18. Jalalvand M, Czél G, Wisnom MR. Numerical modelling of the damage modes in UD thin carbon/glass hybrid laminates. Compos Sci Technol. 2014;94:39–47. doi:10.1016/j.compscitech.2014.01.013. [Google Scholar] [CrossRef]

19. Ge L, Li H, Liu B, Fang D. Multi-scale elastic property prediction of 3D five-directional braided composites considering pore defects. Compos Struct. 2020;244(13):112287. doi:10.1016/j.compstruct.2020.112287. [Google Scholar] [CrossRef]

20. Xiong X, Hua L, Miao M, Shen SZ, Li X, Wan X, et al. Multi-scale constitutive modeling of natural fiber fabric reinforced composites. Compos Part A Appl Sci Manuf. 2018;115(6):383–96. doi:10.1016/j.compositesa.2018.10.016. [Google Scholar] [CrossRef]

21. McWilliams B, Dibelka J, Yen CF. Multi scale modeling and characterization of inelastic deformation mechanisms in continuous fiber and 2D woven fabric reinforced metal matrix composites. Mater Sci Eng A. 2014;618(5):142–52. doi:10.1016/j.msea.2014.08.063. [Google Scholar] [CrossRef]

22. Ha SK, Jin KK, Huang Y. Micro-mechanics of failure (MMF) for continuous fiber reinforced composites. J Compos Mater. 2008;42(18):1873–95. doi:10.1177/0021998308093911; [Google Scholar] [CrossRef]

23. Jin KK, Huang Y, Lee YH, Ha SK. Distribution of micro stresses and interfacial tractions in unidirectional composites. J Compos Mater. 2008;42(18):1825–49. doi:10.1177/0021998308093909; [Google Scholar] [CrossRef]

24. Huang Y, Xu L, Kyu Ha S. Prediction of three-dimensional composite laminate response using micromechanics of failure. J Compos Mater. 2012;46(19–20):2431–42. doi:10.1177/0021998312449888; [Google Scholar] [CrossRef]

25. Huang Y, Jin C, Ha SK. Strength prediction of triaxially loaded composites using a progressive damage model based on micromechanics of failure. J Compos Mater. 2013;47(6–7):777–92. doi:10.1177/0021998312460261; [Google Scholar] [CrossRef]

26. Xu L, Jin CZ, Ha SK. Ultimate strength prediction of braided textile composites using a multi-scale approach. J Compos Mater. 2015;49(4):477–94. doi:10.1177/0021998314521062; [Google Scholar] [CrossRef]

27. Ha SK, Huang Y, Han HH, Jin KK. Micromechanics of failure for ultimate strength predictions of composite laminates. J Compos Mater. 2010;44(20):2347–61. doi:10.1177/0021998310372464; [Google Scholar] [CrossRef]

28. Wang M, Zhang P, Fei Q, Guo F. Modified micro-mechanics based multiscale model for progressive failure prediction of 2D twill woven composites. Chin J Aeronaut. 2020;33(7):2070–87. doi:10.1016/j.cja.2019.10.009. [Google Scholar] [CrossRef]

29. Xia Z, Zhang Y, Ellyin F. A unified periodical boundary conditions for representative volume elements of composites and applications. Int J Solids Struct. 2003;40(8):1907–21. doi:10.1016/S0020-7683(03)00024-6; [Google Scholar] [CrossRef]

30. Jing K, Zhou H, Wang H, Yan H, Xie S. Multiscale damage and low-velocity impact study of three-dimensional woven composites. Thin Walled Struct. 2024;202(6):112132. doi:10.1016/j.tws.2024.112132. [Google Scholar] [CrossRef]

31. Wang D, Huang H, Guo Z, Shan Z, Liu J, Zhang H. Study on the mechanism of Z-directional fiber reinforcement of carbon nanotubes in flexible-oriented three dimensional woven composites. Carbon. 2025;243(13):120584. doi:10.1016/j.carbon.2025.120584. [Google Scholar] [CrossRef]

32. Sun XS, Tan VBC, Tay TE. Micromechanics-based progressive failure analysis of fibre-reinforced composites with non-iterative element-failure method. Comput Struct. 2011;89(11–12):1103–16. doi:10.1016/j.compstruc.2010.12.003. [Google Scholar] [CrossRef]

33. Muthusamy P, Sivakumar SM. A constituent-behavior-motivated model for damage in fiber reinforced composites. Comput Mater Sci. 2014;94(1):163–72. doi:10.1016/j.commatsci.2014.03.048. [Google Scholar] [CrossRef]

34. Wang L, Wang B, Wei S, Hong Y, Zheng C. Prediction of long-term fatigue life of CFRP composite hydrogen storage vessel based on micromechanics of failure. Compos Part B Eng. 2016;97(3):274–81. doi:10.1016/j.compositesb.2016.05.012. [Google Scholar] [CrossRef]

35. Zhu G, Zhang Y, Niu X, Duan C, Wang Z, Zhao X. Novel multiscale modeling strategy for hybrid fiber reinforced composites. Int J Mech Sci. 2025;287(4):109997. doi:10.1016/j.ijmecsci.2025.109997. [Google Scholar] [CrossRef]

36. Hyer MW, Waas AM. Micromechanics of linear elastic continuous fiber composites. In: Comprehensive composite materials. Amsterdam, The Netherlands: Elsevier; 2000. p. 345–75. [Google Scholar]

37. Zheng T, Guo L, Sun R, Wang T, Hong C, Benedictus R, et al. Investigation on the effect of interface properties on compressive failure behavior of 3D woven composites through micromechanics-based multiscale damage model. Compos Struct. 2023;320(2071):117186. doi:10.1016/j.compstruct.2023.117186. [Google Scholar] [CrossRef]

38. Cao XA, Zhu G, Wang Y, Wang Z, Zhao X. Electromechanical behavior and damage monitoring mechanism of CNTs/GFRPs composites. Int J Mech Sci. 2026;312(6348):111273. doi:10.1016/j.ijmecsci.2026.111273. [Google Scholar] [CrossRef]

39. Cao XA, Zhu G, Wang Z, Zhao X. On flexural behavior of 3D-printed continuous hybrid fiber reinforced composites: experimental and multiscale modeling study. Compos Struct. 2025;359(7):119034. doi:10.1016/j.compstruct.2025.119034. [Google Scholar] [CrossRef]

40. Lee KH, Yi JW, Park JS, Park GJ. An optimization algorithm using orthogonal arrays in discrete design space for structures. Finite Elem Anal Des. 2003;40(1):121–35. doi:10.1016/S0168-874X(03)00095-7; [Google Scholar] [CrossRef]

41. Norouzi-Inallu M, Laitinen V, Orue I, Ullakko K. Influence of laser powder bed fusion scanning strategies on microstructure, transformation, and magnetic behavior of Ni-Mn-Ga magnetic shape memory alloys. J Alloys Compd. 2025;1047:185011. doi:10.1016/j.jallcom.2025.185011. [Google Scholar] [CrossRef]