A composite fan blade composite crack propagation simulation method, program, device and storage medium

Abstract

The present application belongs to the technical field of composite material structure mechanics and computational mechanics, and particularly relates to a composite material fan blade composite crack propagation simulation method, program, equipment and storage medium. The present application introduces key rotation and failure criterion based on key deflection angle into the conventional state type near-field dynamic model, which can more accurately capture and simulate the crack propagation behavior of composite materials under shear or mixed mode load. The plate mapping method proposed by the present application successfully extends the improved PD model suitable for plates to the fan blade airfoil section, solving the problem of complex geometric modeling and being suitable for complex geometric structures. The present application adopts the near-field dynamic method, which does not need to preset the crack path or perform grid reconstruction, and can naturally simulate the complex expansion behaviors such as spontaneous bifurcation and merging of cracks; the PD model can finely describe the mechanical behavior and crack evolution of each layer and the interlayer of the laminate, including layer cracking and interlayer delamination.

Description

Technical Field

[0001] This invention belongs to the field of composite material structural mechanics and computational mechanics, specifically relating to a method, program, equipment, and storage medium for simulating the propagation of composite cracks in composite wind turbine blades. Background Technology

[0002] Wind energy is a clean and renewable energy source, and wind turbines are key equipment for converting wind energy. Wind turbine blades, as their core components, directly affect the stable operation and economic benefits of the entire wind power generation system through their structural safety and reliability. Modern large wind turbine blades are generally made of fiber-reinforced composite materials, which have advantages such as high specific strength, high specific modulus, and excellent fatigue resistance. However, they also have characteristics such as anisotropy, complex structure, and diverse fracture modes.

[0003] During operation, wind turbine blades are subjected to complex loads such as wind loads, gravity, and inertial forces. Damage can easily accumulate within the blade structure, potentially leading to crack initiation. Crack propagation can result in decreased blade stiffness, deterioration of aerodynamic performance, and even catastrophic fracture. Blade crack propagation is often not a simple opening type (Type I), but rather a complex propagation mode involving in-plane shear (Type II). Accurately predicting the propagation path and rate of complex cracks in composite blades is crucial for blade design optimization, life assessment, and maintenance strategy development.

[0004] The traditional finite element method (FEM) is well-established for simulating the deformation of continuous media, but it faces challenges in handling discontinuous problems such as cracks. This is especially true for composite laminate structures, where the anisotropy and the interaction of various fracture modes, such as matrix cracking and delamination, make numerical simulation of crack propagation behavior even more challenging.

[0005] Peridynamics (PD) is an emerging method based on the concept of nonlocal interactions to model and describe the mechanical behavior of matter by solving spatial integral equations. It is naturally suitable for simulating discontinuous problems, such as crack propagation. It does not require pre-setting crack paths and can naturally capture the complex behavior of cracks. Existing PD theories have been extensively studied in simulating type I cracks in composite materials. However, for fractures caused by shear deformation, the classical bond-type PD model is limited by Poisson's ratio. While the conventional PD model overcomes the Poisson's ratio limitation, it still has shortcomings in describing shear deformation and corresponding failure modes. Furthermore, its application to large, complex geometries (such as wind turbine blades) of fiber-reinforced composite laminates, and its consideration of intralaminar / interlaminar complex crack propagation, are still insufficient. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of the prior art and provide a method, program, equipment and storage medium for simulating the propagation of complex cracks in composite wind turbine blades. This invention is a numerical method capable of simulating the propagation process of complex cracks in composite wind turbine blades.

[0007] A method for simulating composite crack propagation in composite wind turbine blades includes the following steps:

[0008] Structural analysis was performed on the composite material wind turbine blades to determine the dangerous section under preset load conditions;

[0009] Construct a PD model of a laminated flat plate with a dangerous section, spatially discretize the PD model of the laminated flat plate, and obtain the coordinates and volume of each particle after discretization;

[0010] A three-dimensional discrete model of the dangerous section is constructed, including the upper and lower airfoils; the plate mapping method is used to map the particles in the laminated plate PD model to the upper and lower airfoil surfaces layer by layer;

[0011] An initial crack is set in the three-dimensional discrete model of the critical section of PD, and boundary conditions and loads are applied. In each iteration, for a particle in a certain layer, before deformation, the particle has intralayer matrix bonds with all particles in its in-plane neighborhood, interlayer normal bonds with particles at the projection positions of its adjacent layers, and interlayer shear bonds with all particles in its interlayer neighborhood. Considering the failure criterion of intralayer matrix bond fracture, the force density vector of each particle in the in-plane neighborhood of the particle after deformation is calculated. Considering the failure criterion of interlayer normal bonds, the interlayer normal force density vector of each particle at the projection positions of its adjacent layers after deformation is calculated. Considering the failure criterion of interlayer shear bonds, the interlayer shear force density vector of each particle in the interlayer neighborhood of the particle after deformation is calculated. The crack propagation is simulated by iteratively solving the PD motion equations.

[0012] Furthermore, the composite material of the wind turbine blade includes fiber material and matrix material; the upper and lower airfoils of the dangerous section of the composite wind turbine blade are simulated by N layers of mass points to construct a laminated flat plate PD model with 2N layers.

[0013] Furthermore, the construction of the three-dimensional PD discrete model of the dangerous truncation is specifically as follows:

[0014] For the upper airfoil, the N+1th layer of the laminated flat plate PD model is segmented and mapped to the 1st layer of the upper airfoil. Then, the N+2th to 2Nth layers of the laminated flat plate PD model are mapped to the 2nd to Nth layers of the upper airfoil using the normal extension method.

[0015] For the lower airfoil, the Nth layer of the laminated flat plate PD model is segmented and mapped to the Nth layer of the lower airfoil. Then, the 1st to N-1th layers of the laminated flat plate PD model are mapped to the 1st to N-1th layers of the lower airfoil using the normal extension method.

[0016] Furthermore, the failure criterion considering intralayer matrix bond fracture calculates the force density vector of each particle in its in-plane neighborhood after deformation, specifically as follows:

[0017] For a particle j in the m-th layer, before deformation, particle j has an intra-layer matrix bond with all particles in its in-plane neighborhood. The in-plane neighborhood of particle j is a circular region with particle j as the center and δ as the radius.

[0018] For particle k in the in-plane neighborhood of particle j, before deformation, the position of particle j is: The position of particle k is After deformation, the position of particle j becomes The position of particle k becomes Then the bond elongation of the matrix bond within the layer between particle j and particle k Bond deflection angle for:

[0019]

[0020] in, and Let J be the displacement of particle j and particle k after deformation. n f t is a unit vector parallel to the fiber direction. f It is a unit vector perpendicular to the fiber direction;

[0021] The failure criterion for the breakage of the matrix bonds within the layer is as follows:

[0022] If the intralayer matrix bond between particle j and particle k is along the fiber direction, it will not break;

[0023] If the intralayer matrix bond between particle j and particle k is not along the fiber direction, then determine the bond elongation of the intralayer matrix bond. or key deflection angle Whether the critical value is exceeded will cause failure. The scalar function μ1(j,k) characterizing the matrix bond state within the layer is:

[0024]

[0025] in, Let be the orientation parameter of the intralayer matrix bond between particle j and its in-plane neighboring particle k. If the intralayer matrix bond between particle j and particle k is along the fiber direction, then... If the intralayer matrix bond between particle j and particle k is not along the fiber direction, then s c With γ c These are the critical values ​​for bond elongation and bond deflection angle of the matrix bonds within the layer, respectively.

[0026]

[0027] Among them, K m μ m μ c These are the bulk modulus, Lamé constant, and shear modulus of the matrix material, respectively; G IC The first critical energy release rate;

[0028] Considering the failure criterion of intralayer matrix bond fracture, the force density vector of particle k in the in-plane neighborhood of particle j in the m-th layer after deformation is... for:

[0029]

[0030] in, Q 11 =E 11 / (1-ν 21 ν 12 ), Q 12 =ν 12 E 22 / (1-ν 21 ν 12 ), Q 22 =E 22 / (1-ν 21 ν 12 ),

[0031] Among them, E 11 E represents the elastic modulus of the fibrous material along the fiber direction. 22 G represents the elastic modulus of the fiber material perpendicular to the fiber direction. 12 ν is the in-plane shear modulus of the fiber material. 12 h is the in-plane Poisson's ratio of the fiber material. c The thickness of a single-layer composite material plate; Let k be the set of particles in the in-plane neighborhood of particle k in the m-th layer before deformation. This indicates that particle α lies in the in-plane neighborhood of particle k. Let be the volume of particle α in the m-th layer; Let be the orientation parameter of the matrix bond between particle j and its in-plane neighboring particle k. If the matrix bond between particle j and particle k is perpendicular to the fiber direction, then... If the intralayer matrix bond between particle j and particle k is not perpendicular to the fiber direction, then

[0032] Furthermore, the failure criterion considering interlayer normal bonds calculates the interlayer normal force density vector of the mass point at the projected position of the adjacent layer after deformation, specifically as follows:

[0033] For particle j in the m-th layer, the particle at the projection position of particle j in the (m+1)-th layer is j. + The particle j at the projection position in the (m-1)th layer is called particle j. - ; particle j in layer m and particle j in layer (m+1) + The particle j in the m-th layer and the particle j in the (m-1)-th layer are... - Interlayer normal keys exist;

[0034] Before deformation, the position of particle j is point j + The position is point j - The position is After deformation, the position of particle j becomes point j + The position becomes point j - The position becomes Then particle j and particle j + Bond elongation of interlayer normal bonds particle j and particle j - Bond elongation of interlayer normal bonds for:

[0035]

[0036] The failure criterion for the interlayer normal key is:

[0037] Interlayer normal bonds fail when the bond elongation exceeds a critical value. The scalar function μ2(j,j) characterizes the state of the interlayer normal bonds. + ) and μ2(j,j - )for:

[0038]

[0039] Among them, s d This is the critical value for the bond elongation of the interlayer normal bonds;

[0040]

[0041] Among them, E c The elastic modulus of the matrix material;

[0042] Considering the failure criterion of interlayer normal bonds, the mass j in the (m+1)th layer after deformation + For the interlayer normal force density of particle j in the m-th layer The particle in the (m-1)th layer after deformation is j - For the interlayer normal force density of particle j in the m-th layer for:

[0043]

[0044] in, h m Let be the thickness of the m-th layer.

[0045] Furthermore, the failure criterion considering interlaminar shear bonds calculates the interlaminar shear force density vector of each particle in the interlaminar neighborhood of the deformed particle, specifically as follows:

[0046] For a particle j in the m-th layer, before deformation, particle j and all particles in its interlayer neighborhood have interlayer shear bonds; the interlayer neighborhood of particle j is the projection region of the in-plane neighborhood of particle j in the (m+1)-m-1 layers.

[0047] For particle i in the interlayer neighborhood of particle j, particle i is in the nth layer, n = m + 1 or n = m - 1; before deformation, the position of particle j is... The position of particle i is After deformation, the position of particle j becomes The position of particle i becomes Then the shear angle of the interlaminar shear bond between particle j and particle i for:

[0048]

[0049] The failure criterion for the breakage of the interlaminar shear bonds is as follows:

[0050] Interlayer shear bonds at shear angle It will fail when the critical value is exceeded. The scalar function μ3(j,i) characterizing the interlayer shear bond state is:

[0051]

[0052] Where, φ c The critical shear angle for interlaminar shear bonds;

[0053]

[0054] Among them, G IIC The type II critical energy release rate; G c The shear modulus of the matrix material;

[0055] Considering the failure criterion of interlaminar shear bonds, the interlaminar shear force density of particle i in the nth layer with respect to particle j in the mth layer after deformation. for:

[0056]

[0057] in, Let be the interlayer shear angle between particle i in the nth layer and particle j in the mth layer.

[0058] Furthermore, the motion equation of the PD is:

[0059]

[0060] Where ρ is the density of the composite material wind turbine blade; Let be the set of particles in the in-plane neighborhood of particle j in the m-th layer before deformation; Let be the set of particles in the interlayer neighborhood of particle j in the (m-1)th layer before deformation; Let be the set of particles in the interlayer neighborhood of particle j in the (m+1)th layer before deformation; Let be the mass density of particle j in the m-th layer after deformation.

[0061] A computer device / equipment / system includes a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above-described method for simulating the composite crack propagation of composite wind turbine blades.

[0062] A computer-readable storage medium having a computer program / instructions stored thereon, which, when executed by a processor, implements the steps of the above-described method for simulating the composite crack propagation of composite wind turbine blades.

[0063] A computer program product includes a computer program / instructions that, when executed by a processor, implement the steps of the above-described method for simulating the composite crack propagation of composite wind turbine blades.

[0064] The beneficial effects of this invention are as follows:

[0065] This invention introduces bond rotation and a failure criterion based on bond deflection angle into a conventional near-field dynamics model, enabling more accurate capture and simulation of crack propagation behavior in composite materials under shear or mixed-mode loading. The proposed flat plate mapping method successfully extends the improved PD model, suitable for flat plates, to wind turbine blade airfoil sections, solving the challenge of complex geometric modeling and making it applicable to complex geometric structures. Employing a near-field dynamics method, this invention eliminates the need for pre-setting crack paths or mesh reconstruction, naturally simulating complex propagation behaviors such as spontaneous crack bifurcation and merging. The PD model can meticulously describe the mechanical behavior and crack evolution of each layer and interlayer within the laminate, including intralayer cracking and interlayer delamination. The simulation results of this invention can provide detailed, layer-by-layer crack propagation information and states within the blade section, providing crucial numerical data for refined design, strength verification, life prediction, and maintenance strategy formulation of wind turbine blades. Attached Figure Description

[0066] Figure 1 This is the overall architecture diagram of the present invention.

[0067] Figure 2 This is a schematic diagram of the interactions between particles in a PD model that considers bond rotation effects, illustrating bond elongation and bond rotation.

[0068] Figure 3 This is a schematic diagram of the intralayer bond deflection angle.

[0069] Figure 4 In one embodiment of the present invention, the stress-displacement distribution cloud map of the wind turbine blade under uniform load is determined by the finite element method, wherein (a) is the overall stress-displacement distribution and (b) is the maximum stress section.

[0070] Figure 5 In one embodiment of the present invention, the stress-displacement distribution cloud map of the wind turbine blade under linear load is determined by the finite element method, wherein (a) is the overall stress-displacement distribution and (b) is the maximum stress section.

[0071] Figure 6 This is an embodiment of the present invention, showing the stress distribution cloud diagram of the blade under loads in different directions.

[0072] Figure 7 This is a schematic diagram of the blade segment PD model and initial crack setting used for crack propagation simulation in an embodiment of the present invention.

[0073] Figure 8 This is a diagram showing the damage of each layer of the blade segment in an embodiment of the present invention when only the bond deflection angle is considered as the failure criterion.

[0074] Figure 9This is an embodiment of the present invention, showing the damage of each layer of the blade segment when both bond elongation and deflection angle are considered as failure criteria.

[0075] Figure 10 This is an embodiment of the present invention, showing the damage of each ply of the blade section when a new ply configuration is adopted (the ply configuration of the lower wing is set to [0° / -45° / 90°]2, and the ply configuration of the upper wing is set to [0° / 45° / 90°]2).

[0076] Figure 11 This is a schematic diagram of a blade segment with intersecting initial cracks in an embodiment of the present invention.

[0077] Figure 12 This is a diagram illustrating the propagation of cross cracks in each layer of the blade section in an embodiment of the present invention. Detailed Implementation

[0078] The present invention will now be further described with reference to the accompanying drawings.

[0079] This invention discloses a method for simulating the propagation of complex cracks in composite wind turbine blades, relating to wind turbine blade fracture prediction technology. This invention aims to address the problem that existing methods struggle to accurately simulate the initiation and propagation of shear cracks in composite wind turbine blades under complex loads.

[0080] First, this invention proposes a conventional mode-of-sense near-field dynamics (OSB-PD) method that simultaneously considers bond tension and shear deformation to describe the combined failure mode of intralaminar shear and interlaminar delamination in fiber-reinforced composites. Second, through finite element analysis, a critical section of the wind turbine blade under specific load conditions is identified. For this critical section, the proposed near-field dynamics method is applied to establish a discrete mass model of the actual airfoil surface structure. Finally, an initial crack is set in the PD model, simulated loads are applied, and a failure criterion based on bond fracture is used to simulate the combined crack propagation path and failure evolution process under different layup configurations and different numbers of cracks. This invention can effectively simulate the combined crack propagation of laminated blade structures, providing an accurate numerical simulation method for wind turbine blade structural design optimization and safety assessment.

[0081] A method for simulating composite crack propagation in composite wind turbine blades includes the following steps:

[0082] Step 1: Perform structural analysis on the composite material wind turbine blades to determine the critical section under preset load conditions;

[0083] The composite material includes fiber materials and matrix materials;

[0084] Step 2: Construct a laminated flat plate PD model of the dangerous section, spatially discretize the laminated flat plate PD model, and obtain the coordinates and volume of each particle after discretization;

[0085] The upper and lower airfoils of the dangerous section of the composite wind turbine blade are simulated using N layers of mass points to construct a laminated flat plate PD model with 2N layers.

[0086] Step 3: Construct a three-dimensional discrete PD model of the dangerous section, including the upper and lower airfoils; use the plate mapping method to map the particles in the laminated plate PD model layer by layer to the surfaces of the upper and lower airfoils;

[0087] For the upper airfoil, the N+1th layer of the laminated flat plate PD model is segmented and mapped to the 1st layer of the upper airfoil. Then, the N+2th to 2Nth layers of the laminated flat plate PD model are mapped to the 2nd to Nth layers of the upper airfoil using the normal extension method.

[0088] For the lower airfoil, the Nth layer of the laminated flat plate PD model is segmented and mapped to the Nth layer of the lower airfoil. Then, the 1st to N-1th layers of the laminated flat plate PD model are mapped to the 1st to N-1th layers of the lower airfoil using the normal extension method.

[0089] Step 4: Set an initial crack in the three-dimensional PD discrete model of the critical section, apply boundary conditions and loads, and iteratively solve the PD motion equations according to the failure criterion based on key fracture considering rotation angle to simulate crack propagation.

[0090] Step 4.1: For particle j in the m-th layer, before deformation, particle j has an intra-layer matrix bond with all particles in its in-plane neighborhood. The in-plane neighborhood of particle j is a circular region with particle j as the center and δ as the radius.

[0091] For particle k in the in-plane neighborhood of particle j, before deformation, the position of particle j is: The position of particle k is After deformation, the position of particle j becomes The position of particle k becomes Then the bond elongation of the matrix bond within the layer between particle j and particle k Bond deflection angle for:

[0092]

[0093] in, and Let J be the displacement of particle j and particle k after deformation. n f t is a unit vector parallel to the fiber direction. f It is a unit vector perpendicular to the fiber direction;

[0094] The failure criterion for intralayer matrix bond breakage is:

[0095] If the intralayer matrix bond between particle j and particle k is along the fiber direction, it will not break;

[0096] If the intralayer matrix bond between particle j and particle k is not along the fiber direction, then determine the bond elongation of the intralayer matrix bond. or key deflection angle Whether the critical value is exceeded will cause failure. The scalar function μ1(j,k) characterizing the matrix bond state within the layer is:

[0097]

[0098] in, Let be the orientation parameter of the intralayer matrix bond between particle j and its in-plane neighboring particle k. If the intralayer matrix bond between particle j and particle k is along the fiber direction, then... If the intralayer matrix bond between particle j and particle k is not along the fiber direction, then s c With γ c These are the critical values ​​for bond elongation and bond deflection angle of the matrix bonds within the layer, respectively.

[0099]

[0100] Among them, K m μ m μ c These are the bulk modulus, Lamé constant, and shear modulus of the matrix material, respectively; G IC The first critical energy release rate;

[0101] Considering the failure criterion of intralayer matrix bond fracture, the force density vector of particle k in the in-plane neighborhood of particle j in the m-th layer after deformation is... for:

[0102]

[0103] in, Q 11 =E 11 / (1-ν 21 ν 12 ), Q 12 =ν 12 E 22 / (1-ν 21 ν 12 ), Q 22 =E 22 / (1-ν 21 ν12 ),

[0104] Among them, E 11 E represents the elastic modulus of the fibrous material along the fiber direction. 22 G represents the elastic modulus of the fiber material perpendicular to the fiber direction. 12 ν is the in-plane shear modulus of the fiber material. 12 h is the in-plane Poisson's ratio of the fiber material. c The thickness of a single-layer composite material plate; Let k be the set of particles in the in-plane neighborhood of particle k in the m-th layer before deformation. This indicates that particle α lies in the in-plane neighborhood of particle k. Let be the volume of particle α in the m-th layer; Let be the orientation parameter of the matrix bond between particle j and its in-plane neighboring particle k. If the matrix bond between particle j and particle k is perpendicular to the fiber direction, then... If the intralayer matrix bond between particle j and particle k is not perpendicular to the fiber direction, then

[0105] Step 4.2: For particle j in layer m, the particle at the projected position of particle j in layer m+1 is j + The particle j at the projection position in the (m-1)th layer is called particle j. - ; particle j in layer m and particle j in layer (m+1) + The particle j in the m-th layer and the particle j in the (m-1)-th layer are... - Interlayer normal keys exist;

[0106] Before deformation, the position of particle j is point j + The position is point j - The position is After deformation, the position of particle j becomes point j + The position becomes point j - The position becomes Then particle j and particle j + Bond elongation of interlayer normal bonds particle j and particle j - Bond elongation of interlayer normal bonds for:

[0107]

[0108] The failure criterion for interlayer normal keys is:

[0109] Interlayer normal bonds fail when the bond elongation exceeds a critical value. The scalar function μ2(j,j) characterizes the state of the interlayer normal bonds. + ) and μ2(j,j - )for:

[0110]

[0111] Among them, s d This is the critical value for the bond elongation of the interlayer normal bonds;

[0112]

[0113] Among them, E c The elastic modulus of the matrix material;

[0114] Considering the failure criterion of interlayer normal bonds, the mass j in the (m+1)th layer after deformation + For the interlayer normal force density of particle j in the m-th layer The particle in the (m-1)th layer after deformation is j - For the interlayer normal force density of particle j in the m-th layer for:

[0115]

[0116] in, h m Let m be the thickness of the m-th layer;

[0117] Step 4.3: For particle j in the m-th layer, before deformation, particle j and all particles in its interlayer neighborhood have interlayer shear bonds; the interlayer neighborhood of particle j is the projection region of the in-plane neighborhood of particle j in the (m+1)-m-1 layers.

[0118] For particle i in the interlayer neighborhood of particle j, particle i is in the nth layer, n = m + 1 or n = m - 1; before deformation, the position of particle j is... The position of particle i is After deformation, the position of particle j becomes The position of particle i becomes Then the shear angle of the interlaminar shear bond between particle j and particle i for:

[0119]

[0120] Failure criteria for interlaminar shear bond breakage:

[0121] Interlayer shear bonds at shear angle It will fail when the critical value is exceeded. The scalar function μ3(j,i) characterizing the interlayer shear bond state is:

[0122]

[0123] Where, φ c The critical shear angle for interlaminar shear bonds;

[0124]

[0125] Among them, G IIC The type II critical energy release rate; G c The shear modulus of the matrix material;

[0126] Considering the failure criterion of interlaminar shear bonds, the interlaminar shear force density of particle i in the nth layer with respect to particle j in the mth layer after deformation. for:

[0127]

[0128] in, Let be the interlayer shear angle between particle i in the nth layer and particle j in the mth layer.

[0129] Step 4.4: The motion equation of the PD is:

[0130]

[0131] Where ρ is the density of the composite material wind turbine blade; Let be the set of particles in the in-plane neighborhood of particle j in the m-th layer before deformation; Let be the set of particles in the interlayer neighborhood of particle j in the (m-1)th layer before deformation; Let be the set of particles in the interlayer neighborhood of particle j in the (m+1)th layer before deformation; Let be the mass density of particle j in the m-th layer after deformation.

[0132] Example 1:

[0133] like Figure 1 As shown, this embodiment mainly includes the following steps:

[0134] Step 1: Obtain the geometric model, material properties, and layup design information of the IEA15MW wind turbine blade. Use ABAQUS software to build a three-dimensional finite element model of the blade, setting the blade root to be completely fixed, and apply loads to simulate wind loads. Consider two typical loads:

[0135] (a) Pressure evenly distributed along the blade surface;

[0136] (b) Pressure that increases linearly along the blade's span;

[0137] A dynamic implicit analysis method was used to simulate the blade's response within 0-10 seconds, with a time step of 0.025 seconds. The stress and displacement distributions of the blade were calculated.

[0138] Step 2, analyze the stress contour plot (e.g.) Figure 4 and Figure 5 As shown in the figure, under uniformly distributed load, the maximum stress occurs at a distance of 0.2R from the blade root; under linear load, the maximum stress occurs at a distance of 0.3R from the blade root. These regions are potential danger zones for the blade. Simultaneously, the influence of different load directions on the stress extreme values ​​was analyzed (e.g., Figure 6 As shown in the figure, the stress at the wing spars cap was generally high. Taking all factors into consideration, a segment located 0.25R from the blade root was selected as the object for subsequent PD simulations. The airfoil at this location is FFA-W3-36 with a chord length of 5.7m. A segment length of 1m was selected.

[0139] Step 3, the modeling process for the PD geometric model of the blade segment is as follows:

[0140] (1) Establish a flat plate model. Based on the blade chord length and the blade section length, establish a model with a length L = 5.7m, a width W = 1m, and a thickness h = 8.83 × 10⁻⁶ m. -3 The laminated plate model has a particle spacing Δx = 0.01m and a total of 12 layers. Considering the boundary layer, three more particle layers need to be added in the extension direction. Therefore, the entire laminated plate model is replaced by 570 × (100 + 3) × 12 particles.

[0141] (2) Establish the blade segment model. The upper and lower airfoils are divided into two curves for fitting. The two curves are represented by sixth-order function expressions. First, the particles of the 7th layer of the laminated plate are segmented and mapped to the 1st layer of the upper airfoil. Then, the particles of layers 8-12 of the laminated plate are mapped to layers 2-6 of the upper airfoil using the normal extension method. For the modeling of the lower airfoil, the 6th layer of the laminated plate is selected and segmented and mapped to the 6th layer of the lower airfoil. Then, the particles of layers 1-5 of the laminated plate are mapped to layers 1-5 of the lower airfoil using the same normal extension method. According to this mapping method, the layer numbering of the laminated plate can be kept consistent with the layer numbering of the airfoil.

[0142] (3) Set the blade material properties. To simplify modeling, the material of the blade section is set to laminate, and the material is kept consistent throughout.

[0143] Step 4: A near-field dynamics model considering rotation angle is used to describe the behavior of the composite material of the blade section, such as... Figure 2 As shown.

[0144] Step 5: Establish a bond failure criterion between particles that reflects the composite fracture behavior of composite materials. The core of this criterion is the introduction of a judgment based on the "bond deflection angle." The bond deflection angle is defined as the ratio of the projection of the relative displacement vector after deformation onto the fiber direction to the projection of the relative position of the material point before deformation onto the direction perpendicular to the fiber. For example... Figure 3 As shown.

[0145] Step 6: Set two different numbers of initial cracks and corresponding boundary conditions and loads in the PD model.

[0146] a. The layup is [90° / 0° / 90°]4. The boundary conditions are set to restrict displacement only in the y-direction at one end of the blade section. The crack length is a = 0.2m, the crack is parallel to the x-axis, and the midpoint is located at a chord length of 1.5m. Figure 7 As shown, this is used to compare the simulated crack propagation under two failure criteria: considering only the bond rotation angle and considering both bond elongation and deflection angle. Then, the layup configuration was changed: the lower flange layup was set to [0° / -45° / 90°]², and the upper flange layup was set to [0° / 45° / 90°]², without changing the boundary conditions and initial crack settings, to compare the crack propagation under different layup configurations.

[0147] b. The layup setting remains [90° / 0° / 90°]4. The boundary conditions are set to restrict displacement only in the y-direction at one end of the blade section. The initial cracks are set as two intersecting cracks, each with a length a = 0.2 m, and the intersection center is located at the center point of the airfoil. Figure 11 As shown. Used to compare the effects of different crack modes on crack propagation.

[0148] Step 7: Perform a quasi-static solution using the Adaptive Dynamic Relaxation (ADR) method. Set up loading steps, gradually increasing the external load or applying displacement. In each iteration step, update the mass displacements and calculate the critical tensile elongation s of all matrix bonds. c Critical bond deflection angle γ c Critical elongation s of interlayer normal bonds in Critical shear angle φ of interlaminar shear bonds c The state of the key is determined and updated based on the failure criterion in step 5, and the process continues iterating until the target simulation time is reached.

[0149] Step 8: Output the damage function for each time step. Distribution. Draw damage contour maps for each layer of the blade section, such as... Figure 8As shown, it can be observed that under a 4-ply configuration of [90° / 0° / 90°], cracks in the 90° ply primarily propagate along the fiber direction, exhibiting shear damage characteristics. Under the deflection angle criterion alone, the 0° ply may not propagate. Under the combined criterion, tensile damage occurs along the fiber direction. If the ply configuration is changed, with the lower flange ply set to [0° / -45° / 90°]2 and the upper flange ply set to [0° / 45° / 90°]2, cracks in the ±45° layers propagate along the 45° direction. If intersecting cracks are used, the cracks propagate from the initial crack tip along the favorable direction of their respective plies. The crack propagation paths and damage extent under different failure criteria, different plies, and different numbers of cracks are analyzed to evaluate the crack resistance of the blade section.

[0150] This embodiment details the entire process of applying the method of this invention: The critical section of the IEA 15MW blade was determined using the finite element method; a PD model of this complex section was established using plate mapping technology; and the composite crack propagation behavior under different failure criteria, different layup designs, and different initial crack configurations was successfully simulated on this model. The simulation results revealed the influence of the failure criteria on damage mode judgment, verified the effectiveness of ±45° layup in suppressing shear damage, and confirmed that the crack propagation path is mainly dominated by the fiber direction. These results demonstrate the effectiveness and practical value of the method of this invention, providing a powerful tool for the refined design and safety assessment of composite wind turbine blades.

[0151] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for simulating the propagation of composite cracks in composite wind turbine blades, characterized in that: Structural analysis was performed on the composite material wind turbine blades to determine the dangerous section under preset load conditions; Construct a PD model of a laminated flat plate with a dangerous section, spatially discretize the PD model of the laminated flat plate, and obtain the coordinates and volume of each particle after discretization; A three-dimensional discrete model of the dangerous section is constructed, including the upper and lower airfoils; the plate mapping method is used to map the particles in the laminated plate PD model to the upper and lower airfoil surfaces layer by layer; An initial crack is set in the three-dimensional PD discrete model of the dangerous section, and boundary conditions and loads are applied. In each iteration, for a particle within a certain layer, before deformation, the particle has intralayer matrix bonds with all particles in its in-plane neighborhood, interlayer normal bonds with particles at the projection position of its adjacent layer, and interlayer shear bonds with all particles in its interlayer neighborhood. Considering the failure criterion of intralayer matrix bond breakage, the force density vector of each particle in the in-plane neighborhood of the particle after deformation is calculated. Considering the failure criterion of interlaminar normal bonds, the interlaminar normal force density vector of the particle at the projection position of the adjacent layer of the deformed particle is calculated; considering the failure criterion of interlaminar shear bonds, the interlaminar shear force density vector of each particle in the interlaminar neighborhood of the deformed particle is calculated; crack propagation is simulated by iteratively solving the PD motion equation.

2. The method for simulating composite crack propagation in composite wind turbine blades according to claim 1, characterized in that: The composite material of the wind turbine blade includes fiber material and matrix material; the upper and lower airfoils of the dangerous section of the composite wind turbine blade are simulated by N layers of mass points to construct a laminated flat plate PD model with 2N layers.

3. The method for simulating composite crack propagation in composite wind turbine blades according to claim 2, characterized in that: The construction of the three-dimensional PD discrete model of the dangerous truncation is specifically as follows: For the upper airfoil, the N+1th layer of the laminated flat plate PD model is segmented and mapped to the 1st layer of the upper airfoil. Then, the N+2th to 2Nth layers of the laminated flat plate PD model are mapped to the 2nd to Nth layers of the upper airfoil using the normal extension method. For the lower airfoil, the Nth layer of the laminated flat plate PD model is segmented and mapped to the Nth layer of the lower airfoil. Then, the 1st to N-1th layers of the laminated flat plate PD model are mapped to the 1st to N-1th layers of the lower airfoil using the normal extension method.

4. The method for simulating composite crack propagation in composite wind turbine blades according to claim 2, characterized in that: The failure criterion considering intralayer matrix bond breakage calculates the force density vector of each particle in the in-plane neighborhood of the deformed particle, specifically: For a particle j in the m-th layer, before deformation, particle j has an intra-layer matrix bond with all particles in its in-plane neighborhood. The in-plane neighborhood of particle j is a circular region with particle j as the center and δ as the radius. For particle k in the in-plane neighborhood of particle j, before deformation, the position of particle j is: The position of particle k is After deformation, the position of particle j becomes The position of particle k becomes Then the bond elongation of the matrix bond within the layer between particle j and particle k Bond deflection angle for: in, and Let J be the displacement of particle j and particle k after deformation. n f t is a unit vector parallel to the fiber direction. f It is a unit vector perpendicular to the fiber direction; The failure criterion for the breakage of the matrix bonds within the layer is as follows: If the intralayer matrix bond between particle j and particle k is along the fiber direction, it will not break; If the intralayer matrix bond between particle j and particle k is not along the fiber direction, then determine the bond elongation of the intralayer matrix bond. or key deflection angle Whether the critical value is exceeded will cause failure. The scalar function μ1(j,k) characterizing the matrix bond state within the layer is: in, Let be the orientation parameter of the intralayer matrix bond between particle j and its in-plane neighboring particle k. If the intralayer matrix bond between particle j and particle k is along the fiber direction, then... If the intralayer matrix bond between particle j and particle k is not along the fiber direction, then s c With γ c These are the critical values ​​for bond elongation and bond deflection angle of the matrix bonds within the layer, respectively. Among them, K m μ m μ c These are the bulk modulus, Lamé constant, and shear modulus of the matrix material, respectively; G IC The first critical energy release rate; Considering the failure criterion of intralayer matrix bond fracture, the force density vector of particle k in the in-plane neighborhood of particle j in the m-th layer after deformation is... for: Among them, Q 11 =E 11 / (1-n 21 n 12 ),Q 12 =n 12 E 22 / (1-n 21 n 12 ),Q 22 =E 22 / (1-n 21 n 12 ),Q 66 =G 12 , Among them, E 11 E represents the elastic modulus of the fibrous material along the fiber direction. 22 G represents the elastic modulus of the fiber material perpendicular to the fiber direction. 12 ν is the in-plane shear modulus of the fiber material. 12 h is the in-plane Poisson's ratio of the fiber material. c The thickness of a single-layer composite material plate; Let k be the set of particles in the in-plane neighborhood of particle k in the m-th layer before deformation. This indicates that particle α lies in the in-plane neighborhood of particle k. Let be the volume of particle α in the m-th layer; Let be the orientation parameter of the matrix bond between particle j and its in-plane neighboring particle k. If the matrix bond between particle j and particle k is perpendicular to the fiber direction, then... If the intralayer matrix bond between particle j and particle k is not perpendicular to the fiber direction, then 5. The method for simulating composite crack propagation in composite wind turbine blades according to claim 4, characterized in that: The failure criterion considering interlayer normal keys calculates the interlayer normal force density vector of the mass point at the projected position of the adjacent layer after deformation, specifically as follows: For particle j in the m-th layer, the particle at the projection position of particle j in the (m+1)-th layer is j. + The particle j at the projection position in the (m-1)th layer is called particle j. - ; particle j in layer m and particle j in layer (m+1) + The particle j in the m-th layer and the particle j in the (m-1)-th layer are... - Interlayer normal keys exist; Before deformation, the position of particle j is point j + The position is point j - The position is After deformation, the position of particle j becomes point j + The position becomes point j - The position becomes Then particle j and particle j + Bond elongation of interlayer normal bonds particle j and particle j - Bond elongation of interlayer normal bonds for: The failure criterion for the interlayer normal key is: Interlayer normal bonds fail when the bond elongation exceeds a critical value. The scalar function μ2(j,j) characterizes the state of the interlayer normal bonds. + ) and μ2(j,j - )for: Among them, s d This is the critical value for the bond elongation of the interlayer normal bonds; Among them, E c The elastic modulus of the matrix material; Considering the failure criterion of interlayer normal bonds, the mass j in the (m+1)th layer after deformation + For the interlayer normal force density of particle j in the m-th layer The particle in the (m-1)th layer after deformation is j - For the interlayer normal force density of particle j in the m-th layer for: in, h m Let be the thickness of the m-th layer.

6. The method for simulating composite crack propagation in composite wind turbine blades according to claim 5, characterized in that: The failure criterion considering interlaminar shear bonds calculates the interlaminar shear force density vector of each particle in the interlaminar neighborhood of the deformed particle in the direction of the shear bond. Specifically: For a particle j in the m-th layer, before deformation, particle j and all particles in its interlayer neighborhood have interlayer shear bonds; the interlayer neighborhood of particle j is the projection region of the in-plane neighborhood of particle j in the (m+1)-m-1 layers. For particle i in the interlayer neighborhood of particle j, particle i is in the nth layer, n = m + 1 or n = m - 1; before deformation, the position of particle j is... The position of particle i is After deformation, the position of particle j becomes The position of particle i becomes Then the shear angle of the interlaminar shear bond between particle j and particle i for: The failure criterion for the breakage of the interlaminar shear bonds is as follows: Interlayer shear bonds at shear angle It will fail when the critical value is exceeded. The scalar function μ3(j,i) characterizing the interlayer shear bond state is: Where, φ c The critical shear angle for interlaminar shear bonds; Among them, G IIC The type II critical energy release rate; G c The shear modulus of the matrix material; Considering the failure criterion of interlaminar shear bonds, the interlaminar shear force density of particle i in the nth layer with respect to particle j in the mth layer after deformation. for: in, Let be the interlayer shear angle between particle i in the nth layer and particle j in the mth layer.

7. The method for simulating composite crack propagation in composite wind turbine blades according to claim 6, characterized in that: The motion equation of the PD is: Where ρ is the density of the composite material wind turbine blade; Let be the set of particles in the in-plane neighborhood of particle j in the m-th layer before deformation; Let be the set of particles in the interlayer neighborhood of particle j in the (m-1)th layer before deformation; Let be the set of particles in the interlayer neighborhood of particle j in the (m+1)th layer before deformation; Let be the mass density of particle j in the m-th layer after deformation.

8. A computer device / equipment / system, comprising a memory, a processor, and a computer program stored in the memory, characterized in that: The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program / instructions stored thereon, characterized in that: When the computer program / instructions are executed by the processor, they implement the steps of the method according to any one of claims 1 to 7.

10. A computer program product comprising a computer program / instructions, characterized in that: When the computer program / instructions are executed by the processor, they implement the steps of the method according to any one of claims 1 to 7.

Images