Simulation of the Propeller Disk Inside the Symmetrical Channel

We work with the system of equations describing non-stationary compressible turbulent fluid flow, and we focus on the numerical solution of these equations, and on the boundary conditions. The computational simulation of the propeller disk is a demanding and time-consuming task. Here the propeller disk is represented by the distribution of the vector of velocities along its radius. The main purpose is to describe the special compatible conditions used to simulate the propeller disk on the both its sides. In order to construct these conditions we analyze the equations in the close vicinity of the boundary. We use the analysis of the exact solution of the Riemann problem in order to solve this local boundary problem. The one-side modification of this problem has to be complemented with some other conditions. At the back side of the propeller disk, it is advantageous to use total density and the total pressure distribution, coming from the known distribution of axial velocities on the disk and the total state values at the inlet, and extra added velocities of rotation. At the front side of the disk, it is preferable to use the distribution of the flowing mass, known from the state values computed on the back side of the disk. We analyze the solution of these particular problems. We show the computational results of the flow around such propeller disk, obtained with the own-developed code for the solution of the 3D axis-symmetrical compressible turbulent gas flow.


Introduction
The numerical simulation of the propeller disk is always a difficult, and time-demanding task.For the proper simulation it is usually needed to work in 3D, using some sofisticated meshes.In this paper we work with the simple rotor, represented by the distribution of the velocity.Furthermore we assume the problem to be axis-symmetrical.We simulate this rotor with the use of two connected boundaries.At these boundaries we solve the conservation laws, using the modification of the Riemann problem by the preference of the total quantities at the inlet, and the modification by the preference of the mass flow at the outlet.Numerical example of this approach is shown.In our simulation we consider the Reynolds-Averaged Navier-Stokes equations with the k-ω model of turbulence, rewritten into the axissymmetrical form.

Formulation of the Equations for an Axis-symmetrical Flow
For a symmetrical three dimensional flow we use the following system of the equations where q = ( , u, v, w, E) f (q) = u, u 2 + p, uv, uw, (E + p)u a e-mail: kyncl@vzlu.czb e-mail: pelant@vzlu.czHere p is the pressure, the density, (u, v, w) is the average value vector of velocity: u is the velocity in the direction x, the components v, w are radial and circle velocities, x, y, z denote the cylindrical coordinates: y denote the radius, z the angle of rotation, and t the time.Further, k is the turbulent kinetic energy of flux components of the velocity, ω is the specific turbulent dissipation, P r is laminar and P r T is turbulent Prandtl constant number, μ is the dynamic viscosity coefficient dependent on temperature, μ T = k/ω is the eddy-viscosity coefficient.In the energy equation, E denotes the total energy where ε = p/ (γ − 1) is the internal energy of a unit mass of the fluid where the constant γ > 1.This system of equations ( 1) is an open system for turbulent flow.The system studied (1) can be rewritten into the differential symbolic form and in the integral form it reads where i = 1, 2, 3, 4, 5, Ω is from the space R 3 (t, x, y), and (, ) denotes the scalar product.n is a normal vector to ∂Ω.The positive orientation is given by the outward direction.Here s is the integral measure in the surface ∂Ω.Using the integral form we can study a flow with shock waves, too.

Modification of the k -ω Standard Turbulence Model for Axis-symmetrical Flow
Turbulent model is described by following equations where k the turbulent kinetic energy and ω the turbulent dissipation are functions of time t and space coordinates x, y.The production terms P k and P ω are given by formulas Here functions τ are defined in Sect. 2. for μ = 0, and σ k , β * , β, σ ω , κ are given constants, see [1].This turbulent model k-ω (3), (4) with equations (1) presented the closed system of equations.

1D Riemann Problem for the Euler Equations
We are concerned with one-dimensional Euler equations.
with notation q = ( , u, E) T , f (q) = ( u, u 2 + p, (E + p)u) T .We assume equation of state for ideal gas holds u denotes velocity, density, p pressure.E denotes total energy.The system is hyperbolic.The Riemann problem for hyperbolic system (5) consists in finding its entropy weak solution in Q ∞ , which satisfies the initial condition formed by two known constant states q L , q R q(x, 0) The physical analogue to this problem is so-called shocktube problem.
The general theorem on the solvability of the Riemann problem can be found in [4, page 88].The solution q = q(x, t) of Riemann problem (5),(6),(7) is piecewise smooth and its general form can be seen from figure 1, where a system of half lines is drawn.These half lines define regions, where q is either constant or smooth.Let us define the open sets called wedges (see figure 1): The solution in Ω L , Ω L , Ω R , Ω R is constant (see e.g.[4, page 128]) Exact solution of the Riemann problem has three waves.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 .Wedges Ω R and Ω R are separated by right wave (either 3-shock wave, or 3-rarefaction wave).Therefore there are four possible wave patterns in the solution.
-Contact Discontinuity Pressure p and velocity u don't change across contact discontinuity, which separates Ω L and Ω R region in general, and moving at speed In general, there is a discontinuity in density across half line It is clear, that there must be a jump in temperature, if there is a jump in density, but not in pressure.There is also a jump in entropy.Contact discontinuity is sometimes called entropy wave.In reality, contact discontinuities are smeared out due to diffusive effects, which are neglected and not included in hyperbolic system (5).
It is more convenient to use the vector of primitive variables rather then the vector of conservative variables in solving the Riemann problem.We can denote these vectors of primitive variables in particular regions as Here we show the equations for the primitive variables in Ω L ∪ Ω HT L ∪ Ω L .For the solution of the whole system see [4]. - s 1 denotes speed of the 1-shock wave.Half line x t = s 1 shapes the boundary between Ω L and Ω L .It can be shown (see [3]), that (8) can be rewritten in the form where Here a L = γ p L ρ L is the speed of sound in the Ω L .s T L is speed of the tail of the 1-rarefaction wave.Speed of the head of the rarefaction wave can be expressed Half line x t = s T L shapes the boundary between Ω L and Ω HT L .Pressure positivity in (13) gives condition u < u L + 2 γ−1 a L .Equation (13) can be written in the form State variables in Ω HT L changes continuously, and can be expressed using equations (see [4], (3.1.97),page 118.) We combine both possibilities (11), (17) and we get equation for pressure p using values (21) Combining ( 9) and ( 14) we get equation for density L and pressure EFM 2013 Let us denote s L? the "wave speed" as follows: In case of 1-shock wave it denotes s 1 wave speed, in the case of 1-rarefaction wave it is the speed of the tail of the wave s T L .We combine (10) and ( 15) to obtain relation between s L? and p : There are four unknowns in Ω L to resolve in order to get the solution across left wave and the position of this wave.
Relations of these unknowns L , p , u and the wave speed s L? to left state variables L , u L , p L are given by three equations ( 21),( 22),(23).To get state variables across left wave, i.e. in Ω L ∪Ω S T L ∪Ω L , and position of the wave, we have to add another equation into the system (21),( 22),( 23).There is a jump in density across contact discontinuity.This brings another unknown R to unknowns L , p , u ,s L? .We have to add two properly chosen equations into this system of equations in order to get uniquely solvable system for 5 unknowns sought in

Boundary Condition by the Preference of Total Pressure and Total Density
In this section we will analyze the initial-boundary value problem formed by the Euler equations (5), equipped with the left-hand side initial condition (6) and the complementary conditions given as The equations ( 25),(24) come from the idea, that the total pressure p o and the total temperature θ o are known (prescribed) in the region Ω R .The boundary state q B is the solution of the system (5),( 6),( 24

Solution for the velocity u
The condition (24) has the form where 1) .We will use the analysis in Section 4, and discuss the condition (28) together with the equation for the pressure (21) and the condition (26).We get the velocity u solving the equation where Here a L denotes the speed of sound in Ω L , the functions E 1Ls (u) and E 1Lr (u) are defined in (11),(17).We seek the velocity u in the interval (27).
There is no admissible solution of the system (5),( 6), ( 24), (25 , the extreme value yields the zero pressure p = 0.In practical applications, we may prescribe the closest velocity u = − 2a 2 o (γ−1) + , with > 0 being a small positive constant, though it is not the solution of the system.
Let us analyze the function F(u) for the 1-shock and for the 1-rarefaction wave separately.

1-shock wave
For the 1-shock wave it is u < u L .Using the function for the pressure (11), we write

1-rarefaction wave
γ−1 a L we use the velocity function for the 1-rarefaction wave, defined in (17).Let us find the root of the function For u L > 0 there is no solution with the 1-rarefaction wave.
If F(max{u L , − 2a 2 o (γ−1) }) < 0 and F(min{0, u L + 2 γ−1 a L > 0, then there is a unique root of the function (31), see [3].It is 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.
The problem (29) has a unique solution for the initial data satisfying (33) there is no solution of the system (5),( 6),( 24),( 25) with (26).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 (29) does not have a negative solution, and for the practical applications we choose the velocity c no u = root of (30) Fig. 5. Algorithm for the velocity solution u < 0 of the problem (29).In the case of failure, the problem doesn't have a solution, and the values are chosen.

Solution in Ω
Knowing the velocity u we compute the pressure p using (28), and the density L using (22) uniquely.In the case of the 1-rarefaction wave, the solution in Ω HT L is computed using the equations ( 18), ( 19), (20).At last we are interested in the density R .The equation of state holds in the region Further we use the equation for the total temperature (25) and the equation (28).We write .
Here o = p o Rθ o denotes the total density in Ω R .

Solution in
Then it is u = u R , and R = R , see [3].u = u R , and R = R .This way the prescribed total pressure and total temperature are preferred also in Ω R region.
The complementary conditions (24),( 25),( 35) with (26) yield the well-posed inlet boundary condition for the system (5),(6) only if (33) is satisfied.If (33) is not satisfied, then the total quantities cannot be preferred at the boundary.In the practical application prescribe the closest admissible value for the velocity u and compute the boundary state accordingly, or we use the boundary condition preferring the static pressure p = p o , see [2,3], in this case.

Boundary Condition by the Preference of Mass Flow
In this section, we attempt to solve the incomplete modified Riemann problem for the Euler equations (5),(6), seeking the solution in the general form of the Riemann problem solution, consisting of 4 constant states separated by 3 waves.Variables L , u L , p L are known from the initial condition, also the speed of sound a L is known.We add the boundary condition for the mass flow in Ω L given as where G ≥ 0 is given constant (G = mass f low f ace area ).Here u , L denotes unknown parts of the solution in the region Ω L .Velocity u(x, t) and density (x, t) are constant in this region, as was mentioned in 4.
Here we note, that the system (5),( 6),(36) may have any solution for a given G or it may have multiple solutions for the certain choice of G .In this article we will prefer the solutions with the left rarefaction wave.We will transfer the problem (5),(6),(36) into the boundary problem described in [2] with the velocity u given.It is basically the use of the equations (8),(9),(10), or (13),( 14),(15) for known u .If G = 0, then apparently u = 0, and we may use the analysis for the boundary condition by the preference of velocity shown in [2].
-1-rarefaction wave Here we will try to find the primitive variables L , u , p such that (13),( 14),(36) wih u L ≤ u < u L + 2a L γ−1 is satisfied.From (13),( 14),(36) it is We find the velocity u as the root of the function γ−1 − G on the appropriate interval.Then the unknowns L , p are evaulated from (13), (14).It is is zero at the points Derivative F (u) > 0 in the interval (0, u 0 ) and F (u) < 0 for u ∈ (u 0 , u 1 ) .The maximum of the function F is at the point u 0 .The problem of the existence of the roots depends on the sign of F(max(u L , u 0 )).
Let u L < 0. In the case of 1-shock wave it is u < u L , u L < 0, and (36) cannot be satified.Therefore the solution with 1-rarefaction wave is the only solution possible and u can be found as the solution of the (37).If F(u 0 ) < 0 then the problem doesn't have any solution and one may choose the solution with u = u 0 as this gives the mass flow closest to G for u L < 0.
Let u L ≥ 0. If F(max(u L , u 0 )) ≥ 0 then there exists a solution of the problem (5),( 6), (36) with 1-rarefaction wave.With the velocity u known, we use the analysis for the boundary condition by the preference of velocity shown in [2].If F(max(u L , u 0 )) < 0 then there is no solution with 1rarefaction wave and L u L − G < 0. We will show that for L u L − G < 0, u L < a L there is no solution with 1shock wave either.For u L < a L it is s 1 < 0. The Rankine-Hugoniot condition hold across 1-shock wave (see [4, page 120 (3.1.109)])L (u L −s 1 ) = L (u −s 1 ).It is L < L for 1-shock and therefore L u L −G = s 1 ( L − L ) < 0. This is in contradiction with L u L −G > 0. In the case of u L < a L there is no solution with the 1-rarefaction nor the 1-shock wave.We choose u = max(u 0 , u L ).In the case of u L > a L the solution may not exist, though for an existing solution it is s 1 > 0 or s HL > 0. We choose ( B , u B , p B ) = ( L , u L , p L ) at the boundary in this case.

Examples
Here we present the computational simulation of the propeller disk by a given velocity profile.The computational results are shown in figure 6  back side: the total density and the total pressure distribution, and the given velocity of the rotation (tangential components).Here the total pressure p Do and the total density Do (at the given point) are given (constant in time) as

Conclusion
This paper shows the analysis of the boundary condition for the simulation of the propeller disk by the given velocity profile.This procedure is then used in the finite volume method for the simulation of the propeller placed in the three-dimensional axis-symmetrical channel.

) 2 o
Here θ o > 0, p o > 0 are given constants and a 2 o = γRθ o .Further R denotes the characteristic gas constant, and γ is the adiabatic constant.The equations (24),(25) are considered for the velocity u ∈ − 2a 1) .Using (26), the interval of interest is ),(25) at the half line S B = {( x1 , t); x1 = s B t ; t > 0} , s B = 0.At first, we resolve the velocity u in the Star Region Ω R ∪ Ω L , introduced in Section 4.
and figure 7. The results were obtained with the own-developed codes (the RANS with the k-w model of turbulence).The gas constants for the air are κ = 1.4 and R = 287.04J kg −1 K −1 .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.4K, T re f = 273.15K, μ re f = 1.716e − 05 Pa s, and the heat conduction coefficient k = 0.0211.The figure 7. show the computed result of the flow through the channel with the propeller disk.The following setup was used: EPJ Web of Conferences 02064-p.6 geometry -the given geometry is shown in figure 6. initial condition -constant state in the whole domain, T o = 273.15K, u = (0, 0, 0), p o = 101325 Pa.inlet boundary condition (left) -total quantities and zero tangential velocity (no rotation at the inlet).This boundary condition was described in Sect.5, total pressure p o = 101325 Pa, total temperature T o = 273.15K, and zero tangential velocity (v = 0 m s −1 ).outlet boundary condition (right) -preference of the pressure, shown in [2,3].Given pressure p = 101325 Pa was used.wall boundary condition (up, bottom) -preference of the velocity, shown in [2,3].Zero velocity at the wall was given together with the wall temperature T WALL = 273.15K. the propeller disk:

Here M 2
= v 2 ax /(κRT o ) and v ax is the axial component of the velocity at the given point and p o , T o are the values from the inlet boundary condition.front side: the preference of the mass flow, values computed from the back side of the disk.

Fig. 7 .
Fig. 7. Compressible gas flow through the propeller disk, the pressure isolines, the density isolines, the radial and the axial velocity components in the meridian plane.
2 , and u < u L .-1-rarefaction wave Another possible left wave pattern is rarefaction wave.It forms Ω HT L region.Variables changes smoothly within this wave.Across the 1-rarefaction wave p ≤ p L , u ≥ u L .Following equations are true.