C / C composite material ablation morphology prediction method considering complex microstructure

Abstract

The invention provides a C / C composite material ablation morphology prediction method considering a complex microstructure, and solves the technical problem that the prediction result is not accurate enough due to the fact that the influence of the material microstructure on the ablation morphology is not fully considered in the existing prediction method. A random sequence adsorption algorithm is adopted to create C / C composite material microscopic representative volume cell elements based on a three-dimensional voxel grid, and under a finite element transient heat transfer analysis framework, ablation criteria of all components of the C / C composite material are respectively established according to temperature information and chemical reaction kinetic parameter information; and finally, C / C composite material microscopic ablation morphology prediction comprehensively considering complex microstructure characteristics is realized.

Description

Technical Field

[0001] The invention belongs to the technical field of numerical simulation of ablation behavior of C / C composite materials, and in particular relates to a method for predicting ablation morphology of C / C composite materials taking into account complex microstructural characteristics. Background Art

[0002] Carbon fiber reinforced carbon-based composites (abbreviated as C / C composites) have the characteristics of low density, small thermal expansion coefficient, good high-temperature mechanical properties and ablation resistance. They are the primary candidate materials for hot-end components in the fields of aviation, aerospace and nuclear energy. Under the scouring of high-temperature oxygen-containing gas flow, oxidation ablation will occur on the surface of C / C composites, affecting their service life. In order to predict the ablation behavior of C / C composites, a variety of macroscopic ablation models have been developed. Macroscopic theory believes that C / C composites are homogeneous materials and ignores the influence of microstructural characteristics on the ablation properties of materials. However, in fact, the ablation resistance of carbon fibers and matrices in C / C composites is significantly different. The performance difference between the components will lead to complex micro-morphology on the surface of the material. As the ablation morphology continues to evolve, the aerodynamic thermal load on the surface of the material will also change, thereby affecting the ablation behavior of the material. Further studies have shown that the microstructural characteristics of thermal protection materials, such as pore content and fiber arrangement, also have an important influence on their ablation properties. Therefore, studying the ablation morphology of C / C composites at a microscopic scale is helpful to further reveal the ablation mechanism of the material and lay the foundation for the reliability design and service life assessment of the hot end components of C / C composites.

[0003] At present, the methods for predicting the microscopic ablation morphology of thermal protection materials mainly include the level set method and the fluid volume method. Chinese patent CN113722830B discloses a method for predicting the microscopic ablation morphology of a C / C composite nozzle of a solid rocket engine. This method uses Fick's diffusion theory and additional convection mechanism to construct a dilute species transfer model, and successfully tracks the microscopic ablation morphology of unidirectional C / C composite materials using the level set method. Reference 1 (Acta Mechanica Sinica, 2019, 51(03): 835-844) introduces the solid phase volume fraction and uses the fluid volume method to establish a reaction-diffusion model to describe thermochemical ablation, and solves the microscopic ablation morphology of unidirectional C / C composite materials while ignoring the influence of pores. The above method involves solving complex level set equations or reaction-diffusion coupling equations, and its computational efficiency is limited to a certain extent. It is more suitable for predicting the ablation morphology of thermal protection materials with relatively simple microstructures. However, C / C composites not only include carbon fiber preforms with complex morphology (such as chopped carbon fiber preforms, needle-punched carbon fiber preforms, and three-dimensional woven preforms), but also have a large number of randomly distributed pores inside. These microstructural features will have an important influence on the ablation morphology. Therefore, there is an urgent need for a method to predict the ablation morphology of C / C composites that can comprehensively consider the influence of complex microstructural features. Summary of the invention

[0004] The purpose of the present invention is to solve the technical problem that the existing prediction methods do not fully consider the influence of the material microstructure on the ablation morphology, resulting in inaccurate prediction results, and to provide a C / C composite material ablation morphology prediction method taking into account the complex microstructure.

[0005] To achieve the above purpose, the technical solution provided by the present invention is:

[0006] A method for predicting the ablation morphology of C / C composite materials considering complex microstructures is special in that it includes the following steps:

[0007] Step 1: Establish a transient heat transfer model of C / C composites based on microscale

[0008] Step 1.1: Build the geometry and mesh it

[0009] A cube model with the size of L×L×L was constructed for the C / C composite material;

[0010] The sweep algorithm is used to divide the L×L×L cube model into hexahedral meshes to form NUM tol units to realize the discretization of L×L×L cube geometric model;

[0011] Read the number of each unit and its center point coordinates;

[0012] Step 1.2: Determine the component type to which each unit belongs. Use a random sequence adsorption algorithm to generate the coordinates of the carbon fiber axis endpoints and the pore axis endpoints within the cubic model space. Combine the number of each unit obtained in step 1.1 and its center point coordinates, and calculate the minimum distance between the center point of each unit and all carbon fiber axes in turn. When the minimum distance is less than or equal to the carbon fiber radius, the unit belongs to the carbon fiber component;

[0013] Otherwise, the minimum distance between the unit and all pore axes is calculated. When the minimum distance is less than or equal to the pore radius, the unit belongs to the pore component;

[0014] Otherwise, the unit belongs to the matrix component;

[0015] Step 1.3: According to the component types of each unit determined in step 1.2, a carbon fiber unit set, a pore unit set and a matrix unit set are established, and corresponding material properties are assigned to each unit set. A heat flux boundary condition is applied to the initial ablation surface. Based on this, a transient heat transfer model of C / C composite materials based on microscale is established;

[0016] Step 2: Perform transient heat transfer analysis and determine whether ablation occurs in each unit of the C / C composite cube model

[0017] Initialize the current increment step number K and the unit number N to 1, summarize the numbers of the units corresponding to the initial ablation surface, and form a surface unit number set under the current increment step; wherein the surface unit refers to the unit to which the surface directly contacts the external oxidizing gas belongs;

[0018] Step 2.1: Execute the Kth incremental step of the transient heat transfer model of C / C composite materials based on microscale obtained in step 1, and read the Nth unit; according to the component type of the current unit, determine whether the current incremental step number K is KINC i An integer multiple of (when it is an integer multiple, it proves that the unit has been completely ablated, and the unit can be considered to have disappeared by further assignment; when it is a non-integer multiple, the unit is not completely ablated, so it is necessary to wait for it to be completely ablated in the subsequent incremental step before further processing), if yes, go to step 2.2; otherwise, determine whether N is equal to the total number of units NUM tol ; When N is less than NUM tol When N is equal to NUM tol When K is equal to the total increment number K tol ; When K is equal to K tol When , go to step 3; otherwise, define K = K + 1 and N = 1, update the surface unit number set and re-execute step 2.1;

[0019] Among them, KINC i The number of incremental steps required to completely burn through a cell is calculated by KINC i The values ​​reflect the differences in linear ablation rates of the components of C / C composites;

[0020] Step 2.2: Determine whether the current unit N is a surface unit according to the surface unit number set under the current incremental step number;

[0021] If it is a surface unit, the value of the current unit state variable 2 is defined as 1, otherwise, the value of the current unit state variable 2 is defined as 0;

[0022] Step 2.3: Read the temperature value T of the current unit E , judge T E Is it greater than or equal to the critical ablation temperature T w,i ;

[0023] If T E Greater than or equal to T w,i , indicating that ablation occurs, the value of the current unit state variable 1 is defined as 1; otherwise, indicating that ablation does not occur, the value of the current unit state variable 1 is defined as 0;

[0024] Step 2.4: Determine whether the value of the current unit state variable 1 and the value of the state variable 2 are both 1;

[0025] If the value of the current unit state variable 1 and the value of the state variable 2 are both 1, the number of the current unit is recorded in the set of ablated unit numbers corresponding to the current incremental step number, and step 3 is entered;

[0026] Otherwise, go to step 2.5;

[0027] Step 2.5: Determine whether N is equal to NUM tol ; When N is less than NUM tol When N is equal to NUM tol When K is equal to K tol , when K is equal to K tol , go to step 3; otherwise, define K=K+1 and N=1, update the surface unit number set, and re-execute step 2.1;

[0028] The entire step 2 performs surface judgment and ablation judgment on all units at each incremental step, fully considering the influence of the material microstructure on the ablation morphology and ensuring the accuracy of the model;

[0029] Step 3: Predict the microscopic ablation morphology of C / C composites

[0030] The post-processing function of the finite element analysis software is used to hide the ablated units corresponding to each incremental step obtained in step 2 (i.e., only the non-ablated units are displayed) to obtain the microscopic ablation morphology of the C / C composite material at different times.

[0031] Furthermore, in order to evaluate the anti-ablation performance of the material and verify the accuracy of the ablation morphology prediction results of this method, the following are also included:

[0032] Step 4: Calculate the ablation depth

[0033] In ablation experiments, the maximum height or thickness change of a material within a certain period of time is usually defined as the ablation depth, and the average ablation depth value measured multiple times is used to evaluate the material's anti-ablation performance. Therefore, using the ablation results of the C / C composite material transient heat transfer model generated in step 2, the ablation depth s of the C / C composite material at different times is calculated according to the following formula: t :

[0034]

[0035] Where NUM ab represents the number of cells whose state variable 1 is 1 in the transient heat transfer model of C / C composite materials, EVOL represents the cell volume, and A inirepresents the initial ablated surface area.

[0036] Furthermore, in step 1.2, the distance between the carbon fiber axis / pore axis and the center point of each unit is calculated using the following formula:

[0037]

[0038] In the formula, x p ,y p and z p is the coordinate of the carbon fiber axis or the end point P of the pore axis; q ,y q and z q is the coordinate of the end point Q of the carbon fiber axis or pore axis; c ,y c and z c is the coordinate of the center point C of each unit; the line segment CN intersects the line segment PQ perpendicularly at point N(x n ,y n ,z n ).

[0039] Furthermore, in step 1.3, according to the component type of each unit determined in step 1.2, the corresponding density, thermal conductivity and specific heat capacity are assigned to different units, the analysis step type is defined as transient heat transfer analysis, and a transient heat transfer model of C / C composite materials based on microscale is established.

[0040] Furthermore, in step 2.1, the number of incremental steps required to completely burn through a unit is calculated using the following formula:

[0041]

[0042] In the formula, H ele Indicates the unit thickness; inc is the incremental step size of heat transfer analysis; v D,i is the linear ablation rate of carbon fiber and matrix of C / C composite;

[0043] When the temperature is low, the ablation behavior of the material is affected by the oxidation reaction kinetic parameters, according to the Arrhenius equation:

[0044]

[0045] In the formula, k represents the chemical reaction rate constant, T represents the temperature, A represents the pre-exponential factor, E represents the activation energy of the oxidation reaction, and R represents the ideal gas constant. As the temperature increases, the oxidation ablation rate of the material increases rapidly in an exponential manner. When the temperature exceeds the critical value, the ablation rate of the material remains unchanged due to the limitation of the diffusion rate of the oxidizing gas. Therefore, the ablation behavior of the material below the critical temperature can be ignored. When the temperature exceeds the critical value, it is considered that the material ablates at a constant rate. In summary, the linear ablation rate of carbon fiber and matrix of C / C composite materials is calculated as follows:

[0046]

[0047] Where ρ represents density, C p is the specific heat capacity at constant pressure, represents the relative molecular mass of oxygen, h c represents the convective heat transfer coefficient, M C represents the relative molecular mass of carbon, represents the mass fraction of oxygen, η represents the stoichiometric coefficient of the oxidation reaction equation, and the subscript i=f, m represents the component type of the C / C composite material, f represents the carbon fiber, and m represents the matrix.

[0048] Furthermore, in step 2.4, the critical ablation temperature of the carbon fiber and the matrix of the C / C composite material is calculated using the following formula:

[0049]

[0050] In the formula, ρ represents density, M O represents the relative molecular mass of oxygen atoms, A represents the pre-exponential factor, E represents the activation energy of the oxidation reaction, R represents the ideal gas constant, P W Indicates pressure.

[0051] Furthermore, in step 2, the specific steps of updating the surface unit number set are as follows:

[0052] When K is greater than 1, the set of ablated unit numbers corresponding to the previous incremental step is read;

[0053] Since the unit numbering of the C / C composite material cube model discretized by the sweep algorithm satisfies the following rule: the difference in unit numbers at the corresponding positions of two adjacent layers of units is equal to the total number of units in each layer NUM lay , therefore, the surface unit number of the current incremental step is calculated using the following formula:

[0054]

[0055] In the formula, represents the surface unit number of step K, Indicates the ablated cell number of step K-1.

[0056] Furthermore, step 1.1, step 1.3, step 2 and step 3 all use Abaqus software; wherein step 1.3 uses the Usdfld user subroutine in the Abaqus software;

[0057] Steps 1.2, 4 and updating the surface unit number are all performed using Matlab software.

[0058] Furthermore, in step 1.1, the value range of L is 0.1 mm to 1 mm.

[0059] Furthermore, in step 2.1, the value of inc is 0.01s to 1s.

[0060] Principle of the present invention:

[0061] In order to comprehensively consider the influence of complex microstructural characteristics on the ablation morphology of C / C composites, the random sequence adsorption algorithm was used to create a representative microscopic volume cell of C / C composites based on a three-dimensional voxel grid. Under the framework of finite element transient heat transfer analysis, the ablation criteria of each component of the C / C composite were established according to the temperature information and chemical reaction kinetic parameter information, and finally the prediction of the microscopic ablation morphology of C / C composites with comprehensive consideration of the complex microstructural characteristics was achieved.

[0062] Beneficial effects of the present invention:

[0063] 1. The present invention creates a representative volume cell of C / C composite materials taking into account the preform type and random distribution of pores. By establishing the ablation judgment criteria of each component (carbon fiber and matrix), the microscopic ablation morphology of C / C composite materials is predicted, which provides a reference for modeling the ablation behavior of other thermal protection materials or coatings with complex microstructures.

[0064] 2. The present invention helps to deeply explore the influence of microstructural characteristic parameters on the ablation behavior of C / C composite materials, further reveal the microscopic ablation mechanism of C / C composite materials, and lay the foundation for the local on-demand optimization design of hot end components of C / C composite materials with high ablation resistance. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 Create a flow chart for each component unit set of a representative volume cell of a C / C composite material;

[0066] Figure 2 It is a collection of component units of the microscopic representative volume cell of the chopped C / C composite material; among them, (a) is a collection of carbon fiber units, (b) is a collection of pore units, and (c) is a collection of matrix units;

[0067] Figure 3 The Usdfld user subroutine flow chart for ablation determination of each component of C / C composite material;

[0068] Figure 4 The algorithm flow chart for updating the surface unit number set of C / C composite materials for ablation process;

[0069] Figure 5 Schematic diagram of the unit numbering of the C / C composite material cube model discretized based on the sweeping algorithm;

[0070] Figure 6 This is the ablation depth-time curve of short-cut C / C composite material under oxygen-acetylene flame flow;

[0071] Figure 7 Schematic diagram of the microscopic morphology of short-cut C / C composite materials after ablation at different times (blue represents carbon fiber and green represents matrix); among them, (a) is after 20s, (b) is after 40s, and (c) is after 60s. DETAILED DESCRIPTION

[0072] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments:

[0073] In order to predict the microscopic ablation morphology of C / C composite materials, the present invention proposes a method for predicting the ablation morphology of C / C composite materials taking into account complex microstructures. In order to enable relevant technical personnel to better understand the embodiments of the present invention, the ablation process of short-cut C / C composite materials under oxygen-acetylene flame flow is taken as an example, in which the volume fractions of carbon fibers and pores are 10.6% and 8.5%, respectively. Figure 1-7 The method for predicting the ablation morphology of C / C composite materials considering the complex microstructure characteristics proposed in the present invention comprises the following steps:

[0074] Step 1: Establish a transient heat transfer model for chopped C / C composites based on microscale

[0075] Step 1.1: Create a cube model with a size of 0.4mm×0.4mm×0.4mm for the chopped C / C composite material in the Abaqus / Part module, and discretize the geometric model by the sweeping algorithm (i.e., divide the cube model into hexahedral meshes). After discretization, all unit sizes of the model are 0.005mm×0.005mm×0.005mm, and the total number of units is 512000. Use the getMassProperties function to read the number of each unit and its center point coordinates, and output them to the Unit-CenterCoords.txt file.

[0076] Step 1.2: Assuming that the length and diameter of the carbon fiber are 0.2 mm and 0.05 mm respectively, simplify the pore into a sphere with a diameter of 0.03 mm. Use the random sequence adsorption algorithm in Matlab software to generate the coordinates of all carbon fiber axis endpoints and pore axis endpoints within the cubic model space. Figure 1 As shown, read the generated Unit-CenterCoords.txt file, calculate the minimum distance between each unit center point and all carbon fiber axes, when the minimum distance is less than or equal to the carbon fiber radius, the unit belongs to the carbon fiber component; otherwise, calculate the minimum distance between the unit and all pore axes, when the minimum distance is less than or equal to the pore radius, the unit belongs to the pore component, otherwise it belongs to the matrix component. By traversing all units, the following can be achieved: Figure 2 The unit sets of each component of the C / C composite material shown in (a)-(c) are created. Among them, the distance d between the carbon fiber axis or the pore axis and the center point of each unit is calculated using formulas (1)-(3):

[0077]

[0078] In the formula, x p ,y p and z p is the coordinate of the end point P of the carbon fiber axis or pore axis; q ,y q and z q is the coordinate of the end point Q of the carbon fiber axis or pore axis; c ,y c and z c is the coordinate of the unit center point C; the line segment CN intersects the line segment PQ perpendicularly at point N(x n ,y n ,z n );

[0079] Step 1.3: Define the density, thermal conductivity and specific heat capacity of the chopped C / C composite carbon fiber in the Abaqus / Property module as 1.8 g / cm 3 , 9.4W / m / K and 820J / kg / K; the density, thermal conductivity and specific heat capacity of the chopped C / C composite matrix are defined as 1.5g / cm 3 , 135.68W / m / K and 710J / kg / K; the density, thermal conductivity and specific heat capacity of the pores are defined as 1.29g / cm 3, 0.025W / m / K and 1005J / kg / K. In the Abaqus / Step module, the analysis type is defined as transient heat transfer analysis, the total analysis time is 80s, and the incremental step size is 0.1s. The unit type is set to DC3D8R, the initial temperature field of the short-cut C / C composite material is defined as 25℃, and the heat flux boundary condition is applied to the surface directly acted by the oxygen-acetylene flame flow (i.e., the initial ablation surface), where the heat flux density value is 4.28MW / m 2 Based on this, a transient heat transfer model of C / C composite materials based on microscale is established.

[0080] Step 2: Calculate the linear ablation rate and critical ablation temperature of chopped C / C composite carbon fibers and matrix

[0081] The control mechanisms of oxidation ablation rate of C / C composites under oxygen-acetylene flame flow include: chemical reaction kinetic control mechanism and oxidizing gas diffusion control mechanism. When the temperature is low, the ablation behavior of C / C composites is affected by the oxidation reaction kinetic parameters. According to the Arrhenius formula:

[0082]

[0083] Where k represents the chemical reaction rate constant, T represents temperature, A represents the pre-exponential factor, E represents the activation energy of the oxidation reaction, and R represents the ideal gas constant. As the temperature increases, the oxidation ablation rate of the material increases rapidly in an exponential manner; when the temperature exceeds the critical value, the ablation rate of the material remains unchanged due to the limitation of the diffusion rate of the oxidizing gas. Therefore, the ablation behavior of the material below the critical temperature is ignored, and when the temperature exceeds the critical value, it is assumed that the material ablates at a constant rate. The linear ablation rate and critical ablation temperature of the carbon fiber and matrix of the C / C composite material are calculated according to formulas (5) and (6):

[0084]

[0085] In the formula, ρ represents density, h c represents the convection heat transfer coefficient, C p is the specific heat capacity at constant pressure, represents the mass fraction of oxygen, η represents the stoichiometric coefficient of the oxidation reaction equation, Represents the relative molecular mass of oxygen, M O Represents the relative molecular mass of oxygen atoms, M C Indicates the relative molecular mass of carbon, P W represents pressure, and the subscript i (i = f, m) represents the component type of the C / C composite material, where f is the carbon fiber and m is the matrix.

[0086] During the experiment, the C / C composite material mainly reacted chemically with the oxygen in the oxy-acetylene flame. The pre-exponential factor and activation energy of the reaction between carbon fiber and oxygen were 100 and 1.2×10 5 J / mol, the pre-exponential factor and activation energy of the substrate and oxygen reactions are 1000 and 1.2×10 5 J / mol. After calculation, the linear ablation rate and critical ablation temperature of carbon fiber are 5.9×10 -3 mm / s and 540℃, the linear ablation rate and critical ablation temperature of the substrate are 7.1×10 -3 mm / s and 497°C.

[0087] At the same time, the number of incremental steps required to completely burn through a unit is calculated based on formula (7):

[0088]

[0089] In the formula, H ele represents the unit thickness, inc represents the incremental step size of the heat transfer analysis, and KINC i The values ​​reflect the differences in linear ablation rates of the components of the C / C composite material. The incremental steps corresponding to burning through a single carbon fiber unit, matrix unit and pore unit are finally defined as 9, 7 and 1.

[0090] The purpose of this step is to provide available data for the subsequent step 3. There is no obvious time sequence between this step and step 1. This step can be performed simultaneously with step 1 or before step 1.

[0091] Step 3: If Figure 3 As shown in the figure, a transient heat transfer analysis is performed and the Abaqus / Usdfld user subroutine is used to determine whether ablation occurs in each component unit of the C / C composite material.

[0092] Initialize the current incremental step number K and unit number N to 1, define the unit on the initial heat flux input surface as the surface unit, and output the corresponding unit number to the Surf_K.txt file, where the surface unit refers to the unit belonging to the surface that is in direct contact with the external oxidizing gas.

[0093] Step 3.1: Execute the Kth incremental step of the transient heat transfer model of the C / C composite material based on the micro-scale obtained in step 1, and read the Nth unit. According to the component type of the current unit, determine whether the current incremental step number K is KINC i If it is an integer multiple, then go to step 3.2; otherwise, determine whether N is equal to the total number of units NUM of the model. tol ; When N is less than NUM tol When N is equal to NUMtol When K is equal to the total increment number K tol ; When K is equal to K tol If , then go to step 4, otherwise, define K = K + 1 and N = 1, update the surface unit number set, and re-execute step 3.1;

[0094] like Figure 4 As shown, the specific steps for updating the surface unit number set are as follows:

[0095] When K is equal to 1, there is no need to update the surface unit number set; when K is greater than 1, the Ab_K.txt file generated by the previous incremental step is read based on the Matlab software. Since the model unit numbering implemented by the sweep algorithm for discretization satisfies the following rules ( Figure 5 ): The difference between the unit numbers at the corresponding positions of two adjacent layers is equal to the total number of units NUM in each layer lay . Using the formula:

[0096]

[0097] Calculate the surface unit number of the current incremental step, and output the updated surface unit number set to the Surf_K.txt file; where, represents the surface unit number of step K, Indicates the ablated cell number of step K-1.

[0098] Step 3.2: Determine whether the current unit N is a surface unit according to the surface unit number set under the current incremental step number;

[0099] Based on the function UEXTERNALDB, the generated Surf_K.txt file is read to determine whether the current unit is a surface unit. If it is a surface unit, the state variable STATEV(2) is defined as 1, otherwise the state variable STATEV(2) is defined as 0.

[0100] STATEV(2) is a flag indicating that the current cell belongs to a surface cell set. When STATEV(2)=1, it indicates that the current cell is a surface cell. When STATEV(2)=0, it indicates that the current cell is a non-surface cell.

[0101] Step 3.3: Read the temperature value T of the current unit E , judge T E Is it greater than or equal to the critical ablation temperature T w,i ; Use the function GETVRM to read the temperature value T of the current unit E , judge T E Is it greater than or equal to the critical ablation temperature T w,i , such as T EGreater than or equal to T w,i Then define the value of the state variable STATEV(1) as 1, otherwise define the value of the state variable STATEV(1) as 0.

[0102] STATEV(1) is a flag indicating that the current unit has been ablated. When STATEV(1)=1, it indicates that the current unit has been ablated. When STATEV(1)=0, it indicates that the current unit has not been ablated.

[0103] Step 3.4: Determine whether the STATEV(1) value and STATEV(2) value of the current cell are both 1. If so, write the current cell number into the Ab_K.txt file. Otherwise, go to step 3.5.

[0104] Step 3.5: Determine whether N is equal to NUM tol ; When N is less than NUM tol When N is equal to NUM tol When K is equal to K tol If yes, proceed to step 4, otherwise, define K=K+1 and N=1, update the surface unit number set, and re-execute step 3.1.

[0105] Step 4: Obtain the ablation morphology of C / C composites at different times

[0106] When the incremental step K of the transient heat transfer model reaches the total incremental number K tol Finally, open the transient heat transfer model result file generated in step 3 in the Abaqus / CAE module, select the state variable STATEV(1) in the field variable, and the part where the value of STATEV(1) is equal to 0 is the morphology after ablation.

[0107] Step 5: Calculate the ablation depth

[0108] In the ablation experiment, the ablation resistance of the material is generally evaluated by repeatedly measuring the maximum height or thickness change (i.e., ablation depth) of the material within a certain period of time. Furthermore, the transient heat transfer model result file generated in step 3 is read using a Python script, and the ablation depth s of the C / C composite material at different times is calculated using formula (9): t :

[0109]

[0110] Where NUM ab represents the total number of cells whose state variable STATEV(1) value is 1, i.e., the number of ablated cells in the transient heat transfer model of C / C composite materials; EVOL represents the cell volume; A iniThe surface area where the heat flux boundary condition acts is the initial ablation surface area.

[0111] After calculation, the ablation depth-time curve and microscopic ablation morphology of short-cut C / C composites under oxygen-acetylene flame are shown in Figure 6 and Figure 7 When the test time is equal to 40s, the relative error between the predicted value of ablation depth (0.2775mm) and the experimental value (0.2848mm) is 2.63%. The above conclusions prove the effectiveness of the model.

[0112] The above description is only a specific implementation mode of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and these modifications or substitutions should be included in the protection scope of the present invention.

Claims

1. A method for predicting the ablation morphology of C / C composite materials considering complex microstructures, characterized in that: The following steps are involved: Step 1: Establish a transient heat transfer model of C / C composites based on microscale Step 1.1: Build the geometry and mesh it A cube model with the size of L×L×L was constructed for the C / C composite material; The sweep algorithm is used to divide the L×L×L cube model into hexahedral meshes to form NUM tol Units; Read the number of each unit and its center point coordinates; Step 1.2: Determine the component type of each unit A random sequence adsorption algorithm is used to generate the endpoint coordinates of the carbon fiber axis and the endpoint coordinates of the pore axis within the cubic model space. Combined with the number of each unit obtained in step 1.1 and the coordinates of its center point, the minimum distance between the center point of each unit and all carbon fiber axes is calculated in turn. When the minimum distance is less than or equal to the carbon fiber radius, the unit belongs to the carbon fiber component. Otherwise, the minimum distance between the unit and all pore axes is calculated. When the minimum distance is less than or equal to the pore radius, the unit belongs to the pore component; Otherwise, the unit belongs to the matrix component; Step 1.3: According to the component types of each unit determined in step 1.2, a carbon fiber unit set, a pore unit set and a matrix unit set are established, and corresponding material properties are assigned to each unit set. A heat flux boundary condition is applied to the initial ablation surface. Based on this, a transient heat transfer model of C / C composite materials based on microscale is established; Step 2: Perform transient heat transfer analysis and determine whether ablation occurs in each unit of the C / C composite cube model Initialize the current increment step number K and the unit number N to 1, summarize the numbers of the units corresponding to the initial ablation surface, and form a surface unit number set under the current increment step number; wherein the surface unit refers to the unit to which the surface directly contacts the external oxidizing gas belongs; Step 2.1: Execute the Kth incremental step of the transient heat transfer model of C / C composite materials based on microscale obtained in step 1, and read the Nth unit; according to the component type of the current unit, determine whether the current incremental step number K is KINC i If it is an integer multiple, then go to step 2.2; otherwise, determine whether N is equal to the total number of units NUM tol ; When N is less than NUM tol When N is equal to NUM tol When K is equal to the total number of incremental steps K tol ; When K is equal to K tol When , go to step 3; otherwise, define K = K + 1 and N = 1, update the surface unit number set and re-execute step 2.1; Among them, KINC i The number of incremental steps required to completely burn through a cell; Step 2.2: Determine whether the current unit N is a surface unit according to the surface unit number set under the current incremental step number; If it is a surface unit, the value of the current unit state variable 2 is defined as 1, otherwise, the value of the current unit state variable 2 is defined as 0; Step 2.3: Read the temperature value T of the current unit E , judge T E Is it greater than or equal to the critical ablation temperature T w,i ; If T E Greater than or equal to T w,i , indicating that ablation occurs, the value of the current unit state variable 1 is defined as 1; otherwise, indicating that ablation does not occur, the value of the current unit state variable 1 is defined as 0; Step 2.4: Determine whether the value of the current unit state variable 1 and the value of the state variable 2 are both 1; If the value of the current unit state variable 1 and the value of the state variable 2 are both 1, the number of the current unit is recorded in the set of ablated unit numbers corresponding to the current incremental step number, and step 3 is entered; Otherwise, go to step 2.5; Step 2.5: Determine whether N is equal to NUM tol ; When N is less than NUM tol When N is equal to NUM tol When K is equal to K tol , when K is equal to K tol , go to step 3; otherwise, define K=K+1 and N=1, update the surface unit number set, and re-execute step 2.1; Step 3: Predict the microscopic ablation morphology of C / C composites The post-processing function of the finite element analysis software is used to hide the ablated units corresponding to each incremental step obtained in step 2 to obtain the microscopic ablation morphology of the C / C composite material at different times.

2. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: Also includes: Step 4: Calculate the ablation depth The ablation depth s of C / C composite materials at different times was calculated according to the following formula t : Where NUM ab represents the number of ablated cells in the transient heat transfer model of C / C composite materials, EVOL represents the cell volume, and A ini represents the initial ablated surface area.

3. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: In step 1.2, the distance between the carbon fiber axis / pore axis and the center point of each unit is calculated using the following formula: In the formula, x p ,y p and z p is the coordinate of the end point P of the carbon fiber axis or pore axis; q ,y q and z q is the coordinate of the end point Q of the carbon fiber axis or pore axis; c ,y c and z c is the coordinate of the center point C of each unit; the line segment CN intersects the line segment PQ perpendicularly at point N(x n ,y n ,z n ).

4. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: In step 1.3, according to the component type of each unit determined in step 1.2, the corresponding density, thermal conductivity and specific heat capacity are assigned to different units, the analysis step type is defined as transient heat transfer analysis, and a transient heat transfer model of C / C composite materials based on microscale is established.

5. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: In step 2.1, the number of incremental steps required to completely burn through a cell is calculated using the following formula: In the formula, H ele Indicates the unit thickness; inc is the incremental step size of heat transfer analysis; is the linear ablation rate of carbon fiber and matrix of C / C composite material, where ρ represents density, C p is the specific heat capacity at constant pressure, represents the relative molecular mass of oxygen, h c represents the convective heat transfer coefficient, M C represents the relative molecular mass of carbon, represents the mass fraction of oxygen, η represents the stoichiometric coefficient of the oxidation reaction equation, and the subscript i=f, m represents the component type of the C / C composite material, f represents the carbon fiber, and m represents the matrix.

6. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: In step 2.4, the critical ablation temperature of carbon fiber and matrix of C / C composite material is calculated using the following formula: In the formula, ρ represents density, M O represents the relative molecular mass of oxygen atoms, A represents the pre-exponential factor, E represents the activation energy of the oxidation reaction, R represents the ideal gas constant, P W Indicates pressure.

7. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: In step 2, the specific steps of updating the surface unit number set are as follows: When K is greater than 1, the set of ablated unit numbers corresponding to the previous incremental step is read; The surface element number for the current increment is calculated using the following formula: In the formula, represents the surface unit number of step K, Indicates the number of ablated cells in step K-1, NUM lay Indicates the total number of units in each layer.

8. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to any one of claims 1 to 6, characterized in that: Step 1.1, step 1.3, step 2 and step 3 all use Abaqus software; among them, step 1.3 uses the Usdfld user subroutine in Abaqus software; Steps 1.2, 4 and updating the surface unit number are all performed using Matlab software.

9. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: In step 1.1, the value range of L is 0.1 mm to 1 mm.

10. The method for predicting ablation morphology of C / C composite materials considering complex microstructure according to claim 1, characterized in that: In step 2.1, the value of inc is 0.01s to 1s.

Images