Phase change sweating cooling non-equilibrium flow heat transfer numerical simulation method

Abstract

The invention discloses a phase change sweating cooling non-equilibrium flow heat transfer numerical simulation method which comprises the following steps: after initialization assignment, obtaining physical property parameters of a liquid phase and a gas phase according to a fluid temperature, and calculating a mixing density, a hybrid viscosity, a mixing enthalpy, a convection term coefficient and a diffusion term coefficient; boundary conditions are given to the cooling fluid continuity equation, the cooling fluid momentum equation, the cooling fluid energy equation and the porous framework energy equation; solving a cooling fluid continuity equation and a cooling fluid momentum equation; various items in a cooling fluid energy equation and a porous skeleton energy equation are divided into an unsteady state item, a convection item, a diffusion item and a source item, the four items are expressed by compiling a UDF function and then loaded into a UDS equation in FLUENT, new mixing speed, new mixing pressure intensity, new mixing enthalpy and new solid phase temperature are solved, and the optimal mixing speed, the optimal mixing enthalpy and the optimal solid phase temperature are obtained. The difference value between the initial value and the corresponding initial value in the iteration is judged, if the requirement is met, calculation can be ended, and if not, iteration is conducted again after the solid phase temperature and the fluid saturation are updated.

Description

Technical Field

[0001] The invention relates to a phase change sweat cooling non-equilibrium flow heat transfer numerical simulation method, belonging to the technical field of aircraft thermal protection. Background Art

[0002] Aerospace technology is at the forefront of modern scientific and technological development. In addition to meeting the requirements of flying higher and faster, saving operating costs and the reusability of carriers have also become the inevitable development trend of future aerospace vehicles. Along with this, a series of key technical issues have come along, one of which is the thermal protection problem. The development speed of aerospace vehicles has exceeded the development speed of high-temperature resistant materials. Simply relying on the development of high-temperature resistant materials can no longer meet the design requirements of future aerospace flights. Active cooling technology must be combined to provide higher thermal protection capabilities. Transpiration cooling has excellent thermal protection effects and helps to reduce wall friction resistance, so it is considered to be the most promising cooling technology to solve the thermal protection problem of high-speed aircraft. Since the experimental research on transpiration cooling has high requirements for equipment and it is difficult to observe the internal flow field, numerical simulation of the transpiration cooling process is an important research method. At present, the numerical simulation method of single-phase flow heat transfer for transpiration cooling is relatively mature, while the numerical simulation method for phase change transpiration cooling is relatively small, especially the simulation of the non-equilibrium flow heat transfer process of phase change transpiration cooling cannot be directly completed by existing commercial software. Summary of the invention

[0003] The technical problem solved by the present invention is: to overcome the shortcomings of the prior art, to provide a numerical simulation method for non-equilibrium flow heat transfer of phase change sweating cooling, to carry out secondary development based on the FLUENT platform, to achieve accurate and rapid prediction of non-equilibrium flow heat transfer performance of phase change sweating cooling, and to provide support for the engineering design of phase change sweating cooling.

[0004] The technical solution of the present invention is: a phase change sweat cooling non-equilibrium flow heat transfer numerical simulation method, comprising the following steps:

[0005] Step 1: The mixing velocity u in the cooling fluid continuity equation, the mixing pressure p in the cooling fluid momentum equation, the mixing enthalpy h in the cooling fluid energy equation, and the solid phase temperature T in the porous skeleton energy equation are s Perform initialization assignment;

[0006] Step 2: Fluid saturation s and fluid temperature T f Perform initialization assignment;

[0007] Step 3: According to the fluid temperature T f Obtain the physical properties of liquid and gas phases;

[0008] Step 4. According to the liquid and gas phase physical parameters and fluid saturation s obtained in the previous step, calculate the mixed density ρ, hybrid viscosity μ, convection term coefficient γ and diffusion term coefficient Γ h ;

[0009] Step 5. Assign boundary conditions to the cooling fluid continuity equation, cooling fluid momentum equation, cooling fluid energy equation, and porous skeleton energy equation;

[0010] Step 6. Solve the cooling fluid continuity equation and the cooling fluid momentum equation;

[0011] The cooling fluid energy equation and the porous skeleton energy equation are divided into unsteady-state terms, convection terms, diffusion terms and source terms. The above four terms are expressed by writing UDF functions and loaded into the UDS equation in FLUENT. The cooling fluid energy equation and the porous skeleton energy are solved respectively to obtain the new mixing velocity u new , the new mixing pressure p new 、New mixing enthalpy h new , the new solid phase temperature T s-new ;

[0012] Step 7, determine the new u new 、p new 、h new 、T s-new Compared with the initial u, p, h, T in this iteration s If the difference between the two meets the requirements, the calculation can be ended, otherwise u new 、p new 、h new 、T s-new Assign values ​​to u, p, h, and T respectively s , according to h and T f , the relationship between s is updated T f , s, return to Step3 and iterate again.

[0013] Preferably, a function is written through a UDF program to solve the problem, specifically:

[0014] The initialization function is used to initialize the saturation s and fluid temperature T f ;

[0015] The characteristic definition function is used to define and calculate the fluid mixture density ρ, the hybrid viscosity μ, the convection term coefficient γ and the diffusion term coefficient Γ h ;

[0016] The parameter adjustment function is used to adjust the fluid saturation s and temperature T after each iteration. f , mixture density ρ, hybrid viscosity μ, convection term coefficient γ and diffusion term coefficient Γ h;

[0017] The unsteady-state term definition function defines and solves the unsteady-state terms of the fluid and porous framework energy equations;

[0018] The convection term definition function is used to define and solve the convection term of the fluid energy equation;

[0019] The diffusion term definition function is used to define and solve the diffusion term of the fluid energy equation;

[0020] The heat exchange definition function is used to define and solve the heat exchange between the fluid and the porous framework;

[0021] The boundary definition function is used to define the boundary conditions of the porous skeleton energy equation and the cooling fluid energy equation.

[0022] Preferably, according to the fluid temperature T f The physical properties of the liquid and gas phases are solved by looking up tables or fitting formulas. The physical properties of the liquid and gas phases include: liquid phase density ρ l , gas phase density ρ v , liquid phase dynamic viscosity μ l , gas phase dynamic viscosity μ v , liquid phase saturation enthalpy h l , gas phase saturation enthalpy h v , latent heat of vaporization h fg , liquid specific heat c l Relative heat of gas c v .

[0023] Preferably, the mixed density ρ, the mixed viscosity μ, the convection term coefficient γ and the diffusion term coefficient Γ h The calculation methods are:

[0024] ρ=sρ l +(1-s)ρ v

[0025]

[0026]

[0027]

[0028] Where s is the liquid saturation, ρ l is the liquid density, ρ v is the gas phase density, k rl is the liquid phase relative permeability, μ l is the liquid phase dynamic viscosity, k rv is the gas phase relative permeability, μ v is the gas phase dynamic viscosity, λ is the liquid phase relative fluidity, h l is the liquid phase saturation enthalpy, hv is the gas phase saturation enthalpy, k f,eff is the effective thermal conductivity of the fluid, T f is the fluid temperature, h is the mixing enthalpy, h fg is the latent heat of vaporization, and D(s) is the capillary diffusion coefficient.

[0029] Preferably, boundary conditions are assigned to the cooling fluid continuity equation, the cooling fluid momentum equation, the cooling fluid energy equation, and the porous skeleton energy equation, specifically:

[0030] The velocity and temperature values ​​are given at the cooling fluid inlet, the pressure value is given at the cooling fluid outlet, the cold end of the porous skeleton is adiabatic, and the heat flux density value is given at the hot end of the porous skeleton.

[0031] Preferably, the cooling fluid energy equation is:

[0032]

[0033] in, is the unsteady-state term of the cooling fluid energy equation, is the convection term of the cooling fluid energy equation, is the diffusion term of the cooling fluid energy equation, q sf is the source term, which represents the fluid-solid heat transfer;

[0034] ε is the porosity, t is the time, ρ is the mixing density, h is the mixing enthalpy, γ is the convection term coefficient, u is the hybrid viscosity, Γ h is the diffusion coefficient, represents the Hamiltonian operator.

[0035] Preferably, the porous skeleton energy equation is:

[0036]

[0037] is the non-steady-state term of the porous skeleton energy equation, is the diffusion term of the porous skeleton energy equation, q sf is the source term, which represents the fluid-solid heat transfer;

[0038] ε is the porosity, t is the time, ρ s is the solid phase density, c p,s is the solid specific heat, k f,eff is the solid phase effective thermal conductivity, T s is the solid phase temperature, represents the Hamiltonian operator.

[0039] Preferably, in Step 1, the initialization function of the FLUENT platform is used to initialize the mixing velocity u in the continuity equation of the cooling fluid, the mixing pressure p in the momentum equation of the cooling fluid, the mixing enthalpy h in the energy equation of the cooling fluid, and the solid phase temperature T in the energy equation of the porous skeleton. s Perform initialization assignment;

[0040] In Step 6, boundary conditions are assigned to the cooling fluid continuity equation and the cooling fluid momentum equation through the FLUENT platform, and the cooling fluid continuity equation and the cooling fluid momentum equation are solved through its own flow model.

[0041] Compared with the prior art, the present invention has the following advantages:

[0042] (1) This paper proposes a new simulation model for phase change sweat cooling heat and mass transfer that takes into account non-equilibrium effects, which improves the simulation prediction accuracy under extreme environments and can support the refined design of phase change sweat cooling thermal protection for the next generation of high-speed aircraft;

[0043] (2) This paper develops a numerical simulation method based on the secondary development of the FLUENT platform for a new model of non-equilibrium heat and mass transfer, laying the foundation for rapid iteration in the scheme stage. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 It is a schematic diagram of the technical solution of the present invention;

[0045] Figure 2 The figure is a schematic diagram of the solution process of the present invention. DETAILED DESCRIPTION

[0046] like Figure 1 As shown, the typical phase change sweating cooling physical model is that the cooling fluid enters the porous skeleton region 2 from the inlet 1, and forms a cooling fluid domain 3 in the porous skeleton. After the fluid in the cooling fluid domain 3 completes phase change heat transfer with the porous skeleton, it flows out from the outlet 4;

[0047] The technical solution of the present invention is: the control equations describing the phase change sweating cooling non-equilibrium flow heat transfer process are shown in equations (1) to (4), and the meanings and determination methods of the various parameters in the control equations are shown in Table 1; the present invention solves the continuity equation and momentum equation of the cooling fluid through the existing functions of the FLUENT software; on this basis, a UDF (Userdefine functions, user-defined functions) program is developed to solve the energy equation of the cooling fluid and the energy equation of the porous skeleton; the UDF program includes a DEFINE_INIT function, a DEFINE_PROPERTY function, a DEFINE_ADJUST function, a DEFINE_UDS_USNTEADY function, a DEFINE_UDS_FLUX function, a DEFINE_UDS_DIFFUSIVITY function, a DEFINE_SOURCE function, and a DEFINE_PROFILE function; the DEFINE_INIT function is used to initialize the saturation s and the fluid temperature T f ; The DEFINE_PROPERTY function is used to define and calculate the fluid mixture density, hybrid viscosity, diffusion term coefficient and convection term coefficient; the DEFINE_ADJUST function is used to adjust the fluid saturation s and temperature T after each iteration step. f , mixed density and mixed power viscosity; the DEFINE_UDS_USNTEADY function is used to define and solve the non-steady-state terms of the fluid and porous skeleton energy equations; the DEFINE_UDS_FLUX function is used to define and solve the convection term of the fluid energy equation; the DEFINE_UDS_DIFFUSIVITY function is used to define and solve the diffusion term of the fluid energy equation; the DEFINE_SOURCE function is used to define and solve the source term (heat exchange between the fluid and the porous skeleton); the DEFINE_PROFILE function is used to define the boundary conditions of the cooling fluid energy equation and the porous skeleton energy equation.

[0048] When boundary conditions are assigned to the cooling fluid continuity equation, cooling fluid momentum equation, cooling fluid energy equation, and porous skeleton energy equation: the velocity and temperature values ​​are given at the cooling fluid inlet, the pressure value is given at the cooling fluid outlet, the cold end of the porous skeleton is adiabatic, and the heat flux density value is given at the hot end of the porous skeleton.

[0049] Equations (1) to (4) and the parameter calculations in Table 1 constitute a new simulation model for phase change sweating cooling heat and mass transfer that can take into account non-equilibrium effects. Specifically:

[0050] Cooling fluid continuity equation:

[0051]

[0052] Cooling fluid momentum equation:

[0053]

[0054] Cooling fluid energy equation:

[0055]

[0056] in, is the unsteady-state term of the cooling fluid energy equation, is the cooling fluid energy convection term, is the cooling fluid energy diffusion term, q sf is the source term;

[0057] Porous skeleton energy equation:

[0058]

[0059] in, is the non-steady-state term of the porous skeleton energy equation, is the energy diffusion term of the porous skeleton, q sf is the source term;

[0060] Table 1 Meaning and determination method of parameters in control equation

[0061]

[0062]

[0063]

[0064]

[0065] The workflow of the present invention is as follows Figure 2 As shown, the details are as follows:

[0066] ① Through the initialization function of FLUENT, the mixing velocity u, mixing pressure p, mixing enthalpy h, solid phase temperature T in equations (1) to (4) are s Perform initialization assignment;

[0067] ②Use the DEFINE_INIT function of the UDF program to determine the fluid saturation s and fluid temperature T f Perform initialization assignment;

[0068] ③According to the fluid temperature T f , solve the physical parameters of liquid and gas phases by looking up tables or fitting formulas. The physical parameters of liquid and gas phases include: liquid phase density ρ l , gas phase density ρ v, liquid phase dynamic viscosity μ l , gas phase dynamic viscosity μ v , liquid phase saturation enthalpy h l , gas phase saturation enthalpy h v , latent heat of vaporization h fg , liquid specific heat c l Relative heat of gas c v ;

[0069] ④ According to the liquid and gas phase physical properties and fluid saturation s obtained in the previous step, the mixed density ρ, hybrid viscosity μ, convection term coefficient γ and diffusion term coefficient Γ are calculated through the DEFINE_PROPERTY function in the UDF program. h :

[0070] ρ=sρ l +(1-s)ρ v

[0071]

[0072]

[0073]

[0074] Where s is the liquid saturation, ρ l is the liquid density, ρ v is the gas phase density, k rl is the liquid phase relative permeability, μ l is the liquid phase dynamic viscosity, k rv is the gas phase relative permeability, μ v is the gas phase dynamic viscosity, λ is the liquid phase relative fluidity, h l is the liquid phase saturation enthalpy, h v is the gas phase saturation enthalpy, k f,eff is the effective thermal conductivity of the fluid, T f is the fluid temperature, h is the mixing enthalpy, h fg is the latent heat of vaporization, and D(s) is the capillary diffusion coefficient.

[0075] ⑤ Assign boundary conditions to equation (1) (cooling fluid continuity equation) and equation (2) (cooling fluid momentum equation) through the FLUENT platform, and assign boundary conditions to equation (3) and equation (4) through the DEFINE_PROFILE function;

[0076] ⑥ Use the Boundary Conditions panel in FLUENT to assign boundary conditions to equations (1) and (2), and solve equations (1) and (2) using the built-in flow model;

[0077] Divide the cooling fluid energy equation into the non-steady-state term of the cooling fluid energy equation Convection term of the cooling fluid energy equation Diffusion term in the cooling fluid energy equation and the source term q sf ; Divide the porous skeleton energy equation into the non-steady-state term of the porous skeleton energy equation Diffusion term in the energy equation of porous skeletons and the source term q sf ;

[0078] Solve equations (3) and (4) through the DEFINE_UDS_USNTEADY function, DEFINE_UDS_FLUX function, DEFINE_UDS_DIFFUSIVITY function, and DEFINE_SOURCE function to obtain the new u new , P new 、h new 、T s-new ;

[0079] ⑦ Determine the new mixing speed u new 、New mixed pressure P new 、New mixing enthalpy h new , the new solid phase temperature T s-new Compared with the initial u, P, h, T in this iteration s If the difference between the two meets the requirements, the calculation can be ended, otherwise u new , P new 、h new 、T s-new Assign values ​​to u, P, h, T respectively s , according to h and T f , s, and updates T through the DEFINE_ADJUST function f , s, and then repeat ③~⑦.

[0080] DEFINE_ADJUST function updates T f , s are based on the calculation method in No. 18 and No. 34 in Table 1, specifically:

[0081]

[0082]

[0083] Among them, ρ l is the liquid density, ρ v is the gas phase density, ρ is the mixed density, h is the mixing enthalpy, and h l is the liquid enthalpy, h v is the gas enthalpy, c p,l is the liquid specific heat, c p,vis the relative heat of gas, T sat is the saturation temperature.

[0084] The contents not described in detail in the specification of the present invention belong to the prior art known to the professional and technical personnel in this field.

Claims

1. A numerical simulation method for non-equilibrium flow heat transfer in phase change sweating cooling, characterized in that include: Step 1: The mixing velocity u in the cooling fluid continuity equation, the mixing pressure p in the cooling fluid momentum equation, the mixing enthalpy h in the cooling fluid energy equation, and the solid phase temperature T in the porous skeleton energy equation are s Perform initialization assignment; Step 2: Fluid saturation s and fluid temperature T f Perform initialization assignment; Step 3: According to the fluid temperature T f Obtain the physical properties of liquid and gas phases; Step 4. Calculate the mixed density ρ, mixed flow viscosity μ, convection term coefficient γ and diffusion term coefficient Γ based on the obtained liquid and gas phase physical parameters and fluid saturation s. h ; Step 5. Assign boundary conditions to the cooling fluid continuity equation, cooling fluid momentum equation, cooling fluid energy equation, and porous skeleton energy equation; Step 6. Solve the cooling fluid continuity equation and the cooling fluid momentum equation; The cooling fluid energy equation and the porous skeleton energy equation are divided into unsteady-state terms, convection terms, diffusion terms and source terms. The above four terms are expressed by writing UDF functions and loaded into the UDS equation in FLUENT. The cooling fluid energy equation and the porous skeleton energy are solved respectively to obtain the new mixing velocity u new , the new mixing pressure p new 、New mixing enthalpy h new , the new solid phase temperature T s-new ; Step 7, determine the new u new 、p new 、h new 、T s-new Compared with the initial u, p, h, T in this iteration s If the difference between the two meets the requirements, the calculation can be ended, otherwise u new 、p new 、h new 、T s-new Assign values ​​to u, p, h, and T respectively s , according to h and T f , the relationship between s is updated T f , s, return to Step3 and iterate again.

2. The method for numerical simulation of non-equilibrium flow heat transfer in phase change sweating cooling according to claim 1, characterized in that: Use UDF program to write functions for solving, specifically: The initialization function is used to initialize the saturation s and fluid temperature T f ; The characteristic definition function is used to define and calculate the fluid mixture density ρ, the hybrid viscosity μ, the convection term coefficient γ and the diffusion term coefficient Γ h ; The parameter adjustment function is used to adjust the fluid saturation s and temperature T after each iteration. f , mixture density ρ, hybrid viscosity μ, convection term coefficient γ and diffusion term coefficient Γ h ; The unsteady-state term definition function defines and solves the unsteady-state terms of the fluid and porous framework energy equations; The convection term definition function is used to define and solve the convection term of the fluid energy equation; The diffusion term definition function is used to define and solve the diffusion term of the fluid energy equation; The heat exchange definition function is used to define and solve the heat exchange between the fluid and the porous framework; The boundary definition function is used to define the boundary conditions of the porous skeleton energy equation and the cooling fluid energy equation.

3. The method for numerical simulation of non-equilibrium flow heat transfer in phase change sweating cooling according to claim 1, characterized in that: According to the fluid temperature T f The physical properties of the liquid and gas phases are solved by looking up tables or fitting formulas. The physical properties of the liquid and gas phases include: liquid phase density ρ l , gas phase density ρ v , liquid phase dynamic viscosity μ l , gas phase dynamic viscosity μ v , liquid phase saturation enthalpy h l , gas phase saturation enthalpy h v , latent heat of vaporization h fg , liquid specific heat c l Relative heat of gas c v .

4. The method for numerical simulation of non-equilibrium flow heat transfer in phase change sweating cooling according to claim 1, characterized in that: Mixed density ρ, mixed viscosity μ, convection term coefficient γ and diffusion term coefficient Γ h The calculation methods are: p = s p l +(1-s)ρ v Where s is the liquid saturation, ρ l is the liquid density, ρ v is the gas phase density, k rl is the liquid phase relative permeability, μ l is the liquid phase dynamic viscosity, k rv is the gas phase relative permeability, μ v is the gas phase dynamic viscosity, λ is the liquid phase relative fluidity, h l is the liquid phase saturation enthalpy, h v is the gas phase saturation enthalpy, k f,eff is the effective thermal conductivity of the fluid, T f is the fluid temperature, h is the mixing enthalpy, h fg is the latent heat of vaporization, and D(s) is the capillary diffusion coefficient.

5. The method for numerical simulation of non-equilibrium flow heat transfer in phase change sweating cooling according to claim 1, characterized in that: Boundary conditions are assigned to the cooling fluid continuity equation, cooling fluid momentum equation, cooling fluid energy equation, and porous skeleton energy equation, specifically: The velocity and temperature values ​​are given at the cooling fluid inlet, the pressure value is given at the cooling fluid outlet, the cold end of the porous skeleton is adiabatic, and the heat flux density value is given at the hot end of the porous skeleton.

6. The method for numerical simulation of non-equilibrium flow heat transfer in phase change sweating cooling according to claim 1, characterized in that: The cooling fluid energy equation is: in, is the unsteady-state term of the cooling fluid energy equation, is the convection term of the cooling fluid energy equation, is the diffusion term of the cooling fluid energy equation, q sf is the source term, which represents the fluid-solid heat transfer; ε is the porosity, t is the time, ρ is the mixing density, h is the mixing enthalpy, γ is the convection term coefficient, u is the hybrid viscosity, Γ h is the diffusion coefficient, represents the Hamiltonian operator.

7. The method for numerical simulation of non-equilibrium flow heat transfer in phase change sweating cooling according to claim 1, characterized in that: The energy equation of the porous skeleton is: is the non-steady-state term of the porous skeleton energy equation, is the diffusion term of the porous skeleton energy equation, q sf is the source term, which represents the fluid-solid heat transfer; ε is the porosity, t is the time, ρ s is the solid phase density, c p,s is the solid specific heat, k f,eff is the solid phase effective thermal conductivity, T s is the solid phase temperature, represents the Hamiltonian operator.

8. The method for numerical simulation of non-equilibrium flow heat transfer in phase change sweating cooling according to claim 1, characterized in that: In Step 1, the initialization function of the FLUENT platform is used to initialize the mixing velocity u in the continuity equation of the cooling fluid, the mixing pressure p in the momentum equation of the cooling fluid, the mixing enthalpy h in the energy equation of the cooling fluid, and the solid phase temperature T in the energy equation of the porous skeleton. s Perform initialization assignment; In Step 6, boundary conditions are assigned to the cooling fluid continuity equation and the cooling fluid momentum equation through the FLUENT platform, and the cooling fluid continuity equation and the cooling fluid momentum equation are solved through its own flow model.

Images