# The Finite Element Numerical Study of Transonic Flow of Moist Air with Nonequilibrium Condensation

Authors: Weijun Tian; Liquan Mei
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The finite element numerical simulation of the transonic flow with condensation is studied in this paper. The transformed weak form is numerically discretized by using the classical Galerkin finite element spatial discretization scheme and the backward Euler difference time discretization scheme. Because the coupled Euler equations are convection equations, the results obtained by the classical Galerkin method often have nonphysical numerical oscillations, so the stabilization method is introduced. Another important feature of compressible Euler equation is the existence of shock wave. Because the thickness of shock wave is small and difficult to capture, the shock wave can be better captured by adding artificial viscosity term. Here, the streamline upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability are used to solve the problem of transonic flow with condensation. In the simulation of steady problems, in order to improve the calculation accuracy of the existing area of the shock wave, the adaptive mesh refinement technology is added, which refines the mesh in the region of the shock wave to improve the resolution of the shock wave. Numerical experiments verify the feasibility and stabilization of the numerical method. Finally, the simulation of transonic flow around NACA-0012 and parallel-jet nozzle A1 has obtained the better numerical results.


---

## Body

## 1. Introduction

As for human beings, water is one of the most familiar and needed substances. As we all know, the ocean accounts for three quarters of the global area; in addition, water is a significant component of oceans and lakes. Due to the long-term irradiation of the sun, the water on land evaporates. Meanwhile, water becomes water vapor around us. The wet air mixed by water vapor and air is a component of the troposphere. In the troposphere, the temperature decreases with the increase of the height, and the temperature decreases 6.5°C with the height of 1 km increase. Due to the decrease of temperature in the rise of air, water vapor condenses on the tiny dust particles in the air. The heat released in this process increases the temperature of the surrounding air. Depending on the surrounding environment, these small droplets can continue to increase and form rain clouds to a certain extent. Due to the influence of the Earth's gravity, these small droplets return to the Earth in the form of rain. Nowadays, many famous scholars devote themselves to the research on some possible novel practical applications; for example, Costa et al. [1] proposed the INNOVARE project: innovative plants for distributed poly-generation by residual biomass; Iorio et al. [2] analyzed low enthalpy geothermal systems in structural controlled areas: a sustainability analysis of geothermal resource for heating plant, the Mondragone case in southern Apennines.Cloud formation is a slow condensation process. It occurs when the gas phase reaches saturation or dew point temperature. In fact, it is only suitable for slow condensation phenomena in nature, such as the formation of dew and fog. It is not completely true for all condensation phenomena, but many rapid condensation phenomena occur in nonequilibrium states, for example, the formation of clouds over a tornado. The wind speed around the tornado center is very large, resulting in the air pressure in the center being several times lower than the standard atmospheric pressure, resulting in a sharp drop in the temperature in the center. Therefore, when passing through the water surface, the water sucked up will soon become clouds. We cannot observe a large tornado. However, this phenomenon in a small range also occurs in the vortex generated by the landing of transport aircraft. In the early high-speed wind tunnel, the wet air will condense during the test process. In order to avoid the condensation of wet air in the process of high-speed expansion, experts modified the design of the wind tunnel.In the process of high-speed expansion of wet air, because of the accelerated condensation of condensable water vapor, the original gas changes from unsaturated state to supersaturated metastable state, which occurs in the channel of steam turbine blade and the wing of supersonic aircraft. Condensation will not occur outside the saturated state, and a certain degree of supercooling will be maintained before condensation. This phenomenon can be observed in supersonic nozzle flow and sparse flow in shock tube, but the change of thermodynamic state only occurs in a few milliseconds. The change range of cooling rate spans 10−3‒100 k/μs; due to condensation, the cooling rate is as high as −102 k/μs. In this case, because of the latent heat released by condensation, we can observe a part of small droplet temperature. The pressure and density increase with the occurrence of condensation, and the velocity of the fluid decreases, so that the two-phase nonequilibrium mixture reaches a new equilibrium state. When the heat released into the fluid exceeds the critical amount, the condensation will induce shock waves, such as stable standard shock wave and shock wave of unsteady periodic flow of nozzle flow. The driving part of the shock tube is unstable, and the shock wave propagates in the opposite direction of the fluid.Wet air condensation has been concerned for a long time. Kotake and Class [3] summarized the application fields of wet air condensation, mainly including atmospheric physics, astrophysics, aviation science, nuclear reaction technology, and material science. The condensation of wet air is the formation of clouds, fog, and dew in meteorology and climatology.


### 1.1. The Main Contributions and Thesis Overview

At the end of the 19th century, Boucher [4], Wegener and Mirable [5], Kotake and Class [3], Head [6], and Wegener and Mack [7] expounded the research history of condensation. In 1936, Prandtl [8] proposed the condensation-induced shock wave of moist air supersonic nozzle flow at the Volta Conference held in Rome. So far, the research on the condensation of wet air in high-speed expanding fluid is still in progress [9–11].At first, many authors used non-CFD (computational fluid dynamics) models to study condensation phenomena, but with the deepening of research, CFD models are mostly used to study condensation phenomena in high-speed wind tunnels [12–17]. In 1985, Robinson et al. [18] studied supersonic wet air nozzle flow. In 1990, Schnerr and Dohrmann [19] used a model combining Euler equation and condensation theory to study the transonic wet air flow around the wing. Yamamoto [20] continued Schnerr and Dohrmann’s work and used viscous N–S equation to study the condensation of vortices on delta wings. In 2005, in more in-depth research, Goodhart and Schnerr [21] used the viscous RANS model to study the flow around the wing. At the same time, Yamamoto [22] also went deep into the airfoil and steam turbine.The numerical simulation results of the transonic flow of wet air can make us better understand the flow process around the aircraft wing in high-altitude flight, which is of great significance to the design of airfoil. Transonic flow can be divided into external flow (such as flow around airfoil and fuselage) and internal flow (such as flow at nozzle and cascade), two aspects. The study of transonic flow has wide application value, because most of various types of aircraft are in the transonic range. The cascade flow in aeroengine, the throat flow in rocket engine nozzle, and the flow in various gas flow meters and shut-off valves are often transonic flow. In industry, condensation is an unwanted side effect, and the heat released when wet air condenses will cause the decline of machine performance. Numerical simulation enables us to understand more complex multiphase flow physical phenomena, so as to disperse the huge pressure and shock wave produced by air on the object.In this paper, CFD is mainly used as a tool to predict the condensation of inviscid flow of wet air around the wing. The method is to describe the flow of inviscid flow based on Euler equation and then introduce the condensation of fluid into the equation based on the classical nucleation and droplet expansion theory combined with hill [26]. In this paper, it is assumed that the fluid is stable.The second section mainly introduces the physical background of inviscid flow condensation, as well as the mathematical model and relationship of fluid motion, and describes the properties of many physical quantities, such as surface tension, density, and saturated vapor pressure. Therefore, Euler equation is often used as the governing equation to describe the flow field. For the condensation phenomenon of transonic flow in wet air, the energy equation of Euler equation combines the classical nucleation theory and the Hertz–Knudsen growth model of small droplets, and the hill moment method is used to describe the relationship between various variables of liquid phase. The ordinary differential equation of the change rate of liquid mass fraction is transformed into four equivalent partial differential equations, which are coupled with the previous Euler equation to form the governing equation of moist air transonic flow.In Section3, aiming at the mathematical problem of solving the partial differential equation transformed from this engineering problem, at present, scientists mostly use the finite difference and finite volume method to solve this equation. In this paper, the unstructured grid is used to solve this multi-field coupling problem, and the classical Galerkin finite element method is used to discretize the space to obtain the semi-discrete equation. Then, the backward Euler scheme is used to discretize the time. In order to better solve the steady problem, the local mesh encryption technology is introduced. In fact, the finite element method will produce nonphysical numerical oscillation and cause the instability of the numerical solution. Moreover, in the high-speed expansion flow, the condensation will induce shock waves, but the classical finite element method is not ideal to capture shock waves. Therefore, when solving the problem, the stabilization method is introduced, and the original finite element scheme is modified by using the streamline upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability to find the decoupled equations, so as to achieve better numerical results.In Section4, the method proposed in Section 3 is used for numerical simulation of practical problems, mainly dealing with the transonic inviscid flow field around NACA-0012 airfoil, so as to verify whether the above scheme can obtain satisfactory results in the complex Euler flow model. At the same time, the dry air flow and wet air flow of parallel tail nozzle A1 are also numerically simulated. The feasibility of the numerical method is further verified, and the condensation of wet air in the nozzle is simulated.The fifth section is the conclusion and development trend of the research, summarizes the main contributions of this paper, and the advantages and disadvantages of the method, and looks forward to the development.


## 1.1. The Main Contributions and Thesis Overview

At the end of the 19th century, Boucher [4], Wegener and Mirable [5], Kotake and Class [3], Head [6], and Wegener and Mack [7] expounded the research history of condensation. In 1936, Prandtl [8] proposed the condensation-induced shock wave of moist air supersonic nozzle flow at the Volta Conference held in Rome. So far, the research on the condensation of wet air in high-speed expanding fluid is still in progress [9–11].At first, many authors used non-CFD (computational fluid dynamics) models to study condensation phenomena, but with the deepening of research, CFD models are mostly used to study condensation phenomena in high-speed wind tunnels [12–17]. In 1985, Robinson et al. [18] studied supersonic wet air nozzle flow. In 1990, Schnerr and Dohrmann [19] used a model combining Euler equation and condensation theory to study the transonic wet air flow around the wing. Yamamoto [20] continued Schnerr and Dohrmann’s work and used viscous N–S equation to study the condensation of vortices on delta wings. In 2005, in more in-depth research, Goodhart and Schnerr [21] used the viscous RANS model to study the flow around the wing. At the same time, Yamamoto [22] also went deep into the airfoil and steam turbine.The numerical simulation results of the transonic flow of wet air can make us better understand the flow process around the aircraft wing in high-altitude flight, which is of great significance to the design of airfoil. Transonic flow can be divided into external flow (such as flow around airfoil and fuselage) and internal flow (such as flow at nozzle and cascade), two aspects. The study of transonic flow has wide application value, because most of various types of aircraft are in the transonic range. The cascade flow in aeroengine, the throat flow in rocket engine nozzle, and the flow in various gas flow meters and shut-off valves are often transonic flow. In industry, condensation is an unwanted side effect, and the heat released when wet air condenses will cause the decline of machine performance. Numerical simulation enables us to understand more complex multiphase flow physical phenomena, so as to disperse the huge pressure and shock wave produced by air on the object.In this paper, CFD is mainly used as a tool to predict the condensation of inviscid flow of wet air around the wing. The method is to describe the flow of inviscid flow based on Euler equation and then introduce the condensation of fluid into the equation based on the classical nucleation and droplet expansion theory combined with hill [26]. In this paper, it is assumed that the fluid is stable.The second section mainly introduces the physical background of inviscid flow condensation, as well as the mathematical model and relationship of fluid motion, and describes the properties of many physical quantities, such as surface tension, density, and saturated vapor pressure. Therefore, Euler equation is often used as the governing equation to describe the flow field. For the condensation phenomenon of transonic flow in wet air, the energy equation of Euler equation combines the classical nucleation theory and the Hertz–Knudsen growth model of small droplets, and the hill moment method is used to describe the relationship between various variables of liquid phase. The ordinary differential equation of the change rate of liquid mass fraction is transformed into four equivalent partial differential equations, which are coupled with the previous Euler equation to form the governing equation of moist air transonic flow.In Section3, aiming at the mathematical problem of solving the partial differential equation transformed from this engineering problem, at present, scientists mostly use the finite difference and finite volume method to solve this equation. In this paper, the unstructured grid is used to solve this multi-field coupling problem, and the classical Galerkin finite element method is used to discretize the space to obtain the semi-discrete equation. Then, the backward Euler scheme is used to discretize the time. In order to better solve the steady problem, the local mesh encryption technology is introduced. In fact, the finite element method will produce nonphysical numerical oscillation and cause the instability of the numerical solution. Moreover, in the high-speed expansion flow, the condensation will induce shock waves, but the classical finite element method is not ideal to capture shock waves. Therefore, when solving the problem, the stabilization method is introduced, and the original finite element scheme is modified by using the streamline upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability to find the decoupled equations, so as to achieve better numerical results.In Section4, the method proposed in Section 3 is used for numerical simulation of practical problems, mainly dealing with the transonic inviscid flow field around NACA-0012 airfoil, so as to verify whether the above scheme can obtain satisfactory results in the complex Euler flow model. At the same time, the dry air flow and wet air flow of parallel tail nozzle A1 are also numerically simulated. The feasibility of the numerical method is further verified, and the condensation of wet air in the nozzle is simulated.The fifth section is the conclusion and development trend of the research, summarizes the main contributions of this paper, and the advantages and disadvantages of the method, and looks forward to the development.


## 2. Condensation Flow Model

Solid, liquid, and gas are the three most common phases. The so-called phase refers to the part with completely uniform physical properties and chemical components in the system. Many substances can exist in more than one material state. From the perspective of molecular dynamics, the emergence of different material states is due to the different orders of molecules in various states, and substances can exist in single phase. It can also coexist in equilibrium in the form of two or three phases. The thermodynamic phase diagram of pure substances is shown in Figure1. In the diagram, it, respectively, represents the coexistence areas of solid-liquid, liquid-gas, and solid-gas phases, which can become the dissolution surface, vaporization surface, and sublimation surface. Either side represents a single-phase area of solid, liquid, and gas, and the intersection of the three surfaces represents the equilibrium coexistence of solid, liquid, and gas. It is called a three-phase line (TPL). For a certain substance, the temperature and pressure of its three-phase point are certain. For example, the three-phase point temperature of water is 273.15K and the pressure is 0.6112 kPa. Another important point is the critical point, indicating the limit state of phase transition. At this time, the interface between liquid and gas disappears, and the vaporization process becomes a point at the critical pressure, indicating that the vaporization is completed in an instant. There is no longer a wet steam state in which vapor-liquid two phases coexist, and the difference between liquid and gas disappears above the critical temperature. See Table 1 for the critical point parameters of some substances.Figure 1
Thermodynamic phase diagram of pure substance [23].Table 1
Critical point parameters of some substances [24].
MaterialCritical point temperatureTc/KCritical point pressurep-/(MPa)Helium (He)5.10.234Hydrogen (H2)33.31.30Air132.53.77Water (H2O)647.322.09Due to the high-speed expansion, the equilibrium state of the fluid is destroyed and the fluid reaches the supersaturation state. If there is no condensation nodule in the wet air or the condensation nucleus is too small, even if the steam pressure exceeds the saturated steam pressure at this temperature, or even several times higher, the droplets cannot form or grow up, which is called the supersaturation phenomenon. The degree of nonequilibrium state can be expressed by saturationΦ. It is defined as the ratio of the actual water vapor pressure Pv to the water vapor pressure Pv,eqT at equilibrium.(1)Φ=PvPv,eqT.IfΦ>1, it is called supersaturation. Supersaturation is a metastable state. With slight external interference, such as the addition of dust, impurities, or charged particles, it can be used as a sufficiently large condensate, so that there can be growing droplets in the supersaturated vapor, which changes from metastable state to stable gas-liquid equilibrium state.Therefore, the first stage of condensation is called nucleation process. Different kinds of substances affect the initial stage of droplet formation, such as small particles and aerosols on the molecular surface (suspended particles in gas, such as smoke and fog). However, due to rapid expansion, the nucleation of the same kind of substances far exceeds the number of particles formed by the nucleation of different substances. The second stage of condensation is called the growth process of small droplets. In this stage, the thermodynamic equilibrium is maintained, and the vapor molecules on the equilibrium condensation core continue to combine, and the latent heat generated by the combination will increase the temperature and pressure of the fluid and slow down the growth rate of liquid droplets. Only those small droplets that reach the critical radius can continue to grow. Small droplets that do not reach the critical radius cannot condense because they are easy to evaporate due to instability.The same kind of nucleation requires many gas molecules to form metastable molecular clusters in supersaturated state, which plays an important role in condensation. Nucleation is the first step in the formation of liquid droplets. The supersaturation rate in nucleation theory is given by equation (1). The change value of free energy when a molecular cluster is composed of molecules is given by Wisman [25],(2)ΔGn=σa0n2/3−nkBTlnΦ,kB is a universal constant of micro-particle behavior, called the Boltz constant. kB=1.38066×10−23J/K, and the appearance of kB characterizes the system related to heat. σ is the surface tension of the flat liquid surface, which represents the surface area of a molecule, defined as(3)a0=36π1/3vl2/3,vl represents the volume of liquid molecules unit:m3/cluster.The first part at the right end of equation (2) represents the volume strain energy caused by the shape mismatch between the nuclear embryo (liquid) and the matrix (gas), and the second part represents the reduced free energy for the formation of nuclear embryo. We can see that when Φ<1, ΔGn is a monotonically increasing function with respect to n. Therefore, the phase transition is avoided by reducing the value of n.WhenΦ>1, there is always a value of n, which makes the wet steam prone to phase transition. The critical value of n∗ is determined by ∂ΔGn/∂nn∗=0:(4)n∗=2θ3lnΦ3,θ is the dimensionless surface tension, defined as(5)θ=σa0kBT.The critical free energy is obtained by substitutingn∗ into equation (2).(6)ΔGn∗=427kBTθ3ln2Φ=43πr∗2σ.Then, usen∗υl=4/3πr∗3 to get the Kelvin relation:(7)r∗=2σρl0RvTlnΦ,In the formula, ρl0 is condensation density and Rv is gas constant of water vapor.The expression of critical radiusr∗ is given by Kelvin relation, and r∗ is also the minimum radius of the molecular cluster when the supersaturated state reaches equilibrium. In other words, only those molecular clusters that reach the critical radius can continue to grow, while those less than the critical radius decompose due to instability.The Kelvin relation has the following relations:vl=M/NAρl0, Mkg/mol is the molar mass of the gas, and NA=6.02×1026/mol is the Avogadro constant.Classical nucleation theory is widely used in the study of transonic condensation flow. It is widely used in the study of transonic condensation flow.Nucleation rate is defined as the number of cores formed per unit volume in the parent phase matrix without transformation or interaction between components, recorded asJNumber of cores/m3⋅s, and the definition formula is(8)J=Ke−ΔGn∗/kBT.ΔGn∗ is given by equation (6), and K is a constant.The condensation rate is given by the classical nucleation theory as follows:(9)JCNT=ρv2ρl02σπm3e−4/27θ3/ln2Φ.In the formula, ρv is density of water vapor and m is the mass of a water vapor molecule.The growth rate of small droplets is expressed by the growth rate of small droplets. The growth rate of small droplets in this paper mainly adopts the famous Hertz–Knudsen growth model [26]. Hertz–Knudsen growth theorem is established based on the balance between the condensation of gas molecules on the surface of small droplets and the evaporation of gas molecules in small droplets,(10)r˙=αρl0pv2πRvT−pv,r2πRvTd.α is a free parameter, generally let α=0.2, Td is the temperature of small droplets, pv is the pressure of water vapor, and pv,r is the supersaturated vapor pressure with a small droplet radius of r, defined as(11)pv,r=pv,eqe2σ/rρl0RvTd.We usually assume that the temperature of small droplets is the same as that of the surrounding air, that is,T=Td.The Euler equation is extended by Luijten [26], which well describes the flow of high-speed expansion condensation flow. We compare the numerical results with the experimental data and come to the conclusion that the viscous effect of the fluid is negligible, and the viscous term may not be included in the equation.When the composition of the gas is constant, the Euler equation can be obtained from the conservation of mass, momentum, and energy.(12)∂U∂t+∂F1∂x+∂F2∂y=0.In the formula, U is vector consisting of conserved variables; Fi is vector consisting of convective flux i=1,2; and U, Fi is shown below:(13)U=ρρuρvρet,(14)F1=ρuρu2+pρuvρuht,F2=ρvρuvρv2+pρvht.ρ is density. p is pressure. u,v are the velocity components of the direction of x,y. The total energy et is the sum of internal energy e and kinetic energy 1/2u2+v2,(15)et=e+12u2+v2.And the total enthalpyht is defined as below:(16)ht=et+pρ.The Euler equation needs to be solved in combination with the equation of state. The equation of state in thermodynamics will be given in the following discussion.Consider a very small liquid containing a certain number of small droplets with a volume ofV mixed gas, and the mixed gas consists of inert gas and steam, for example, a mixture of dry air and water vapor. Inert gas refers to a gas that does not condense, and its mass is recorded as Ma. The mass of steam that can condense is recorded as Mv. The mass of the condensed liquid is recorded as Ml; thus, volume V of total mass of mixed gas is M:(17)M=Ma+Mv+Ml.When the total massM is constant, because the inert gas does not condense, the sum of the mass of the condensable steam and the mass of the condensed liquid is a constant, recorded as Mv0,(18)Mv0=Mv+Ml.The ratio of the mass of the condensed liquid to the total mass is called the mass fraction of the liquid and is recorded asg,(19)g=MlM,0≤g≤gmax.The mass fraction of liquid reaches the maximum when all steam condenses, which is called the maximum mass fraction of liquid, i.e.,(20)gmax=Mv0M.At a certain volume, the density of the mixture is defined as(21)ρ≡MV.Similarly, the density of inert gas and steam is(22)ρa=MaV=1−gmaxρ,(23)ρv=MvV=gmax−gρ.The density of the condensed liquidρl (not condensation density ρl0) similar can be expressed as(24)ρl=MlV=gρ.From equations (23) and (24), the maximum mass fraction of the liquid can be obtained, which can also be expressed as(25)gmax=ρv0ρ.Andρv0 for the sum of ρv and ρl,(26)ρv0=ρv+ρl.The first equation of equation (12) is the mass conservation equation, which is defined u⟶=u,v as the velocity of inert gas and steam mixed gas, so the mass conservation equation of each component is(27)∂ρa∂t+∂ρau∂x+∂ρav∂y=0,(28)∂ρv0∂t+∂ρv0u∂x+∂ρv0v∂y=0.Note that the precondition for the establishment of equations (27) and (28) is that the velocity of gas and condensed liquid is the same. Therefore, we assume that there is no slip between gas phase and liquid phase. Bring equation (25) into equation (28),(29)∂∂tρgmax+∂∂xρgmaxu+∂∂yρgmaxv=0.Expand equation (29),(30)gmax∂ρ∂t+∂∂xρu+∂∂yρv+ρ∂gmax∂t+u∂gmax∂x+v∂gmax∂y=0.The part in braces of the first term represents the mass conservation equation of the mixed gas, and its value is 0. Becauseρ≠0, equation (30) is simplified as(31)∂gmax∂t+u∂gmax∂x+v∂gmax∂y=0.It can also be written as(32)DgmaxDt=0.Thus, it can be obtained thatgmax is a constant along a certain trajectory. From equation (32), we can get an important conclusion: if gmax is a constant on the inlet boundary, gmax is a constant on the whole flow field, so gmax can be regarded as a global constant.We can obtain the mass conservation equation of steam and liquid by substituting equation (26) into equation (28),(33)∂ρv∂t+∂ρvu∂x+∂ρvv∂y+∂ρl∂t+∂ρlu∂x+∂ρlv∂y=0.Then, combine formula (24) to obtain(34)∂ρg∂t+∂∂xρgu+∂∂yρgv=b,(35)∂ρv∂t+∂∂xρvu+∂∂yρvv=−b.In the formula,b is energy provided by phase transition.The mass conservation equation of steam can be obtained by substituting equation (23) into equation (35). Therefore, the mass conservation equations of mixed gas can be transformed into(36)∂ρ∂t+∂ρu∂x+∂ρv∂y=0,∂ρg∂t+∂ρgu∂x+∂ρgv∂y=b.The equations are composed of the mass conservation equation of mixed gas and the mass conservation equation of condensed liquid. However, the phase transition energy represented in the right termb is not known. This equation describes the condensation process and will be discussed in the following sections.In the general dynamic equation, the condensation and decomposition of small droplets follow the conservation condition. The phase space isx⟶,r=x,y,r, where ris the radius of small droplets. GDE can be referred to Seinfeld [27]:(37)∂f∂t+∂∂rr˙f+∂fu∂x+∂fv∂x=δr−r∗J,where f is distribution function; δ is Dirac delta function; J is condensation rate; and f is defined as the partial derivative of the number density. nr,x⟶,t of small droplets to the radius r is(38)fr,x⟶,t≡∂nr,x⟶,t∂r,(39)ornr,x⟶,t=∫0rfr,x⟶,tdr.In order to simplify the equation, the following independent variablesf are omitted.Both ends of equation (37) are multiplied by rk and integrated r from 0 to ∞, and(40)∂∂t∫0∞rkfdr+∫0∞rk∂∂rr˙fdr+∂∂x∫0∞rkfdru+∂∂y∫0∞rkfdrv=∫0∞rkδr−r∗drJ.From the properties of Dirac delta function,(41)∫0∞rkδr−r∗dr=r∗k.Introducingk−order distribution function,(42)μk≡∫0∞rkfdr,k=0,1,….According to equations (41) and (42), equation (40) can be simplified as(43)∂μk∂t+∂μku∂x+∂μkv∂y=r∗kJ−∫0∞rk∂∂rr˙fdr.The integral can be simplified, and(44)∫0∞rk∂∂rr˙fdr=rkr˙f|0∞−k∫0∞rk−1r˙fdr.The first term at the right end of the equation is 0 because no small droplet radius reaches infinity. In order to put forward the growth rate of small droplet radiusr˙ from the integral, Hill [26] introduced the average small droplet radius:(45)rH=μ2μ0.Therefore, the radius growth rater˙ of small droplets no longer depends on the radius rH of droplets, but depends on the “Hill” radius r, that is, r˙=r˙rH,x⟶,t; formula (44) becomes(46)∫0∞rk∂∂rr˙fdr=−kr˙∫0∞rk−1fdr=−kr˙μk−1.Therefore, formula (43) is transformed into(47)∂μk∂t+∂μku∂x+∂μkv∂y=r∗kJ+kr˙μk−1.This equation is aboutμk cycle. If k=0,…,3, partial differential equations can be obtained:(48)k=0:∂μ0∂t+∂μ0u∂x+∂μ0v∂y=J,k=1:∂μ1∂t+∂μ1u∂x+∂μ1v∂y=r∗J+r˙μ0,k=2:∂μ2∂t+∂μ2u∂x+∂μ2v∂y=r∗2J+2r˙μ1,k=3:∂μ3∂t+∂μ3u∂x+∂μ3v∂y=r∗3J+3r˙μ2.Introduce equation(49)μk≡ρQk.Qk is the so-called liquid moment, and equation (48) can be transformed into(50)∂ρQ0∂t+∂ρQ0u∂x+∂ρQ0v∂y=J,∂ρQ1∂t+∂ρQ1u∂x+∂ρQ1v∂y=r∗J+r˙ρQ0,∂ρQ2∂t+∂ρQ2u∂x+∂ρQ2v∂y=r∗2J+2r˙ρQ1,∂ρQ3∂t+∂ρQ3u∂x+∂ρQ3v∂y=r∗3J+3r˙ρQ2.The densityρl of condensed liquid can be expressed as the following distribution function:(51)ρl=∫0∞43πρl0r3fdr.We assume that when the material is constant, the condensation densityρl0 is constant, which is obtained from formulas (42) and (51),(52)ρl=43πρl0∫0∞r3fdr=43πρl0μ3.From equations (24) and (49),(53)g=43πρl0Q3.Therefore, equation (50) can be expressed as the equation of liquid mass fraction, and the following four equations can be obtained from equation (50), which is called “Hill moment method”(54)∂ρQ0∂t+∂ρQ0u∂x+∂ρQ0v∂y=J,∂ρQ1∂t+∂ρQ1u∂x+∂ρQ1v∂y=r∗J+r˙ρQ0,∂ρQ2∂t+∂ρQ2u∂x+∂ρQ2v∂y=r∗2J+2r˙ρQ1,∂ρg∂t+∂ρgu∂x+∂ρgv∂y=43πρl0r∗3J+3r˙ρQ2.The average radius of Hill is obtained by replacing equation (45) with liquid moment Qk,(55)rH=Q2Q0.In order to make Hill moment better describe the liquid mass fraction in PDE, we make the following assumptions:(1)
There is no slip between gas phase and liquid phase.(2)
Droplet growth rate ofr˙ is independent of droplet radius.(3)
The condensation densityρl0 is constant.The mass fractiong of condensed liquid in Hill's method of moments can also be expressed by the following expression (Hill [26]):(56)gt=∫0t43πρl0r3t,τJτρτdτ.Compared with formula (53) in the previous section, this expression is based on physical properties instead of being expressed by the third moment Q3. This expression is tenable on the assumption that when t=0. There are no small droplets, condensation, and decomposition at time. When the liquid density ρ10 does not depend on time, the derivative of g to time in formula (56) is(57)dgdt=43πρl0r∗3Jtρt+4πρl0∫0tJτρτr˙r2t,τdτ.Similar to formula (55), Hill radius rH is introduced, and since the growth rate r˙ of small droplets can be moved out of the integral, the following auxiliary variables can be introduced:(58)Q2=∫0tr2t,τJτρτdτ,Q1=∫0trt,τJτρτdτ,Q0=∫0tJτρτdτ,rH=Q2Q0.The equations can be gained,(59)∂ρQ0∂t+∂ρQ0u∂x+∂ρQ0v∂y=J,∂ρQ1∂t+∂ρQ1u∂x+∂ρQ1v∂y=r∗J+r˙ρQ0,∂ρQ2∂t+∂ρQ2u∂x+∂ρQ2v∂y=r∗2J+2r˙ρQ1,∂ρg∂t+∂ρgu∂x+∂ρgv∂y=43πρl0r∗3J+3r˙ρQ2.The liquid moment in the equations here is an integral about time, while the liquid moment in equation (50) is an integral about the droplet radius. Hagmeijer [28] has proved that the liquid moment obtained by these two integrals is equivalent. Therefore, equations (54) and (59) are the same.The governing equation of wet air condensation flow is composed of Euler equation and the equation obtained by Hill’s method of moments. It also needs the state equation of gas. The state equation is an expression composed of pressurep and temperature T, which is used to calculate flux and condensation terms.The static pressure of the gas mixture is the sum of the pressurepv of the steam and the pressure pa of the inert gas, i.e.,(60)p=pv+pa.Here, both inert gas and steam are regarded as ideal gases. Given that the two gas constants areRa and Rv, respectively, the thermodynamic equation of state of gas can be obtained:(61)pa=ρaRaTa,pv=ρvRvTv.Assuming thermodynamic equilibrium, the temperature of inert gas and steam is the same, i.e.,(62)Ta=Tv=T.The equation of state of the gas mixture is(63)p=ρRT.R is the gas constant of the gas mixture,(64)R=R0−gRv.R0 is the gas constant at which the mass of the condensed liquid is defined as(65)R0=1−gmaxRa+gmaxRv.The enthalpyhρ of the mixed gas is the sum of the enthalpies of the components, i.e.,(66)hρ=haρa+hvρv+hlρl.Combining equations (20), (22), and (23), the expression of specific enthalpy h (enthalpy per unit mass) can be obtained from equation (66):(67)h=1−gmaxha+gmax−ghv+ghl.The heat of vaporizationL is defined as the difference between the specific enthalpy hv of steam and the specific enthalpy hl of condensed liquid:(68)L=hv−hl,L>0.Substitute equation (68) into equation (67) to obtain(69)h=1−gmaxha+gmaxhv−gL.Similar to equation (67), the expression of specific internal energy e can be obtained:(70)e=1−gmaxea+gmax−gev+gel.According to the definition of formulah=e+p/ρ and formula (68), we can get(71)e=1−gmaxea+gmaxev−gL+gpv/ρv−plρl0.Assuming that the partial pressure of the condensed liquid is equal to the pressure of the mixed gas, i.e.,p=pl, and then using the equation of state of the ideal gas, equation (71) can be transformed into(72)e=1−gmaxea+gmaxev−gL+gRv−ρρl0RT.Complete the gas according to the calorimetry, and obtain(73)e=cvT.cv represents the specific heat under a certain volume, which is defined as(74)cv≡∂e∂Tv.Assuming that both inert gas and steam are ideal gases, formula (72) can be converted into(75)e=cv0T−gL+gRv−ρρl0RT.cv0 represents the specific heat of the liquid at a certain volume when the mass fraction of the liquid is, i.e.,(76)cv0=1−gmaxcva+gmaxcvv.If the specific internal energye, density ρ, and mass fraction of the liquid are known, the temperature of the mixed gas can be obtained from formula (75) and the pressure of the mixed gas can be obtained from formula (63).Becauseρ≪ρl0, the last term ρ/ρl0R at the right end of equation (75) can be ignored, i.e.,(77)e=cv0T+gRvT−L.According to the discussion of latent heat in section 2.3.4, when the mixed gas is a mixture of dry air and water vapor, the linear expression of latent heat is formula (99), and the temperature T can be obtained by substituting it into formula (77):(78)T=e+gL0cv0+gRv−L1.From the discussion in the previous section, the final form of the governing equations is obtained by combining the Euler equation with the “Hill moment” method(79)∂U∂t+∂F1∂x+∂F2∂y=B,U represents a vector composed of conserved variables:(80)UT=ρ,ρu,ρv,ρet,ρQ0,ρQ1,ρQ2,ρgT.Fii=1,2 represents the composition of convective flux:(81)F1=ρuρu2+pρuvρuhtρuQ0ρuQ1ρuQ2ρug,F2=ρvρuvρv2+pρvhtρvQ0ρvQ1ρvQ2ρvg.B represents a vector consisting of condensation source terms:(82)B=0000Jr∗J+r˙ρQ0r∗2J+2r˙ρQ143πρl0r∗3J+3r˙ρQ2.Sound velocity is a very important physical quantity, which determines the velocity of the fluid in the minimal disturbance propagation medium. Because the composition of the mixture is not invariable when condensation occurs, we use frozen sound velocity instead of sound velocity. Frozen sound velocityaf represents the velocity of a fixed mixture in the minimal disturbance propagation medium, which is defined as(83)af2≡∂p∂ρs,g.s is specific entropy.In general, the equation of state is given as(84)p=pρ,e,g.When the condensation mass fractiong is constant, the growth rate of pressure p is as follows:(85)dp=∂p∂ρe,gdρ+∂p∂eρ,gde.The first law of thermodynamics is(86)de=Tds+pρ2dρ.Substitute equation (86) into equation (85)(87)dp=∂p∂ρe,g+pρ2∂p∂eρ,gdρ+T∂p∂eρ,gds.We definep as a function of ρ, s, and g. When g is constant, the growth rate of p is(88)dp=∂p∂ρs,gdρ+∂p∂sρ,gds.Compare the increments fordρ of equations (87) and (88), and get(89)∂p∂ρs,g=∂p∂ρe,g+pρ2∂p∂eρ,g.Equation (89) is the definition of frozen sound velocity af. Substitute equations (63) and (78) into equation (89) to obtain(90)af2=RT1+Rcv0+gRv−L1.For mixed gas, when the mass fraction of condensed liquid is 0, the temperature isT0, the pressure is p0, the mass fraction of maximum condensed liquid is gmax, and the saturation is Φ0, it can be expressed as follows:(91)Φ0=pv0pv,eqT0=ρv0RvTvpv,eqT0=gmaxρ0RvT0pv,eqT0=gmaxp0Rv/R0pv,eqT0.Because saturationΦ0 is known, gmax can be obtained(92)gmax=11+Rv/Rap0/pv,eqTΦ0−1.In the expressions of condensation rateJ, droplet radius growth rate r˙, and equation of state discussed in the previous sections, there are some properties and constants of substances. For water, these properties and constants are given in the following sections.When the fluid expands at a high speed, the temperature of the condensed matter will soon drop below the three-phase point temperature. In this case, it is difficult for us to judge whether the condensed matter is liquid or solid. However, for the steady flow of wet air, the calculated results of Dohrmann et al. [29] are consistent with the experimental results under the assumption that there is only liquid condensation. Therefore, we assume that the fluid only condenses into liquid.The surface tensionσTN/m of water adopts the expression in Dohrmann et al. [29], as follows:(93)σT=19.861TTc4−42.150TTc5−5.6464×10−6,T<249.39,0.118−0.100TTc,T≥249.39.where T is Kelvin temperature, critical temperature Tc=647.3K of water, and three-phase point temperature Ttr=273.15K.The densityρl0kg/m3 of water is a function of temperature. The expression is given by Pruppacher and Klett [30]. The form at Kelvin temperature is(94)ρl0=ρ1∑i=05aiT/Ttri∑i=01biT/Ttri,T≥Ttr,∑i=02ciTTtri,T<Ttr.ρ1=999.84kg/m3 is the density of water at the three-phase point temperature, and the constants in formula (94) are given in Table 2.Table 2
The constants in formula (94).
a0=−2.0079a3=−10.447b0=−3.9603c0=0.17058a1=−3.5552a4=3.8244b1=4.9603c1=1.6353a2=13.784a5=−0.59812c2=−0.80593The pressurepv,eqN/m2 of steam in equilibrium is given by Sonntag and Heinze [31]: (95)pv,eqT=ea6+a7T+a8T2+b2lnT+c3T−1.and, a6=21.125, a7=−2.7246×10−2, a8=1.6853×10−5, b2=2.4576, and c3=−6094.4642.The latent heatLJ/kg is given by Clausius–Clapeyron equation:(96)dpv,eqdT=Lpv,eqRvT2.The expression of latent heat is as follows:(97)L=RvT21pv,eqdpv,eqdT.The derivation of equation (95) shows that the latent heat of water is(98)LT=Rva7T2+2a8T3+b2T−c3.When the temperature is between200−300K, the latent heat of equation (98) can be approximately regarded as a linear function of T:(99)LT=L0+L1T,and L0=3105913.39J/kg, and L1=−2212.97J/kg⋅K.The gas constantsR of dry air and specific heat capacity cp and cv, and specific heat γ are shown in Table 3.Table 3
Some constants of air.
RJ/kg⋅KcpJ/kg⋅KcvJ/kg⋅KγDry air287.041004.0716.961.4004Water vapor461.521859.01397.51.3302The mass of a water moleculem=2.991×10−26kg.


## 3. Numerical Methods

In Section2, the final control equation of the condensation flow model is obtained by combining the Euler equation with the classical nucleation theory and using the Hill moment method (79). In this section, the numerical model of the condensation flow is given, and the equation with the required solution is obtained by quasilinearization of the control equation. The boundary conditions and initial conditions of the equation are given. The numerical methods adopted in this paper are the classical Galerkin finite element spatial discretization scheme and the backward Euler difference time discretization scheme. For the steady problem of unstructured grid, the adaptive grid technique is also introduced to solve the problem because the Euler equation is a convection equation. The results obtained by the classical Galerkin method often have nonphysical numerical oscillations, so the stabilization method is introduced. Another important feature of the compressible Euler equation is the existence of shock waves. Because the thickness of shock waves is small and difficult to capture, artificial viscosity terms are added to better capture shock waves. The stream upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability are used to solve the decoupling problem.


### 3.1. Numerical Model

The final form of the governing equation is given in Section2. Equation (79) is the conservative form of the governing equation:(100)∂U∂t+∂F1∂x+∂F2∂y=B.Combined with the equation of state of gas proposed in Section2, each variable is as follows:(101)U=ρρuρvρetρQ0ρQ1ρQ2ρg,F1=ρuρu2+pρuvρuhtρuQ0ρuQ1ρuQ2ρug,F2=ρvρuvρv2+pρvhtρvQ0ρvQ1ρvQ2ρvg,B=0000Jr∗J+r˙ρQ0r∗2J+2r˙ρQ143πρl0r∗3J+3r˙ρQ2.Refer to Section2 for the meaning of each physical quantity.Next, equation (100) is transformed into another form from the chain relationship of derivative:(102)∂U∂t=∂U∂W∂W∂t,(103)∂F1∂x=∂F1∂U∂U∂W∂W∂x,(104)∂F2∂y=∂F2∂U∂U∂W∂W∂y.Substituting the above three equations into equation (100) can obtain(105)∂U∂W∂W∂t+∂F1∂U∂U∂W∂W∂x+∂F2∂U∂U∂W∂W∂y=B,i.e.,(106)∂W∂t+A1∂W∂x+A2∂W∂y=G,W, G is a vector, and A1=∂F1/∂U, A2=∂F2/∂U is a matrix, which can be obtained by combining the relationship between internal energy et, enthalpy ht, and pressure p introduced in Section 2,(107)W=ρuvpQ0Q1Q2g,A1=uρ0000000u01ρ000000u000000γ0p0u00000000u00000000u00000000u00000000u,A2=v0ρ000000v00000000v1ρ000000γ0pv00000000v00000000v00000000v00000000v,G=000γ0−1QJρr∗Jρ+r˙Q0r∗2Jρ+2r˙Q043πρl0r∗3Jρ+3r˙Q0.Here, γ0=1+R0−gRv/cv0+gRv+gL1, Q=−4/3πρl0cv0Rv+RvR0−L1R0p−L0ρR0−gRv2/R0−gRv2r∗3J/ρ+3r˙Q2 or recorded as (108)Wt+A⋅∇W=G,Here A=A1,A2, called Jacobian matrix.The selection of boundary conditions of Euler equation is a difficulty. Only by selecting appropriate boundary conditions can the existence and uniqueness of Euler equation solution be guaranteed. Similarly, the selection of boundary conditions is also a difficulty for the control equation of wet air condensation flow. At present, the boundary conditions are mainly divided into three types: inlet and outlet boundary, far-field boundary condition, and surface boundary conditions.Using one-dimensional fluid to discuss the inletΓin and outlet Γout boundary conditions, two-dimensional fluid is similar. Combined with the important relationship: sound velocity a=γ0p/ρ, the governing equation of one-dimensional fluid is as follows:(109)∂∂tρupQ0Q1Q2g+uρ000000u1ρ00000ρa2u0000000u0000000u0000000u0000000u∂∂xρupQ0Q1Q2g=00G4G5G6G7G8,Gii=5,6,7,8 given by (107). The characteristic polynomial obtained the Jacobian matrix of 7×7 from equation (109):(110)u−λ5u−λ2−a2=0,Seven real eigenvalues are obtained(111)λ1=u,λ2=u+a,λ3=u−a,λi=ui=4,5,6,7.It can be seen from the eigenvalues that the direction of five eigenvalues is consistent with the direction of fluid velocity streamline, and the direction of the other two eigenvalues is related to the direction of sound wave.Boundary conditions are divided into physical boundary conditions and numerical boundary conditions. Physical boundary conditions can be determined by physical data such as experimental values, and numerical boundary conditions can be obtained by numerical extrapolation in the internal flow region. The number of physical boundaries is determined by whether the fluid is supersonic or subsonic [32].For the inflow boundary conditions of one-dimensional subsonic flow0<u<a, there are six positive eigenvalues, that is, u,u+a,u,u,u,u; it indicates that the information enters the internal region from the outside, and the physical boundary conditions must be given. Another negative eigenvalue u−a indicates that the propagation of information is from inside to outside, that is, given the numerical boundary conditions. Therefore, six physical boundary conditions and one numerical boundary condition must be given. In the actual numerical calculation, the method of giving the physical boundary conditions gives the liquid moment ρQ0,ρQ1,ρQ2,ρg, total pressure p0, and total temperature T0. After giving the p0 and T0, the density ρ0 can be obtained from the gas equation of state.For supersonic flowu>a, the seven eigenvalues are positive numbers, indicating that the propagation direction of information is from outside to inside, and the physical boundary conditions must be given. Therefore, we must give the values of all fluid variables at the boundary.For the exit boundary condition of one-dimensional subsonic flow0<u<a, there are six positive eigenvalues u,u+a,u,u,u,u for all eigenvalues, the information propagates from inside to outside, and the numerical boundary condition is given. In a negative eigenvalue u−a, the information propagates from the outside to the inside, and the physical boundary conditions are given. Therefore, only one boundary condition is given at the outlet boundary conditions, and the pressure pout is given in the calculation in this paper.Similarly, for supersonic flowu>a, all characteristic directions are from inside to outside, so there is no need to give boundary conditions at the outlet boundary.At the far-fieldΓ∞ boundary, the viscous effect is negligible. In the actual numerical calculation, the far-field boundary is used to limit the size of the calculation area. Therefore, the nonreflective boundary condition needs to be introduced at the far-field boundary. In the numerical calculation in this paper, because the normal vector of velocity at the far-field boundary is very small, all eigenvalues at the inlet boundary are positive. Therefore, the values of all variables of the fluid must be given at the inlet boundary of the far-field boundary. For the far-field outlet boundary condition of subsonic velocity, the treatment method of the outlet boundary given in (1) is used.It is assumed that there is no condensation in the far field, so there is condensation in the far field.(112)g=Q2=Q1=Q0=0,x¯⟶−∞.For inviscid flow at the object surfaceΓb, the nonpenetration condition is generally used in numerical calculation; that is, the normal velocity component: u⟶⋅n=0 on the object surface is 0, where u⟶=u,v is the velocity vector and n is the external normal vector.The initial condition is the state that the field variables should satisfy at the initial time, that is, the value of each variable at the time oft=0, W|t=0=W0. When solving the steady flow problem, there is no need to give the initial condition. When solving unsteady flow problems, the initial condition is the steady solution of a given state.


### 3.2. Numerical Methods

The mathematical model of inviscid condensation flow problem is given earlier. In order to solve the practical problem, it is necessary to solve this mathematical problem and then simulate the actual condensation problem. For the partial differential equations (100), the classical Galerkin method is used to discretize them in space and the backward Euler scheme is used to discretize them in time. The grid is encrypted by grid adaptive technology.According to the discussion earlier, the mathematical model for solving the problem is obtained:(113)Wt+W=G,Ω×0,T,Win⋅n=H,Γ∞×0,T,u⟶⋅n=0,Γb×0,T,Wx⟶,0=W0x⟶onΩ,whent=0.The variables in equation (113) refer to equation (106), Ω is the calculation area, and u⟶=u,v, x⟶=x,y. Γ is the boundary of Ω, and Γ=Γb∪Γ∞, Γ=Γb∪Γ∞=ϕ.The classical Galerkin finite element method is used to discretize equation (113) in space, and the test function space is introduced(114)V=U|Ux⟶,t∈H1Ω,t∈0,T.TakeV the weight function space of ρ,u,v,p,Q0,Q1,Q2,g, respectively, multiply ∀U∈V8 at both ends of equation (113), and integrate it on the Ω,(115)∫ΩUWtdΩ+∫ΩUA⋅∇WdΩ=∫ΩUGdΩ.The second term of the above formula is obtained by using Green's formula and substituting the boundary conditions.(116)∫ΩUWtdΩ−∫Ω∇UAWdΩ=∫ΩUGdΩ−∫Γ∞UHdΓ∞.Therefore, equation (113) is transformed into a weak form problem: for ∀t∈0,T, find Wx⟶,t∈V8 such that Wx⟶,0=W0x⟶, and(117)U,Wt−∇U,AW=U,G−U,HΓ∞,∀U∈V8.Next, the space is discretized to obtain the semi-discrete scheme of the problem. The bounded regionΩ is approximated by a finite polygon Ωh, and Ωh is divided into several triangular elements Ωe, including Ω¯=∪e=1NeΩ¯e, ∩e=1NeΩe=ϕ, where Ne is the total number of elements of Ωh divided into Ωe. VΩe is the local finite element space on Ωe, which is taken as the finite element space k≥0, ∀t∈0,T of k-order polynomial.ForΩh establishing finite dimensional discrete space,(118)Vh=U∈H1Ω,U|Ωe∈VΩe,∀Ωe∈Ωh.The spaceVh is discretized with the weight function of ρ,u,v,p,Q0,Q1,Q2,g, respectively, and the Galerkin finite element space discretization of equation (117) is obtained: for time t∈0,T, find Wx⟶,t∈Vh8 such that Wx⟶,0=W0x⟶.(119)Uh,Wt−∇Uh,AW=Uh,G−Uh,HΓ∞∀Uh∈Vh8.The vectorU uses the shape function Nh of the standard finite element.In (119), only the space is discretized to obtain the semi-discrete scheme, and then the time is discretized to obtain the fully discrete scheme.The time discretization of (119) adopts the backward Euler scheme. Firstly, the notation Qn=Ω×In is introduced, where In=tn,tn+1,n=0,1,…nst−1. Thus, the unit element is transformed into Ωen=Ωe×In,e=1,2,…,Ne, and the test function space is defined as follows:(120)Vhn=U∈H1Ω,U|Ωen∈VΩen,∀Ωen∈Qn.The problem to be solved by backward Euler scheme is as follows: for∀n=0,1,…,nst−1, find Whx⟶,t∈Vhn8 such that Whx⟶,0=W0x⟶, and(121)∬QnUhWhn+1−WhnΔtdΩdt−∬Qn∇UhAWhn+1dΩdt=∬QnUhGn+1dΩdt−∬UhHn+1dΓ∞dt∀Uh∈Vhn8.where Δt=tn+1−tn, Whn+1=Whtn+1:=Wht for ∀t∈tn,tn+1.


### 3.3. Stabilization Techniques

Shock wave is a strong compression wave in the moving gas. The weak disturbance in the gas propagates around at the local sound speed. When the object moves at subsonic speed, the disturbance propagation speed is greater than that of the object, so the disturbance cannot be concentrated. At this time, the distribution of flow parameters (including velocity and pressure) on the whole flow field is continuous. When the object moves at supersonic speed, there is no time for the disturbance to reach the front of the object. As a result, the gas in front is suddenly compressed by the moving object to form a concentrated strong disturbance wave. At this time, an interface in the compression process is called shock wave.As we all know, there may be discontinuities in the solution of nonlinear hyperbolic partial differential equations. When solving aerodynamic problems with Euler equations, even if the initial distribution of the solution is smooth, there may be discontinuities in the solution with the passage of time. When solving such physical problems with second-order or higher-order accuracy schemes, in fact, when the solution of the differential equation is smooth before the discontinuity is formed, the numerical solutions of some schemes will also have nonphysical oscillations.Due to these characteristics of shock wave, many stabilization methods [34–42] are introduced. A key feature of these techniques in dealing with shock waves and equations with discontinuous solutions is to add an artificial viscosity term to the equations. The purpose of adding this term is to widen the thickness of shock waves in some elements of the grid, so as to smooth discontinuities and eliminate nonphysical shocks before shock waves occur.For general differential equations,(122)LW=f,Ω,W=h,ΓN.where L is differential operator. The stability term is added based on the standard Galerkin finite element method:(123)aW,v+∑e∫ΩePvτRWdΩ=v,f+v,hΓN︸stabilitytermv∈Vh,where aW,v is bilinear form and RW=LW−f represents the residual.WhenPv takes different expressions, it is a different stabilization method: streamline upwind/Petrov–Galerkin method.For the numerical model formula (113) proposed in this paper, the differential operator l=∂/∂t+A⋅∇ can be obtained. If Pv=A⋅∇v is taken in equation (123), only the artificial diffusion term is added in the streamline direction. Since the numerical method uses the finite element method to discretize the equation, it provides a method to add the artificial viscosity term without affecting the original equation. This method is called SUPG (streamline upwind/Petrov–Galerkin) method, which is based on the weak form of differential equation.The variational form of Galerkin method is obtained from equation (119):(124)Uh,Wt−∇Uh,AW=Uh,G−Uh,HΓ∞∀Uh∈Vh8.The SUPG method is added to the above formula to obtain the variational form after adding the stabilization term: fort∈0,T, find Wx⟶,t∈Vh8, such that Wx⟶,0=W0x⟶ and(125)Uh,Wt−∇Uh,AW+∑eA⋅∇Uh,τRWΩe=Uh,G−Uh,HΓ∞∀Uh∈Vh8,the stable term is(126)A⋅∇Uh,τA⋅∇WΩe.That is, the artificial viscosity term is added in the direction of the characteristic line. For the SUPG method of compressible flow problem, the stabilization coefficientτ is very important. The stabilization coefficient makes the artificial viscosity term added only in the direction of the characteristic line. In general, it is not easy to determine the structure of the coefficient τ, because for the two-dimensional Euler equation, the Jacobian matrix Aii=1,2 will not be diagonalized at the same time.In this paper, we take the stabilization coefficient as(127)τ=h2a,where a is the element order.When the diffusion term is small or does not exist, the convection dominated equation needs to be modified to reduce the nonphysical oscillation because the role of diffusion is becoming smaller and smaller.In this method, the artificial diffusion coefficient is added to the equation,(128)τ=hmaxa1,a22,where a1,a2 is the spectral radius of matrix A1,A2 (the maximum of the absolute value of eigenvalues).According to this method, we do not solve the problem of the original formula (113), but the modified equation:(129)Wt+A⋅∇W+∇⋅−τ∇W=GΩ×0,T.The variational form is as follows: fort∈0,T, find Wx⟶,t∈Vh8 such that Wx⟶,0=W0x⟶, and(130)Uh,Wt−∇Uh,AW+∑e∇Uh,τ∇WΩe=Uh,G−Uh,HΓ∞∀Uh∈Vh8.The added artificial diffusion term can reduce the numerical oscillation and the diffusion of oscillation to other areas.According to the previous discussion of numerical method and stabilization method, combined with the governing equation of condensation flow, the weak form of the modified differential equation is as follows: for timet∈0,T, find Whx⟶,t∈Vh8 such that Whx⟶,0=W0x⟶, and(131)∫Ω∂Wh∂t+A1∂Wh∂x+A2∂Wh∂y⋅NhdΩ+∑e=1Ne∫ΩτSUPG⋅∂Wh∂t+A1∂Wh∂x+A2∂Wh∂y−G⋅A1∂Nh∂x+A2∂Nh∂ydΩ+∑e=1Ne∫Ωτiso⋅∇Wh⋅∇NhdΩ=∑e=1Ne∫ΩGNhdΩ−∫ΓWhHdΓ∀Nh∈Vh8.In the above formula, the first term is the classical Galerkin finite element discretization scheme, the second term is the perturbation term of SUPG method in the streamline direction, and the third term is the artificial viscosity term added by various isotropic diffusion methods in all directions, where(132)τSUPG=h2a,τiso=hmaxa1,a22,where a is the element order and a1,a2 is the spectral radius of matrix A1,A2 (the maximum of the absolute value of eigenvalues).


### 3.4. Grid Adaptive Technology

In numerical calculation, some special regions, such as large deformation, boundary layer, and shock surface, have large singularity, and the solution results are often poor. In order to improve the calculation accuracy of these regions, local mesh encryption technology is introduced to continuously adjust and refine the mesh of these special regions in the iterative process. It is called grid adaptive technology to strive to match the grid density with the physical solution, so as to improve the accuracy of numerical solution [33]. Especially in the field of aeronautics, the shock surface will be generated due to the high-speed expansion flow of air. In order to improve the resolution of shock wave, the grid adaptive technology is introduced in the calculation of steady condensation flow. The grid in the area where the shock surface exists is locally refined, and the sparse grid is used in the area where the solution changes relatively slowly. Experimental results show that the method of introducing grid adaptive technology is better than that without this method.


### 3.5. Section

In this chapter, the mathematical model and numerical method for solving the transonic flow of wet air are proposed, which are considered from the following aspects:According to the discussion in Chapter 2, the mathematical model of the practical problem is obtained, and the boundary conditions and initial conditions of the problem are given, which completely transforms the practical problem into the solution of a system of partial differential equations.For the solution of partial differential equations, the classical Galerkin finite element method is used to discretize the equations in space, and the semi-discrete scheme is obtained.The weak form is discretized in space, and then the backward Euler scheme is used to discretize time to obtain the fully discrete scheme.In the process of solving Euler equations, due to the existence of shock, the numerical solution produces nonphysical shock. In order to better capture the shock and obtain better numerical results, the streamline upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability are introduced here.For the steady problem, the mesh adaptive technology is also introduced to encrypt the local mesh, so that the numerical solution can be better understood.


## 3.1. Numerical Model

The final form of the governing equation is given in Section2. Equation (79) is the conservative form of the governing equation:(100)∂U∂t+∂F1∂x+∂F2∂y=B.Combined with the equation of state of gas proposed in Section2, each variable is as follows:(101)U=ρρuρvρetρQ0ρQ1ρQ2ρg,F1=ρuρu2+pρuvρuhtρuQ0ρuQ1ρuQ2ρug,F2=ρvρuvρv2+pρvhtρvQ0ρvQ1ρvQ2ρvg,B=0000Jr∗J+r˙ρQ0r∗2J+2r˙ρQ143πρl0r∗3J+3r˙ρQ2.Refer to Section2 for the meaning of each physical quantity.Next, equation (100) is transformed into another form from the chain relationship of derivative:(102)∂U∂t=∂U∂W∂W∂t,(103)∂F1∂x=∂F1∂U∂U∂W∂W∂x,(104)∂F2∂y=∂F2∂U∂U∂W∂W∂y.Substituting the above three equations into equation (100) can obtain(105)∂U∂W∂W∂t+∂F1∂U∂U∂W∂W∂x+∂F2∂U∂U∂W∂W∂y=B,i.e.,(106)∂W∂t+A1∂W∂x+A2∂W∂y=G,W, G is a vector, and A1=∂F1/∂U, A2=∂F2/∂U is a matrix, which can be obtained by combining the relationship between internal energy et, enthalpy ht, and pressure p introduced in Section 2,(107)W=ρuvpQ0Q1Q2g,A1=uρ0000000u01ρ000000u000000γ0p0u00000000u00000000u00000000u00000000u,A2=v0ρ000000v00000000v1ρ000000γ0pv00000000v00000000v00000000v00000000v,G=000γ0−1QJρr∗Jρ+r˙Q0r∗2Jρ+2r˙Q043πρl0r∗3Jρ+3r˙Q0.Here, γ0=1+R0−gRv/cv0+gRv+gL1, Q=−4/3πρl0cv0Rv+RvR0−L1R0p−L0ρR0−gRv2/R0−gRv2r∗3J/ρ+3r˙Q2 or recorded as (108)Wt+A⋅∇W=G,Here A=A1,A2, called Jacobian matrix.The selection of boundary conditions of Euler equation is a difficulty. Only by selecting appropriate boundary conditions can the existence and uniqueness of Euler equation solution be guaranteed. Similarly, the selection of boundary conditions is also a difficulty for the control equation of wet air condensation flow. At present, the boundary conditions are mainly divided into three types: inlet and outlet boundary, far-field boundary condition, and surface boundary conditions.Using one-dimensional fluid to discuss the inletΓin and outlet Γout boundary conditions, two-dimensional fluid is similar. Combined with the important relationship: sound velocity a=γ0p/ρ, the governing equation of one-dimensional fluid is as follows:(109)∂∂tρupQ0Q1Q2g+uρ000000u1ρ00000ρa2u0000000u0000000u0000000u0000000u∂∂xρupQ0Q1Q2g=00G4G5G6G7G8,Gii=5,6,7,8 given by (107). The characteristic polynomial obtained the Jacobian matrix of 7×7 from equation (109):(110)u−λ5u−λ2−a2=0,Seven real eigenvalues are obtained(111)λ1=u,λ2=u+a,λ3=u−a,λi=ui=4,5,6,7.It can be seen from the eigenvalues that the direction of five eigenvalues is consistent with the direction of fluid velocity streamline, and the direction of the other two eigenvalues is related to the direction of sound wave.Boundary conditions are divided into physical boundary conditions and numerical boundary conditions. Physical boundary conditions can be determined by physical data such as experimental values, and numerical boundary conditions can be obtained by numerical extrapolation in the internal flow region. The number of physical boundaries is determined by whether the fluid is supersonic or subsonic [32].For the inflow boundary conditions of one-dimensional subsonic flow0<u<a, there are six positive eigenvalues, that is, u,u+a,u,u,u,u; it indicates that the information enters the internal region from the outside, and the physical boundary conditions must be given. Another negative eigenvalue u−a indicates that the propagation of information is from inside to outside, that is, given the numerical boundary conditions. Therefore, six physical boundary conditions and one numerical boundary condition must be given. In the actual numerical calculation, the method of giving the physical boundary conditions gives the liquid moment ρQ0,ρQ1,ρQ2,ρg, total pressure p0, and total temperature T0. After giving the p0 and T0, the density ρ0 can be obtained from the gas equation of state.For supersonic flowu>a, the seven eigenvalues are positive numbers, indicating that the propagation direction of information is from outside to inside, and the physical boundary conditions must be given. Therefore, we must give the values of all fluid variables at the boundary.For the exit boundary condition of one-dimensional subsonic flow0<u<a, there are six positive eigenvalues u,u+a,u,u,u,u for all eigenvalues, the information propagates from inside to outside, and the numerical boundary condition is given. In a negative eigenvalue u−a, the information propagates from the outside to the inside, and the physical boundary conditions are given. Therefore, only one boundary condition is given at the outlet boundary conditions, and the pressure pout is given in the calculation in this paper.Similarly, for supersonic flowu>a, all characteristic directions are from inside to outside, so there is no need to give boundary conditions at the outlet boundary.At the far-fieldΓ∞ boundary, the viscous effect is negligible. In the actual numerical calculation, the far-field boundary is used to limit the size of the calculation area. Therefore, the nonreflective boundary condition needs to be introduced at the far-field boundary. In the numerical calculation in this paper, because the normal vector of velocity at the far-field boundary is very small, all eigenvalues at the inlet boundary are positive. Therefore, the values of all variables of the fluid must be given at the inlet boundary of the far-field boundary. For the far-field outlet boundary condition of subsonic velocity, the treatment method of the outlet boundary given in (1) is used.It is assumed that there is no condensation in the far field, so there is condensation in the far field.(112)g=Q2=Q1=Q0=0,x¯⟶−∞.For inviscid flow at the object surfaceΓb, the nonpenetration condition is generally used in numerical calculation; that is, the normal velocity component: u⟶⋅n=0 on the object surface is 0, where u⟶=u,v is the velocity vector and n is the external normal vector.The initial condition is the state that the field variables should satisfy at the initial time, that is, the value of each variable at the time oft=0, W|t=0=W0. When solving the steady flow problem, there is no need to give the initial condition. When solving unsteady flow problems, the initial condition is the steady solution of a given state.


## 3.2. Numerical Methods

The mathematical model of inviscid condensation flow problem is given earlier. In order to solve the practical problem, it is necessary to solve this mathematical problem and then simulate the actual condensation problem. For the partial differential equations (100), the classical Galerkin method is used to discretize them in space and the backward Euler scheme is used to discretize them in time. The grid is encrypted by grid adaptive technology.According to the discussion earlier, the mathematical model for solving the problem is obtained:(113)Wt+W=G,Ω×0,T,Win⋅n=H,Γ∞×0,T,u⟶⋅n=0,Γb×0,T,Wx⟶,0=W0x⟶onΩ,whent=0.The variables in equation (113) refer to equation (106), Ω is the calculation area, and u⟶=u,v, x⟶=x,y. Γ is the boundary of Ω, and Γ=Γb∪Γ∞, Γ=Γb∪Γ∞=ϕ.The classical Galerkin finite element method is used to discretize equation (113) in space, and the test function space is introduced(114)V=U|Ux⟶,t∈H1Ω,t∈0,T.TakeV the weight function space of ρ,u,v,p,Q0,Q1,Q2,g, respectively, multiply ∀U∈V8 at both ends of equation (113), and integrate it on the Ω,(115)∫ΩUWtdΩ+∫ΩUA⋅∇WdΩ=∫ΩUGdΩ.The second term of the above formula is obtained by using Green's formula and substituting the boundary conditions.(116)∫ΩUWtdΩ−∫Ω∇UAWdΩ=∫ΩUGdΩ−∫Γ∞UHdΓ∞.Therefore, equation (113) is transformed into a weak form problem: for ∀t∈0,T, find Wx⟶,t∈V8 such that Wx⟶,0=W0x⟶, and(117)U,Wt−∇U,AW=U,G−U,HΓ∞,∀U∈V8.Next, the space is discretized to obtain the semi-discrete scheme of the problem. The bounded regionΩ is approximated by a finite polygon Ωh, and Ωh is divided into several triangular elements Ωe, including Ω¯=∪e=1NeΩ¯e, ∩e=1NeΩe=ϕ, where Ne is the total number of elements of Ωh divided into Ωe. VΩe is the local finite element space on Ωe, which is taken as the finite element space k≥0, ∀t∈0,T of k-order polynomial.ForΩh establishing finite dimensional discrete space,(118)Vh=U∈H1Ω,U|Ωe∈VΩe,∀Ωe∈Ωh.The spaceVh is discretized with the weight function of ρ,u,v,p,Q0,Q1,Q2,g, respectively, and the Galerkin finite element space discretization of equation (117) is obtained: for time t∈0,T, find Wx⟶,t∈Vh8 such that Wx⟶,0=W0x⟶.(119)Uh,Wt−∇Uh,AW=Uh,G−Uh,HΓ∞∀Uh∈Vh8.The vectorU uses the shape function Nh of the standard finite element.In (119), only the space is discretized to obtain the semi-discrete scheme, and then the time is discretized to obtain the fully discrete scheme.The time discretization of (119) adopts the backward Euler scheme. Firstly, the notation Qn=Ω×In is introduced, where In=tn,tn+1,n=0,1,…nst−1. Thus, the unit element is transformed into Ωen=Ωe×In,e=1,2,…,Ne, and the test function space is defined as follows:(120)Vhn=U∈H1Ω,U|Ωen∈VΩen,∀Ωen∈Qn.The problem to be solved by backward Euler scheme is as follows: for∀n=0,1,…,nst−1, find Whx⟶,t∈Vhn8 such that Whx⟶,0=W0x⟶, and(121)∬QnUhWhn+1−WhnΔtdΩdt−∬Qn∇UhAWhn+1dΩdt=∬QnUhGn+1dΩdt−∬UhHn+1dΓ∞dt∀Uh∈Vhn8.where Δt=tn+1−tn, Whn+1=Whtn+1:=Wht for ∀t∈tn,tn+1.


## 3.3. Stabilization Techniques

Shock wave is a strong compression wave in the moving gas. The weak disturbance in the gas propagates around at the local sound speed. When the object moves at subsonic speed, the disturbance propagation speed is greater than that of the object, so the disturbance cannot be concentrated. At this time, the distribution of flow parameters (including velocity and pressure) on the whole flow field is continuous. When the object moves at supersonic speed, there is no time for the disturbance to reach the front of the object. As a result, the gas in front is suddenly compressed by the moving object to form a concentrated strong disturbance wave. At this time, an interface in the compression process is called shock wave.As we all know, there may be discontinuities in the solution of nonlinear hyperbolic partial differential equations. When solving aerodynamic problems with Euler equations, even if the initial distribution of the solution is smooth, there may be discontinuities in the solution with the passage of time. When solving such physical problems with second-order or higher-order accuracy schemes, in fact, when the solution of the differential equation is smooth before the discontinuity is formed, the numerical solutions of some schemes will also have nonphysical oscillations.Due to these characteristics of shock wave, many stabilization methods [34–42] are introduced. A key feature of these techniques in dealing with shock waves and equations with discontinuous solutions is to add an artificial viscosity term to the equations. The purpose of adding this term is to widen the thickness of shock waves in some elements of the grid, so as to smooth discontinuities and eliminate nonphysical shocks before shock waves occur.For general differential equations,(122)LW=f,Ω,W=h,ΓN.where L is differential operator. The stability term is added based on the standard Galerkin finite element method:(123)aW,v+∑e∫ΩePvτRWdΩ=v,f+v,hΓN︸stabilitytermv∈Vh,where aW,v is bilinear form and RW=LW−f represents the residual.WhenPv takes different expressions, it is a different stabilization method: streamline upwind/Petrov–Galerkin method.For the numerical model formula (113) proposed in this paper, the differential operator l=∂/∂t+A⋅∇ can be obtained. If Pv=A⋅∇v is taken in equation (123), only the artificial diffusion term is added in the streamline direction. Since the numerical method uses the finite element method to discretize the equation, it provides a method to add the artificial viscosity term without affecting the original equation. This method is called SUPG (streamline upwind/Petrov–Galerkin) method, which is based on the weak form of differential equation.The variational form of Galerkin method is obtained from equation (119):(124)Uh,Wt−∇Uh,AW=Uh,G−Uh,HΓ∞∀Uh∈Vh8.The SUPG method is added to the above formula to obtain the variational form after adding the stabilization term: fort∈0,T, find Wx⟶,t∈Vh8, such that Wx⟶,0=W0x⟶ and(125)Uh,Wt−∇Uh,AW+∑eA⋅∇Uh,τRWΩe=Uh,G−Uh,HΓ∞∀Uh∈Vh8,the stable term is(126)A⋅∇Uh,τA⋅∇WΩe.That is, the artificial viscosity term is added in the direction of the characteristic line. For the SUPG method of compressible flow problem, the stabilization coefficientτ is very important. The stabilization coefficient makes the artificial viscosity term added only in the direction of the characteristic line. In general, it is not easy to determine the structure of the coefficient τ, because for the two-dimensional Euler equation, the Jacobian matrix Aii=1,2 will not be diagonalized at the same time.In this paper, we take the stabilization coefficient as(127)τ=h2a,where a is the element order.When the diffusion term is small or does not exist, the convection dominated equation needs to be modified to reduce the nonphysical oscillation because the role of diffusion is becoming smaller and smaller.In this method, the artificial diffusion coefficient is added to the equation,(128)τ=hmaxa1,a22,where a1,a2 is the spectral radius of matrix A1,A2 (the maximum of the absolute value of eigenvalues).According to this method, we do not solve the problem of the original formula (113), but the modified equation:(129)Wt+A⋅∇W+∇⋅−τ∇W=GΩ×0,T.The variational form is as follows: fort∈0,T, find Wx⟶,t∈Vh8 such that Wx⟶,0=W0x⟶, and(130)Uh,Wt−∇Uh,AW+∑e∇Uh,τ∇WΩe=Uh,G−Uh,HΓ∞∀Uh∈Vh8.The added artificial diffusion term can reduce the numerical oscillation and the diffusion of oscillation to other areas.According to the previous discussion of numerical method and stabilization method, combined with the governing equation of condensation flow, the weak form of the modified differential equation is as follows: for timet∈0,T, find Whx⟶,t∈Vh8 such that Whx⟶,0=W0x⟶, and(131)∫Ω∂Wh∂t+A1∂Wh∂x+A2∂Wh∂y⋅NhdΩ+∑e=1Ne∫ΩτSUPG⋅∂Wh∂t+A1∂Wh∂x+A2∂Wh∂y−G⋅A1∂Nh∂x+A2∂Nh∂ydΩ+∑e=1Ne∫Ωτiso⋅∇Wh⋅∇NhdΩ=∑e=1Ne∫ΩGNhdΩ−∫ΓWhHdΓ∀Nh∈Vh8.In the above formula, the first term is the classical Galerkin finite element discretization scheme, the second term is the perturbation term of SUPG method in the streamline direction, and the third term is the artificial viscosity term added by various isotropic diffusion methods in all directions, where(132)τSUPG=h2a,τiso=hmaxa1,a22,where a is the element order and a1,a2 is the spectral radius of matrix A1,A2 (the maximum of the absolute value of eigenvalues).


## 3.4. Grid Adaptive Technology

In numerical calculation, some special regions, such as large deformation, boundary layer, and shock surface, have large singularity, and the solution results are often poor. In order to improve the calculation accuracy of these regions, local mesh encryption technology is introduced to continuously adjust and refine the mesh of these special regions in the iterative process. It is called grid adaptive technology to strive to match the grid density with the physical solution, so as to improve the accuracy of numerical solution [33]. Especially in the field of aeronautics, the shock surface will be generated due to the high-speed expansion flow of air. In order to improve the resolution of shock wave, the grid adaptive technology is introduced in the calculation of steady condensation flow. The grid in the area where the shock surface exists is locally refined, and the sparse grid is used in the area where the solution changes relatively slowly. Experimental results show that the method of introducing grid adaptive technology is better than that without this method.


## 3.5. Section

In this chapter, the mathematical model and numerical method for solving the transonic flow of wet air are proposed, which are considered from the following aspects:According to the discussion in Chapter 2, the mathematical model of the practical problem is obtained, and the boundary conditions and initial conditions of the problem are given, which completely transforms the practical problem into the solution of a system of partial differential equations.For the solution of partial differential equations, the classical Galerkin finite element method is used to discretize the equations in space, and the semi-discrete scheme is obtained.The weak form is discretized in space, and then the backward Euler scheme is used to discretize time to obtain the fully discrete scheme.In the process of solving Euler equations, due to the existence of shock, the numerical solution produces nonphysical shock. In order to better capture the shock and obtain better numerical results, the streamline upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability are introduced here.For the steady problem, the mesh adaptive technology is also introduced to encrypt the local mesh, so that the numerical solution can be better understood.


## 4. Numerical Simulation

In Chapter 3, a numerical method for solving the coupled equations is proposed, which requires numerical simulation to verify the feasibility and stability of the method. At the same time, the numerical simulation of dry air flow and wet air flow in parallel nozzle A1 further verifies the feasibility of the numerical method and the condensation of wet air in the nozzle.


### 4.1. NACA-0012 Airfoil

For NACA-0012 airfoil, its ordinate can be expressed as(133)y=±0.60.2969x−0.126x−0.3516x2+0.2843x3−0.1015x4,where x is the distance to the leading edge of the wing, 1 is chord length c, and y is the dimensionless thickness. It can be seen from the above discussion that NACA-0012 is shown in Figure 2.Figure 2
NACA-0012 airfoil.The following will be a numerical simulation of the air flow around NACA-0012 wet air.Example 1.
The airfoil is NACA-0012, and the free flow conditions given areM∞=1.4, p∞=61060Pa, T∞=253.8K, and relative humidity Φ0=50%. When dealing with the numerical simulation of steady wet air flow field, adaptive grid technology is added, and the numerical results are shown in the figure below.
Figures3–5 are the contour maps of Mach number M, pressure p, and density ρ obtained without adaptive grid technology and with adaptive grid technology, respectively. Figure 6 shows the grids without adaptive grid technology and those with adaptive grid technology. From the results, it can be seen that the adaptive grid technology encrypts the grid in the area where the shock exists, while the sparse grid is used in the area where the solution is relatively flat. The original numerical results are modified, and better results are obtained.Figure 3
Contour map of Mach numberM obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 4
Contour map of pressurep obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 5
Contour map of densityρ obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 6
Original grid (a) and grid of grid adaptive technology (b).
(a)(b)Example 2.
The airfoil is NACA-0012, and the free flow condition is given asp∞=65600Pa, T∞=259K, and Φ0=50%. The numerical results obtained when m is taken as 0.8 and 1.4, respectively, are given below.
Figures7–9 are the contour maps of Mach number M, pressure p, and density ρ with Mach numbers M∞=0.8 and M∞=1.4. It can be seen from the figure that the position of the shock wave moves downstream with the increase of Mach number, and the stability of the method is also verified.Figure 7
Contour map of Mach numberM of M∞=0.8 (a) and M∞=1.4 (b).
(a)(b)Figure 8
Contour map of pressurep of M∞=0.8 (a) and M∞=1.4 (b).
(a)(b)Figure 9
Contour map of densityρ of M∞=0.8 (a) and M∞=1.4 (b).
(a)(b)Example 3.
the airfoil is NACA-0012, and the free flow condition is given asM∞=2, p∞=65600Pa, T∞=259K, and Φ0=50%; the numerical results obtained when the velocity angle of attack α is 0° and 1.25°, respectively, are given below.
Figures10–12 are the contour maps of Mach number M, pressure p, and density ρ of angles of attack α = 0° and α = 1.25°, respectively. From the results in the figure, it can be seen that the change of angle of attack has a greater impact on the upstream of the convection field and less impact on the downstream of the flow field.Figure 10
Contour map of Mach numberM of α = 0° (a) and α = 1.25° (b).
(a)(b)Figure 11
Contour map of pressurep of α = 0° (a) and α = 1.25° (b).
(a)(b)Figure 12
Contour map of densityρ of α = 0° (a) and α = 1.25° (b).
(a)(b)Example 4.
Moist air flows through NACA-0012 airfoil, and the free flow condition isM∞=0.8, p∞=65600Pa, and T∞=259K. When the relative humidity is Φ0=0%, Φ0=50%, and Φ0=60%, the numerical results are as follows.
Figure13 shows a comparison diagram of differential pressure resistance cp obtained when the relative humidity is Φ0=0%, Φ0=50%, and Φ0=60%, respectively, and gives the results obtained by adding stabilization technology according to the classical finite element method proposed in Chapter 3. It can be seen from the figure that the pressure increases in the supersonic region. With the increase of humidity, the position of the first shock wave moves forward. When the pressure decreases, the shock force decreases. The relative humidity is 50% mainly discussed, and the following results are obtained.
Figure14 shows the isoline diagram of condensation mass fraction. It can be seen from the figure that condensation only occurs at the trailing edge of the airfoil. It is consistent with the region where the shock wave is generated when the relative humidity is 50% in Figure 13. The shock wave is generally generated in the supersonic region, so the condensation of wet air around the airfoil also occurs in the supersonic region.Figure 13
Differential pressure resistancecp under different humidity conditions.Figure 14
Isoline ofg/gmax when relative humidity is 50%.


### 4.2. Parallel-Jet Nozzle A1

There is little difference between the downstream part of parallel tail nozzle A1 and the throat part, and the nozzle has a straight outlet part. Therefore, when the fluid passes through the nozzle, the fluid is very sensitive to the change of pressure and the composition of incoming wet steam. The geometric model of the nozzle is given in Figure15.Figure 15
Geometric model of parallel tail nozzle A1.The nozzle has a straight outflow part with a length of about104.5mm and the throat part of the nozzle at the downstream 87.6mm, its height is 2y∗=90mm, and the shape of the arc is a circular arc with a radius R∗=300mm. The height of the outflow part of the nozzle is 2y∗≈92.69mm, so the Mach number of isentropic outflow is Mex=1.2. See Adam [43] for details.In the wind tunnel experiment, the parallel outflow part of the nozzle is connected with the adjustable diffuser, which can control the outflow pressure. The experimental device is shown in Figure16. In this way, the fluid passing through the nozzle completely depends on the geometry of the nozzle, not the pressure at the parallel outflow boundary of the nozzle. In the numerical simulation, the geometric model we use is shown in Figure 17.Figure 16
Structural diagram of supersonic wind tunnel [43].Figure 17
Geometric model of nozzle in numerical simulation.The nozzle flow of dry air and wet air will be numerically simulated below. The same boundary conditions are given for dry air and wet air, and the incoming boundary isT0=298.7K and p0=1.004bar. For the initial saturation (relative humidity) Φ0=35.6% of wet air, the maximum mass fraction can be obtained gmax=7.253g/kg from equation (92). For dry air, the forced boundary condition p=0.6339bar is given at the outlet boundary.Example 4-5: The forced boundary conditionp=0.6339bar is given at the outlet boundary of dry air. According to the given boundary conditions, the nozzle flow of dry air generates shock wave due to the different pressures of inlet and outlet. Moreover, since the given inlet and outlet boundary conditions are subsonic boundary conditions, the Mach number M=0.7 is taken.In the process of solving, different grid numbers, degrees of freedom used, and CPU time consumed are given in Table4.Table 4
Comparison of degrees of freedom and CPU when taking different grids.
Grid numberDegrees of freedomCPU (s)516123255146432308220644522114585612310612825617294726According to the above different grid divisions, the change of pressurep/p0 is shown in Figure 18.Figure 18
Contour map of Mach numberM (a) and pressure p/p0 (b). When the grid is 516.
(a)(b)As can be seen from Figure19, with the increase of the number of grids, the initial position of the shock wave moves backward, but changes little, which shows the stability of the numerical method.Figure 19
Values of pressurep/p0 obtained with different grids.When the grid numbers are 516, 2064, and 5856, the obtained values of Mach numberM and pressure p/p0 are shown in Figures 18, 20, and 21.Figure 20
Contour map of Mach numberM (a) and pressure p/p0 (b). When the grid is 2064.
(a)(b)Figure 21
Contour map of Mach numberM (a) and pressure p/p0 (b). When the grid number is 5856.
(a)(b)The comparison of the values ofp/p0 on the right also verifies the conclusion in Figure 19. With the increase of the number of grids, the compression wave moves to the right, while the comparison of the values of M on the left shows that the mesh encryption makes the solution more and more smooth.Example 6.
For the steady problem of wet air, add adaptive grid technology and take Mach numberM=4, Φ0=30%. The results are shown in the following figure.
Figures22–24 are the contour maps of Mach number M, pressure p, and density ρ obtained without adaptive grid technology and with adaptive grid technology, respectively. Figure 25 shows the grid without adaptive grid technology and the grid with adaptive grid technology. It can be seen from the comparison diagram of the above results that the results obtained by using the original grid to solve the wet air nozzle flow are not smooth, but the grid obtained by modifying the original grid with adaptive grid technology. The obtained solution is smooth, and the existence of shock wave is more obvious. This shows that the mesh adaptive technology is an effective numerical method in solving unsteady problems. It not only makes the solution smooth, but also can encrypt the local mesh in the region where the shock wave exists, which is more helpful to capture the shock wave.Figure 22
Contour map of Mach numberM obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 23
Contour map of densityρ obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 24
Contour map of pressurep obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 25
Original grid (a) and grid of grid adaptive technology (b).
(a)(b)Example 7.
Take the result of dry air pressure contraction wave as the initial condition of wet air nozzle flow, saturationΦ0=35.6%, the boundary condition is supersonic boundary condition, and the Mach numbers M=1.4 and M=2.5 of the incoming flow are taken. The numerical results are shown in the figure below.
Figures26–28 show the contour maps of Mach number M, pressure p, and density ρ with Mach numbers M=1.4 and M=2.5. It can be seen from the results that the Mach number of the whole flow field changes greatly with the increase of Mach number, while the contour maps of condensation mass fraction and condensation rate change little in the whole flow field. The condensation of wet air occurs in the throat of the nozzle. The linear comparison diagram of g/gmax and log10 J at different Mach numbers is obtained from the above contour map.
Figure29 shows g/gmax linear comparison diagram when the Mach number is 1.4 and 2.5, respectively. It can be seen that when the Mach number is 1.4, condensation occurs at the throat of the nozzle, while when the Mach number is 2.5, condensation exists in the whole flow field, but it is still obvious at the throat. When the fluid reaches the throat, condensation occurs. Next, due to the shock force at the throat, some droplets evaporate, making the condensation mass fraction continuously decrease. The subsequent expansion increases g/gmax.
Figure30 shows log10J linear comparison diagram when the Mach number is 1.4 and 2.5, respectively, indicating that the condensation rate increases with the increase of Mach number. The change of condensation rate and condensation mass fraction is synchronous. When the condensation rate reaches the maximum, the condensation mass fraction also reaches the maximum. The condensation rate is reduced to zero due to the impact of shock wave.Figure 26
Contour map of Mach numberM for M=1.4 (a) and M=2.5 (b).
(a)(b)Figure 27
Contour map ofg/gmax for M=1.4 (a) and M=2.5 (b).
(a)(b)Figure 28
Contour map of log10J for M=1.4 (a) and M=2.5 (b).
(a)(b)Figure 29
Linear comparison diagram ofg/gmax with different Mach numbers.Figure 30
Linear comparison diagram of log10J with different Mach numbers.


### 4.3. Conclusions

In this chapter, the stability and feasibility of the numerical method proposed in Chapter 3 are verified by the numerical simulation experiments of transonic flow and supersonic flow around the wing and parallel tail nozzle. The experimental results show that the stabilization method proposed in Chapter 3 has good stability in solving practical problems.The mesh adaptive technology is added to the steady problem, which better modifies the numerical solution of the original method and has a good ability to capture the shock wave.The numerical simulation of wet air flow around the wing shows that condensation occurs in the supersonic region and the shock wave will affect the condensation.The numerical simulation of wet air parallel wake nozzle flow shows that condensation occurs in the throat of the nozzle, which has certain application value.


## 4.1. NACA-0012 Airfoil

For NACA-0012 airfoil, its ordinate can be expressed as(133)y=±0.60.2969x−0.126x−0.3516x2+0.2843x3−0.1015x4,where x is the distance to the leading edge of the wing, 1 is chord length c, and y is the dimensionless thickness. It can be seen from the above discussion that NACA-0012 is shown in Figure 2.Figure 2
NACA-0012 airfoil.The following will be a numerical simulation of the air flow around NACA-0012 wet air.Example 1.
The airfoil is NACA-0012, and the free flow conditions given areM∞=1.4, p∞=61060Pa, T∞=253.8K, and relative humidity Φ0=50%. When dealing with the numerical simulation of steady wet air flow field, adaptive grid technology is added, and the numerical results are shown in the figure below.
Figures3–5 are the contour maps of Mach number M, pressure p, and density ρ obtained without adaptive grid technology and with adaptive grid technology, respectively. Figure 6 shows the grids without adaptive grid technology and those with adaptive grid technology. From the results, it can be seen that the adaptive grid technology encrypts the grid in the area where the shock exists, while the sparse grid is used in the area where the solution is relatively flat. The original numerical results are modified, and better results are obtained.Figure 3
Contour map of Mach numberM obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 4
Contour map of pressurep obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 5
Contour map of densityρ obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 6
Original grid (a) and grid of grid adaptive technology (b).
(a)(b)Example 2.
The airfoil is NACA-0012, and the free flow condition is given asp∞=65600Pa, T∞=259K, and Φ0=50%. The numerical results obtained when m is taken as 0.8 and 1.4, respectively, are given below.
Figures7–9 are the contour maps of Mach number M, pressure p, and density ρ with Mach numbers M∞=0.8 and M∞=1.4. It can be seen from the figure that the position of the shock wave moves downstream with the increase of Mach number, and the stability of the method is also verified.Figure 7
Contour map of Mach numberM of M∞=0.8 (a) and M∞=1.4 (b).
(a)(b)Figure 8
Contour map of pressurep of M∞=0.8 (a) and M∞=1.4 (b).
(a)(b)Figure 9
Contour map of densityρ of M∞=0.8 (a) and M∞=1.4 (b).
(a)(b)Example 3.
the airfoil is NACA-0012, and the free flow condition is given asM∞=2, p∞=65600Pa, T∞=259K, and Φ0=50%; the numerical results obtained when the velocity angle of attack α is 0° and 1.25°, respectively, are given below.
Figures10–12 are the contour maps of Mach number M, pressure p, and density ρ of angles of attack α = 0° and α = 1.25°, respectively. From the results in the figure, it can be seen that the change of angle of attack has a greater impact on the upstream of the convection field and less impact on the downstream of the flow field.Figure 10
Contour map of Mach numberM of α = 0° (a) and α = 1.25° (b).
(a)(b)Figure 11
Contour map of pressurep of α = 0° (a) and α = 1.25° (b).
(a)(b)Figure 12
Contour map of densityρ of α = 0° (a) and α = 1.25° (b).
(a)(b)Example 4.
Moist air flows through NACA-0012 airfoil, and the free flow condition isM∞=0.8, p∞=65600Pa, and T∞=259K. When the relative humidity is Φ0=0%, Φ0=50%, and Φ0=60%, the numerical results are as follows.
Figure13 shows a comparison diagram of differential pressure resistance cp obtained when the relative humidity is Φ0=0%, Φ0=50%, and Φ0=60%, respectively, and gives the results obtained by adding stabilization technology according to the classical finite element method proposed in Chapter 3. It can be seen from the figure that the pressure increases in the supersonic region. With the increase of humidity, the position of the first shock wave moves forward. When the pressure decreases, the shock force decreases. The relative humidity is 50% mainly discussed, and the following results are obtained.
Figure14 shows the isoline diagram of condensation mass fraction. It can be seen from the figure that condensation only occurs at the trailing edge of the airfoil. It is consistent with the region where the shock wave is generated when the relative humidity is 50% in Figure 13. The shock wave is generally generated in the supersonic region, so the condensation of wet air around the airfoil also occurs in the supersonic region.Figure 13
Differential pressure resistancecp under different humidity conditions.Figure 14
Isoline ofg/gmax when relative humidity is 50%.


## 4.2. Parallel-Jet Nozzle A1

There is little difference between the downstream part of parallel tail nozzle A1 and the throat part, and the nozzle has a straight outlet part. Therefore, when the fluid passes through the nozzle, the fluid is very sensitive to the change of pressure and the composition of incoming wet steam. The geometric model of the nozzle is given in Figure15.Figure 15
Geometric model of parallel tail nozzle A1.The nozzle has a straight outflow part with a length of about104.5mm and the throat part of the nozzle at the downstream 87.6mm, its height is 2y∗=90mm, and the shape of the arc is a circular arc with a radius R∗=300mm. The height of the outflow part of the nozzle is 2y∗≈92.69mm, so the Mach number of isentropic outflow is Mex=1.2. See Adam [43] for details.In the wind tunnel experiment, the parallel outflow part of the nozzle is connected with the adjustable diffuser, which can control the outflow pressure. The experimental device is shown in Figure16. In this way, the fluid passing through the nozzle completely depends on the geometry of the nozzle, not the pressure at the parallel outflow boundary of the nozzle. In the numerical simulation, the geometric model we use is shown in Figure 17.Figure 16
Structural diagram of supersonic wind tunnel [43].Figure 17
Geometric model of nozzle in numerical simulation.The nozzle flow of dry air and wet air will be numerically simulated below. The same boundary conditions are given for dry air and wet air, and the incoming boundary isT0=298.7K and p0=1.004bar. For the initial saturation (relative humidity) Φ0=35.6% of wet air, the maximum mass fraction can be obtained gmax=7.253g/kg from equation (92). For dry air, the forced boundary condition p=0.6339bar is given at the outlet boundary.Example 4-5: The forced boundary conditionp=0.6339bar is given at the outlet boundary of dry air. According to the given boundary conditions, the nozzle flow of dry air generates shock wave due to the different pressures of inlet and outlet. Moreover, since the given inlet and outlet boundary conditions are subsonic boundary conditions, the Mach number M=0.7 is taken.In the process of solving, different grid numbers, degrees of freedom used, and CPU time consumed are given in Table4.Table 4
Comparison of degrees of freedom and CPU when taking different grids.
Grid numberDegrees of freedomCPU (s)516123255146432308220644522114585612310612825617294726According to the above different grid divisions, the change of pressurep/p0 is shown in Figure 18.Figure 18
Contour map of Mach numberM (a) and pressure p/p0 (b). When the grid is 516.
(a)(b)As can be seen from Figure19, with the increase of the number of grids, the initial position of the shock wave moves backward, but changes little, which shows the stability of the numerical method.Figure 19
Values of pressurep/p0 obtained with different grids.When the grid numbers are 516, 2064, and 5856, the obtained values of Mach numberM and pressure p/p0 are shown in Figures 18, 20, and 21.Figure 20
Contour map of Mach numberM (a) and pressure p/p0 (b). When the grid is 2064.
(a)(b)Figure 21
Contour map of Mach numberM (a) and pressure p/p0 (b). When the grid number is 5856.
(a)(b)The comparison of the values ofp/p0 on the right also verifies the conclusion in Figure 19. With the increase of the number of grids, the compression wave moves to the right, while the comparison of the values of M on the left shows that the mesh encryption makes the solution more and more smooth.Example 6.
For the steady problem of wet air, add adaptive grid technology and take Mach numberM=4, Φ0=30%. The results are shown in the following figure.
Figures22–24 are the contour maps of Mach number M, pressure p, and density ρ obtained without adaptive grid technology and with adaptive grid technology, respectively. Figure 25 shows the grid without adaptive grid technology and the grid with adaptive grid technology. It can be seen from the comparison diagram of the above results that the results obtained by using the original grid to solve the wet air nozzle flow are not smooth, but the grid obtained by modifying the original grid with adaptive grid technology. The obtained solution is smooth, and the existence of shock wave is more obvious. This shows that the mesh adaptive technology is an effective numerical method in solving unsteady problems. It not only makes the solution smooth, but also can encrypt the local mesh in the region where the shock wave exists, which is more helpful to capture the shock wave.Figure 22
Contour map of Mach numberM obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 23
Contour map of densityρ obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 24
Contour map of pressurep obtained without adding grid adaptive technology (a) and adding grid adaptive technology (b).
(a)(b)Figure 25
Original grid (a) and grid of grid adaptive technology (b).
(a)(b)Example 7.
Take the result of dry air pressure contraction wave as the initial condition of wet air nozzle flow, saturationΦ0=35.6%, the boundary condition is supersonic boundary condition, and the Mach numbers M=1.4 and M=2.5 of the incoming flow are taken. The numerical results are shown in the figure below.
Figures26–28 show the contour maps of Mach number M, pressure p, and density ρ with Mach numbers M=1.4 and M=2.5. It can be seen from the results that the Mach number of the whole flow field changes greatly with the increase of Mach number, while the contour maps of condensation mass fraction and condensation rate change little in the whole flow field. The condensation of wet air occurs in the throat of the nozzle. The linear comparison diagram of g/gmax and log10 J at different Mach numbers is obtained from the above contour map.
Figure29 shows g/gmax linear comparison diagram when the Mach number is 1.4 and 2.5, respectively. It can be seen that when the Mach number is 1.4, condensation occurs at the throat of the nozzle, while when the Mach number is 2.5, condensation exists in the whole flow field, but it is still obvious at the throat. When the fluid reaches the throat, condensation occurs. Next, due to the shock force at the throat, some droplets evaporate, making the condensation mass fraction continuously decrease. The subsequent expansion increases g/gmax.
Figure30 shows log10J linear comparison diagram when the Mach number is 1.4 and 2.5, respectively, indicating that the condensation rate increases with the increase of Mach number. The change of condensation rate and condensation mass fraction is synchronous. When the condensation rate reaches the maximum, the condensation mass fraction also reaches the maximum. The condensation rate is reduced to zero due to the impact of shock wave.Figure 26
Contour map of Mach numberM for M=1.4 (a) and M=2.5 (b).
(a)(b)Figure 27
Contour map ofg/gmax for M=1.4 (a) and M=2.5 (b).
(a)(b)Figure 28
Contour map of log10J for M=1.4 (a) and M=2.5 (b).
(a)(b)Figure 29
Linear comparison diagram ofg/gmax with different Mach numbers.Figure 30
Linear comparison diagram of log10J with different Mach numbers.


## 4.3. Conclusions

In this chapter, the stability and feasibility of the numerical method proposed in Chapter 3 are verified by the numerical simulation experiments of transonic flow and supersonic flow around the wing and parallel tail nozzle. The experimental results show that the stabilization method proposed in Chapter 3 has good stability in solving practical problems.The mesh adaptive technology is added to the steady problem, which better modifies the numerical solution of the original method and has a good ability to capture the shock wave.The numerical simulation of wet air flow around the wing shows that condensation occurs in the supersonic region and the shock wave will affect the condensation.The numerical simulation of wet air parallel wake nozzle flow shows that condensation occurs in the throat of the nozzle, which has certain application value.


## 5. Conclusions and Suggestions

### 5.1. Conclusions

In this paper, a set of partial differential equations and some algebraic relations are used to describe the condensation of compressible gas. The components are the liquid correlation system described by Euler equation and Hill moment method. The most important assumption in the solution process is that the fluid is inviscid flow and there is no slip between gas and liquid. The work done in this paper is as follows.When solving the coupled equations, most of the previous articles used the finite volume method to solve the transonic condensation flow problem, while this paper used the finite element method to solve this practical problem.Because the coupled Euler equations are convection equations, the results obtained by the classical Galerkin method often have nonphysical numerical oscillations, so the stabilization method is introduced. Another important feature of the compressible Euler equation is the existence of shock waves. Because the thickness of shock waves is small and difficult to capture, artificial viscosity terms are added to better capture shock waves. The streamline upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability are used to solve the transonic condensation flow problem.In dealing with the steady Euler equations, this paper also adds the adaptive mesh technology to mesh the region where the shock surface exists.The transonic flow around NACA-0012 airfoil and the flow field of parallel nozzle are numerically simulated by the method proposed in this paper.The conclusions of this paper are as follows:Hill moment method is derived from general dynamic equation.For single gas or mixed gas, the growth rate of liquid droplets follows the Hertz–Knudsen growth model and satisfies the state equation of ideal gas. The 2-D coupled equations with 8 variables are obtained by analysis and simplification.Through the numerical simulation of transonic flow around NACA-0012 airfoil and parallel nozzle, good numerical results are obtained, which verifies the stability and feasibility of this method.The mesh adaptive technology is added to the steady problem, which better modifies the numerical solution of the original method and has a good ability to capture the shock wave.The numerical simulation of wet air flow around the wing shows that condensation occurs in the supersonic region and the shock wave will affect the condensation.The numerical simulation of wet air parallel wake nozzle flow shows that condensation occurs in the throat of the nozzle.


### 5.2. Suggestions

Based on the research of this paper, we have the following plans for further research work:The numerical method is extended to the updated geometric model, such as ONERA M6 and F6 airfoils; the 2-D model can also be extended to 3-D model.When solving the problem, a new stabilization method can be added to the numerical method.The equation of state of gas can be extended from the ideal gas equation of state to the actual gas equation of state, so that the numerical results are closer to the actual results.The Euler equation can be replaced by N-S equation, and the inviscid flow model can be transformed into viscous flow model, which can be applied to a wider range of practical problems.


## 5.1. Conclusions

In this paper, a set of partial differential equations and some algebraic relations are used to describe the condensation of compressible gas. The components are the liquid correlation system described by Euler equation and Hill moment method. The most important assumption in the solution process is that the fluid is inviscid flow and there is no slip between gas and liquid. The work done in this paper is as follows.When solving the coupled equations, most of the previous articles used the finite volume method to solve the transonic condensation flow problem, while this paper used the finite element method to solve this practical problem.Because the coupled Euler equations are convection equations, the results obtained by the classical Galerkin method often have nonphysical numerical oscillations, so the stabilization method is introduced. Another important feature of the compressible Euler equation is the existence of shock waves. Because the thickness of shock waves is small and difficult to capture, artificial viscosity terms are added to better capture shock waves. The streamline upwind/Petrov–Galerkin (SUPG) finite element method with good stability and the isotropic diffusion method with shock capture ability are used to solve the transonic condensation flow problem.In dealing with the steady Euler equations, this paper also adds the adaptive mesh technology to mesh the region where the shock surface exists.The transonic flow around NACA-0012 airfoil and the flow field of parallel nozzle are numerically simulated by the method proposed in this paper.The conclusions of this paper are as follows:Hill moment method is derived from general dynamic equation.For single gas or mixed gas, the growth rate of liquid droplets follows the Hertz–Knudsen growth model and satisfies the state equation of ideal gas. The 2-D coupled equations with 8 variables are obtained by analysis and simplification.Through the numerical simulation of transonic flow around NACA-0012 airfoil and parallel nozzle, good numerical results are obtained, which verifies the stability and feasibility of this method.The mesh adaptive technology is added to the steady problem, which better modifies the numerical solution of the original method and has a good ability to capture the shock wave.The numerical simulation of wet air flow around the wing shows that condensation occurs in the supersonic region and the shock wave will affect the condensation.The numerical simulation of wet air parallel wake nozzle flow shows that condensation occurs in the throat of the nozzle.


## 5.2. Suggestions

Based on the research of this paper, we have the following plans for further research work:The numerical method is extended to the updated geometric model, such as ONERA M6 and F6 airfoils; the 2-D model can also be extended to 3-D model.When solving the problem, a new stabilization method can be added to the numerical method.The equation of state of gas can be extended from the ideal gas equation of state to the actual gas equation of state, so that the numerical results are closer to the actual results.The Euler equation can be replaced by N-S equation, and the inviscid flow model can be transformed into viscous flow model, which can be applied to a wider range of practical problems.



---

*Source: 9110254-2022-08-30.xml*