Numerical simulation of the passive gas mixture ﬂow

. The aim of this paper is the numerical solution of the equations describing the non-stationary compressible turbulent multicomponent ﬂow in gravitational ﬁeld. The mixture of perfect inert gases is assumed. We work with the RANS equations equipped with the k-omega and the EARSM turbulence models. For the simulation of the wall roughness we use the modiﬁcation of the speciﬁc turbulent dissipation. The ﬁnite volume method is used, with thermodynamic constants being functions in time and space. In order to compute the ﬂuxes through the boundary faces we use the modiﬁcation of the Riemann solver, which is the original result. We present the computational results, computed with the own-developed code (C, FORTRAN, multiprocessor, unstructured meshes in general).


Introduction
The aim of this work is to simulate the complicated behaviour of the perfect gas mixture.In this contribution we consider the Reynolds-Averaged Navier-Stokes equations with the k-ω model of turbulence.This system is equipped with the equation of state in more general form, and with the mass conservation of the additional gas specie.We choose the well-known finite volume method to discretize the analytical problem, represented by the system of the equations in generalized (integral) form.In order to apply this method we split the area of the interest into the elements, and we construct a piecewise constant solution in time.The crucial problem of this method lies in the evaluation of the so-called fluxes through the edges/faces of the particular elements.We use the analysis of the exact solution of the Riemann problem for the discretization of the fluxes through the boundary edges/faces.Own algorithms for the solution of the boundary problem were coded, and used in the numerical examples.

Formulation of the Equations
We consider the conservation laws for viscous compressible turbulent flow of ideal gas with the zero heat sources in a domain Ω ∈ I R N , and time interval (0, T ), with T > 0. The system of the Reynolds-Averaged Navier-Stokes equations in 3D has the form ∂I R s (w, ∇w) ∂x s +S in Q T = Ω×(0, T ).
(1) Here x 1 , x 2 , x 3 are the space coordinates, t the time, w = w(x, t) = ( , v 1 , v 2 , v 3 , E) T is the state vector, f s = a e-mail: kyncl@vzlu.czb e-mail: pelant@vzlu.cz( v s , v s v 1 + δ s1 p, v s v 2 + δ s2 p, v s v 3 + δ s3 p, (E + p) v s ) T are the inviscid fluxes, I R s = (0, τ s1 , τ s2 , τ s3 , 3 r=1 τ sr v r + C k ∂θ/∂x s ) T are the viscous fluxes, S are additional sources.u = (v 1 , v 2 , v 3 ) T denotes the velocity vector, is the density, p the pressure, θ the absolute temperature, E = e + 1 2 u 2 + k the total energy.Further τ i j = (μ + μ T )S i j , i j (μ + μ T )S i j − 2  3 k, i = j , with S 11 = 2 3 2 ∂v 1 ∂x 1 − ∂v 2 ∂x 2 − ∂v 3 ∂x 3 , S 12 = ∂v 1 ∂x 2 + ∂v 2 ∂x 1 , S 13 = ∂v 1 ∂x 3 + ∂v 3 ∂x 1 , S 21 = S 12 , S 22 = 2 3 − ∂v 1 ∂x 1 + 2 ∂v 2 ∂x 2 − ∂v 3 ∂x 3 , S 23 = ∂v 2 ∂x 3 + ∂v 3 ∂x 2 , S 31 = S 13 , S 32 = S 23 , S 33 = 2  3 − ∂v 1 ∂x 1 − ∂v 2 ∂x 2 + 2 ∂v 3 ∂x 3 , where μ is the dynamic viscosity coefficient dependent on temperature, μ T is the eddy-viscosity coefficient.For the specific internal energy e = c v θ we assume the caloric equation of state e = p/ (γ − 1), c v is the specific heat at constant volume, γ > 1 is called the Poisson adiabatic constant.The constant C k denotes the heat conduction coefficient C k = ( μ P r + μ T P r T )c v γ, and P r is laminar and P r T is turbulent Prandtl constant number.In our application of flow in the gravitational field we set the source terms to S = (0, g 1 , g 2 , g 3 , g • u), where g = (g 1 , g 2 , g 3 ) is the gravity vector.For the gas mixture with two species we use the Dalton's law for the total mixture pressure where p 1 and p 2 are the partial pressures of the first and second component gas.Let 1 and 2 be the mass density of these components.Then the total mass density of mixture is Temperature θ is same for all gases in the mixture, and the equation of state holds where R g = 8.3144621 is universal gas constant, and m i denotes the mollar mass of the ith specie.We can introduce The thermodynamic constants of the mixture satisfy (using the decomposition of the internal specific energy and enthalpy) then the adiabatic constant γ, needed in the solution of (1), can be written as The system (1) is then extended with the conservation law of the mass for one gas component (specie) Here σ C is diffusion coefficient.The mass conservation for the second specie is automatically satisfied via the system (1).
Here we assume the system (1),( 2) equipped with the two-equation turbulent model k−ω (Kok), described in [1].The effective turbulent viscosity is where k the specific turbulent kinetic energy and ω the turbulent dissipation are functions of time t and space coordinates x 1 , x 2 , x 3 .The production terms P k and P ω are given by formulas where functions τ are defined in (1) where σ d = 0.5 is constant.

Numerical Method
For the discretization of the system we proceed as described in [6].We use either explicit or implicit finite volume method to solve the systems sequentionally.Here we present the discretization of the system (1).By Ω h let us the denote the polyhedral approximation of Ω.The system of the closed polyhedrons with mutually disjoint interiors D h = {D i } i∈J , where J ⊂ Z + = {0, 1, . ..} is an index set and h > 0, will be called a finite volume mesh.This system D h approximates the domain Ω, we write Ω h = i∈J D i .The elements D i ∈ D h are called the finite volumes.For two neighboring elements D i , D j we set Γ i j = ∂D i ∩ ∂D j = Γ ji .Similarly, using the negative index j we may denote the boundary faces.Here we will work with the so-called regular meshes, i.e. the intersection of two arbitrary (different) elements is either empty or it consists of a common vertex or a common edge or a common face (in 3D).The boundary ∂D i of each element D i is Here the set Using the Green's theorem on D i it is Here n = (n 1 , n 2 , n 3 ) is the unit outer normal to ∂D i .Further we use (5), and we rewrite (6) We define a finite volume approximate solution of the system studied (1) as a piecewise constant vector-valued functions w k h , k = 0, 1, . . ., where w k h is constant on each element D i , and t k is the time instant.By w k i we denote the value of the approximate solution on D i at time t k .We approximate the integral over the element Further we proceed with the approximation of the fluxes.
Usually the flux with w l i , w l j denoting the approximate solution on the elements adjacent to the edge Γ i j at the time instant t l .In the case of a boundary face the vector w l j has to be specified.Here we show the numerical flux based on the solution of the Riemann problem for the split Euler equations.By w l Γ i j let us denote the state vector w at the center of the edge Γ i j at the time instant t l , and let us suppose w l Γ i j is known.Evaluation of these values will be a question of the further analysis, here we use them to approximate the integrals with the one-point rule (10) Here ∇w l Γ i j denotes the ∇w at the center of the edge Γ i j at time instant t l .Now it is possible to approximate the system (8) by the following explicit finite volume scheme With this finite volume formula one computes the values of the approximate solution at the time instant t k+1 , using the values from the time instant t k , and by evaluating the values w k Γ i j at the faces Γ i j .In order to achieve the stability of the used method, the time step τ k must be restricted by the so-called CFL condition, see [4].The crucial problem of this discretization lies with the evaluation of the edge values w k Γ i j .Or one deals with the problem of finding the face fluxes H(w k i , w k j , n i j ).It is also possible to use the implicit scheme The crucial problem of this discretization lies with the evaluation of the face fluxes H(w k+1 i , w k+1 j , n i j ).One possibility is to use the linearization via the Taylor expansion of the vector function H(wI, wJ, n), this was shown in [6].One possibility of the face flux evaluation is to approximate the face values w k Γ i j and then compute the numerical flux To approximate the face values w k Γ i j at time instant t k we solve the simplified system (14) in the vicinity of the face Γ i j in time with the initial condition formed by the state vectors w k i and w k j .Using the rotational invariance of the Euler equations, the system is expressed in a new Cartesian coordinate system x1 , x2 , x3 with the origin at the center of the gravity of Γ i j and with the new axis x1 in the direction of n = (n 1 , n 2 , n 3 ), given by the face normal n = n i j .
x 1 Then the derivatives with respect to x2 , x3 are neglected and we get the so-called split 3D Euler equations, see [4, page 138]: The values w k i and w k j adjacent to the face Γ i j are known, forming the initial conditions The transformation matrix Q is defined in [6].The problem ( 14), ( 15), ( 16) has a unique "solution" in (−∞, ∞) × (0, ∞), the analysis can be found in [4, page 138].Let q RS (q L , q R , x1 , t) denote the solution of this problem at the point ( x1 , t).We are interested in the solution of this local problem at the time axis, which is the sought solution in the local coordinates q Γ i j = q RS (q L , q R , 0, t).The backward transformation of the state vector q Γ i j into the global coordinates is The constructed numerical flux can be written as Let Γ i j be the face of the element D i laying at the boundary of the computational area.To approximate the face values w k Γ i j at time instant t k we solve the incomplete simplified system (19) in the vicinity of the face Γ i j in time with the initial condition (20) determined by the state vector w k i .By adding properly chosen equations into the system (19),(20) it is possible to reconstruct the boundary state q B such that the system (19),(20) has a unique solution at the boundary, see [3].We will refer to these added equations as to complementary conditions.Several choices of the complementary conditions were discussed in [3], [5].

The Riemann Problem for the Split Euler Equations
For many numerical methods dealing with the two or three dimensional equations, describing the compressible flow, it is useful to solve the Riemann problem for the 3D split Euler equations.We search the solution of the system of partial differential equations in time t and space (x, y, z) equipped with the initial conditions Vector u = (u, v, w) denotes the velocity, density, p pressure, E = e + |u| 2 is the total energy, with e denoting the specific internal energy.We assume the equation of state for the calorically ideal gas e = p (γ − 1) .
'Split' means here that we still have 5 equations in 3D, but the dependence on y, z (space coordinates) is neglected, and we deal with the system for one space variable x.The system (19) is considered in the set The solution of this problem is fundamentally the same as the solution of the Riemann problem for the 1D Euler equations, see [4, page 138].In fact, the solution for the pressure, the first component of the velocity, and the density is exactly the same as in one-dimensional case.It is a characteristic feature of the hyperbolic equations, that there is a possible raise of discontinuities in solutions, even in the case when the initial conditions are smooth, see [8, page 390].Here the concept of the classical solution is too restrictive, and therefore we seek a weak solution of this problem.To distinguish physically admissible solutions from nonphysical ones, entropy condition must be introduced, see [8, page 396].By the solution of the problem (19),(20),(21) we mean the weak entropy solution of this problem in Q ∞ .The analysis to the solution of this problem can be found in many books, i.e. [4], [8], [9].The general theorem on the solvability of the Riemann problem can be found in [4, page 88].
The solution is piecewise smooth and its general form can be seen in Fig. 2, where the system of half lines is drawn.These half lines define regions, where solution is either constant or given by a smooth function.Let us define the open sets called wedges (see Fig. 2): Using the theory in [4], [8], [9], we can write the solution for the primitive variables as The folowing relations for these variables hold: Here a L = γp L / L , a R = γp R / R , and γ denotes the adiabatic constant.Further s 1 T L denotes "unknown left wave speed", s 3  T R "unknown right wave speed".Note, that the system (22) -(27) is the system of 6 equations for 6 unknowns p ,u , L R ,s 1 T L ,s 3 T R .

Solution for the Pressure p
The combination of the equations ( 22), (25) gives the implicit equation for the unknown pressure p with This is a nonlinear algebraic equation, and one cannot express the analytical solution of this problem in a closed form.The solution p can be found as the root of the function F(p) defined as The analysis of this function can be found in [4].Here we state, that F(p) is monotone and concave.Also the analytic expression for its derivative is available.For a positive solution for the pressure F(0) < 0 is required.This gives the pressure positivity condition The Newton method can be used to find the root of F(p) = 0.With the pressure p known, we use the equations ( 22)-( 27) to obtain the whole solution.

Remarks
-Once the pressure p is known, the solution on the lefthand side of the contact discontinuity depends only on the left-hand side initial condition (20).And similarily, with p known, only the right-hand side initial condition ( 21) is used to compute the solution on the righthand side of the contact discontinuity.-The solution in Ω L ∪Ω HT L ∪Ω L (across 1 wave) (solvability in general case) The solution components in Ω L ∪ Ω HT L ∪ Ω L region must satisfy the system of equations ( 22)-(24).It is a system of three equations for four unknowns ( L , p , u , s 1 T L ).We have to add another equation in order to get the uniquely solvable system in Ω L ∪ Ω HT L ∪ Ω L .This is the key problem for the outlet boundary condition.
-The equation ( 22) can be written as the equation for pressure, see [3] p = E 1 (u ), (30) In our work we use the analysis of this problem also for the solution of the initial-boundary value at the boundary faces.

Boundary Condition by Preference of Pressure
Using this boundary condition, we try to prescribe given static pressure at the face.This boundary condition corresponds to the real-world problem, when we deal with the experimentally obtained pressure distribution at the boundary.The conservation laws must be satisfied in close vicnity to the boundary face.We use the analysis of the incomplete Riemann problem to construct the values for the density and the pressure at each boundary face.
Here p denotes the pressure in Ω L ∪Ω R part of the solution of the Riemann problem for the split Euler equations, shown in Section 4. This way we have prescribed the desired pressure value whenever it is possible (see the general form of the solution in Section 4).Now it is possible to use ( 22),( 23),(24) Further we must discuss all the possibilities of the left (shock, rarefaction) and the middle (discontinuity) wave in the solution, which is shown in figure 3.In the case of outlet (u ≥ 0), the transformation of the velocity into the global coordinates is Here u B denotes the computed velocity in the normal direction.In the case of an inlet (u < 0) we must prescribe another conditions in order to obtain a uniquely solvable system.Here we may choose (for example) arbitrary density GIVEN and the direction of the velocity d.
Then it is The analysis of this problem was shown also in [3]. 02064-p.5

Boundary Condition by Preference of Velocity
Here we show the boundary condition prefering the prescribed given velocity at the face.The conservation laws must be satisfied in close vicnity to the boundary face.
Here u denotes the normal component of velocity in Ω L ∪ Ω R part of the solution of the Riemann problem for the split Euler equations, shown in Section 4. This way we have prescribed the desired normal component of velocity whenever it is possible (see the general form of the solution in Section 4).Now it is possible to use (30),( 23),(24).Further we must discuss all the possibilities of the left (shock, rarefaction) and the middle (discontinuity) wave in the solution, which is shown in figure 4. The tangential velocities are conserved (same as "left" values) in the case of outlet (u ≥ 0), and the velocity vector in the global coordinates is Here u B denotes the computed velocity in the normal direction.In the case of an inlet (u < 0) we must prescribe another conditions in order to obtain a uniquely solvable system.Here we may choose (for example) arbitrary density GIVEN and the whole vector of velocity For example, it is possible to choose the density R using the given entropy index p o o γ (p o , o are some reference values for the pressure and the density) EFM 2015 Another possibility is to conserve (choose) the total temperature θ 0 for the inlet case.Then, using thermodynamic relations (energy equation for a perfect gas), it is The algorithm for the construction of the primitive variables B , u B , p B at the boundary face is shown in Fig. 4. The complete analysis of this problem is shown in [3].
u L ≤ 0 (1-rarefaction wave) The velocity in local coordinates is (0, v L , w L ).The backward transformation of the computed values back into global coordinates (transformation of the velocity components) gives the sought velocity vector at the boundary Remark.It is possible to solve this boundary problem by prescribing the special right-hand side initial condition to the original Riemann problem, described in Section 4., i.e.
Non-slip wall For the non-slip wall (viscous flow) we set the velocity vector u w = 0 at the boundary.The turbulent kinetic energy at such non-slip wall is k = 0.For the evaulation of the specific dissipation at the boundary ω we use the statements from [1].Let us use the following notation and relations θ w static temperature at the wall μ w dynamic viscosity at the wall, using Sutherland's law unit outer normal to face force (vector) acting on surface surface shear stress, friction tension (force magnitude in surface direction) friction velocity (velocity scale) In the real application it is necessary to prescribe the static temperature of the surface θ w .
The law of wall, [1, pages 13-18] The law of the wall is empirically-determined relation between the streamwise velocity U and distance from surface y, see [ The coefficient S R is defined with the use of the roughness height k S = z S e 8.5κ .
Here k S represents the roughness height.S R is a continuous function of k + S .

The Flat Plate Example
Here we show example of the compressible viscous flow along the flat plate.The velocity regime u = 15 m s −1 .The total pressure p o = 101325 Pa, total temperature θ o = 273.15K, and properties R = 287.04J K −1 mol −1 , γ = 1.4,gas flows from the left to the right side.Initial condition for the computation is given by isentropic approximation (and the use of energy equation for a perfect gas) The initial conditions for the free stream turbulent kinetic energy k and the specific turbulent dissipaton ω are given by the turbulent intensity I = 0.1 and μ T /μ = 0.01:

Numerical Tests on Shortened Domains
In order to achieve faster computations it is desired to use smaller computation areas.In such case, the boundary conditions must be chosen appropriately.Here we present numerical experiments with the boundary conditions on smaller domains.At first we assume the compressible viscous flow in the long channel (x-coordinate from -13 to 4), with only small area of interest (x-coordinate from -1 to 2).The velocity regime u = 15 m s −1 .The total pressure p o = 101325 Pa, total temperature θ o = 273.15K, and the gas constants set to R = 287.04J K −1 mol −1 , γ = 1.4,gas flows from the left to the right side.Initial condition for the computation is given by isentropic approximation (and the use of energy equation for a perfect gas).The initial conditions for the free stream turbulent kinetic energy k and the specific turbulent dissipaton ω are given by the turbulent intensity I = 0.1 and μ T /μ = 0.01, see also section 7.1.The results are depicted in figure 7.
The computed results were used in another series of computations on shortened domain (with y-coordinate -1 to 2), in order to test various boundary conditions.The aim was to test various boundary and initial conditions in order to match previously computed data on the original (longer) channel.Here we present the results using the boundary conditions given by the preference of pressure at the outlet, and with the preference of total quantities and the direction of velocity at the inlet.The input values for these boundary conditions were given from the previous computation on larger domain (different values at each boundary face).Results are shown in figures 8 -10.
Further we show the simulation of emission of another gas specie into the area.We choose the additional gas with  the properties R = 461 J K −1 mol −1 , c V = 1450 J K −1 mol −1 , which enters the area.This emission is realized by introducing the fixed source located at the given cell (with cell center coordinates [0.092,0.0997,0.5]).The mass fraction of additional gas component at this cell was fixed to Y 1 = 1, and density was set accordingly to equation of state, using the computed pressure and temperature.Figure 11.
shows the propagation of mixture in time.

Conclusion
This paper shows the formulation of the equations describing the mixture of two inert perfect gases in 3D, it is further focused on the boundary conditions and simulation of 02064-p.9   the wall roughness.The numerical method (finite volume method) is applied for the solution of these equations.Own software was programmed.The modification of the Riemann problem is used at the boundaries.

Fig. 3 .
Fig. 3. Algorithm for the solution of the problem (19),(20),(31),(32),(33) at the half line S B = {(0, t); t > 0}.The value of the pressure p is prescribed.Possible situations are illustrated by the pictures showing Ω L ∪ Ω HT L ∪ Ω L region.The sought boundary state located at the time axis, marked red.

Fig. 5 .
Fig. 5. Computational simulation, the flow along the flat plate and along the rough wall (z S = 0.0002) simulation, regime 15 m s −1 , isolines of k turbulent kinetic energy.The Y + − U + graph (log scale) at chosen line cuts (marked by the lines in the central picture).The dotted line shows U+ given by (39), full line shows the computed results.

Fig. 6 .
Fig. 6.Computational simulation, the flow along the flat plate and along the rough wall (z S = 0.0002) simulation, regime 15 m s −1 , isolines of the turbulent kinetic energy k and dissipation ω.The wall distancek (green) and ω (red) graphs at the chosen line cuts (marked by the lines in the central picture).

Fig. 7 .
Fig. 7.The flow inside the simple channel with bump, regime 15 m s −1 , isolines of u velocity component and k turbulent kinetic energy.The y − , u, p graphs at chosen line cuts (shown red).
forms the boundary ∂D i .By n i j let us denote the unit outer normal to ∂D i on Γ i j .Let us construct a partition 0 = t 0 < t 1 < . . . of the time interval [0, T ] and denote the time steps τ k = t k+1 − t k .We integrate the system (1) over the set D i × (t k , t k+1 ).With the integral form of the equations we can study a flow with discontinuities, such as shock waves, too.
The following notation is used (some values are used only for an INLET case)

7 Boundary Condition for the Imermeable Wall
Slip wall (Eulerian, inviscid, symmetric boundary) Using this boundary condition, we set the normal component velocity at the wall to zero.The conservation laws must be satisfied in close vicnity to the boundary face.We use the analysis of the incomplete Riemann problem to construct the values for the density and the pressure at each boundary face.The following notation is used At the vicinity of the boundary we solve the following problem (modified Riemann problem for the split Euler equations), with the left-hand side initial condition and the complementary condition (zero normal velocity at the boundary)