The simulation of the gas flow through the porous media and fences

Abstract. The paper is focused on the numerical simulation of the compressible gas flow through the porous media and fences. We work with the the non-stationary viscous compressible fluid flow, described by the RANS equations. The flow through the porous media is characterized by the loss of momentum. It is possible to use various methods for the simulation of such flow. Here we present the approach with the modification of the source term, and other possibilities using the modification of the face flux. The original approach was presented recently by the authors analysing the modification of the Riemann problem with one-side initial condition, complemented with the Darcy’s law and added inertial loss. Another aim of this paper is the evaluation and estimate of the forces acting on the diffusible barrier (fence) with given parameters. The presented examples were obtained with the own-developed code for the solution of the compressible gas flow.


Introduction
The physical theory of the compressible fluid motion is based on the principles of conservation laws of mass, momentum, and energy. The mathematical equations describing these fundamental conservation laws form a system of partial differential equations (the Euler equations, the Navier-Stokes equations, the Navier Stokes equations with turbulent models). We focus on the flow through the porous media and throug the diffusible barriers (fences). We choose the well-known finite volume method to discretize the analytical problem, represented by the system of the equations in generalized (integral) form. 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. In order to compute these fluxes, various methods can be used. For example it is possible to use the solution of the Riemann problem for the construction of the face fluxes. Recently, in [1][2][3] we have shown the use of the analysis of this exact solution also for the discretization of the fluxes through the boundary edges/faces and on the edges/faces simulating the diffusible barrier (fences). Here we present other more simple method for the construction of the flux through the face representing the diffusible barrier. This method was implemented into own computational code, and used in the numerical examples. Another method for the numerical simulation of the porous media is described in [4]. * e-mail: kyncl@vzlu.cz * * e-mail: pelant@vzlu.cz

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 (1) Here x 1 , x 2 , x 3 are the space coordinates, t the time, w is the state vector, f s are the inviscid fluxes, I R s 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.
Here δ i j denotes Cronecker's delta, and 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 = S g = (0, g 1 , g 2 , g 3 , g · u), where g = (g 1 , g 2 , g 3 ) is the gravitational acceleration vector.

Porous media source term
The porous media is simulated using the modification of the system of equations (1). The simple porous media can be simulated via the new source term, written as Here α is the permeability coefficient, C o is the pressure gradient coefficient. For the simulation of the thin fence with small viscous effect we choose large value of the koefficient alpha = 1e7, which leads to neligible viscous term µ α . Based on the work [5] we may estimate the parameter C 0 based on the fence porosity as 0.01 C 0 = k r = 0.52(1 − η 2 )/η 2 . Here η is the porosity parameter of the barrier. For the 30% barrier we choose η = 0.3. For the barriers 15%,30%,42%,50%,70% we may estimate C 0 = 2300, 525, 250, 150, 55.

Model of turbulence
Here we assume the system (1) equipped with the twoequation turbulent model k − ω (Kok), described in [6]. The effective turbulent viscosity is µ T = k/ω.
∂ ω ∂t 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 σ d = 0.5 is constant.

Numerical method
For the discretization of the system we proceed as described in [7,8]. We use either explicit or implicit finite volume method (FVM) to solve the systems sequentionally. Other possible discretizations were shown in [16,17]. 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 polyhe- . .} 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 ∂D i of each element D i is Here the set Γ D i = {Γ i j ; Γ i j ⊂ ∂D i } 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.
∂I R s (w, ∇w) ∂x s dx dt (6) 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 s=1 f s (w)(n i j ) s | Γ i j is being approximated by a numerical flux at suitable time instant t l 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 [9]. 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 [9][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 [7].

Porous media discretization
Simpliest way to include the porous media within the FVM is to choose the elements (porous area), where the source term (2) is applied. At the face Γ i j simulating the diffusible barrier (fence) we may use the modified face flux as Here |Γ i j | d 2 is the volume of the non-existing (artifitial) element adjacent to the face, representing the diffusible barrier (fence) with diameter d, and are the approximated values of density and velocity vector at the face.

Example
Let us suppose the gas flow with the total temperature θ 0 = 293.15K, total pressure p 0 = 101325 Pa, and given velocity magnitude u ∞ (freestream values). Let us place the fence with the known parameter C 0 into such flow. Our aim is to compute/ estimate the force acting on such fence. The force is computed by integration of the source term S PM at the barrier. Figure 1 shows the graph of the force component F X depending on the barrier parameter C 0 and the velocity regime u ∞ . Based on the CFD simulations, figure 1, it is possible to use rough approximation of the force acting on the barrier as To be more precise, the fluid density should be also included in these simplified relations The estimated values for the coefficients for the considered barriers are shown in table 1.
Further, let us suppose the gas flow above the plate with the 0.1m fence. Let us form the initial and the boundary conditions using the flat plate simulation with the regime θ 0 = 293.15K, p 0 = 101325, u ∞ = 15m.s −1 , stationary state, cut at 70m from the plate start. The regime is shown in figure 2. The resulting velocity is shown in figure 3. The comparison of the CFD computed force and simplified estimation is shown in figure 4.
Further examples involve the body placed behing the fence in various distances. The aim was to estimate the reduction of the force acting on the body behind the barrier with the studied porosities. The figures 5, 6 show the reduction of the horizontal force acting on the body.

Conclusion
This paper is focused on the viscous compressible flow, with the focus on the porous media and the diffusible barrier. The boundary condition for the diffusible barrier basen on the modification of the Riemann Problem was shown in [1][2][3], here we have shown other possibility. We have shown the simple estimate of the force acting on the fence, and further we used onw CFD for the estimation of the forces acting on the body placed behind the barrier. All codes were implemented into the own-developed software.
The numerical examples were presented.
The results originated with the support of Ministry of the Interior of the Czech Republic, project SCENT, Grant MSM 0001066902 of the Ministry of Education of the Czech Republic, and Ministry of Industry and Trade of the Czech Republic for the long-term strategic development of the research organization. The authors acknowledge this support.     6. Regime 1570, the force acting on the body in the distance 5H,10H,15H,20H behind the barrier with porosity 50%, 42%,30%,15% and height H.