NUMERICAL SIMULATION OF THE FLOW AROUND DIFFUSIBLE BARRIERS

We work with the system of equations describing the non-stationary compressible fluid flow, and we focus on the numerical solution of these equations, and on the boundary conditions. In order to construct the boundary conditions we analyze the equations in the close vicinity of the boundary. Many boundary conditions (i.e. fixed, linearized) do not solve the local problem in an appropriate way, which may bring non-physical errors into the solution, slow down the convergent process, or even ruin the solution in the whole domain. We use the analysis of the exact solution of the Riemann problem in order to solve this local boundary problem. In the case of the diffusible barrier we solve the modified Riemann problem, with the one-side initial condition, complemented with the Darcy’s law and added inertial loss. We analyze the solution of this problem. We show the computational results of the flow through the diffusible barrier, 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 numerical solution of these equations.The correct design of the boundary conditions plays also an important role in the numerical modeling of the processes involved.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.One of the most accurate method (and perhaps the most accurate method) is based on the solution of the so-called Riemann problem for the 2D/3D split Euler equations.Unfortunately, the exact solution of this problem cannot be expressed in a closed form, and has to be computed by an iterative process (to given accuracy).Therefore this method is also one of the most demanding.Nevertheless, on account of the accuracy of the Riemann solver, 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.The incomplete local Riemann problem at the diffusible barrier is complemented with the Darcy's law and added inertial loss.In this work we analyze the modified local problem.We construct own algorithm for the solution of the boundary problem at the a e-mail: kyncl@vzlu.czb e-mail: pelant@vzlu.czdiffusible barrier, and we use it in the numerical examples.

The Navier-Stokes Equations
Here we show the so-called Navier-Stokes equations.We consider the conservation laws for viscous compressible flow with the zero outer volume forces and heat sources in a domain Ω ∈ I R N , and time interval (0, T ), with T > 0. Using the notation in [2, page 376], we write ∂I R s (w, ∇w) ∂x s in Q T = Ω × (0, T ).
(1) Here N denotes the dimension (N = 2 or 3 for 2D or 3D flow, respectively), w = w(x, t) is the state vector, x ∈ Ω, t denotes the time, Q T is called a space-time cylinder, f s are the inviscid (Euler) fluxes, I R s are viscous fluxes: with u = (v 1 , . . ., v N ) T denoting the velocity vector, being the density, p the pressure, θ the absolute temperature, E the total energy, k the heat conduction coefficient, and µ, λ denoting the viscosity coefficients.In this work we assume the relation 3λ + 2µ = 0 is valid.We add the thermodynamical relations to the system (1) with the γ denoting the Poisson adiabatic constant and c v the specific heat at constant volume.Within the numerical solution of this problem we usually consider the equations in the more general integral form, which allows the discontinuities in the solution.It is well known that the classical Navier-Stokes equation for three-dimensional incompressible viscous flow is invariant under the Galilean transformations, see [1, page 69].We can choose a new Cartesian coordinate system ( x1 , x2 , x3 ) such that the origin is at the point σ, and the axis x1 is in the direction of a chosen unit vector n = (n 1 , n 2 , n 3 ), |n| = 1.This transformation can be written Here Q 0 is the rotation matrix, the elements of the rotation matrix are not all independent.The vectors n, o, p must have the properties We can choose the vector o And then the vector p = (p 1 , p 2 , p 3 ) is determined by p = n × o.The axis x2 points in the direction of the vector o, and the axis x3 in the direction of the vector p.
The matrix Q 0 is orthogonal, and the vectors are orthogonal unit vectors.The inverse of the orthogonal matrix Q 0 is The transformation of the state vector w yields the state vector q in the new coordinates.It is Using this transformation of the coordinates, the Navier-Stokes equations keep the form

The Euler Equations
Ommiting the vicsous terms in (1) we get the system of the Euler equations The system of the Euler equations is rotationally invariant (see e.g.[2, page 108]).

FVM Discretization of the Problem
Here we describe the so-called finite volume discretization of the problem (1) in the domain Ω.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 .
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 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 H(w l i , w l j , n i j ) at suitable time instant t l .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 in (5) with the one-point rule.It is possible to approximate the system (5) by the following explicit finite volume scheme Here ∇w l Γ i j denotes the ∇w at the center of the edge Γ i j at time instant t l .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 EFM 2012 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 [2].The approximation of the values ∇w l Γ i j can be done via the dual mesh, see [3].The crucial problem of our work lies with the evaluation of the edge values w k Γ i j .

Inner Faces
To approximate the values w k Γ i j at the inner faces/edges we solve the simplified system (8) 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 (4) 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 .Then the derivatives with respect to x2 , x3 are neglected and we get the so-called split 3D Euler equations, see [2, page 138]: The values w k i and w k j adjacent to the face Γ i j are known, forming the initial conditions (10) The transformation matrix Q is defined in (3).In this work we will refer to these initial conditions as to the left-hand side initial condition (9) and the right-hand side initial condition (10).The problem (8), ( 9), (10) has a unique "solution" in (−∞, ∞) × (0, ∞), the analysis can be found in [2, page 138].The solution of this problem is fundamentally the same as the solution of the Riemann problem for the 1D Euler equations, see [2, page 138].In fact, the solution for the pressure, the first component of the velocity (in the direction of the axis x1 ), and the density is exactly the same as in one-dimensional case.The remaining components of the velocity change only across the middle wave (described in Section 4) in general.Let q RS (q L , q R , x1 , t) denote the solution of the problem (8), (9), (10) 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 (3) of the state vector q Γ i j into the global coordinates is The described process of finding the face values w k Γ i j is independent on the choice of the vectors o, p, determining the transformation matrix Q defined in (3), see [3].

Boundary Faces
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 (8) in the vicinity of the face Γ i j in time with the initial condition (9) determined by the state vector w k i .Similarily to the situation at the inner face we use the rotational invariance of the Euler equations, and we express the system (4) in a new Cartesian coordinate system x1 , x2 , x3 with the origin at the center of gravity of Γ i j and with the new axis x1 in the direction of vector n = (n 1 , n 2 , n 3 ).The vector n corresponds to the face normal vector n i j pointing out of the volume D i , however for certain boundary conditions it is convenient to define n in a different way (still pointing out of the volume D i ).Then the derivatives with respect to x2 , x3 are neglected and we get the so-called split 3D Euler equations (8).The state vector w k i on the element D i adjacent to the face Γ i j is known, forming the initial condition for the local problem (8),(9).We seek the boundary state q Γ i j = q(0, t).This problem is not uniquely solvable and it must be further specified with some chosen complementary conditions, see [3].The boundary state q B is then found as a solution of the problem (8),(9) and the complementary conditions.In order to construct such complementary conditions we use the analysis of the 1D Riemann problem in Section 4, see [6], [3].For example, it is possible to prefer the given pressure, or the given velocity, or the given temperature, or the total quantities and the direction of the velocity, or the given mass flow at the boundary.In this work we focus on the boundary condition at the diffusible barrier.Section 5 shows the solution of this problem for the pressure, density, and the first component of velocity at the barrier (in local coordinates), the other components of velocity are chosen from the known initial condition.Once the local state vector q Γ i j is known, we use transformation (3) to evaluate the face state vector

1D Riemann Problem for the Euler Equations
We are concerned with one-dimensional Euler equations with the notation q = ( , u, E) T , f (q) = ( u, u 2 + p, (E + p)u) T , u denotes the velocity, the density, p the 01054-p.3 pressure, E denotes the total energy.We assume that the equation of state for ideal gas holds p = (γ − 1) Here γ is the Poisson adiabatic constant.
The system is hyperbolic.The Riemann problem for hyperbolic system (12) consists of finding its entropy weak solution in Q ∞ , which satisfies the initial condition formed by two constant states q L , q R : Here The physical analogue to this problem is so-called shock-tube problem.It is more convenient to use the vector of primitive variables ( , u, p) rather then the vector of conservative variables ( , u, E) in solving the Riemann problem.The general theorem on the solvability of the Riemann problem can be found in [2, page 88].The solution q = q( x1 , t) of the Riemann problem (12),( 13),( 14) is piecewise smooth and its general form can be seen from Fig. 3, where a system of half lines is drawn.These half The solution in Exact solution of the Riemann problem has three waves in general.The wedges Ω L and Ω L are separated by left wave (either 1-shock wave, or 1-rarefaction wave).There is a contact discontinuity in between regions Ω L and Ω R , it is formed by the half line ( x1 , t); x1 t = u ; t > 0 .Wedges Ω R and Ω R are separated by right wave (either 3-shock wave, or 3-rarefaction wave).For special state q R there is no right wave and it is q R = q R .The pressure p and the velocity u don't change across the contact discontinuity, while the density changes across this discontinuity in general.We can denote the vectors of primitive variables in particular regions as Here we show the relations between primitive variables in Ω L ∪ Ω HT L ∪ Ω L (across the left wave), and the set of relations for primitive variables in Ω R ∪Ω R ∪Ω HT R (across the right wave).
-1-shock wave One of the possible wave patterns connecting Ω L and Ω L is a shock wave.Region Ω HT L degenerates into single half-line.Primitive variables , u, p "jump" accross this wave.Inviscid shock jump relations can be derived, we call them Rankine-Hugoniot relations.These leads us to the following relations across the 1-shock wave (see [5], or [2, page 125]) where a L = γ p L ρ L is the speed of sound in the Ω L , s 1 denotes speed of the 1-shock wave.Half line x1 t = s 1 shapes the boundary between Ω L and Ω L .With the (13) known, (15),( 16),(17) form the system of three equations for four unknowns L , p , u , s 1 .It is u < u L , p > p L , L > L across this wave, see [ Here s T L is speed of the tail of the 1-rarefaction wave.
From the speed of sound definition a = γp and (19) we get another useful relation across the 1-rarefaction wave It is a L ≤ a L for 1-rarefaction wave.Speed of the head of the rarefaction wave can be expressed s HL = u L −a L .Half line x1 t = s T L forms the boundary between Ω L and Ω HT L .Primitive variables in Ω HT L change continuously, see [2], (3.1.97),page 118.Pressure positivity for p in (18) gives the condition u < u L + 2 γ−1 a L . 01054-p.4 -3-shock wave One of the possible right wave pattern is a 3-shock wave.Region Ω HT R degenerates into a single half-line.
The primitive variables , u, p are discontinues across this half line.Following relations for the primitive variables hold Here s 3 denotes the velocity of the 3-shock wave,

rarefaction wave
Another possible wave pattern forming the region Ω HT R is the rarefaction wave.Primitive variables change smoothly across this wave.It holds that In the last relation s T R denotes the speed of the tail of the 3-rarefaction wave.Half line S HR = {( x1 , t); x1 = s HR t; t > 0} is the boundary between Ω HT R and Ω R , and half line S T R = {( x1 , t); x1 = s T R t; t > 0} is the boundary between Ω R and Ω HT R .The velocity s HR can be expressed It is p ≤ p R , u ≤ u R across the 3-shock wave.Primitive variables in Ω HT R change smoothly, see [2], [7].
Now, the combination of ( 15), ( 18), ( 22) and ( 25) gives the implicit equation for the unknown pressure p , the analysis can be found in [2].This is a nonlinear algebraic equation, and one cannot express the analytical solution of this problem in a closed form.The Newton method can be used to find the solution for the pressure p .With the pressure p known, it is possible to use the equations ( 15)-( 27) to compute the remaining unknowns in the solution of the Riemann problem (12),( 13),(14).

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 (31) 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 4 we have shown the unique solution (entropy weak) of the Riemann problem (31),(32).Adding another condition (30) 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, U, P of the system (31), equipped with the one-side initial condition 1 , u 1 , p 1 and the boundary condition (30).-Problem 2. Find the solution R, U, P of the system (31), equipped with the one-side initial condition 2 , u 2 , p 2 and the boundary condition (30).
R, U, P R, U, P Further we will show the solution of these two problems.
Remark 1.The condition (30) is not rotationally invariant!Let us think of a Problem 1. with the given initial state 1 , u 1 , p 1 = ˜ , ũ, p and p 2 = p X , and let R 1 , U 1 , P 1 be the solution satisfying ).This situation is symmetric to the Problem 2. with the given 2 , u 2 , p 2 = ˜ , −ũ, p and p 1 = p X .Let R 2 , U 2 , P 2 be the solution of such problem satisfying . This solution should be symmetric (to the solution of the previous Problem 1.): R 2 , U 2 , 01054-p.5 Remark 2. Non-uniqueness of the solution: There is a possibility of multiple solutions to the Problem 1.Let u 1 > 0, and let us choose the velocity U < 0 arbitrarily.Then there is a solution with the 1-shock wave, and the sought boundary state is R, U, P where . This relation, shown in [3], is rewritten relation (15).Further (due to the possible density jump across the contact discontinuity) we choose density R = (d(p 1 −p 2 )−CU)/U 2 .
As a consequence of the remarks above, we limit the Problem 1. in order to achieve the uniqueness of the solution.We suppose the initial velocity u 1 ≥ 0, and we seek the solution with U > 0. According to the analysis of the Riemann problem we use the equations for the left wave and the complementary condition (30).There are two possible patterns of the left wave.
1-shock wave: u 1 ≥ U ≥ 0 In this case we use the formulae (15),( 16) for the shock wave, valid for P ≥ p 1 .The velocity of the shock wave s 1 = u 1 − a 1 (P)/ 1 , see (17).We require the boundary value condition (30) to be applicable.So it must be s 1 ≤ 0. Using (17), this is satisfied only if where c 1 = γp 1 / 1 is the speed of sound, and M denotes the Mach number.If u 1 < c 1 , then s 1 ≤ 0 for every P ≥ p 1 .
The function U(P) is decreasing, let us denote P o the pressure satisfying U(P o ) = 0, it is P o ≥ p 1 .Now we define a function l(P) = R(P) U 2 (P) + C U(P) in the interval p 1 , P o , using the formula (33),(34).We have got l a (P) = (a 1 (P)u 1 − (P − p 1 )) 2 P(γ − 1)/2 + p 1 (γ + 1)/2 , l b = C u 1 − P − p 1 a 1 (P) where l(p 1 ) = 1 u 2 1 +Cu 1 , l(P o ) = 0. We want to prove that the function l(P) is monotonous in the interval p 1 , P o .The function l b (P) is decreasing.We want to prove, that the function l a (P) has the same property in p 1 , P o : It is difficult to guess the sign of the first derivative of this function except l a (p 1 ) = u 1 1 γp 1 (u 1 − 2c 1 ) and l a (P o ) = 0.
The added inequality u 1 < c 1 implies l a (p 1 ) < 0. To study the monotony of l a (P) we use the following substitutions h 2 (P) = l a (P), l a (P) = 2 h(P) h (P), where h(P) > 0 in p 1 , P o ), and h(P o ) = 0. We use following trick to find the sign of h (P).We use the next substitution to the formulae for h (P) It yields 1-rarefaction (expansion) wave: u 1 ≤ U < u 1 + 2 γ−1 c 1 In this case we use the formulas for the 1-rarefaction wave (18), (19).Using these, the relation for the density R can be rewritten Let us define the function s(U) = U 2 R(U)+CU dependent on the velocity U: ).In order to apply the boundary value condition (30) it is necessary to study the velocities 01054-p.6The results have been proved to be unique.

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

Conclusion
This paper shows the analysis of the boundary condition for the diffusible barrier.This boundary condition is then used in the finite volume method for the simulation of the compressible gas flow through the diffusible barrier.The numerical examples of the two-dimensional flow around diffusible barrier were presented.

RU 2 +Fig. 7 .
Fig. 7. Diffusible barrier, R, U, P is the sought state at the boundary, 2 , u 2 , p 2 , p 1 present the given initial condition, situation with the rarefaction wave.
4 and R = 287.04[J/kgK].We use the viscosity coefficient (depending on the local temperature) computed with the use of the Sutherland's law µ = µ re f T T re f 3/2 T re f +S T +S , with S = 110.4[K],T re f = 273.15[K],µ re f = 1.716e − 05[Pa s], and the heat conduction coefficient k = 0.0211.The Figs. 9. -11.show the computed result of the flow through the diffusible barrier enclosed in the channel.The following setup was used: geometry -the given geometry is shown in Fig. 8. initial condition -constant state in the whole domain, T o = 273.15[K],u = (0, 0), p o = 101325[Pa].inlet boundary condition (left) -total quantities and velocity direction.This boundary condition was described in [4], [3], total pressure p o = 101325[Pa], total temperature T o = 273.15[K],and zero tangential velocity.outlet boundary condition (right) -preference of the pressure, see [3].Given pressure p = 95009.2[Pa]and averaging technique was used.wall boundary condition (up, bottom) -preference of the velocity, see [3].Zero velocity at the wall was given together with the wall temperature T WALL = 273.15[K].diffusible barrier -boundary condition described in Section 5. with constants C 0 = 10 2 , ∆m = 0.01, α = 10 −8 .
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 .