Modification of the Riemann problem and the application for the boundary conditions in computational fluid dynamics

We work with the system of partial differential equations describing the non-stationary compressible turbulent fluid flow. 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. The fundamental problem in this area is the solution of the so-called Riemann problem for the split Euler equations. It is the elementary problem of the one-dimensional conservation laws with the given initial conditions (LIC left-hand side, and RIC right-hand side). The solution of this problem is required in many numerical methods dealing with the 2D/3D fluid flow. The exact (entropy weak) solution of this hyperbolical problem cannot be expressed in a closed form, and has to be computed by an iterative process (to given accuracy), therefore various approximations of this solution are being used. The complicated Riemann problem has to be further modified at the close vicinity of boundary, where the LIC is given, while the RIC is not known. Usually, this boundary problem is being linearized, or roughly approximated. The inaccuracies implied by these simplifications may be small, but these have a huge impact on the solution in the whole studied area, especially for the non-stationary flow. Using the thorough analysis of the Riemann problem we show, 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. Algorithms for the solution of the modified Riemann problems were coded and used within our own developed code for the solution of the compressible gas flow (the Euler, the Navier-Stokes, and the RANS equations). Numerical examples show superior behaviour of the suggested boundary conditions. Constructed boundary conditions are robust and accelerate the convergence of the method. The original result of our work is the analysis of various modifications of the Riemann problem and its applications.


Introduction
Undoubtedly, the boundary conditions play important role in the computational fluid dynamics (CFD).We work with the system of equations describing non-stationary compressible turbulent fluid flow, i.e. the Reynolds-Averaged Navier-Stokes (RANS) equations in 2D and 3D.We focus on the numerical solution of these equations, and on the boundary conditions.We suggest the original approach to the boundary conditions.The aim is to satisfy the conservation laws in the close vicinity of the boundary.Usually, the boundary problem is being linearized, or roughly aproximated.The inaccuracies implied by these simplifications may be small, but it has a huge impact on the solution in the whole studied area, especially for the non-stationary flow.In our approach we try to be as exact as possible.Therefore we use the analysis of the so-called Riemann problem for the 2D/3D split Euler equations in order to construct the boundary values (for the density, velocity, pressure).a e-mail: kyncl@vzlu.czb e-mail: pelant@vzlu.cz We solve the boundary modifications of this initial-value problem.The crucial problem of many numerical methods for the solution of the compressible flow lies in the evaluation of the so-called fluxes through the edges/faces Γ i j of the particular elements.The state values in the vicinity of the edge Γ i j are known at time instant t k , and these values form the initial conditions (LIC -left-hand side, and RIC -right-hand side) for the so-called Riemann problem for the 2D/3D split Euler equations.The exact (entropy weak) solution of this problem cannot be expressed in a closed form, and has to be computed by an iterative process (to given accuracy).Therefore various approximations of this solution are usually analyzed.At the boundary faces we deal with the local modified Riemann problem, where the LIC is given, while the RIC is not known.In some cases (far field boundary) it is wise to choose the RIC here as the solution of the local Riemann problem with given far field values, which gives better results than the solution of the linearized Riemann problem, see [3].Another boundary condition based on the exact Riemann problem solution, simulating the impermeable wall on move, was shown in [13, pages 221-225], where the RIC is constructed in a special way, in order to obtain the desired solution.Using the analysis of the Riemann problem we show, that the RIC for the local problem can be partially replaced by the suitable complementary condition.Some of the suggested boundary conditions were shown in [8] (by preference of pressure, temperature, velocity,...).We complement the boundary problem suitably, and we show the analysis of the resulting uniquely-solvable modified Riemann problem.We construct own algorithm for the solution of this boundary problem.It can be used within various methods in CFD.The algorithm was coded and used within our own developed code for the solution of the compressible gas flow (the Euler, NS, and RANS equations).

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 (x, t) = L , u(x, t) = u L , p(x, t) = p L , x < 0, (x, t) = R , u(x, t) = u R , p(x, t) = p R , x > 0. ( 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 (1) 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 [9, 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 [12, 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 [12, page 396].By the solution of the problem (1),( 2),(3) 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. [9], [12], [13].The general theorem on the solvability of the Riemann problem can be found in [9, page 88].
The solution is piecewise smooth and its general form can be seen in Fig. 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 Fig. 1): Using the theory in [9], [12], [13], 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 (4) -( 9) 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 ( 4), (7) 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 [9].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 ( 4)-( 9) 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 (2).And similarily, with p known, only the right-hand side initial condition (3) 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 ( 4)- (6).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 ( 4) can be written as the equation for pressure, see [8] p = P(u ), ( 12) 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.At the vicinity of the boundary we solve the modified Riemann problem for the split Euler equations (1), introduced in Section 2, with the left-hand side initial condition (2) and the complementary conditions (prescribed pressure).
Here p denotes the pressure in Ω L ∪Ω R part of the solution of the Riemann problem for the split Euler equations, shown in Section 2. This way we have prescribed the desired pressure value whenever it is possible (see the general form of the solution in Section 2).Now it is possible to use (4),( 5),( 6) 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 2. 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 [8].
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 2. This way we have prescribed the desired normal component of velocity whenever it is possible (see the general form of the solution in Section 2).Now it is possible to use ( 12),( 5), (6).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. 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) 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. 3.The complete analysis of this problem is shown in [8].

Boundary Condition by Preference of Temperature
Using this boundary condition, we try to prescribe given static temperature at the face, comlete analysis was shown in [8].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, pressure and velocity.The following notation is used (some values are used only for an INLET  At the vicinity of the boundary we solve the modified Riemann problem for the split Euler equations (1), introduced in Section 2, with the left-hand side initial condition (2) and the complementary conditions (prescribed pressure).
Here θ L denotes the static temperature in Ω L part of the solution of the Riemann problem for the split Euler equations, shown in Section 2. This way we have prescribed the desired temperature value whenever it is possible (see the general form of the solution in Section 2).Now it is possible to use (4),( 5),( 6) together with the equation of state 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. 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 use the condition (21), which ensures that p R = p , u R = u , R = R = L .Further it is possible to prescribe the direction of the velocity d Then it is given value for the total temperature (=desired total temperature at the face) At the vicinity of the boundary we solve the modified Riemann problem for the split Euler equations (1), introduced in Section 2, with the left-hand side initial condition (2) and the complementary conditions (prescribed direction of velocity and given total quantities).24),(25) come from the idea, that the total pressure p o and the total temperature θ o are known (prescribed) in the region Ω R .
According to analysis shown in [5], there are two possibillities for the 1-wave, either there is a 1-shock wave and the sought velocity u can be found as a root of the function considered for u < u L , − 2a 2 o (γ−1) < u < 0, and with P S (u) defined in (12).
Or there is a 1-rarefaction wave, and the velocity u can be computed combining (24),( 25), (12) as Summarizing the above possibilities of the 1-shock and the 1-rarefaction wave we can construct the algorithm in figure 5. for the computation of the velocity u in the star region.
c no u = root of (27) Fig. 5. Algorithm for the velocity solution u < 0.In the case of failure, the problem doesn't have a solution, and the values are chosen.
Knowing the velocity u we compute the pressure p using (25), and the density R = p Rθ R .At last we set the boundary values as The problem (1),( 2), (25), (24) has a unique solution for the initial data satisfying u there is no solution.In this case we prescribe the velocity u = − 2a 2 o (γ−1) + , with > 0 being a small positive constant.If F(min{0, u L + 2 γ−1 a L }) ≤ 0 then the problem does not have a negative solution, and for the practical applications we choose the velocity u = min{0, u L + 2 γ−1 a L }, or we use the boundary condition preferring the static pressure p = p o , see [2,8], in this case.
Let p0 , θ0 be the given total quantities together with the given tangential velocity components : v t and w t are the velocity components in the direction of o and p.Here o•n = 0, |o| = 1, and p = n × o.Then we set and we use the same algorithm (figure 5.) for the computation of the velocity u .Then it is p B = p , B = R , and given value for the total temperature (=desired total temperature at the face) At the vicinity of the boundary we solve the modified Riemann problem for the split Euler equations (1), introduced in Section 2, with the left-hand side initial condition (2) and the complementary conditions prescribing the mass flow, direction of velocity, and total temperature in Ω R given as where G ≤ 0 is given constant (G = mass f low f ace area ).Here γ and R are gas constants.The equation ( 31) is considered for the velocity u ∈ − 2a 2 o (γ−1) , 0 .In general, there are two possibilities of the wave pattern, which may interest us.We seek the unknown velocity u as a root of the function with P(u) defined in (12), and The figure 6. shows example of such function for the chosen values.The function F(u) is monotone, and it has a   Let θ0 be the given total temperature together with the given tangential velocity components : v t and w t are the velocity components in the direction of o and p.Here o • n = 0, |o| = 1, and p = n × o.Then we set and we use the same algorithm (figure 7.) for the computation of the velocity u .Then it is p B = p , B = R , and

Outlet Boundary Condition by Preference of Massflow
Using this boundary condition, we try to prescribe given total quantities at the face, comlete analysis was shown in [6,7].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, pressure and velocity.The where G ≥ 0 is given constant (G = mass f low f ace area ).We are interested in the solution with u > 0, and s L? < 0, which guarantees the possibility of the values to be prescribed at the boundary.In general, there are two possibilities of the wave pattern, which may interest us.
-1-shock wave Here, we are interested in the solution with s 1 < 0, u > 0. Let us construct the function F S (u), using the relations (4),( 5) where and P S (u) defined in (12).The sought velocity u is the root of this function F S (u).It is [8].The first derivative is F S > 0 for u < min{u L , u X }.The problem has a solution with the 1-shock wave if F S (min{u L , u X }) > 0. -

1-rarefaction wave
Let us assume the solution with the 1-rarefaction wave.Using (4),( 5),(34), the velocity u is the root of the function Derivative F R (u) > 0 in the interval (0, u 0 ) and F (u) < 0 for u ∈ (u 0 , u 1 ) .The maximum of the function F R is at the point u 0 .We are interested in the solution with s T L < 0. This is satisfied only fo u < u 0 .y

Algorithm
Here we present the possible algorithm for the solution of the sought values at the boundary.

The Shortened Domain Example
The

The Outlet Mass Flow Example
The developed software with presented boundary conditions was used for the simulation of the compressible turbulent flow in the 3D axis-symmetrical channel.Axis

Conclusion
This paper is focused on the boundary conditions for 3D compressible flow, and summarizes the work presented in [2], [8], [5], [6], [7], and other publications.The algorithms for the construction of 3D boundary conditions are presented.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.The solution of such problems is shown.

Fig. 1 .
Fig. 1.Structure of the solution of the Riemann problem (1),(2),(3) ) u < 0.(26) Here θ o > 0, p o > 0 are given constants, R denotes the characteristic gas constant, and γ is the adiabatic constant.The equations (24),(25) are considered for the velocity u ∈ 1) .The equations ( condition, we try to prescribe given total quantities at the face, comlete analysis was shown in[7].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, pressure and velocity.The following notation is used (some values are used only for an INLET case) B static density at the boundary (unknown) p B static pressure at the boundary (unknown) u B velocity vector at the boundary (in global coordinates) (unknown) u B normal component of velocity at the boundary (local coordinates) (unknown) n unit outer normal to face n = (n 1 , n 2 , n 3 ) L static density near the wall, time 0 p L static pressure near the wall, time 0 u L velocity vector, state near the wall at time 0 u L normal component of velocity, state near the wall at time 0, u L = u L • n G given inlet massflow (=desired massflow at the face) θ 0

2a 2 0γ− 1 .
Once the velocity u < 0 is known, it is possible to compute the remaining valuesB = G /u , p B = p(u ), u B = du d • n.

Fig. 7 .
Fig.7.Algorithm for the velocity solution u < 0.In the case of failure, the problem doesn't have a solution, and the values are chosen.
following notation is used (some values are used only for an INLET case) B static density at the boundary (unknown) p B static pressure at the boundary (unknown) u B velocity vector at the boundary (in global coordinates) (unknown) u B normal component of velocity at the boundary (local coordinates) (unknown) n unit outer normal to face n = (n 1 , n 2 , n 3 ) L static density near the wall, time 0 p L static pressure near the wall, time 0 u L velocity vector, state near the wall at time 0 u L normal component of velocity, state near the wall at time 0, u L = u L • n G given inlet massflow (=desired massflow at the face) At the vicinity of the boundary we solve the modified Riemann problem for the split Euler equations (1), introduced in Section 2, with the left-hand side initial condition (2) and the complementary conditions prescribing the mass flow u R = G ,
. e l s e i f u L > u 0 t h e n u = u L e l s e i f ( u 1 < u 0 ) u = u 1 − ε e l s e i f F(u 0 ) < 0 then u = u 0 e l s e u = r o o t o f F R (u), u ∈ (u L , u 0 ), (36) .end e n d i f end end Compute p B = p , B = L u s i n g (12), (5).

9 EXAMPLES 9 . 1
The Flow Around the ObstacleHere we present a computational result of the 2D nonstationary inviscid channel flow at Mach number M = 0.67.A body immersed in the flowing fluid establishes a certain wave pattern which evolves in time and eventually exits the channel.At figure1we show, that the fixed (values are fixed at the boundary) and linearized (as described in [9]) boundary conditions do not give the expected result in time.The inlet is located left, outlet right, other boundaries are considered as wall.The fixed boundary conditions give incorrect results near boundaries.The linearized boundary condition reflects the waves into the domain, leading to the oscillations in the solution.The new suggested boundary conditions do not suffer from these drawbacks.The residual behavior (shown right) demonstrates this result.

Fig. 9 .
Fig. 9. Incompressible flow, body in the channel.Comparison of the various boundary conditions.
following numerical example shows superior behavior of the suggested boundary condition.This is the test case of the the inviscid flow through the Double Circular Arc (DCA) blade cascade DCA08.The blades of this cascade are composed of two circular arcs with the relative thickness 8%.At the inlet we use the boundary condition conserving the total pressure p o = 101325 Pa, the total temperature θ o = 273.15K, and the direction of the velocity α IN = 5.2 • .At the outlet, the outlet boundary condition with the averaging technique described in [8], preferring the pressure p = 45722.351Pa in average.The part of the inlet flow is supersonic.It can be seen in figure10(left), that the used inlet boundary condition does not reflect the shock waves into the computational area.The right picture shows that this boundary condition can be used on the shortened computational domain with the similar result.Constructed boundary condition is robust and accelerates the convergence of the method, it can be also used on the shortened domains.This is the original result of our work.

Fig. 10 .
Fig. 10.The compressible gas flow, transonic regime.Two computations with the same boundary data.The bottom picture shows the computation with the shortened inlet part.Inlet located left.Pressure isolines, density and Mach number isolines, highlighted Mach number 1.

β y 2 s , c ω 6 β
x is the axis of symmetry.Example is shown in figure 11.At the inlet we used the boundary condition shown in section 6, conserving the total pressure, total temperature, and zero tangential velocity, with θ o = 273.15K , p o = 101325 Pa.The boundary condition shown in section 8 was used at the outlet, with the given massflow G = 4.0 kg•s −1 •m −2 .At each face, the value G was computed (in each iteration) in order to match the massflow across the whole boundary.The boundary condition preffering the zero normal velocity was used in the case of the inviscid flow.For the viscous turbulent flow, this condition was modified by the zero velocity at the wall, and wall temperature θ WALL = 273.15K was set.Further k WALL = 0 m 2 •s −2 and ω WALL = c ω 6μ = 120.Here by y s we mean the distance between the face and the center of the neighbouring element.

Table 1 .
The initial data ),(15) 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.L = 1.00, u L , p L = 1.00, p and the solution B , u B , p B for the test problems.

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.Therefore we use the analysis of the incomplete Riemann problem to construct the values for the density and the velocity at each boundary face.The following notation is used (some values are used only for an INLET case) near the wall, time 0 p L static pressure near the wall, time 0 u L velocity vector, state near the wall at time 0 u L normal component of velocity, state near the wall at time 0, u L = u L • n u GIVEN given velocity vector (=desired value at the face) u GIVEN given value for the normal vcomponent of velocity (=desired normal velocity at the face) Lstatic density case) GIVEN given value for the density (=desired INLET density at the face) d given INLET velocity direction d = (d 1 , d 2 , d 3 )

Table 2 .
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.The initial data u Algorithm for the solution of the problem (1),(2),(20),(21),(22) at the half line S B = {(0, t); t > L , θ L and the solution B , u B , p B for test problems with L = 1.25, p L = 100000.

Boundary Condition by Preference of Total Quantities and Direction of Velocity Using
[5]s boundary condition, we try to prescribe given total quantities at the face, comlete analysis was shown in[8],[5].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, pressure and velocity.The following notation is used (some values are used only for an INLET case)

Table 3 .
n 2 u +o 2 v t +p 2 w t , n 3 u +o 3 v t +p 3 w t , ).The initial data u L , p L , G and the solution B , u , p for the test problems with L = 1.25, θ 0 = 273.15,d • n = 1.

Table 4 .
Outlet boundary condition by preference of massflow.
The input data u L , G and the solution (rounded) B , u B , p B , G B = B u B for the test problems with L = 1.25, p L = 100000.