ABAQUS-based composite material lightning stroke dynamic ablation simulation method

Abstract

The invention discloses a composite material lightning stroke dynamic ablation simulation method based on ABAQUS. The method comprises the specific steps that S1, an anisotropic composite material three-dimensional model is established; s2, applying lightning stroke current and heat flow load, setting a variable initial value at normal temperature, and calculating the pyrolysis degree, potential and temperature field of each node; s3, giving a grid deformation displacement load in a static ablation depth direction, and driving the nodes to move to simulate deformation; and S4, constructing an information transfer matrix, mapping a temperature update grid, and taking the update grid as a next electrothermal analysis initial material domain. According to the ABAQUS-based composite material lightning stroke dynamic ablation simulation method provided by the invention, flowing of a dynamic grid between different layers of an anisotropic composite material laminated plate is realized, and parallel coupling calculation of composite material temperature rise and surface ablation depression is realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of lightning ablation damage simulation methods, and in particular relates to an ABAQUS-based lightning ablation simulation method for composite materials. Background Art

[0002] Carbon fiber-reinforced resin-based composites (CFRPs) are increasingly being used in civil and military aircraft due to their exceptional properties of lightweight, high strength, and high modulus. However, when a composite material is struck by lightning, its extremely low electrical conductivity makes it difficult to quickly release the energy generated by the lightning strike, resulting in high temperatures ranging from several hundred to several thousand Kelvin in the area where the lightning strike was applied. When the temperature reaches 700 K, the resin completely pyrolyzes, and as the temperature rises further, the resin vaporizes. At 3589 K, the carbon fiber sublimates. As the phase transition occurs, the composite surface gradually ablates, and the area where the lightning strike was applied continues to expand and deform, forming dynamic concave damage. To simulate lightning ablation damage in CFRPs, ABAQUS provides user-defined lightning arc load models and material pyrolysis kinetic models. Through solid-state electro-thermal-chemical coupled calculations, it solves for the temperature, degree of pyrolysis, and displacement of the composite material along the fiber, matrix, and depth directions, thereby accurately simulating the dynamic ablation damage of CFRPs during a lightning strike.

[0003] When simulating the dynamic ablation deformation of anisotropic composite materials, the UMESHMOTION user subroutine provided by ABAQUS can only simulate the deformation and damage of isotropic materials using a dynamic mesh. It does not allow the mesh to flow between different layers of anisotropic composite laminates during the material concave deformation process. This is because when the anisotropic material degenerates from one layer to another, the severely distorted mesh at the interface between the different layers will terminate the dynamic mesh movement analysis simulated by the UMESHMOTION user subroutine. Therefore, ABAQUS cannot yet effectively simulate the dynamic ablation of anisotropic composite materials. Summary of the Invention

[0004] The purpose of the present invention is to provide an ABAQUS-based simulation method for dynamic ablation of composite materials caused by lightning strikes, which uses a manual grid movement method to realize the flow of dynamic grids between different layers of anisotropic composite laminates.

[0005] The technical solution adopted by the present invention is a method for simulating dynamic ablation of composite materials by lightning strike based on ABAQUS, which specifically includes the following steps: S1, establish a three-dimensional geometric model of anisotropic composite materials, divide the mesh, define material properties related to pyrolysis degree and set boundary conditions; S2, applies lightning current and heat flux loads, sets the initial values ​​of variables to those at room temperature, and calculates the thermal decomposition degree, electric potential, and temperature field of each node based on electrothermal analysis. When the temperature exceeds 700K, the static ablation rate and depth are calculated. If the cumulative depth exceeds the threshold, the electrothermal analysis is terminated. S3 sets the direction of the static ablation depth and converts it into a displacement load on the mesh nodes, driving the nodes to move and simulate mesh deformation; S4, construct the information transfer matrix based on RBF interpolation theory, calculate the grid domain node temperature in combination with the material domain temperature, and use the updated grid as the initial material domain for the next electrothermal analysis until the total time reaches 30s and the analysis is completed.

[0006] The present invention is also characterized in that: The material properties related to the degree of thermal decomposition described in S1 include electrical conductivity, thermal conductivity, density and specific heat capacity; the boundary conditions are that the four sides and bottom of the composite material are grounded, the top surface is electrically insulated, the top surface and the four sides radiate heat, and the bottom surface is insulated.

[0007] The expression of the lightning current load described in S2 is: (1) Where, t For time.

[0008] The expression of the lightning heat flux load described in S2 is: (2) Where, t It's time, is the radial coordinate, is the radius of the lightning channel.

[0009] The calculation of the pyrolysis degree field, electric potential field and temperature field described in S2 is specifically as follows: the pyrolysis degree field, temperature field and electric potential field are calculated through the pyrolysis reaction kinetic equation, the electric field and temperature field control equations; the pyrolysis reaction heat source or latent heat source is introduced into the temperature field control equation, when the temperature is lower than 3589K, it is a pyrolysis reaction heat source, and when the temperature is higher than 3589K, it is a latent heat source.

[0010] The calculation of the ablation rate and depth described in S2 is specifically as follows: when the node temperature is higher than 700K, the static ablation rate and depth are calculated based on the Hertz-Knudsen equation and the trapezoidal rule. The specific formula is: (6) Where, is the vaporization coefficient, m is the atomic mass of the solid material, is the Boltzmann constant, is the density, is the latent heat of vaporization, and are the boiling temperature and boiling pressure at atmospheric pressure, is the temperature of the corresponding grid at the current time step; go through n The cumulative ablation depth per time increment is calculated based on the static ablation velocity using the trapezoidal rule as follows: (7) Where, t For time, i =1, 2, 3…, v is the calculated static ablation rate, For the k +1 increment time step.

[0011] The S3 is specifically as follows: using the DISP subroutine to convert the static ablation depth setting direction into a displacement load, and driving the mesh domain nodes to move to realize mesh deformation.

[0012] The RBF interpolation basis function described in S4 is a volume spline function: (8) Where, is the Euclidean distance between nodes.

[0013] Compared with the prior art, the present invention has the following beneficial effects: (1) The present invention provides a method for simulating dynamic ablation of composite materials by lightning strike based on ABAQUS. The method calculates the material domain node temperature and thermal decomposition degree of the composite material through electro-thermal analysis, calculates the static ablation rate and depth of the material surface according to the Hertz-Knudsen equation and the trapezoidal rule, and processes the node movement and progressive dynamic removal of the material in the grid domain of the lightning attachment surface during the depression process based on the manual grid movement method of ABAQUS. The RBF interpolation theory is used to achieve the corresponding matching of the material domain and grid domain temperatures after dynamic ablation of the composite material. Finally, the parallel coupling calculation of the composite material heating process and the material surface depression removal process is realized, which provides a feasible method for evaluating the dynamic ablation of composite materials by lightning strike.

[0014] (2) The present invention provides an ABAQUS-based simulation method for dynamic ablation of composite materials caused by lightning strikes. Under the action of a lightning arc, the resin and carbon fibers of the composite material rapidly heat up, undergo pyrolysis and sublimation, and cause ablation on the material surface. Based on electro-thermal analysis, the present invention calculates the material domain node temperature and pyrolysis degree of the composite material caused by heat conduction, Joule heating, and pyrolysis reaction. Combined with the Hertz-Knudsen equation and the trapezoidal rule, the static ablation rate and depth of the material surface are calculated, realizing the simulation of independent heating of the composite material in the material domain and static ablation of the material surface.

[0015] (3) The present invention provides a method for simulating dynamic ablation of composite materials by lightning strikes based on ABAQUS. With the continuous injection of lightning arc energy, the thermal decomposition and sublimation effects of the composite material cause ablation depression along the fiber, matrix, and depth directions on the surface where the lightning strike is attached. In order to simulate the dynamic depression process of the composite material, the present invention uses a manual mesh movement method based on dynamic analysis to calculate the movement of the composite material surface nodes during the depression process. The depression direction is set by the static ablation depth of the material surface obtained by electro-thermal analysis. The mesh domain nodes are driven to move in a step-by-step manner according to the manual mesh movement algorithm. The mesh flow between different layers of the anisotropic composite laminate during the material depression deformation process is simulated, thereby realizing the progressive dynamic removal of the ablated material in the mesh domain.

[0016] (4) The present invention provides a method for simulating dynamic ablation of composite materials by lightning strikes based on ABAQUS. As the composite material continues to heat up and the surface of the lightning strike gradually sinks, the surface material is gradually removed, the ablation sink area continues to deform, and the node positions of the material domain and the grid domain no longer match. To achieve the corresponding matching of the material domain and the grid domain after dynamic ablation of the composite material, the present invention uses RBF interpolation theory to construct an information transfer matrix from the material domain to the grid domain based on the relative position relationship between the nodes in the material domain and the grid domain, and realizes the interpolation mapping of the node temperature in the material domain to the node temperature in the grid domain. As the temperature in the material domain is continuously transferred to the grid domain, the parallel coupling of the composite material heating and the surface sinking removal process is achieved. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is a flow chart of the ABAQUS-based simulation method for dynamic ablation of composite materials caused by lightning strikes.

[0018] Figure 2 This is a schematic diagram of the geometric model and boundary conditions in Example 7 of the ABAQUS-based simulation method for dynamic ablation of composite materials struck by lightning.

[0019] Figure 3 It is a temperature distribution diagram of the simulation results in Example 7 of the ABAQUS-based simulation method for dynamic ablation of composite materials struck by lightning of the present invention.

[0020] Figure 4 This is a damage distribution diagram of the simulation results along the depth direction section in Example 7 of the ABAQUS-based composite material lightning strike dynamic ablation simulation method of the present invention. DETAILED DESCRIPTION

[0021] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. The embodiments described are only some embodiments of the present invention, not all embodiments.

[0022] Example 1 The present invention provides a composite material lightning strike dynamic ablation simulation method based on ABAQUS, such as Figure 1 As shown, the following steps are included: S1. Establish a three-dimensional geometric model of anisotropic composite materials, divide the mesh and define the element type as C3D8R, define the material properties related to the thermal decomposition degree and set the boundary conditions.

[0023] The material properties related to the degree of thermal decomposition in S1 include electrical conductivity, thermal conductivity, density and specific heat capacity; the boundary conditions are that the four sides and bottom of the composite material are grounded, the top surface is electrically insulated, the top surface and the four sides radiate heat, and the bottom surface is insulated.

[0024] S2, apply lightning current and heat flux loads, set the initial values ​​of all variables to the values ​​at room temperature, and calculate the thermal field, electric potential field and temperature field of each grid node in each time analysis step in the electro-thermal analysis module; read the temperature field calculated by the electro-thermal analysis and judge the temperature value. When the node temperature is higher than 700K, calculate the static ablation rate and depth based on the Hertz-Knudsen equation and trapezoidal rule. If the cumulative ablation depth exceeds the set tolerance threshold, terminate the electro-thermal analysis and transmit the calculated static ablation depth to the dynamic analysis module.

[0025] S3 sets the direction of the static ablation depth and converts it into a displacement load on the mesh nodes, driving the nodes to move and simulate mesh deformation.

[0026] S4, based on the node position relationship between the material domain and the grid domain, the RBF interpolation theory is used to construct the information transfer matrix between the material domain and the grid domain. The temperature of the grid domain node is calculated based on the material domain node temperature combined with the information transfer matrix to achieve the corresponding matching of the material domain and grid domain temperatures after dynamic ablation. The updated grid is then used as the initial material domain for the next analysis step of electrothermal analysis, and the analysis is terminated until the total analysis step time reaches 30s.

[0027] Example 2 Based on Example 1, Figure 1 As shown in the figure, S1 is specifically as follows: establish a three-dimensional geometric model of the composite material, according to the unit node definition rules in the ABAQUS three-dimensional finite element model, use [45° / 0° / -45° / 90°] 2s In a symmetric order, define 16 composite plies, each containing 2400 elements. Assign material orientation to each ply using Orientation and Solid Section. Then, set the material properties at different pyrolysis states. Create the assembly and apply boundary conditions and initial conditions at room temperature. The boundary conditions include: radiative heat dissipation from the top and sides, a thermal emissivity of 0.9, an adiabatic bottom, and zero ground potential on the bottom and sides.

[0028] Example 3 Based on Example 2, S2 specifically includes applying lightning current and heat flow loads to the composite material model in step 1 according to the following calculation formula: The lightning current load expression is: (1) Where, t For time.

[0029] The lightning heat flux load expression is: (2) Where, t It's time, is the radial coordinate, is the radius of the lightning channel.

[0030] Then, the electric potential field, temperature field and pyrolysis degree field at each time analysis step are calculated according to the electric field and temperature field control equations and the pyrolysis reaction kinetics equations. The control equations are as follows: The governing equations for the electric field are: (3) Where, V is the grid cell volume; is the electric potential; is the conductivity matrix; S is the grid cell surface; , is the current density, n It's the surface S The normal vector of the surface S The current density entering the grid cell; It is the volume current source per unit volume.

[0031] Temperature field control equation: (4) Where, c 、 、 k Respectively represent the specific heat capacity, density and thermal conductivity of the composite material, v is the ablation rate, t For time, x 、 y and z is the three-dimensional direction, 、 、 They are x 、 y and z Thermal conductivity in three dimensions, is Joule heat, It is a pyrolysis reaction heat source or latent heat source. When the temperature is lower than 3589K, it is a pyrolysis reaction heat source. When the temperature is higher than 3589K, it is a latent heat source. As the radiation heat transfer source, is the heat flux density.

[0032] Pyrolysis reaction kinetic equation: (5) Where, is the thermal decomposition field of the previous time step, n is the reaction order, A is the pre-exponential factor, is the activation energy, R is the gas constant, represents the thermal decomposition field of the current time step, is the temperature of the corresponding grid at the previous time step, is the time step.

[0033] The temperature field of the current analysis step is obtained to determine the temperature of each node. If the node temperature is lower than 700K, the electro-thermal analysis is continued. When the node temperature is higher than 700K, the static ablation rate and depth of the composite material are calculated based on the Hertz-Knudsen equation and the trapezoidal rule. The electro-thermal analysis is terminated when the cumulative ablation depth of the surface nodes exceeds the set tolerance threshold.

[0034] Static ablation speed: (6) Where, is the vaporization coefficient, m is the atomic mass of the solid material, is the Boltzmann constant, is the density, is the latent heat of vaporization, and are the boiling temperature and boiling pressure at atmospheric pressure, is the temperature of the corresponding grid at the current time step.

[0035] go through n time increment (i.e., from the time increment i To time increment n+i- 1, where i =1,2,3…), and use the trapezoidal rule to calculate the cumulative ablation depth: (7) Where, t For time, i =1, 2, 3…, v is the calculated static ablation rate, For the k+ A time step of 1 increment.

[0036] When the cumulative ablation depth of the surface nodes exceeds the set tolerance threshold, the electro-thermal analysis is terminated. At this point, the static ablation depth calculated independently in the electro-thermal analysis is transferred to the dynamic analysis module to provide displacement loads for the subsequent simulation of mesh deformation caused by ablation depression and material removal in the mesh deformation module.

[0037] Example 4 Based on Example 3, S3 is specifically as follows: the static ablation depth calculated in S2 is used as the displacement load for moving the corresponding nodes in the grid domain through the DISP subroutine, and the direction is set at the same time to drive the corresponding nodes in the grid domain to move to the new position. Manual grid movement is achieved by changing the node position, simulating the deformation of the grid.

[0038] Example 5 On the basis of Example 4, Python is called to obtain the deformed mesh from S3 to read the new coordinates of the mesh domain nodes, and the temperature field of the material domain is obtained at the same time. Then, the temperature is mapped to the deformed mesh in Matlab to realize the temperature mapping from the material domain to the mesh domain. In order to improve the computational efficiency, the grid node temperature in the non-ablated area with a static ablation rate and a depth of 0 is directly transferred from the corresponding node temperature in the material domain. For the ablation area with a calculated static ablation rate and a depth greater than 0, due to the ablation removal caused by the increase in the material temperature in the area, the mesh undergoes concave deformation and the node position changes. The grid node temperature in the ablation area is calculated using the volume spline interpolation function method: The RBF interpolation basis function is a volume spline function: (8) Where, is the Euclidean distance between nodes.

[0039] First, construct a linear equation system based on the material domain node information Calculate weight coefficient : Construct volume spline interpolation basis functions based on the spatial coordinates of the material domain nodes The basis function matrix M With the polynomial matrix P , combined into an extended matrix At the same time, the right-hand matrix of the linear equations is constructed according to the temperature field of the material domain node , as follows: Basis function matrix M : (8) in, N S is the number of mesh nodes in the material domain, Any two points in the material domain mesh node i 、 j (i , j =1,2,… N S ), the basis function of volume spline interpolation is , is the Euclidean distance between two points, i 、 j The two nodes (space coordinates are ( x i , y i , z i )、( x j , y j , z j The Euclidean distance calculation formula is: (9) Polynomial Matrix P : (10) Right matrix : (11) in, is the temperature of the material domain mesh node: (12) The weight coefficients can be calculated based on the linear equations Obtain: (13) Then construct the evaluation matrix based on the coordinates of the material domain nodes and the grid domain ablation nodes : (14) in, N f is the number of mesh nodes in the ablation area within the mesh domain.

[0040] According to the obtained weight coefficient and evaluation matrix , calculate the information transfer matrix from the material domain to the mesh domain H : (15) Calculate the information transfer matrix between the two domains, where the transferable information includes temperature, so the information transfer matrix can be used to calculate the information transfer matrix. H The temperature of the mesh domain nodes is calculated using the material domain mesh node temperature: (16) After completing the mapping of the temperatures of the nodes in the ablated and unablated areas from the material domain to the mesh domain, the obtained dynamic mesh results are passed to the electro-thermal analysis module as the initial step for the next time step. Each time the electro-thermal analysis starts, the field variables are calculated and mapped to the dynamic mesh temperature, and the results are passed to the electro-thermal module again, which constitutes a cycle until the total analysis step time reaches 30s and the analysis is terminated.

[0041] Example 6 The ABAQUS-based composite material lightning strike dynamic ablation simulation method provided in Example 5 is used, and the specific process is as follows: S1. Divide the mesh in ABAQUS and set the element type to C3D8R. The composite material model consists of 16 layers with 2400 elements per layer, and the mesh is [45° / 0° / -45° / 90°]. 2s The layers are laid in a symmetrical order and then the electric and temperature field boundary conditions are defined.

[0042] S2, set the initial value of the variable to the value at room temperature, and apply lightning current and heat flow loads in ABAQUS: Lightning current load: (1) Where, t For time.

[0043] Lightning heat flux load: (2) Where, t It's time, is the radial coordinate, is the radius of the lightning channel.

[0044] Composite materials heat up rapidly under the action of lightning current load and heat flow load, undergoing pyrolysis and sublimation, resulting in ablation of the material surface. The static ablation rate and depth of the composite material surface are calculated by the user-defined USDFLD subroutine in the finite element software ABAQUS: the pyrolysis degree and temperature of the material domain mesh node after electro-thermal analysis are read. When the node temperature is lower than 700K, the static ablation rate is 0. When the node temperature is higher than 700K, the static ablation rate is calculated based on the node temperature, expressed as: (6) Where, is the vaporization coefficient, m is the atomic mass of the solid material, is the Boltzmann constant, is the density, is the latent heat of vaporization, and are the boiling temperature and boiling pressure at atmospheric pressure, is the temperature of the corresponding grid at the current time step.

[0045] go through n time increment (i.e., from the time increment i To time increment n+i -1, where i =1,2,3…), the cumulative ablation depth calculated by the USDFLD subroutine is: (7) Where, t For time, i =1, 2, 3…, v is the calculated static ablation rate, For the k +1 increment time step.

[0046] When the cumulative ablation depth reaches the specified tolerance threshold ( h >Tol), terminate the electro-thermal analysis using the LSTOP flag in the URDFIL subroutine.

[0047] To achieve temperature matching between the material and mesh domains after dynamic ablation, temperature mapping of the deformed mesh is required. The temperature of mesh nodes that have not undergone ablation remains unchanged, while the temperature of mesh nodes that have undergone ablation is calculated using RBF interpolation. The results of the temperature mapping at each time step are used to initiate the electrothermal module analysis in the next time step, and the analysis continues until the total analysis step time reaches 30 seconds.

[0048] S3, after the electrothermal analysis is terminated, the static ablation depth is applied as a displacement with the same mesh using the DISP subroutine and the direction is set to drive the surface nodes in the mesh domain to move to the new position, simulating the deformation of the mesh by changing the coordinate position of the mesh domain nodes.

[0049] After mesh deformation, to ensure a consistent match between the mesh and material domains after dynamic ablation, RBF interpolation is used in Matlab to perform temperature mapping from the material domain to the mesh domain. First, node information for the deformed mesh domain is obtained. Ablated and unablated regions are distinguished based on whether the material domain nodes have ablated. Node temperatures in unablated regions are directly transferred from the material domain to the corresponding nodes in the mesh domain. Node temperatures in ablated regions are calculated using RBF interpolation. An information transfer matrix is ​​constructed between the material and mesh domains based on their relative positions. The mesh domain node temperatures are then inferred from the material domain node temperatures. The resulting temperature mapping is then output as the material domain, serving as the initial step for the next electrothermal analysis.

[0050] Example 7 A three-dimensional geometric model of the composite material was established according to the unit node definition rules in the ABAQUS three-dimensional finite element model, and the material orientation of each ply was specified through Orientation and Solid Section. Subsequently, the performance parameters of the material under different pyrolysis states were set. Next, an assembly was established and boundary conditions and an initial temperature field at room temperature were applied. The boundary conditions were: heat dissipation by radiation from the upper surface and side surfaces, a thermal radiation coefficient of 0.9, a bottom surface insulated, and the bottom and side surfaces were grounded, the electric potential was 0, and the initial conditions were set at room temperature of 25°C.

[0051] The model established in ABAQUS is as follows Figure 2 As shown, the model size is 150mm×100mm×2mm, with a total of 16 layers, each layer contains 2400 units, each layer is 0.125mm thick, and the layers are arranged in [45° / 0° / -45° / 90°] 2s The layers are laid in a symmetrical order. After modeling is completed, the material properties are set under the initial conditions, including the electrical conductivity, thermal conductivity, density and specific heat capacity of the composite material. These parameters change with the pyrolysis state. The reference performance parameters are shown in Table 1.

[0052] Table 1

[0053] The lightning current component D is introduced during the load application process through surface current. The current load is set using the Analytical field and Amplitude functions, and the time-varying heat flux load is defined using the user-defined DFLUX subroutine interface. Under the analytical calculations of the electric and temperature fields, the custom subroutine HETVAL is used to introduce a reaction heat source term into the temperature field governing equations to simulate the heat released by the pyrolysis reaction during the lightning strike. After the pyrolysis reaction is complete, a latent heat source is introduced to prevent further temperature rise in the material.

[0054] Initial field variables are updated, and the static ablation rate and depth are calculated using the user-defined subroutine USDFLD. In the electro-thermal analysis module, the user-defined subroutine USDFLD reads the electro-thermal analysis temperature field. The temperature of each mesh node is then determined. If the node temperature is below 700K, the ablation rate is zero. When the node temperature exceeds 700K, the static ablation rate and depth of the composite material are calculated using the Hertz-Knudsen equation and the trapezoidal rule. The electro-thermal analysis is terminated when the cumulative ablation depth at the surface nodes reaches the set tolerance threshold. The static ablation depth calculated from the electro-thermal analysis is transferred to the dynamic analysis module, where manual mesh movement is used to simulate progressive dynamic material removal. Finally, RBF interpolation is used to map the temperature from the material domain to the mesh domain, enabling parallel coupled calculations of the composite material heating and surface dent removal processes. After the temperature mapping from the material domain to the mesh domain is completed, the meshed results are transferred to the electro-thermal analysis module to be used as the material domain for the next step. The analysis process from the electro-thermal analysis of the material domain to the update of the material domain after the dynamic mesh is a complete analysis process. The calculation and analysis are repeated until the total analysis step time reaches 30s and the analysis is terminated.

[0055] like Figure 3 、 4 As shown in the figure, the temperature distribution diagram of the simulation results in the heat dissipation stage and the damage distribution diagram along the depth direction section are respectively. As the heat dissipation stage accumulates in time, the temperature in the composite laminate gradually decreases, which is manifested as the reduction of the high temperature area and the expansion of the low temperature area. The temperature in the plate gradually approaches the starting temperature of resin ablation, and the damage in the surface and depth tends to be stable. Figure 3 The temperature distribution diagram can intuitively show the basically stable temperature distribution of the composite laminate damage, which can be compared with the experiment to verify the accuracy of the simulation model; Figure 4 As shown in the figure, from the damage distribution diagram along the depth direction, it can be clearly observed that the grid deformation caused by the ablation damage produces a depression in the center of the composite laminate, making the simulation of dynamic ablation more consistent with the real process.

Claims

1. The ABAQUS-based simulation method for dynamic ablation of composite materials caused by lightning strikes is characterized by: The specific steps are: S1, establish a three-dimensional geometric model of anisotropic composite materials, divide the mesh, define material properties related to pyrolysis degree and set boundary conditions; S2, applies lightning current and heat flux loads, sets the initial values ​​of variables to those at room temperature, calculates the thermal decomposition degree, electric potential, and temperature field of each node, and calculates the ablation rate and depth when the temperature is higher than 700K. The electrothermal analysis is terminated when the cumulative depth exceeds the threshold. S3 sets the direction of the static ablation depth and converts it into a displacement load on the mesh nodes, driving the nodes to move and simulate mesh deformation; S4, construct the information transfer matrix based on RBF interpolation theory, combine the temperature to calculate the grid domain node temperature, and use the updated grid as the initial material domain for the next electrothermal analysis. The analysis ends when the total time reaches 30s.

2. The ABAQUS-based composite material lightning strike dynamic ablation simulation method according to claim 1 is characterized in that: The material properties related to the degree of thermal decomposition described in S1 include electrical conductivity, thermal conductivity, density and specific heat capacity; the boundary conditions are that the four sides and bottom of the composite material are grounded, the top surface is electrically insulated, the top surface and the four sides radiate heat, and the bottom surface is insulated.

3. The ABAQUS-based composite material lightning strike dynamic ablation simulation method according to claim 1, characterized in that: The expression of the lightning current load described in S2 is: （1） Where, t For time.

4. The ABAQUS-based composite material lightning strike dynamic ablation simulation method according to claim 1, characterized in that: The expression of the lightning heat flux load described in S2 is: （2） Where, t It's time, is the radial coordinate, is the radius of the lightning channel.

5. The ABAQUS-based composite material lightning strike dynamic ablation simulation method according to claim 1, characterized in that: The calculation of the pyrolysis degree field, electric potential field and temperature field described in S2 is specifically as follows: the pyrolysis degree field, electric potential field and temperature field are calculated through the pyrolysis reaction kinetic equation, the electric field and temperature field control equations; the pyrolysis reaction heat source or latent heat source is introduced into the temperature field control equation, when the temperature is lower than 3589K, it is a pyrolysis reaction heat source, and when the temperature is higher than 3589K, it is a latent heat source.

6. The ABAQUS-based composite material lightning strike dynamic ablation simulation method according to claim 1, characterized in that: The calculation of the static ablation rate and depth described in S2 is as follows: When the node temperature is higher than 700K, the static ablation rate and depth are calculated based on the Hertz-Knudsen equation and the trapezoidal rule. The specific formula is: （6） Where, is the vaporization coefficient, m is the atomic mass of the solid material, is the Boltzmann constant, is the density, is the latent heat of vaporization, and are the boiling temperature and boiling pressure at atmospheric pressure, is the temperature of the corresponding grid at the current time step; go through n The cumulative ablation depth per time increment is calculated based on the ablation velocity using the trapezoidal rule as follows: （7） Where, t For time, i =1, 2, 3…, v is the calculated static ablation rate, For the k +1 increment time step.

7. The ABAQUS-based composite material lightning strike dynamic ablation simulation method according to claim 1, characterized in that: The S3 is specifically as follows: using the DISP subroutine to convert the static ablation depth setting direction into a displacement load, and driving the mesh domain nodes to move to realize mesh deformation.

8. The ABAQUS-based composite material lightning strike dynamic ablation simulation method according to claim 1, characterized in that: The RBF interpolation basis function described in S4 is a volume spline function: （8） Where, is the Euclidean distance between nodes.

Images