Thermal-mechanical coupling simulation method and system for fiber reinforced polymer-based flexible thermal protection materials

Abstract

The application discloses a kind of force-thermal coupling simulation method and system for fiber reinforced polymer-based flexible thermal protection material, and relates to material multi-physical field numerical simulation technical field.The method comprises the following steps: constructing an ablation sample model and calibrating parameters;Full coupling model is constructed, including heat conduction, pyrolysis dynamics, large deformation mechanics equilibrium and force-thermal coupling constitutive, the two-way dynamic correlation between material stiffness and thermal expansion coefficient and pyrolysis conversion rate and temperature is established;Initial and loaded boundary conditions are set;Multi-physical field iterative solution is carried out by full coupling solver according to time step, and temperature, conversion rate and stress displacement and other physical fields are updated;Data is extracted to calculate ablation and deformation indexes for evaluation.The system includes the module for executing the above method.The application realizes the two-way full closed loop coupling of force, heat and chemical multi-field, accurately characterizes the non-monotonic mechanical evolution behavior of the material, and provides reliable support for the formulation optimization and structure design of thermal protection materials.

Description

Technical Field

[0001] This invention relates to the field of multiphysics numerical simulation technology for materials, and more specifically, to a mechanical-thermal-chemical coupling simulation method and system for fiber-reinforced polymer-based flexible thermal protection materials. Background Technology

[0002] Fiber-reinforced polymer-based flexible thermal protection materials are core materials for thermal protection systems in aerospace and other fields. Their service process in high-temperature and high-heat-flow environments involves complex multi-physics coupling processes such as heat conduction, polymer pyrolysis and carbonization chemical reactions, gas diffusion in porous media, evolution of material mechanical properties, and large deformation stress-strain response.

[0003] To assess and predict the service condition of such materials, numerical simulation is currently the commonly used method in the industry. However, traditional simulation methods often focus on single thermal response analysis or simplified thermo-chemical coupling analysis, which has the following technical limitations: First, existing material constitutive models often neglect the non-monotonic mechanical behavior of the polymer matrix during pyrolysis, characterized by "softening-carbonization stiffness recovery." Furthermore, existing methods do not adequately consider the inverse coupling effect of stress and strain states on heat transfer boundaries and macroscopic material properties under large structural deformations, making it difficult to accurately reflect the complex failure modes such as deformation, delamination, and cracking during ablation.

[0004] Second, in traditional simulation frameworks, data on pyrolysis reactions, temperature fields, and mechanical fields are mostly transmitted in a unidirectional manner, failing to achieve a two-way, fully closed-loop coupling of "temperature-pyrolysis conversion rate-material mechanical properties-deformation / stress-heat transfer / boundary conditions". Due to the lack of real-time mutual correction between physical fields, simulation results under long-term or extreme thermal loads deviate significantly from actual service behavior.

[0005] Third, for fiber-reinforced polymer-based materials, there is currently a lack of a comprehensive meso-macro equivalent modeling framework. Existing macro equivalent methods are insufficient to accurately characterize the composite effect between fibers and the matrix, thus failing to directly support the practical engineering needs for formulation optimization, structural design, and thermal protection performance evaluation of such flexible thermal protection materials. Summary of the Invention

[0006] The present invention aims to solve at least one of the aforementioned technical problems existing in the prior art.

[0007] Therefore, the first aspect of the present invention provides a mechanical-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials.

[0008] The second aspect of the present invention provides a mechanical-thermal-chemical coupling simulation system for fiber-reinforced polymer-based flexible thermal protection materials.

[0009] This invention provides a mechanical-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials, including: S1. Construct a macroscopic ablation specimen geometric model of fiber-reinforced polymer-based flexible thermal protection material and mesh it, and calibrate and obtain the material's basic stiffness tensor, basic thermal expansion coefficient and pyrolysis kinetic parameters based on the material's basic experimental data. S2. Construct a fully coupled mechanical-thermal-chemical model, which includes the heat conduction equation, the pyrolysis reaction kinetic equation, the large deformation quasi-static mechanical equilibrium equation, the thermo-elastic constitutive equation, and the coupled mechanical-thermal-chemical constitutive model. The coupled mechanical-thermal-chemical constitutive model establishes a correlation between the material's stiffness tensor and the pyrolysis conversion rate and temperature to characterize the non-monotonic evolution of the material's stiffness at different pyrolysis stages, and simultaneously establishes a correlation between the material's thermal expansion coefficient and the pyrolysis conversion rate and temperature. S3. Initialize the temperature field, pyrolysis conversion rate field, stress field and strain field, and apply heat flow load boundary conditions and displacement constraint boundary conditions corresponding to the actual ablation service scenario to the macroscopic ablation specimen geometric model. S4. Perform multiphysics coupling iterative solution according to a preset time step. During the iteration process of the current time step, the following steps are executed sequentially: solve the heat conduction equation based on the stress state and material properties of the previous time step to update the current temperature field; solve the pyrolysis reaction kinetic equation based on the updated current temperature field to update the current pyrolysis conversion rate field; use the mechanized thermo-coupled constitutive model to update the current stiffness tensor and thermal expansion coefficient based on the current temperature field and pyrolysis conversion rate field; solve the large deformation quasi-static mechanical equilibrium equation and thermoelastic constitutive equation based on the updated stiffness tensor and thermal expansion coefficient to obtain the current stress field, strain field and displacement field; use the obtained current stress field to update the thermo-mechanical coupling term in the heat conduction equation and proceed to the next time step iteration until the preset simulation end time is reached. S5. Extract simulation data at each time step, calculate the linear ablation rate, mass ablation rate and deformation engineering indicators, and evaluate the performance of the fiber-reinforced polymer-based flexible thermal protection material.

[0010] The mechano-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials according to the above-described technical solution of the present invention may also have the following additional technical features: In the above technical solution, the specific form of the heat conduction equation in step S2 is as follows:

[0011] in, The density of the material; Specific heat capacity at constant pressure; It is the product of the material density and the specific heat capacity at constant pressure, which changes dynamically with the pyrolysis conversion rate; For temperature; For time; For gradient operators; It is a divergence operator; It is the equivalent thermal conductivity, and satisfies ; The thermal conductivity of the original material; The thermal conductivity of the carbonized layer; This refers to the pyrolysis conversion rate, i.e., the carbonization rate. This is the internal heat source term, which includes the enthalpy change of the pyrolysis reaction; For reference temperature; It is the stress vector; This is a vector of thermal expansion coefficients.

[0012] In the above technical solution, in step S2, the pyrolysis reaction kinetic equation includes an equation describing the evolution of material density, and its specific form is as follows:

[0013] in, The density of the material; For time; This represents the total number of parallel pyrolysis reactions. These represent different reaction stages, components, or reaction pathways. For the first The pre-exponential factor of a reaction; For the first The activation energy of each reaction; It is the gas constant; For the first The temperature at which the reaction occurs; For the first The density of substances in a reaction; For the first The current density of the relevant substances in the reaction; For the first Density of residual carbon in each reaction; It is a natural constant; denoted as the reaction order.

[0014] In the above technical solution, the fully coupled mechanical-thermal model also includes a porous medium mass transport equation to describe the diffusion and pressure evolution of pyrolysis gas in the porous layer, the specific form of which is:

[0015] in, Indicates porosity; It is the ideal gas constant; For temperature; The pressure of the pyrolysis gas within the pores of the porous medium; For time; For gradient operators; Let be the Darcy velocity of the gas; Permeability of porous media; The instantaneous density of the solid skeleton; This refers to the carbonization layer density; The density of the original material; Porosity of the original material; The porosity is the value of the carbonized layer.

[0016] In the above technical solution, the specific form of the large deformation quasi-static mechanical equilibrium equation in step S2 is as follows:

[0017] in, It is the stress vector; The divergence of stress; Represents the volume force vector; The specific form of the thermoelastic constitutive equation is as follows:

[0018] in, It is the stress vector; For stiffness tensor; For total strain; For thermal expansion strain term, where This is a vector of thermal expansion coefficients. This represents the change in temperature.

[0019] In the above technical solution, in step S2, the specific form of the calculation equation for the stiffness tensor in the mechanically coupled constitutive model is as follows:

[0020]

[0021] in, Let be the stiffness tensor under the current state; This represents the basic stiffness tensor of the original fiber-reinforced polymer matrix material at room temperature. This is a piecewise function for the stiffness evolution during the pyrolysis stage; This is a temperature-thermal coupling stiffness correction term; This refers to the pyrolysis conversion rate, i.e., the carbonization rate. For temperature; The critical carbonization rate for pyrolysis initiation; The critical carbonization rate at the end of rapid pyrolysis softening; The critical carbonization rate at which the carbonization stiffness recovers is reached; This represents the stiffness correlation coefficient for the first stage. This is the stiffness correlation coefficient for the second stage.

[0022] In the above technical solution, the specific calculation equation for the temperature-thermal coupling stiffness correction term is as follows:

[0023]

[0024] in, The calculated value of the temperature-thermal coupling stiffness correction term; exp is the natural exponential function; The attenuation coefficient; For reference temperature; The basic thermal degradation coefficient; This is the temperature-carbonization coupling coefficient.

[0025] In the above technical solution, in step S2, the specific form of the equation for calculating the thermal expansion coefficient in the mechano-thermal coupling constitutive model is as follows:

[0026] in, The current coefficient of thermal expansion; The coefficient of thermal expansion of the original material at room temperature; This is a temperature correction term; To support the coefficients of piecewise functions; This refers to the change in temperature. This is a correction term for pyrolysis conversion rate; This is a correction factor; This refers to the pyrolysis conversion rate, also known as the carbonization rate.

[0027] In the above technical solution, step S4 involves using a globally discrete equation set based on the fully coupled solver of MOOSE (Multiphysics Object-Oriented Simulation Environment) for unified iterative solution. The specific matrix form of the globally discrete equation set is as follows:

[0028] in, This represents the coefficient matrix of the master equation for the temperature field. This represents the coefficient matrix of the master equation for the pyrolysis field; This represents the coefficient matrix of the master equation for the pore pressure field; Here is the coefficient matrix of the master equation for the displacement field; Represents the thermal coupling term matrix; Represents the thermo-mechanical coupling term matrix; This represents the Newton iteration increment of the temperature field; Represents the Newton iteration increment of the pyrolysis field; This represents the Newton iteration increment of the pore pressure field; This represents the Newton iteration increment of the displacement field; Represents the residual vector of the temperature field; Represents the residual vector of the thermal decomposition field; Represents the residual vector of the pore pressure field; This represents the residual vector of the displacement field.

[0029] This invention provides a mechano-thermo-chemical coupling simulation system for fiber-reinforced polymer-based flexible thermal protection materials, used to execute the mechano-thermo-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials as described in any of the above technical solutions. The system includes: The data calibration and modeling preprocessing module is used to construct a macroscopic ablation specimen geometric model of fiber-reinforced polymer-based flexible thermal protection material and divide it into meshes, as well as to obtain the basic stiffness tensor, basic thermal expansion coefficient and pyrolysis kinetic parameters of the material based on the material's basic experimental data. The fully coupled mathematical model construction module is used to construct a fully coupled mechanical-thermal model, including the heat conduction equation, the pyrolysis reaction kinetic equation, the large deformation quasi-static mechanical equilibrium equation, the thermo-elastic constitutive equation, and the mechanical-thermal coupled constitutive model. The boundary conditions and initialization module is used to initialize the temperature field, pyrolysis conversion rate field, stress field and strain field, and to apply heat flow load boundary conditions and displacement constraint boundary conditions corresponding to the actual ablation service scenario to the geometric model of the macroscopic ablation specimen. The multiphysics fully coupled solution module is used to perform iterative solutions by time step using a fully coupled solver. During the iteration of the current time step, the temperature field, pyrolysis conversion rate field, current stiffness tensor and thermal expansion coefficient, current stress field, strain field and displacement field are updated sequentially until the simulation ends. The post-processing and evaluation module is used to extract simulation data at each time step and calculate engineering indicators to complete the performance and failure risk assessment of the material.

[0030] In summary, due to the adoption of the above-mentioned technical features, the beneficial effects of the present invention are: This invention is the first to construct a mechano-thermal coupled constitutive model capable of accurately characterizing the non-monotonic mechanical behavior of materials during the "softening-carbonization stiffness recovery" process. Traditional simulation methods often neglect or simplify the complex mechanical evolution of polymer matrices at high temperatures. However, this invention introduces a piecewise function for stiffness evolution during the pyrolysis stage and a temperature-pyrolysis coupled stiffness correction term, establishing a direct two-way mathematical relationship between the material's stiffness tensor and the pyrolysis conversion rate and temperature field. This accurately reproduces the four stages of the material from pyrolysis initiation, rapid pyrolysis softening, carbonization stiffness recovery, to complete carbonization stability, thus realistically reflecting the complex failure modes of materials under extreme thermal loads, such as deformation, delamination, and cracking.

[0031] This invention achieves true bidirectional, fully closed-loop coupling of force, thermal, and chemical fields, enabling accurate solutions for large deformations. Traditional models often involve unidirectional data transfer between physical fields. This invention, based on a global discrete equation set and time-step iteration, not only drives the degradation and reconstruction of material mechanical properties through the temperature and pyrolysis fields in each calculation, but also applies the stress and strain fields induced by large deformations inversely to the thermal conductivity boundaries and thermo-mechanical coupling source terms. This multi-field closed-loop mechanism breaks down the barriers of single-response analysis, significantly improving the accuracy of numerical calculations under extreme thermal environments. Verification shows that the temperature and displacement errors predicted in its simulations can be strictly controlled within 3%.

[0032] The simulation method of this invention possesses excellent material versatility and engineering application value. This core coupled constitutive model, through directional correction of the underlying stiffness tensor and thermal expansion coefficient tensor, can not only accurately simulate the ablation response of isotropic materials, but also precisely capture and solve the complex strain concentration and redistribution behavior of anisotropic materials (such as textured materials like carbon fiber composites) caused by directional coupling. Combined with an automated iterative solution architecture, this invention significantly reduces the high cost of physical testing, providing reliable digital support for the formulation optimization, structural design, and quantitative assessment of failure risks of flexible thermal protection materials.

[0033] Additional aspects and advantages of the invention will become apparent in the following description or may be learned by practice of the invention. Attached Figure Description

[0034] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which: Figure 1 The flowchart illustrates the calculation steps of a mechanical-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials provided in this embodiment of the invention. Figure 2 The embodiment of the present invention provides a 2MW / m 2 Comparison of back temperature and surface temperature between ablation experimental data and numerical model calculations at heat flux density; Figure 3 The embodiment of the present invention provides a 2MW / m 2 A comparison of carbonization rate (conversion rate) between ablation experimental data and numerical model calculations at heat flux density; Figure 4 This is a comparison chart of numerically calculated displacement results and ablation experimental data provided in the embodiments of the present invention; Figure 5 Comparison cloud maps of strain distribution of materials in various directions provided for embodiments of the present invention. Detailed Implementation

[0035] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0036] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0037] The following reference Figures 1 to 5 This invention describes a mechanical-thermal-chemical coupling simulation method and system for fiber-reinforced polymer-based flexible thermal protection materials, provided by some embodiments of the present invention.

[0038] Some embodiments of this application provide a mechanical-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials.

[0039] To address the shortcomings of existing simulation methods in predicting the service behavior of fiber-reinforced polymer-based flexible thermal protection materials, such as insufficient characterization of non-monotonic mechanical evolution and unidirectional decoupling of multi-physics fields, this invention constructs a fully closed-loop coupled mechanical-thermal-mechanical numerical simulation framework. By introducing a piecewise function for stiffness evolution during the pyrolysis stage and a temperature-pyrolysis coupling correction mechanism, a mechanical-thermal-coupled constitutive model is established that directly correlates the material stiffness tensor and thermal expansion coefficient with the temperature field and pyrolysis conversion rate. This objectively characterizes the complex dynamic mechanical behavior of the polymer matrix at high temperatures, from "pyrolysis initiation - rapid pyrolysis softening - carbonization stiffness recovery - complete carbonization stabilization." Simultaneously, based on the fully coupled equation set and time-step iterative solution mechanism, real-time data interaction and reverse correction between heat conduction, pyrolysis reaction, porous media parameters, and large deformation mechanical response are achieved. This allows for accurate prediction of the material's multi-field response and failure risk in complex thermal load environments on a macroscopic scale.

[0040] Specifically, in combination Figure 1As shown, the present invention provides a mechanical-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials, which mainly includes the following steps S1-S5: S1. Construct a macroscopic ablation specimen geometric model of fiber-reinforced polymer-based flexible thermal protection material and mesh it, and calibrate and obtain the material's basic stiffness tensor, basic thermal expansion coefficient and pyrolysis kinetic parameters based on the material's basic experimental data. S2. Construct a fully coupled mechanical-thermal-chemical model, which includes the heat conduction equation, the pyrolysis reaction kinetic equation, the large deformation quasi-static mechanical equilibrium equation, the thermo-elastic constitutive equation, and the coupled mechanical-thermal-chemical constitutive model. The coupled mechanical-thermal-chemical constitutive model establishes a correlation between the material's stiffness tensor and the pyrolysis conversion rate and temperature to characterize the non-monotonic evolution of the material's stiffness at different pyrolysis stages, and simultaneously establishes a correlation between the material's thermal expansion coefficient and the pyrolysis conversion rate and temperature. S3. Initialize the temperature field, pyrolysis conversion rate field, stress field and strain field, and apply heat flow load boundary conditions and displacement constraint boundary conditions corresponding to the actual ablation service scenario to the macroscopic ablation specimen geometric model. S4. Perform multiphysics coupling iterative solution according to a preset time step. During the iteration process of the current time step, the following steps are executed sequentially: solve the heat conduction equation based on the stress state and material properties of the previous time step to update the current temperature field; solve the pyrolysis reaction kinetic equation based on the updated current temperature field to update the current pyrolysis conversion rate field; use the mechanized thermo-coupled constitutive model to update the current stiffness tensor and thermal expansion coefficient based on the current temperature field and pyrolysis conversion rate field; solve the large deformation quasi-static mechanical equilibrium equation and thermoelastic constitutive equation based on the updated stiffness tensor and thermal expansion coefficient to obtain the current stress field, strain field and displacement field; use the obtained current stress field to update the thermo-mechanical coupling term in the heat conduction equation and proceed to the next time step iteration until the preset simulation end time is reached. S5. Extract simulation data at each time step, calculate the linear ablation rate, mass ablation rate and deformation engineering indicators, and evaluate the performance of the fiber-reinforced polymer-based flexible thermal protection material.

[0041] Next, we will combine Figure 1 The implementation details of each step of the method provided by the present invention will be described in detail.

[0042] like Figure 1 As shown, step S1 is first performed to construct a macroscopic ablation specimen geometric model of fiber-reinforced polymer-based flexible thermal protection material and divide it into meshes, and to calibrate and obtain the basic stiffness tensor, basic thermal expansion coefficient and pyrolysis kinetic parameters of the material based on the material's basic experimental data.

[0043] In some embodiments, to accurately obtain the constitutive parameters of materials under different temperatures and physical fields, a series of basic material experiments are required to complete parameter calibration. Specifically, thermogravimetric (TG) analysis experiments are conducted to test the mass change of the material at different heating rates, thereby calibrating and extracting the pyrolysis kinetic parameters of the material (e.g., pre-exponential factor, activation energy, etc.); dynamic mechanical analysis (DMA) and high-temperature mechanical tests are conducted to test the stress-strain response of the material at different temperature nodes and degrees of pyrolysis, in order to calibrate the basic stiffness tensor at room temperature and the mechanical property evolution parameters required in subsequent constitutive models (e.g., stiffness evolution rate constant, etc.); simultaneously, thermal expansion tests are conducted to calibrate the basic thermal expansion coefficient of the material and its correction parameters as a function of temperature and pyrolysis conversion rate. These calibration data provide reliable physical input for subsequent fully coupled mechanical-thermal-chemical models.

[0044] In one specific embodiment, considering that fiber-reinforced polymer-based flexible thermal protection materials (such as carbon fiber-reinforced rubber matrix composites) are typically composed of reinforcing fibers and a polymer matrix, and have complex microstructures, directly assigning parameters at the macroscopic scale may lead to deviations. Therefore, before constructing the macroscopic geometric model, a microscopic equivalent model of the fiber-reinforced polymer matrix material can be established. Since the specific equivalent algorithm is not the core limitation of this invention, those skilled in the art can use well-known homogenization theories (such as the representative volume element RVE method or the mixing rule) to equivalently calculate the microscopic properties of the carbon fiber and rubber matrix as macroscopic model parameters, and then apply these equivalent macroscopic parameters to the macroscopic simulation calculations of this invention to further improve the accuracy of obtaining parameters for complex composite materials.

[0045] After obtaining the above parameters, a macroscopic ablation specimen geometric model (e.g., a cylindrical specimen model with a specific diameter and height) matching the actual ablation service scenario is constructed using a multiphysics object-oriented simulation environment (MOOSE) or other finite element analysis frameworks with equivalent computational capabilities. Subsequently, a finite element mesh is generated for this geometric model, and the aforementioned calibrated or equivalently obtained material properties (basic stiffness tensor, basic thermal expansion coefficient, etc.) are defined and assigned to the mesh elements, thereby completing the preprocessing preparation work for the simulation model.

[0046] After completing the model preprocessing and parameter calibration in step S1, step S2 is executed to construct a fully coupled mechanical-thermal-chemical model. This model aims to accurately describe the transient response and multi-field interaction mechanism of materials under high-temperature thermal loads through mathematical physics equations, specifically including the heat conduction equation, the pyrolysis reaction kinetic equation, the large deformation quasi-static mechanical equilibrium equation, the thermo-elastic constitutive equation, and the mechanical-thermal-chemical coupled constitutive model.

[0047] In one specific embodiment, for the thermal response part, the present invention adopts a heat conduction equation that includes an internal heat source and a mechanical-thermal coupling term, the specific form of which is:

[0048] in, The density of the material; Specific heat capacity at constant pressure; It is the product of material density and specific heat capacity at constant pressure, which changes dynamically with the increase of pyrolysis conversion rate, reflecting the change of material heat capacity property by the pyrolysis process; For temperature; For time; For gradient operators; It is a divergence operator; The equivalent thermal conductivity is used to accurately reflect the thermal conductivity characteristics at different pyrolysis stages. The equivalent thermal conductivity must satisfy the following conditions: ; The thermal conductivity of the original material; The thermal conductivity of the carbonized layer; , where is the pyrolysis conversion rate, i.e., the carbonization rate; Q is the internal heat source term, which includes the enthalpy change of the pyrolysis reaction; This is a thermo-mechanical coupling term used to describe the effect of thermoelastic deformation on the temperature field. For reference temperature; It is the stress vector; This is a vector of thermal expansion coefficients.

[0049] In some embodiments, a pyrolysis reaction kinetic model was constructed to characterize the chemical degradation process of the polymer matrix at high temperatures. This model describes the multi-order decay process of material density based on the Arrhenius equation, which is as follows:

[0050] This equation represents the total density decay caused by It is composed of multiple parallel pyrolysis reactions superimposed on each other. The density of the material; For time; This represents the total number of parallel pyrolysis reactions. These represent different reaction stages, components, or reaction pathways. For the first The pre-exponential factor of a reaction; For the first The activation energy of each reaction; It is the gas constant; For the first The temperature at which the reaction occurs; For the first The density of the substances (which may be volatile substances, etc.) in the reaction. For the first The current density of the relevant substances in the reaction, For the first The residual carbon density in each reaction, together with the terms they form, constitute a density-related reaction rate correction term; It is a natural constant; The reaction order reflects the sensitivity of the reaction rate to changes in density; different reaction orders reflect differences in the reaction mechanism.

[0051] Based on this, in order to quantitatively characterize the degree of pyrolysis of materials, this invention clarifies the definition of carbonization rate (i.e., pyrolysis conversion rate), and its calculation equation is as follows:

[0052] Among them, pyrolysis conversion rate It is a core indicator of the degree of material pyrolysis, representing the extent of the reaction, and usually indicates the proportion of material that has undergone pyrolysis; These represent the original mass of the material, its current mass, and its mass after complete carbonization, respectively.

[0053] Meanwhile, the pyrolysis rate (i.e., the change in carbonization rate over time) is described by the following equation:

[0054] This equation further refines the dynamic evolution of the pyrolysis process with temperature. Among other things, and These are the pyrolysis kinetic parameters (i.e., pre-exponential factor and activation energy) calibrated through the aforementioned basic experiments. Let be the ideal gas constant. For reaction-related parameters (such as heating rate). This is the reaction mechanism function.

[0055] Furthermore, considering that a porous medium will be formed during the ablation and pyrolysis of the material, the diffusion and pressure evolution of pyrolysis gases in the porous carbon layer will significantly affect the internal pressure and deformation of the material. Therefore, in a specific embodiment, a porous medium mass transport equation is also introduced into the fully coupled mechanical-thermal-chemical model:

[0056] in, Indicates porosity; It is the ideal gas constant; For temperature; The pressure of the pyrolysis gas within the pores of the porous medium; For time; For gradient operators; Let be the Darcy velocity of the gas; Permeability of porous media; The instantaneous density of the solid skeleton; This refers to the carbonization layer density; The density of the original material; Porosity of the original material; The porosity is the value of the carbonized layer.

[0057] In the equation, the first term on the right side of the equation represents the Darcy diffusion behavior of gas controlled by the permeability of the porous medium, the second term represents the change in gas density caused by temperature change, and the third term represents the source term for gas production from solid pyrolysis.

[0058] Ultimately, to obtain the coupling relationship between porous media, chemical reactions, and heat transfer, it is also necessary to use the energy transport equation to solve it:

[0059] in, Indicates the density of the pyrolysis gas; This represents the specific heat capacity of a gas at constant pressure. This represents the specific heat capacity of a solid at constant pressure. Indicates the Darcy velocity of the gas; Indicates the specific enthalpy of a gas; Indicates the specific enthalpy of a solid; This represents the enthalpy of the pyrolysis reaction. This equation is one of the core coupled equations of the entire model, taking into account the internal energy changes of the solid framework and porous gas, heat conduction, gas convection heat transfer, and the enthalpy change and internal heat source terms of the pyrolysis reaction.

[0060] For the mechanical response part, in order to simulate the large geometric deformation characteristics of materials during service, this invention adopts the large deformation quasi-static mechanical equilibrium equation and the thermoelastic constitutive equation.

[0061] The specific form of the large deformation quasi-static mechanical equilibrium equation is as follows:

[0062] in, It is the stress vector; The divergence of stress; Represents the volume force vector; The specific form of the thermoelastic constitutive equation is as follows:

[0063] in, It is the stress vector; For stiffness tensor; For total strain; For thermal expansion strain term, where This is a vector of thermal expansion coefficients. This represents the change in temperature.

[0064] As one of the core innovations of this invention, in order to organically combine the aforementioned relatively independent physical fields and overcome the defect in traditional models where the mechanical properties of materials only weaken monotonically with temperature, a mechanical-thermal coupled constitutive model is innovatively constructed in step S2. This constitutive model establishes a direct correlation between the material's stiffness tensor and the pyrolysis conversion rate and temperature, and the calculation formula is as follows:

[0065]

[0066] in, Let be the stiffness tensor under the current state; This is the basic stiffness tensor of the original fiber-reinforced polymer matrix material at room temperature. The relevant parameters can be obtained through standard experimental measurements. It is a piecewise function for stiffness evolution during the pyrolysis stage, used to accurately describe the four-stage mechanical behavior of the material: "pyrolysis not started - rapid pyrolysis softening - carbonization stiffness recovery - complete carbonization stabilization". This is a temperature-thermal coupling stiffness correction term; This refers to the pyrolysis conversion rate, i.e., the carbonization rate. For temperature; The critical carbonization rate for pyrolysis initiation; The critical carbonization rate at the end of rapid pyrolysis softening; The critical carbonization rate at which the carbonization stiffness recovers is reached; This is the stiffness correlation coefficient for the first stage (rapid pyrolysis stage); This is the stiffness correlation coefficient for the second stage (the stage of stiffness recovery due to carbonization).

[0067] Meanwhile, the temperature-thermal coupling stiffness correction term introduced in this coupled constitutive model describes the further weakening effect of the high-temperature thermal field on the material stiffness. Its specific form is:

[0068]

[0069] in, The calculated value of the temperature-thermal coupling stiffness correction term; exp is the natural exponential function; The attenuation coefficient; For reference temperature; The basic thermal degradation coefficient; Here, represents the temperature-carbonization coupling coefficient. By introducing a decay coefficient that depends on the pyrolysis conversion rate, the high-temperature softening effect in this model can be dynamically adjusted according to the degree of carbonization of the material itself.

[0070] In addition to the deep coupling of the stiffness tensor, this core constitutive model also achieves multi-field corrections to the thermal expansion behavior. The thermal expansion coefficient of the material is controlled by both temperature changes and carbonization reactions, and the calculation equation is as follows:

[0071] in, The current coefficient of thermal expansion; The coefficient of thermal expansion of the original material at room temperature; This is a temperature correction term; To support the coefficients of piecewise functions; This refers to the change in temperature. This is a correction term for pyrolysis conversion rate, reflecting the change in the thermal expansion behavior of the material before and after carbonization; This is a correction factor; This refers to the pyrolysis conversion rate, also known as the carbonization rate.

[0072] The fully coupled force-thermal-chemical model constructed through the above implementation method breaks the isolation between the force field, thermal field, and chemical field at the underlying governing equation framework. The evolution of any physical parameter (such as temperature or stress) can instantly trigger the linkage response of the material constitutive parameters, thus laying the mathematical foundation for the fully coupled closed-loop solution of large deformations in subsequent steps.

[0073] It should be noted that this thermo-mechanical coupling model is applicable to both isotropic and anisotropic materials.

[0074] After constructing the above-mentioned fully coupled force-thermal-chemical model, the next step is to execute step S3, initialize the temperature field, pyrolysis conversion rate field, stress field and strain field, and apply heat flow load boundary conditions and displacement constraint boundary conditions corresponding to the actual ablation service scenario to the geometric model of the macroscopic ablation sample.

[0075] In some embodiments, to perform multiphysics transient calculations, it is first necessary to assign initial state values ​​to each physical field within the simulation computational domain. Specifically, the initial temperature field of the model is set as the ambient reference temperature. Since the material has not yet undergone thermal degradation in its initial state, the initial value of the pyrolysis conversion rate (carbonization rate) field is set to 0 globally. At the same time, assuming that the material is in a stress-free and strain-free free state at the initial moment of calculation (or ignoring manufacturing residual stress), both the initial stress field and the initial strain field are set to 0.

[0076] In one specific embodiment, to accurately simulate the thermal state of materials under actual ablation tests (e.g., typical oxyacetylene ablation tests) or real service scenarios, it is necessary to apply engineered heat flux loads on specific boundaries of the macroscopic ablation specimen's geometric model. Specifically, the thermal boundary conditions can be set as follows: applying a dynamically varying heat flux load that changes with spatial location and time to the heated surface of the specimen. The input data for this heat flux load can be aerodynamic heating data obtained in advance through computational fluid dynamics (CFD) numerical calculations based on the oxyacetylene combustion reaction, or it can be a spatiotemporal distribution curve obtained by actual measurement and fitting using experimental instruments (such as calorimeters) to reflect the non-uniform and transient heating characteristics in the real ablation environment.

[0077] In some embodiments, in addition to thermal boundary conditions, displacement constraint boundary conditions also need to be applied to simulate the mechanical clamping and fixing state of the material during experimental testing or actual service. Taking the geometric model of a cylindrical macroscopic ablation specimen as an example, the mechanical boundary conditions can be specifically set as follows: applying axial displacement constraints on the bottom surface of the specimen, i.e. Apply radial displacement constraint to the side of the specimen, i.e. The combination of these engineered boundary conditions is crucial, as it accurately reflects the complex thermal stress distribution induced by external mechanical constraints and uneven internal temperature gradients during the thermal expansion and rapid carbonization of raw materials.

[0078] After initializing all the physical fields and applying boundary conditions, you can further set the global control parameters of the solver in the finite element simulation software (such as the aforementioned MOOSE framework), including the total simulation duration and the initial time step. The maximum number of iterations and the nonlinear solution convergence tolerance are established, thus creating a complete solution environment for subsequent fully coupled time-step iterative calculations.

[0079] After completing the physics initialization and engineering boundary condition settings, the core computational stage of this method will begin. Step S4 involves performing multiphysics coupling iterative solutions at preset time steps. During the iteration process of the current time step, each physics field is updated sequentially and the coupling terms are processed until the preset simulation end time is reached.

[0080] In some embodiments, the multiphysics coupling iterative process involved in step S4 is implemented using a fully coupled solver in the Multiphysics Object-Oriented Simulation Environment (MOOSE). Due to the strong nonlinear interactions between heat transfer, chemical reactions, and large deformation during the ablation process, this invention employs a time-step-based iterative calculation method.

[0081] In one specific embodiment, within each current time step, the system sequentially performs the following operations according to the logic of "thermal-chemical-mechanical" bidirectional feedback: First, the solver solves the heat conduction equation constructed in step S2 based on the stress state and material properties (such as the thermal conductivity of the previous time step) to update the temperature field of the current time step; second, based on the updated current temperature field, it solves the pyrolysis reaction kinetic equation to update the current pyrolysis conversion rate field distribution (i.e., carbonization rate distribution); simultaneously, it solves the porous medium mass transport equation based on the updated temperature field to update the pore gas pressure field within the porous medium; third, using the mechanical-chemical-thermal coupled constitutive model determined in step S2, it substitutes the latest temperature field and pyrolysis conversion rate field into the piecewise function and correction term to update the stiffness tensor and thermal expansion coefficient within the current material model mesh element in real time; then, based on the updated stiffness tensor and thermal expansion coefficient, it solves the large deformation quasi-static mechanical equilibrium equation and the thermoelastic constitutive equation to calculate the stress field, strain field, and nodal displacement field of the current time step. After obtaining all field variables, the thermo-mechanical coupling terms in the heat conduction equation (i.e., the change in internal energy caused by mechanical deformation) are updated in reverse using the latest stress field, as well as the squeezing effect of large deformation on the porosity of the porous medium. Then, the convergence of the overall residual is evaluated. If it does not converge, nonlinear iteration continues in the current time step. If it converges, the calculation of this step is completed and the next time step is entered until the simulation physical time ends.

[0082] To achieve unified discretization and fully closed-loop solution of the aforementioned physical fields at the underlying mathematical architecture, and to ensure computational stability and accuracy, in some embodiments, the system employs a globally discrete equation set based on a fully coupled solver for Newtonian unified iterative solution. Specifically, the large sparse matrices of heat transfer, pyrolysis, porous media fluids, and solid mechanics are assembled into the following globally discrete equation set:

[0083] in, This represents the coefficient matrix of the master equation for the temperature field (including heat capacity, mass, and thermal stiffness). This represents the coefficient matrix of the master equation of the pyrolysis field (corresponding to the lumped mass of the pyrolysis field). This represents the coefficient matrix of the master equation for the pore pressure field (including gas mass and permeability stiffness). The coefficient matrix of the displacement field master equation is composed of the material nonlinear stiffness matrix and the geometric large deformation nonlinear stiffness matrix, used to handle the mechanical response of large deformation; the off-diagonal terms of this matrix reflect a strong inter-field physical coupling mechanism. This represents the thermal-fluid coupling term matrix (reflecting the driving force of the temperature gradient on diffusion). This represents the thermo-mechanical coupling term matrix (reflecting the thermal effects generated by work done during large deformations). This represents the Newton iteration increment of the temperature field; Represents the Newton iteration increment of the pyrolysis field; This represents the Newton iteration increment of the pore pressure field; This represents the Newton iteration increment of the displacement field; Represents the residual vector of the temperature field; Represents the residual vector of the thermal decomposition field; Represents the residual vector of the pore pressure field; This represents the residual vector of the displacement field.

[0084] In the specific solution process, the solver discretizes the spatial domain using finite element discretization, employs an implicit time integration scheme to ensure numerical stability, and continuously corrects itself through iteration. This continues until the norm of each residual vector falls below the user-defined convergence tolerance. This fully coupled nonlinear equation system architecture completely breaks the limitation of "one-way data interpolation" from the temperature field to the mechanical field in traditional simulation software, ensuring a strict mathematical closed loop for the evolution of non-monotonic material properties and feedback from complex boundary conditions during ablation.

[0085] In a more specific embodiment, to ensure that the global discrete equations can be accurately executed by the computer program, a finite element spatial discretization scheme based on the Galerkin method and a backward Euler time integration scheme are used to construct the coefficient matrices and residual vectors of each physical field in the global discrete equations. Specifically, in the element-level discrete equations, the superscript e indicates the element level, and the superscripts n and n+1 represent the current and the time step to be solved, respectively. Furthermore, the meanings of the parameters involved in the following discrete equations are shown in Table 1: Table 1. Schematic diagram of finite element discretization parameters

[0086] The discrete equation systems are constructed as follows: (a) Discretization of the heat conduction equation

[0087] in, Represents the heat capacity mass matrix. This represents the temperature at time step n+1. This represents the temperature at the nth time step. Represents the thermal conductivity stiffness matrix. This represents the thermal load vector.

[0088] (b) Heat capacity mass matrix

[0089] in, Represents the elements of the heat capacity mass matrix; This indicates the computational domain under the reference configuration, with the superscript 'e' indicating a finite element. and Represents a shape function.

[0090] (c) Thermal conductivity stiffness matrix

[0091] in, Represents the elements of the thermal conductivity stiffness matrix; It is the equivalent thermal conductivity.

[0092] (d) Thermal load vector

[0093] in, This indicates the internal heat source of the enthalpy change in the pyrolysis reaction. This represents the thermo-coupling heat source term (the change in internal energy caused by mechanical deformation). Indicates the boundary thermal load (corresponding to the heat of reaction of oxyacetylene combustion). It represents the boundary of the computational domain under the reference configuration.

[0094] (e) Discrete pyrolysis reaction kinetic equations

[0095] in, Represents the lumped mass matrix of the pyrolysis. This represents the source term vector of the pyrolysis reaction.

[0096] (f) Pyrolysis lumped mass matrix

[0097] in, represents the elements of the pyrolysis mass matrix; diag(·) represents the diagonal operator.

[0098] (g) Pyrolysis reaction source term vector

[0099] in, This represents the rate function of the pyrolysis reaction.

[0100] (h) Discretization of gas transport equations in porous media

[0101] in, Represents the gas mass matrix. Represents the gas pressure matrix. Represents the permeability stiffness matrix. This represents the pressure load matrix.

[0102] (i) Gas mass matrix

[0103] in, Represents the elements of the gas mass matrix; Indicates porosity. This represents the pyrolysis conversion rate over n+1 time steps. Let be the ideal gas constant. This represents the average molar mass of the pyrolysis gas.

[0104] (j) Permeability stiffness matrix

[0105] in, Represents the elements of the permeability stiffness matrix. Indicates the thermal conductivity of a gas. This indicates the dynamic viscosity of the pyrolysis gas.

[0106] (k) Pressure load vector

[0107] in, This represents the pyrolysis gas production source term at time step n+1. This represents the boundary gas mass flux at time step n+1.

[0108] (l) Discretization of the equilibrium equations for large deformation mechanics

[0109] in, Represents the material stiffness matrix. Indicates the displacement increment. This represents the displacement field residual vector.

[0110] (m) Material tangent stiffness matrix (material nonlinearity + geometric nonlinearity)

[0111]

[0112]

[0113] in, Represents the nonlinear stiffness matrix of the material. Represents the geometric nonlinear stiffness matrix. Represents the elements of the material's nonlinear stiffness matrix. Represents the elements of the geometric nonlinear stiffness matrix. Represents the coupling stiffness matrix. Represents the strain-displacement matrix for large deformation. express The transpose of the matrix, Let represent the second type of Piola-Kirchhoff stress tensor at time step (n+1) and iteration (nl).

[0114] (n) Mechanical residual vector

[0115] in, Represents the surface force at time step n+1. Force per unit mass (volume force) This represents the initial density.

[0116] Through the specific assembly and construction of the aforementioned subfield matrix, it is possible to accurately realize the data transfer and nonlinear coupling solution of multiple fields of heat, force, and chemistry in macroscopic materials at each computational grid node.

[0117] After completing the time-series iterative solution of the multiphysics field full coupling, step S5 is executed to extract the simulation data at each time step and calculate the linear ablation rate, mass ablation rate and deformation engineering indicators to evaluate the performance of the fiber-reinforced polymer-based flexible thermal protection material.

[0118] In some embodiments, the post-processing module of the finite element simulation software is first used to extract the spatiotemporal distribution data of the temperature field, pyrolysis conversion rate field, stress field, strain field, and displacement field at each time step throughout the entire simulation calculation cycle. Based on the extracted data, key quantitative indicators for engineering evaluation can be further calculated. Specifically, these can be achieved using formulas. Calculate the line ablation rate, where, Linear ablation rate, The original thickness of the material. The current remaining thickness is given by t, where t is the elapsed time; similarly, the formula is used... Calculate the mass ablation rate, where, For mass ablation rate, For original quality, This represents the current remaining mass. Furthermore, mechanical response indicators such as maximum stress, maximum strain, and macroscopic deformation in various directions of the material can be directly extracted from the field data.

[0119] After obtaining the aforementioned engineering indicators, the structural reliability and thermal insulation performance of the material can be comprehensively evaluated. In one specific embodiment, the failure risk area inside the material can be accurately located by analyzing the spatiotemporal evolution of stress concentration. For example, due to the significant differences in thermodynamic properties between the carbonized layer, the pyrolysis layer, and the original material layer, stress concentration due to thermal expansion mismatch usually occurs at the pyrolysis interface. By combining anisotropic strain field data, failure modes such as delamination, cracking, or spalling that may occur in the material under extreme thermal loads can be quantitatively predicted, thereby providing digital support for the formulation optimization and structural design of thermal protection materials.

[0120] In some embodiments, a typical silicone rubber-based fiber-reinforced flexible thermal protection material (carbon fiber-reinforced rubber-based composite material) is used as the research object. A cylindrical sample model (30 mm in diameter and 10 mm in height) is established, simulating a surface heat source heat flux density of 2 MW / m². 2 The ablation process was simulated for 30 seconds, and the boundary conditions and model parameters were set as described above.

[0121] pass Figure 2 and Figure 3 Comparative analysis shows that this numerical model exhibits extremely high computational precision and accuracy, with the overall simulation error strictly controlled within 3%. For example... Figure 2 The numerical model accurately captured the rapid rise in surface temperature (reaching nearly 1900 K within 5 seconds and then leveling off) and the slow, nonlinear increase in back surface temperature over 30 seconds. The simulated curves closely matched the experimental data points, verifying the correctness of the heat transfer boundary and thermophysical parameters. Figure 3 Demonstrates the carbonization rate of the material The variation with temperature is shown. The typical S-shaped conversion rate curve obtained from the numerical simulation perfectly matches the experimental discrete points, accurately reflecting the nonlinear degradation evolution process of the material in the range of 400K to 1100K. In summary, the numerical model has been fully validated at the thermo-chemical coupling level and can reliably and quantitatively predict the complex behavior of materials under extreme thermal loads.

[0122] Furthermore, combined Figure 4 The numerical displacement results shown are compared with the ablation experimental data. By extracting the displacement field and processing the data, the structural evolution process of the carbonized layer, pyrolysis layer and original material layer formed sequentially from the surface to the interior of the material can be intuitively reproduced. The physical phenomenon of the pyrolysis interface advancing into the interior over time is clearly visible, and the macroscopic ablation displacement at the center point can be quantitatively compared (i.e., the displacement vector modulus disp_vector shown in the figure).

[0123] Specifically, Figure 4Simulation results show that the material forms a carbonized layer, a pyrolysis layer, and a pristine material layer sequentially from the surface to the interior, with the pyrolysis interface advancing inward over time. Quantitative data indicate that the measured displacement induced by ablation at the material's center at a power density of 2 MW / m³... 2 The value is 0.96 mm at 5 MW / m 2 The actual value was 2.74 mm; while the model's calculated prediction was 1.01 mm (2MW / m²). 2 ) and 2.39 mm (5MW / m 2 The corresponding quantization errors for the two models are 5.2% and 12.7%, respectively, both within the acceptable deviation range for engineering requirements. Besides the errors in experimental operations and numerical model calculations, the model errors may also be influenced by mechanical damage.

[0124] To more clearly demonstrate that the model can calculate and predict ablation behavior not only under isotropic conditions but also under anisotropic conditions, this paper establishes a rectangular material model and calculates its ablation behavior under both isotropic and anisotropic conditions. Figure 5As shown, the strain contour plot clearly demonstrates the constitutive model's ability to calculate the ablation mechanical behavior of both isotropic and anisotropic materials. The first three figures (a), (b), and (c) represent the isotropic calculation results, while the last three figures (d), (e), and (f) represent the anisotropic calculation results. Figures (a), (b), and (c) exhibit high geometric symmetry. Under isotropic conditions, the elastic modulus and Poisson's ratio of the material are identical in all directions. Therefore, the strain fields of strain_xx(a) (representing the normal strain in the x-direction) and strain_zz(c) (representing the normal strain in the z-direction) exhibit a nearly perfect, orthogonally mirror-image circular symmetrical distribution. The central region (red area) of strain_yy(b) (the normal strain in the y-direction) shows a regular circular high-strain concentration area, uniformly decreasing outwards, which is consistent with the typical response of isotropic materials under centrosymmetric loading. Symmetry breaking and waveform changes in Figures (d), (e), and (f): When anisotropy is introduced (Figures (d), (e), and (f)), the distribution characteristics of the strain field change significantly. The strain contours of strain_xx(d) and strain_zz(f) are no longer regular rings, but exhibit obvious direction selectivity. The stripes and wave characteristics reflect the differences in stiffness in different crystal orientations or material principal axes. Strain reversal and redistribution: The most significant change is reflected in strain_yy(e). Compared to the high strain at the center (red region) in isotropic conditions, the central strain in anisotropic conditions is reduced to the lowest (blue region), and the maximum strain shifts to the outer ring region. This strain redistribution phenomenon strongly depends on the anisotropic components of the elastic tensor of the anisotropic matrix. The comparison chart strongly demonstrates the model's versatility: when an isotropic matrix is ​​input, the model can accurately degenerate and ensure the disappearance of anisotropic features, outputting a perfectly symmetric field; when directional stiffness (such as the constitutive properties of carbon fiber composites or textured materials) is taken into account, the model can accurately capture and solve the complex strain concentration and wave behavior caused by directional coupling.

[0125] The above examples fully demonstrate that the method of the present invention can accurately solve complex mechanical behaviors caused by directional coupling, and has excellent engineering versatility and application value.

[0126] Based on the same inventive concept as the aforementioned mechano-thermal coupling simulation method, other embodiments of the present invention also provide a mechano-thermal coupling simulation system for fiber-reinforced polymer-based flexible thermal protection materials. This system can be implemented through software, hardware, or a combination of both, for example, by developing a dedicated mechano-thermal coupling calculation program deployed on a computer device to achieve fully automated solution processing. Specifically, the simulation system mainly includes the following functional modules: The data calibration and modeling preprocessing module is mainly used to perform the aforementioned step S1. Specifically, this module provides a model preprocessing interface for constructing a macroscopic ablation specimen geometric model of fiber-reinforced polymer-based flexible thermal protection material and generating a finite element mesh. Simultaneously, this module receives calibration data based on fundamental material experiments such as thermogravimetric analysis, dynamic mechanical analysis, high-temperature mechanical tests, and thermal expansion tests to obtain and store the material's fundamental stiffness tensor, fundamental coefficient of thermal expansion, and pyrolysis kinetic parameters, and assigns these attribute parameters to the corresponding mesh elements of the model.

[0127] The fully coupled mathematical model construction module is mainly used to execute the aforementioned step S2. This module is the core configuration unit of the entire simulation system, used to establish a complete set of physical control equations in the underlying architecture of the simulation environment, including heat conduction equations, pyrolysis reaction kinetic equations, porous medium mass transport equations, large deformation quasi-static mechanical equilibrium equations, and thermoelastic constitutive equations. In particular, this module is used to write the control logic of the mechanically coupled constitutive model into the system, so that the stiffness tensor and thermal expansion coefficient of the material can be strictly correlated with the pyrolysis conversion rate and the local temperature field according to the set piecewise function and dynamic correction term.

[0128] The boundary conditions and initialization module is mainly used to execute the aforementioned step S3. This module is responsible for initializing all state variables, such as the temperature field, pyrolysis conversion rate field, stress field, and strain field, within the system's computational domain before the transient time-series calculation starts. At the same time, this module supports the binding and application of engineering boundary conditions, allowing users to apply dynamically changing thermal flux load boundary conditions that vary with space and time on specific boundaries of the specimen's geometric model, as well as bottom and side displacement constraint boundary conditions used to simulate the actual clamping state of the specimen.

[0129] The multiphysics fully coupled solution module is mainly used to execute the aforementioned step S4. In a specific embodiment, this module relies on the fully coupled solver of the multiphysics underlying framework (such as the MOOSE framework) and uses the finite element spatial discretization and time integration algorithm to perform Newtonian nonlinear iterative calculations based on the global discrete equation system. Within each simulation time step, the module automatically coordinates the data transfer of each subphysics field, sequentially solves and updates the current temperature field, pyrolysis conversion rate field, material stiffness and thermal expansion coefficient, stress field, strain field, and displacement field, and corrects the force-thermal coupling source term in the heat transfer equation by extracting the current stress field for back feedback. The module continuously iterates within the current time step until the residual converges, and then advances to the next time step until the set physical ablation cycle is completed.

[0130] The post-processing and evaluation module is primarily used to execute the aforementioned step S5. After the fully coupled iterative solution is completed, this module extracts the spatiotemporal distribution data of field variables at each time step from the system and automatically calculates and outputs key engineering evaluation indicators such as line ablation rate, mass ablation rate, and macroscopic deformation. Furthermore, this module also features data comparison and visualization analysis capabilities, rendering strain distribution cloud maps in various directions. This assists engineers in accurately locating stress concentration phenomena in areas such as the interface between the carbonized layer and the pyrolysis layer, thereby completing a quantitative analysis of material failure risks and providing reliable digital support for the formulation optimization and structural design of flexible thermal protection materials.

[0131] In some embodiments, the mechano-thermo-chemical coupling simulation system for fiber-reinforced polymer-based flexible thermal protection materials provided by this invention can be specifically deployed or integrated into an electronic computing device (such as a computer terminal or server cluster) containing at least one processor and at least one memory. The memory is used to store computer program instructions for implementing the various functional modules described above, and the processor is used to read and execute the computer program instructions, thereby implementing the method flow described in this invention. Those skilled in the art should understand that the above division of system modules is only a logical functional division; in actual software code implementation, each module can be integrated into a main program or distributed in different subroutines or dynamic link libraries, all without departing from the protection scope of this invention.

[0132] In this specification, the illustrative expressions of the terms used do not necessarily refer to the same embodiments or examples. Moreover, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0133] Any modifications, equivalent substitutions, or improvements made within the spirit and principles of this invention shall be included within the scope of protection of this invention.

Claims

1. A mechano-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials, characterized in that, include: S1. Construct a macroscopic ablation specimen geometric model of fiber-reinforced polymer-based flexible thermal protection material and mesh it, and calibrate and obtain the material's basic stiffness tensor, basic thermal expansion coefficient and pyrolysis kinetic parameters based on the material's basic experimental data. S2. Construct a fully coupled mechanical-thermal-chemical model, which includes the heat conduction equation, the pyrolysis reaction kinetic equation, the large deformation quasi-static mechanical equilibrium equation, the thermo-elastic constitutive equation, and the coupled mechanical-thermal-chemical constitutive model. The coupled mechanical-thermal-chemical constitutive model establishes a correlation between the material's stiffness tensor and the pyrolysis conversion rate and temperature to characterize the non-monotonic evolution of the material's stiffness at different pyrolysis stages, and simultaneously establishes a correlation between the material's thermal expansion coefficient and the pyrolysis conversion rate and temperature. S3. Initialize the temperature field, pyrolysis conversion rate field, stress field and strain field, and apply heat flow load boundary conditions and displacement constraint boundary conditions corresponding to the actual ablation service scenario to the macroscopic ablation specimen geometric model. S4. Perform multiphysics coupling iterative solution according to a preset time step. During the iteration process of the current time step, the following steps are executed sequentially: solve the heat conduction equation based on the stress state and material properties of the previous time step to update the current temperature field; solve the pyrolysis reaction kinetic equation based on the updated current temperature field to update the current pyrolysis conversion rate field; use the mechanized thermo-coupled constitutive model to update the current stiffness tensor and thermal expansion coefficient based on the current temperature field and pyrolysis conversion rate field; solve the large deformation quasi-static mechanical equilibrium equation and thermoelastic constitutive equation based on the updated stiffness tensor and thermal expansion coefficient to obtain the current stress field, strain field and displacement field; use the obtained current stress field to update the thermo-mechanical coupling term in the heat conduction equation and proceed to the next time step iteration until the preset simulation end time is reached. S5. Extract simulation data at each time step, calculate the linear ablation rate, mass ablation rate and deformation engineering indicators to evaluate the performance of the fiber-reinforced polymer-based flexible thermal protection material. In step S2, the specific form of the equation for calculating the stiffness tensor in the mechanically coupled constitutive model is as follows: in, Let be the stiffness tensor under the current state; This represents the basic stiffness tensor of the original fiber-reinforced polymer matrix material at room temperature. This is a piecewise function for the stiffness evolution during the pyrolysis stage; This is a temperature-thermal coupling stiffness correction term; This refers to the pyrolysis conversion rate, i.e., the carbonization rate. For temperature; The critical carbonization rate for pyrolysis initiation; The critical carbonization rate at the end of rapid pyrolysis softening; The critical carbonization rate at which the carbonization stiffness recovers is reached; This represents the stiffness correlation coefficient for the first stage. This represents the stiffness correlation coefficient for the second stage. The specific calculation equation for the temperature-thermal decoupling stiffness correction term is as follows: in, The calculated value of the temperature-thermal coupling stiffness correction term; exp is the natural exponential function; The attenuation coefficient; For reference temperature; The basic thermal degradation coefficient; This is the temperature-carbonization coupling coefficient; In step S2, the specific form of the equation for calculating the thermal expansion coefficient in the mechano-thermal coupled constitutive model is as follows: in, The current coefficient of thermal expansion; The coefficient of thermal expansion of the original material at room temperature; This is a temperature correction term; To support the coefficients of piecewise functions; This refers to the change in temperature. This is a correction term for pyrolysis conversion rate; This is a correction factor; This refers to the pyrolysis conversion rate, also known as the carbonization rate.

2. The mechano-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials according to claim 1, characterized in that, In step S2, the specific form of the heat conduction equation is: in, The density of the material; Specific heat capacity at constant pressure; It is the product of the material density and the specific heat capacity at constant pressure, which changes dynamically with the pyrolysis conversion rate; For temperature; For time; For gradient operators; For divergence operators; It is the equivalent thermal conductivity, and satisfies ; The thermal conductivity of the original material; The thermal conductivity of the carbonized layer; This refers to the pyrolysis conversion rate, i.e., the carbonization rate. This is the internal heat source term, which includes the enthalpy change of the pyrolysis reaction; For reference temperature; It is the stress vector; This is a vector of thermal expansion coefficients.

3. The mechano-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials according to claim 1, characterized in that, In step S2, the pyrolysis reaction kinetic equation includes an equation describing the evolution of material density, which is in the following specific form: in, The density of the material; For time; This represents the total number of parallel pyrolysis reactions. These represent different reaction stages, components, or reaction pathways. For the first The pre-exponential factor of a reaction; For the first The activation energy of each reaction; It is the gas constant; For the first The temperature at which the reaction occurs; For the first The density of substances in a reaction; For the first The current density of the relevant substances in the reaction; For the first Density of residual carbon in each reaction; It is a natural constant; denoted as the reaction order.

4. The mechano-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials according to claim 1, characterized in that, The aforementioned fully coupled mechanical-thermal model also includes a porous medium mass transport equation to describe the diffusion and pressure evolution of pyrolysis gases in the porous layer. Its specific form is as follows: in, Indicates porosity; It is the ideal gas constant; For temperature; The pressure of the pyrolysis gas within the pores of the porous medium; For time; For gradient operators; The Darcy velocity of the gas; Permeability of porous media; The instantaneous density of the solid skeleton; This refers to the carbonization layer density; The density of the original material; Porosity of the original material; The porosity is the value of the carbonized layer.

5. The mechano-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials according to claim 1, characterized in that, In step S2, the specific form of the large deformation quasi-static mechanical equilibrium equation is as follows: in, It is the stress vector; The divergence of stress; Represents the volume force vector; The specific form of the thermoelastic constitutive equation is as follows: in, It is the stress vector; For stiffness tensor; For total strain; For thermal expansion strain term, where This is a vector of thermal expansion coefficients. This represents the change in temperature.

6. The mechano-thermal-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials according to claim 1, characterized in that, In step S4, a unified iterative solution is performed using a global discrete equation system based on the MOOSE fully coupled solver. The specific matrix form of the global discrete equation system is as follows: in, This represents the coefficient matrix of the master equation for the temperature field. This represents the coefficient matrix of the master equation for the pyrolysis field; This represents the coefficient matrix of the master equation for the pore pressure field; Here is the coefficient matrix of the master equation for the displacement field; Represents the thermal coupling term matrix; Represents the thermo-mechanical coupling term matrix; This represents the Newton iteration increment of the temperature field; Represents the Newton iteration increment of the pyrolysis field; This represents the Newton iteration increment of the pore pressure field; This represents the Newton iteration increment of the displacement field; Represents the residual vector of the temperature field; Represents the residual vector of the thermal decomposition field; Represents the residual vector of the pore pressure field; This represents the residual vector of the displacement field.

7. A mechano-thermo-chemical coupling simulation system for fiber-reinforced polymer-based flexible thermal protection materials, characterized in that, A system for performing a mechano-thermo-chemical coupling simulation method for fiber-reinforced polymer-based flexible thermal protection materials as described in any one of claims 1 to 6, the system comprising: The data calibration and modeling preprocessing module is used to construct a macroscopic ablation specimen geometric model of fiber-reinforced polymer-based flexible thermal protection material and mesh it, as well as to obtain the basic stiffness tensor, basic thermal expansion coefficient and pyrolysis kinetic parameters of the material based on the material's basic experimental data. The fully coupled mathematical model construction module is used to construct a fully coupled mechanical-thermal model, including the heat conduction equation, the pyrolysis reaction kinetic equation, the large deformation quasi-static mechanical equilibrium equation, the thermo-elastic constitutive equation, and the mechanical-thermal coupled constitutive model. The boundary conditions and initialization module is used to initialize the temperature field, pyrolysis conversion rate field, stress field and strain field, and to apply heat flow load boundary conditions and displacement constraint boundary conditions corresponding to the actual ablation service scenario to the geometric model of the macroscopic ablation specimen. The multiphysics fully coupled solution module is used to perform iterative solutions by time step using a fully coupled solver. During the iteration of the current time step, the temperature field, pyrolysis conversion rate field, current stiffness tensor and thermal expansion coefficient, current stress field, strain field and displacement field are updated sequentially until the simulation ends. The post-processing and evaluation module is used to extract simulation data at each time step and calculate engineering indicators to complete the performance and failure risk assessment of the material.

Images