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

e nite element numerical simulation of the transonic ow with condensation is studied in this paper. e transformed weak form is numerically discretized by using the classical Galerkin nite element spatial discretization scheme and the backward Euler dierence 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 dicult to capture, the shock wave can be better captured by adding articial viscosity term. Here, the streamline upwind/ Petrov–Galerkin (SUPG) nite element method with good stability and the isotropic diusion method with shock capture ability are used to solve the problem of transonic ow 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 renement technology is added, which renes 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 ow around NACA-0012 and parallel-jet nozzle A1 has obtained the better numerical results.


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 signi cant component of oceans and lakes. Due to the longterm irradiation of the sun, the water on land evaporates. Meanwhile, water becomes water vapor around us. e 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. e 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 in uence 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. e 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. erefore, 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 highspeed 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.
is 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. e change range of cooling rate spans 10 − 3 -10 0 k/μs; due to condensation, the cooling rate is as high as − 10 2 k/μs. In this case, because of the latent heat released by condensation, we can observe a part of small droplet temperature. e 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. e 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. e condensation of wet air is the formation of clouds, fog, and dew in meteorology and climatology.

e Main Contributions and esis 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][10][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][13][14][15][16][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. e 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. e study of transonic flow has wide application value, because most of various types of aircraft are in the transonic range. e 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. e 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. e 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. erefore, 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. e 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 Section 3, 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. en, 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. erefore, 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 Section 4, 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. e feasibility of the numerical method is further verified, and the condensation of wet air in the nozzle is simulated. e 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.

Condensation Flow Model
Solid, liquid, and gas are the three most common phases. e 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. e thermodynamic phase diagram of pure substances is shown in Figure 1. In the diagram, it, respectively, represents the coexistence areas of solid-liquid, liquidgas, 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. ere 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.
Due 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. e degree of nonequilibrium state can be expressed by saturation Φ. It is defined as the ratio of the actual water vapor pressure P v to the water vapor pressure P v,eq (T) at equilibrium.
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. erefore, 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. e 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. e 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. e supersaturation rate in nucleation theory is given by equation (1). e change value of free energy when a molecular cluster is composed of molecules is given by Wisman [25], k B is a universal constant of micro-particle behavior, called the Boltz constant. k B � 1.38066 × 10 − 23 J/K, and the appearance of k B 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 v l represents the volume of liquid molecules (unit: m 3 /cluster). e 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, ΔG n is a monotonically increasing function with respect to n. erefore, 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. e critical value of n * is determined by (zΔG n /zn) n * � 0: θ is the dimensionless surface tension, defined as e critical free energy is obtained by substituting n * into equation (2).
en, use n * υ l � 4/3πr * 3 to get the Kelvin relation: In the formula, ρ l0 is condensation density and R v is gas constant of water vapor. e expression of critical radius r * 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.
e Kelvin relation has the following relations: v l � M/N A ρ l0 , M(kg/mol) is the molar mass of the gas, and N A � 6.02 × 10 26 /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 as J(Number of cores/m 3 · s), and the definition formula is ΔG n * is given by equation (6), and K is a constant. e condensation rate is given by the classical nucleation theory as follows: In the formula, ρ v is density of water vapor and m is the mass of a water vapor molecule. e growth rate of small droplets is expressed by the growth rate of small droplets. e 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, α is a free parameter, generally let α � 0.2, T d is the temperature of small droplets, p v is the pressure of water vapor, and p v,r is the supersaturated vapor pressure with a small droplet radius of r, defined as We usually assume that the temperature of small droplets is the same as that of the surrounding air, that is, e 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.
In the formula, U is vector consisting of conserved variables; F i is vector consisting of convective flux (i � 1, 2); and U, F i is shown below: ρ is density. p is pressure. u, v are the velocity components of the direction of x, y. e total energy e t is the sum of internal energy e and kinetic energy 1/2(u 2 + v 2 ), And the total enthalpy h t is defined as below: e Euler equation needs to be solved in combination with the equation of state. e 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 of V 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 M a . e mass of steam that can condense is recorded as M v . e mass of the condensed liquid is recorded as M l ; thus, volume V of total mass of mixed gas is M: When the total mass M 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 M v0 , e ratio of the mass of the condensed liquid to the total mass is called the mass fraction of the liquid and is recorded as g, e mass fraction of liquid reaches the maximum when all steam condenses, which is called the maximum mass fraction of liquid, i.e., At a certain volume, the density of the mixture is defined as Similarly, the density of inert gas and steam is e density of the condensed liquid ρ l (not condensation density ρ l0 ) similar can be expressed as From equations (23) and (24), the maximum mass fraction of the liquid can be obtained, which can also be expressed as And ρ v0 for the sum of ρ v and ρ l , e 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 Note that the precondition for the establishment of equations (27) and (28) is that the velocity of gas and condensed liquid is the same. erefore, we assume that there is no slip between gas phase and liquid phase. Bring equation (25) into equation (28), Expand equation (29), Mathematical Problems in Engineering e 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 It can also be written as us, it can be obtained that g max is a constant along a certain trajectory. From equation (32), we can get an important conclusion: if g max is a constant on the inlet boundary, g max is a constant on the whole flow field, so g max 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), en, combine formula (24) to obtain In the formula, b is energy provided by phase transition. e mass conservation equation of steam can be obtained by substituting equation (23) into equation (35). erefore, the mass conservation equations of mixed gas can be transformed into e 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 term b is not known. is 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. e phase space is ( x → , r) � (x, y, r), where ris the radius of small droplets. GDE can be referred to Seinfeld [27]: 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. n(r, x → , t) of small droplets to the radius r is In order to simplify the equation, the following independent variables f are omitted.
Both ends of equation (37) are multiplied by r k and integrated r from 0 to ∞, and From the properties of Dirac delta function, Introducing k − order distribution function, According to equations (41) and (42), equation (40) can be simplified as e integral can be simplified, and e 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 radius _ r from the integral, Hill [26] introduced the average small droplet radius: erefore, the radius growth rate _ r of small droplets no longer depends on the radius r H of droplets, but depends on the "Hill" radius r, that is, erefore, formula (43) is transformed into is equation is about μ k cycle. If k � 0, . . . , 3, partial differential equations can be obtained: 6 Mathematical Problems in Engineering Introduce equation Q k is the so-called liquid moment, and equation (48) can be transformed into e density ρ l of condensed liquid can be expressed as the following distribution function: We assume that when the material is constant, the condensation density ρ l0 is constant, which is obtained from formulas (42) and (51), From equations (24) and (49), erefore, 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" e average radius of Hill is obtained by replacing equation (45) with liquid moment Q k , In order to make Hill moment better describe the liquid mass fraction in PDE, we make the following assumptions: (1) ere is no slip between gas phase and liquid phase.
(2) Droplet growth rate of _ r is independent of droplet radius.
(3) e condensation density ρ l0 is constant. e mass fraction g of condensed liquid in Hill's method of moments can also be expressed by the following expression (Hill [26]): Compared with formula (53) in the previous section, this expression is based on physical properties instead of being expressed by the third moment Q 3 . is expression is tenable on the assumption that when t � 0.
ere 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 Mathematical Problems in Engineering Similar to formula (55), Hill radius r H 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: e equations can be gained, (59) e 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. erefore, equations (54)  e state equation is an expression composed of pressure p and temperature T, which is used to calculate flux and condensation terms. e static pressure of the gas mixture is the sum of the pressure p v of the steam and the pressure p a of the inert gas, i.e., Here, both inert gas and steam are regarded as ideal gases. Given that the two gas constants are R a and R v , respectively, the thermodynamic equation of state of gas can be obtained: Assuming thermodynamic equilibrium, the temperature of inert gas and steam is the same, i.e., (62) e equation of state of the gas mixture is R is the gas constant of the gas mixture, R 0 is the gas constant at which the mass of the condensed liquid is defined as e enthalpy hρ of the mixed gas is the sum of the enthalpies of the components, i.e., Combining equations (20), (22), and (23), the expression of specific enthalpy h (enthalpy per unit mass) can be obtained from equation (66): e heat of vaporization L is defined as the difference between the specific enthalpy h v of steam and the specific enthalpy h l of condensed liquid: Substitute equation (68) into equation (67) to obtain Similar to equation (67), the expression of specific internal energy e can be obtained: According to the definition of formula h � e + p/ρ and formula (68), we can get Assuming that the partial pressure of the condensed liquid is equal to the pressure of the mixed gas, i.e., p � p l , and then using the equation of state of the ideal gas, equation (71) can be transformed into Complete the gas according to the calorimetry, and obtain c v represents the specific heat under a certain volume, which is defined as

Mathematical Problems in Engineering
Assuming that both inert gas and steam are ideal gases, formula (72) can be converted into c v 0 represents the specific heat of the liquid at a certain volume when the mass fraction of the liquid is, i.e., If the specific internal energy e, 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 ρ/ρ l0 R at the right end of equation (75) can be ignored, i.e., 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): 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 U represents a vector composed of conserved variables: F i (i � 1, 2) represents the composition of convective flux: B represents a vector consisting of condensation source terms: 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 velocity a f represents the velocity of a fixed mixture in the minimal disturbance propagation medium, which is defined as s is specific entropy. In general, the equation of state is given as When the condensation mass fraction g is constant, the growth rate of pressure p is as follows: e first law of thermodynamics is Substitute equation (86) into equation (85) Mathematical Problems in Engineering dp � zp zρ e,g + p We define p as a function of ρ, s, and g. When g is constant, the growth rate of p is dp � zp zρ s,g dρ + zp zs ρ,g ds.
Compare the increments for dρ of equations (87) and (88), and get Equation (89) is the definition of frozen sound velocity a f . Substitute equations (63) and (78) into equation (89) to obtain For mixed gas, when the mass fraction of condensed liquid is 0, the temperature is T 0 , the pressure is p 0 , the mass fraction of maximum condensed liquid is g max , and the saturation is Φ 0 , it can be expressed as follows: Because saturation Φ 0 is known, g max can be obtained In the expressions of condensation rate J, 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 threephase 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.
erefore, we assume that the fluid only condenses into liquid. e surface tension σ(T)(N/m) of water adopts the expression in Dohrmann et al. [29], as follows: where T is Kelvin temperature, critical temperature T c � 647.3 K of water, and three-phase point temperature T tr � 273.15 K. e density ρ l0 (kg/m 3 ) of water is a function of temperature. e expression is given by Pruppacher and Klett [30]. e form at Kelvin temperature is ρ 1 � 999.84 kg/m 3 is the density of water at the three-phase point temperature, and the constants in formula (94) are given in Table 2. e pressure p v,eq (N/m 2 ) of steam in equilibrium is given by Sonntag and Heinze [31]: p v,eq (T) � e a 6 +a 7 T+a 8 T 2 +b 2 ln(T)+c 3 e expression of latent heat is as follows: e derivation of equation (95) shows that the latent heat of water is When the temperature is between 200 − 300 K, the latent heat of equation (98) can be approximately regarded as a linear function of T: and L 0 � 3105913.39(J/kg), and L 1 � − 2212.97(J/kg · K).  e gas constants R of dry air and specific heat capacity c p and c v , and specific heat c are shown in Table 3.

Numerical Methods
In Section 2, 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.
e boundary conditions and initial conditions of the equation are given. e 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. e 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. e 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.

Numerical Model.
e final form of the governing equation is given in Section 2. Equation (79) is the conservative form of the governing equation: Combined with the equation of state of gas proposed in Section 2, each variable is as follows: Refer to Section 2 for the meaning of each physical quantity.
Next, equation (100) is transformed into another form from the chain relationship of derivative: i.e., zW zt W, G is a vector, and A 1 � zF 1 /zU, A 2 � zF 2 /zU is a matrix, which can be obtained by combining the relationship between internal energy e t , enthalpy h t , and pressure p introduced in Section 2, Here, Here A � (A 1 , A 2 ), called Jacobian matrix. e 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 � ����� c 0 p/ρ, the governing equation of one-dimensional fluid is as follows: G i (i � 5, 6, 7, 8) given by (107). e characteristic polynomial obtained the Jacobian matrix of (7 × 7) from equation (109): Seven real eigenvalues are obtained

5, 6, 7). (111)
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. e number of physical boundaries is determined by whether the fluid is supersonic or subsonic [32].
For the inflow boundary conditions of one-dimensional subsonic flow (0 < 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. erefore, 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 (ρQ 0 , ρQ 1 , ρQ 2 , ρg), total pressure p 0 , and total temperature T 0 . After giving the p 0 and T 0 , the density ρ 0 can be obtained from the gas equation of state.
For supersonic flow (u > 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. erefore, we must give the values of all fluid variables at the boundary.
For the exit boundary condition of one-dimensional subsonic flow (0 < 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. erefore, only one boundary condition is given at the outlet boundary conditions, and the pressure p out is given in the calculation in this paper.
Similarly, for supersonic flow (u > 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. erefore, 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. erefore, 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.
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. e 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 of t � 0, W| t�0 � W 0 . 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.

Numerical Methods.
e 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. e grid is encrypted by grid adaptive technology.
According to the discussion earlier, the mathematical model for solving the problem is obtained: e variables in equation (113) refer to equation (106), Ω is the calculation area, and u e classical Galerkin finite element method is used to discretize equation (113) in space, and the test function space is introduced Take V the weight function space of ρ, u, v, p, Q 0 , Q 1 , Q 2 , g, respectively, multiply ∀U ∈ V 8 at both ends of equation (113), and integrate it on the Ω, Mathematical Problems in Engineering e second term of the above formula is obtained by using Green's formula and substituting the boundary conditions.
Next, the space is discretized to obtain the semi-discrete scheme of the problem. e bounded region Ω is approximated by a finite polygon Ω h , and Ω h is divided into several triangular elements Ω e , including Ω � ∪ where N e 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, e space V h is discretized with the weight function of ρ, u, v, p, Q 0 , Q 1 , Q 2 , g, respectively, and the Galerkin finite element space discretization of equation (117) is obtained: (119) e vector U uses the shape function N h of the standard finite element.
In (119), only the space is discretized to obtain the semidiscrete scheme, and then the time is discretized to obtain the fully discrete scheme. e time discretization of (119) adopts the backward Euler scheme. Firstly, the notation Q n � Ω × I n is introduced, where I n � [t n , t n+1 ], n � 0, 1, . . . n st − 1.
us, the unit element is transformed into Ω n e � Ω e × I n , e � 1, 2, . . . , N e , and the test function space is defined as follows: e problem to be solved by backward Euler scheme is as follows: for ∀n � 0, 1, . . . ,  where

Stabilization Techniques.
Shock wave is a strong compression wave in the moving gas. e 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][35][36][37][38][39][40][41][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. e 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, where L is differential operator. e stability term is added based on the standard Galerkin finite element method: where a(W, v) is bilinear form and R(W) � L(W) − f represents the residual.
When P(v) 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 � z/zt + (A · ∇) can be obtained. If P(v) � (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. is method is called SUPG (streamline upwind/ Petrov-Galerkin) method, which is based on the weak form of differential equation. e variational form of Galerkin method is obtained from equation (119): e SUPG method is added to the above formula to obtain the variational form after adding the stabilization term: the stable term is at 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. e 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 A i (i � 1, 2) will not be diagonalized at the same time.
In this paper, we take the stabilization coefficient as 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, where a 1 , a 2 is the spectral radius of matrix A 1 , A 2 (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: e 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 time 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 where a is the element order and a 1 , a 2 is the spectral radius of matrix A 1 , A 2 (the maximum of the absolute value of eigenvalues).

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. e 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.

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.
e 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.

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.

NACA-0012
Airfoil. For NACA-0012 airfoil, its ordinate can be expressed as 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. e following will be a numerical simulation of the air flow around NACA-0012 wet air.

Example 2.
e airfoil is NACA-0012, and the free flow condition is given as p ∞ � 65600 Pa, T ∞ � 259 K, and Φ 0 � 50%. e numerical results obtained when m is taken as 0.8 and 1.4, respectively, are given below. Example 3. the airfoil is NACA-0012, and the free flow condition is given as M ∞ � 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. Figures 10-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. Example 4. Moist air flows through NACA-0012 airfoil, and the free flow condition is M ∞ � 0.8, p ∞ � 65600 Pa, and T ∞ � 259 K. When the relative humidity is Φ 0 � 0%, Φ 0 � 50%, and Φ 0 � 60%, the numerical results are as follows. Figure 13 shows a comparison diagram of differential pressure resistance c p 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. e relative humidity is 50% mainly discussed, and the following results are obtained. Figure 14 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. e shock wave is generally generated in the supersonic region, so the condensation of wet air around the airfoil also occurs in the supersonic region.

Parallel-Jet Nozzle A1.
ere is little difference between the downstream part of parallel tail nozzle A1 and the throat part, and the nozzle has a straight outlet part. erefore, when the fluid passes through the nozzle, the fluid is very

24
Mathematical Problems in Engineering sensitive to the change of pressure and the composition of incoming wet steam. e geometric model of the nozzle is given in Figure 15. e nozzle has a straight outflow part with a length of about 104.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. e height of the outflow part of the nozzle is 2y * ≈ 92.69mm, so the Mach number of isentropic outflow is M ex � 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. e experimental device is shown in Figure 16. 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.
e nozzle flow of dry air and wet air will be numerically simulated below. e same boundary conditions are given for dry air and wet air, and the incoming boundary is T 0 � 298.7K and p 0 � 1.004bar. For the initial saturation (relative humidity) Φ 0 � 35.6% of wet air, the maximum mass fraction can be obtained g max � 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: e forced boundary condition p � 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 Table 4.
According to the above different grid divisions, the change of pressure p/p 0 is shown in Figure 18.
As can be seen from Figure 19, 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.
When the grid numbers are 516, 2064, and 5856, the obtained values of Mach number M and pressure p/p 0 are shown in Figures 18, 20, and 21. e comparison of the values of p/p 0 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.  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. e obtained solution is smooth, and the existence of shock wave is more obvious. is 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. 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. e condensation of wet air occurs in the throat of the nozzle. e linear comparison diagram of g/g max and log 10 J at different Mach numbers is obtained from the above contour map. Figure 29 shows g/g max 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. e subsequent expansion increases g/g max . Figure 30 shows log 10 J 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. e 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. e condensation rate is reduced to zero due to the impact of shock wave.

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.
e experimental results show that the stabilization method proposed in Chapter 3 has good stability in solving practical problems. e 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. e 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. e numerical simulation of wet air parallel wake nozzle flow shows that condensation occurs in the throat of the nozzle, which has certain application value.

Conclusions.
In this paper, a set of partial differential equations and some algebraic relations are used to describe the condensation of compressible gas. e components are the liquid correlation system described by Euler equation and Hill moment method. e most important assumption in the solution process is that the fluid is inviscid flow and there is no slip between gas and liquid. e 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. e 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. e transonic flow around NACA-0012 airfoil and the flow field of parallel nozzle are numerically simulated by the method proposed in this paper. e 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. e 2-D coupled equations with 8 variables are obtained by analysis and simplification.
rough 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. e 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. e 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. e numerical simulation of wet air parallel wake nozzle flow shows that condensation occurs in the throat of the nozzle.

Suggestions.
Based on the research of this paper, we have the following plans for further research work: e 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.
e 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.
e 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.

Mathematical Problems in Engineering
Data Availability e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.