chapter. Governing Equations 第2章 控制方程[cfd-0003]

If we restrict our attention to single-phase fluids, the law of mass conservation expresses the fact that mass cannot be created in such a fluid system, nor it can disappear. There is also no diffusive flux contribution to the continuity equation, since for a fluid at rest, any variation of mass would imply a displacement of the fluid particles.

In order to derive the continuity equation, consider the model of a finite control volume fixed in space, as sketched in Fig. 2.1. At a point on the control surface, the flow velocity is \(\vec{v}\), the unit normal vector is \(\vec{n}\), and \(dS\) denotes an elemental surface area. The conserved quantity in this case is the density \(\rho\). For the time rate of change of the total mass inside the finite volume \(\Omega\) we have

The mass flow of a fluid through some surface fixed in space equals to the product of (density) × (surface area) × (velocity component perpendicular to the surface). Therefore, the contribution from the convective flux across each surface element \(dS\) becomes

Since by convection \(\vec{n}\) always points out of the control volume, we speak of inflow if the product \(\left(\vec{v}\cdot\vec{n}\right)\) is negative, and of outflow if it is positive and hence the mass leaves the control volume.

As stated above, there are no volume or surface sources present. Thus, by taking into account the general formulation of Eq. (2.1), we can write

This represents the integral form of the continuity equation - the conservation law of mass.

We may start the derivation of the momentum equation by recalling the particular form of Newton's second law which states that the variation of momentum is caused by the net force acting on an mass element. For the momentum of an infinitesimally small portion of the control volume \(\Omega\) (see Fig. 2.1) we have

The variation in time of momentum within the control volume equals

Hence, the conserved quantity is here the product of the density and the velocity, i.e.,

The convective flux tensor, which describes the transfer of momentum across the boundary of the control volume, consists in the Cartesian coordinate system of the following three components

The contribution of the convective flux tensor to the conservation of momentum is then given by

The diffusive flux is zero since there is no diffusion of momentum possible for a fluid at rest. Thus, the remaining question is now, what are the forces the fluid element is exposed to? We can identify two kinds of forces acting on the control volume:

• External volume or body forces, which act directly on the mass of the volume. These are for example gravitational, buoyancy, Coriolis or centrifugal forces. In some cases, there can be electromagnetic forces present as well.
• Surface forces, which act directly on the surface of the control volume. They result from only two sources:
  • the pressure distribution, imposed by the outside fluid surrounding the volume,
  • the shear and normal stresses, resulting from the friction between the fluid and the surface of the volume.

From the above, we can see that the body force per unit volume, further denoted as \(\rho\vec{f}_e\), corresponds to the volume sources in Eq. (2.2). Thus, the contribution of the body (external) force to the momentum conservation is

The surface sources consist then of two parts - of an isotropic pressure component and of a viscous stress tensor \(\overline{\overline{\tau}}\), i.e.,

with \(\overline{\overline{I}}\) being the unit tensor (for tensors see, e.g., [2]). The effect of the surface sources on the control volume is sketched in Fig. 2.2. In Section 2.3, we shall elaborate the form of the stress tensor in more detail, and in particular show how the normal and the shear stresses are connected to the flow velocity.

Hence, if we now sum up all the above contributions according to the general conservation law (Eq. (2.2)), we finally obtain the expression

for the momentum conservation inside an arbitrary control volume \(\Omega\) which is fixed in space.

The underlying principle that we will apply in the derivation of the energy equation, is the first law of thermodynamics. Applied to the control volume displayed in Fig. 2.1, it states that any changes in time of the total energy inside the volume are caused by the rate of work of forces acting on the volume and by the net heat flux into it. The total energy per unit mass \(E\) of a fluid is obtained by adding its internal energy per unit mass, \(e\), to its kinetic energy per unit mass \(|\vec{v}|^{2}/2\). Thus, we can write for the total energy

The conserved quantity is in this case the total energy per unit volume, i.e., \(\rho E\). Its variation in time within the volume \(\Omega\) can be expressed as

Following the discussion in course of the derivation of the general conservation law (Eq. (2.1)), we can readily specify the contribution of the convective flux as

In contrast to the continuity and the momentum equation, there is now a diffusive flux. As we already stated, it is proportional to the gradient of the conserved quantity per unit mass (Fick's law). Since the diffusive flux \(\vec{F}_D\) is defined for a fluid at rest, only the internal energy becomes effective and we obtain

In the above, \(\gamma = c_p/c_v\) is the ratio of specific heat coefficients, and \(\kappa\) denotes the thermal diffusivity coefficient. The diffusion flux represents one part of the heat flux into the control volume, namely the diffusion of heat due to molecular thermal conduction - heat transfer due to temperature gradients. Therefore, Equation (2.7) is in general written in the form of Fourier's law of heat conduction, i.e.,

with \(k\) standing for the thermal conductivity coefficient and \(T\) for the absolute static temperature.

The other part of the net heat flux into the finite control volume consists of volumetric heating due to the absorption or emission of radiation, or due to chemical reactions. We will denote the heat sources - the time rate of heat transfer per unit mass - as \(\dot{q}_h\). Together with the rate of work done by the body forces \(\vec{f}_e\), which we have introduced for the momentum equation, it completes the volume sources

The last contribution to the conservation of energy, which we have yet to determine, are the surface sources \(Q_S\). They correspond to the time rate of work done by the pressure as well as the shear and normal stresses on the fluid element (see Fig. 2.2), i.e.,

Sorting now all the above contributions and terms, we obtain for the energy conservation equation the expression

The energy equation (2.11) is usually written in a slightly different form. For that purpose, we will utilise the following general relation between the total enthalpy, the total energy and the pressure

When we now gather the convective (\(\rho E\vec{v}\)) and the pressure term (\(p\vec{v}\)) in the energy conservation law (2.11), and apply the formula (2.12), we can finally write the energy equation in the form

Herewith, we have derived integral formulations of the three conservation laws: the conservation of mass (2.3), of momentum (2.5), and of energy (2.13). In the next section, we shall work out the formulation of the normal and the shear stresses in more detail.

In pure aerodynamics, it is generally reasonable to assume that the working fluid behaves like a calorically perfect gas, for which the equation of state assumes the form [13], [14]

where \(R\) denotes the specific gas constant. The enthalpy results from

It is convenient to express the pressure in terms of the conservative variables. For that purpose, we have to combine Equation (2.12), relating the total enthalpy to the total energy, together with the equation of state (2.26). Substituting expression (2.27) for the enthalpy and using the definitions

we finally obtain for the pressure

The temperature is then calculated with the aid of the relationship Eq. (2.26). The coefficient of the dynamic viscosity \(\mu\) is, for a perfect gas, strongly dependent on temperature but only weakly dependent on pressure. The so-called Sutherland formula is frequently used. The result for air is (in SI units)

where the temperature \(T\) is in degree Kelvin (K). Thus, at \(T = 288\,\mathrm{K}\) one obtains \(\mu = 1.78\cdot 10^{-5}\,\mathrm{kg/ms}\). The temperature dependence of the thermal conductivity coefficient \(k\) resembles that of \(\mu\) in the case of gases. By contrast, \(k\) is virtually constant in the case of liquids. For this reason, the relationship

is generally used for air. In addition, it is commonly assumed that the Prandtl number \(Pr\) is constant in the entire flow field. For air, the Prandtl number takes the value \(Pr = 0.72\).

The matter becomes more complicated when one has to deal with a real gas. The reason is that now we have to model a thermodynamic process and chemical reactions in addition to the fluid dynamics. Examples for a real gas flow are the simulation of combustion, the hypersonic flow past a re-entry vehicle, or the flow in a steam turbine. In principle, two different methods can be pursued to solve the problem. The first methodology is applicable in cases, where the gas is in chemical and in thermodynamical equilibrium. This implies that there is a unique equation of state. Then, the governing equations (2.19) remain unchanged. Only the values of pressure, temperature, viscosity, etc. are interpolated from lookup tables using curve fits [15]-[17]. But in practice, the gas is more often in chemical and/or thermodynamical non-equilibrium and has to be treated correspondingly.

Let us for illustration consider a gas mixture consisting of \(N\) different species. For a finite Damköhler number, which is defined as the ratio of flow-residence time to chemical-reaction time, we have to include finite-rate chemistry into our model. It has to describe the generation/destruction of species due to chemical reactions. In what follows, we will furthermore assume that the temporal and the spatial scales of fluid dynamics and of chemical reactions are much larger compared to those of thermodynamics. Thus, we suppose the gas is thermodynamically in equilibrium but chemically in non-equilibrium. In order to simulate the behaviour of such a gas mixture, the Navier-Stokes equations have to be augmented by \((N-1)\) additional transport equations for the \(N\) species [18]-[23]. Hence, we obtain formally the same system like Eq. (2.19), but now with the vectors of the conservative variables \(\vec{W}\), the flux vectors \(\vec{F}_c\) and \(\vec{F}_v\), as well as with the source term \(\vec{Q}\) extended by \((N-1)\) species equations. Recalling the expressions (2.20) to (2.25), the vector of the conservative variables reads now

The convective and the viscous flux vectors transform into

where

Finally, the source term now becomes

In the above expressions Eqs. (2.32)-(2.35), \(Y_m\) denotes the mass fraction, \(h_m\) the enthalpy, and \(D_m\) the effective binary diffusivity of species \(m\), respectively. Furthermore, \(\dot{s}_m\) is the rate of change of species \(m\) due to chemical reactions. Note that the total density \(\rho\) of the mixture is equal to the sum of the densities of the species \(\rho Y_m\). Therefore, since the total density is regarded as an independent quantity, there are only \((N-1)\) independent densities \(\rho Y_m\) left. The remaining mass fraction \(Y_N\) is obtained from

In order to find an expression for the pressure \(p\), we first assume that the individual species behave like ideal gases, i.e.,

with \(R_u\) denoting the universal gas constant and \(W_m\) being the molecular weight, respectively. Together with Dalton's law,

we can write

It is important to notice that because the gas is in thermodynamical equilibrium, all species possess the same temperature \(T\). The temperature has to be calculated iteratively from the expression [21], [24]

The internal energy of the gas mixture \(e\) is obtained from Eq. (2.6). The quantities \(h^{0}_{f,m}\), \(c_{p,m}\), and \(T_{ref}\) in Eq. (2.40) denote the heat of formation, the specific heat at constant pressure, and the reference temperature of the \(m\)-th species, respectively. Values of the above quantities as well as of the thermal conductivity \(k\) and of the dynamic viscosity \(\mu\) of the species are determined from curve fits [19], [21], [23].

The last part, which remains to be modelled, is the chemical source term \(\dot{s}_m\) in Eq. (2.35). The rate equations for a set of \(N_R\) elementary reactions involving \(N\) species can be written in the general form

In the above Eq. (2.41), \(\nu'_{lm}\) and \(\nu''_{lm}\) are the stoichiometric coefficients for species \(m\) in the \(l\)-th forward and backward reaction, respectively. Furthermore, \(C_m\) stands for the molar concentration of species \(m\) (\(C_m = \rho Y_m/W_m\)), and finally \(K_{fl}\) and \(K_{bl}\), respectively, denote the forward and the backward reaction rate constants for the \(l\)-th reaction step. They are given by the empirical Arrhenius formulae

where \(A_f\) and \(A_b\) are the Arrhenius coefficients, \(E_f\) and \(E_b\) represent the activation energies, and \(B_f\) as well as \(B_b\) are constants, respectively. The rate of change of molar concentration of species \(m\) by the \(l\)-th reaction is given by

Hence, together with Eq. (2.43) we can calculate the total rate of change of species \(m\) from

More details can be found in the references cited above. A detailed overview of the equations governing a chemically reacting flow, together with the Jacobian matrices of the fluxes and their eigenvalues, can also be found in [24].

Another practical example of real gas is the simulation of steam or, which is more demanding, of wet steam in turbomachinery applications [25]-[32]. In the later case, where the steam is mixed with water droplets, so that we speak of multiphase flow, it is either possible to solve an additional set of transport equations, or to trace the water droplets along a number of streamlines. These simulations have very important applications in the design of modern steam turbine cascades. The analysis of flow past turbine blades can for instance help to understand the occurrence of supercritical shocks by condensation and of flow instabilities, responsible for an additional dynamic load on the bladings and resulting in a loss of the efficiency.

In the following, we shall consider three common simplifications to the Navier-Stokes equations (2.19). We shall restrict our attention here to the physical reasoning behind each of the approximations. The equations for the first two simplified forms of the Navier-Stokes equations are provided in the Appendix.

When simulating flows around bodies for high Reynolds numbers (i.e., when the boundary layer is thin with respect to a characteristic dimension), the Navier-Stokes equations (2.19) can be simplified. One necessary condition is that there is no large area of separated boundary layer. It can then be anticipated that

only the gradients of the flow quantities in the normal direction to the surface of the body (\(\eta\)-direction in Fig. 2.4) contribute to the viscous stresses [33], [34]. On the other hand, the gradients in the other coordinate directions (\(\xi\) in Fig. 2.4) are neglected in the evaluation of the shear stress tensor (Eqs. (2.14, 2.15)). We speak here of the so-called Thin Shear Layer (TSL) approximation of the Navier-Stokes equations. The motivation for the TSL modification is that the numerical evaluation of the viscous terms becomes computationally less expensive, but, within the assumptions, the solution remains sufficiently accurate. The TSL approximation can also be justified from a practical point of view. In the case of high Reynolds number flows, the grid has to be very fine in the wall normal direction in order to resolve the boundary layer properly. Because of the limited computer memory and speed, much coarser grid has to be generated in the other directions. This in turn results in significantly lower numerical accuracy of the gradient evaluation compared to the normal direction. The TSL equations are for completeness presented in the Appendix (A.6). Due to the fact that secondary flow (e.g., like in a blade row) cannot be resolved appropriately, the TSL simplification is usually applied only in external aerodynamics.

In cases, where the following three conditions are fulfilled:

• the flow is steady (i.e. \(\partial\vec{W}/\partial t = 0\)),
• the fluid moves predominantly in one main direction (e.g., there must be no boundary layer separation),
• the cross-flow components are negligible,

the governing equations (2.19)) can be simplified to a form called the Parabolised Navier-Stokes (PNS) equations [8], [35]-[37]. The above conditions allow us to set the derivatives of \(u\), \(v\), and \(w\) with respect to the streamwise direction to zero in the viscous stress terms (Eq. (2.15)). Furthermore, the components of the viscous stress tensor \(\overline{\overline{\tau}}\), of the work performed by it (\(\overline{\overline{\tau}}\cdot\vec{v}\)), and of the heat conduction \(k\nabla T\) in the streamwise direction are dropped from the viscous flux vector in Eq. (2.23). The continuity equation, as well as the convective fluxes (Eq. (2.21)) remain unchanged. For details, the reader is referred to the Appendix (A.7). Considering the situation sketched in Fig. 2.5, where the main flow direction coincides with the \(x\) coordinate, it can be shown that the PNS approximation leads to a mixed set of parabolic / elliptic equations. Namely, the momentum equation in the flow direction becomes parabolic together with the energy equation, and hence they can be solved by marching in the \(x\)-direction. The momentum equations in the \(y\)- and in the \(z\)-direction are elliptic and they have to be solved iteratively in each \(x\)-plane. Thus, the main benefit of the PNS approach is in the largely reduced complexity of the flow solution - from a complete 3-D field to a sequence of 2-D problems. A typical application of the parabolised Navier-Stokes equations is the calculation of internal flows in ducts and in pipes, and also the simulation of steady supersonic flows using the space-marching method [38]-[41].

As we have seen, the Navier-Stokes equations describe the behaviour of a viscous fluid. In many instances, it is a valid approximation to neglect the viscous effects completely, like for example for high Reynolds-number flows, where the boundary layer is very thin compared to the dimensions of the body. In such cases, we can simply omit the vector of viscous fluxes, \(\vec{F}_v\), from the Equations (2.19). Thus, we are left with

The remaining terms are given by the same relations (2.20)-(2.22) and Eq. (2.25) as before. This simplified form of the governing equations is called the Euler equations. They describe the pure convection of flow quantities in an inviscid fluid. If the Euler equations are formulated in conservative way (like above), they allow for accurate representation of such important phenomena like shocks, expansion waves and vortices over delta wings (with sharp leading edges). Furthermore, the Euler equations served in the past - and still do - as the basis for the development of discretisation methods and boundary conditions.

However, it should be noted that today, due to the computational power of even personal computers and due to the increased demands on the quality of the simulations, the Euler equations are only relatively seldom employed for flow computations.