The simulation of the atmospheric flow with emmisions

The aim of this work is to simulate the flow of the passive gas mixture and estimate the resulting concentration of the potentionally dangerous pollutant emmisions. We assume non-stationary compressible gas flow in the gravitational field described by the system of the RANS equations, equipped with the equation for the concentration of additional specie. Due to gravitational force effect we modify the boundary conditions and data resulting admissible situations. The real simulations of the atmospheric flow require the use of the wall functions for the boundary condition simulating the ground surface with given roughness. We show the modification of such wall functions. The computational results are obtained with the own-developed CFD code for the compressible turbulent mixture flow. The originality of this work lies with the special handling of the boundary conditions, and own computational code.


Introduction
The main topic of this work is to numerically simulate 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. The expected result is to construct a method for the estimates of the pollutant concentrations in the case of a sudden leakage of the substances hazardous to health. This contribution follows previous work [1][2][3].

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 = ( 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 the total energy. Further τ i j = (µ+µ T )S i j −δ i j 2 3 k, 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 + * e-mail: kyncl@vzlu.cz * * e-mail: pelant@vzlu.cz 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 species mass fractions Y 1 , 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).

Numerical method
For the discretization of the system we proceed as described in [8,9]. We use either explicit or implicit finite volume method (FVM) to solve the systems sequentionally. Other possible discretizations were shown in [17,18]. Here we present the discretization of the system (1) using FVM. By Ω h let us the denote the polyhedral approximation of Ω. The system of the closed polyhedrons with mutually disjoint interiors . .} 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 ele- 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 Here the set With the integral form of the equations we can study a flow with discontinuities, such as shock waves, too.
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 D i Further we proceed with the approximation of the fluxes. Usually the flux 3 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 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 [10]. 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 Newton method. We may proceed as shown in [10][page 216]. At each time level t k+1 the following iterative scheme is applied: Here DΦ(w) Dw is the Jacobian matrix of the mapping Φ. Using such process we obtain the sequence w k+1,r i converging to w k+1 i as r → ∞. 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 [8].

The gas mixture equations discretization
As noted before, we will solve the system (1) and the equation (2) sequentially at each time instant. In order to determine the unknown mass fraction Y 1 we consider the following system Using the FVM discretization, shown above, we obtain the following implicit scheme for the unknown ∆w = w| n+1 − w| n at each time instant t n Here H| n i j is the approximation of the fluxes 3 s=1 f s (w)n i j s | n i j through the face Γ i j at time instant t k . Using the Vijayasundaram flux we can approximate these fluxes as follows: 3 s=1 f s (w)n i j s n i j = u | n i j · n i j w| n i j ≈ H| n i j H| n i j = H(w n i , w n j , n i j ) = a + w n i + a − w n j , where a + = max(u| n i j · n i j , 0), a − = min(u| n i j · n i j , 0) are computed constants. The values u| n i j are either known (from the Riemann flux across the face Γ i j ), or we can choose u| n i j = (u| D i + u| D j )/2. The Jacobians of such mapping H(wI, wJ, n i j ) are The symbol R denotes the approximation of the viscous flux 3 s=1 I R s (w, ∇w)(n i j ) s . By ∂R DwI n i j , ∂R ∂wJ n i j we mean the Jacobian of the mapping R(wI, wJ, ...) dependent on the states wI, wJ at the neighbouring elements over the face Here µ T is the approximation (known constant) of the viscosity µ T | n i j at the face Γ i j at time instant t n . The gradient ∇Y 1 | n i j at the face Γ i j can be approximated as ∇Y 1 , using the m-point rule.
The system (16) can be written in a matrix form with and A Y is the matrix with We use the GMRES method with the ILUK preconditioning in order to solve this system.

Boundary conditions
In order to to get the state variables at each boundary face we solve the local boundary problem with the use the original analysis of exact solution of the Riemann problem. This approach was shown and described in [4,16] and analyzed also in [9], [13]. Using the thorough analysis of the Riemann problem we have shown, that the RIC for the local problem can be partially replaced by the suitable complementary conditions. We suggest such complementary conditions accordingly to the desired preference. This way it is possible to construct the boundary conditions by the preference of total values, by preference of pressure, velocity, mass flow, temperature. Further, using the suitable complementary conditions, it is possible to simulate the flow in the vicinity of the diffusible barrier. On the contrary to the initial-value Riemann problem, the solution of such modified problems can be written in the closed form for some cases. Moreover, using such construction, the local conservation laws are not violated.

Wall functions
In the case of the coarse mesh with the first cell in the log layer we use the so-called wall functions for the values inside the first cell.

Standard wall functions
We estimate the friction velocity u τ , assuming the law of the wall (see [3,7]) to be valid at the closest volume, i.e. solve the implicit equation The wall friction τ w = µ ∂U ∂y | w is then estimated as τ w = w u 2 τ . Measurements indicate C ≈ 5.0 for smooth surfaces, κ ≈ 0.41 (Karman's constant), y P denotes the distance from surface (of the first cell). For the rough surface, with elements of average height k S , experiments show that the C is a function of k S , Here k + S is dimensionless roughness height. The momentum equations are solved with the modified effective viscosity µ e at the wall: µ e U P y P = τ w . The values for the k, ω at the first cell are set using log layer equations The main drawback of this simple wall function are the cases with separated flow where τ w → 0, consequently giving k p → 0.  Launder-Spalding/Chieng-Launder methodology Estimate the friction velocity using the known value k P at the first cell: (u * τ ) 2 = k P √ β * , then estimate shear stress τ w using modified log layer equation The momentum equations are solved with the modified effective viscosity µ e at the wall: µ e U P y P = τ w , i.e. µ e = w u * τ y P 1 κ ln(u * τ y P w /µ w ) + C .
Further, the production term P k = τ w ∂U ∂y (in the log layer) can be estimated as Integrating over the log layer in the first cell (height y N ), the average value at this cell is then estimated Here y V is the sublayer thickness, typically set such that y V √ k P w /µ w = 20. We use other estimate y V = 2µ w κ The average dissipation term β * ωk at the near-wall cell is then estimated using viscous and log layer approximations for k, ω Then solve the turbulent kinetic energy with modified integrated production terms. At last set ω at the first cell as Another modification is to set average omega production P ω at the first cell as P ω = α ω P k ω/k.

Examples
All the results and mesh were obtained with the owndeveloped codes. The aim of the following computations was to achieve relatively quickly some estimate of the gas dispersion in the case of some possibly dangerous pollution. Let the surface data be given by a set of point coordinates (left picture) and scaled. Then the surface mesh is constructed, and volumes created using own softwar. Computation run with the initial condition with the total quantities θ 0 = 273.15 K, p 0 = 101325 Pa, and regime velocity v 1 = 5 m s −1 , and fixed emission source located at the chosen point. The figure 4 shows the simulation of the gas mixture flow over such simple terrain.

Conclusion
This paper is focused on the numerical simulation of the mixture of two inert perfect gases in 3D. The finite volume method is applied for the solution of the system of equations. The modification of the Riemann problem and its solution was used at the boundaries. All codes were implemented into the own-developed software. The numerical examples were presented. The original intention was to use the modified software for the quick estimation of the gaseous pollution of the air. This may be critical in the case of the sudden leakage of the substances hazardous to health. The estimates computed with this software are undoubtedly more precise than a simple set of concentric circles. Further development and comparison with the ex-perimental data (from the wind tunnel) is still in the process.