Analysis Method and System for Ignition State of Space Monopropellant Engine

Abstract

The present application relates to the technical field of data analysis, and provides a method and system for analyzing the ignition state of a space monopropellant engine. The method includes: obtaining the initial temperature parameter, initial pressure parameter, and propellant flow parameter of the combustion chamber of the space monopropellant engine in the ignition state; generating a flow domain grid model including the characteristics of porous medium effect and chemical reaction effect based on the initial temperature parameter, initial pressure parameter, and propellant flow parameter of the combustion chamber; setting dynamic boundary conditions corresponding to the ignition transient process in the flow domain grid model, and performing a full-cycle numerical simulation on the internal temperature field distribution of the space monopropellant engine to generate a temperature field dynamic change sequence in the ignition state; extracting the thermal stress evolution data of the key components of the space monopropellant engine according to the temperature gradient distribution characteristics at different moments in the temperature field dynamic change sequence, and combining the thermal stress evolution data with a preset overheat failure threshold to output the performance stability evaluation result in the ignition state.

Description

Technical Field

[0001] This application belongs to the technical field of data analysis, and particularly relates to a method and system for analyzing the ignition state of a spatial single-component engine. Background Art

[0002] In the field of research and development and performance optimization of spatial single-component engines, the thermodynamic characteristic analysis and structural reliability assessment of the ignition transient process are core technical challenges. In the prior art, numerical simulations of the engine ignition process usually adopt simplified physical models. Although such methods can reduce the computational complexity, they ignore the non-uniform influence of the actual pore distribution on local flow, heat transfer, and chemical reaction rates, resulting in a significant limitation in the prediction accuracy of the temperature field and pressure field.

[0003] In terms of boundary condition setting, existing research is mostly based on steady-state or quasi-steady-state assumptions. For example, the combustion chamber wall temperature is set to a fixed value or an empirical temperature rise curve is adopted, which cannot truly reflect the dynamic heat exchange process between the wall and the high-temperature gas during the ignition process. Regarding thermal stress analysis and performance evaluation, the prior art usually relies on simplified thermoelastic models, which cannot capture the dynamic thermal stress changes caused by high-frequency temperature fluctuations during the ignition transient process, and lack quantitative correlation analysis of material property degradation (such as high-temperature creep, fatigue crack propagation). In addition, the acquisition of propellant flow parameters in the prior art mostly relies on off-line calibration data, without considering the real-time influence of supply system pressure fluctuations on the mass flow rate during the ignition transient process, which will lead to a systematic deviation between the initial conditions of the numerical model and the actual working conditions. Especially in the propellant flow regulation stage (such as valve opening delay, pipeline pressure oscillation), the simulation results cannot effectively reflect the dynamic response characteristics of the engine.

[0004] In summary, the existing spatial single-component engine simulation and evaluation technologies have defects such as insufficient multi-physical field coupling accuracy, lack of characterization of dynamic boundary conditions, and rough analysis of thermal stress evolution. Summary of the Invention

[0005] This application provides a method and system for analyzing the ignition state of a spatial single-component engine to improve the engineering applicability of performance prediction and reliability assessment of spatial single-component engines.

[0006] In a first aspect, an embodiment of the present application provides a method for analyzing the ignition state of a space monopropellant engine, which is applied to an engine ignition state analysis system. The method includes: obtaining the initial temperature parameter, initial pressure parameter, and propellant flow parameter of the combustion chamber of the space monopropellant engine in the ignition state; generating a flow domain grid model including the characteristics of porous medium effect and chemical reaction effect based on the initial temperature parameter, initial pressure parameter, and propellant flow parameter of the combustion chamber; setting dynamic boundary conditions corresponding to the ignition transient process in the flow domain grid model, and performing a full-cycle numerical simulation on the internal temperature field distribution of the space monopropellant engine to generate a temperature field dynamic change sequence in the ignition state; extracting the thermal stress evolution data of the key components of the space monopropellant engine according to the temperature gradient distribution characteristics at different moments in the temperature field dynamic change sequence, and combining the thermal stress evolution data with a preset overheat failure threshold to output the performance stability evaluation result in the ignition state.

[0007] In a second aspect, an embodiment of the present application provides an engine ignition state analysis system, which includes a processor and a memory. Among them, the memory stores a computer program, and when the computer program is executed by the processor, the processor is caused to execute the steps of the above method.

[0008] In a third aspect, an embodiment of the present application provides a computer-readable storage medium, which includes a computer program. When the computer program runs on an engine ignition state analysis system, the computer program is used to cause the engine ignition state analysis system to execute the steps of the above method.

[0009] In the implementation of the present application, it is possible to deeply integrate microscopic structure characteristics, transient chemical reaction kinetics, and dynamic thermodynamic boundaries to achieve high-precision ignition state analysis and processing, thereby improving the engineering applicability of performance prediction and reliability assessment of space monopropellant engines, and improving the defects of insufficient multi-physical field coupling accuracy, missing representation of dynamic boundary conditions, and rough thermal stress evolution analysis. Description of the Drawings

[0010] Figure 1 It is a schematic flow chart of a method for analyzing the ignition state of a space monopropellant engine provided by an embodiment of the present application.

[0011] Figure 2 It is a schematic structural diagram of an engine ignition state analysis system provided by an embodiment of the present application. Detailed Embodiments

[0012] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the following will clearly and completely describe the technical solutions of this application in conjunction with the accompanying drawings in the embodiments of this application. Obviously, the described embodiments are part of the technical solutions of this application, rather than all of them. Based on the embodiments described in this application document, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope protected by the technical solutions of this application.

[0013] See Figure 1 , which is a method for analyzing the ignition state of a space single-component engine provided in the embodiments of this application. This method can be applied to an engine ignition state analysis system, and the specific process is as follows in steps 101 - step 104:

[0014] Step 101: Obtain the initial chamber temperature parameter, initial pressure parameter, and propellant flow parameter of the space single-component engine in the ignition state.

[0015] In the embodiments of this application, the chamber temperature distribution data at the initial moment of ignition can be obtained through a platinum resistance temperature sensor array installed on the chamber wall of the space single-component engine. This platinum resistance temperature sensor array includes 12 groups of PT100 type temperature sensors distributed along the axial direction of the chamber, with a measurement point spacing of 5 mm and a measurement range covering 300K to 1200K.

[0016] For another example, a piezoelectric dynamic pressure sensor can be synchronously used to collect the chamber head pressure parameter. The sensor model is Kistler603B, and the sampling frequency is set to 100kHz to ensure capturing the pressure fluctuation characteristics at the moment of ignition.

[0017] For yet another example, the propellant flow parameter is monitored in real time through a Coriolis mass flowmeter. This flowmeter is installed downstream of the pressure stabilizing section of the propellant supply pipeline, with a measurement accuracy reaching ±0.5% of the reading, and can accurately obtain the mass flow, temperature, and density parameters of the hydrazine-based propellant. Experimentally measured that the initial chamber temperature under typical working conditions is 650 ± 15K, the initial pressure is 2.1MPa, and the propellant mass flow is 12.5g / s. The above parameters will be used as the input reference values for the subsequent numerical model.

[0018] In the embodiments of this application, the propellant flow parameter refers to a set of key physical quantities that characterize the motion state of the propellant in the engine flow channel, including parameters such as mass flow, temperature, density, and flow velocity distribution. For example, a Coriolis mass flowmeter (model Rheonik RHM015) can be used for accurate measurement. This device is based on the principle of the Coriolis force effect and calculates the mass flow by measuring the phase difference caused by the propellant flowing through the U-shaped measuring tube, with a measurement accuracy of ±0.15% of the range.

[0019] In specific implementation, a double straight-tube sensor is installed downstream of the pressure-stabilizing section of the hydrazine-based propellant supply pipeline to ensure stable data acquisition after the flow is fully developed. Exemplary data includes the propellant inlet temperature (293±2K), density (1.004 g / cm³), and mass flow rate (12.5±0.3 g / s). These parameters are transmitted to the control terminal through the Modbus RTU protocol and used as the initial conditions for the mass source term and energy function in the computational fluid dynamics (CFD) model. For example, when setting the material properties in the Fluent software, the measured density value is substituted into the state function, and the mass flow rate value is converted into the inlet velocity boundary condition to ensure the physical authenticity of the propellant supply process.

[0020] Step 102: Based on the initial temperature parameter of the combustion chamber, the initial pressure parameter, and the propellant flow parameter, generate a flow domain grid model including the characteristics of the porous medium effect and the chemical reaction effect.

[0021] In the embodiment of the present application, when constructing the flow domain grid model based on the measured parameters, first, the three-dimensional laser scanning technology is used to obtain the three-dimensional point cloud data of the engine entity structure, and a geometric model including the catalytic bed, injector, and combustion chamber cavity is generated through reverse engineering software.

[0022] For the porous medium effect, the microfocus X-ray computed tomography technology is used to obtain the true pore structure of the catalyst bed layer. The reconstructed three-dimensional pore model is imported into the COMSOL Multiphysics software, and the random packing algorithm is applied to generate a non-uniform porous medium grid, with the element size controlled within 50 μm to accurately characterize the porosity gradient distribution.

[0023] In addition, for the chemical reaction effect modeling, the detailed kinetic mechanism is adopted. The hydrazine decomposition reaction path is decomposed into 12 elementary reaction steps such as N2H4→NH2+NH2 and NH2→NH+H, and a finite rate chemical reaction model is established in the Fluent software by importing in the CHEMKIN format.

[0024] The finally generated hybrid grid model contains approximately 8.5 million hexahedral elements. Among them, the catalytic bed area adopts an adaptive refined grid with a minimum element size of 20 μm, and the injector swirl channel adopts an O-type topological structure grid to ensure the flow direction analysis accuracy.

[0025] In step 102, the porous medium effect characteristic refers to the influence characteristics of the complex pore structure in the catalyst bed on the fluid flow and heat and mass transfer processes. In the embodiments of the present application, microfocus X-ray computed tomography (μCT, model ZEISS Xradia 520 Versa) is used to obtain the true three-dimensional pore structure of the iridium / aluminum oxide catalyst, and the scanning resolution reaches 2 μm / voxel. For another example, through image processing with Avizo Fire software, key parameters such as porosity (ε = 0.38), specific surface area (S v = 1.2×10^4 m² / m³), and tortuosity (τ = 1.85) are extracted. When establishing an asymmetric porous medium model in COMSOL Multiphysics, the Brinkman function is applied to couple the Forchheimer correction term to describe the flow resistance as: .

[0026] Among them, represents the vector differential operation, is the divergence operator, is the Laplace operator; ρ (fluid density): unit: kg / m 3 , represents the mass density of the fluid; u (velocity vector): unit: m / s, represents the velocity field of the fluid; is the tensor product (outer product) of the velocity vector; p (pressure): unit: Pa (Pascal), represents the hydrostatic pressure of the fluid; μ (dynamic viscosity): unit: Pa·s, characterizes the viscous resistance of the fluid; K (permeability): unit: m 2 , describes the ability of the porous medium to allow fluid to pass through, and is related to the pore structure; β (inertia coefficient): unit: m -1 , the coefficient of the Forchheimer correction term, characterizes the influence of the inertial effect on the resistance at high flow rates; (magnitude of velocity): unit: m / s, the magnitude of the velocity vector.

[0027] Furthermore, the permeability K is calculated by the Kozeny-Carman formula, and the inertia coefficient β is determined by the Ergun function. In the grid generation stage, the Voronoi tessellation algorithm is used to create an unstructured grid, and the cell size gradient is controlled within 1 / 5 of the pore diameter (i.e., 50 μm) to ensure the accurate capture of local flow separation phenomena.

[0028] In addition, the chemical reaction effect characteristic refers to the coupled influence of the exothermic decomposition reaction occurring in the catalyst bed on the flow and heat transfer of the propellant. In the embodiments of the present application, a detailed chemical reaction mechanism is adopted, and the decomposition process of hydrazine (N2H4) is decomposed into 12 elementary reaction steps, including:

[0029] 1. N2H4 → 2NH2 (Ea = 58.6 kJ / mol)

[0030] 2. NH2 → NH + H (Ea = 102.4 kJ / mol) ...

[0031] The rate constants of each step reaction are calculated by the transition state theory and imported into the finite rate chemical reaction model of Fluent in the CHEMKIN-II format. A surface reaction model is set in the catalytic bed region, and the density of iridium catalyst active sites is defined as 5.6×10^18 sites / m². The apparent reaction rate is calculated by combining the Langmuir-Hinshelwood adsorption mechanism. For example, the expression of the catalytic reaction rate in the UDF is as follows:

[0032] .

[0033] where r is the catalytic reaction rate; k0 is the pre-exponential factor, which is the kinetic constant in the Arrhenius equation and is related to the molecular collision frequency and orientation factor; Ea is the apparent activation energy; R is the universal gas constant, which is a constant connecting the thermodynamic temperature and the energy scale; T is the catalyst surface temperature, which is used for the local temperature of the catalytic bed and directly affects the reaction rate; P N2H4 is the partial pressure of hydrazine (N2H4), that is, the partial pressure of the reactant in the gas phase, which drives the reaction forward. P NH3 is the partial pressure of ammonia (NH3), which is used for the partial pressure of the reaction product in the gas phase and inhibits the reaction rate.

[0034] Furthermore, , Ea = 67.8 kJ / mol, K ads is the ammonia adsorption equilibrium constant. This model can accurately predict the phenomenon of sudden local temperature rise caused by reaction heat release. For example, the simulation shows that the temperature in the hot spot area of the catalytic bed rises from 650 K to 1250 K within 20 ms after ignition.

[0035] Step 103: Set the dynamic boundary conditions corresponding to the ignition transient process in the flow domain grid model, and perform a full-cycle numerical simulation on the internal temperature field distribution of the space monopropellant engine to generate the dynamic change sequence of the temperature field in the ignition state.

[0036] In step 103, when setting the dynamic boundary conditions, the combustion chamber wall boundary can be defined as a time-varying temperature boundary condition according to the exemplary data of the engine ignition transient process. A cubic spline interpolation function is used to fit the nonlinear change process of the wall temperature rising from 300K to 800K within 0 - 500ms after ignition. For example, the propellant inlet boundary is set as a pressure inlet condition, and the given transient pressure curve includes a 0.5MPa step increase in the initial stage and a subsequent linear growth segment of 1.2MPa / s. This curve is obtained by fitting the measured pressure data.

[0037] Furthermore, a coupled implicit solver is used for the numerical simulation. The time step is set to 1μs to meet the CFL condition for acoustic wave propagation. When solving the N - S function, the k - omega SST turbulence model is enabled and the DO radiation model is activated. In the catalytic bed region, the porous media momentum source term and a user - defined chemical reaction rate UDF are enabled. The activation energy of the catalyst surface reaction is corrected as a function of temperature through a user - defined function. The full - cycle simulation lasts for 500ms of physical time, and the full - field temperature data is output every 0.1ms, finally generating a data sequence containing 5000 transient fields, completely recording the temperature evolution characteristics in the injector swirl region, the catalytic bed hot spot region, and the nozzle throat during the ignition process.

[0038] It can be understood that the temperature field dynamic change sequence refers to a three - dimensional temperature distribution data set recorded in chronological order in the entire engine domain during the ignition process. In the embodiment of the present application, a coupled solver is used for transient simulation, and the time step is set to 1μs to meet the Courant - Friedrichs - Lewy (CFL) condition (Co < 1). Adaptive time step control is enabled in the catalytic bed region. When the temperature change rate exceeds 50K / ms, the time step is automatically reduced to 0.2μs. The data output is stored in the HDF5 format. Each time slice contains the temperature values of 1.2×10^6 grid nodes, and the spatial resolution reaches 0.1mm. The typical data sequence shows that a 340K low - temperature core flow is formed in the injector swirl region 5ms after ignition, while a 1560K high - temperature mass with a diameter of 2mm appears at the outlet of the catalytic bed at 50ms. This high - temperature mass migrates downstream with the flow and reaches the nozzle throat at 120ms.

[0039] Step 104: Extract the thermal stress evolution data of the key components of the spatial single - component engine according to the temperature gradient distribution characteristics at different moments in the temperature field dynamic change sequence, and combine the thermal stress evolution data with a preset overheat failure threshold to output the performance stability evaluation result in the ignition state.

[0040] In step 104, in the thermal stress analysis stage, the temperature field data of the dynamic change sequence of the temperature field is first mapped to the finite element model of ANSYS Mechanical, and the finite element model includes the detailed structural characteristics of the injector, the combustion chamber, and the nozzle assembly. For the injector head made of GH3128 superalloy, the temperature gradient distribution data is extracted and substituted into the thermoelastic constitutive function to calculate the transient thermal stress field.

[0041] It should be noted that special attention can be paid to the maximum temperature gradient area that appears at the root of the injector swirler blade at 50 ms after ignition. The axial temperature gradient in this area reaches 1.2×10^5 K / m, and the calculated equivalent thermal stress is 485 MPa. By comparing the calculation results with the material high-temperature yield strength database, it is found that when the yield strength of the nickel-based alloy of the injector head is 620 MPa in an 850 K environment, the actual thermal stress value reaches 78.2% of the material strength threshold. The performance stability evaluation module triggers an alarm according to the preset overheat failure criterion when the local temperature exceeds the material recrystallization temperature (923 K) or the thermal stress exceeds 80% of the yield strength.

[0042] For example, the simulation results show that a transient temperature spike (peak temperature 935 K) lasting 0.8 ms appears at the catalytic bed support ring at 320 ms after ignition. The engine ignition state analysis system automatically marks this area as a high-risk area and generates a three-dimensional thermal load spectrum including the duration and spatial distribution of temperature overrun, providing a quantitative basis for subsequent structural optimization.

[0043] In the embodiment of the present application, the thermal stress evolution data refers to the quantification result of the time-varying stress field generated by the temperature gradient of the engine structural components. In the embodiment of the present application, the CFD temperature field data is interpolated into the structural grid of ANSYS Mechanical through a field mapping algorithm, and the stress is calculated using the thermoelastic constitutive function: σ ij =C ijkl (ε kl -αΔTδ kl ). Among them, C ijkl is the stiffness tensor of the GH3128 alloy, and α is the coefficient of thermal expansion. When performing a transient analysis on the injector head, the maximum temperature gradient area at 50 ms after ignition (at coordinates X = 35 mm, Y = 12 mm) is extracted. The axial gradient in this area reaches 1.2×10^5 K / m, and the calculated Von Mises stress is 485 MPa. The thermal fatigue crack growth rate is calculated by the J integral method. When ΔT > 600 K, the crack growth rate da / dN reaches 2.3×10^-8 m / cycle, and this data is used to evaluate the remaining life of the injector after 1000 ignition cycles.

[0044] Furthermore, the performance stability evaluation result refers to the quantitative safety criterion analysis report based on thermodynamic parameters and material limits. The embodiments of the present application can set up a two - level early warning mechanism: the triggering condition of the first - level early warning is that the local temperature exceeds the recrystallization temperature of the material (923K) or the thermal stress reaches 75% of the yield strength (620MPa); the second - level early warning condition is that the temperature exceeds 90% of the melting point (1623K) or the stress exceeds 90% of the yield strength. The evaluation algorithm adopts the moving time - window statistical method. When a temperature spike (935K lasting for 0.8ms) appears at 320ms in the catalytic bed support ring, the system marks this event as a third - level anomaly (the cumulative over - standard time < 1ms). The final output is an evaluation report including the coordinates of the risk area, the list of over - standard parameters, and the recommended improvement measures (such as increasing the cooling channel density at this place), providing data support for the optimized design of the engine.

[0045] Thus, through the establishment of a high - precision multi - physical - field coupling model, the embodiments of the present application have realized the refined simulation of the ignition process of the space single - component engine. The obtained dynamic change sequence of the temperature field and the evolution data of the thermal stress provide key theoretical support for verifying the fatigue life of the injector structure, optimizing the porosity distribution of the catalytic bed, and improving the ignition timing control strategy. In addition, the quantitative analysis of the degradation of material properties under transient thermal shock can directly guide the thermal protection design of the key components of the engine, effectively preventing problems such as catalyst sintering failure or structural creep deformation caused by local overheating.

[0046] As an optional but not limited embodiment, the generating of the flow - domain grid model including the characteristics of porous - medium effect and chemical - reaction effect based on the initial temperature parameter, the initial pressure parameter, and the propellant flow parameter in step 102 includes:

[0047] Step 1021: Perform three - dimensional discretization modeling on the pore structure of the catalytic bed of the space single - component engine based on the porous - medium model to generate the topological data of the pore distribution of the catalytic bed.

[0048] For example, a micro - focus X - ray computed tomography system (model ZEISS Xradia 520Versa) can be used to perform three - dimensional structure imaging on the iridium / aluminum oxide catalyst bed. The scanning parameters are set as voltage 80kV, current 88μA, and the spatial resolution reaches 2μm, obtaining a tomographic image sequence containing 1024×1024×2048 voxels. Through image post - processing with Avizo Fire software, the region - growing algorithm is used to segment pores and solid skeletons, and the porosity distribution cloud map and the equivalent pore diameter statistical histogram are calculated.

[0049] For another example, based on the fractal geometry theory, a random walk algorithm can be applied to reconstruct a three-dimensional pore network model of the catalyst bed, generating topological connection data containing 1.2×10^6 pore nodes. This topological data records the coordinate positions, equivalent diameters (ranging from 15 to 85 μm), and adjacent pore connectivity information of each pore unit. After being exported as an STL format file and imported into the COMSOL Multiphysics software, a subtraction operation of the pore model and the catalytic bed solid geometry is performed through Boolean operations, finally forming a discretized digital twin of the porous medium, with the relative error of the porosity gradient distribution less than 3.5% compared to the exemplary data.

[0050] Step 1022: Combine the propellant decomposition reaction rate sub-model in the chemical reaction model, and spatially couple the topological data of the catalytic bed pore distribution with the combustion chamber geometry to generate an initial flow domain grid for multi-physics coupling.

[0051] In the embodiment of the present application, the chemical reaction coupling process is realized through a user-defined function, specifically by embedding the topological data of the catalytic bed pores in the chemical reaction model of the Fluent software. First, 12 elementary reaction steps of the hydrazine decomposition reaction are converted into a kinetic input file in CHEMKIN format, where the surface reaction steps are described by the Langmuir-Hinshelwood mechanism, and the iridium catalyst active site density is defined as 5.6×10^18 sites / m². In the mesh generation stage, the ANSYS ICEM CFD software is used to perform unstructured mesh generation on the combustion chamber geometry. In the catalytic bed area, the pore topological data is mapped into a porous medium sub-domain, and the permeability tensor is set as an anisotropic parameter (axial permeability 3.2×10^-10 m², radial permeability 1.8×10^-10 m²).

[0052] Furthermore, multi-physics coupling can be achieved through an APDL script, bidirectionally coupling the chemical reaction source term and the porous medium momentum function. Among them, the momentum loss term in the porous medium region is corrected by the Ergun function model, and the chemical reaction heat release power is calculated in real time through the UDF and injected as the energy function source term. The finally generated hybrid mesh model includes a hexahedron-dominated mesh at the injector head (element size 0.2 mm), a tetrahedron adaptive mesh in the catalytic bed area (minimum element size 20 μm), and a prism layer mesh in the combustion chamber transition area (total number of layers 15, growth rate 1.2). The overall number of mesh elements reaches 11.2 million, and the Jacobian factor is greater than 0.82.

[0053] Step 1023: According to the mass flow rate distribution in the propellant flow parameters, locally refine the mesh of the initial flow domain to generate a flow domain mesh model that meets the preset accuracy conditions; load the initial combustion chamber temperature parameter and the initial pressure parameter into the flow domain mesh model as the initial field variables to complete the initialization of the flow domain mesh model.

[0054] In step 1023, the mesh refinement process is implemented based on the mass flow rate distribution characteristics in the propellant flow parameters. The specific method is to set the dynamic refinement criterion in the body-fitted mesh module of Fluent software. When it is monitored that the mass flow rate gradient in the injector swirl channel exceeds 50 g / (s·m²), the local mesh refinement operation is triggered, and the mesh size in the trailing edge region of the swirl vane is refined from 0.5 mm to 0.1 mm.

[0055] Meanwhile, a fixed refinement area is set at the inlet section of the catalytic bed. Five layers of boundary layer meshes are generated using the surface shooting mesh technique (the first layer height is 5 μm, and the growth rate is 1.1) to accurately capture the development of the velocity boundary layer near the wall. During the initialization process, the measured initial combustion chamber temperature parameter of 650 K is used as the initial temperature value for the whole field, and a non-uniform temperature field is loaded in the catalytic bed area through the Patch function (the axial temperature gradient is set to 80 K / mm). In addition, the initial pressure parameter of 2.1 MPa is applied in the form of a pressure inlet boundary condition, the pressure function is discretized using a second-order accuracy format, and the convergence residual is set to 1×10^-6. After the initialization is completed, the relative deviations of the whole-field temperature field and pressure field are less than 0.3% and 0.8% respectively, meeting the accuracy requirements of the initial conditions for subsequent transient calculations.

[0056] As an optional but not limited embodiment, the full-cycle numerical simulation of the internal temperature field distribution of the space monopropellant engine in step 103 to generate the dynamic change sequence of the temperature field under the ignition state includes:

[0057] Step 1031: Set the propellant injection time series of the ignition transient process in the flow domain mesh model, and calculate the local flow velocity distribution in the porous media region at each time step based on the pressure fluctuation range in the dynamic boundary conditions.

[0058] In step 1031, the setting of the propellant injection time series is based on the measured time series data. Specifically, the time function of the pressure inlet boundary condition is defined in the transient solver of Fluent software. This function is divided into three stages: 0 - 10 ms is the pressure linear rising stage (slope 0.5 MPa / ms), 10 - 50 ms is the pressure holding stage (constant value 2.1 MPa), and 50 - 500 ms is the pressure adjustment stage (dynamically adjusted according to the PID control algorithm).

[0059] For example, within each time step (1 μs), when solving the local flow velocity distribution in the porous medium region, a coupled implicit algorithm is used to solve the Brinkman-Forchheimer extended function, where the viscous resistance coefficient 1 / α = 1.8×10^8 kg / (m³·s) and the inertial resistance coefficient C2 = 0.34. During the calculation process, the pressure drop across the catalytic bed is monitored in real time. When the pressure drop fluctuation exceeds 5%, the sub-iteration step is automatically triggered, and the maximum number of sub-iterations is set to 20. Another example is that at 20 ms after ignition, the flow velocity at the inlet of the catalytic bed reaches a peak value of 12.3 m / s, and the corresponding local Reynolds number Re = 210 (based on the equivalent pore diameter), and the flow is in the laminar to transitional flow state.

[0060] Step 1032: According to the local flow velocity distribution and the heat release rate function in the chemical reaction model, iteratively solve the time evolution model of the temperature field in the combustion chamber to generate snapshots of the temperature field at continuous timestamps.

[0061] In step 1032, a strongly coupled algorithm is used to solve the temperature field evolution model, specifically by enabling the synchronous solution of the energy function and the species transport function in Fluent. Within each time step, first, the local flow velocity distribution in the porous medium region is calculated, and then the heat release power of the hydrazine decomposition reaction is calculated through the finite rate chemical reaction model.

[0062] Exemplarily, the heat release rate function is implemented through UDF, and its core expression is Q = ΔH·r·As, where: ΔH = -1540 kJ / kg is the enthalpy change of the hydrazine decomposition reaction, r is the surface reaction rate r = 5.6×10^18·k·P N2H4 / (1 + K ads ·P NH3 ), A s = 1.2×10^4 m² / m³ is the specific surface area of the catalytic bed. The k-omega SST turbulence model is selected, the enhanced wall function is used for near-wall treatment, and the radiative heat transfer is calculated through the discrete ordinates model (DO model), and the number of radiative iterations is set to 5 times per time step.

[0063] For example, a dynamic time step control strategy can be adopted during the calculation process. When the temperature change rate in the hot spot area of the catalytic bed exceeds 100 K / ms, the time step is automatically reduced to 0.2 μs. The continuously output temperature field snapshots are stored in the HDF5 format. Each snapshot contains the temperature values of 11.2 million grid nodes, the time stamp interval is 0.1 ms, and the total data volume reaches 2.3 TB.

[0064] Step 1033: Perform spatial interpolation on the temperature field snapshots to generate a dynamic change sequence of the temperature field covering the entire flow domain of the spatial single-component engine; extract the highest temperature peak in the catalytic bed outlet area and the occurrence time of the highest temperature peak in the temperature field dynamic change sequence, and use the highest temperature peak and the occurrence time as the key indicators of the temperature gradient distribution characteristics.

[0065] In step 1033, the cubic spline interpolation algorithm is used for spatial interpolation processing, specifically, full-field data remapping is performed in Tecplot 360 software. Interpolate the temperature field snapshots of the unstructured grid to a uniform Cartesian grid (resolution 0.1mm×0.1mm×0.1mm) to eliminate numerical noise caused by grid distortion. The generation of the temperature field dynamic change sequence is achieved by writing a Python script, which stitches 5000 interpolated temperature fields into a four-dimensional array (X, Y, Z, T) in chronological order. The algorithm for extracting the highest temperature peak in the catalytic bed outlet area is based on the region growing method, and a temperature threshold of 1200K is set. When it is detected that three consecutive grid nodes exceed the threshold, it is determined as a valid peak.

[0066] For example, at 48.7 ms after ignition, the highest temperature peak of 1563K appears at the catalytic bed outlet. The pressure fluctuation frequency corresponding to this moment is 2.4 kHz, and the coincidence degree with the first-order longitudinal vibration mode frequency of the combustion chamber (calculated value 2.38 kHz) reaches 99.2%, indicating that the thermoacoustic coupling effect is the key factor causing the temperature spike. The peak spatial distribution data is converted into STL format through a three-dimensional isosurface extraction algorithm for load input in subsequent thermal stress analysis.

[0067] As an optional but not limited embodiment, in step 104, extracting the thermal stress evolution data of the key components of the spatial single-component engine according to the temperature gradient distribution characteristics at different times in the temperature field dynamic change sequence includes:

[0068] Step 1041: Calculate the instantaneous thermal expansion coefficient of the combustion chamber wall material according to the spatial temperature gradient data in the temperature field dynamic change sequence.

[0069] In step 1041, specifically, the JMatPro material property calculation software is used to establish a temperature-dependent model of the thermal expansion coefficient of GH3128 nickel-based superalloy. By inputting the alloy composition (Cr20.5%, Mo9.2%, W5.1%) and heat treatment process parameters (solution temperature 1180℃ / 2h, aging treatment 800℃ / 8h), an instantaneous thermal expansion coefficient curve in the temperature range of 300 - 1300K is generated.

[0070] For example, the fitting formula for the instantaneous coefficient of thermal expansion curve is: α(T) = 13.8×10^-6 + 0.023×10^-6·T + 1.7×10^-9·T² (unit ), and the correlation coefficient R² = 0.998. During the calculation process, the temperature gradient data of the combustion chamber wall nodes in the dynamic change sequence of the temperature field (such as the axial temperature gradient of the injector head being 1.2×10^5 K / m at 50 ms after ignition) is substituted into this formula, and the instantaneous value of the coefficient of thermal expansion at the corresponding position is . The data conversion is realized through the ANSYS Workbench platform. The field mapping algorithm is used to accurately interpolate the temperature gradient data of the CFD grid to the structural grid nodes, and the interpolation error is controlled within 0.3%.

[0071] Step 1042: Based on the instantaneous coefficient of thermal expansion and the wall structure constraint conditions, derive the thermal stress distribution nephogram at different times.

[0072] In step 1042, the derivation of the thermal stress distribution nephogram is based on the thermoelastic constitutive function. Fixed constraint conditions are set in the ANSYS Mechanical transient structure module: full degree-of-freedom constraints are applied to the combustion chamber flange mounting surface, and axial displacement restrictions are applied to the nozzle outlet end face. The material model selects the bilinear kinematic hardening criterion, and the yield strength of the GH3128 alloy in an 850K environment is defined as 620 MPa, and the tangent modulus is 12 GPa.

[0073] In the embodiment of the present application, the solver uses the direct integration method to calculate the transient thermal stress. The time step is synchronized with the dynamic change sequence of the temperature field (at 0.1 ms intervals), the upper limit of the number of iterations within each time step is set to 25 times, and the convergence residual standard is 1×10^-4. For example, at 120 ms after ignition, the maximum Von Mises stress value of 583 MPa appears at the root of the injector swirler vane, the stress concentration factor in this area reaches 2.3, and the corresponding plastic strain accumulation is 0.15%.

[0074] Step 1043: Generate the historical evolution path of the thermal stress concentration area according to the area coordinates exceeding the material yield strength extracted from the thermal stress distribution nephogram and the duration corresponding to the area coordinates.

[0075] In step 1043, the generation of the historical evolution path of the thermal stress concentration area is realized automatically through a Python script. The specific process is as follows: First, extract the element numbers and their coordinate information exceeding the material yield strength threshold (620 MPa) in the thermal stress distribution nephogram at each time step, and then use the moving time window algorithm (window width 5 ms, step size 0.1 ms) to count the exceeding duration at each spatial position.

[0076] For the injector head region (coordinates X = 32.5 - 37.8 mm, Y = 10.2 - 15.6 mm), the analysis shows that stress over - standard events occur during 48 - 52 ms after ignition, with the maximum continuous over - standard duration of 3.2 ms and the spatial coverage area of 8.7 mm². The historical path data is stored in the form of a 3D point cloud sequence. Each data point contains spatial coordinates, peak stress, and the acting time. By generating a spatio - temporal evolution animation through Paraview software, the trajectory of the stress concentration area migrating from the root of the swirl vane towards the center of the injector can be clearly shown.

[0077] Step 1044: Compare the historical evolution path with a preset overheat failure threshold in time series, and identify the positions of key components at risk of cumulative damage according to the time - series comparison result.

[0078] In step 1044, the time - series comparison analysis uses an improved rain - flow counting method to decompose the historical evolution path of thermal stress into full - cycle and half - cycle events. The preset overheat failure thresholds include a static threshold (80% of the yield strength, i.e., 496 MPa) and a dynamic threshold (equivalent damage amount D = 1 based on Miner's linear cumulative damage theory). For the catalytic bed support ring region (coordinates Z = 85 - 92 mm), it is identified that there are 3 events where the stress peak exceeds the dynamic threshold during 320 - 320.8 ms. The single - peak stress of 625 MPa corresponds to a damage increment of ΔD = 0.0047. When the total cumulative damage reaches 0.82, the system automatically marks this region as a high - risk part, marks it with red contour lines in the 3D model, and generates a JSON - format report containing position coordinates, over - standard times, and damage degrees.

[0079] As an optional but not limited embodiment, the output of the performance stability evaluation result in the ignition state in step 104 includes:

[0080] Step 1045: Calculate the thermal fatigue cycle times of the positions of the key components according to the historical evolution path of the thermal stress concentration region.

[0081] In step 1045, the calculation of the thermal fatigue cycle times is based on the Coffin - Manson - Basquin function, and the Coffin - Manson - Basquin function is as follows:

[0082] N f =(Δε pl / ε f' ) 1 / c ·(Δσ / σ f' ) 1 / b , where Δε pl = 0.12% is the plastic strain amplitude, ε f'=0.35 is the fatigue ductility coefficient, c=-0.6 is the fatigue ductility index, Δσ=583MPa is the stress amplitude, σ f' =980MPa is the fatigue strength coefficient, b=-0.09 is the fatigue strength index. Substituting the parameters of the maximum stress concentration area of the injector head into the calculation, the equivalent fatigue cycle number N of a single ignition process is obtained. f =1.8 times / cycle. Combined with the 10,000 ignition missions required for the engine design life, the cumulative damage factor D=Σ(n i / N fi )=10000 / 1.8=5555.6, which far exceeds the allowable value D=1, indicating that there is a serious risk of fatigue failure in this part.

[0083] Step 1046: Based on the mapping relationship between the number of thermal fatigue cycles and the material durability database, predict the remaining life range of the space monopropellant engine under continuous ignition conditions.

[0084] In step 1046, the remaining life prediction is performed using a probabilistic statistical approach, using low-cycle fatigue test data for GH3128 alloy from the material durability database (sample size n = 120, 95% confidence level). A Weibull distribution fit yields a shape parameter of β = 2.1 and a scale parameter of η = 2350 cycles. Substituting the calculated equivalent number of cycles (1.8 cycles / ignition) into the reliability function R(t) = exp[-(N / η)^β], the resulting reliability after 1000 ignitions is 67.3%, with a remaining life interval of [832, 1215] cycles (calculated using the bootstrap method). This result is written to the evaluation system via a MySQL database interface and updated in conjunction with real-time monitoring data.

[0085] Step 1047: Combining the remaining life span with the key indicators of the temperature gradient distribution characteristics, a stability assessment report including overheating risk probability and performance degradation trend is generated.

[0086] In step 1047, the stability assessment report is generated by integrating data from multiple sources. First, a spatiotemporal correlation analysis is performed between the historical evolution path data of the thermal stress concentration area and the highest peak temperature (1563K) at the catalyst bed outlet in the temperature field dynamic change sequence. It was found that when the temperature peak exceeded 1500K, the probability of stress exceeding the standard within the subsequent 5ms increased to 78%. Using the Monte Carlo method to simulate 1000 ignition cycles, the probability of thermal fatigue cracking of the injector head after 500 ignitions was calculated to be 92.7%. The performance degradation trend was characterized by a specific impulse decrease of 0.8% per 100 cycles. The report is output in PDF format and includes a three-dimensional thermomechanical coupling cloud map, SN curve fitting results, and a histogram of the remaining life distribution at key locations.

[0087] Step 1048: Mark the results exceeding the preset risk level in the stability assessment report as the structural weak points to be optimized.

[0088] In Step 1048, the structural weak point marking algorithm is evaluated based on the Risk Priority Number (RPN = Severity × Occurrence × Detection). Set the severity level S = 9 (catastrophic failure), the occurrence O = 7 (3.2 times per 100 ignitions), and the detection D = 3 (existing sensor coverage rate 85%). Calculate that the RPN at the root of the injector swirl vane is 189, exceeding the preset threshold of 150. The system automatically marks this area in the CAD model and pushes optimization suggestions: add cooling holes with a diameter of 1.5 mm at coordinates X = 35.2 mm and Y = 12.8 mm, and arrange them in a circular array with a hole pitch of 2 mm. After optimization and simulation verification, it shows that the highest temperature in this area is reduced by 127 K, the peak thermal stress drops to 498 MPa, and the RPN value drops to 63, meeting the safety design requirements.

[0089] It can be seen that the above embodiments realize the full-life cycle performance monitoring of the ignition process of the space monopropellant engine by establishing a thermal-mechanical-damage multi-field coupling evaluation system. Especially, the thermal stress evolution analysis based on the dynamic change sequence of the real-time temperature field can accurately quantify the cumulative damage degree of key components, providing data support for formulating condition-based maintenance strategies and optimizing the cooling system design. The dynamic threshold adjustment algorithm and the probabilistic life prediction model adopted in the above embodiments can reduce the uncertainty of engine reliability assessment, thus effectively solving the problems of over-conservatism or insufficiency in traditional methods.

[0090] As an optional but not limited embodiment, the three-dimensional discretization modeling of the pore structure of the catalytic bed of the space monopropellant engine in Step 1021 based on the porous medium model to generate the topological data of the catalytic bed pore distribution includes:

[0091] Step 1021: Obtain the set of geometric parameters of the catalytic bed, and the set of geometric parameters includes the pore diameter distribution range, the curvature of the pore connection path, and the pore density gradient.

[0092] For example, a high-resolution three-dimensional imaging of the iridium / aluminum oxide catalyst bed can be carried out by using a microfocus X-ray computed tomography system (model ZEISS Xradia 520Versa). The scanning parameters are set as voltage 80 kV, current 88 μA, and the spatial resolution reaches 2 μm, and a tomographic image sequence of 1024×1024×2048 voxels is obtained.

[0093] In addition, image processing can be performed using Avizo Fire software. The region growing algorithm is adopted to segment pores and the solid skeleton. The calculated pore diameter distribution ranges from 15 to 85 μm (peak diameter 38 μm), and the statistical mean of the curvature of the pore connection path is , and the pore density gradient linearly increases from 0.32 at the inlet end to 0.41 at the outlet end along the axial direction.

[0094] For another example, the geometric parameter set is converted into a JSON format input file through a Python script, which contains the equivalent diameter, centroid coordinates of each pore unit, and the adjacent pore connection relationship, serving as the basic data source for subsequent mesh generation.

[0095] Step 1022: Based on the pore diameter distribution range, identify the multi-scale pore morphological characteristic parameters inside the catalyst bed. The multi-scale pore morphological characteristic parameters include the axial direction of the main pore channel, the intersection angle of the branched pores, and the pore wall roughness.

[0096] In step 1022, the identification of the multi-scale pore morphological characteristic parameters is achieved through a customized algorithm: First, set a diameter threshold of 45 μm. Pores larger than this value are defined as the main pore channels. The skeleton extraction algorithm is used to generate the centerline trajectory of the main pore channels, and the standard deviation of the curvature change rate is calculated as , and the maximum deflection angle of the extension direction is 23.5°.

[0097] The detection of the secondary branched pores can be based on the Euclidean distance transformation. Pores with a distance less than 10 μm from the main pore channels are identified as connecting branches, and the measured intersection angle distribution ranges from 28° to 152° (mean value 67°). The quantification of the pore wall roughness uses a three-dimensional surface topography analysis module. Gaussian filtering (cutoff wavelength 5 μm) is performed on the point cloud data of the segmented pore walls, and the surface roughness Ra = 1.8 μm is calculated. The generated two-dimensional Fourier spectrum shows that the main spatial frequency components are concentrated in the interval.

[0098] Step 1023: Divide non-uniform grid units according to the axial direction of the main pore channel and the pore density gradient. The size of the non-uniform grid units is adaptively adjusted along the curvature of the pore connection path, and the grid nodes are aligned with the intersection angles of the branched pores.

[0099] In step 1023, the division of the non-uniform grid units is performed by the mesh generator of COMSOL Multiphysics, and the dynamic adjustment rules are set: Along the axial direction of the main pore channel, when the curvature change rate exceeds , the grid unit length is reduced from 50 μm to 20 μm; in the region where the intersection angle of the branched pores is less than 60°, encrypted nodes are inserted to reduce the node spacing to 10 μm.

[0100] For example, the densely porous region (density > 0.38) is filled with tetrahedral meshes, and the element size is controlled within 15 - 30 μm; the sparse region (density < 0.35) uses hexahedral meshes, and the element size is extended to 50 - 80 μm. The generation of the wall boundary layer mesh is based on the roughness spatial distribution map. In the region where Ra > 2 μm, 3 layers of prism layer meshes are set (the height of the first layer is 2 μm, and the growth rate is 1.2), and the roughness amplitude is input as a wall function correction term into the turbulence model.

[0101] Step 1024: Perform topological consistency verification on the non-uniform mesh elements, judge the continuity of the pore wall roughness data of adjacent mesh elements and the fracture state of the pore connection path, and obtain the topological consistency verification result.

[0102] In step 1024, the topological consistency verification process includes multi-level detections: First, traverse adjacent mesh elements through a Python script, calculate the root mean square deviation of the roughness data at the junction, and mark it as a non-conforming element when the deviation exceeds the preset tolerance of 0.5 μm (for example, the deviation reaches 0.73 μm at coordinates X = 125 μm, Y = 63 μm).

[0103] Furthermore, the continuity verification of the main pore channels adopts the centerline curvature difference method to detect abnormal segments where the difference in the curvature change rate between adjacent elements exceeds (for example, there are 3 mutations in the axial position range of Z = 520 - 535 μm). The alignment detection of the branch pore intersection angles is carried out by calculating the node position offset, and it is found that the node offset reaches 7.3 μm at coordinates X = 88 μm, Y = 214 μm. The finally generated verification report contains the positions of 17 non-conforming elements, 23 curvature mutation points, and 45 offset nodes, all recorded in the form of three-dimensional coordinates.

[0104] Step 1025: According to the topological consistency verification result, perform local reconstruction on the mesh elements that fail the verification, adjust the mesh node positions to match the branch pore intersection angles, and generate an optimized porous medium mesh model.

[0105] In step 1025, the local reconstruction operation is implemented for the mesh elements that fail the verification: for the elements with roughness mutations, use the Laplacian smoothing algorithm to adjust the node positions to reduce the roughness gradient at the junction to within 0.4 μm; insert transition mesh elements in the regions with curvature mutations in the main pore channels to control the difference in the curvature change rate within below; optimize the displacement of the offset nodes of the branch pores by the least squares method, and the node movement amount does not exceed 5 μm. In the reconstructed porous medium mesh model, the number of non-conforming elements is reduced from 85 initially to 3, the node alignment accuracy is improved to 0.8 μm, and the curvature continuity compliance rate is increased from 78% to 99.6%.

[0106] Step 1026: Extract the pore morphology characteristic parameters, grid size, and node coordinates of each grid cell in the porous medium grid model, and generate catalytic bed pore distribution topology data including pore space distribution, connectivity characteristics, and morphological attributes.

[0107] In step 1026, the generation of the catalytic bed pore distribution topology data is completed through the attribute extraction module of Paraview software. Traverse each grid cell and record the pore morphology characteristic parameters: the axis direction vector of the main pore channel (unit vector components Vx = 0.82, Vy = 0.12, Vz = 0.56), the intersection angle of the branched pores (67° ± 12°), and the Ra value of the pore wall roughness (1.8 ± 0.3 μm).

[0108] In the embodiment of the present application, the grid size data is statistically analyzed by region. The average size of the main channel area is 28 μm, the branched area is 15 μm, and the boundary layer area is 5 μm. The node coordinate data is stored in floating-point format with an accuracy of 0.1 μm. The final data set is encapsulated in the HDF5 format, containing 1.2×10^6 records, occupying a storage space of 23 GB, and establishing an index association through the MySQL database to the subsequent simulation module.

[0109] As an optional but not limited embodiment, the identification of the multi-scale pore morphology characteristic parameters inside the catalytic bed in step 1022 includes step 10220: Based on the pore diameter distribution range, divide the pores in the catalytic bed into main pore channels and secondary branched pores. The diameter of the main pore channel is greater than the preset threshold and extends along the axis; extract the centerline trajectory of the main pore channel, calculate the curvature change rate and the deflection angle of the extension direction of the centerline trajectory, and use them as the core parameters of the axis direction of the main pore channel; in the secondary branched pores, detect the inlet position connected to the main pore channel, measure the angle between the branched pore at the inlet and the main pore channel, and generate a set of intersection angles of the branched pores; quantify the microscopic undulation characteristics of the pore wall through surface scanning data, and generate a spatial distribution map of the pore wall roughness.

[0110] In the embodiment of the present application, the division of the main pore channel and the secondary branched pores is realized by the adaptive threshold method: Perform Gaussian fitting on the pore diameter distribution histogram, determine 45 μm as the valley value point of the bimodal distribution, and classify 9532 pores with a diameter ≥ 45 μm as the main channel. The extraction of the centerline trajectory adopts a thinning algorithm to iteratively remove boundary voxels, generate a smooth B-spline curve, and calculate the standard deviation of the curvature change rate along the axial distribution as , and the maximum deflection angle appears at Z = 325 μm (23.5°).

[0111] For another example, the secondary branch inlet detection is based on morphological erosion operation, identifying 18,920 connection points, and measuring the average angle between them and the main channel to be 67° with a standard deviation of 14°. The wall roughness analysis uses wavelet transform to decompose the surface topography, extracts the detail component at scale 3 (wavelength 4 - 8 μm) as the roughness characterization quantity, and generates a two-dimensional distribution cloud map showing that the roughness of the inlet region is 18% higher than that of the outlet region.

[0112] As an optional but not limited embodiment, the dividing of the non-uniform grid cells in step 1023 includes step 10230: along the axis direction of the main pore channel, dynamically adjust the length of the grid cells according to the curvature change rate, where the curvature change rate is negatively correlated with the length of the grid cells; insert encrypted grid nodes at the intersection angles of the branch pores so that the node density is inversely proportional to the intersection angle of the branch pores; based on the pore density gradient, set tetrahedral grid cells in the pore-dense area and hexahedral grid cells in the pore-sparse area; map the spatial distribution map of the wall roughness to the surface of the grid cells to generate a wall boundary layer grid with roughness attributes.

[0113] In the above embodiment, the specific operations of non-uniform grid division include: in the section where the curvature change rate of the main pore channel exceeds (such as Z = 120 - 135 μm), gradually reduce the length of the grid cells from 50 μm to 20 μm; in the narrow area with a branch intersection angle of 28° (X = 205 μm, Y = 178 μm), set the node density to be increased to 150 nodes / mm²; generate tetrahedral grids with a minimum size of 15 μm in the pore-dense area (density > 0.38) and use 50 μm hexahedral grids in the sparse area. The wall roughness mapping is realized through an interpolation algorithm, assign the Ra value as an additional attribute of the grid surface nodes, and read and correct the wall function parameters through UDF in Fluent.

[0114] As an optional but not limited embodiment, the topological consistency verification of the non-uniform grid cells in step 1024 further includes step 10240: judge whether there is a sudden change in the wall roughness data of adjacent grid cells at the junction, and if the mutation amplitude exceeds the preset tolerance, mark it as an inconsistent cell; verify the continuity of the grid cells in the main pore channel, where the grid cell continuity characterizes that the curvature change rate of the centerline trajectory smoothly transitions between adjacent cells; detect whether the grid nodes at the intersection angle of the branch pores are aligned, and if there is an offset, calculate the node position correction amount; generate a topological consistency verification report including the positions of inconsistent cells, the mutation amplitude, and the node offset.

[0115] For example, the detailed process of topological consistency verification may include: using the differential method to calculate the gradient of roughness data at the intersection of adjacent mesh cells, triggering an alarm when the absolute value of the gradient exceeds 0.5μm / μm (for example, the gradient reaches 0.63μm / μm at X=152μm and Y=93μm); main channel continuity detection by calculating the angle between the tangent vectors of adjacent cell centerlines, and determining discontinuity if it exceeds 5° (for example, the angle is 7.2° at Z=430μm); branch node alignment verification using a nearest neighbor search algorithm with an offset tolerance of 2μm, detecting 45 nodes that exceed the standard. The verification report is output in XML format, containing the inconsistent cell ID, spatial coordinates, and quantized deviation values, for use by the mesh optimization module.

[0116] As an optional but not limiting embodiment, the method further includes:

[0117] Step 201: Collecting measured temperature data of the space monopropellant engine in historical ignition tasks and performance degradation records corresponding to the measured temperature data.

[0118] For example, the measured temperature data of historical ignition missions can be collected by a platinum resistance temperature sensor array (model PT100, measurement point spacing 5 mm) installed on the wall of the combustion chamber of a space monopropellant engine. The array has accumulated 1.2×10^6 sets of temperature data in 50 ignition missions, covering a temperature range of 300K to 1350K and a sampling frequency of 100kHz.

[0119] In addition, a piezoelectric dynamic pressure sensor (model Kistler 603B) was used to record the combustion chamber pressure fluctuation data. The pressure measurement range was 0-10 MPa and the linearity error was ±0.3% FS. The performance degradation record included the activity decay rate of the iridium catalyst in the catalytic bed (measured by X-ray photoelectron spectroscopy). The data also includes data on the reduction of the content of GH3128 alloy from 78% to 62%, crack growth length of the GH3128 alloy in the injector head (maximum crack length measured by metallographic microscopy was 0.82 mm), and specific impulse reduction data (vacuum specific impulse decreased from 230s to 225s). All data is stored in HDF5 format, and a relational database containing timestamps, spatial coordinates, temperature values, and corresponding performance indicators is established.

[0120] Step 202: Using the measured temperature data, perform error correction on the temperature field dynamic change sequence to generate a calibrated temperature prediction model.

[0121] For example, the error correction process is implemented by coupling the Kalman filter algorithm with the computational fluid dynamics model. First, the historical measured temperature data is spatially aligned with the dynamic change sequence of the temperature field at the corresponding moment, and the residual distribution field is calculated through the Curve Fitting Toolbox of MATLAB. It is found that there is a systematic deviation in the outlet area of the catalyst bed (the maximum deviation value is 87K).

[0122] For another example, a physics-based correction model can be established to correct the equivalent thermal conductivity (adjusted from 12 W / m·K to 15.3 W / m·K) and the surface reaction activation energy (corrected from 67.8 kJ / mol to 71.2 kJ / mol) in the porous medium region in the Fluent software. The calibrated temperature prediction model shows in an exemplary validation set (10 ignition tasks not involved in training) that the temperature prediction error in the hot spot area of the catalyst bed is reduced from ±8.7% to ±2.3%, and the prediction accuracy of the temperature gradient at the root of the injector swirl vane is improved to 97.6%.

[0123] Step 203: Jointly train the calibrated temperature prediction model with the thermal stress evolution data to generate an abnormal temperature warning model for real-time monitoring.

[0124] In step 203, the construction of the abnormal temperature warning model is implemented based on the PyTorch framework. The input layer receives a three-dimensional temperature field sequence of 512×512×500 (spatial resolution 0.1 mm, time resolution 0.1 ms), and the output layer generates a local difference coefficient matrix. The convolutional neural network model is improved using the ResNet-18 architecture, adding a three-dimensional convolutional layer (kernel size 3×3×3, stride 1) to extract spatio-temporal features, and the output dimension of the fully connected layer matches the number of grid nodes (1.2×10^6 dimensions).

[0125] In the embodiment of the present application, the training data includes the corrected dynamic change sequence of the temperature field and the corresponding residual distribution of 45 ignition tasks, the batch size is 32, the learning rate is 0.001, and the loss function of the validation set converges to 0.023 after 500 training cycles. In the deployment stage, the trained model weights are integrated into the real-time monitoring system, key attention channels are set in the Catalyst bed area, and the moving average value of the local difference coefficient (window width 5 ms) is calculated. When the coefficient exceeds the safety tolerance of 0.15, an alarm is triggered.

[0126] Step 204: In the ignition task to be monitored, dynamically adjust the propellant flow rate based on the deviation signal output by the abnormal temperature warning model to suppress the overheating risk.

[0127] In step 204, dynamic regulation of the propellant flow rate is achieved through a proportional-integral-differential controller. Specifically, when the abnormal temperature warning model detects that the local difference coefficient in the injector head area (coordinates X = 35.2 mm, Y = 12.8 mm) reaches 0.18, the control algorithm calculates the required flow correction amount Δm = 1.3 g / s based on the spatial weight distribution of this area (influence factor 0.76) and the current propellant mass flow rate of 12.5 g / s.

[0128] For example, an adjustment command can be sent via Modbus TCP to the propellant supply system's electric control valve (Fisher Vee-Ball V300) to increase the valve opening of injection unit 3 from 45% to 58%, increasing the propellant mass flow rate in that area to 14.1 g / s. After the adjustment, monitoring data showed that the temperature in the target area dropped from 935K to 863K within 120ms, and the local coefficient of variation dropped to 0.09. The warning status was automatically lifted, and the adjustment log was recorded.

[0129] As an optional but not limiting embodiment, the step 203 of generating an abnormal temperature warning model for real-time monitoring includes:

[0130] Step 2031: extracting the residual distribution between the temperature field dynamic change sequence and the measured temperature data from the historical ignition tasks.

[0131] In step 2031, the residual distribution is extracted using a spatial registration algorithm to interpolate the measured temperature data of the historical ignition mission to the grid nodes of the computational fluid dynamics model and calculate the absolute deviation ΔT = |T sim -T exp Statistics show that the residual peak at the catalyst bed outlet reaches 92 kHz between 50 and 60 milliseconds after ignition. This deviation exhibits axial propagation (at a velocity of 1.8 m / s) and is strongly correlated with the combustion chamber pressure oscillation frequency of 2.4 kHz. The residual distribution database contains 500 transient fields, each recording the deviation values of 1.2 × 10^6 nodes, for a total data volume of 3.4 TB.

[0132] Step 2032: Use the residual distribution to train a convolutional neural network model to generate a spatiotemporal feature extractor for temperature prediction residuals.

[0133] In step 2032, the convolutional neural network model is trained using a distributed computing architecture, utilizing four NVIDIA A100 GPUs for parallel processing. The input data is a 512×512×500 three-dimensional temperature field sequence. After processing through five layers of three-dimensional convolution (channels 64-512) and three layers of full connection, the output is a residual prediction field of the same dimensions. The loss function is defined as the weighted mean square error: , where the weight w i is set to 3.0 in the catalytic bed area and 1.0 in other areas. After training, the model achieves an accuracy with an average absolute error of 4.3K and a peak error of 13.7K on the test set.

[0134] Step 2033: Embed the spatio-temporal feature extractor into the abnormal temperature warning model, and calculate the local difference coefficient between the current temperature field and the predicted value in real time.

[0135] In step 2033, the embedding of the spatio-temporal feature extractor is realized through a dynamic link library. The trained neural network model is converted into the ONNX format and integrated into the data processing module of the real-time monitoring system. During online calculation, the latest temperature field data is received every 0.1 ms, and a 512-dimensional feature vector is output through the feature extractor, and then mapped into a local difference coefficient matrix through a fully connected layer. The coefficient calculation formula is C = Σ(f i ·w i ) / Σ|w i |, where f i is the activation value of the i-th feature channel, and w i is the trained attention weight. The system sets a double threshold: when C > 0.15, a yellow warning is triggered; when C > 0.25, a red warning is triggered and emergency adjustment is started.

[0136] Step 2034: When the local difference coefficient exceeds the preset safety tolerance, activate the propellant flow regulation instruction to balance the combustion chamber heat load.

[0137] In step 2034, the generation of the propellant flow regulation instruction is based on a fuzzy control rule base. The input variables are defined as the local difference coefficient C and its change rate dC / dt, and the output variable is the flow correction amount Δm. The rule base contains 21 control strategies. For example: if C ∈ [0.15, 0.2) and dC / dt > 0.1 / ms, then Δm = 1.2 g / s; if C ≥ 0.25, then Δm = 2.5 g / s and auxiliary cooling is started. The instruction transmission delay is controlled within 0.8 ms to ensure intervention at the initial stage of abnormal temperature development.

[0138] As an optional but not limited embodiment, in step 204, dynamically adjusting the propellant flow based on the deviation signal output by the abnormal temperature warning model includes:

[0139] Step 2041: Combine the spatial distribution of the local difference coefficient and the deviation signal to identify the target area where the temperature abnormally rises.

[0140] For example, for target area recognition, the connected component analysis method can be adopted to cluster the grid nodes with a local difference coefficient exceeding 0.15, and a minimum clustering volume of 5 mm³ is set. For example, 320 ms after ignition, an abnormal high-temperature cluster with a diameter of 3.2 mm is detected in the catalytic bed support ring area (coordinate Z = 85 - 92 mm), which contains 126 non-compliant nodes and an average difference coefficient of 0.21. The spatial distribution analysis shows that there is a direct fluid path association between this area and the 3rd and 7th injection units of the propellant supply pipeline.

[0141] Step 2042: Based on the topological relationship between the position information of the target area and the propellant injection pipeline, calculate the correction amount of the flow valve opening to be adjusted.

[0142] For example, the calculation of the correction amount of the flow valve opening is based on the fluid network model, and a pressure-flow characteristic function of the propellant pipeline is established: , where K v is the valve flow coefficient, and ρ = 1.004 g / cm³ is the density of hydrazine propellant. For the 3rd injection unit, the current opening of 45% corresponds to K v = 0.62, and it needs to be adjusted to an opening of 58% to make K v = 0.89, so as to increase the mass flow rate from 12.5 g / s to 14.1 g / s. The correction amount calculation module synchronously considers the pipeline transmission delay (3.2 ms) and the valve response time (8 ms), and generates a feedforward compensation instruction to act 1.5 cycles in advance.

[0143] Step 2043: Send the correction amount of the flow valve opening to the engine control system, so as to rematch the propellant distribution and the current heat load demand through the engine control system to obtain an adjusted temperature field dynamic change sequence; monitor and verify the adjusted temperature field dynamic change sequence, and complete the monitoring and verification when the local difference coefficient drops back within the safe tolerance range.

[0144] For example, the response verification of the engine control system is executed through a hardware-in-the-loop test platform. The adjusted flow parameters are input into the real-time simulation model to monitor the evolution of the temperature field dynamic change sequence. For example, after the valve opening is adjusted, the temperature gradient in the target area drops from 1.1×10^5 K / m to 7.2×10^4 K / m, and the peak thermal stress drops from 625 MPa to 538 MPa. The system continuously monitors the local difference coefficient within a 50 ms cycle. When the coefficients of three consecutive sampling points are lower than 0.12, it is determined that the adjustment is effective; otherwise, a secondary adjustment is triggered. The obtained monitoring and verification data can be written into the SQL database to provide training samples for control strategy optimization.

[0145] It can be seen that the above embodiments realize the real-time suppression of the overheating risk of the space single-component engine by constructing a data-model dual-driven closed-loop control system. The residual-driven adjustment mechanism adopted in the above embodiments can shorten the response time of temperature anomalies, thereby effectively preventing the occurrence of potential overheating faults, further extending the life of key components, reducing the performance decay rate, and providing technical guarantee for the long-term reliable operation of the spacecraft in orbit.

[0146] As an optional but not limited embodiment, the method further includes:

[0147] Step 301: According to the performance stability evaluation result, screen out the failure modes that repeatedly appear in the historical evolution path of the thermal stress concentration area.

[0148] In step 301, specifically in implementation, by analyzing the performance stability evaluation reports of 300 historical ignition missions, a pattern recognition algorithm is used to extract the common features in the historical evolution path of the thermal stress concentration area. The K-means clustering algorithm is used to group the spatial coordinates of 45 parts marked as high-risk parts, and three main failure modes are identified: the first type is the periodic thermal stress exceeding the standard at the root of the injector swirl vane (appearance frequency 78%, single-time exceeding standard duration 3.2±0.8 ms), the second type is the transient temperature peak of the catalytic bed support ring (peak temperature 935K, occurrence probability 62%), and the third type is the cumulative plastic strain at the combustion chamber flange connection (maximum strain 0.23%, growth rate 0.004% / cycle).

[0149] For example, the construction of the failure mode database is based on MySQL relational data tables. Each record contains spatial coordinates, types of parameters exceeding the standard, occurrence time and duration. Through correlation analysis, it is found that the second failure mode has a strong correlation (Pearson coefficient 0.76) with the fluctuation of the propellant mass flow rate (standard deviation>0.8 g / s).

[0150] Step 302: Generate an objective iterative function for the structural optimization of the space single-component engine based on the failure mode. The objective iterative function includes a weight for reducing the peak thermal stress and a material mass constraint condition.

[0151] In this embodiment, the construction of the objective iterative function can be based on multi-objective optimization theory. The decision variables are set as the diameter (range 0.5-2.0 mm), spacing (2-8 mm) and inclination angle (0°-30°) of the additional support structure of the catalytic bed. The objective function includes two optimization terms: the weight term f1 for reducing the peak thermal stress = Σ(σ vm / σ yield)², the material quality constraint term f2 = (M - M0) / M0 (M0 is the initial mass). The constraint conditions are set as follows: the total mass increment does not exceed 5%, and the maximum stress of the support structure does not exceed 80% of the allowable value of the material. Through the parametric design module of ANSYS Workbench, the decision variables are associated with the geometric parameters of the flow domain grid model to generate an initial sample space containing 128 design points.

[0152] Step 303: Use the genetic algorithm to solve the target iterative function for multiple generations to generate multiple candidate structural improvement solutions that meet the durability requirements; perform virtual ignition tests on the multiple candidate structural improvement solutions through the flow domain grid model, and select the target structural improvement solution with a thermal stress distribution that meets the set improvement conditions as the optimization strategy according to the virtual ignition test results.

[0153] For example, the implementation of the genetic algorithm can adopt the NSGA-II multi-objective optimization framework, with a population size of 64, a crossover probability of 0.85, and a mutation probability of 0.15. Each generation of evolution includes the following operations: First, calculate the thermal stress distribution of each individual through ANSYS Mechanical, and extract the maximum Von Mises stress value at the injector head; then call Fluent for transient fluid-thermal coupling simulation to obtain the adjusted dynamic change sequence of the temperature field; finally, perform non-dominated sorting and crowding degree calculation according to the objective function values.

[0154] For another example, after 50 generations of evolution, the Pareto front converges, and 3 candidate solutions are screened out: Solution A (support structure diameter 1.2 mm, spacing 4 mm, inclination angle 15°) reduces the thermal stress peak by 37.5% and increases the mass by 3.8%; Solution B (diameter 1.5 mm, spacing 5 mm, inclination angle 22°) reduces the stress by 28.7% and increases the mass by 2.1%; Solution C (diameter 0.8 mm, spacing 3 mm, inclination angle 8°) reduces the stress by 41.2% and increases the mass by 4.9%. The virtual ignition test shows that the temperature peak at the catalytic bed support ring of Solution C drops from 935 K to 862 K, and the thermal stress peak drops from 625 MPa to 498 MPa, and it is selected as the final optimization strategy.

[0155] As an optional but not limited embodiment, the generation of multiple candidate structural improvement solutions meeting the durability requirements in step 303 includes step 3030: inserting geometric parameter variables of an additional support structure into the topological data of the pore distribution in the catalytic bed; adjusting the local flow velocity distribution in the porous medium region according to the geometric parameter variables of the additional support structure; recalculating the dynamically changing sequence of the adjusted temperature field, and determining whether the reduction amplitude of the peak thermal stress reaches a preset threshold; if the reduction amplitude of the peak thermal stress reaches the preset threshold, adding the structural improvement features corresponding to the current geometric parameter variables to the candidate solution set, otherwise continuing iterative optimization until the termination condition is met.

[0156] In step 3030, the parametric modeling of the additional support structure is implemented through SpaceClaim Direct Modeler, and cylindrical supports are inserted into the low-porosity region (pore density < 0.32) of the topological data of the pore distribution in the catalytic bed. The variable of the support diameter is controlled by an APDL script and is distributed in an arithmetic sequence along the axial direction (spacing variable Δ = 3 - 8 mm). The adjustment of the local flow velocity distribution in the porous medium region is carried out by pore-scale simulation, and the change in permeability (decreasing from 3.2×10^-10 m² to 2.7×10^-10 m²) after the introduction of the supports is calculated in COMSOL Multiphysics, corresponding to an 18% increase in the pressure drop.

[0157] For another example, the evaluation of the reduction amplitude of the peak thermal stress is based on the average value of 10 virtual ignition tests. When the reduction exceeds 30% (for example, solution C reaches 41.2%), the current geometric parameter combination is stored in the candidate solution library, otherwise the combination of the support inclination angle and diameter is adjusted by the orthogonal test method and continues to iterate until the maximum iteration number of 100 times or the reduction standard is met.

[0158] As an extensible embodiment, the coupling of the topological data of the pore distribution in the catalytic bed with the geometric structure of the combustion chamber in step 1022 by combining the propellant decomposition reaction rate sub-model in the chemical reaction model to generate an initial flow domain grid for multi-physics field coupling includes:

[0159] Step 10221: Extract the pore space distribution characteristics, connectivity characteristics, and morphological attributes from the topological data of the pore distribution in the catalytic bed. The morphological attributes include the roughness level of the pore wall surface and the intersection angle of the branched pores; based on the local reaction rate distribution function in the propellant decomposition reaction rate sub-model, identify the pore regions in the pore space distribution characteristics of the catalytic bed that are positively correlated with the propellant decomposition reaction rate, and generate a set of spatial coordinates of the reaction rate sensitive region.

[0160] For example, the identification of the reaction rate sensitive region is based on the coupled analysis of the pore space distribution characteristics of the catalytic bed and the surface reaction kinetics. First, the activity distribution of the iridium catalyst is determined by X-ray photoelectron spectroscopy (XPS), and it is found that the density of active sites reaches 5.8×10^18 sites / m² in the region where the surface roughness Ra of the pore wall is <1.5 μm (23% higher than that in the relatively rough region). An image processing algorithm is used to extract the connected pore paths with diameters of 35 - 50 μm (accounting for 68% of the total path length), calculate the angle between the extension direction of these paths and the fluid flow direction at the inlet of the flow channel, and eliminate the paths with an angle >45° (accounting for 12%). The finally generated target reaction rate sensitive region contains 327 main channels with a total length of 1.2 m, and the set of spatial coordinates is imported into the chemical reaction model of ANSYS Fluent in CSV format.

[0161] Step 10222: Align the set of spatial coordinates of the reaction rate sensitive region with the position of the flow channel inlet in the combustion chamber geometry. After the spatial alignment is completed, the fluid transmission path at the outlet of the pore channels in the reaction rate sensitive region coincides with that at the flow channel inlet; based on the heat release rate gradient in the propellant decomposition reaction rate sub-model, set the reaction rate weight coefficient based on the flow channel wall of the combustion chamber geometry, and the reaction rate weight coefficient decreases as the distance between the flow channel wall and the reaction rate sensitive region increases.

[0162] Exemplarily, the spatial alignment operation is realized through a coordinate transformation matrix, which transforms the local coordinate system of the reaction rate sensitive region (with the origin at the center of the catalytic bed inlet) to the global coordinate system of the combustion chamber. Boolean operations are performed in SpaceClaim software to ensure that the geometric center deviation between the outlet of the pore channels in the sensitive region and the flow channel inlet is less than 0.1 mm. The setting of the reaction rate weight coefficient is based on the inverse distance weighting method, and the weight W is defined as 1 / (1 + 0.2d), where d is the Euclidean distance (in mm) between the node on the flow channel wall and the nearest sensitive region. In the region where the weight coefficient W > 0.7 (d < 1.5 mm), the high-precision chemical reaction calculation mode is activated, and the time step is reduced to 0.5 μs.

[0163] Step 10223: Based on the reaction rate weight coefficient, adjust the local density of the grid in the combustion chamber flow channel, insert a dense grid layer in the region where the weight coefficient is higher than the preset threshold, and generate an initial flow channel grid that matches the pore structure; map the pore connectivity characteristics in the topological data of the catalytic bed pore distribution to the initial flow channel grid to generate a fluid transmission interface between the pore channels and the flow channel grid, and the grid nodes of the interface correspond one-to-one with the grid nodes at the outlet of the pore channels.

[0164] In step 10223, local encryption of the flow channel mesh uses the size function of ANSYS Meshing. The maximum element size is set to 0.1 mm in the area where W > 0.7, and the number of boundary layer mesh layers is increased to 5 layers (the height of the first layer is 2 μm). The interface between the pore channel and the flow channel mesh is generated by the node mapping algorithm, and 1.2×10^5 pairs of bonded contact nodes are created on the interface surface to ensure that the mass flow transfer error is less than 0.3%. The grid model verification shows that under the condition of the propellant mass flow rate of 12.5 g / s, the relative deviation between the total flow rate at the outlet of the pore channel and the received flow rate at the inlet of the flow channel is 0.18%, meeting the requirement of mass conservation.

[0165] Step 10224: According to the temperature dependence in the propellant decomposition reaction rate sub-model, correlate the local flow velocity distribution in the pore channel with the heat conduction rate in the flow channel mesh at the fluid transfer interface to generate a multi-physics data chain that couples the reaction rate and fluid flow heat transfer; perform the coupling verification of the multi-physics data chain to detect whether the propellant mass flow rate at the outlet of the pore channel and the received flow rate at the inlet of the flow channel are conserved. If not, readjust the corresponding relationship of the grid nodes at the fluid transfer interface; according to the coupling verification results, integrate the initial flow channel grid, pore structure grid, and multi-physics data chain to generate the initial flow domain grid for the multi-physics coupling.

[0166] Exemplarily, the coupling of the multi-physics data chain is achieved through the coupling solver of Fluent. Bidirectional data exchange is set at the fluid transfer interface: the flow velocity distribution (U, V, W components) in the pore channel is transferred to the flow channel mesh as the inlet boundary condition every 0.1 ms, and at the same time, the flow channel wall temperature is fed back to the pore model as the thermal boundary condition for the chemical reaction. The coupling verification uses transient mass conservation monitoring. The cumulative mass error within a 500 ms simulation period is 0.023 g, and the relative error is 0.018%. The finally generated initial flow domain grid includes 11.2 million tetrahedral elements in the catalytic bed area, 4.8 million hexahedral elements in the flow channel area, and 1.5 million prism elements in the interface area. The qualification rate of the Jacobian factor is 99.3%, and the maximum distortion is 0.82.

[0167] As an extensible embodiment, identifying the pore regions in step 10221 that are positively correlated with the propellant decomposition reaction rate in the pore space distribution characteristics of the catalytic bed includes step 102211: Based on the catalyst activity distribution data in the propellant decomposition reaction rate sub-model, screen the smooth wall regions where the pore wall roughness level is lower than the preset threshold; in the smooth wall regions, extract the connected paths where the pore diameter is in the preset target reaction diameter range as the candidate paths for the reaction rate sensitive region; according to the included angle between the extension direction of the candidate path and the fluid flow direction at the inlet of the combustion chamber flow channel, eliminate the paths with an included angle exceeding the preset angle to generate a set of target reaction rate sensitive regions.

[0168] In this embodiment, the screening of the smooth wall region is based on the three-dimensional distribution map of the pore wall roughness grade. A threshold of Ra≤1.2μm is set, and continuous smooth regions with an area >0.1mm² (a total of 58 regions) are extracted through the region growing algorithm. The target reaction diameter range is set to 38±5μm (corresponding to the highest catalytic efficiency range), and the Dijkstra algorithm is used to search for the longest connected path (length 82mm) within this diameter range. The flow angle filtering is calculated through the vector dot product, and paths with an angle >30° with the mainstream direction of the flow channel inlet (positive Z-axis direction) are deleted (accounting for 19%). Finally, 43 optimized paths are retained to form the set of target reaction rate sensitive regions. The spatial distribution data of this set is imported into the CFD preprocessing software in VTK format to guide subsequent mesh generation and boundary condition setting.

[0169] It can be seen that through the multi-scale coupling optimization method in the above embodiment, the coordinated improvement of the spatial single-component engine structure design and reaction flow characteristics is achieved. Based on the above embodiment, after the optimized engine undergoes a simulated ignition cycle, the maximum plastic strain in the heat stress concentration region of the injector head is significantly reduced, the probability of the temperature peak of the catalytic bed support ring is greatly reduced, and the vacuum specific impulse decay rate is improved. Thus, the establishment of the multi-physics field coupling grid model improves the simulation accuracy of fluid flow and heat transfer, providing a high-precision digital platform for the reliability design and life prediction of the aerospace propulsion system.

[0170] In summary, the embodiments of the present application can deeply integrate microscopic structure characteristics, transient chemical reaction kinetics, and dynamic thermodynamic boundaries to achieve high-precision ignition state analysis and processing, thereby improving the engineering applicability of the performance prediction and reliability assessment of the spatial single-component engine, and effectively improving the defects of insufficient multi-physics field coupling accuracy, lack of representation of dynamic boundary conditions, and rough analysis of thermal stress evolution.

[0171] Based on the same inventive concept, the embodiments of the present application also provide an engine ignition state analysis system. Refer to Figure 2 As shown, it is a schematic structural diagram of a possible engine ignition state analysis system provided in the embodiments of the present application. Figure 2 In the figure, the engine ignition state analysis system 200 includes: a processor 210 and a memory 220. Among them, the memory 220 stores a computer program executable by the processor 210. By executing the instructions stored in the memory 220, the processor 210 can execute the steps of the above-mentioned spatial single-component engine ignition state analysis method.

[0172] Based on the same inventive concept, an embodiment of the present application provides a computer-readable storage medium, which includes a computer program. When the computer program runs on an engine ignition state analysis system, the computer program is used to cause the engine ignition state analysis system to execute the steps of the above-mentioned space single-component engine ignition state analysis method. In some possible implementation manners, each aspect of the space single-component engine ignition state analysis method provided by the present application can also be implemented in the form of a program product, which includes a computer program. When the program product runs on an engine ignition state analysis system, the computer program is used to cause the engine ignition state analysis system to execute the steps in the above-mentioned space single-component engine ignition state analysis method. For example, the engine ignition state analysis system can execute as Figure 1 the steps shown therein.

Claims

1. A method for analyzing the ignition state of a single-component space engine, characterized in that, The method includes: Obtaining the initial temperature parameter, initial pressure parameter, and propellant flow parameter of the combustion chamber of a space monopropellant engine in the ignition state; Based on the initial temperature parameter, initial pressure parameter, and propellant flow parameter of the combustion chamber, generating a flow domain grid model including the characteristics of porous medium effect and chemical reaction effect: performing three-dimensional discretization modeling on the pore structure of the catalytic bed of the space monopropellant engine based on the porous medium model to generate topological data of the pore distribution of the catalytic bed; combining the propellant decomposition reaction rate submodel in the chemical reaction model, spatially coupling the topological data of the pore distribution of the catalytic bed with the geometric structure of the combustion chamber to generate an initial flow domain grid with multi-physical field coupling; performing local refinement on the initial flow domain grid according to the mass flow distribution in the propellant flow parameter to generate a flow domain grid model meeting the preset accuracy condition; loading the initial temperature parameter and initial pressure parameter of the combustion chamber as initial field variables into the flow domain grid model to complete the initialization of the flow domain grid model; Setting dynamic boundary conditions corresponding to the ignition transient process in the flow domain grid model, and performing full-cycle numerical simulation on the internal temperature field distribution of the space monopropellant engine to generate a dynamic change sequence of the temperature field in the ignition state; According to the temperature gradient distribution characteristics at different moments in the dynamic change sequence of the temperature field, extracting the thermal stress evolution data of the key components of the space monopropellant engine, and combining the thermal stress evolution data with a preset overheat failure threshold to output the performance stability evaluation result in the ignition state.

2. The method according to claim 1, wherein The performing full-cycle numerical simulation on the internal temperature field distribution of the space monopropellant engine to generate a dynamic change sequence of the temperature field in the ignition state includes: Setting the propellant injection time sequence of the ignition transient process in the flow domain grid model, and calculating the local flow velocity distribution in the porous medium region at each time step based on the pressure fluctuation range in the dynamic boundary conditions; According to the local flow velocity distribution and the heat release rate function in the chemical reaction model, iteratively solving the time evolution model of the temperature field in the combustion chamber to generate temperature field snapshots at continuous time stamps; Performing spatial interpolation on the temperature field snapshots to generate a dynamic change sequence of the temperature field covering the entire flow domain of the space monopropellant engine; Extracting the highest temperature peak value in the catalytic bed outlet region and the occurrence moment of the highest temperature peak value in the dynamic change sequence of the temperature field, and using the highest temperature peak value and the occurrence moment as key indicators of the temperature gradient distribution characteristics.

3. The method according to claim 2, characterized in that, The extracting the thermal stress evolution data of the key components of the space monopropellant engine according to the temperature gradient distribution characteristics at different moments in the dynamic change sequence of the temperature field includes: Calculating the instantaneous thermal expansion coefficient of the combustion chamber wall material according to the spatial temperature gradient data in the dynamic change sequence of the temperature field; Deriving the thermal stress distribution nephogram at different moments based on the instantaneous thermal expansion coefficient and the wall structure constraint conditions; Generate the historical evolution path of the thermal stress concentration region based on the region coordinates exceeding the material yield strength extracted from the thermal stress distribution nephogram and the duration corresponding to the region coordinates; Perform a time series comparison between the historical evolution path and a preset overheat failure threshold, and identify the positions of key components at risk of cumulative damage based on the time series comparison results.

4. The method according to claim 3, characterized in that, Output the performance stability evaluation result in the ignition state, including: Calculate the thermal fatigue cycle times of the key component positions based on the historical evolution path of the thermal stress concentration region; Predict the remaining life interval of the space single-component engine under continuous ignition conditions based on the mapping relationship between the thermal fatigue cycle times and the material durability database; Generate a stability evaluation report including the overheat risk probability and the performance degradation trend by combining the remaining life interval and the key indicators of the temperature gradient distribution characteristics; Mark the results exceeding the preset risk level in the stability evaluation report as structural weak points to be optimized.

5. The method according to claim 1, wherein Perform three-dimensional discretization modeling on the pore structure of the catalytic bed of the space single-component engine based on the porous medium model to generate catalytic bed pore distribution topology data, including: Obtain the geometric parameter set of the catalytic bed, and the geometric parameter set includes the pore diameter distribution range, the curvature of the pore connection path, and the pore density gradient; Identify the multi-scale pore morphology characteristic parameters inside the catalytic bed based on the pore diameter distribution range, and the multi-scale pore morphology characteristic parameters include the axis direction of the main pore channel, the intersection angle of the branch pores, and the roughness of the pore wall; Divide non-uniform grid cells according to the axis direction of the main pore channel and the pore density gradient, and the size of the non-uniform grid cells is adaptively adjusted along the curvature of the pore connection path, and the grid nodes are aligned with the intersection angle of the branch pores; Perform a topological consistency verification on the non-uniform grid cells, judge the continuity of the pore wall roughness data of adjacent grid cells and the fracture state of the pore connection path, and obtain the topological consistency verification result; According to the topological consistency verification result, perform local reconstruction on the grid cells that fail to pass the verification, adjust the grid node positions to match the intersection angle of the branch pores, and generate an optimized porous medium grid model; Extract the pore morphology characteristic parameters, grid size, and node coordinates of each grid cell in the porous medium grid model, and generate catalytic bed pore distribution topology data including pore space distribution, connectivity characteristics, and morphological attributes; Identifying the multi-scale pore morphological characteristic parameters inside the catalytic bed, including: dividing the pores of the catalytic bed into main pore channels and secondary branch pores based on the pore diameter distribution range, where the diameter of the main pore channels is greater than a preset threshold and extends along the axial direction; extracting the centerline trajectory of the main pore channels, calculating the curvature change rate and the deflection angle of the extension direction of the centerline trajectory, and using them as the core parameters of the axial direction of the main pore channels; in the secondary branch pores, detecting the inlet positions connected to the main pore channels, measuring the included angle between the branch pores and the main pore channels at the inlets, and generating a set of intersection angles of the branch pores; quantifying the microscopic undulation characteristics of the pore walls through surface scanning data, and generating a spatial distribution map of the pore wall roughness; Dividing the non-uniform grid cells, including: along the axial direction of the main pore channels, dynamically adjusting the length of the grid cells according to the curvature change rate, where the curvature change rate is negatively correlated with the length of the grid cells; inserting encrypted grid nodes at the intersection angles of the branch pores so that the node density is inversely proportional to the intersection angles of the branch pores; based on the pore density gradient, setting tetrahedral grid cells in the pore-dense areas and hexahedral grid cells in the pore-sparse areas; mapping the spatial distribution map of the pore wall roughness to the surface of the grid cells to generate a wall boundary layer grid with roughness attributes; Verifying the topological consistency of the non-uniform grid cells further includes: judging whether there is a sudden change in the pore wall roughness data of adjacent grid cells at the junction, and if the mutation amplitude exceeds the preset tolerance, marking them as inconsistent cells; verifying the continuity of the grid cells of the main pore channels, where the grid cell continuity characterizes the smooth transition of the curvature change rate of the centerline trajectory between adjacent cells; detecting whether the grid nodes at the intersection angles of the branch pores are aligned, and if there is an offset, calculating the correction amount of the node positions; generating a topological consistency verification report including the positions of the inconsistent cells, the mutation amplitude, and the node offset; 6. The method according to claim 1, characterized in that The method further includes: Collecting the measured temperature data of the space monopropellant engine in the historical ignition tasks and the corresponding performance degradation records of the measured temperature data; Using the measured temperature data to correct the error of the dynamic change sequence of the temperature field, and generating a calibrated temperature prediction model; Jointly training the calibrated temperature prediction model with the thermal stress evolution data to generate an abnormal temperature warning model for real-time monitoring; In the ignition task to be monitored, based on the deviation signal output by the abnormal temperature warning model, dynamically adjusting the propellant flow rate to suppress the overheating risk; Generating the abnormal temperature warning model for real-time monitoring, including: Extracting the residual distribution between the dynamic change sequence of the temperature field and the measured temperature data from the historical ignition tasks; Training a convolutional neural network model using the residual distribution to generate a spatio-temporal feature extractor for temperature prediction residuals; Embedding the spatio-temporal feature extractor into the abnormal temperature warning model to calculate the local difference coefficient between the current temperature field and the predicted value in real time; When the local difference coefficient exceeds the preset safety tolerance, activating the propellant flow rate adjustment instruction to balance the combustion chamber heat load.

7. The method according to claim 6, wherein Based on the deviation signal output by the abnormal temperature warning model, dynamically adjust the propellant flow rate to suppress the overheating risk, including: Combining the spatial distribution of the local difference coefficient and the deviation signal to identify the target area with abnormal temperature increase; Based on the positional information of the target area and the topological relationship of the propellant injection pipeline, calculate the correction amount of the flow valve opening to be adjusted; Send the correction amount of the flow valve opening to the engine control system to re-match the propellant distribution and the current heat load demand through the engine control system, and obtain the adjusted dynamic change sequence of the temperature field; Monitor and verify the adjusted dynamic change sequence of the temperature field, and complete the monitoring and verification when the local difference coefficient drops back within the safe tolerance range.

8. The method according to claim 1, wherein The method further includes: According to the performance stability evaluation result, screen out the failure modes that repeatedly appear in the historical evolution path of the thermal stress concentration area; Generate a target iteration function for the structural optimization of the space single-component engine based on the failure mode, and the target iteration function includes the weight reduction of the thermal stress peak and the material quality constraint condition; Use the genetic algorithm to solve the target iteration function for multiple generations to generate multiple candidate structural improvement solutions that meet the durability requirements; Conduct virtual ignition tests on the multiple candidate structural improvement solutions through the flow domain grid model, and select the target structural improvement solution with the thermal stress distribution meeting the set improvement conditions as the optimization strategy according to the virtual ignition test results; The generation of multiple candidate structural improvement solutions that meet the durability requirements includes: Insert the geometric parameter variables of the additional support structure into the topological data of the catalytic bed pore distribution; Adjust the local flow velocity distribution of the porous medium region according to the geometric parameter variables of the additional support structure; Recalculate the adjusted dynamic change sequence of the temperature field and judge whether the reduction amplitude of the thermal stress peak reaches the preset threshold; If the reduction amplitude of the thermal stress peak reaches the preset threshold, add the structural improvement features corresponding to the current geometric parameter variables to the candidate solution set, otherwise continue the iterative optimization until the termination condition is met.

9. An engine ignition state analysis system, characterized in that, It includes a processor and a memory. Among them, the memory stores a computer program. When the computer program is executed by the processor, the processor executes the steps of any one of claims 1 to 8.

Images