The simulation of the flow through the porous media and diffusible barriers

Here we work with the RANS equations describing the non-stationary viscous compressible fluid flow. We focus on the numerical simulation of the flow through the porous media, characterized by the loss of momentum. Further we simulate the flow through the set of diffusible barriers. Here we analyze the modification of the Riemann problem with one-side initial condition, complemented with the Darcy’s law and added inertial loss. We show the computational results 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.We decided to use the analysis of the exact solution also for the discretization of the fluxes through the boundary edges/faces and on the edges/faces simulating the diffusible barrier, as was shown in [1].We construct own algorithm for the solution of the boundary problem at the diffusible barrier, and we use it 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 a e-mail: kyncl@vzlu.cz 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.
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 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 gravitational acceleration vector.
EPJ Web of Conferences 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.

Porous media simulation
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.

Numerical method
For the discretization of the system we proceed as described in [3,4].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 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.
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 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 [5].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 [5][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 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 [3].

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 (15) 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 [5, 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 [6, 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 [6, page 396].By the solution of the problem (15),( 16),(17) 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. [5], [6], [7].The general theorem on the solvability of the Riemann problem can be found in [5, page 88].
The solution is piecewise smooth and its general form can be seen in figure 1, 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 figure 1): Using the theory in [5], [6], [7], 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 (18) -( 23) 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 ( 18), (21) 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 [5].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 ( 18)-( 23) to obtain the whole solution.
In our work we use the analysis of this problem also for the solution of the initial-boundary value at the boundary faces, see also [8][9][10][11].

Boundary condition for the diffusible barrier
Here we present the diffusible barrier condition as a combination of Darcy's law with the additional innertial losses. here We are interested in the boundary values R, U, P at the barrier.Further we require that the conservation laws (1D Euler equations) are satisfied in the close vicinity of the barrier (28) Here the axis x1 is perpendicular to the barrier, ( x1 , t) denotes the density, p( x1 , t) is the pressure, u( x1 , t) is the velocity (with the direction perpendicular to the barrier), E( x1 , t) denotes the total energy: E = u 2 /2 + p/(γ − 1).The initial condition is formed by the two states near the barrier.Let us denote these states 1 , u 1 , p 1 and 2 , u 2 , p 2 : In Section 5 we have shown the unique solution (entropy weak) of the Riemann problem (28),(29).Adding another condition (27) into this complete system doesn't make sense.Therefore we think of the barrier problem as of two boundary problems with the two particular solutions (different in general).
-Problem 1. Find the solution R 1 , U 1 , P 1 of the system (28), equipped with the one-side initial condition 1 , u 1 , p 1 and the boundary condition ( 27).-Problem 2. Find the solution R 2 , U 2 , P 2 of the system (28), equipped with the one-side initial condition 2 , u 2 , p 2 and the boundary condition (27).
Further we will show the solution of these two problems.The analysis of these problems was shown in [1], here we write the resulting algorithm.We suppose the initial velocity u 1 ≥ 0, and we seek the solution with U > 0.
INLET to the diffusible barrier

Examples
Here we show the computational results of the compressible gas flow, obtained with the own-developed codes (the RANS with the k-ω model of turbulence).The gas constants for the air are γ = 1.

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 (analyzed by authors in [1]) is shown together with the computational algorithm.The work is based on the analysis of the Riemann Problem for the split Euler equations and the modifications of this problem.Here the left hand side initial condition is replaced by various complementary conditions.All codes were implemented into the own-developed software.The numerical examples were presented.

Fig. 2 .
Fig. 2. Algorithm for the solution of the Problem 1 (INLET to the barrier).Both posibilities are taken into account.