Simulation of wind-blown sand erosion damage of CFRP laminate based on discrete element method

Abstract

To investigate the wind-blown sand erosion resistance of Carbon Fiber Reinforced Polymer (CFRP) laminates, this study conducted simulated wind-blown sand erosion experiments using erosion velocity as a variable and tested the mechanical properties of the materials after erosion. Meanwhile, Furthermore, a VUMAT subroutine integrating the Discrete Element Method (DEM) with the 3D-Hashin criterion was developed in ABAQUS. This approach aims to establish a quantitative relationship between stochastic multi-particle erosion and the anisotropic micro-damage of CFRP laminate, thereby overcoming the limitations of traditional methods in quantifying erosion mechanisms. The model was first validated for accuracy using experimental results obtained at different erosion velocities. Subsequently, to further examine the model’s generalizability and prediction accuracy, supplementary experiments under different sand flow rates were conducted, and the simulation results were compared with the measured data. The results indicate that during the validation phase (different erosion velocities), the relative errors between the simulated and experimental values for the tensile strength and elastic modulus of CFRP laminate were 1.37%–4.18% and 3.56%–6.48%, respectively; whereas during the prediction validation phase (different sand flow rates), these relative error ranges were 0.83%–6.53% and 0.44%–1.4%, respectively. These results not only validate the model’s accuracy but also demonstrate its robust predictive potential for untested scenarios. The model can reliably simulate the mechanical degradation and damage evolution of CFRP laminate under real-world wind-blown sand environments, thereby providing an effective numerical tool for the durability assessment of CFRP structures.

Keywords

CFRP laminate wind-blown sand erosion discrete element method VUMAT damage evolution

Introduction

Carbon Fiber Reinforced Polymer (CFRP) is a high-performance material consisting of fibers embedded in a resin matrix, characterized by excellent properties such as high tensile strength, lightweight, durability, and corrosion resistance. CFRP can significantly enhance the durability of concrete structures in harsh environments and extend their service life; consequently, it has been widely applied in engineering practices^1,2 However, when such structures are situated in wind-sand environments, the continuous impact of high-velocity sand particles on the CFRP surface causes surface damage, significantly reducing the load-bearing capacity and service life of the material. Due to the instantaneous and violent nature of the erosion process, along with complex damage modes, it is difficult to reveal the continuous mechanism from particle impact to material failure relying solely on experimental observations^3,4 Therefore, a deep understanding of the damage mechanism of CFRP under wind-blown sand erosion, along with the establishment of effective simulation and prediction methods, is of great significance for ensuring structural safety and stability as well as extending service life^5,6

Laboratory particle erosion testing serves as a primary method for simulating the process of high-velocity particle flows impacting material surfaces. Based on the acceleration mechanisms, these tests are categorized into jet, rotary, and centrifugal erosion^7,8 with jet erosion being the most widely applied. The principle involves utilizing high-pressure gas to accelerate abrasive particles through either a Laval nozzle^9,10 or a straight pipe nozzle, generating a high-velocity jet to impinge upon the target material. Owing to its advantages—such as high particle velocity, controllable impact direction, and adjustable flow rate—this method has established itself as a standard approach for investigating erosion mechanisms and comparing material performance. In light of this, this study adopts the jet erosion method to conduct wind-blown sand erosion experiments.

Existing research indicates that erosion outcomes are influenced by erosion conditions, target material properties, and particle characteristics, among which erosion conditions exert the most significant influence. Choudhary et al.¹¹ discovered that erosion velocity governs the failure mode; specifically, FRP predominantly exhibits delamination cracking under low-velocity erosion, whereas penetration is likely to occur under high-velocity conditions. Regarding impact angle, the erosion rate of ductile materials initially increases and then decreases with increasing angle (peaking at 15°–30°), while that of brittle materials increases with the angle (peaking near 90°)^12,13 In contrast, the maximum erosion rate for semi-ductile FRP occurs at 45°–60°.¹⁴ In terms of particle parameters, Chellaganesh et al.¹⁵ noted that the contribution of particle flux to the erosion rate can reach up to 20%. Arnold et al.¹⁶ further observed that under low-flow conditions, the erosion rate increases as the flow rate decreases, a phenomenon attributed to inter-particle interactions.¹⁷ Furthermore, scholars have proposed the concept of a “critical particle size”,¹⁸ suggesting that the erosion rate remains constant once the particle size exceeds a certain threshold. Zhang et al.¹⁹ indicated that, at a constant concentration, an increase in particle size reduces the number of particles, thereby lowering the impact probability. However, for FRP materials characterized by significant anisotropy and multiphase structures, the erosion response mechanism exhibits more complex features. On the one hand, the behavior at the fiber-matrix interface undergoes nonlinear evolution under impact²⁰; on the other hand, the voids and weaving paths inherent in the fiber weave structure further enhance the non-uniformity and directionality of damage development.²¹ Consequently, traditional erosion models based on the assumption of homogeneous materials struggle to accurately describe the damage evolution patterns of FRP materials in sandstorm environments^20,22

Although the aforementioned experimental studies have revealed macroscopic laws, limitations persist in capturing microscopic mechanisms and conducting parametric analysis. Currently, research on the erosion resistance of CFRP relies primarily on methods such as sand blasting and wind tunnel tests^23–25 While these methods can provide macroscopic damage data, they are often costly and time-consuming. Furthermore, they face difficulties in capturing the microscopic evolutionary process at the instant of particle impact in real time. In terms of numerical simulation, the conventional Finite Element Method (FEM) faces significant limitations in handling the complex interactions between numerous discrete particles and continuous media.²⁶ Although meshless methods such as Smoothed Particle Hydrodynamics (SPH) have been applied to simulate material damage, their applicability to large-scale particle impact problems remains to be verified^27,28 Furthermore, although scholars have investigated the influence of various failure criteria—such as Puck, Hashin, and Chang-Chang—on the dynamic progressive failure behavior of CFRP laminate^29–32 there is a notable scarcity of research that combines efficient particle generation algorithms with refined damage models. Although the general-purpose finite element software ABAQUS is widely applied in structural analysis^33,34 the application of its built-in particle generator module remains underexplored regarding the simulation of realistic wind-sand flow fields and the refined damage mechanisms of composites under particle impact. While the FEM provides an effective means to characterize the erosion process,³⁵ existing models predominantly rely on simplified material removal criteria or empirical erosion functions, failing to simultaneously capture the anisotropic characteristics of the material and the progressive accumulation of erosion damage^36,37 Furthermore, current erosion simulation studies primarily concentrate on metallic materials or adopt simplified impact models^38–40 Research that specifically addresses the anisotropic nature of CFRP by utilizing a particle generator to simulate random impacts coupled with complex damage models remains particularly scarce.

To address the aforementioned issues, this paper proposes a numerical simulation method based on the ABAQUS platform and its particle emitter module to efficiently and accurately simulate the erosion damage process of wind-sand particles on CFRP laminate. To overcome the limitations of existing studies that are confined to macroscopic mass loss descriptions or oversimplified single-projectile models, this study constructs a multi-scale predictive model bridging random particle fields, microscopic anisotropic damage, and macroscopic performance degradation. Specifically, the DEM within ABAQUS/Explicit is employed to achieve a highly efficient equivalent simulation of real wind-blown sand flow fields. Furthermore, a VUMAT user subroutine is developed to enable the real-time capture of transient damage evolution in CFRP laminate during the erosion process. The accuracy of the model was validated by comparing the simulated damage micromorphology and mechanical properties of the CFRP laminate with experimental results under different erosion velocities. Furthermore, the model was utilized to predict damage morphology and mechanical properties under varying sand flow rates, and the reliability of these predictions was verified through actual wind-blown sand erosion tests. This model intuitively illustrates the impact process and reveals microscopic damage mechanisms, thereby providing a theoretical basis and numerical tool for the anti-erosion design and life assessment of CFRP structures in wind-sand environments.

Finite element model

3D model creation of CFRP laminate subjected to particle erosion

Figure 1 shows the finite element model of a carbon fiber reinforced polymer (CFRP) laminate subjected to particle impact. In this model, the CFRP laminate is modeled using the reduced-integration element C3D8R, with hourglass control activated to enhance computational stability.

Figure 1.

Finite element model of CFRP laminate: (a) Particle Emission surface; (b) Material orientation; (c) Meshing; (d) Boundary conditions.

As shown in Figure 1(a), a circular particle emission surface is positioned above the model. This surface is defined as a rigid body and discretized using SFM3D4R 4-node quadrilateral surface elements. The emitted discrete particles are modeled as PD3D elements with a circular shape, and their size distribution is set to range from 0.2 mm to 0.6 mm. The dimensions of the single-layer CFRP laminate are 450 mm × 75 mm × 0.167 mm. To simplify the model and reduce computational cost, the end tab regions on both sides are omitted during modeling. Regarding the contact definition, the surfaces of all incident particles were first grouped into a surface set. A surface-to-surface contact algorithm was then employed to define the interaction between the particles and the CFRP laminate. Subsequently, the contact properties in the INP file were manually modified to accurately establish the final contact definition between the two components. As illustrated in Figure 1(b), the X-axis is defined as the fiber direction and the Y-axis as the matrix direction. Furthermore, as shown in Figure 1(c), the mesh of the CFRP laminate is generated using a biased seeding strategy, with local refinement applied in the particle impact zone to improve simulation accuracy in that region. Figure 1(d) presents the boundary condition setup: reference points RP-1 and RP-2 are created at the two end faces of the specimen and are coupled to their respective faces. Constraints are applied around the specimen perimeter and at the reference points to restrict displacement in the Y-direction, preventing lateral drift of the specimen during the erosion process. Specifically, RP-1 is treated as the fixed end and is fully clamped; RP-2 is the loading end, with only the translational degree of freedom along the X-axis released, and a uniform tensile load is applied via an amplitude curve. The material constitutive relations are defined through user-defined material properties, and the specific parameters are listed in Table 1.

Table 1.

Constitutive parameters of the CFRP laminate model.

Material properties		Mechanical properties
E₁ (MPa)	234,000	X_t (MPa)	4323
E₂ (MPa)	2400	X_c (MPa)	20
E₃ (MPa)	2400	Y_t (MPa)	38
Nu₁₂	0.3	Y_c(MPa)	20
Nu₁₃	0.3	S₁₂ (MPa)	46.9
Nu₂₃	0.4	S₁₃ (MPa)	30
G₁₂ (MPa)	4000	S₂₃ (MPa)	30
G₁₃ (MPa)	4000
G₂₃ (MPa)	2500
Ρ(t/mm³)	1.6 × 10^-9

Constitutive relation of CFRP laminate

In this paper, ABAQUS software is used to simulate the residual mechanical properties of CFRP laminate subjected to wind-blown sand erosion. The 3D-Hashin criterion is adopted as the damage criterion for the composite material. The elastic constitutive equations employed^29,41 are expressed as follows:

{\begin{cases} σ_{11} \\ σ_{22} \\ σ_{33} \\ σ_{12} \\ σ_{23} \\ σ_{31} \end{cases}} = [\begin{array}{l} C_{11}^{0} & C_{12}^{0} & C_{13}^{0} & 0 & 0 & 0 \\ C_{12}^{0} & C_{22}^{0} & C_{23}^{0} & 0 & 0 & 0 \\ C_{13}^{0} & C_{23}^{0} & C_{33}^{0} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44}^{0} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55}^{0} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66}^{0} \end{array}] {\begin{cases} ε_{11} \\ ε_{22} \\ ε_{33} \\ ε_{12} \\ ε_{23} \\ ε_{31} \end{cases}}

(1)

where:

\begin{array}{l} C_{11}^{0} = E_{1} (1 - ν_{23} ν_{32}) Δ; C_{12}^{0} = E_{1} (ν_{12} + ν_{31} ν_{23}) Δ; \\ C_{13}^{0} = E_{3} (ν_{31} + ν_{21} ν_{32}) Δ; C_{22}^{0} = E_{2} (1 - ν_{13} ν_{31}) Δ; \\ C_{23}^{0} = E_{2} (ν_{32} + ν_{12} ν_{31}) Δ; C_{33}^{0} = E_{3} (1 - ν_{12} ν_{21}) Δ; \\ C_{44}^{0} = G_{12}; C_{55}^{0} = G_{23}; C_{66}^{0} = G_{13} \end{array}

(2)

Δ = 1 / (1 - ν_{12} ν_{21} - ν_{23} ν_{32} - ν_{13} ν_{31} - 2 ν_{13} ν_{21} ν_{32})

(3)

where

E_{i}

represents the Young’s modulus in the

i

-direction;

ν_{ij}

denotes the Poisson’s ratio; and

G_{ij}

is the shear modulus.

After an element undergoes damage failure, material property degradation and stiffness reduction are required. Consequently, the stiffness matrix element $C_{i j}$ varies with the damage state $C_{i j}^{0}$ . Meanwhile, damage state variables are introduced to simulate the continuous damage process of the model more accurately⁴²; the specific form is as follows:

{\begin{cases} σ_{11} \\ σ_{22} \\ σ_{33} \\ σ_{12} \\ σ_{23} \\ σ_{31} \end{cases}} = [\begin{array}{l} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{13} & C_{23} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 C_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & 2 C_{66} \end{array}] {\begin{cases} ε_{11} \\ ε_{22} \\ ε_{33} \\ ε_{12} \\ ε_{23} \\ ε_{31} \end{cases}}

(4)

\begin{array}{l} C_{11} = (1 - d_{f}) C_{11}^{0}; C_{12} = (1 - d_{f}) (1 - d_{m}) C_{12}^{0}; C_{13} = (1 - d_{f}) (1 - d_{m}) C_{13}; \\ C_{22} = (1 - d_{f}) (1 - d_{m}) C_{22}^{0}; C_{23} = (1 - d_{f}) (1 - d_{m}) C_{23}^{0}; \\ C_{33} = (1 - d_{f}) (1 - d_{m}) C_{33}^{0}; C_{44} = (1 - d_{f}) (1 - s_{m t} d_{m t}) (1 - s_{m c} d_{m c}) C_{44}^{0}; \\ C_{55} = (1 - d_{f}) (1 - s_{m t} d_{m t}) (1 - s_{m c} d_{m c}) C_{55}^{0}; C_{66} = (1 - d_{f}) (1 - s_{m t} d_{m t}) (1 - s_{m c} d_{m c}) C_{66}^{0} \end{array}

(5)

d_{f} = 1 - (1 - d_{f t}) (1 - d_{f c})

(6)

d_{m} = 1 - (1 - d_{m t}) (1 - d_{m c})

(7)

Where

d_{f}

represents the global damage variable for the fiber, and

d_{m}

represents the global damage variable for the matrix.

d_{s} = 1 - (1 - d_{f t}) (1 - d_{f c}) (1 - d_{m t}) (1 - d_{m c})

(8)

In the above expression,

d_{f t}

and

d_{f c}

represent the fiber tensile and compressive damage variables, respectively, while

d_{m t}

and

d_{m c}

correspond to the matrix tensile and compressive damage variables, respectively. These four variables in equation (8) are determined by equation (13) in the following section.

Damage initiation

To determine the onset of material damage, the 3D-Hashin damage criterion^43,44 is introduced. This criterion accounts for the failure modes and the combined influence of individual stress components under interactive stress conditions. Its specific expressions are as follows:

Fiber tensile mode ( $σ_{11} \geq 0$ ):

F_{f t} = {(\frac{σ_{11}}{X_{t}})}^{2} + {(\frac{σ_{12}}{S_{12}})}^{2} + {(\frac{σ_{31}}{S_{13}})}^{2}

(9)

Fiber compression mode ( $σ_{11} < 0$ ):

F_{f c} = \frac{| σ_{11} |}{X_{c}}

(10)

Matrix tensile mode ( $σ_{22} + σ_{33} \geq 0$ ):

F_{m t} = \frac{1}{Y_{t}^{2}} {(σ_{22} + σ_{33})}^{2} + \frac{1}{S_{23}^{2}} (σ_{23}^{2} - σ_{22} σ_{33}) + \frac{1}{S_{12}^{2}} (σ_{12}^{2} + σ_{31}^{2})

(11)

Matrix compression mode ( $σ_{22} + σ_{33} < 0$ ):

F_{m c} = \frac{1}{Y_{c}} [{(\frac{Y_{c}}{2 S_{23}})}^{2} - 1] (σ_{23} + σ_{33}) + \frac{1}{4 S_{23}^{2}} {(σ_{23} + σ_{33})}^{2} + \frac{1}{S_{23}^{2}} (σ_{23}^{2} - σ_{22} σ_{33}) + \frac{1}{S_{12}^{2}} (σ_{12}^{2} + σ_{31}^{2})

(12)

Where

F_{f t}, F_{f c}, F_{m t}, F_{m c}

represent the damage parameters for different failure modes;

S_{12}, S_{13}, S_{23}

are the shear strengths on the 12, 13, and 23 planes, respectively;

X_{t}

and

X_{c}

denote the tensile strength and compressive strength in the fiber direction;

Y_{t}

and

Y_{c}

are the tensile strength and compressive strength in the matrix direction, respectively.

Damage evolution

To characterize the progressive degradation of material stiffness, this paper adopts a damage evolution law based on the equivalent displacement damage coefficient.⁴⁵ In this law, the evolution of the damage variable from 0 (damage-free) to 1 (complete failure) achieves a continuous quantification of damage accumulation. By correlating damage evolution with the material fracture energy, this method ensures the validity of energy dissipation. Consequently, it effectively overcomes the mesh dependence problem in finite element simulations, rendering the computational results more objective and stable. The specific formulas are as follows:

d_{i} = \frac{δ_{e q}^{f, i} (δ_{e q}^{i} - δ_{e q}^{0, i})}{δ_{e q}^{i} (δ_{e q}^{f, i} - δ_{e q}^{0, i})}, (i = f_{t}, f_{c}, m t, m c)

(13)

Where

δ_{e q}^{0, i}

represents the equivalent strain at the onset of damage, and

δ_{e q}^{f, i}

represents the equivalent strain at complete failure. The superscripts denote four distinct failure modes: fiber tension (

f_{t}

), fiber compression (

f_{c}

), matrix tension (

m t

), and matrix compression (

m c

). Furthermore, the evolution of any local damage variable leads to strain softening along the path indicated by line segment A-B in Figure 2. Strain softening introduces strong mesh dependency due to strain localization, meaning the dissipated energy decreases as the mesh is refined. This occurs because the energy consumed during damage evolution is proportional to the volume of the failing element, not to the area of the fracture surface. Consequently, this dissipated energy approaches zero as the mesh is continuously refined. The calculation method for linear stiffness degradation⁴⁶ is shown in equation (14) through equation (17). Meanwhile, the equivalent strain under three-dimensional stress states in equation (13) is also calculated based on equation (14) through equation (17).

Figure 2.

Linear stiffness degradation curve.

Fiber tensile mode:

δ_{e q}^{f t} = \sqrt{{〈 e_{11} 〉}^{2} + e_{12}^{2} + e_{13}^{2}}; δ_{e q}^{o, f t} = \frac{X_{t}}{E_{1}}; δ_{e q}^{f, f t} = \frac{2 G_{c}^{f t}}{X_{t} L_{c}}

(14)

Fiber compression mode:

δ_{e q}^{f c} = 〈 - e_{11} 〉; δ_{e q}^{o, f c} = \frac{X_{C}}{E_{1}}; δ_{e q}^{f, f c} = \frac{2 G_{c}^{f c}}{X_{C} L_{c}}

(15)

Matrix tensile mode:

δ_{e q}^{m t} = \sqrt{{〈 e_{22} 〉}^{2} + {〈 e_{33} 〉}^{2} + e_{12}^{2} + e_{13}^{2}}; δ_{e q}^{o, m t} = \frac{Y_{t}}{E_{2}}; δ_{e q}^{f, m t} = \frac{2 G_{c}^{m t}}{Y_{t} L_{c}}

(16)

Matrix compression mode:

δ_{e q}^{m c} = \sqrt{{〈 - e_{22} 〉}^{2} + {〈 - e_{33} 〉}^{2} + e_{12}^{2} + e_{13}^{2}}; δ_{e q}^{o, m c} = \frac{Y_{c}}{E_{2}}; δ_{e q}^{f, m c} = \frac{2 G_{c}^{m c}}{Y_{c} L_{c}}

(17)

where

〈 〉

is an operator symbol, meaning

〈 x 〉 = (x + \frac{| x |}{2})

Experimental study

Specimen design

The carbon fiber fabric used in this study was supplied by a well-established Chinese manufacturer. The tensile strengths of the carbon fiber fabric and the epoxy adhesive were 4323 MPa and 38 MPa, respectively, and their elastic moduli were 234 GPa and 2.4 GPa, respectively. As shown in Figure 3, five electrical resistance strain gauges were bonded to the specimen surface to measure the strain distribution during the erosion process. To prevent premature failure at the gripping ends, end tabs were bonded to both ends of the carbon fiber fabric specimens to distribute stresses, and steel plates were used to protect against clamping-induced crushing, thereby ensuring that tensile failure occurred within the central region of the epoxy adhesive system. The designed dimensions of the test specimens were 750 mm in length and 75 mm in width. After fabrication, the specimens were cured for 12 h under standard laboratory conditions at 20°C, after which the wind-blown sand erosion tests were conducted.

Figure 3.

Macro-detail view of the specimen.

Test setup

The experiment employed the jet erosion method, which effectively simulates the characteristics of actual wind-sand environments and allows for accurate and efficient control of the erosion angle and sand flow rate, thereby facilitating the successful completion of the erosion tests. The experimental setup is illustrated in Figure 4 and consists of three parts: a high-pressure gas supply system, a sand feeding system, and an erosion chamber. Specifically, the sand feeding system was used to precisely regulate the sand delivery rate per unit time, while the erosion angle was controlled by rotating and positioning the specimen stage. Inside the erosion chamber, the nozzle diameter was set to 15 mm, and the distance from the nozzle to the specimen surface (stand-off distance) was fixed at 150 mm. Under these parameters, the diameter of the erosion-damaged area on the specimen surface closely matched the specimen width. This ensured that the central region of the specimen was subjected to uniform erosion, thereby enhancing the consistency and comparability of the experiments.

Figure 4.

Experimental system of blown-sand erosion test.

The remaining mechanical properties of the specimens subjected to wind-blown sand erosion were evaluated using an electronic universal testing machine. During the tests, the deformation of the CFRP laminate was monitored in real time by combining Digital Image Correlation (DIC) with a strain acquisition system. A high-speed camera was used to continuously capture deformation images, which were processed by the associated analysis system to generate full-field strain contours, thereby enabling precise quantification of both local and global deformation behavior of the material. The specific arrangement of the experimental setup is shown in Figure 5.

Figure 5.

Setup of residual mechanical properties test.

Experimental results of wind-blown sand erosion on CFRP laminate

Mechanical properties of CFRP laminate after wind-blown sand erosion

Taking into account the environmental characteristics such as wind speed and dust concentration during sandstorms in the central and western regions of Inner Mongolia, and based on existing research and preliminary test results, the experimental conditions shown in Table 2 were determined. To ensure the accuracy and statistical reliability of the test data, three parallel duplicate specimens were prepared for each condition.

Table 2.

Test conditions.

Erosion velocity/(m·s^-1)	Erosion time/min	Erosion angle/°	Sand flow rate/(g·min^-1)
16	40	90	55
19
23
26
31

Figure 6 illustrates the evolutionary patterns of the tensile strength and elastic modulus of CFRP laminate following wind-blown sand erosion. It is evident that wind-blown sand erosion significantly weakens the tensile properties of the material. At an erosion angle of 90°, the tensile strength of the CFRP laminate exhibited a clear downward trend with increasing erosion velocity; the maximum strength loss occurred at 31 m/s, with a reduction of approximately 23.16%. Meanwhile, while the elastic modulus fluctuated slightly with increasing velocity, it remained generally stable overall. Notably, the most significant decrease in elastic modulus was observed at 26 m/s, with a reduction of 7.2%. As an intrinsic property, the elastic modulus is primarily determined by the material’s composition and internal structure. Although wind-blown sand erosion, as a surface physical damage process, may induce surface fiber breakage or resin spalling in the CFRP laminate, it does not significantly alter the overall internal structure of the material. Consequently, within the range of erosion parameters examined in this study, the elastic modulus of the CFRP laminate showed no drastic changes, demonstrating excellent stability. In summary, wind-blown sand erosion has a minor impact on the elastic modulus of CFRP laminate, but it significantly degrades their tensile strength.

Figure 6.

Mechanical properties of CFRP laminate after wind-blown sand erosion.

Damage morphology of CFRP laminate after wind-blown sand erosion

As shown in Figure 7, significant damage occurred on the surface of the CFRP laminate after wind-blown sand erosion. At an erosion velocity of 16 m/s, slight damage was visible on the specimen surface, accompanied by localized slight fragmentation of the epoxy resin matrix. As the erosion velocity increased to 19 m/s and 23 m/s, the degree of surface damage intensified significantly, and spalling of the surface resin matrix occurred. When the velocity further increased to 26 m/s, the spalling of the resin matrix aggravated, internal fibers became exposed, and partial fiber bundles fractured. At an erosion velocity of 31 m/s, the surface damage was the most severe, characterized by large-scale spalling of the resin matrix and extensive exposure of internal fibers, which exhibited obvious fracture and damage under the continuous impact of wind-blown sand particles.

Figure 7.

Damage morphology of CFRP laminate after wind-blown sand erosion: (a) 16 m/s; (b) 19 m/s; (c) 23 m/s; (d) 23 m/s; (e) 31 m/s.

Damage mechanism of CFRP laminate under wind-blown sand erosion

The damage to CFRP fabrics caused by wind-blown sand erosion is a continuous process evolving from the accumulation of microscopic defects to the deterioration of macroscopic properties. As shown in Figure 8, this process is primarily driven by the kinetic energy carried by impacting sand particles, which is predominantly converted into localized compressive stress at the impact points. Specifically, high-speed sand particles striking the surface of the CFRP fabric initially induce elastoplastic deformation and denudation of the epoxy resin matrix. As erosion progresses, the surface layer of the matrix gradually wears away, thereby exposing the internal fibers. These exposed fibers are then subjected to direct impact, leading to fracture or fragmentation under erosive action. Ultimately, these microscopic structural damages directly cause macroscopic performance degradation: the decrease in tensile strength primarily stems from fiber fracture, while the reduction in stiffness is attributed to the combined effects of interfacial failure and microcracking. Consequently, the damage of CFRP fabrics by wind-blown sand erosion is essentially a progressive physical destruction process that is driven by kinetic energy, proceeds from the surface to the interior, and extends from the interface to the whole structure.

Figure 8.

Damage mechanism of CFRP laminate under wind-blown sand erosion.

Model validation

Single element verification

The specimen parameters obtained through experiments are used to construct the finite element model to perform numerical prediction and simulation analysis of the experimental process. This model couples the material property parameters via a VUMAT user material subroutine written in Fortran, thereby achieving a more accurate simulation of the tensile mechanical behavior after erosion. The stress update algorithm and the specific calculation flow within the VUMAT subroutine are illustrated in Figure 9. In this study, a dynamic explicit algorithm is employed. Compared to static algorithms, the dynamic explicit algorithm offers higher computational accuracy when dealing with large deformation nonlinear problems and demonstrates significant advantages in the efficient solution of impact-type simulations.

Figure 9.

Flowchart of VUMAT workflow.

To verify the accuracy of the written VUMAT subroutine, and considering that the actual thickness of the CFRP laminate is relatively thin (0.167 mm), it was found during the simulation process that if a single mesh (element) were used for verification, the elements in the thickness direction would be prone to severe distortion. Therefore, this study employs a coarse mesh model of the overall structure for verification, rather than targeting a single mesh element. In the result analysis phase, one of the typical element meshes is selected for comparative verification. To validate the accuracy of the model, a uniaxial tensile simulation along the fiber direction was conducted, and the model was assigned the actual material parameters of the CFRP laminate. Correspondingly, a 3D solid element model was established with clearly defined material orientations, as shown in Figure 10, and the analysis steps, and boundary conditions were set as shown in Figure 11.

Figure 10.

Material orientation.

Figure 11.

Boundary conditions.

A comparison between the simulation results of the material’s mechanical properties and the experimental data is presented in Figure 12. The bar chart indicates that the simulated tensile strength and elastic modulus agree well with the experimental values, with relative errors of 2.7% for tensile strength and 3.8% for elastic modulus. Figure 13 shows the stress-strain response curve obtained from the simulation. At a strain of 0.01,351, the simulated ultimate stress is 3486.2 MPa, which is close to the experimentally measured ultimate stress of 3586.5 MPa, further confirming the reliability of the established model. Figure 14 illustrates the variations in the fiber initial damage coefficient $F_{f t}$ , and the damage evolution coefficient $d_{f t}$ during axial tension. As the strain increases, $F_{f t}$ rises nonlinearly, while $d_{f t}$ remains 0 before damage initiation. When the strain reaches 0.01,351 and the stress is 3486.2 MPa, $F_{f t}$ attains 1.0, at which point $d_{f t}$ begins to increase, rising rapidly to 1.0 and then stabilizing. This indicates that when $F_{f t}$ reaches 1.0, the material enters the linear damage evolution stage; when $d_{f t}$ reaches 1.0, the corresponding material element is fully damaged, leading to a decline in the overall load-bearing capacity and rapid failure. The damage evolution process simulated by this model is consistent with the characteristics of CFRP as a brittle material; that is, damage develops over a very short period, with strength degrading linearly until final failure.

Figure 12.

Comparison of the experimental results with the simulation results.

Figure 13.

Stress-strain response curve.

Figure 14.

Damage variable - strain response curve.

Mesh sensitivity analysis

To ensure the reliability and accuracy of numerical simulation results, it is crucial to reduce the mesh dependency of the model and guarantee its convergence. This study evaluates the convergence behavior of the model through a mesh sensitivity analysis, specifically by employing a method of gradually reducing the mesh size. In the CFRP laminate model, reference planes are used for partitioning, and mesh refinement is applied to the erosion zone. Four sets of mesh sizes (5 mm, 3 mm, 1 mm, and 0.5 mm) were defined for the sensitivity analysis of this region. As shown in Figure 15, the axial tensile simulation results indicate that the load-displacement curves for different mesh sizes overlap significantly before reaching the peak load, suggesting that the response in this stage is essentially independent of the mesh; however, the post-peak stage exhibits slight mesh dependency. Considering computational efficiency and convergence stability comprehensively, the model demonstrates good convergence and stable results with a 1 mm mesh size. Consequently, the mesh size for the erosion zone is determined to be 1 mm. To balance computational costs, the non-erosion region utilizes a bias-based gradient mesh where the element size gradually increases from 1 mm to 15 mm. This mesh strategy significantly improves computational efficiency while ensuring computational accuracy and convergence.

Figure 15.

Load-displacement responses of elements of different sizes under tensile simulation.

Comparison of mechanical properties after erosion

In this section, the mechanical property degradation of the CFRP laminate model under discrete particle erosion was numerically simulated and compared with experimental results, as shown in Figure 16. As shown in Figure 16(a), with the continuous increase of particle impact velocity, the simulated tensile strength of the CFRP laminate gradually decreases. This degradation trend is consistent with the law observed in the experiments where strength declines with rising velocity. The error range between the simulated and experimental values under various velocity conditions is between 1.37% and 4.18%, indicating that the model possesses high prediction accuracy. Figure 16(b) illustrates the comparison between the simulated and experimental values of the elastic modulus of the CFRP laminate after erosion. The simulation results indicate that the elastic modulus does not change significantly with the increase of erosion velocity, which is basically consistent with the experimental conclusions. Meanwhile, the simulated values are close to the experimental ones, with an error range of 3.56%–6.48%. It is worth noting that, as seen from the bar chart, the experimental values exhibit a slight decrease at 23 m/s, 26 m/s, and 31 m/s, whereas the simulated values remain stable. This discrepancy may stem from the relatively ideal simulation environment, whereas in actual experiments, factors such as fiber dispersion and manufacturing errors can induce slight fluctuations in material properties, thereby resulting in a slight reduction in the experimental values of the elastic modulus.

Figure 16.

Mechanical properties of CFRP laminate after wind-blown sand erosion: (a) Tensile strength; (b) Elastic modulus.

In summary, the simulation results of the model agree well with the experimental results, and the errors fall within an acceptable range for engineering applications. Therefore, the model can be used for further analysis and prediction of the mechanical properties of CFRP laminate under Wind-Blown Sand Erosion conditions.

Damage morphology after particle impact

Under particle erosion at different initial velocities (16 m/s, 19 m/s, 23 m/s, 26 m/s, and 31 m/s), the surfaces of CFRP laminate specimens exhibited varying degrees of damage morphologies. Similarly, in numerical simulation, the material surface is primarily subjected to compressive stress from the particles upon impact, a stress shared by both the fiber bundles and the matrix. Therefore, this study adopts the fiber compression initial damage coefficient ( $F_{f c}$ ), the fiber compression damage evolution coefficient ( $d_{f c}$ ), the matrix compression initial damage coefficient ( $F_{m c}$ ), and the matrix damage evolution coefficient ( $d_{m c}$ ) to determine the initiation and evolution of damage during the impact process. When $F_{f c}$ or $F_{m c}$ equals 1, damage is considered to have initiated in the material; if $d_{f c}$ or $d_{m c}$ reaches or exceeds 1, the corresponding material element is deemed to have failed and is subsequently deleted from the model.

At the relatively low impact velocity of 16 m/s, the kinetic energy of the particles is relatively low, yet the initial damage coefficients $F_{f c}$ and $F_{m c}$ have already reached 1. As shown in the contours of SDV19 and SDV21 in Figure 17, elements with colors tending towards red have values greater than or equal to 1.0, indicating that these regions have entered the damage initiation stage. However, it is clearly evident from the contours of SDV2 and SDV4 that although the values of the damage evolution coefficients $d_{f c}$ and $d_{m c}$ approach 1.0, neither reaches the failure threshold of 1.0. This indicates that only micro-cracks have developed in the material without element failure; consequently, no element deletion is observed on the model surface. At this stage, the impact of erosion on the macroscopic mechanical properties of the material is relatively limited. Although micro-cracks have been generated, the material retains its load-bearing capacity, and the stress in the fiber direction shows no significant decline. This simulation result is consistent with the observed phenomenon of minor material strength loss under low-velocity conditions in Wind-Blown Sand Erosion experiments.

Figure 17.

Erosion damage diagram of 16 m/s.

When the initial particle velocity is increased to 19 m/s and 23 m/s, the damage behavior of the model undergoes significant changes. As shown in Figures 18 and 19, the contours of SDV19 and SDV21 indicate a marked increase in the number of elements whose damage coefficients have reached the threshold compared to the 16 m/s case. Meanwhile, the number of elements undergoing strength degradation in the damage evolution coefficient contours (SDV2 and SDV4) also rises significantly. Distinct holes can be observed on the model surface, indicating that some elements have reached their ultimate strength and have been deleted. It is worth noting that although elements approaching 1.0 are observed in the SDV2 and SDV4 contours, values equal to or greater than 1.0 are not visible. This is because once the damage evolution coefficient reaches 1.0, the corresponding element is immediately deleted and therefore no longer appears in the contours. The stress contour in the fiber direction further reveals a significant increase in dark blue spots, reflecting that internal damage within the impacted elements has caused them to lose their tensile load-bearing capacity. Furthermore, the appearance of dark red spots in the figures indicates the occurrence of local stress concentration due to element deletion. As the external loading continues, these stress concentration zones will undergo fracture first, ultimately triggering the failure of the overall material.

Figure 18.

Erosion damage diagram of 19 m/s.

Figure 19.

Erosion damage diagram of 23 m/s.

Figure 20 illustrates a significant increase in element deletion on the specimen surface when the erosion velocity is raised to 26 m/s. The damage coefficient contour SDV19 reveals that the damage coefficients of a large number of elements have exceeded 1.0, indicating significant damage to the fibers under impact loading. Meanwhile, the SDV21 contour reflects that the matrix damage areas are relatively concentrated; however, no elements with values greater than or equal to 1 appear in this contour. Combined with the SDV4 damage evolution coefficient contour, it can be inferred that elements in SDV21 that originally exceeded 1.0 have been deleted due to reaching the ultimate strength, and thus are no longer displayed in the SDV21 contour. Observing the damaged surface of the model, it is evident that the number of deleted elements has increased significantly. The overall proliferation of blue elements, which represent the loss of load-bearing capacity, indicates that the material damage has become severe, ultimately leading to a decrease in the ultimate strength of the overall model.

Figure 20.

Erosion damage diagram of 26 m/s.

As shown in Figure 21, when the velocity is further increased to 31 m/s, the number of deleted elements on the specimen surface rises sharply. The damage coefficient contour SDV19 indicates significant fiber damage, while the SDV21 contour shows that the matrix damage coefficients generally exceed 1.0, with more concentrated damage zones. The corresponding damage evolution coefficient contours (SDV4) also exhibit a dense distribution of high-value elements, reflecting a highly concentrated spatial pattern of damage evolution. Multiple large-area holes are visible in the model surface morphology contour, indicating widespread generation of failed elements, and impact damage has extended significantly from the surface into the interior. At this stage, the CFRP laminate model is in a severely damaged state, with markedly degraded macroscopic mechanical properties. This simulated morphology is consistent with the severe damage characteristics of the material observed in high-speed Wind-Blown Sand Erosion experiments.

Figure 21.

Erosion damage diagram of 31 m/s.

Prediction of mechanical properties

Case design

Based on the comparison between experimental data and simulation results, the accuracy of the model in predicting the mechanical properties and damage morphologies of CFRP laminate under various erosion velocities has been validated. To further evaluate the model’s predictive capability in real-world wind-sand environments, this study simulated the mechanical response of CFRP laminate under different sand flow rates by varying the number of particles in the model. Specifically, the conversion between the DEM particle count and the sand flow rate is strictly based on the mass conservation principle. The corresponding simulation cases are summarized in Table 3

Table 3.

Case design.

Erosion velocity/(m/s)	Erosion angle/°	Number of particles
26	90	100
		200
		300
		400
		500

Mechanical properties of CFRP laminate model post-erosion

Figure 22 illustrates the evolution of the mechanical properties of CFRP laminate under erosion conditions with different numbers of particles, as predicted by the ABAQUS model. It can be observed that when the number of particles is 100, 200, and 300, the predicted tensile strength of the CFRP laminate model does not exhibit a significant decrease. This may be attributed to the relatively small number of particles; the cumulative energy upon impact on the CFRP laminate model is low and insufficient to cause significant damage to the surface elements, thus the overall mechanical properties remain relatively stable. As the number of particles increases to 400, the tensile strength of the CFRP laminate exhibits a significant decline, decreasing by approximately 16.71% compared to the 300-particle case. At this stage, the total energy exerted by particle impact on the CFRP laminate model is substantial, causing noticeable damage to the material. When the number of particles further increases to 500, the tensile strength continues to decrease, dropping by approximately 29.26% relative to the 300-particle case. The aforementioned results indicate that as the number of particles increases, the damage sustained by the CFRP laminate model intensifies. Consequently, the number of elements deleted due to reaching the ultimate strength increases correspondingly, leading to a distinct deterioration trend in the tensile strength.

Figure 22.

Mechanical properties of CFRP laminate model after erosion with different numbers of particles.

It is worth noting that the elastic modulus of the CFRP laminate model does not exhibit significant fluctuations under any of the particle number cases, remaining relatively stable throughout. This indicates that, although particle erosion in the predictive model leads to a decrease in macroscopic strength, it does not alter the fundamental elastic properties of the material.

Damage morphology of CFRP laminate model post-erosion

The erosion damage morphology on the surface of the CFRP laminate model under different numbers of particles is shown in the figure. From the damage morphology images in Figure 23(a)–(c), it can be seen that at an initial velocity of 26 m/s and with particle numbers of 100, 200, and 300, no obvious damage appears on the surface of the CFRP laminate model. The dark blue spots in the damage contour indicate that the elements at these locations have met the damage initiation criterion; microcracks may have initiated or softening may have occurred internally, but since they have not yet reached the ultimate strength, these elements are retained and can continue to bear tensile stress. As the number of particles increases to 400, voids begin to appear on the model surface, as shown in Figure 23(d). This is because the increased particle number leads to a rise in cumulative impact energy; some surface elements reach their strength limit with a damage evolution coefficient greater than or equal to 1.0, and are therefore deleted. These deleted elements no longer participate in load carrying, resulting in a decrease in the ultimate tensile strength of the model. When the number of particles further increases to 500, a large number of dense voids appear on the surface of the CFRP laminate model. The increase in particle number significantly aggravates material damage; more surface elements are deleted due to reaching the failure condition, thereby severely weakening the overall tensile strength of the model.

Figure 23.

Surface damage morphology of CFRP laminate model after erosion with different numbers of particles: (a) 100 particles; (b) 200 particles; (c) 300 particles; (d) 400 particles; (e) 500 particles.

Analysis of prediction accuracy

To verify the accuracy of the model predictions, corresponding supplementary experiments were conducted. By controlling the sand flow rate in the supplementary tests, the effects of different sand concentrations under extreme weather conditions on the erosion damage of CFRP laminate were simulated. The correspondence between the sand flow rate and the number of particles in the supplementary experimental conditions is shown in Table 4.

Table 4.

Correspondence between supplementary experimental conditions and particle numbers.

Erosion velocity/(m·s^-1)	Erosion time/min	Erosion angle/°	Sand flow rate/(g·min^-1)	Number of particles
26	40	90	15	100
			25	200
			35	300
			45	400
			55	500

Through a comparative analysis of the mechanical property test results of CFRP laminate after erosion at different sand flow rates and the model prediction results, the following conclusions can be drawn. As shown in Table 5, when the sand flow rate is low (15 g/min, 25 g/min, and 35 g/min), the tensile strength of the material does not exhibit significant changes, which is consistent with the model predictions for a smaller number of particles (100, 200, and 300). When the sand flow rate increases to 45 g/min, the tensile strength of the CFRP laminate decreases significantly, aligning with the predicted trend when the number of particles increases to 400. As the sand flow rate continues to increase, the tensile strength further decreases, a trend that is also consistent with the model predictions. The error between the model-predicted and experimentally measured tensile strengths ranges from 0.83% to 6.53%, indicating that the model possesses high accuracy in predicting strength degradation. Meanwhile, as shown in

Table 5.

Tensile strength comparison.

Sand flow rate/(g·min^-1)	Tensile strength/MPa	Number of particles	Tensile strength/MPa	Relative error/%
15	3669.5	100	3638.68	0.83
25	3783.5	200	3536.19	6.53
35	3764	300	3529.26	6.23
45	2914.5	400	2939.23	0.84
55	2647.5	500	2483.63	6.18

Table 6, with the increase in sand flow rate and particle number, neither the predicted nor experimental values of the elastic modulus show significant fluctuations. The error between the two lies between 0.44% and 1.4%, further verifying the reliability of the model. This also demonstrates that both the wind-blown sand erosion test and discrete particle impact are physical damage processes and do not alter the material’s intrinsic constitutive properties. In summary, through the comparison with experimental results in terms of tensile strength and elastic modulus, this model demonstrates high prediction accuracy and can be used to simulate the evolution of the mechanical properties of CFRP laminate in real wind-sand environments.

Table 6.

Elastic modulus comparison.

Sand flow rate/(g·min^-1)	Elastic modulus/GPa	Number of particles	Elastic modulus/GPa	Relative error/%
15	283.6	100	280.5	1.09
25	284.6	200	286.1	0.52
35	285.7	300	283.4	0.78
45	271.5	400	272.7	0.44
55	266.7	500	270.5	1.40

Figure 24 illustrates the evolution of the surface morphology of CFRP laminate after wind-blown sand erosion at different sand flow rates. At lower sand flow rates (15, 25, and 35 g/min), no significant macroscopic damage was observed on the material surface, with only slight spalling of the epoxy resin matrix (Figure 24(a)–(c)). When the sand flow rate increased to 45 g/min, the degree of damage significantly intensified, characterized by large-area spalling of the resin matrix, fiber exposure, and extensive fracture (Figure 24(d)). Upon further increasing the sand flow rate to 55 g/min, distinct erosion pits formed on the surface, and the damage further propagated (Figure 24(e)). The damage modes observed in the aforementioned experiments are highly consistent with the model predictions shown in Figure 11, both demonstrating a progressive damage process ranging from surface resin abrasion to fiber fracture and finally to macroscopic pitting. In summary, the erosion model established in this study exhibits high prediction accuracy and can reliably simulate the degradation of mechanical properties and the damage evolution process of CFRP laminate under real wind-sand environments, thereby providing an analytical basis for the durability assessment of CFRR structures in such conditions.

Figure 24.

Wind-blown sand erosion damage morphology of CFRP laminate under different sand flow rates: (a) 15 g/min; (b) 25 g/min; (c) 35 g/min; (d)45 g/min; (e) 55 g/min.

It should be objectively noted that, constrained by the ultra-thin thickness of the CFRP laminate (0.167 mm) and the size effect induced by mesh discretization in the numerical model, the comparison of damage morphology in this study primarily focuses on the qualitative agreement of evolutionary patterns and trends, rather than directly extracting damage area or depth for quantitative measurement. However, the high consistency of the macroscopic tensile strength degradation data (with an error of less than 5%) indirectly corroborates the reliability of the model in assessing the extent of damage evolution from a mechanical response perspective. This provides robust support for the qualitative morphological comparison, thereby validating the accuracy of the proposed model based on damage morphological features.

Conclusion

In this study, the erosion process of CFRP laminate by wind-sand particles was simulated using ABAQUS software, where the dynamic erosion of particles was achieved via a built-in particle generator. A VUMAT subroutine was developed based on the 3D Hashin failure criterion to characterize the material’s damage initiation, evolution, and stiffness degradation behavior during the erosion process. The subroutine was validated at the element scale, and a mesh sensitivity analysis was performed on the computational model. On this basis, the evolution of damage morphology of CFRP laminate under different erosion velocities was systematically analyzed, and the simulation results were verified by comparison with experimentally measured variations in mechanical properties. Finally, the model was utilized to predict the mechanical property degradation and damage morphology characteristics of CFRP laminate under different sand flow rates. The conclusions are drawn as follows:

(1) In the wind-sand erosion tests, the tensile strength of the CFRP laminate exhibited a significant downward trend with increasing erosion velocity. At an erosion velocity of 31 m/s, the tensile strength experienced its maximum reduction of approximately 23.16%. However, the increased erosion velocity had no noticeable impact on the elastic modulus of the CFRP laminate. The pronounced degradation in tensile strength is attributed to the fracture failure of locally load-bearing fibers driven by the kinetic energy of the sand particles. In contrast, the stability of the elastic modulus is primarily due to the fact that wind-sand erosion is essentially a surface physical damage mechanism; it mainly induces local microstructural destruction and strength degradation, while minimally affecting the inherent constitutive properties that determine the overall stiffness of the material.

(2) To verify the accuracy of the developed VUMAT subroutine, this study employed coarse mesh modeling and selected a representative element for validation. The simulated tensile strength and elastic modulus agreed well with the experimental values, with relative errors of 2.7% and 3.8%, respectively. The results indicate that the subroutine can reliably characterize the mechanical behavior of CFRP laminate under tensile loads.

(3) Based on the ABAQUS platform, this study utilized its built-in particle generator to simulate the erosion process of sand particles on CFRP laminate. The results indicate that the degree of damage in the model significantly intensifies as the initial velocity of the particles increases. When the velocity reaches 31 m/s, the number of surface elements deleted due to reaching the damage limit is the highest, resulting in the most significant impact on the macroscopic mechanical properties of the material. The simulated damage morphology is consistent with experimental observations, and the prediction errors for tensile strength and elastic modulus are maintained within 1.37%–4.18% and 3.56%–6.48%, respectively. This verifies the feasibility and reliability of the simulation method for investigating wind-blown sand erosion behavior.

(4) A wind-blown sand erosion model was established on the ABAQUS platform, utilizing the built-in particle generator to predict the degradation of mechanical properties of CFRP laminate under different sand flow rate conditions. Supplementary experiments were conducted to compare the experimental data with the model predictions. The results indicated that the tensile strength of the CFRP laminate exhibited a consistent downward trend in both the experiments and predictions, with an error range of 0.83%–6.53% between them. The elastic modulus remained stable in both cases, with errors ranging from 0.44% to 1.4%. The comparison of damage morphologies further revealed that the surface damage evolution patterns observed in the experiments were highly consistent with the simulation predictions.

(5) Although the proposed model demonstrates sound reliability in predicting the macroscopic mechanical degradation trends of CFRP laminate, due to computational efficiency constraints, it incorporates several idealized assumptions. These include neglecting the irregular morphologies, impact fragmentation, and inter-particle interactions of sand particles, as well as modeling the CFRP laminate as a macroscopic homogeneous medium while overlooking the inherent randomness of the fiber materials. Furthermore, the current model focuses solely on purely mechanical kinetic impacts, without incorporating the multi-physics effects present in real-world wind-blown sand environments. Consequently, the proposed method is presently more applicable for evaluating the performance degradation and damage evolution trends of CFRP laminate under idealized erosion conditions. Future research could expand its applicability by introducing mesoscale characterizations and multi-physics coupled models.

Footnotes

Acknowledgments

The authors are deeply grateful for the financial support from National Natural Science Foundation of China (Grant no. 52508560), Natural Science Foundation of Inner Mongolia Autonomous Region of China (Grant No. 2024QN05010), Research Program of Science and Technology at Universities of Inner Mongolia Autonomous Region (Grant No. NJZZ23077), the Basic Scientific Research Expenses Program of Universities directly under Inner Mongolia Autonomous Region (Grant No. JY20250089).

ORCID iD

Jinpeng Zhang

Author contributions

Jinpeng Zhang: Methodology, writing-original draft; Siha A: Supervision, conceptualization; funding acquisition; validation; writing-review and editing; Shuai Hao: Software and writing-review.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The financial support from National Natural Science Foundation of China (Grant no. 52508560), Natural Science Foundation of Inner Mongolia Autonomous Region of China (Grant No. 2024QN05010), Research Program of Science and Technology at Universities of Inner Mongolia Autonomous Region (Grant No. NJZZ23077), the Basic Scientific Research Expenses Program of Universities directly under Inner Mongolia Autonomous Region (Grant No. JY20250089).

Declaration of conflicting interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data Availability Statement

Data will be made available on request. .⁴⁷

References

Solahuddin

Yahaya

. A narrative review on strengthening of reinforced concrete beams using carbon fibre reinforced polymer composite material through experimental investigation and numerical modelling. Structures 2023; 52: 666–710. https://doi.org/10.1016/j.istruc.2023.03.168

Frhaan

WKM

Abu Bakar

Hilal

, et al. CFRP for strengthening and repairing reinforced concrete: a review. Innov Infrastruct Solut 2021; 6(2): 49. https://doi.org/10.1007/s41062-020-00417-5

Gietzen

Karimi

Goel

, et al. Experimental investigation of low velocity and high temperature solid particle impact erosion wear. Wear 2022; 506-507: 204441. https://doi.org/10.1016/j.wear.2022.204441

Zhao

Zhu

, et al. Erosion mechanism and damage behavior of CFRP in plastic abrasive jet machining. Int J Adv Manuf Technol 2022; 121(7): 4905–4918. https://doi.org/10.1007/s00170-022-09667-2

Zhang

Liu

, et al. Erosion mechanism and mathematical model of CFRPs in wind-grit flow environments. Constr Build Mater 2025; 458: 139489. https://doi.org/10.1016/j.conbuildmat.2024.139489

Chen

Yin

Z-Y

, et al. A micromechanical-based study on the tribological and creep-relaxation behavior of sand-FRP composite interfaces. Compos Struct 2021; 275: 114423. https://doi.org/10.1016/j.compstruct.2021.114423

Tewari

Harsha

Häger

, et al. Solid particle erosion of carbon Fibre– and glass fibre–epoxy composites. Compos Sci Technol 2003; 63(3): 549–557. https://doi.org/10.1016/s0266-3538(02)00210-5

Wang

Cao

, et al. Damage mechanism and residual tensile strength of CFRP laminates subjected to high-velocity sand erosion. Compos B Eng 2025; 300: 112459. https://doi.org/10.1016/j.compositesb.2025.112459

Zhang

Liu

, et al. Coupled modeling and numerical simulation of gas flows laden with solid particles in de Laval nozzles. Shock Waves 2022; 32(2): 213–230. https://doi.org/10.1007/s00193-021-01063-1

10.

Yang

, et al. Impact of Laval nozzle structure on the flow characteristics of supersonic gas-solid two-phase flow. Powder Technol 2024; 439: 119657. https://doi.org/10.1016/j.powtec.2024.119657

11.

Choudhary

Sharma

Agarwal

, et al. Optimization of solid particle erosion behaviour of waste marble dust filled glass fiber polymer composite using Taguchi approach. Mater Today Proc 2021; 44: 4908–4912. https://doi.org/10.1016/j.matpr.2020.11.884

12.

Desale

Gandhi

Jain

. Effect of erodent properties on erosion wear of ductile type materials. Wear 2006; 261(7): 914–921. https://doi.org/10.1016/j.wear.2006.01.035

13.

Harsha

Thakre

. Investigation on solid particle erosion behaviour of polyetherimide and its composites. Wear 2007; 262(7): 807–818. https://doi.org/10.1016/j.wear.2006.08.012

14.

Miyazaki

. Solid particle erosion of composite materials: a critical review. J Compos Mater 2016; 50(13): 3175–3217. https://doi.org/10.1177/0021998315617818

15.

Chellaganesh

Adam Khan

Mohamed

, et al. Hybrid composite material and solid particle erosion studies. IOP Conf Ser Mater Sci Eng 2018; 346(1): 012014. https://doi.org/10.1088/1757-899x/346/1/012014

16.

Arnold

Hutchings

. Flux rate effects in the erosive wear of elastomers. J Mater Sci 1989; 24(3): 833–839. https://doi.org/10.1007/bf01148765

17.

Shipway

Hutchings

. A method for optimizing the particle flux in erosion testing with a gas-blast apparatus. Wear 1994; 174(1): 169–175. https://doi.org/10.1016/0043-1648(94)90099-x

18.

Liang

Yuan

Wang

, et al. Current situation and prospect of erosion wear. J Phys Conf 2020; 1600(1): 012015. https://doi.org/10.1088/1742-6596/1600/1/012015

19.

Zhang

Deguo

Weihua

. Investigation of abrasive erosion of polymers. J Mater Sci 1995; 30(18): 4561–4566. https://doi.org/10.1007/bf01153063

20.

Zhu

, et al. The effect of interfacial imperfection on the nonlinear dynamic behavior of composite laminates under low-velocity impact. Nonlinear Dyn 2024; 112(14): 11869–11886. https://doi.org/10.1007/s11071-024-09702-9

21.

Zheng

, et al. Two-scale damage failure analysis of 3D braided composites considering pore defects. Compos Struct 2021; 260: 113556. https://doi.org/10.1016/j.compstruct.2021.113556

22.

Zhao

Zhu

, et al. Development of erosion equations for removal of organic coating on carbon fiber reinforced polymer surface by plastic abrasive. J Coating Technol Res 2023; 20(2): 671–685.

23.

Hao

Liu

, et al. Energy consumption and erosion mechanism of polyester fiber reinforced cement composite in wind-blown sand environments. J Wuhan Univ Technol -Materials Sci Ed 2022; 37(4): 666–676. https://doi.org/10.1007/s11595-022-2581-4

24.

Srivastava

Pawar

. Solid particle erosion of glass fibre reinforced flyash filled epoxy resin composites. Compos Sci Technol 2006; 66(15): 3021–3028. https://doi.org/10.1016/j.compscitech.2006.02.004

25.

Zhang

Feng

, et al. Concrete erosion mechanism and model research under the action of wind gravel flow in Gobi gale area. Constr Build Mater 2023; 409: 133940. https://doi.org/10.1016/j.conbuildmat.2023.133940

26.

Jaensson

Hulsen

Anderson

. A comparison between the XFEM and a boundary-fitted mesh method for the simulation of rigid particles in cahn–hilliard fluids. Computers & Fluids 2017; 148: 121–136. https://doi.org/10.1016/j.compfluid.2017.02.003

27.

Leblanc

Boyer

Joslin

. Modelling and animation of impact and damage with smoothed particle hydrodynamics. Vis Comput 2014; 30(6): 909–917. https://doi.org/10.1007/s00371-014-0981-y

28.

Nekahi

Villasenor Vazquez

Papini

. Numerical modeling of the erosion of alumina particulate reinforced epoxy-matrix composites: material removal mechanisms including reinforcement fracture. Wear 2023; 522: 204710. https://doi.org/10.1016/j.wear.2023.204710

29.

Liu

Liao

Jia

, et al. Finite element analysis of dynamic progressive failure of carbon fiber composite laminates under low velocity impact. Compos Struct 2016; 149: 408–422. https://doi.org/10.1016/j.compstruct.2016.04.012

30.

Petrolo

Tortorelli

Saputo

, et al. Tensile damage analysis of composite structures based on 3D hashin failure criteria and 2D higher-order structural theories. Compos Commun 2025; 56: 2452. https://doi.org/10.1016/j.coco.2025.102395

31.

Zouggar

Guerraiche

Rabouh

, et al. Numerical damage assessment in T700/epoxy composite laminate under low and high velocity impacts using a modified hashin-puck criterion. Compos Struct 2025; 372: 0263–8223. https://doi.org/10.1016/j.compstruct.2025.119578

32.

Sha

Xiao

, et al. Enhancing compressive failure prediction in repaired composites: a comparative study of 2D and 3D hashin criteria with experimental validation. Thin-Walled Struct 2025; 217: 0263–8231. https://doi.org/10.1016/j.tws.2025.113911

33.

Hung

P-Y

Lau

K-T

Cheng

L-K

, et al. Impact response of hybrid carbon/glass fibre reinforced polymer composites designed for engineering applications. Compos B Eng 2018; 133: 86–90. https://doi.org/10.1016/j.compositesb.2017.09.026

34.

Hou

Tie

, et al. Low-velocity impact behaviors of repaired CFRP laminates: effect of impact location and external patch configurations. Compos B Eng 2019; 163: 669–680. https://doi.org/10.1016/j.compositesb.2018.12.153

35.

Magliaro

Papini

. Streamlined numerical modeling of solid particle impacts and erosion using a fully eulerian technique: an experimentally validated comparison to Lagrangian and meshfree methods. Wear 2025; 570: 205918. https://doi.org/10.1016/j.wear.2025.205918

36.

Dong

Cao

, et al. A new mechanistic model for abrasive erosion using discrete element method. Powder Technol 2021; 380: 486–496. https://doi.org/10.1016/j.powtec.2020.11.017

37.

Zhang

. DEM-FEM based numerical analysis on mechanical responses of sandy soil and pipeline to seepage erosion. Eng Geol 2022; 310: 106868. https://doi.org/10.1016/j.enggeo.2022.106868

38.

Zhang

Jia

Y-F

Sun

X-W

, et al. A model of erosion rate prediction for component with complex geometry based on numerical simulation. Wear 2024; 546-547: 205328. https://doi.org/10.1016/j.wear.2024.205328

39.

Izard

Ben Hamouda

Vande

. High-stress impact–abrasion test by discrete element modeling. Computational Particle Mechanics 2021; 8(5): 1061–1073. https://doi.org/10.1007/s40571-020-00377-8

40.

Zhou

Liu

Wang

. Finite element analysis on impact response and damage mechanism of composite laminates under single and repeated low-velocity impact. Aero Sci Technol 2022; 129: 107810. https://doi.org/10.1016/j.ast.2022.107810

41.

Zhang

Qiao

Zhang

, et al. Prediction of deformation and failure behavior of continuous fiber reinforced composite fabricated by additive manufacturing. Compos Struct 2021; 265: 113738. https://doi.org/10.1016/j.compstruct.2021.113738

42.

Wang

, et al. Tensile properties and damage evolution of composite joints with lap and stitched structure: experimental and numerical investigation. Thin-Walled Struct 2025; 215: 113531. https://doi.org/10.1016/j.tws.2025.113531

43.

Hashin

. Failure criteria for unidirectional fiber composites. J Appl Mech 1980; 47(2): 329–334. https://doi.org/10.1115/1.3153664

44.

Han

Lei

Liu

, et al. Fatigue damage analysis and life prediction of 3D braided SiC/SiC composites with multi-scale approach. Compos Struct 2024; 349-350: 118447. https://doi.org/10.1016/j.compstruct.2024.118447

45.

Zhang

Yang

Tong

, et al. Multi-scale modeling and tensile failure prediction of 3D needled C/C–SiC composites considering real microstructure. J Mater Res Technol 2023; 25: 120–137. https://doi.org/10.1016/j.jmrt.2023.05.171

46.

Divse

Marla

Joshi

. 3D progressive damage modeling of fiber reinforced plastics laminates including drilling-induced damage. Compos Appl Sci Manuf 2022; 163: 107230. https://doi.org/10.1016/j.compositesa.2022.107230

47.

Zhang

. CFRP FEA data. figshare. Dataset. https://doi.org/10.6084/m9.figshare.31867579.v1.