Three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupling prediction method

Abstract

The application relates to the technical field of computer-aided design, and particularly discloses a three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupling prediction method, which comprises the following steps: constructing a porous medium geometric model of carbon / phenolic composite under three-dimensional pore scale; setting a calculation domain, generating a structured grid covering the fluid domain and the solid domain of the porous medium geometric model according to the grid node distribution set in the calculation domain; according to the spatial distribution information of the porous medium geometric model, performing attribute identification on each grid unit in the structured grid, marking the solid grid unit and the fluid grid unit, and assigning the solid phase material attribute to the solid grid unit; based on the marked grid units and the assigned solid phase material attribute, performing multi-physical field coupling iterative solving on the fluid domain and the solid domain of the carbon / phenolic composite, obtaining a coupling calculation result, and calculating the ablation recession distance and the ablation recession rate of the carbon / phenolic composite according to the coupling calculation result.

Description

Technical Field

[0001] This application relates to the field of computer-aided design technology, and in particular to a multi-field coupled prediction method for pyrolysis ablation of three-dimensional porous carbon / phenolic composites. Background Technology

[0002] Hypersonic vehicles, flying at high speeds within the atmosphere, endure extremely harsh aerodynamic heating environments, with surface temperatures reaching thousands of degrees Celsius, placing extremely high demands on thermal protection systems. Carbon / phenolic composite materials, as a typical passive ablation thermal protection material, have been widely used in high-speed vehicles such as ballistic missiles, reentry capsules / entry vehicles, and near-space vehicles due to their excellent high-temperature resistance, good ablation resistance, and lightweight properties. When the material is in a high-enthalpy incoming flow environment, complex physicochemical changes occur: the phenolic resin matrix undergoes pyrolysis upon heating, absorbing a large amount of heat and generating pyrolysis gases; the carbon fiber reinforcement undergoes oxidation with oxidizing components in the incoming flow, further consuming the solid phase material. These processes involve multi-physics coupling effects, and the material's pore structure continuously evolves during ablation, making the ablation mechanism extremely complex and directly related to the reliability assessment and refined design of thermal protection systems.

[0003] Existing simulation methods for the ablation behavior of carbon / phenolic composites have significant limitations. Macro-scale numerical simulations typically employ homogenized models, simplifying porous media into continuous media, which fails to analyze the flow details and the actual evolution of the solid-phase material at the pore scale. In terms of experimental research, while high-temperature ablation experiments can characterize the surface morphology and overall retreat of the material after ablation, they struggle to reveal the dynamic evolution mechanism of multi-field coupling during ablation and suffer from difficulties such as expensive equipment, complex implementation, and significant material loss. Some pore-scale simulation methods often use simplified two-dimensional models or only consider the consumption of a single solid-phase material, failing to distinguish between the two different reaction mechanisms of carbon fiber oxidation and phenolic resin pyrolysis. Summary of the Invention

[0004] In view of this, this application provides a three-dimensional porous-scale carbon / phenolic composite pyrolysis ablation multi-field coupled prediction method and device, storage medium, and computer equipment. By constructing a porous medium geometric model containing spatial distribution information of carbon fiber and phenolic resin matrix, and generating a structured mesh covering the fluid and solid domains based on the OpenFOAM platform, it overcomes the limitation of macroscopic homogenized models that cannot resolve pore-scale flow details and the real evolution process of the solid phase, and achieves accurate characterization of the material's microstructure. By identifying the attributes of the mesh elements and assigning corresponding solid phase material attributes to the carbon fiber and phenolic resin matrix respectively, it solves the technical blind spot in existing porous-scale simulation methods that only consider the consumption of a single solid phase material and fail to distinguish between the two different reaction mechanisms of carbon fiber oxidation and phenolic resin pyrolysis, thus enabling a true reflection of the differentiated ablation behavior of different components in the composite material. Based on the labeled mesh elements, the ThermalFOAM solver configured for multi-physics field coupled iterative solution is used, which effectively overcomes the shortcomings of traditional experimental methods that are difficult to reveal the dynamic evolution mechanism of flow-solid-thermal-chemical processes during ablation, thus facilitating the acquisition of high-precision three-dimensional dynamic prediction results.

[0005] According to one aspect of this application, a multi-field coupled prediction method for pyrolysis ablation of three-dimensional porous carbon / phenolic composites is provided, comprising: A porous media geometric model of carbon / phenolic composite material at the three-dimensional pore scale is constructed. The porous media geometric model includes the spatial distribution information of the solid domain composed of carbon fiber and phenolic resin matrix, and the spatial distribution information of the fluid domain composed of pores. In the OpenFOAM environment, a computational domain containing the porous medium geometric model is set, and a grid node distribution is set within the computational domain. Based on the grid node distribution, a structured grid covering the fluid and solid domains of the porous medium geometric model is generated. Based on the spatial distribution information of the porous medium geometric model, attribute identification is performed on each grid cell in the structured grid. Solid grid cells belonging to carbon fiber and phenolic resin matrix are marked, as well as fluid grid cells belonging to pores. The solid grid cells are then assigned corresponding solid material properties. Based on the marked mesh elements and the assigned solid material properties, the ThermalFOAM solver, configured to perform multiphysics coupled iterative solution, is used in the OpenFOAM environment to perform multiphysics coupled iterative solution on the fluid and solid domains of the carbon / phenolic composite. The coupled calculation results are obtained, and based on the coupled calculation results, the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite are calculated.

[0006] According to another aspect of this application, a three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupled prediction device is provided, comprising: The geometric model construction module is used to construct a porous medium geometric model of carbon / phenolic composite at the three-dimensional pore scale. The porous medium geometric model includes spatial distribution information of the solid domain composed of carbon fiber and phenolic resin matrix, and spatial distribution information of the fluid domain composed of pores. The mesh generation module is used to set the computational domain containing the porous medium geometric model in the OpenFOAM environment, set the mesh node distribution in the computational domain, and generate a structured mesh covering the fluid domain and solid domain of the porous medium geometric model according to the mesh node distribution. The grid marking module is used to identify the attributes of each grid cell in the structured grid according to the spatial distribution information of the porous medium geometric model, mark the solid grid cells belonging to the carbon fiber and phenolic resin matrix, and the fluid grid cells belonging to the pores, and assign the corresponding solid phase material properties to the solid grid cells. The coupled solution module is used to perform multiphysics coupled iterative solution on the fluid and solid domains of the carbon / phenolic composite material in the OpenFOAM environment based on the marked mesh elements and the assigned solid material properties, using the ThermalFOAM solver configured to perform multiphysics coupled iterative solution. The coupled calculation results are obtained, and based on the coupled calculation results, the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite material are calculated.

[0007] According to another aspect of this application, a storage medium is provided that stores a computer program thereon, which, when executed by a processor, implements the above-described three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupling prediction method.

[0008] According to another aspect of this application, a computer device is provided, including a storage medium, a processor, and a computer program stored on the storage medium and executable on the processor, wherein the processor executes the program to implement the above-described three-dimensional porous-scale carbon / phenolic composite pyrolysis ablation multi-field coupling prediction method.

[0009] Using the above technical solution, this application provides a three-dimensional porous-scale multi-field coupled prediction method and apparatus for pyrolysis ablation of carbon / phenolic composites, a storage medium, and a computer device. First, a porous medium geometric model of the carbon / phenolic composite material at the three-dimensional porous scale is constructed. Then, in the OpenFOAM environment, a computational domain encompassing this porous medium geometric model is defined, and a structured mesh covering the entire porous medium geometric model is generated. After mesh generation, each mesh cell can be attribute-identified and labeled to determine whether it is a solid mesh cell or a fluid mesh cell. For mesh cells labeled as solid, corresponding solid-phase material properties can be assigned. Based on the above labeling and assignment results, a specially developed ThermalFOAM solver is used for multi-physics coupled iterative solution, finally outputting the coupled calculation results. Based on the coupled calculation results, the ablation retreat distance and ablation retreat rate can be calculated. This application's embodiments construct a porous media geometric model containing spatial distribution information of carbon fiber and phenolic resin matrix, and generate a structured mesh covering both fluid and solid domains based on the OpenFOAM platform. This overcomes the limitations of macroscopic homogenized models in failing to analyze pore-scale flow details and the actual evolution of the solid phase, achieving accurate characterization of the material's microstructure. By identifying the properties of the mesh elements and assigning corresponding solid-phase material properties to the carbon fiber and phenolic resin matrix, the technical blind spots commonly found in existing pore-scale simulation methods—which only consider the consumption of a single solid-phase material and fail to distinguish between the two different reaction mechanisms of carbon fiber oxidation and phenolic resin pyrolysis—are addressed. This allows for a true reflection of the differentiated ablation behavior of different components in the composite material. Based on the labeled mesh elements, a ThermalFOAM solver configured for multi-physics coupled iterative solution is used, effectively overcoming the shortcomings of traditional experimental methods in revealing the dynamic evolution mechanism of flow-solid-thermal-chemical processes during ablation, thus facilitating the acquisition of high-precision three-dimensional dynamic prediction results.

[0010] The above description is only an overview of the technical solution of this application. In order to better understand the technical means of this application and to implement it in accordance with the contents of the specification, and to make the above and other objects, features and advantages of this application more obvious and understandable, the following are specific embodiments of this application. Attached Figure Description

[0011] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings: Figure 1 The diagram shows a flowchart of a three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupling prediction method provided in an embodiment of this application. Figure 2This paper presents a schematic diagram showing the change curve of ablation retreat rate over time according to an embodiment of this application. Figure 3 This diagram illustrates the solid-phase evolution of a carbon fiber and phenolic resin matrix according to an embodiment of this application. Figure 4 This illustration shows a schematic diagram of a three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupling prediction device provided in an embodiment of this application. Figure 5 A schematic diagram of the device structure of a computer device provided in an embodiment of this application is shown. Detailed Implementation

[0012] The present application will be described in detail below with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in the embodiments of the present application can be combined with each other.

[0013] This embodiment provides a three-dimensional pore-scale multi-field coupled prediction method for pyrolysis ablation of carbon / phenolic composites, such as... Figure 1 As shown, the method includes: Step 101: Construct a porous media geometric model of carbon / phenolic composite at the three-dimensional pore scale. The porous media geometric model includes spatial distribution information of the solid domain composed of carbon fiber and phenolic resin matrix, and spatial distribution information of the fluid domain composed of pores.

[0014] Step 102: In the OpenFOAM environment, set the computational domain containing the porous medium geometric model, and set the grid node distribution in the computational domain. Based on the grid node distribution, generate a structured grid covering the fluid domain and solid domain of the porous medium geometric model.

[0015] Step 103: Based on the spatial distribution information of the porous medium geometric model, identify the attributes of each grid cell in the structured grid, mark the solid grid cells belonging to the carbon fiber and phenolic resin matrix, and the fluid grid cells belonging to the pores, and assign the corresponding solid material properties to the solid grid cells.

[0016] Step 104: Based on the marked mesh elements and the assigned solid material properties, the ThermalFOAM solver configured to perform multiphysics coupled iterative solution is used in the OpenFOAM environment to perform multiphysics coupled iterative solution on the fluid and solid domains of the carbon / phenolic composite, and the coupled calculation results are obtained. Based on the coupled calculation results, the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite are calculated.

[0017] This application provides a three-dimensional pore-scale multi-field coupled prediction method for the pyrolysis ablation of carbon / phenolic composites. First, a porous medium geometric model of the carbon / phenolic composite material at the three-dimensional pore scale is constructed. Since real materials possess complex microstructures, including the distribution morphology of carbon fibers, the filling state of the phenolic resin matrix, and the connectivity of pores, these microscopic features directly influence subsequent flow and reaction behavior. Therefore, the core of this step lies in establishing a digital model that can accurately characterize these microstructures: the solid domain part uses spatial distribution information to clarify the diameter, length, and orientation of carbon fibers, as well as the spatial location of the phenolic resin matrix; the fluid domain part accurately describes the geometry and connectivity of pores using spatial distribution information, providing a realistic physical spatial basis for subsequent flow calculations.

[0018] After establishing the porous medium geometric model, a computational domain encompassing this model can be defined within the OpenFOAM environment, generating a structured mesh covering the entire model. OpenFOAM, as an open-source computational fluid dynamics platform, provides powerful mesh generation and numerical solution capabilities. The computational domain defined here reserves some flow development space upstream and downstream of the porous medium geometric model to ensure that the inlet and outlet boundary conditions do not affect the realism of the flow within the material. The fineness of the mesh node distribution determines the computational resolution; denser nodes allow for richer resolution of pore details, but also incur higher computational costs. The final generated structured mesh can consist of hexahedral elements, completely covering both the fluid and solid regions of the porous medium geometric model, discretizing the continuous physical space into a finite number of elements that a computer can process.

[0019] After the mesh is generated, each mesh cell can be attributed and labeled. Specifically, by comparing the center coordinates of the mesh cell with the spatial distribution information in the porous media geometric model, it can be accurately determined whether the mesh cell is located inside the carbon fiber, inside the phenolic resin, or in the pores, thus labeling it as the corresponding solid mesh cell or fluid mesh cell. For mesh cells labeled as solid, values ​​can be assigned to them; for example, carbon fiber cells can be assigned the density and reactivity characteristics of carbon fibers, and phenolic resin matrix cells can be assigned the density and pyrolysis characteristics of phenolic resin.

[0020] Based on the above labeling and assignment results, a specially developed ThermalFOAM solver is used for multiphysics coupled iterative solution. On each grid cell, the continuity equation, momentum equation, component transport equation, and energy equation are solved, corresponding to physical laws such as mass conservation, force balance, component diffusion, and energy conservation, respectively. Furthermore, on the solid grid cells, the remaining mass and solid volume fraction of the two materials are tracked according to the assigned material properties, thereby analyzing the solid-phase evolution results. In addition, the solver can establish a bidirectional data exchange mechanism at the fluid-solid interface. Specifically, the temperature and oxygen concentration of the flow field drive the solid-phase reaction, while the changes in pore structure caused by solid-phase consumption (such as changes in permeability and specific surface area) can be fed back to the flow field calculation in real time, affecting the flow and heat transfer in the next time step. After repeated iterations until convergence, the coupled calculation results are finally output. Based on the coupled calculation results, the ablation retreat distance and ablation retreat rate can be calculated, thus providing a quantitative evaluation index for the ablation resistance performance of the material.

[0021] By applying the technical solution of this embodiment, firstly, a porous medium geometric model of carbon / phenolic composite material at a three-dimensional pore scale is constructed. Then, in the OpenFOAM environment, a computational domain encompassing this porous medium geometric model is defined, and a structured mesh covering the entire porous medium geometric model is generated. After mesh generation, each mesh cell can be attribute-identified and labeled to determine whether it is a solid mesh cell or a fluid mesh cell. For mesh cells labeled as solid, corresponding solid-phase material properties can be assigned. Based on the above labeling and assignment results, a specially developed ThermalFOAM solver is used for multiphysics coupled iterative solution, ultimately outputting the coupled calculation results. Based on the coupled calculation results, the ablation retreat distance and ablation retreat rate can be calculated. This application's embodiments construct a porous media geometric model containing spatial distribution information of carbon fiber and phenolic resin matrix, and generate a structured mesh covering both fluid and solid domains based on the OpenFOAM platform. This overcomes the limitations of macroscopic homogenized models in failing to analyze pore-scale flow details and the actual evolution of the solid phase, achieving accurate characterization of the material's microstructure. By identifying the properties of the mesh elements and assigning corresponding solid-phase material properties to the carbon fiber and phenolic resin matrix, the technical blind spots commonly found in existing pore-scale simulation methods—which only consider the consumption of a single solid-phase material and fail to distinguish between the two different reaction mechanisms of carbon fiber oxidation and phenolic resin pyrolysis—are addressed. This allows for a true reflection of the differentiated ablation behavior of different components in the composite material. Based on the labeled mesh elements, a ThermalFOAM solver configured for multi-physics coupled iterative solution is used, effectively overcoming the shortcomings of traditional experimental methods in revealing the dynamic evolution mechanism of flow-solid-thermal-chemical processes during ablation, thus facilitating the acquisition of high-precision three-dimensional dynamic prediction results.

[0022] In this embodiment of the application, optionally, step 104, "based on the marked grid cells and the assigned solid material properties, using the ThermalFOAM solver configured to perform multiphysics coupled iterative solution in the OpenFOAM environment to perform multiphysics coupled iterative solution on the fluid domain and solid domain of the carbon / phenolic composite, and obtaining the coupled calculation results," includes: giving initial and boundary conditions for the flow field, and initializing the flow field parameters and temperature field according to the initial and boundary conditions; based on the initialized flow field parameters and temperature field, solving the fluid control equations on each grid cell to obtain the flow field parameters and temperature field at the current time step, wherein the fluid control equations include the continuity equation, component transport equation, energy equation, and momentum equation; determining the values ​​of each component in the fluid control equations at the current time step. If the computational residual of the equation is less than a preset convergence threshold, and the computational residual is greater than or equal to the preset convergence threshold, return to the step of solving the fluid control equations based on the flow field parameters and temperature field of the current time step; if the computational residual is less than the preset convergence threshold, update the solid material properties of each solid mesh element based on the flow field parameters and temperature field of the current time step to obtain the solid evolution result, and determine whether the current computation time is less than the preset total computation time; if yes, update the time step, and update the geometric parameters of the fluid control equations based on the solid evolution result, and return to the step of solving the fluid control equations based on the latest flow field parameters, temperature field, and the updated fluid control equations; if no, end the computation, and output the solid evolution results of all time steps as the final coupled computation result.

[0023] In this embodiment, firstly, initial and boundary conditions of the flow field are given, and the flow field parameters and temperature field are initialized accordingly. Here, the initial conditions may include the distribution of physical quantities such as velocity, pressure, temperature, and component concentration at the initial moment. Boundary conditions may define the physical behavior at the inlet, outlet, and wall of the computational domain, such as the incoming flow velocity, temperature, and oxygen concentration at the inlet, the zero-gradient condition at the outlet, and the no-slip or adiabatic condition at the wall. The initialization operation assigns these conditions to each grid cell within the computational domain, serving as the starting point for time-progressive calculations and providing initial flow field parameters and temperature fields for subsequent iterative solutions.

[0024] Next, based on the initialized flow field parameters and temperature field, the fluid control equations are solved on each grid cell to obtain the flow field parameters and temperature field at the current time step. Here, the fluid control equations consist of four core equations: the continuity equation describes mass conservation, ensuring that the net mass flow rate into and out of the grid cell is balanced with the internal mass change; the momentum equation describes the force balance relationship, including pressure gradient, viscous forces, and the Darcy drag term specific to porous media; the component transport equation describes the convection, diffusion, and consumption or generation of gaseous components (such as oxygen) due to solid-phase reactions; and the energy equation describes heat transfer, including convective heat transfer, heat conduction, and heat released or absorbed by solid-phase chemical reactions. These four equations constitute a complete mathematical framework describing the flow, mass transfer, and heat transfer behavior within the pores. Solving these equations yields the velocity field, pressure field, oxygen concentration field, and temperature field at the current time step.

[0025] After obtaining the solution result for the current time step, it is determined whether the computational residuals of each equation in the fluid control equations are less than a preset convergence threshold. The residual is a numerical indicator measuring the degree to which each equation is satisfied in the current iteration step; specifically, it refers to the difference between the left and right sides of the equation after substituting the currently calculated velocity, pressure, and temperature fields into the equation. If the residual is small, it indicates that the current solution satisfies the fluid control equations well; if the residual is large, it means that the solution is not yet accurate enough. The preset convergence threshold can be a very small positive number (e.g., 1 × 10⁻⁶). -9 This threshold serves as the standard for determining whether the solution has reached the required accuracy. Only when the computational residuals of all equations are less than this preset convergence threshold is the flow field calculation at the current time step considered to have converged.

[0026] If the calculated residuals are greater than or equal to the preset convergence threshold, it indicates that the flow field solution at the current time step has not yet met the accuracy requirements. Therefore, it is necessary to return and resolve the fluid control equations based on the flow field parameters and temperature field at the current time step. Here, returning to resolve represents an iterative process, where the velocity, pressure, oxygen concentration, and temperature fields calculated in the current iteration step are used as new initial values ​​and substituted back into the fluid control equations for solving. This process can be repeated until the calculated residuals of all equations are reduced below the preset convergence threshold. This internal iterative mechanism is crucial for ensuring the accuracy and stability of the nonlinear equation solution.

[0027] When the calculated residual is less than the preset convergence threshold, it indicates that the flow field calculation for the current time step has converged. At this point, based on the converged flow field parameters and temperature field, the solid material properties of each solid mesh element are updated to obtain the solid evolution results. Updating the solid material properties involves calculating the oxidation reaction rate of carbon fiber and the pyrolysis reaction rate of phenolic resin based on the current temperature field and oxygen concentration field, and then updating the remaining mass and solid volume fraction of both materials at the current time step. These solid evolution results reflect the degree of material consumption and microstructural changes during the ablation process. After completing the solid-phase update, determine whether the current calculation time is less than the preset total calculation time: if the preset total calculation time has not been reached, update the time step and update the geometric parameters (such as porosity, permeability, and specific surface area) in the fluid control equations based on the solid-phase evolution results just obtained. Then, based on the latest flow field parameters, temperature field, and updated geometric parameters, return to solve the fluid control equations for the next time step; if the preset total calculation time has been reached, end the calculation and output the solid-phase evolution results for all time steps as the final coupled calculation result.

[0028] This application's embodiments ensure the accuracy of the flow field solution within each time step through an internal iteration mechanism, and achieve bidirectional dynamic coupling between the flow field and solid-state evolution through external iteration, forming a hierarchical and logically rigorous solution framework. It clearly distinguishes the solid-state update step after flow field convergence, as well as the real-time update of geometric parameters based on solid-state evolution results, truly reflecting the feedback influence of pore structure changes on subsequent flow during ablation. Through dual control of preset convergence threshold and total computation time, it not only ensures the computational accuracy of each time step, but also realizes the time-progress simulation of the entire ablation process, and finally outputs the solid-state evolution results for all time steps, providing a complete data foundation for subsequent ablation performance evaluation.

[0029] Optionally, in this embodiment, the step of "solving the fluid control equations on each grid cell based on the initialized flow field parameters and temperature field to obtain the flow field parameters and temperature field at the current time step" includes: based on the initialized flow field parameters and temperature field, using a pressure-velocity coupling algorithm, iteratively solving the coupled equations of the continuity equation and the momentum equation on each grid cell until a preset convergence condition is met, thereby obtaining the converged velocity field and pressure field at the current time step. The momentum equation includes a Darcy term characterizing the flow resistance of the porous medium, and the coefficient of the Darcy term is updated in real time based on the solid-phase evolution results of the solid grid cells. Based on the initialized flow field parameters and temperature field, and the velocity field, the component transport equations are solved to obtain... The oxygen concentration field at the current time step, wherein the component transport equation includes an interface mass transfer source term, which is calculated based on the specific surface area of ​​the current solid grid cell and the solid-phase reaction rate, and the solid-phase reaction rate is calculated based on the temperature field and the oxygen concentration field in the initialized flow field parameters; based on the velocity field, oxygen concentration field and initialized temperature field at the current time step, the energy equation is solved to obtain the temperature field at the current time step, and the velocity field, pressure field and oxygen concentration field at the current time step are used as the flow field parameters at the current time step, wherein the energy equation includes a reaction heat source term, which acts on the fluid domain and solid domain at the fluid-solid interface to transfer heat bidirectionally between the fluid and solid.

[0030] In this embodiment, firstly, based on the initialized flow field parameters and temperature field, a pressure-velocity coupling algorithm is used to iteratively solve the coupled equations of the continuity equation and the momentum equation on each grid cell. When the iterative calculation meets the preset convergence condition, the velocity field and pressure field at this time are taken as the velocity field and pressure field of the current time step convergence. The pressure-velocity coupling algorithm is a standard method in computational fluid dynamics for handling incompressible or low Mach number flows. Since the continuity equation does not explicitly contain a pressure term, while the momentum equation has a pressure gradient term but no independent pressure equation, this embodiment uses this coupling algorithm to match the pressure field and the velocity field. Here, the momentum equation specifically introduces a Darcy term to characterize the flow resistance of porous media, which reflects the hindering effect of the pore structure on fluid flow. It is important to note that the coefficients of the Darcy term are not fixed, but are updated in real time based on the solid-phase evolution results of the solid grid cells. This is because as carbon fiber oxidation and phenolic resin pyrolysis proceed, the solid-phase material is continuously consumed, the pore structure changes, and the permeability changes accordingly. This real-time update mechanism ensures that the momentum equation can accurately reflect the dynamic influence of pore structure evolution on flow resistance during ablation.

[0031] After obtaining the converged velocity field, the component transport equation is further solved based on the initialized flow field parameters, temperature field, and the velocity field just solved to obtain the oxygen concentration field at the current time step. The component transport equation describes the convection, diffusion, and consumption or generation processes of gaseous components in the pores due to chemical reactions. The convection term depends on the velocity field, which is why the velocity field obtained in the previous step needs to be used as input. This equation contains a key interfacial mass transfer source term, which represents the consumption or generation of gaseous components by solid-phase chemical reactions, i.e., the oxygen consumption of the carbon fiber oxidation reaction and the mass exchange of gaseous products generated by the pyrolysis of phenolic resin at the fluid-solid interface. The calculation of this source term is based on the specific surface area of ​​the current solid mesh cells and the solid-phase reaction rate, where the specific surface area defines the fluid-solid interface area per unit volume, and the solid-phase reaction rate determines the reaction intensity per unit area. The calculation of the solid-phase reaction rate depends on the current temperature field and the oxygen concentration field in the initialized flow field parameters. Here, the initialized oxygen concentration field is used instead of the oxygen concentration field being solved at the current time step, which cleverly avoids circular dependence and reflects the explicit processing strategy in numerical solution.

[0032] After solving the component transport equations, the energy equations are solved based on the velocity field, oxygen concentration field, and initialized temperature field at the current time step to obtain the temperature field at the current time step. The energy equations describe the heat transfer process between the fluid and solid, including convective heat transfer caused by fluid motion, heat conduction due to temperature gradients, and heat release or absorption from chemical reactions. The convection term in the equations depends on the velocity field, while the reaction heat source term depends on the oxygen concentration field and temperature field. The oxygen concentration determines the intensity of the reaction, while the temperature affects the reaction rate constant. It is important to note that the reaction heat source term is set to act simultaneously on both the fluid and solid domains at the fluid-solid interface, indicating that the heat released or absorbed by solid-phase chemical reactions not only affects the temperature of the solid itself but can also be transferred to the surrounding fluid through convection and conduction; conversely, changes in fluid temperature can also affect the reaction rate of the solid phase. By simultaneously solving the energy equations over the global computational domain, bidirectional heat transfer and dynamic equilibrium between the fluid and solid are achieved without the need for additional artificial coupling at the interface. After solving the energy equation, the velocity field, pressure field, and oxygen concentration field of the current time step are used together as the flow field parameters for the current time step, providing complete input data for subsequent solid-phase evolution calculations and the advancement of the next time step.

[0033] In a specific embodiment, the continuity equation for the gas phase can be expressed as follows: ; in, t Indicates time, The porosity is the proportion of the fluid domain space within a single grid cell. For fluid density, For fluid velocity, This refers to the fluid / solid mass transfer rate that is non-zero only at the interface due to heterogeneous chemical reactions.

[0034] The momentum equation can be expressed as follows: ; in, It is the fluid pressure gradient. It is fluid viscosity. It is the permeability of the grid cells. This refers to Darcy's term.

[0035] As can be seen from the equations above, the continuity equation lacks a pressure term, and while the momentum equation contains a pressure gradient term, it lacks an independent pressure term. Therefore, the key to solving the pressure-velocity coupling problem lies in finding a pressure field such that the velocity field calculated from the momentum equation precisely satisfies the continuity equation.

[0036] Therefore, in the ThermalFOAM solver of OpenFOAM, the pressure-velocity coupled solution process in each time step is as follows: (1) Momentum prediction input: using the pressure field of the previous iteration step (or initial conditions). and velocity field Operation: [The operation involves] transferring the known... Substituting the pressure gradient term into the right-hand side of the momentum equation, discretizing and solving the equation yields a new velocity field. The solution at this point is obtained under an inaccurate old pressure field, so although it satisfies the momentum equation, it usually does not satisfy the continuity equation. (2) The solution obtained in the previous step Substituting this into the continuity equation, we find that the velocity divergence requires certain conditions to be met. Since the current velocity does not meet these conditions, a mass imbalance can occur. Through mathematical transformation, this mass imbalance can be linked to changes in pressure, thus deriving a pressure correction factor. The equation (i.e., the pressure Poisson equation). (3) Solution: Solve the pressure Poisson equation to obtain the pressure correction. Update the pressure field: Update the velocity field: Based on the new pressure gradient, explicitly correct the velocity field using the discrete form of the momentum equation, and obtain... At this point, we get and It is a pair of approximate solutions that simultaneously satisfy the continuity equation and the momentum equation. (4) Check convergence: Judgment: Check whether the changes (residuals) of the velocity field and pressure field before and after this iteration are less than the preset convergence condition. Loop: If it does not converge, update the... and As new initial values and Return to step (1) and start the next iteration; if convergence is achieved, it is considered that the velocity field and pressure field of the current time step have been accurately coupled, the inner iteration ends, and the calculation of the next physical field (such as the component transport equation) can begin.

[0037] It is important to note that during the pressure-velocity coupled calculation, i.e., the internal iteration process, only the pressure and velocity fields are continuously updated, while the other parameters remain unchanged. Furthermore, although the temperature field is not a directly solved variable in the pressure-velocity coupled process, it plays a crucial role in solving the velocity and pressure fields by influencing the physical properties of the fluid, the source term, and the structural parameters of the porous medium; therefore, it also needs to be included as input.

[0038] The above equation is obtained by introducing Porosity offers significant advantages: when calculating complex porous media structures, structured meshes can be used, and porosity values ​​can be assigned within the mesh to distinguish between the fluid and solid domains. This eliminates the need for traditional large numbers of unstructured meshes to describe complex porous structures, improving computational efficiency and accuracy. For example, it is effective across the entire computational domain for the fluid domain (i.e., the pore region), the solid domain, and the gas / solid interface. Specifically, the initial solid volume fraction in solid mesh cells can be set to 0.999 (not "1" to avoid singular value problems during calculation), while the initial solid volume fraction in fluid mesh cells can be set to 0.001. This allows the fluid and solid domains to be calculated using the same set of fluid governing equations.

[0039] In a specific embodiment, the component transport equation can be expressed as follows: ; in, Indicates the concentration of the reactant (i.e., oxygen). Indicates the effective diffusion coefficient. Represents the specific surface area in a grid cell. R This represents the reaction rate constant (which depends on temperature).

[0040] This step will use the known velocity field Substituting the velocity field (i.e., the output from the previous step) into the component transport equation, we can solve for the oxygen concentration field. In the above formula This represents the interfacial mass transfer source term, where the solid-phase reaction rate refers to the rate at which gas is consumed or generated per unit reaction interface due to a chemical reaction, and is related to... Relevant, specifically It is positively correlated with the solid-phase reaction rate.

[0041] solid-phase reaction rate The calculation formula can be expressed as: ; in, It is a function of the temperature field. It is a function of the oxygen concentration field.

[0042] It is important to note that in the above solution process, apart from the new velocity field, the oxygen concentration field and temperature field in the initialized flow field parameters being used as inputs, and the new oxygen concentration field being used as the output, the other parameters are fixed. Specifically, the oxygen concentration field in the initialized flow field parameters is substituted into the right side of the component transport equation, and the new oxygen concentration field can be obtained by solving the oxygen concentration field on the left side of the component transport equation.

[0043] In a specific embodiment, the energy equation can be expressed as follows: ; in, Indicates the mixed density. , This indicates the density of the fluid (gas in the pores). Indicates the solid density (weighted average or individual density of carbon fiber and phenolic resin matrix). The proportion of solid domain space within a grid cell. . T Represents the temperature field. It is a mixed heat capacity. , This represents the specific heat capacity at constant pressure of a fluid. This represents the specific heat capacity at constant pressure of a solid. It is the effective thermal conductivity. , Indicates the thermal conductivity of a fluid. The term represents the thermal conductivity of a solid; the term representing the heat source of the reaction. Used to describe the heat release and absorption resulting from heterogeneous chemical reactions.

[0044] This step will use the known velocity field Substituting the velocity field (i.e., the output from the previous step) into the energy equation, and simultaneously the oxygen concentration field... And the initialized temperature field is substituted into the energy equation to calculate the reaction heat source term. Thus, the temperature field T at the current time step can be obtained.

[0045] It should be noted that in the above solution process, the velocity field, oxygen concentration field, and old temperature field are the inputs, the new temperature field is the output, and the other parameters are fixed.

[0046] In this application, the iterative solution of the pressure-velocity coupling algorithm ensures the decoupling accuracy of the continuity equation and the momentum equation. The real-time update mechanism of the Darcy term coefficients enables the momentum equation to dynamically respond to the changes in pore structure caused by solid-phase evolution, realizing the initial coupling between flow field calculation and solid-phase evolution. In solving the component transport equation, the influence of solid-phase response on gaseous components is introduced in the form of a source term through the interface mass transfer source term, and the solid-phase reaction rate is calculated using the initialized oxygen concentration field, which avoids cyclic dependence and ensures the accuracy of the physical process. In solving the energy equation, the bidirectional action of the reaction heat source term at the fluid-solid interface realizes the natural transfer and dynamic balance of heat between the fluid and the solid, without the need for additional artificial coupling processing.

[0047] Optionally, in this embodiment, the step of "returning to solve the fluid control equations based on the flow field parameters and temperature field of the current time step" includes: updating the flow field parameters and temperature field used in the previous iteration based on the velocity field, pressure field, oxygen concentration field, and temperature field calculated at the current time step, to obtain updated flow field parameters and temperature field; recalculating the fluid density, fluid viscosity, and solid-phase reaction rate based on the updated temperature field and oxygen concentration field, and updating the physical property parameters of each equation in the fluid control equations according to the calculation results; substituting the updated flow field parameters and temperature field into the updated fluid control equations, and resolving the updated fluid control equations.

[0048] In this embodiment, when the calculated residual does not meet the preset convergence requirements, the fluid control equations need to be resolved to provide more accurate input data for the next iteration. Specifically, firstly, based on the velocity field, pressure field, oxygen concentration field, and temperature field obtained at the current time step, the flow field parameters and temperature field used in the previous iteration are updated to obtain the updated flow field parameters and temperature field. Here, the update is not simply discarding the old values, but using the new results calculated in this iteration as the initial guess values ​​for the next iteration. That is, in the iterative solution, each round of calculation is based on the output of the previous round. Therefore, using the result that is closer to the true solution in this round as the starting point for the next round allows the iteration process to gradually move towards convergence. This operation of replacing old values ​​with new values ​​reflects the essential characteristic of gradual approximation in iterative algorithms.

[0049] After updating the flow field parameters and temperature field, the fluid density, fluid viscosity, and solid-phase reaction rate are recalculated based on the updated temperature and oxygen concentration fields. The physical property parameters of each equation in the fluid control equation set are then updated according to the calculation results. Fluid density and viscosity are functions of temperature and composition; changes in the temperature field directly affect the values ​​of these physical property parameters. The solid-phase reaction rate also depends on the temperature and oxygen concentration fields; an increase in temperature typically leads to an exponential increase in the reaction rate. If these physical property parameters are not updated promptly after updating the temperature field, the coefficients in the momentum equation, component transport equation, and energy equation will not match the current physical state, causing the iterative process to fail to converge accurately. Therefore, this step of updating physical property parameters is essential to ensure that the physical processes described by the equations remain consistent with the current state.

[0050] Finally, the updated flow field parameters and temperature field are substituted into the fluid control equations containing the updated physical property parameters, and these equations are solved again. Here, substitution means using the updated velocity field, pressure field, oxygen concentration field, and temperature field as initial values, and reconstructing the coefficient matrix of the fluid control equations based on the updated fluid density, fluid viscosity, and solid-phase reaction rate, and then solving them again. Through this cycle of "parameter update—substitution—resolving," each iteration brings the solution one step closer to the true convergence value. This process can be repeated until the computational residuals of all equations are reduced to below the preset convergence threshold, at which point the flow field calculation for the current time step is considered to have converged, and the solid-phase update stage can begin.

[0051] This application's embodiments, by updating the flow field parameters, temperature field, and physical property parameters, not only ensure the evolution of the solution vector itself during the iteration process, but also ensure the consistency of the coefficients of the fluid control equations with the current physical state; by incorporating the update of the solid-phase reaction rate into the scope of physical property parameters, the chemical reaction, a strongly nonlinear factor, can be adjusted in a timely manner during the iteration process, avoiding convergence difficulties caused by reaction rate lag.

[0052] In an embodiment of this application, optionally, the step of "updating the solid material properties of each solid grid cell based on the flow field parameters and temperature field of the current time step to obtain solid phase evolution results" includes: for each solid grid cell, calculating the oxidation reaction rate of carbon fiber and the pyrolysis reaction rate of phenolic resin matrix based on the oxygen concentration field included in the flow field parameters of the current time step and the temperature field of the current time step; updating the remaining mass of carbon fiber in the current time step according to the oxidation reaction rate of carbon fiber; updating the remaining mass of phenolic resin matrix in the current time step according to the pyrolysis reaction rate of phenolic resin matrix; recalculating the solid volume fraction of each solid grid cell based on the updated remaining mass of carbon fiber and phenolic resin matrix, and using the remaining mass and the solid volume fraction as the solid phase evolution results; correspondingly, the step of "updating the geometric parameters of the fluid control equation set based on the solid phase evolution results" includes: recalculating the porosity, permeability and specific surface area of ​​each grid cell based on the latest solid phase evolution results of the current time step, and updating the geometric parameters of each equation in the fluid control equation set according to the calculation results.

[0053] In this embodiment, for each grid cell labeled as solid, the oxidation rate of carbon fiber and the pyrolysis rate of the phenolic resin matrix are calculated separately based on the oxygen concentration field and temperature field in the flow field parameters of the current time step. This separate calculation is necessary because carbon fiber and phenolic resin have distinctly different reaction mechanisms. Carbon fiber oxidation is a surface reaction, its rate controlled by both temperature and oxygen concentration, typically following an Arrhenius-form kinetic equation. In contrast, phenolic resin pyrolysis is a volumetric reaction, primarily driven by temperature, and its rate depends on the current temperature and the remaining resin content. This differentiated approach allows for precise tracking of the ablation behavior of different components in the composite material.

[0054] After obtaining the reaction rates of the two materials, the remaining mass of the carbon fiber at the current time step is updated based on the oxidation reaction rate of the carbon fiber. Specifically, in each time step, the mass consumed by the carbon fiber due to the oxidation reaction is equal to the reaction rate multiplied by the reaction area within that grid cell, and then multiplied by the time step length. This consumed mass is subtracted from the remaining mass of the previous time step to obtain the remaining mass of the current time step. This process intuitively reflects the manifestation of the fundamental physical law of mass conservation in solid-phase materials.

[0055] Similarly, the remaining mass of the phenolic resin matrix at the current time step is updated based on the pyrolysis reaction rate of the phenolic resin matrix. The pyrolysis process of phenolic resin differs from carbon fiber oxidation; it involves the decomposition of organic polymers into gaseous products and residual carbon at high temperatures, and its mass loss rate can be controlled by a pyrolysis reaction kinetic model. Through this step, the remaining mass of the phenolic resin matrix at the current time step is accurately updated, providing fundamental data for subsequent calculations of the solids volume fraction.

[0056] After updating the remaining mass of the two materials, the solid volume fraction of each solid grid cell is recalculated based on the updated remaining mass of the carbon fiber and phenolic resin matrix. The solid volume fraction is defined as the proportion of solid material within a grid cell. It is calculated by dividing the remaining mass of the carbon fiber and phenolic resin matrix by their respective material densities to obtain the volume occupied by each material within the current grid cell. These two volumes are then added together and divided by the total volume of the grid cell to obtain the solid volume fraction. Both the remaining mass and the solid volume fraction are output as solid-phase evolution results because these two quantities describe the state of the solid phase from different dimensions. Specifically, the remaining mass reflects the absolute consumption of material, while the solid volume fraction reflects the relative occupancy of material within the grid cell. Combining the two provides a complete description of the degree of solid-phase evolution.

[0057] Accordingly, after obtaining the solid-phase evolution results, it is necessary to recalculate the porosity, permeability, and specific surface area of ​​each grid cell based on the latest solid-phase evolution results at the current time step, and update the geometric parameters of each equation in the fluid control equation set according to the calculation results. Porosity is defined as the volume proportion of the fluid domain within a grid cell, directly derived from the solid volume fraction; specifically, porosity equals 1 minus the solid volume fraction. Permeability is a parameter describing the conductivity of porous media to fluids, usually positively correlated with porosity, and calculated from porosity and pore structure characteristics using empirical formulas (such as the Kozeny-Carman equation). Specific surface area is the fluid-solid interface area per unit volume, also dependent on the distribution pattern and remaining mass of the solid phase. These three geometric parameters will play a crucial role in the flow field calculation at the next time step. Porosity appears in the continuity and energy equations, permeability directly affects the Darcy drag term in the momentum equation, and specific surface area is used to calculate the interface mass transfer source term in the component transport equation. Through this update mechanism, the influence of solid-phase evolution on the flow field can be fed back in real time, forming a complete two-way coupled closed loop.

[0058] This application's embodiments achieve parallel calculations of two differentiated reaction mechanisms—carbon fiber oxidation and phenolic resin pyrolysis—within the same framework. This accurately reflects the ablation behavior of different components in the composite material, avoiding calculation errors caused by simplifying the composite material into a single material. Through a progressive update path of "remaining mass → solid phase evolution results → porosity / permeability / specific surface area," the microscopic results of solid phase evolution are gradually transformed into macroscopic geometric parameters that affect the flow field calculation. Through two levels of solid phase update and geometric parameter update, both the driving effect of the flow field on the solid phase (positive coupling) and the feedback effect of solid phase evolution on the flow field (reverse coupling) are reflected, thus achieving bidirectional dynamic coupling.

[0059] In this embodiment of the application, optionally, the computational domain has a fluid inlet boundary, and the coupled computational results include the solid-state evolution results of each time step; the step 104, "calculating the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite based on the coupled computational results," includes: for each time step, based on the solid-state evolution results of the time step, extracting the spatial coordinates of all grid nodes on the ablation surface, calculating the minimum and maximum distances between the ablation surface and the fluid inlet boundary according to the spatial coordinates, and taking the average of the maximum and minimum distances as the ablation retreat distance of the time step; performing numerical differentiation on the ablation retreat distances of different time steps to obtain the curve of the ablation retreat rate changing with time.

[0060] In this embodiment, the final step of the entire calculation process is to calculate the ablation retreat distance and ablation retreat rate based on the coupled calculation results, transforming the massive amount of simulated data into intuitive engineering evaluation indicators. Before starting the calculation, two prerequisites can be clarified: first, the computational domain has a fluid inlet boundary, indicating that the position of the ablated surface is measured relative to a fixed reference surface (i.e., the fluid inlet), thus ensuring the comparability of calculation results at different time steps; second, the coupled calculation results include the solid-state evolution results at each time step, which serves as the input basis for subsequent calculations, as the position information of the ablated surface is extracted from the solid-state evolution results. These two prerequisites lay the necessary physical benchmark and data foundation for the calculation of ablation indicators.

[0061] For each time step where calculations have been completed, the spatial coordinates of all mesh nodes on the ablation surface are extracted based on the solid-state evolution results of that time step. Here, the ablation surface refers to the interface between the solid material and the fluid domain, that is, the new boundary formed after the carbon fiber and phenolic resin matrix are consumed. In the three-dimensional structured mesh, these nodes are located adjacent to the solid mesh cells and the fluid mesh cells. After obtaining the coordinates of these nodes, the minimum and maximum distances from the ablation surface to the fluid inlet boundary are further calculated. The minimum distance reflects the position of the ablation surface closest to the inlet, and the maximum distance reflects the position of the ablation surface furthest from the inlet. These two extreme values ​​together characterize the spatial distribution range of the ablation surface. The average value between the maximum and minimum distances is taken as the ablation retreat distance for that time step, which is a direct indicator for evaluating the degree of material ablation.

[0062] After obtaining the ablation retreat distance for all time steps, the ablation retreat rate is further calculated. By numerically differentiating the ablation retreat distance for different time steps, the curve of the ablation retreat rate versus time can be obtained. Numerical differentiation allows estimation of the derivative from discrete data points; specifically, the change in ablation retreat distance between adjacent time steps can be calculated and divided by the time step size to obtain the instantaneous retreat rate at each moment. This curve visually reflects the dynamic characteristics of the material ablation process, such as whether the ablation retreat rate gradually increases or tends to level off, and at what moments abrupt changes occur. This information is of great value for understanding the ablation mechanism and evaluating material properties. Figure 2 As shown, a curve illustrating the change in ablation retreat rate over time is provided in an embodiment of this application.

[0063] This application's embodiments eliminate computational uncertainties caused by the choice of reference frame by using a fixed fluid inlet boundary as a reference benchmark, thus providing a unified metric for ablation retreat distance at different time steps. By incorporating the spatial coordinates of all nodes on the ablated surface into the calculation, and considering both minimum and maximum distances, it reflects the non-uniform evolution characteristics of the ablated surface and avoids the bias that may result from single-location measurements. Furthermore, by numerically differentiating the ablation retreat distance at discrete time steps, the static ablation retreat distance is transformed into a dynamic ablation retreat rate curve, providing a dynamic index that changes over time for evaluating the ablation resistance performance of materials.

[0064] In one specific embodiment, after determining the final coupling calculation results, the Iso-Surfaces function can be used to plot the solid-phase evolution results of the carbon fiber and phenolic resin matrix, such as... Figure 3As shown, streamlines are plotted using Streamtraces. Furthermore, the coupled calculation results can also include the velocity field, pressure field, temperature field, and oxygen concentration field at each time step. Therefore, the changes in the velocity field, pressure field, temperature field, and oxygen concentration field within the porous medium during ablation can be analyzed over time using the Contour function box.

[0065] Optionally, in this embodiment, step 104, "using a ThermalFOAM solver configured to perform multiphysics coupled iterative solutions in the OpenFOAM environment to perform multiphysics coupled iterative solutions on the fluid and solid domains of the carbon / phenolic composite," includes: spatially partitioning the computational domain using the ThermalFOAM solver to divide the structured mesh into multiple sub-regions, and assigning an independent computational process to each sub-region; establishing processor boundaries at adjacent interfaces between sub-regions using the ThermalFOAM solver, and establishing logical connections between sub-regions through the processor boundaries; for each time step, performing coupled iterative solutions on the fluid and solid domains within the corresponding sub-regions through each computational process, and after the solution is completed, exchanging field variable information at the interfaces of adjacent sub-regions at the processor boundaries through the communication mechanism built into the ThermalFOAM solver, wherein the field variable information includes velocity field, pressure field, temperature field, and oxygen concentration field.

[0066] In this embodiment, the computational domain is first spatially partitioned using the ThermalFOAM solver, dividing the originally complete structured mesh into multiple sub-regions, and assigning an independent computational process to each sub-region. Here, spatial partitioning refers to decomposing the entire computational domain into several blocks according to a certain strategy (such as uniform division along the coordinate axes) based on the distribution of mesh cells; each block is called a sub-region. Assigning an independent computational process to each sub-region allows the computational tasks of these sub-regions to be distributed to different CPU cores for simultaneous execution. This enables the massive computational load that would otherwise require a single processor to handle to be distributed across multiple processors for collaborative completion, thereby significantly improving computational efficiency.

[0067] After spatial partitioning, the ThermalFOAM solver establishes processor boundaries at the interfaces between adjacent sub-regions and uses these boundaries to establish logical connections between them. Here, processor boundaries are a logical concept, not corresponding to physical geometric boundaries. They refer to the mesh interfaces between different sub-regions that are spatially adjacent but belong to different processors in the computation process. When the computational domain is partitioned, previously adjacent mesh cells may be assigned to different sub-regions; the interfaces between these sub-regions are the processor boundaries. By establishing logical connections at these boundaries, the solver can identify which information needs to be exchanged between different computation processes, preparing for subsequent data communication.

[0068] For each time step, each computational process executes a coupled iterative solution of the fluid and solid domains within its corresponding sub-region. After completion, the process exchanges field variable information at the interface between adjacent sub-regions at the processor boundary via the built-in communication mechanism of the ThermalFOAM solver. Here, "executed separately" means that each computational process independently performs a complete coupled iterative solution process for its assigned sub-region, including flow field calculation, residual determination, and solid phase update. Since boundaries exist between sub-regions, mesh cells near these boundaries rely on information from corresponding cells in adjacent sub-regions (e.g., velocity, pressure, temperature). Therefore, after each time step, data exchange can be performed at the processor boundary via a built-in communication mechanism (e.g., based on the MPI message passing interface), ensuring that each computational process obtains the latest field variable information at the interface between adjacent sub-regions. The exchanged information includes the velocity field, pressure field, temperature field, and oxygen concentration field; these are core solution variables in the fluid control equations, ensuring the continuity and consistency of the solution across the global computational domain.

[0069] The embodiments of this application distribute large-scale computing tasks across multiple processors for parallel execution through domain decomposition and computational process allocation, significantly shortening the computation time for three-dimensional pore-scale simulation and making it feasible to perform fine simulations that were originally difficult to achieve due to excessive computational load. By establishing processor boundaries and communication mechanisms, the accuracy and efficiency of data exchange between sub-regions are guaranteed, ensuring that the results of parallel computation are completely consistent with those of serial computation and do not affect the solution accuracy.

[0070] Furthermore, as Figure 1 In terms of specific implementation, this application provides a three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupled prediction device, such as... Figure 4 As shown, the device includes: The geometric model construction module is used to construct a porous medium geometric model of carbon / phenolic composite at the three-dimensional pore scale. The porous medium geometric model includes spatial distribution information of the solid domain composed of carbon fiber and phenolic resin matrix, and spatial distribution information of the fluid domain composed of pores. The mesh generation module is used to set the computational domain containing the porous medium geometric model in the OpenFOAM environment, set the mesh node distribution in the computational domain, and generate a structured mesh covering the fluid domain and solid domain of the porous medium geometric model according to the mesh node distribution. The grid marking module is used to identify the attributes of each grid cell in the structured grid according to the spatial distribution information of the porous medium geometric model, mark the solid grid cells belonging to the carbon fiber and phenolic resin matrix, and the fluid grid cells belonging to the pores, and assign the corresponding solid phase material properties to the solid grid cells. The coupled solution module is used to perform multiphysics coupled iterative solution on the fluid and solid domains of the carbon / phenolic composite material in the OpenFOAM environment based on the marked mesh elements and the assigned solid material properties, using the ThermalFOAM solver configured to perform multiphysics coupled iterative solution. The coupled calculation results are obtained, and based on the coupled calculation results, the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite material are calculated.

[0071] Optionally, the coupled solution module is used for: Given the initial and boundary conditions of the flow field, the flow field parameters and temperature field are initialized according to the initial and boundary conditions; Based on the initialized flow field parameters and temperature field, the fluid control equations are solved on each grid cell to obtain the flow field parameters and temperature field at the current time step. The fluid control equations include the continuity equation, component transport equation, energy equation, and momentum equation. Determine whether the computational residuals of each equation in the fluid control equation set at the current time step are less than a preset convergence threshold; If the calculated residual is greater than or equal to the preset convergence threshold, the process returns to the step of solving the fluid control equations based on the flow field parameters and temperature field at the current time step. If the calculated residual is less than the preset convergence threshold, the solid material properties of each solid mesh element are updated based on the flow field parameters and temperature field at the current time step to obtain the solid evolution result. It is then determined whether the current calculation time is less than the preset total calculation time. If so, the time step is updated, and the geometric parameters of the fluid control equations are updated based on the solid evolution result. The solution to the fluid control equations is then returned based on the latest flow field parameters, temperature field, and the updated fluid control equations. If not, the calculation ends, and the solid evolution results for all time steps are output as the final coupled calculation result.

[0072] Optionally, the coupled solution module is further configured to: Based on the initialized flow field parameters and temperature field, a pressure-velocity coupled algorithm is used to iteratively solve the coupled equations of the continuity equation and the momentum equation on each grid cell until the preset convergence condition is met, and the velocity field and pressure field converged at the current time step are obtained. The momentum equation includes a Darcy term to characterize the flow resistance of porous media, and the coefficient of the Darcy term is updated in real time based on the solid phase evolution results of the solid grid cell. Based on the initialized flow field parameters, temperature field, and velocity field, the component transport equation is solved to obtain the oxygen concentration field at the current time step. The component transport equation includes an interface mass transfer source term, which is calculated based on the specific surface area of ​​the current solid mesh cell and the solid phase reaction rate. The solid phase reaction rate is calculated based on the temperature field and the oxygen concentration field in the initialized flow field parameters. Based on the velocity field, oxygen concentration field, and initialized temperature field of the current time step, the energy equation is solved to obtain the temperature field of the current time step. The velocity field, pressure field, and oxygen concentration field of the current time step are used as the flow field parameters of the current time step. The energy equation includes a reaction heat source term, which acts on the fluid and solid domains at the fluid-solid interface to transfer heat bidirectionally between the fluid and solid.

[0073] Optionally, the coupled solution module is further configured to: Based on the velocity field, pressure field, oxygen concentration field and temperature field obtained from the current time step, the flow field parameters and temperature field used in the previous iteration are updated to obtain the updated flow field parameters and temperature field. Based on the updated temperature field and oxygen concentration field, the fluid density, fluid viscosity and solid-phase reaction rate are recalculated, and the physical property parameters of each equation in the fluid control equation set are updated according to the calculation results. Substitute the updated flow field parameters and temperature field into the updated fluid control equations and solve the updated fluid control equations again.

[0074] Optionally, the coupled solution module is further configured to: For each solid mesh cell, the oxidation reaction rate of the carbon fiber and the pyrolysis reaction rate of the phenolic resin matrix are calculated based on the oxygen concentration field and the temperature field in the flow field parameters of the current time step. Update the remaining mass of the carbon fiber at the current time step based on the oxidation reaction rate of the carbon fiber. The remaining mass of the phenolic resin matrix at the current time step is updated based on the pyrolysis reaction rate of the phenolic resin matrix. Based on the remaining mass of the updated carbon fiber and phenolic resin matrix, the solid volume fraction of each solid grid cell is recalculated, and the remaining mass and the solid volume fraction are used as the solid phase evolution result. Accordingly, the coupled solution module is further configured to: Based on the latest solid-phase evolution results at the current time step, the porosity, permeability, and specific surface area of ​​each grid cell are recalculated, and the geometric parameters of each equation in the fluid control equation set are updated according to the calculation results.

[0075] Optionally, the computational domain has a fluid inlet boundary, and the coupled computational results include the solid-phase evolution results at each time step; the coupled solution module is further used for: For each time step, based on the solid-state evolution results of the time step, the spatial coordinates of all grid nodes on the ablation surface are extracted. According to the spatial coordinates, the minimum and maximum distances between the ablation surface and the fluid inlet boundary are calculated, and the average of the maximum and minimum distances is taken as the ablation retreat distance of the time step. Numerical differentiation was performed on the ablation retreat distance at different time steps to obtain the curve of ablation retreat rate as a function of time.

[0076] Optionally, the coupled solution module is further configured to: The ThermalFOAM solver spatially partitions the computational domain to divide the structured mesh into multiple sub-regions, and assigns an independent computational process to each sub-region. The ThermalFOAM solver establishes processor boundaries at adjacent interfaces between sub-regions, and establishes logical connections between sub-regions through these processor boundaries. For each time step, the coupled iterative solution of the fluid domain and solid domain in the corresponding sub-region is performed through each calculation process. After the solution is completed, the field variable information of the interface between adjacent sub-regions is exchanged at the processor boundary through the communication mechanism built into the ThermalFOAM solver. The field variable information includes velocity field, pressure field, temperature field and oxygen concentration field.

[0077] It should be noted that for other corresponding descriptions of the functional units involved in the three-dimensional pore-scale carbon / phenolic composite pyrolysis ablation multi-field coupling prediction device provided in the embodiments of this application, please refer to... Figures 1 to 3 The corresponding descriptions in the method will not be repeated here.

[0078] This application also provides a computer device, which may specifically be a personal computer, a server, a network device, etc. Figure 5 As shown, the computer device includes a bus, a processor, memory, and a communication interface, and may also include an input / output interface and a display device. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores an operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The database stores location information. The network interface allows communication with external terminals via a network connection. When the computer program is executed by the processor, it implements the steps in the various method embodiments.

[0079] Those skilled in the art will understand that Figure 5 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0080] In one embodiment, a computer-readable storage medium is provided, which may be non-volatile or volatile, having stored thereon a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0081] In one embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0082] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties.

[0083] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM). The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc., and are not limited to these.

[0084] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0085] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of this patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this application should be determined by the appended claims.

Claims

1. A three-dimensional pore-scale multi-field coupled prediction method for pyrolysis ablation of carbon / phenolic composites, characterized in that, include: A porous media geometric model of carbon / phenolic composite material at the three-dimensional pore scale is constructed. The porous media geometric model includes the spatial distribution information of the solid domain composed of carbon fiber and phenolic resin matrix, and the spatial distribution information of the fluid domain composed of pores. In the OpenFOAM environment, a computational domain containing the porous medium geometric model is set, and a grid node distribution is set within the computational domain. Based on the grid node distribution, a structured grid covering the fluid and solid domains of the porous medium geometric model is generated. Based on the spatial distribution information of the porous medium geometric model, attribute identification is performed on each grid cell in the structured grid. Solid grid cells belonging to carbon fiber and phenolic resin matrix are marked, as well as fluid grid cells belonging to pores. The solid grid cells are then assigned corresponding solid material properties. Based on the marked mesh elements and the assigned solid material properties, the ThermalFOAM solver, configured to perform multiphysics coupled iterative solution, is used in the OpenFOAM environment to perform multiphysics coupled iterative solution on the fluid and solid domains of the carbon / phenolic composite. The coupled calculation results are obtained, and based on the coupled calculation results, the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite are calculated.

2. The method according to claim 1, characterized in that, Based on the labeled mesh elements and the assigned solid material properties, the ThermalFOAM solver, configured to perform multiphysics coupled iterative solutions, is used in the OpenFOAM environment to perform multiphysics coupled iterative solutions on the fluid and solid domains of the carbon / phenolic composite, yielding coupled calculation results, including: Given the initial and boundary conditions of the flow field, the flow field parameters and temperature field are initialized according to the initial and boundary conditions; Based on the initialized flow field parameters and temperature field, the fluid control equations are solved on each grid cell to obtain the flow field parameters and temperature field at the current time step. The fluid control equations include the continuity equation, component transport equation, energy equation, and momentum equation. Determine whether the computational residuals of each equation in the fluid control equation set at the current time step are less than a preset convergence threshold; If the calculated residual is greater than or equal to the preset convergence threshold, the process returns to the step of solving the fluid control equations based on the flow field parameters and temperature field at the current time step. If the calculated residual is less than the preset convergence threshold, the solid material properties of each solid mesh element are updated based on the flow field parameters and temperature field at the current time step to obtain the solid evolution result. It is then determined whether the current calculation time is less than the preset total calculation time. If so, the time step is updated, and the geometric parameters of the fluid control equations are updated based on the solid evolution result. The solution to the fluid control equations is then returned based on the latest flow field parameters, temperature field, and the updated fluid control equations. If not, the calculation ends, and the solid evolution results for all time steps are output as the final coupled calculation result.

3. The method according to claim 2, characterized in that, Based on the initialized flow field parameters and temperature field, the fluid control equations are solved on each grid cell to obtain the flow field parameters and temperature field at the current time step, including: Based on the initialized flow field parameters and temperature field, a pressure-velocity coupled algorithm is used to iteratively solve the coupled equations of the continuity equation and the momentum equation on each grid cell until the preset convergence condition is met, and the velocity field and pressure field converged at the current time step are obtained. The momentum equation includes a Darcy term to characterize the flow resistance of porous media, and the coefficient of the Darcy term is updated in real time based on the solid phase evolution results of the solid grid cell. Based on the initialized flow field parameters, temperature field, and velocity field, the component transport equation is solved to obtain the oxygen concentration field at the current time step. The component transport equation includes an interface mass transfer source term, which is calculated based on the specific surface area of ​​the current solid mesh cell and the solid phase reaction rate. The solid phase reaction rate is calculated based on the temperature field and the oxygen concentration field in the initialized flow field parameters. Based on the velocity field, oxygen concentration field, and initialized temperature field of the current time step, the energy equation is solved to obtain the temperature field of the current time step. The velocity field, pressure field, and oxygen concentration field of the current time step are used as the flow field parameters of the current time step. The energy equation includes a reaction heat source term, which acts on the fluid and solid domains at the fluid-solid interface to transfer heat bidirectionally between the fluid and solid.

4. The method according to claim 3, characterized in that, The step of returning to solve the fluid control equations based on the flow field parameters and temperature field at the current time step includes: Based on the velocity field, pressure field, oxygen concentration field and temperature field obtained from the current time step, the flow field parameters and temperature field used in the previous iteration are updated to obtain the updated flow field parameters and temperature field. Based on the updated temperature field and oxygen concentration field, the fluid density, fluid viscosity and solid-phase reaction rate are recalculated, and the physical property parameters of each equation in the fluid control equation set are updated according to the calculation results. Substitute the updated flow field parameters and temperature field into the updated fluid control equations and solve the updated fluid control equations again.

5. The method according to claim 3 or 4, characterized in that, The process of updating the solid material properties of each solid mesh element based on the flow field parameters and temperature field at the current time step to obtain solid evolution results includes: For each solid mesh cell, the oxidation reaction rate of the carbon fiber and the pyrolysis reaction rate of the phenolic resin matrix are calculated based on the oxygen concentration field and the temperature field in the flow field parameters of the current time step. Update the remaining mass of the carbon fiber at the current time step based on the oxidation reaction rate of the carbon fiber. The remaining mass of the phenolic resin matrix at the current time step is updated based on the pyrolysis reaction rate of the phenolic resin matrix. Based on the remaining mass of the updated carbon fiber and phenolic resin matrix, the solid volume fraction of each solid grid cell is recalculated, and the remaining mass and the solid volume fraction are used as the solid phase evolution result. Accordingly, updating the geometric parameters of the fluid control equations based on the solid-phase evolution results includes: Based on the latest solid-phase evolution results at the current time step, the porosity, permeability, and specific surface area of ​​each grid cell are recalculated, and the geometric parameters of each equation in the fluid control equation set are updated according to the calculation results.

6. The method according to claim 1, characterized in that, The computational domain has a fluid inlet boundary, and the coupled computation results include the solid-state evolution results at each time step; the calculation of the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite based on the coupled computation results includes: For each time step, based on the solid-state evolution results of the time step, the spatial coordinates of all grid nodes on the ablation surface are extracted. According to the spatial coordinates, the minimum and maximum distances between the ablation surface and the fluid inlet boundary are calculated, and the average of the maximum and minimum distances is taken as the ablation retreat distance of the time step. Numerical differentiation was performed on the ablation retreat distance at different time steps to obtain the curve of ablation retreat rate as a function of time.

7. The method according to claim 1, characterized in that, The process of using the ThermalFOAM solver configured to perform multiphysics coupled iterative solutions in the OpenFOAM environment to solve the fluid and solid domains of the carbon / phenolic composite includes: The ThermalFOAM solver spatially partitions the computational domain to divide the structured mesh into multiple sub-regions, and assigns an independent computational process to each sub-region. The ThermalFOAM solver establishes processor boundaries at adjacent interfaces between sub-regions, and establishes logical connections between sub-regions through these processor boundaries. For each time step, the coupled iterative solution of the fluid domain and solid domain in the corresponding sub-region is performed through each calculation process. After the solution is completed, the field variable information of the interface between adjacent sub-regions is exchanged at the processor boundary through the communication mechanism built into the ThermalFOAM solver. The field variable information includes velocity field, pressure field, temperature field and oxygen concentration field.

8. A three-dimensional porous-scale carbon / phenolic composite pyrolysis ablation multi-field coupled prediction device, characterized in that, include: The geometric model construction module is used to construct a porous medium geometric model of carbon / phenolic composite at the three-dimensional pore scale. The porous medium geometric model includes spatial distribution information of the solid domain composed of carbon fiber and phenolic resin matrix, and spatial distribution information of the fluid domain composed of pores. The mesh generation module is used to set the computational domain containing the porous medium geometric model in the OpenFOAM environment, set the mesh node distribution in the computational domain, and generate a structured mesh covering the fluid domain and solid domain of the porous medium geometric model according to the mesh node distribution. The grid marking module is used to identify the attributes of each grid cell in the structured grid according to the spatial distribution information of the porous medium geometric model, mark the solid grid cells belonging to the carbon fiber and phenolic resin matrix, and the fluid grid cells belonging to the pores, and assign the corresponding solid phase material properties to the solid grid cells. The coupled solution module is used to perform multiphysics coupled iterative solution on the fluid and solid domains of the carbon / phenolic composite material in the OpenFOAM environment based on the marked mesh elements and the assigned solid material properties, using the ThermalFOAM solver configured to perform multiphysics coupled iterative solution. The coupled calculation results are obtained, and based on the coupled calculation results, the ablation retreat distance and ablation retreat rate of the carbon / phenolic composite material are calculated.

9. A storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 7.

10. A computer device, comprising a storage medium, a processor, and a computer program stored on the storage medium and executable on the processor, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 7.

Images