Method and system for analyzing ignition state of spatial monopropellant engine

Abstract

The invention relates to the technical field of data analysis, and provides a spatial monopropellant engine ignition state analysis method and system, and the method comprises the steps: obtaining an initial temperature parameter, an initial pressure parameter and a propellant flow parameter of a combustion chamber of a spatial monopropellant engine in an ignition state; generating a flow domain grid model containing porous medium effect characteristics and chemical reaction effect characteristics based on the initial temperature parameter, the initial pressure parameter and the propellant flow parameter of the combustion chamber; setting a dynamic boundary condition corresponding to an ignition transient process in the flow domain grid model, and performing full-period numerical simulation on internal temperature field distribution of the spatial single-component engine to generate a temperature field dynamic change sequence in an ignition state; and according to the temperature gradient distribution characteristics at different moments in the temperature field dynamic change sequence, extracting thermal stress evolution data of key parts of the spatial single-component engine, and outputting a performance stability evaluation result in an ignition state in combination with the thermal stress evolution data and a preset overheating failure threshold value.

Description

Technical Field

[0001] The present application belongs to the field of data analysis technology, and specifically relates to a method and system for analyzing the ignition state of a space monopropellant engine. Background Art

[0002] In the field of research and development and performance optimization of space monopropellant engines, the analysis of thermodynamic characteristics of ignition transient processes and structural reliability assessment are core technical challenges. In existing technologies, numerical simulations of engine ignition processes usually use simplified physical models. Although such methods can reduce computational complexity, they ignore the effects of the actual pore distribution on the non-uniformity of local flow, heat transfer, and chemical reaction rates, resulting in significantly limited prediction accuracy of temperature and pressure fields.

[0003] In terms of boundary condition setting, existing studies are mostly based on steady-state or quasi-steady-state assumptions, such as setting the combustion chamber wall temperature to a fixed value or using an empirical temperature rise curve, which cannot truly reflect the dynamic heat exchange process between the wall and the high-temperature combustion gas during the ignition process. Regarding thermal stress analysis and performance evaluation, existing technologies usually rely 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 performance degradation (such as high-temperature creep and fatigue crack extension). In addition, the acquisition of propellant flow parameters in existing technologies mostly relies on offline calibration data, and does not consider the real-time impact of supply system pressure fluctuations on mass flow during the ignition transient process, which will lead to systematic deviations between the initial conditions of the numerical model and the actual operating 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 space monopropellant engine simulation and evaluation technology has the defects of insufficient multi-physical field coupling accuracy, lack of dynamic boundary condition characterization, and rough analysis of thermal stress evolution. Summary of the invention

[0005] The present application provides a method and system for analyzing the ignition state of a space monopropellant engine, so as to improve the engineering applicability of performance prediction and reliability evaluation of a space monopropellant engine.

[0006] In the 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, and the method includes: obtaining the initial temperature parameters, initial pressure parameters and propellant flow parameters of the combustion chamber of the space monopropellant engine in the ignition state; based on the initial temperature parameters of the combustion chamber, the initial pressure parameters and the propellant flow parameters, generating a flow domain grid model including porous medium effect characteristics and chemical reaction effect characteristics; setting dynamic boundary conditions corresponding to the ignition transient process in the flow domain grid model, performing a full-cycle numerical simulation of the internal temperature field distribution of the space monopropellant engine, and generating a dynamic change sequence of the temperature field under 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 times in the dynamic change sequence of the temperature field, combining the thermal stress evolution data with a preset overheating failure threshold, and outputting the performance stability evaluation result under 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, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes 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 is run on an engine ignition state analysis system, the computer program is used to enable the engine ignition state analysis system to perform the steps of the above method.

[0009] In the implementation of this application, it is possible to deeply integrate microstructural 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, lack of dynamic boundary condition characterization, and rough thermal stress evolution analysis. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0011] Figure 2 A schematic diagram of the structure of an engine ignition state analysis system provided in an embodiment of the present application. DETAILED DESCRIPTION

[0012] In order to make the purpose, technical solution and advantages of the embodiments of the present application clearer, the technical solution of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the technical solution of the present application, rather than all of the embodiments. Based on the embodiments recorded in the application documents, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the technical solution of the present application.

[0013] See also Figure 1 , which is a space monopropellant engine ignition state analysis method provided in an embodiment of the present application. The method can be applied to an engine ignition state analysis system. The specific process is as follows: Step 101 to Step 104:

[0014] Step 101: Obtaining the initial temperature parameters, initial pressure parameters and propellant flow parameters of the combustion chamber of the space monopropellant engine in the ignition state.

[0015] In an embodiment of the present application, the combustion chamber temperature distribution data at the initial moment of ignition can be obtained by installing a platinum resistance temperature sensor array on the combustion chamber wall of a space monopropellant engine. The platinum resistance temperature sensor array includes 12 groups of PT100 temperature sensors distributed along the axial direction of the combustion chamber. The measurement point spacing is 5 mm, and the measurement range covers 300K to 1200K.

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

[0017] For another example, the propellant flow parameters are monitored in real time by a Coriolis mass flow meter installed downstream of the pressure stabilization section of the propellant supply pipeline. The measurement accuracy reaches ±0.5% of the reading, and the mass flow, temperature and density parameters of the hydrazine-based propellant can be accurately obtained. The experiment measured that the initial temperature of the combustion chamber 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 ​​of the subsequent numerical model.

[0018] In the embodiment of the present application, the propellant flow parameters refer 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 RheonikRHM015) can be used for accurate measurement. The 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. The measurement accuracy reaches ±0.15% of the range.

[0019] In the specific implementation, a dual straight tube sensor is installed downstream of the pressure stabilization section of the hydrazine-based propellant supply pipeline to ensure that stable data is obtained after the flow is fully developed. Exemplary data include propellant inlet temperature (293±2K), density (1.004g / cm³) and mass flow rate (12.5±0.3g / s), which are transmitted to the control terminal via the ModbusRTU protocol as the initial conditions of the mass source term and energy function in the computational fluid dynamics (CFD) model. For example, when setting material properties in the Fluent software, the measured density value is substituted into the state function, and the mass flow 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 combustion chamber initial temperature parameter, the initial pressure parameter and the propellant flow parameter, a flow domain grid model including porous medium effect characteristics and chemical reaction effect characteristics is generated.

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

[0022] In view of the porous medium effect, microfocus X-ray computed tomography technology was used to obtain the real pore structure of the catalyst bed. The reconstructed three-dimensional pore model was imported into COMSOL Multiphysics software, and the random stacking algorithm was applied to generate a non-uniform porous medium grid. The unit size was controlled at 50μm to accurately characterize the porosity gradient distribution.

[0023] In addition, the chemical reaction effect modeling adopts a detailed kinetic mechanism to decompose the hydrazine decomposition reaction path into 12 elementary reaction steps, such as N2H4→NH2+NH2, NH2→NH+H, etc., and imports it into Fluent software through the CHEMKIN format to establish a finite rate chemical reaction model.

[0024] The final generated hybrid grid model contains about 8.5 million hexahedral units, of which the catalyst bed area uses an adaptive encrypted grid with a minimum unit size of 20μm, and the injector swirl channel uses an O-type topology grid to ensure the accuracy of flow direction analysis.

[0025] In step 102, the porous media effect characteristics refer to the influence of the complex pore structure in the catalyst bed on the fluid flow and heat and mass transfer process. In the embodiment of the present application, micro-focus X-ray computed tomography (μCT, model ZEISS Xradia 520 Versa) is used to obtain the real three-dimensional pore structure of the iridium / aluminum oxide catalyst, and the scanning resolution reaches 2μm / voxel. For another example, the porosity (ε=0.38), specific surface area (S v =1.2×10^4m² / m³) and tortuosity (τ=1.85). When establishing an asymmetric porous media model in COMSOL Multiphysics, the Brinkman function is coupled with the Forchheimer correction term to describe the flow resistance: .

[0026] in, represents 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; It is the tensor product (outer product) of the velocity vector; p (pressure): unit: Pa (Pascal), represents the static pressure of the fluid; μ (dynamic viscosity): unit: Pa·s, represents the viscous resistance of the fluid; K (permeability): unit: m 2 , describes the ability of porous media to allow fluid to pass through, which is related to the pore structure; β (Inertia coefficient): Unit: m -1 , the coefficient of the Forchheimer correction term, characterizing the influence of inertial effect on resistance at high flow rate; (Velocity modulus): 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. The Voronoitessellation algorithm is used to create an unstructured grid in the grid generation stage, and the cell size gradient is controlled at 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 coupling effect of the exothermic decomposition reaction of the propellant in the catalytic bed on the flow heat transfer. The embodiment of the present application adopts a detailed chemical reaction mechanism to decompose the decomposition process of hydrazine (N2H4) into 12 elementary reaction steps, including:

[0029] 1. N2H4→2NH2(Ea=58.6kJ / mol) 2. NH2→NH+H (Ea=102.4kJ / mol) ... The reaction rate constants of each step are calculated by transition state theory and imported into Fluent's finite rate chemical reaction model using the CHEMKIN-II format. A surface reaction model is set up in the catalyst bed area, and the active site density of the iridium catalyst is defined as 5.6×10^18sites / m². The apparent reaction rate is calculated in combination with the Langmuir-Hinshelwood adsorption mechanism. For example, the catalytic reaction rate expression in UDF is:

[0030] .

[0031] Where r is the catalytic reaction rate; k0 is the prefactor, which is the kinetic constant in the Arrhenius formula 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 thermodynamic temperature and energy scale; T is the catalyst surface temperature, which is used for the local temperature of the catalyst bed and directly affects the reaction rate; P N2H4 is the partial pressure of hydrazine (N2H4), that is, the partial pressure of the reactants in the gas phase, driving the reaction forward. NH3 is the partial pressure of ammonia (NH3), which is used to control the reaction rate by the partial pressure of the reaction products in the gas phase.

[0032] Further, , Ea=67.8kJ / mol, K ads is the ammonia adsorption equilibrium constant. The model can accurately predict the local temperature rise phenomenon caused by the reaction exotherm. For example, the simulation shows that the temperature of the hot spot area of ​​the catalyst bed rises from 650K to 1250K within 20ms after ignition.

[0033] Step 103: setting dynamic boundary conditions corresponding to the ignition transient process in the flow domain grid model, performing full-cycle numerical simulation on the internal temperature field distribution of the spatial monopropellant engine, and generating a dynamic change sequence of the temperature field under the ignition state.

[0034] In step 103, when setting the dynamic boundary conditions, the combustion chamber wall boundary can be defined as a time-varying temperature boundary condition based on exemplary data of the engine ignition transient process, and 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 a given transient pressure curve includes a 0.5MPa step rise in the initial stage and a subsequent 1.2MPa / s linear growth segment, which is obtained by fitting the measured pressure data.

[0035] Furthermore, the numerical simulation uses a coupled implicit solver, and the time step is set to 1μs to meet the CFL condition of sound wave propagation. The k-omegaSST turbulence model is enabled and the DO radiation model is activated when solving the NS function. The porous medium momentum source term and the custom chemical reaction rate UDF are enabled in the catalyst bed area, and the catalyst surface reaction activation energy is corrected as a function of temperature through a user-defined function. The full-cycle simulation lasts 500ms of physical time, and the full-field temperature data is output every 0.1ms. Finally, a data sequence containing 5000 transient fields is generated, which fully records the temperature evolution characteristics of the injector swirl area, the catalyst bed hot spot area and the nozzle throat during the ignition process.

[0036] It can be understood that the dynamic change sequence of the temperature field refers to a three-dimensional temperature distribution data set recorded in chronological order throughout the engine during the ignition process. The embodiment of the present application uses a coupled solver 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 catalyst bed area, and the step size is automatically reduced to 0.2μs when the temperature change rate exceeds 50K / ms. The data output is stored in HDF5 format, and each time slice contains the temperature values ​​of 1.2×10^6 grid nodes, with a spatial resolution of 0.1mm. The typical data sequence shows that the injector swirl zone forms a 340K low-temperature core flow 5ms after ignition, and a 1560K high-temperature group with a diameter of 2mm appears at the catalyst bed outlet at 50ms. The high-temperature group migrates downstream with the flow and reaches the nozzle throat at 120ms.

[0037] Step 104: Extract thermal stress evolution data of key components of the space monopropellant engine based on the temperature gradient distribution characteristics at different moments in the temperature field dynamic change sequence, combine the thermal stress evolution data with a preset overheating failure threshold, and output the performance stability evaluation result under the ignition state.

[0038] In step 104, the thermal stress analysis phase first maps the temperature field data of the temperature field dynamic change sequence to the finite element model of ANSYS Mechanical, which contains the detailed structural features of the injector, combustion chamber and nozzle assembly. For the injector head made of GH3128 high-temperature alloy, its temperature gradient distribution data is extracted and substituted into the thermoelastic constitutive function to calculate the transient thermal stress field.

[0039] It should be noted that special attention should be paid to the maximum temperature gradient zone at the root of the injector swirl vane 50ms after ignition. The axial temperature gradient in this area reaches 1.2×10^5K / m, and the equivalent thermal stress generated is calculated to be 485MPa. 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 in an environment of 850K is 620MPa, the actual thermal stress value reaches 78.2% of the material strength threshold. The performance stability evaluation module triggers an early warning when the local temperature exceeds the material recrystallization temperature (923K) or the thermal stress exceeds 80% of the yield strength according to the preset overheating failure criteria.

[0040] For example, the simulation results show that 320ms after ignition, a transient temperature spike (peak temperature 935K) lasting 0.8ms appears at the catalyst bed support ring. 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 the temperature excess, providing a quantitative basis for subsequent structural optimization.

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

[0042] Furthermore, the performance stability evaluation result refers to a quantitative safety criterion analysis report based on thermodynamic parameters and material limits. The embodiment of the present application can set up a two-level early warning mechanism: the first-level early warning trigger condition is that the local temperature exceeds the material recrystallization temperature (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 a moving time window statistical method. When the catalyst bed support ring is at 320ms and a temperature spike occurs (935K for 0.8ms), the system marks the event as a third-level abnormality (cumulative exceeding time <1ms). The final output includes an evaluation report containing the coordinates of the risk area, a list of exceeding parameters, and recommended improvement measures (such as increasing the density of the cooling channel at this location), providing data support for engine optimization design.

[0043] It can be seen that the embodiment of the present application has achieved a refined simulation of the ignition process of a space monopropellant engine by establishing a high-precision multi-physics field coupling model. The obtained temperature field dynamic change sequence and thermal stress evolution data 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 material performance degradation under transient thermal shock can directly guide the thermal protection design of key engine components and effectively prevent problems such as catalyst sintering failure or structural creep deformation caused by local overheating.

[0044] As an optional but not limiting embodiment, the step 102 of generating a flow domain grid model including porous medium effect characteristics and chemical reaction effect characteristics based on the combustion chamber initial temperature parameter, the initial pressure parameter and the propellant flow parameter includes: Step 1021: Perform three-dimensional discretization modeling on the pore structure of the catalyst bed of the spatial monopropellant engine based on the porous medium model to generate catalyst bed pore distribution topology data.

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

[0046] For example, based on fractal geometry theory, a random walk algorithm can be used to reconstruct the three-dimensional pore network model of the catalyst bed and generate topological connection data containing 1.2×10^6 pore nodes. The topological data records the coordinate position, equivalent diameter (range 15-85μm) and adjacent pore connectivity information of each pore unit. After being exported as an STL format file, it is imported into the COMSOL Multiphysics software. The pore model and the solid geometry of the catalyst bed are subtracted through Boolean operations to finally form a discretized porous medium digital twin, whose porosity gradient distribution has a relative error of less than 3.5% compared with the exemplary data.

[0047] Step 1022: In combination with the propellant decomposition reaction rate sub-model in the chemical reaction model, the catalytic bed pore distribution topology data is spatially coupled with the combustion chamber geometry to generate an initial flow domain mesh for multi-physics field coupling.

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

[0049] Furthermore, multi-physics coupling can be achieved through APDL scripts, and the chemical reaction source term can be bidirectionally coupled with the porous media momentum function. The momentum loss term in the porous media region adopts the Ergun function correction model, and the chemical reaction exothermic power is calculated in real time by UDF and injected as the energy function source term. The final generated hybrid grid model includes the hexahedral dominant grid of the injector head (unit size 0.2mm), the tetrahedral adaptive grid of the catalytic bed area (minimum unit size 20μm) and the prismatic layer grid of the transition zone of the combustion chamber (total number of layers 15, growth rate 1.2), the total number of grids reaches 11.2 million units, and the Jacobian factor is greater than 0.82.

[0050] Step 1023: Based on the mass flow distribution in the propellant flow parameters, the initial flow domain grid is locally encrypted to generate a flow domain grid model that meets the preset accuracy conditions; the initial temperature parameters and the initial pressure parameters of the combustion chamber are loaded into the flow domain grid model as initial field variables to complete the initialization of the flow domain grid model.

[0051] In step 1023, the mesh encryption process is implemented based on the mass flow distribution characteristics in the propellant flow parameters. The specific method is to set the dynamic encryption criteria in the adapter mesh module of the Fluent software. When the mass flow gradient in the injector swirl channel is monitored to exceed 50g / (s·m²), the local mesh refinement operation is triggered, and the mesh size of the trailing edge area of ​​the swirl blade is encrypted from 0.5mm to 0.1mm.

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

[0053] As an optional but not limiting embodiment, the full-cycle numerical simulation of the internal temperature field distribution of the spatial monopropellant engine in step 103 to generate a dynamic change sequence of the temperature field under the ignition state includes: Step 1031: setting the propellant injection time series of the ignition transient process in the flow domain grid model, and calculating the local flow velocity distribution of the porous medium region in each time step based on the pressure fluctuation range in the dynamic boundary conditions.

[0054] 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 the Fluent software. The function is divided into three stages: 0-10ms is the pressure linear rise stage (slope 0.5MPa / ms), 10-50ms is the pressure holding stage (constant value 2.1MPa), and 50-500ms is the pressure regulation stage (dynamically adjusted according to the PID control algorithm).

[0055] For example, in each time step (1μs), when solving the local velocity distribution in the porous media area, the coupled implicit algorithm is used to solve the Brinkman-Forchheimer expansion 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 of the catalyst 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 times. For another example, 20ms after ignition, the velocity at the inlet of the catalyst bed reaches a peak of 12.3m / s, and the corresponding local Reynolds number Re=210 (based on the equivalent pore diameter), and the flow is in a laminar to transitional state.

[0056] Step 1032: Iteratively solve the time evolution model of the temperature field in the combustion chamber according to the local flow velocity distribution and the heat release rate function in the chemical reaction model, and generate snapshots of the temperature field under continuous time stamps.

[0057] In step 1032, the temperature field evolution model is solved by using a strong coupling algorithm, specifically, enabling the simultaneous solution of the energy function and the component transport function in Fluent. In each time step, the local velocity distribution in the porous medium region is first calculated, and then the exothermic power of the hydrazine decomposition reaction is calculated by a finite rate chemical reaction model.

[0058] Exemplarily, the heat release rate function is implemented by UDF, and its core expression is Q=ΔH·r·As, where: ΔH=-1540kJ / kg is the enthalpy change of 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 catalyst bed. The k-omega SST format is used for the turbulence model, the enhanced wall function is used for near-wall treatment, the radiation heat transfer is calculated by the discrete coordinate model (DO model), and the number of radiation iterations is set to 5 times / time step.

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

[0060] Step 1033: perform spatial interpolation processing on the temperature field snapshot to generate a temperature field dynamic change sequence covering the entire flow domain of the spatial monocomponent engine; extract the highest temperature peak value in the catalyst bed outlet area and the time of occurrence of the highest temperature peak value in the temperature field dynamic change sequence, and use the highest temperature peak value and the time of occurrence as key indicators of the temperature gradient distribution characteristics.

[0061] In step 1033, the spatial interpolation process uses a cubic spline interpolation algorithm, specifically, full-field data remapping is performed in the Tecplot 360 software. The temperature field snapshot of the unstructured grid is interpolated to a uniform Cartesian grid (resolution 0.1mm×0.1mm×0.1mm) to eliminate the numerical noise caused by grid distortion. The generation of the dynamic change sequence of the temperature field is achieved by writing a Python script, which splices 5000 interpolated temperature fields into a four-dimensional array (X, Y, Z, T) in chronological order. The highest temperature peak extraction algorithm in the catalyst bed outlet area is based on the regional growth method, and the temperature threshold is set to 1200K. When three consecutive grid nodes are monitored to exceed the threshold, it is determined to be a valid peak.

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

[0063] As an optional but not limiting embodiment, the step 104 of 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 includes: Step 1041: Calculate the instantaneous thermal expansion coefficient of the combustion chamber wall material based on the spatial temperature gradient data in the temperature field dynamic change sequence.

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

[0065] For example, the fitting formula of the instantaneous thermal expansion coefficient curve is: α(T)=13.8×10^-6+0.023×10^-6·T+1.7×10^-9·T² (unit: ), with a correlation coefficient of R²=0.998. During the calculation process, the temperature gradient data of the combustion chamber wall node in the temperature field dynamic change sequence (such as the axial temperature gradient of the injector head 1.2×10^5K / m 50ms after ignition) is substituted into the formula to obtain the instantaneous value of the thermal expansion coefficient at the corresponding position: The data conversion is realized through the ANSYS Workbench platform, and 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%.

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

[0067] In step 1042, the derivation of the thermal stress distribution cloud map is based on the thermoelastic constitutive function, and fixed constraints are set in the ANSYS Mechanical transient structure module: the combustion chamber flange mounting surface is subject to full degree of freedom constraints, and the nozzle outlet end face is subject to axial displacement restrictions. The material model uses the bilinear kinematic hardening criterion, and defines the yield strength of the GH3128 alloy at 850K as 620MPa and the tangent modulus as 12GPa.

[0068] 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 (0.1ms interval), the upper limit of the number of iterations in each time step is set to 25 times, and the convergence residual standard is 1×10^-4. For example, 120ms after ignition, the maximum Von Mises stress value of 583MPa appears at the root of the injector swirl plate, the stress concentration factor in this area reaches 2.3, and the corresponding plastic strain accumulation is 0.15%.

[0069] Step 1043: Generate a 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 cloud map and the duration corresponding to the area coordinates.

[0070] In step 1043, the generation of the historical evolution path of the thermal stress concentration area is automated through a Python script. The specific process is as follows: first, the unit number and coordinate information of the unit exceeding the material yield strength threshold (620 MPa) in the thermal stress distribution cloud map at each time step are extracted, and then the moving time window algorithm (window width 5 ms, step length 0.1 ms) is used to count the duration of exceeding the standard at each spatial position.

[0071] For the injector head area (coordinates X=32.5-37.8mm, Y=10.2-15.6mm), analysis shows that stress exceeds the limit during 48-52ms after ignition, with the maximum continuous exceeding limit duration of 3.2ms and a spatial coverage area of ​​8.7mm². The historical path data is stored in the form of a three-dimensional point cloud sequence. Each data point contains spatial coordinates, stress peak value and action time. The spatiotemporal evolution animation generated by Paraview software can clearly show the trajectory of the stress concentration area migrating from the root of the swirl plate to the center of the injector.

[0072] Step 1044: Perform a time series comparison between the historical evolution path and a preset overheat failure threshold, and identify the positions of key components with cumulative damage risks based on the time series comparison results.

[0073] In step 1044, the time series comparison analysis uses an improved rain flow counting method to decompose the thermal stress history evolution path into full-cycle and half-cycle events. The preset overheating 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 catalyst bed support ring area (coordinate Z=85-92 mm), three events in which the stress peak exceeded the dynamic threshold during the period of 320-320.8 ms were identified, and the single peak stress of 625 MPa corresponded to a damage increment of ΔD=0.0047. When the total cumulative damage reaches 0.82, the system automatically marks the area as a high-risk area and marks it with a red contour line in the three-dimensional model, and generates a JSON format report containing the location coordinates, the number of times the standard is exceeded, and the degree of damage.

[0074] As an optional but not limiting embodiment, the outputting of the performance stability evaluation result under the ignition state in step 104 includes: Step 1045: Calculate the number of thermal fatigue cycles at the key component position according to the historical evolution path of the thermal stress concentration area.

[0075] In step 1045, the calculation of the number of thermal fatigue cycles is implemented based on the Coffin-Manson-Basquin function, which is as follows: 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 is far beyond the allowable value D=1, indicating that there is a serious risk of fatigue failure in this part.

[0076] 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.

[0077] In step 1046, the remaining life prediction adopts the probability statistics method, calling the GH3128 alloy low-cycle fatigue test data in the material durability database (sample size n=120, confidence level 95%). The shape parameter β=2.1 and the scale parameter η=2350cycles are obtained by fitting the Weibull distribution. The calculated equivalent cycle number 1.8 times / ignition is substituted into the reliability function R(t)=exp[-(N / η)^β], and the reliability after running 1000 ignitions is 67.3%, and the remaining life range is [832, 1215] times (Bootstrap method calculation). The result is written into the evaluation system through the MySQL database interface and updated in conjunction with the real-time monitoring data.

[0078] 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.

[0079] In step 1047, the generation of the stability assessment report integrates multi-source data: first, the historical evolution path data of the thermal stress concentration area is analyzed in time and space with the highest temperature peak (1563K) at the catalyst bed outlet in the temperature field dynamic change sequence. It is found that when the temperature peak exceeds 1500K, the probability of stress exceeding the standard within the subsequent 5ms increases to 78%. The Monte Carlo method is used to simulate the 1000 ignition processes, and the probability of thermal fatigue cracking of the injector head after 500 ignitions is calculated to be 92.7%, and the performance attenuation trend is a specific impulse decrease of 0.8% / 100 cycles. The report output is in PDF format, including a three-dimensional thermal-mechanical coupling cloud map, SN curve fitting results, and a histogram of the remaining life distribution at key locations.

[0080] Step 1048: Mark the results in the stability assessment report that exceed a preset risk level as structural weaknesses that need to be optimized.

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

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

[0083] As an optional but not limiting embodiment, the three-dimensional discretization modeling of the pore structure of the catalytic bed of the spatial monopropellant engine based on the porous medium model in step 1021 to generate catalytic bed pore distribution topology data includes: Step 1021: Obtain a set of geometric parameters of the catalytic bed, wherein the set of geometric parameters includes a pore diameter distribution range, a pore connection path curvature, and a pore density gradient.

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

[0085] In addition, the image can be processed by Avizo Fire software, and the pores and solid skeleton can be segmented by region growing algorithm. The pore diameter distribution range is calculated to be 15-85 μm (peak diameter 38 μm), and the statistical mean of the curvature of the pore connection path is The pore density gradient increases linearly along the axial direction from 0.32 at the inlet to 0.41 at the outlet.

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

[0087] Step 1022: Based on the pore diameter distribution range, identify multi-scale pore morphology characteristic parameters inside the catalyst bed, wherein the multi-scale pore morphology characteristic parameters include the main pore channel axis direction, branch pore intersection angle and pore wall roughness.

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

[0089] The detection of secondary branch pores can be based on Euclidean distance transformation, identifying pores with a distance less than 10μm from the main pore channel as connecting branches, and measuring their intersection angle distribution range is 28°-152° (mean 67°). The quantification of pore wall roughness uses a three-dimensional surface morphology analysis module to perform Gaussian filtering (cutoff wavelength 5μm) on the segmented pore wall point cloud data, calculate the surface roughness Ra=1.8μm, and generate a two-dimensional Fourier spectrum showing that the main spatial frequency components are concentrated in interval.

[0090] Step 1023: Divide the 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 branch pores.

[0091] In step 1023, the division of the non-uniform grid cells is performed by the grid generator of COMSOL Multiphysics, and a dynamic adjustment rule is set: along the axis direction of the main pore channel, when the curvature change rate exceeds When the grid unit length is reduced from 50μm to 20μm, in the area where the branch pore intersection angle is less than 60°, the node spacing is reduced to 10μm by inserting encrypted nodes.

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

[0093] Step 1024: Perform topological consistency verification on the non-uniform grid cells, determine the continuity of the pore wall roughness data of adjacent grid cells and the fracture state of the pore connection path, and obtain a topological consistency verification result.

[0094] In step 1024, the topology consistency verification process includes multi-level detection: first, the adjacent grid cells are traversed through the Python script to calculate the root mean square deviation of the roughness data at the junction. When the deviation exceeds the preset tolerance of 0.5μm, it is marked as an inconsistent cell (such as the deviation at coordinates X=125μm, Y=63μm is 0.73μm).

[0095] Furthermore, the continuity of the main pore channel is verified by using the centerline curvature difference method to detect the difference in curvature change rate between adjacent units exceeding The abnormal sections (such as 3 mutations in the axial position Z=520-535μm interval). The branch pore intersection angle alignment test calculated the node position offset and found that the node offset at coordinates X=88μm and Y=214μm was 7.3μm. The final verification report contains 17 inconsistent unit positions, 23 curvature mutation points and 45 offset nodes, all recorded in the form of three-dimensional coordinates.

[0096] Step 1025: Based on the topological consistency verification result, locally reconstruct the grid units that have not passed the verification, adjust the grid node positions to match the branch pore intersection angles, and generate an optimized porous medium grid model.

[0097] In step 1025, local reconstruction operations are performed on the grid cells that have not passed the verification: for the cells with sudden roughness changes, the Laplacian smoothing algorithm is used to adjust the node positions so that the roughness gradient at the junction is reduced to within 0.4 μm; transition grid cells are inserted into the region where the curvature of the main pore channel suddenly changes, and the difference in the curvature change rate is controlled within The following figure shows that the displacement of the branch pore offset node is optimized by the least square method, and the node movement does not exceed 5μm. In the reconstructed porous media grid model, the number of inconsistent units is reduced from the initial 85 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%.

[0098] Step 1026: extract the pore morphological characteristic parameters, grid size and node coordinates of each grid unit in the porous medium grid model, and generate catalytic bed pore distribution topology data including pore spatial distribution, connectivity characteristics and morphological attributes.

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

[0100] In the embodiment of the present application, the grid size data is counted by region, with an average size of 28 μm in the main channel area, 15 μm in the branch area, and 5 μm in the boundary layer area. The node coordinate data is stored in floating point format with an accuracy of 0.1 μm. The final data set is packaged in HDF5 format, contains 1.2×10^6 records, occupies 23GB of storage space, and is indexed and associated to subsequent simulation modules through the MySQL database.

[0101] As an optional but not limiting embodiment, the identification of multi-scale pore morphological characteristic parameters inside the catalytic bed in step 1022 includes step 10220: based on the pore diameter distribution range, dividing the catalytic bed pores into main pore channels and secondary branch pores, the diameter of the main pore channel is greater than a preset threshold and extends axially; extracting the centerline trajectory of the main pore channel, calculating the curvature change rate and the deflection angle of the extension direction of the centerline trajectory as the core parameters of the axial direction of the main pore channel; in the secondary branch pores, detecting the entrance position connected to the main pore channel, measuring the angle between the branch pores at the entrance and the main pore channel, and generating a set of branch pore intersection angles; quantifying the microscopic undulation characteristics of the pore wall through surface scanning data, and generating a spatial distribution map of the pore wall roughness.

[0102] In the embodiment of the present application, the division of the main pore channel and the secondary branch pores is achieved by the adaptive threshold method: Gaussian fitting is performed on the pore diameter distribution histogram, 45 μm is determined as the bimodal distribution valley point, and 9532 pores with a diameter ≥ 45 μm are classified as the main channel. The centerline trajectory extraction uses 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: , the maximum deflection angle occurs at Z = 325 μm (23.5°).

[0103] For another example, the secondary branch entrance detection was based on morphological corrosion operations, identifying 18,920 connection points, and the average angle between them and the main channel was 67°, with a standard deviation of 14°. The wall roughness analysis used wavelet transform to decompose the surface morphology, extracting the detail component of scale 3 (wavelength 4-8μm) as the roughness characterization, and generating a two-dimensional distribution cloud map showing that the roughness of the entrance area is 18% higher than that of the exit area.

[0104] As an optional but not limiting embodiment, the division of non-uniform grid units in step 1023 includes step 10230: dynamically adjusting the grid unit length according to the curvature change rate along the axis of the main pore channel, and the curvature change rate is negatively correlated with the grid unit length; inserting encrypted grid nodes at the branch pore intersection angle so that the node density is inversely proportional to the branch pore intersection angle; based on the pore density gradient, setting tetrahedral grid units in the pore-dense area and setting hexahedral grid units in the pore-sparse area; mapping the spatial distribution map of the pore wall roughness to the grid unit surface to generate a wall boundary layer grid with roughness attributes.

[0105] In the above embodiment, the specific operation of non-uniform meshing includes: when the curvature change rate of the main pore channel exceeds In the section with narrow branch intersection angle of 28° (such as Z=120-135μm), the length of the grid unit is gradually reduced from 50μm to 20μm; in the narrow area with branch intersection angle of 28° (X=205μm, Y=178μm), the node density is set to 150 nodes / mm²; the pore-dense area (density>0.38) generates a tetrahedral mesh with a minimum size of 15μm, and the sparse area uses a 50μm hexahedral mesh. The wall roughness mapping is realized by the interpolation algorithm, and the Ra value is assigned as an additional attribute of the mesh surface node. The wall function parameters are read and corrected by UDF in Fluent.

[0106] As an optional but not limiting embodiment, the topological consistency verification of the non-uniform grid unit in step 1024 also includes step 10240: determining whether there is a mutation in the pore wall roughness data of adjacent grid units at the junction, and marking the unit as inconsistent if the mutation amplitude exceeds a preset tolerance; verifying the grid unit continuity of the main pore channel, and the grid unit continuity characterizes the smooth transition of the curvature change rate of the centerline trajectory between adjacent units; detecting whether the grid nodes at the intersection angle of the branch pores are aligned, and calculating the node position correction if there is an offset; generating a topological consistency verification report including inconsistent unit positions, mutation amplitudes and node offsets.

[0107] For example, the detailed process of topological consistency verification can include: using the differential method to calculate the gradient of the roughness data at the junction of adjacent grid units, triggering an alarm when the absolute value of the gradient exceeds 0.5μm / μm (such as the gradient of 0.63μm / μm at X=152μm and Y=93μm); the main channel continuity detection calculates the angle of the tangent vectors of the center lines of adjacent units, and it is judged as discontinuous if it exceeds 5° (such as the angle of 7.2° at Z=430μm); the branch node alignment verification uses the nearest neighbor search algorithm, the offset tolerance is set to 2μm, and 45 nodes that exceed the standard are detected. The verification report is output in XML format, including inconsistent unit IDs, spatial coordinates, and quantized deviation values ​​for the grid optimization module to call.

[0108] As an optional but not limiting embodiment, the method further includes: Step 201: Collecting the measured temperature data of the space monopropellant engine in the historical ignition mission and the performance degradation records corresponding to the measured temperature data.

[0109] For example, the measured temperature data of historical ignition missions can be collected by installing a platinum resistance temperature sensor array (model PT100, measurement point spacing 5 mm) 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.

[0110] 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-10MPa and the linearity error was ±0.3%FS. The performance degradation record included the activity decay rate of the iridium catalyst in the catalyst bed (measured by X-ray photoelectron spectroscopy). The content of GH3128 alloy was reduced from 78% to 62%), the crack extension length of the injector head GH3128 alloy (the maximum crack length measured by metallographic microscope was 0.82mm), and the specific impulse reduction data (the vacuum specific impulse was reduced from 230s to 225s). All data are stored in HDF5 format, and an associated database containing timestamps, spatial coordinates, temperature values ​​and corresponding performance indicators is established.

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

[0112] 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 by MATLAB's Curve Fitting Toolbox. It is found that there is a systematic deviation in the catalyst bed outlet area (the maximum deviation value is 87K).

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

[0114] Step 203: jointly train the calibrated temperature prediction model and the thermal stress evolution data to generate an abnormal temperature early warning model for real-time monitoring.

[0115] In step 203, the abnormal temperature warning model is constructed based on the PyTorch framework. The input layer receives a 512×512×500 three-dimensional temperature field sequence (spatial resolution 0.1mm, time resolution 0.1ms), 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 convolution layer (kernel size 3×3×3, step size 1) to extract spatiotemporal features, and the output dimension of the fully connected layer matches the number of grid nodes (1.2×10^6 dimensions).

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

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

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

[0119] For example, the adjustment command can be sent to the electric control valve (model Fisher Vee-Ball V300) of the propellant supply system through the Modbus TCP protocol to increase the valve opening of the No. 3 injection unit from 45% to 58%, so that the mass flow rate of propellant in this area increases to 14.1g / s. After adjustment, the monitoring data shows that the temperature in the target area drops from 935K to 863K within 120ms, and the local difference coefficient drops back to 0.09. The warning state is automatically lifted and the adjustment log is recorded.

[0120] As an optional but not limiting embodiment, the step 203 of generating an abnormal temperature warning model for real-time monitoring includes: Step 2031: extracting the residual distribution between the temperature field dynamic change sequence and the measured temperature data from the historical ignition tasks.

[0121] In step 2031, the residual distribution extraction uses a spatial registration algorithm to interpolate the measured temperature data of the historical ignition task to the grid nodes of the computational fluid dynamics model and calculate the absolute deviation ΔT=|T sim -T exp |. Statistics show that the peak value of the residual in the catalyst bed outlet area during the period of 50-60ms after ignition reaches 92K. The deviation shows the characteristics of axial propagation (propagation speed 1.8m / s) and has a strong correlation with the combustion chamber pressure oscillation frequency of 2.4kHz. The residual distribution database contains 500 transient fields, each field records the deviation values ​​of 1.2×10^6 nodes, and the total data volume is 3.4TB.

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

[0123] In step 2032, the training process of the convolutional neural network model adopts a distributed computing architecture and uses 4 NVIDIA A100 GPUs for parallel processing. The input data is a 512×512×500 three-dimensional temperature field sequence. After 5 layers of three-dimensional convolution (channel number 64-512) and 3 layers of full connection processing, the residual prediction field of the same dimension is output. The loss function is defined as the weighted mean square error: , where the weight w i It is set to 3.0 in the catalytic bed area and 1.0 in other areas. After training, the model achieves an average absolute error of 4.3K and a peak error of 13.7K on the test set.

[0124] Step 2033: embed the spatiotemporal 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.

[0125] In step 2033, the spatiotemporal feature extractor is embedded through a dynamic link library, and the trained neural network model is converted into 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.1ms, and a 512-dimensional feature vector is output through the feature extractor, which is then mapped to 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, w i is the attention weight obtained through training. The system sets dual thresholds: when C>0.15, a yellow warning is triggered, and when C>0.25, a red warning is triggered and emergency adjustment is initiated.

[0126] Step 2034: When the local difference coefficient exceeds a preset safety tolerance, the propellant flow adjustment instruction is activated to balance the thermal load of the combustion chamber.

[0127] In step 2034, the generation of the propellant flow adjustment command is based on the fuzzy control rule base, which defines the input variable as the local difference coefficient C and its change rate dC / dt, and the output variable as the flow correction Δ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.2g / s; if C≥0.25, then Δm=2.5g / s and start auxiliary cooling. The command transmission delay is controlled within 0.8ms to ensure intervention in the early stage of temperature anomaly development.

[0128] As an optional but not limiting embodiment, the step 204 of dynamically adjusting the propellant flow rate to suppress the overheating risk based on the deviation signal output by the abnormal temperature warning model includes: Step 2041: Combining the spatial distribution of the local difference coefficient and the deviation signal, identify the target area with abnormally high temperature.

[0129] For example, the target area identification can use the connected domain analysis method to cluster the grid nodes with a local difference coefficient exceeding 0.15, and set the minimum cluster volume to 5 mm³. For example, 320ms after ignition, an abnormally high temperature cluster with a diameter of 3.2mm was detected in the catalyst bed support ring area (coordinate Z=85-92mm), containing 126 excessive nodes with an average difference coefficient of 0.21. Spatial distribution analysis shows that this area has a direct fluid path association with the No. 3 and No. 7 injection units of the propellant supply pipeline.

[0130] Step 2042: Calculate the flow valve opening correction amount to be adjusted based on the position information of the target area and the topological relationship of the propellant injection pipeline.

[0131] For example, the calculation of the flow valve opening correction is based on the fluid network model, and the propellant pipeline pressure-flow characteristic function is established: , where K v is the valve flow coefficient, ρ=1.004g / cm³ is the density of hydrazine propellant. For injection unit No. 3, the current opening of 45% corresponds to K v =0.62, need to adjust to 58% opening to make K v =0.89, thereby increasing the mass flow rate from 12.5g / s to 14.1g / s. The correction calculation module simultaneously considers the pipeline transmission delay (3.2ms) and valve response time (8ms), and generates a feedforward compensation command to act 1.5 cycles in advance.

[0132] Step 2043: Send the flow valve opening correction amount 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 falls back to a safe tolerance range.

[0133] For example, the response verification of the engine control system is performed through the hardware-in-the-loop test platform, and the adjusted flow parameters are input into the real-time simulation model to monitor the evolution of the dynamic change sequence of the temperature field. For example, after the valve opening is adjusted, the temperature gradient in the target area drops from 1.1×10^5K / m to 7.2×10^4K / m, and the peak thermal stress drops from 625MPa to 538MPa. The system continuously monitors the local difference coefficient within a 50ms period. When the coefficient of three consecutive sampling points is less than 0.12, the adjustment is determined to be effective, otherwise the 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.

[0134] It can be seen that the above embodiment realizes the real-time suppression of the overheating risk of the space monopropellant engine by constructing a closed-loop control system driven by data and model. The residual drive adjustment mechanism adopted in the above embodiment can shorten the response time of temperature anomalies, thereby effectively preventing the occurrence of potential overheating failures, thereby extending the life of key components, reducing the performance degradation rate, and providing technical guarantee for the long-term reliable operation of spacecraft in orbit.

[0135] As an optional but not limiting embodiment, the method further includes: Step 301: According to the performance stability evaluation result, repeatedly appearing failure modes in the historical evolution path of the thermal stress concentration area are screened out.

[0136] In step 301, the performance stability evaluation reports of 300 historical ignition tasks were analyzed, and the common features in the historical evolution path of the thermal stress concentration area were extracted using a pattern recognition algorithm. The spatial coordinates of 45 high-risk parts were grouped using the K-means clustering algorithm, and three main failure modes were identified: the first type was the periodic thermal stress exceeding the standard at the root of the injector swirl plate (occurrence frequency 78%, single exceeding duration 3.2±0.8ms), the second type was the transient temperature spike of the catalyst bed support ring (peak temperature 935K, occurrence probability 62%), and the third type was the accumulated plastic strain at the combustion chamber flange connection (maximum strain 0.23%, growth rate 0.004% / cycle).

[0137] For example, the construction of the failure mode database is based on a MySQL relational data table. Each record contains spatial coordinates, type of exceeded parameter, time of occurrence and duration. Through correlation analysis, it was found that the second type of failure mode was strongly correlated with the propellant mass flow fluctuation (standard deviation > 0.8g / s) (Pearson coefficient 0.76).

[0138] Step 302: Generate a target iteration function for structural optimization of a space monopropellant engine based on the failure mode, wherein the target iteration function includes a thermal stress peak reduction weight and a material quality constraint condition.

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

[0140] Step 303: using a genetic algorithm to perform multiple generations of solving the target iterative function to generate multiple candidate structural improvement schemes that meet the durability requirements; performing a virtual ignition test on the multiple candidate structural improvement schemes through the flow domain grid model, and selecting a target structural improvement scheme whose thermal stress distribution meets the set improvement conditions as an optimization strategy based on the virtual ignition test results.

[0141] For example, the implementation of the genetic algorithm can adopt the NSGA-II multi-objective optimization framework, set the population size to 64, the crossover probability to 0.85, and the mutation probability to 0.15. Each generation of evolution includes the following operations: first, the thermal stress distribution of each individual is calculated by ANSYS Mechanical, and the maximum Von Mises stress value of the injector head is extracted; then Fluent is called to perform transient fluid-thermal coupling simulation to obtain the adjusted temperature field dynamic change sequence; finally, non-dominated sorting and crowding calculation are performed according to the objective function value.

[0142] For another example, after 50 generations of evolution, the Pareto front converged and three candidate solutions were selected: Solution A (support structure diameter 1.2mm, spacing 4mm, inclination 15°) reduced the peak thermal stress by 37.5% and increased the mass by 3.8%; Solution B (diameter 1.5mm, spacing 5mm, inclination 22°) reduced stress by 28.7% and increased mass by 2.1%; Solution C (diameter 0.8mm, spacing 3mm, inclination 8°) reduced stress by 41.2% and increased mass by 4.9%. Virtual ignition tests showed that the temperature peak of Solution C at the catalyst bed support ring dropped from 935K to 862K, and the thermal stress peak dropped from 625MPa to 498MPa, and was selected as the final optimization strategy.

[0143] As an optional but not limiting embodiment, the generation of multiple candidate structural improvement schemes that meet the durability requirements in step 303 includes step 3030: inserting geometric parameter variables of the additional support structure into the pore distribution topology data of the catalytic bed; adjusting the local flow velocity distribution in the porous medium area according to the geometric parameter variables of the additional support structure; recalculating the adjusted temperature field dynamic change sequence to determine whether the reduction in the thermal stress peak reaches a preset threshold; if the reduction in the thermal stress peak reaches the preset threshold, adding the structural improvement features corresponding to the current geometric parameter variables to the candidate scheme set, otherwise continuing the iterative optimization until the termination condition is met.

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

[0145] For example, the evaluation of the reduction in thermal stress peak is based on the average of 10 virtual ignition tests. When the reduction exceeds 30% (such as 41.2% for Scheme C), the current geometric parameter combination is stored in the candidate solution library. Otherwise, the support body inclination angle and diameter combination are adjusted through the orthogonal test method to continue iteration until the maximum number of iterations is 100 or the reduction meets the standard.

[0146] As an expandable embodiment, the propellant decomposition reaction rate sub-model in the combined chemical reaction model in step 1022 spatially couples the catalytic bed pore distribution topology data with the combustion chamber geometry to generate an initial flow domain mesh coupled with multiple physical fields, including: Step 10221: extracting pore space distribution characteristics, connectivity characteristics and morphological attributes from the pore distribution topology data of the catalytic bed, wherein the morphological attributes include the roughness grade of the pore wall and the intersection angle of branch pores; based on the local reaction rate distribution function in the propellant decomposition reaction rate sub-model, identifying the pore area in the pore space distribution characteristics of the catalytic bed that is positively correlated with the propellant decomposition reaction rate, and generating a spatial coordinate set of the reaction rate sensitive area.

[0147] For example, the identification of reaction rate sensitive areas is based on the coupling analysis of the spatial distribution characteristics of the pores of the catalytic bed and the surface reaction dynamics. First, the activity distribution of the iridium catalyst was determined by X-ray photoelectron spectroscopy (XPS), and the density of active sites in the area with pore wall roughness Ra<1.5μm was determined to reach 5.8×10^18sites / m² (23% higher than the rough area). The image processing algorithm was used to extract the connected pore paths with a diameter of 35-50μm (accounting for 68% of the total path length), and the angle between its extension direction and the flow direction of the fluid at the flow channel inlet was calculated, and the paths with angles >45° (accounting for 12%) were eliminated. The target reaction rate sensitive area finally generated contains 327 main channels with a total length of 1.2m, and the spatial coordinate set is imported into the chemical reaction model of ANSYS Fluent in CSV format.

[0148] Step 10222: spatially align the spatial coordinate set of the reaction rate sensitive area with the flow channel inlet position in the combustion chamber geometry, wherein after the spatial alignment is completed, the pore channel outlet of the reaction rate sensitive area coincides with the fluid transmission path of the flow channel inlet; according to the heat release rate gradient in the propellant decomposition reaction rate sub-model, a reaction rate weight coefficient is set 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 area increases.

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

[0150] Step 10223: Based on the reaction rate weight coefficient, the local density of the combustion chamber flow channel grid is adjusted, and an encrypted grid layer is inserted in the area where the weight coefficient is higher than a preset threshold value to generate an initial flow channel grid that matches the pore structure; the pore connectivity characteristics in the pore distribution topology data of the catalytic bed are mapped to the initial flow channel grid to generate a fluid transmission interface between the pore channel and the flow channel grid, and the grid nodes of the interface correspond one-to-one to the grid nodes of the pore channel outlet.

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

[0152] Step 10224: According to the temperature dependence in the propellant decomposition reaction rate sub-model, the local flow velocity distribution in the porous channel and the heat conduction rate in the flow channel grid are associated at the fluid transmission interface to generate a multi-physics data chain of coupled reaction rate and flow heat transfer; the coupling verification of the multi-physics data chain is performed to detect whether the propellant mass flow rate at the porous channel outlet and the received flow rate at the flow channel inlet are conserved. If not, the corresponding relationship of the grid nodes of the fluid transmission interface is readjusted; according to the coupling verification result, the initial flow channel grid, the pore structure grid and the multi-physics data chain are integrated to generate the initial flow domain grid of the multi-physics coupling.

[0153] Exemplarily, the coupling of the multi-physics data chain is realized through Fluent's coupling solver, and a two-way data exchange is set at the fluid transmission interface: the velocity distribution (U, V, W components) in the pore channel is transferred to the flow channel grid as the inlet boundary condition every 0.1ms, and the flow channel wall temperature is fed back to the pore model as the thermal boundary condition of the chemical reaction. The coupling verification uses transient mass conservation monitoring, and the cumulative mass error is 0.023g and the relative error is 0.018% within a 500ms simulation cycle. The final generated initial flow domain mesh contains 11.2 million tetrahedral units in the catalyst bed area, 4.8 million hexahedral units in the flow channel area, and 1.5 million prismatic units in the interface area. The Jacobian factor pass rate is 99.3% and the maximum distortion is 0.82.

[0154] As an expandable embodiment, the identification of pore regions in the pore space distribution characteristics of the catalytic bed that are positively correlated with the propellant decomposition reaction rate in step 10221 includes step 102211: based on the catalyst activity distribution data in the propellant decomposition reaction rate sub-model, screening smooth wall regions with pore wall roughness levels below a preset threshold; in the smooth wall regions, extracting connected paths with pore diameters in a preset target reaction diameter interval as candidate paths for reaction rate sensitive areas; based on the angle between the extension direction of the candidate paths and the fluid flow direction at the inlet of the combustion chamber flow channel, eliminating paths with angles exceeding a preset angle to generate a set of target reaction rate sensitive areas.

[0155] In this embodiment, the screening of smooth wall areas is based on the three-dimensional distribution map of the roughness level of the pore wall, and the threshold value Ra≤1.2μm is set. The continuous smooth area with an area of ​​>0.1mm² is extracted by the regional growing algorithm (a total of 58 areas). The target reaction diameter interval is set to 38±5μm (corresponding to the highest catalytic efficiency interval), 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 by vector dot product, and the paths (accounting for 19%) with an angle of >30° with the mainstream direction (positive direction of the Z axis) of the flow channel inlet are deleted. Finally, 43 optimized paths are retained to constitute the target reaction rate sensitive area set. The spatial distribution data of this set is imported into the CFD pre-processing software in VTK format to guide subsequent meshing and boundary condition setting.

[0156] It can be seen that the above embodiment achieves the coordinated improvement of the structural design of the space monopropellant engine and the reaction flow characteristics through the multi-scale coupling optimization method. Based on the above embodiment, after the optimized engine simulates the ignition cycle, the maximum plastic strain in the thermal stress concentration area of ​​the injector head is significantly reduced, the probability of the temperature spike of the catalyst bed support ring is greatly reduced, and the vacuum specific impulse attenuation rate is improved. Therefore, the establishment of the multi-physics field coupling grid model improves the accuracy of flow heat transfer simulation and provides a high-precision digital platform for the reliability design and life prediction of the aerospace propulsion system.

[0157] In summary, the embodiments of the present application can deeply integrate microstructural 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 effectively improving the defects of insufficient multi-physical field coupling accuracy, lack of dynamic boundary condition characterization, and rough thermal stress evolution analysis.

[0158] Based on the same inventive concept, the present application also provides an engine ignition state analysis system. Figure 2 As shown, it is a structural schematic diagram of a possible engine ignition state analysis system provided in an embodiment of the present application, Figure 2 In the embodiment, the engine ignition state analysis system 200 includes: a processor 210 and a memory 220. The memory 220 stores a computer program executable by the processor 210, and the processor 210 can execute the steps of the above-mentioned spatial monopropellant engine ignition state analysis method by executing the instructions stored in the memory 220.

[0159] 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 enable the engine ignition state analysis system to execute the steps of the above-mentioned spatial monopropellant engine ignition state analysis method. In some possible implementations, various aspects of the spatial monopropellant 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 enable the engine ignition state analysis system to execute the steps in the above-mentioned spatial monopropellant engine ignition state analysis method. For example, the engine ignition state analysis system can execute the following steps: Figure 1 Follow the steps shown in .

Claims

1. A method for analyzing the ignition state of a space monopropellant engine, characterized in that: The method comprises: Obtain the initial temperature parameters, initial pressure parameters and propellant flow parameters of the combustion chamber of the space monopropellant engine under ignition state; Based on the combustion chamber initial temperature parameter, the initial pressure parameter and the propellant flow parameter, a flow domain grid model including porous medium effect characteristics and chemical reaction effect characteristics is generated; Setting dynamic boundary conditions corresponding to the ignition transient process in the flow domain grid model, performing full-cycle numerical simulation on the internal temperature field distribution of the spatial monopropellant engine, and generating a dynamic change sequence of the temperature field under the ignition state; According to the temperature gradient distribution characteristics at different moments in the dynamic change sequence of the temperature field, the thermal stress evolution data of the key components of the space monopropellant engine are extracted, and the performance stability evaluation results under the ignition state are output by combining the thermal stress evolution data with the preset overheating failure threshold.

2. The method according to claim 1, characterized in that The generating of a flow domain grid model including porous medium effect characteristics and chemical reaction effect characteristics based on the combustion chamber initial temperature parameter, the initial pressure parameter and the propellant flow parameter comprises: Based on the porous medium model, a three-dimensional discretization model is performed on the pore structure of the catalyst bed of the spatial monopropellant engine to generate pore distribution topology data of the catalyst bed; Combined with the propellant decomposition reaction rate sub-model in the chemical reaction model, the catalytic bed pore distribution topology data is spatially coupled with the combustion chamber geometric structure to generate an initial flow domain grid for multi-physics field coupling; According to the mass flow distribution in the propellant flow parameters, the initial flow domain grid is locally encrypted to generate a flow domain grid model that meets a preset accuracy condition; The combustion chamber initial temperature parameter and the initial pressure parameter are loaded into the flow domain grid model as initial field variables to complete the initialization of the flow domain grid model.

3. The method according to claim 2, characterized in that The full-cycle numerical simulation of the internal temperature field distribution of the space monopropellant engine is performed to generate a dynamic change sequence of the temperature field under the ignition state, including: Setting the propellant injection time series of the ignition transient process in the flow domain grid model, and calculating the local flow velocity distribution of the porous medium region in each time step based on the pressure fluctuation range in the dynamic boundary condition; Iteratively solving a time evolution model of the temperature field in the combustion chamber according to the local velocity distribution and the heat release rate function in the chemical reaction model to generate snapshots of the temperature field at continuous time stamps; Performing spatial interpolation processing on the temperature field snapshots to generate a temperature field dynamic change sequence covering the entire flow domain of the spatial monopropellant engine; The highest temperature peak value in the catalyst bed outlet area and the occurrence time of the highest temperature peak value in the temperature field dynamic change sequence are extracted, and the highest temperature peak value and the occurrence time are used as key indicators of the temperature gradient distribution characteristics.

4. The method according to claim 3, characterized in that The step of extracting thermal stress evolution data of key components of the space monopropellant engine according to the temperature gradient distribution characteristics at different moments in the temperature field dynamic change sequence comprises: Calculating 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; Based on the instantaneous thermal expansion coefficient and the wall structure constraint conditions, deriving the thermal stress distribution cloud diagram at different times; Generate a historical evolution path of the thermal stress concentration area according to the coordinates of the area exceeding the material yield strength extracted from the thermal stress distribution cloud map and the duration corresponding to the coordinates of the area; The historical evolution path is compared with a preset overheat failure threshold in time series, and the positions of key components with cumulative damage risks are identified based on the time series comparison results.

5. The method according to claim 4, characterized in that Outputting the performance stability evaluation result under the ignition state includes: Calculating the number of thermal fatigue cycles at the key component position according to the historical evolution path of the thermal stress concentration area; Based on the mapping relationship between the number of thermal fatigue cycles and the material durability database, predict the remaining life interval of the space monopropellant engine under continuous ignition conditions; Combining the remaining life span with the key indicators of the temperature gradient distribution characteristics, generating a stability assessment report including overheating risk probability and performance degradation trend; Results in the stability assessment report that exceed a preset risk level are marked as structural weaknesses that need to be optimized.

6. The method according to claim 2, characterized in that The method of performing three-dimensional discretization modeling on the pore structure of the catalytic bed of the spatial monopropellant engine based on the porous medium model to generate catalytic bed pore distribution topology data includes: Acquiring a set of geometric parameters of the catalytic bed, wherein the set of geometric parameters includes a pore diameter distribution range, a pore communication path curvature, and a pore density gradient; Based on the pore diameter distribution range, identifying multi-scale pore morphology characteristic parameters inside the catalyst bed, the multi-scale pore morphology characteristic parameters including the main pore channel axis direction, the branch pore intersection angle and the pore wall roughness; According to the axis direction of the main pore channel and the pore density gradient, non-uniform grid units are divided, the size of the non-uniform grid units is adaptively adjusted along the curvature of the pore communication path, and the grid nodes are aligned with the intersection angles of the branch pores; Performing topological consistency verification on the non-uniform grid units, determining the continuity of pore wall roughness data of adjacent grid units and the fracture state of pore connection paths, and obtaining a topological consistency verification result; According to the topology consistency verification result, the grid units that failed the verification are locally reconstructed, the grid node positions are adjusted to match the branch pore intersection angles, and an optimized porous medium grid model is generated; Extracting pore morphological characteristic parameters, grid size and node coordinates of each grid unit in the porous medium grid model, and generating catalytic bed pore distribution topology data including pore spatial distribution, connectivity characteristics and morphological attributes; The method for identifying multi-scale pore morphological characteristic parameters inside the catalytic bed includes: dividing the pores of the catalytic bed into main pore channels and secondary branch pores based on the pore diameter distribution range, wherein the diameter of the main pore channel is greater than a preset threshold and extends axially; extracting the centerline trajectory of the main pore channel, calculating the curvature change rate and the extension direction deflection angle of the centerline trajectory as the core parameters of the axial direction of the main pore channel; in the secondary branch pores, detecting the entrance position connected to the main pore channel, measuring the angle between the branch pore at the entrance and the main pore channel, and generating a set of branch pore intersection angles; quantifying the microscopic undulation characteristics of the pore wall through surface scanning data, and generating a spatial distribution map of the pore wall roughness; The non-uniform grid unit division includes: dynamically adjusting the grid unit length according to the curvature change rate along the axis direction of the main pore channel, wherein the curvature change rate is negatively correlated with the grid unit length; inserting encrypted grid nodes at the intersection angle of the branch pores so that the node density is inversely proportional to the intersection angle of the branch pores; setting tetrahedral grid units in the pore-dense area and setting hexahedral grid units in the pore-sparse area based on the pore density gradient; mapping the spatial distribution map of the pore wall roughness to the grid unit surface to generate a wall boundary layer grid with roughness attributes; The topological consistency verification of the non-uniform grid unit also includes: determining whether there is a mutation in the pore wall roughness data of adjacent grid units at the junction, and marking the unit as inconsistent if the mutation amplitude exceeds a preset tolerance; verifying the grid unit continuity of the main pore channel, and the grid unit continuity characterizes the smooth transition of the curvature change rate of the centerline trajectory between adjacent units; detecting whether the grid nodes at the intersection angle of the branch pores are aligned, and calculating the node position correction if there is an offset; and generating a topological consistency verification report including inconsistent unit positions, mutation amplitudes and node offsets.

7. The method according to claim 1, characterized in that The method further comprises: Collecting measured temperature data of the space monopropellant engine in historical ignition missions and performance degradation records corresponding to the measured temperature data; Using the measured temperature data to perform error correction on the temperature field dynamic change sequence to generate a calibrated temperature prediction model; Jointly training the calibrated temperature prediction model with the thermal stress evolution data to generate an abnormal temperature early 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, the propellant flow rate is dynamically adjusted to suppress the overheating risk; The generating of the abnormal temperature early warning model for real-time monitoring comprises: Extracting the residual distribution between the temperature field dynamic change sequence and the measured temperature data from the historical ignition tasks; Using the residual distribution to train a convolutional neural network model to generate a spatiotemporal feature extractor for temperature prediction residuals; The spatiotemporal feature extractor is embedded into the abnormal temperature early 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 a preset safety margin, a propellant flow adjustment instruction is activated to balance the thermal load of the combustion chamber.

8. The method according to claim 7, characterized in that The method of dynamically adjusting the propellant flow rate to suppress the overheating risk based on the deviation signal output by the abnormal temperature warning model includes: Combining the spatial distribution of the local difference coefficient and the deviation signal, identifying a target area with abnormally elevated temperature; Calculating the flow valve opening correction amount to be adjusted based on the position information of the target area and the topological relationship of the propellant injection pipeline; The flow valve opening correction amount is sent to the engine control system, so that the propellant distribution and the current heat load demand are re-matched through the engine control system to obtain an adjusted temperature field dynamic change sequence; The adjusted temperature field dynamic change sequence is monitored and verified, and the monitoring and verification is completed when the local difference coefficient falls back to a safe tolerance range.

9. The method according to claim 1, characterized in that: The method further comprises: Based on the performance stability evaluation results, recurring failure modes in the historical evolution path of the thermal stress concentration area are screened out; generating a target iteration function for structural optimization of a space monopropellant engine based on the failure mode, wherein the target iteration function includes a thermal stress peak reduction weight and a material quality constraint condition; Using a genetic algorithm to solve the target iterative function for multiple generations, generating multiple candidate structural improvement solutions that meet durability requirements; Performing a virtual ignition test on the plurality of candidate structural improvement schemes through the flow domain grid model, and selecting a target structural improvement scheme whose thermal stress distribution meets the set improvement conditions as an optimization strategy according to the virtual ignition test results; The generating of multiple candidate structural improvement solutions that meet the durability requirements includes: Inserting geometric parameter variables of additional support structures into the pore distribution topology data of 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; Recalculate the adjusted temperature field dynamic change sequence to determine whether the reduction in the thermal stress peak value reaches a preset threshold; If the reduction amplitude of the thermal stress peak reaches the preset threshold, the structural improvement feature corresponding to the current geometric parameter variable is added to the candidate solution set, otherwise the iterative optimization is continued until the termination condition is met.

10. An engine ignition state analysis system, characterized in that: It comprises a processor and a memory, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of any one of the methods of claims 1 to 9.

Images