A fluid equivalent method for calculating complex two-phase flow phase change processes

Abstract

This invention discloses a fluid equivalent method for calculating complex two-phase flow phase change processes. Addressing the challenges of Jacobian matrix singularity and computational divergence arising from the fully coupled solution of multiphase flow models and the kinetic source terms of strongly endothermic / exothermic reactions, this method isolates the computational domain of the real gas-liquid two-phase flow and replaces it with an equivalent single-phase virtual fluid. Based on the latent heat characteristics of the fluid phase change, a specific heat amplification factor M is dynamically set, and a maximum equivalent thermal conductivity is set in conjunction with the system's series thermal resistance model, accurately replicating the isothermal source boundary properties of the gas-liquid phase change process at the physical level. Furthermore, the reduced-dimensional equivalent heat transfer boundary conditions are extracted and mapped to the inner tube of the reaction, completing multi-field coupled solutions and forward optimization of the flow channel configuration. This invention completely eliminates the numerical rigidity of complex phase change calculations, achieving unconditional stability of the underlying algorithm and providing accurate simulation technology support for the design of high-efficiency phase change heat transfer reactors and internal enhanced heat transfer structures.

Description

Technical Field

[0001] This invention relates to the fields of computational fluid dynamics and thermal engineering, specifically to a fluid equivalent simulation method for complex two-phase flow phase change processes, and a phase change heating reactor design method based on the equivalent boundary conditions for forward structural optimization. Background Technology

[0002] In advanced aerospace propulsion systems, particularly in the highly efficient hybrid architecture of integrated solid oxide fuel cells (SOFCs), the front-end aviation kerosene endothermic reforming reactor is the core fuel conversion unit that determines the overall power-to-weight ratio and thermal management efficiency. To meet extremely high endothermic demands, the reactor often adopts a shell-and-tube structure: the outer tube uses the latent heat of phase change of gaseous metals (such as sodium heat pipes) for heating, while the inner tube is the catalytic reforming endothermic zone. Performing three-dimensional numerical simulations of such multi-field coupled systems to guide the design of the physical structure presents significant challenges in underlying algorithms. Existing technologies in fluid dynamics calculation software often directly call multiphase flow models (such as VOF or Mixture models) to calculate the external phase change and solve the internal reforming reaction kinetic source terms in a fully coupled manner. This computational framework suffers from a fatal numerical rigidity defect: the system internally deploys a rigid set of partial differential equations with strong nonlinearity and multiple physics fields (phase change latent heat field, component transport field, and strong endothermic reaction heat field). This coupling causes a sharp deterioration in the Jacobian matrix condition number when the solver simultaneously updates the latent heat of vaporization and component mass source terms, leading to residual inflation and floating-point anomalies, ultimately causing the energy equation to diverge. Computational instability severely restricts the efficient optimization of flow channel configurations within complex reactors. Summary of the Invention

[0003] The purpose of this invention is to overcome the numerical divergence defects of existing fully coupled calculation methods when dealing with high-temperature phase transitions and strong endothermic reactions, and to provide a fluid equivalent simulation method with stable calculation and accurate boundary mapping. Based on this method, a reactor entity structure that balances high fuel conversion rate and low fluid pressure drop is generated.

[0004] To achieve the above objectives, this invention provides a simulation and structural optimization design method for phase change heating reactors based on fluid equivalence, comprising the following steps:

[0005] Step 1: Decouple the model establishment from the computational domain and establish an initial geometric model that includes the endothermic reaction inner tube, the solid tube wall, and the phase change heating outer tube; strip the actual gas-liquid two-phase flow computational domain inside the phase change heating outer tube and replace it with an equivalent single-phase virtual fluid computational domain;

[0006] Step 2: Equivalent reconstruction of physical property parameters. Set the equivalent specific heat capacity of the single-phase virtual fluid. ,in The specific heat amplification factor M is set based on the ratio of the latent heat of phase change of the fluid to the specific heat capacity at constant pressure in the gaseous state, and the equivalent thermal conductivity is set. Greater than a first preset threshold, to characterize the limiting convective heat transfer state of the fluid condensation phase change process;

[0007] Step 3: Dimensional reduction boundary condition extraction. Perform fully coupled conjugate heat transfer calculations on the single-phase virtual fluid computational domain, solid pipe wall and inner pipe single-phase fluid. After convergence, extract the temperature distribution matrix of the outer wall of the inner pipe and use it as the equivalent heat transfer boundary condition.

[0008] Preferably, the conjugate heat transfer solution domain in step 3 includes only the single-phase virtual fluid computation domain, the solid tube wall domain coupled with it, and the inner tube single-phase fluid domain, excluding the simulation of the reforming endothermic reaction in the inner tube, so as to achieve computational decoupling between the phase change heating side boundary extraction process and the inner tube side reaction solution process.

[0009] Furthermore, the equivalent heat transfer boundary condition extracted in step 3 is preferably the interface temperature distribution on the solid tube wall that corresponds one-to-one with the outer surface of the heat-absorbing reaction inner tube in step 4 at the axial and circumferential positions. The temperature distribution can be represented as a discrete temperature matrix, a set of nodal temperatures, or a continuous temperature field after interpolation reconstruction.

[0010] Extract the interface temperature distribution on the solid tube wall that corresponds geometrically to the outer surface of the endothermic reaction inner tube in step 4, and map it onto the outer surface of the endothermic reaction inner tube.

[0011] Step 4: Multi-field coupling solution and forward structure optimization. The equivalent heat transfer boundary conditions are mapped to the outer surface of the endothermic reaction inner tube. An endothermic chemical kinetic model is loaded in the inner tube to solve the flow field and concentration field. Based on the component concentration or temperature field distribution characteristics of the target product in the tube, the forward optimization design of the internal flow channel configuration or heat transfer enhancement structure of the reactor is carried out based on the equivalent heat transfer boundary conditions until the target fluid parameters at the reactor outlet reach the preset threshold.

[0012] Preferably, the endothermic reaction inner tube model in step 4 includes the inner tube solid wall and the reaction fluid domain inside the inner tube. The equivalent heat transfer boundary conditions extracted in step 3 are applied to the outer surface of the inner tube solid wall in a one-way mapping manner, serving as the external thermal boundary conditions for the coupled solution of flow, heat transfer, component transport and chemical reaction kinetics on the inner tube side.

[0013] In a preferred embodiment, the temperature field, heat flux density, or reaction heat effect obtained from the inner tube side are not fed back to the single-phase virtual fluid computational domain described in step 3 for iterative updates, thereby maintaining the computational decoupling between the phase change heating side and the reaction side and avoiding the reintroduction of rigid closed-loop coupling between the real two-phase flow and the strong reaction source term.

[0014] The true gas-liquid two-phase flow equivalently substituted inside the phase change heat supply outer pipe covers low-temperature, medium-temperature or high-temperature phase change heat transfer working fluids; including but not limited to the following fluid working fluids, low-temperature phase change heat transfer working fluids: conventional refrigerants (such as R134a, R410A, ammonia), cryogenic fluids (such as liquid nitrogen , liquid oxygen, liquid hydrogen); medium-temperature phase change heat transfer working fluids: water / steam, heat-conducting oil (such as some boilable phase change diphenyl ether-based Therminol biphenyl mixtures), alcohols, etc.; high-temperature phase change heat transfer working fluids: gaseous sodium (Na), gaseous potassium (K), gaseous lithium (Li), etc. The specific heat amplification coefficient M is valued according to the thermodynamic characteristics of the specific heat transfer working fluid, and its theoretical calculation satisfies the formula:

[0015] In the formula, is the molar phase change latent heat of the heat transfer working fluid, is the molar mass of this working fluid, is the equivalent minimum temperature gradient, and the value range is set to be from 0.1 K to 2 K; is the specific heat at constant pressure of the gaseous state of the working fluid. The specific value range includes: (1) When the heat supply working fluid is a low-temperature phase change fluid and the applicable temperature range is from -200 °C to 100 °C, the specific heat amplification coefficient is set to 50 < M < 2000; (2) When the heat supply working fluid is a medium-temperature phase change fluid and the applicable temperature range is from 100 °C to 400 °C, the specific heat amplification coefficient is set to 200 < M < 5000; (3) When the heat supply working fluid is a high-temperature phase change fluid and the applicable temperature range is from 600 °C to 1200 °C, the specific heat amplification coefficient is set to 500 < M < 10000.

[0016] For the high-temperature phase change fluid, the determination condition for the lower bound of the value of the specific heat amplification coefficient M is: Based on the energy equivalence principle, when the equivalent thermal inertia corresponding to the M value is not sufficient to absorb the temperature剧变generated by the internal endothermic reaction source term, the axial temperature drop gradient of the single-phase virtual fluid in the flow direction exceeds the latent heat constant temperature physical threshold, resulting in the failure of the equivalent heat transfer boundary condition; the lower bound of the value of the equivalent thermal conductivity is determined based on the series heat transfer resistance model of the casing reactor system; the series heat transfer resistance model satisfies the formula:

[0017]

[0018] In the formula, K is the total heat transfer coefficient of the system, is the external convective heat transfer resistance on the side of the single-phase virtual fluid, is the thermal conductivity resistance of the solid pipe wall, is the internal convective heat transfer resistance on the side of the endothermic reaction inner pipe; for the high-temperature phase change fluid, when the specific heat amplification coefficient is set to M > 500, the equivalent thermal conductivity The value of must be greater than 1000 W / (m·K) to forcefully smooth the spatial distribution of the temperature gradient through the algorithm, so that the external convective heat transfer thermal resistance in the formula approaches zero, thereby replicating the latent heat constant temperature boundary at the physical level.

[0019] The endothermic chemical reaction is the reaction of aviation kerosene reforming to produce hydrogen, and the heat supply working medium is gaseous sodium; the target product is hydrogen.

[0020] The specific numerical discretization and control methods for the multi-field coupling solution in step 4 include: using the laminar finite rate model to calculate the kinetic source term of the reforming reaction; enabling the pseudo-transient solution framework to suppress the singularity of the Jacobian matrix caused by the strong endothermic reaction, where the time scale control method is set to automatic and the length scale method is set to conservative; the spatial discretization of the momentum equation, energy equation, and component transport equation all adopts the second-order upwind scheme.

[0021] The beneficial effects of the present invention are as follows: First, eliminate numerical stiffness and achieve unconditional stability of the algorithm. The present invention constructs a fluid equivalent model based on the mechanism of heat transfer thermal resistance series, and limits the specific heat amplification coefficient within a reasonable range (such as setting it in the range of 500 < M < 10000 for high-temperature working fluids). This parameter domain not only avoids the singularity of the Jacobian matrix caused by the phase interface tracking equation but also ensures that there is no significant sensible heat temperature drop gradient change in the virtual fluid; combined with the extreme value setting effectively eliminates the external convective heat transfer thermal resistance and accurately replicates the constant temperature source boundary property during the latent heat release process of gaseous metal liquefaction at the physical level. Further combined with the pseudo-transient solution framework and implicit relaxation mechanism, the divergence problem of fuel reforming multiphase flow coupling simulation is completely overcome. Second, break through the conversion rate bottleneck and balance the low-pressure drop requirement at the system level. Based on the spiral finned inner tube structure designed with equivalent heat transfer boundary conditions, the three-dimensional spiral flow disturbance effect is used to destroy the near-wall thermal boundary layer. Data shows that under the constraint of an equivalent constant wall temperature of 1180 K, the hydrogen mole fraction at the outlet of the design group jumps from 0.3 of the control group to 0.5, and the conversion rate increases by more than 60%. Description of the Drawings

[0022] Figure 1 It is a curve graph of the temperature distribution on the outer wall surface of the inner tube under the equivalent heat transfer boundary conditions extracted in the embodiment of the present invention.

[0023] Figure 2 It is a global temperature distribution nephogram of the conjugate heat transfer model established in the embodiment of the present invention.

[0024] Figure 3 It is a nephogram of the internal temperature distribution of the design group model (with spiral fins) in the embodiment of the present invention.

[0025] Figure 4This is a cloud map showing the internal temperature distribution of the control group model (finless light tube) in the comparative example.

[0026] Figure 5 To design a cloud map of the hydrogen mole fraction distribution within the model group.

[0027] Figure 6 This is a cloud map showing the hydrogen mole fraction distribution within the control group model.

[0028] Figure 7 The hydrogen mole fraction curve along the axis of the design group model is shown.

[0029] Figure 8 The graph shows the hydrogen mole fraction along the axis of the control group model.

[0030] Figure 9 To design the reactants within the group model ( Mole fraction distribution cloud map.

[0031] Figure 10 The reactants in the control model ( Mole fraction distribution cloud map. Detailed Implementation

[0032] This invention discloses a fluid equivalent method for calculating complex two-phase flow phase change processes. Addressing the challenges of Jacobian matrix singularity and computational divergence arising from the fully coupled solution of multiphase flow models and the kinetic source terms of strongly endothermic / exothermic reactions, this method isolates the computational domain of the real gas-liquid two-phase flow and replaces it with an equivalent single-phase virtual fluid. Based on the latent heat characteristics of the fluid phase change, a specific heat amplification factor M is dynamically set, and a maximum equivalent thermal conductivity is set in conjunction with the system's series thermal resistance model, accurately replicating the isothermal source boundary properties of the gas-liquid phase change process at the physical level. Furthermore, the reduced-dimensional equivalent heat transfer boundary conditions are extracted and mapped to the inner tube of the reaction, completing multi-field coupled solutions and forward optimization of the flow channel configuration. This invention completely eliminates the numerical rigidity of complex phase change calculations, achieving unconditional stability of the underlying algorithm and providing accurate simulation technology support for the design of high-efficiency phase change heat transfer reactors and internal enhanced heat transfer structures.

[0033] The present invention will be further described in detail below with reference to the accompanying drawings. It should be understood that this example is used to illustrate the physical logic and computational framework of the present invention and does not constitute a limitation on the scope of protection of the claims.

[0034] The main technical solutions in this patent are as follows:

[0035] Step 1: Model Establishment and Decoupling of Computational Domain: Establish an initial geometric model that includes the endothermic reaction inner tube, solid tube wall, and phase change heating outer tube; strip away the actual gas-liquid two-phase flow computational domain inside the phase change heating outer tube and replace it with an equivalent single-phase virtual fluid computational domain;

[0036] Step 2 Physical property parameter equivalent reconstruction: Set the equivalent specific heat capacity of the single-phase virtual fluid , and its core idea is to equivalently convert the latent heat near the phase change temperature into sensible heat within an extremely small temperature gradient t. Dynamically set the value range of the specific heat amplification coefficient M to be limited within (such as the interval of 50 < M < 10000); Set the equivalent thermal conductivity to be greater than the first preset threshold value;

[0037] Step 3 Dimensionality reduction boundary condition extraction: Conduct conjugate heat transfer calculations on the single-phase virtual fluid and the solid tube wall, extract the steady-state temperature distribution matrix of the inner wall surface of the solid tube, and use it as the equivalent heat transfer boundary condition;

[0038] Step 4 Multi-field coupling solution and forward structure optimization: Map the equivalent heat transfer boundary condition to the outer surface of the inner tube of the endothermic reaction, and load the chemical reaction kinetics model; According to the component concentration or temperature field distribution characteristics of the target product in the tube, perform forward optimization design of the internal flow channel configuration of the reactor based on the equivalent heat transfer boundary condition until indicators such as the conversion rate of the reactor outlet product reach the standard.

[0039] As a preferred technical solution, the specific numerical discretization and control method for the multi-field coupling solution includes: Using the laminar finite rate model to calculate the kinetic source term; Enabling the pseudo-transient solution framework (with an automatic time scale and a conservative length scale) to suppress rigid divergence; All physical quantity transport equations are spatially discretized using the second-order upwind scheme.

[0040] Specifically:

[0041] Core mechanism deduction of the embodiment: According to the heat transfer principle, the total heat transfer thermal resistance of the double-pipe reactor system can be systematically expressed as the following formula . Among them, the external convective heat transfer thermal resistance 1 / h. Due to the extremely large convective heat transfer coefficient of the condensation phase change of gaseous sodium vapor, tends to zero. Based on this, it is inferred that the temperatures of the inner and outer tubes should tend to be the same as the external steam temperature.

[0042] Example Verification Step 1: Extraction and Verification of Dimensional Reduction Boundary Conditions. A steady-state conjugate heat transfer model is constructed, consisting of a "virtual single-phase fluid domain + solid tube wall domain + inner tube single-phase fluid domain". In this example, the solution model in the dimensionality reduction boundary extraction stage only retains the outer single-phase virtual fluid domain, the solid tube wall domain in direct contact with it, and the inner tube single-phase fluid model. The endothermic chemical reaction kinetic source term of the inner tube reaction fluid domain is not introduced in this stage. After convergence, the interface temperature field data on the solid tube wall corresponding to the geometry of the outer surface of the subsequent endothermic reaction inner tube model are exported, and a temperature distribution matrix is ​​formed according to axial and circumferential coordinates, serving as the equivalent heat transfer boundary condition applied in subsequent steps. The specific heat amplification factor M is set to 1000. The model is solved, and the temperature field data of the outer wall of the inner tube is extracted. Combined with... Figure 1 and Figure 2 It can be seen that the outer wall of the inner tube remains stable at around 1180K, with minimal temperature gradient along the axial direction. This result verifies that when M>500, the fluid equivalent parameters successfully replicate the isothermal physical characteristics of the real latent heat environment of phase change.

[0043] Example Verification Step Two: Forward Structural Design Based on Equivalent Boundaries. As a specific embodiment of this equivalent method in the design of a reforming reactor for an aero-engine system, two sets of reactor inner tube models were established for comparison and verification. The control group model was a bare tube without fins on the inner wall, with dimensions of 120mm length, 17mm outer diameter, and 15mm inner diameter. The design group model, based on the hydrogen mole fraction distribution characteristics inside the tube, underwent forward optimization based on the aforementioned equivalent heat transfer boundary conditions. Its specific enhanced heat transfer structure was designed by adding helical fins to the inner wall of the control group model, with specific geometric parameters of 20mm pitch, 3.5 continuous turns, and 3mm fin height. To ensure stable flow field at the inlet and outlet, straight tube buffer zones of 15mm and 35mm were reserved at the front and rear, respectively.

[0044] Example Verification Step 3: Solver Control and Discretization Settings. To ensure accurate coupling between the complex component transport, the strong endothermic reaction source term, and the helical turbulence characteristics within the inner tube, rigid control rules are established: a coupled pressure-velocity solver is used. The component transport equations are coupled with a laminar finite-rate model. For the rigid partial differential equations induced by this model, a pseudo-transient solution framework is enabled, with the time scale set to automatic and the length scale set to conservative. An implicit pseudo-time step is introduced to suppress residual divergence. The spatial discretization of momentum, energy, and each component transport equation is uniformly set to a second-order upwind scheme to accurately capture abrupt changes in the cross-sectional concentration field.

[0045] Within the pseudo-transient solution framework, the fluid equivalent simulation method of this invention strictly adheres to the following set of partial differential governing equations reconstructed through parametric equivalence within the computational domain. By discretizing and solving the following equations, the rigidity of the multi-field coupled system is eliminated:

[0046] (1) Continuity equation (mass conservation):

[0047]

[0048] In the formula, ρ is the fluid density, v is the velocity vector, and t is time.

[0049] (2) Momentum conservation equation

[0050]

[0051] In the formula, p is the static pressure, ρg and F are the gravitational body force and the external body force, respectively; As an effective stress tensor, it includes not only the viscous stress caused by molecular viscosity, but also the Reynolds stress term, which is used to precisely close the three-dimensional swirling flow induced by the helical fins and the boundary layer separation.

[0052] (3) Energy conservation equation for equivalent reconstruction:

[0053]

[0054] In the formula, it can always .

[0055] This invention, through equivalent reconstruction of physical property parameters, provides strong thermal inertial damping in the transient and convection terms on the left side of the energy equation with a very large specific heat capacity, absorbing / releasing heat from the source term without producing a significant temperature gradient (T); simultaneously, the heat conduction term is endowed with an extreme equivalent thermal conductivity. This forcibly smooths out the spatial distribution of the temperature gradient.

[0056] (4) Component transport and chemical reaction kinetic equations:

[0057]

[0058] In the formula, Let i be the mass fraction of the i-th component. Diffusion flux. Net formation rate of the reaction. Calculated directly from the laminar finite rate model, it follows Arrhenius's law:

[0059]

[0060] In the prior art, the Arrhenius index is extremely likely to cause the divergence of the component equation residuals. In the present invention, it is strongly coupled with an equivalent energy equation with extremely large thermal inertia damping, and the source term is realized within the pseudo-time step and the stable iteration of the latent heat constant temperature boundary, solving the divergence problem of the fully coupled calculation from the underlying partial differential equation level.

[0061] Step 4 of the embodiment verification: Multi-field coupling solution and result comparison.

[0062] The extracted equivalent wall temperature distribution of about 1180 K is unidirectionally loaded onto the outer surfaces of two sets of endothermic reaction inner tube models according to the spatial correspondence relationship. In this embodiment, both sets of endothermic reaction inner tube models include an inner tube solid wall and an inner tube inner reaction fluid domain, and the equivalent wall temperature distribution serves as the thermal boundary condition for the outer surface of the inner tube solid wall; subsequently, the coupled processes of inner tube wall heat conduction, internal flow, component transport, and reforming reaction are solved respectively. To maintain the decoupling relationship between the outer side phase change heat supply solution and the inner side reaction solution, the heat flux information and temperature information obtained from the inner tube side are not reversely transmitted back to the aforementioned virtual single-phase fluid domain for re-iteration. The inlet flow velocity is 0.01 m / s, the temperature is 776 K, and the gas phase molar fraction ratio : =0.027:0.973.

[0063] Comparison Figure 3 (design group) and Figure 4 (control group), the radial temperature gradient of the internal fluid in the design group is significantly reduced, and the spiral turbulent flow forces the high-temperature fluid to mix into the cold source central area. Figures 5 to 8 It shows that when the fluid flows through the fin area of the design group, the hydrogen molar fraction climbs non-linearly and accelerates. Figure 9 And Figure 10 It shows that the design group achieved the complete decomposition of the core component within an extremely short axial distance. At the final outlet section, the hydrogen molar fraction of the control group was about 0.3, and that of the design group reached about 0.5. Through the physical layout guided by the algorithm of the present invention, the conversion rate of the reforming reaction achieved an engineering benefit jump of more than 60%.

[0064] To verify the physical accuracy of the fluid equivalence method described in the present invention, parameter reconstruction calculations were carried out for typical working fluids in different temperature ranges. The equivalent temperature gradient t = 1 K was set, and the obtained equivalent physical properties are shown in Table 1. For metallic sodium commonly used in aeroengine thermal management, the calculated equivalent specific heat capacity is about 3871304 J / (kg·K), and the specific heat amplification coefficient M is about 4300 (calculated based on ≈900 J / (kg·K)), which falls within the protection range of 500 < M < 10000 set in claim 2 of the present invention, proving the physical accuracy of the equivalence method and the protection numerical range.

[0065] Table 1. Partial Equivalent Physical Property Parameters

[0066] Material molecular weight Phase transition temperature (K) Latent heat of phase transition (kJ / mol) Equivalent specific heat (J / (Kg·K)) sodium 23 1154.55 89.04 3871304.348 Potassium 39 1029.65 77.33 1982820.513 water 18 373.2 39.5 2194444.444 Ethylene glycol 62 470.5 52.49 846612.9032 Propylene glycol 76 487.6 57.86 761315.7895 ammonia 17 239.7 22.77 1339411.765 nitrogen 28 77.35 5.68 202857.1429 oxygen 32 90.17 6.74 210625 hydrogen 2 20.39 0.46 230000

Images