In the second test case a ramp with 20 degrees slope angle is located in a flow channel. The air enters into a channel with Mach number 2.85. Before the ramp an oblique shockwave develops. The geometry and Mach number distribution can be found in Fig. 3. The comparison of the measured and simulated results shows similar shock wave pattern, however they can not be compared with each other directly due to slight difference between the inlet conditions and slope values. The reason why the presented condition is used in the simulation is the available measured quantitative parameters found in the Gerolymos' publication (Gerolymos et al., 2003). The locations and directions of the coordinate systems,

software can be accepted for this benchmark.

layer at upstream and the ω.

yEXP [m]

0

0.01

0.02

0.03

along which velocity distributions are measured is shown in Fig. 4.

Fig. 4. Locations and directions of the coordinate systems, along which the velocity distributions of the measurements are compared with the results of simulation

yEXP is the distance from the wall. The velocity profiles of the recent viscous flow solver with k-ω model and the experiments are found in Figs. 5-7. for comparison. The velocity profiles in Fig. 5 and 6. are even quantitatively agreed with each other, but the results shown in Fig. 7. are slightly far from the experiments. The reason of that can be caused by the fact, that the error, which is generated by the velocity profile at the beginning of the computational domain is growing along with the flow and so the small difference becomes larger at downstream. The other problem can be the free-stream sensitivity of the k-ω model described in Subchapter 1.2. The solution could be improved by further adjusting boundary

s = -0.0508 m

experiment k-ω model

Fig. 5. Validation of the viscous flow solver. Velocity profiles at s = -0.0508 m (see Fig. 4.)

The compressible viscous flow solver presented in Chapter 2. requires relative high computational time due to the significant number of equations and fine mesh especially in the boundary layer. Hence, this approach can not be used economically for coupling with optimization methods in the explicit time marching manner. However, assuming a frictionless and convection dominated problems, the NS equations can be reduced to the Euler equations, which are the highest level approximation of inviscid flows. The Euler equations are valid for modelling compressible high speed flows outside of the boundary layer without separation. As most of the industrial process under the interest, as well as significant number of flow situations encountered in nature, are dominated by convective effects, and therefore, they are well approximated by the Euler equations as it will be seen in the validation also (Manna, 1992). Furthermore, the number of equations and the desired cell number are significantly reduced compared with NS based solver, hence it is more suitable for applying them in the optimization methods.

Although the mathematical characteristics of the presented numerical method for the Euler equations are investigated in many articles (e.g. Barth et al., 2004), a measurement of the

Application of Finite Volume Method in Fluid Dynamics and Inverse Design Based Optimization 21

Fig. 9. Lift coefficient in the function of angle of attack at the profile NACA 65-410

*stat*,1 *p* are the dynamic and static pressure respectively at wall surface position 1.

*P*

dispersions near to the leading edge due to the geometrical inaccuracy (sharp edge).

numerical tool is reasonable, it can be used for further applications.

.1 is the solidity, 7.12 *cmc* is the chord and the *g* is the tangential spacing. The

direction and chord (angle of attack) respectively in following. The boundary conditions are set to provide the same Reynolds number as in the experiment. The total inlet pressure is *ptot,in*=101750 [Pa], the total inlet temperature is *Ttot,in*=293.15 [K] and the static outlet pressure is *pstat,out*=101325 [Pa] over the H-type mesh (110×60). The pressure coefficients (38) along the profile are considered in the validation, where *ptot* is the upstream total pressure, 1 *q* and

represents the angle between flow direction and the rotation of axis and the flow

,1

(38)

1

C *tot stat*

*p p q*

The measured parameters and the results of the simulation are compared with each other and the quantitative parameters of the pressure coefficients are shown in Fig. 10. The investigated variables of the calculation, at different angle of attack, are in a good correlation with the measurements (Emery et al., 1958). The difference between them is under the limit of the acceptance, the overall deviation is less then 8 %. The *Cp* values show higher

Although the mathematical aspects of the applied methods as consistency, stability and convergence characteristics are strongly investigated and published in many articles (e.g. Barth et al., 2004), the validation of the Euler based CFD solver was completed in the present subchapter. The results of the analysis and measurements are compared with each other, for a 2D NACA 65-410 wing profile and its low speed cascade. The resembling shows acceptable agreement in engineering point of view. The average deviation between the real tests and the analyses is less then 8 % in both investigated cases, the accuracy of the

where 5

1 and

wing profile NACA 65-410 has been used in the first case to check the accuracy of the calculation. The measurements are performed by Abbott et al. in 1945 and they include the experimental analysis of the lift, drag, pitching moment, etc. of the NACA 6 series airfoils.

Fig. 8. Mach number and pressure distribution over the profile NACA 65-410 at α=4 degrees angle of attack

Most of the data on airfoil section characteristics were obtained in the Langley twodimensional low-turbulence pressure tunnel with a rectangular test section (0.9144 meters wide and 2.286 meters high), in which usually 0.6 meters chord models were tested. The test models completely span the width of the tunnel has a maximum speed of about 70 m/s. More information about the experiments is found in (Abbott, 1945). The Mach number and pressure distribution around the profile is found in Fig. 8. as a result of the computation at α=4 degrees angle of attack (angle between the up stream flow and chord). The boundary conditions are the following: inlet total pressure: *ptot,in*=112800 [Pa]; inlet total temperature: *Ttot,in*=293.15 [K]; outlet static pressure: *pstat,out*=101325 [Pa]. The mesh size is 87×120. The result of the analysis and measurements are compared with each other and they are shown in Fig. 9. in the plot of the lift coefficient in the function of angle of attack. The results of the analysis are accepted in engineering point of view, the overall deviation is less then 5 % at the investigated range. However, it must be considered that the results depend on the mesh resolution and the mesh sensitivity analyses are indispensable to have.

The second test case for the validation is a compressor cascade analysis2 based on the technical report by (Emery et al., 1958). The cascade has been constructed by using NACA 65-410 profile also with:

$$
\sigma = \frac{c}{\mathcal{S}} \tag{37}
$$

2 The original and full version of the present investigation is found in (Veress et al., 2010)

wing profile NACA 65-410 has been used in the first case to check the accuracy of the calculation. The measurements are performed by Abbott et al. in 1945 and they include the experimental analysis of the lift, drag, pitching moment, etc. of the NACA 6 series airfoils.

angle of attack

65-410 profile also with:

Fig. 8. Mach number and pressure distribution over the profile NACA 65-410 at α=4 degrees

Most of the data on airfoil section characteristics were obtained in the Langley twodimensional low-turbulence pressure tunnel with a rectangular test section (0.9144 meters wide and 2.286 meters high), in which usually 0.6 meters chord models were tested. The test models completely span the width of the tunnel has a maximum speed of about 70 m/s. More information about the experiments is found in (Abbott, 1945). The Mach number and pressure distribution around the profile is found in Fig. 8. as a result of the computation at α=4 degrees angle of attack (angle between the up stream flow and chord). The boundary conditions are the following: inlet total pressure: *ptot,in*=112800 [Pa]; inlet total temperature: *Ttot,in*=293.15 [K]; outlet static pressure: *pstat,out*=101325 [Pa]. The mesh size is 87×120. The result of the analysis and measurements are compared with each other and they are shown in Fig. 9. in the plot of the lift coefficient in the function of angle of attack. The results of the analysis are accepted in engineering point of view, the overall deviation is less then 5 % at the investigated range. However, it must be considered that the results depend on the mesh

The second test case for the validation is a compressor cascade analysis2 based on the technical report by (Emery et al., 1958). The cascade has been constructed by using NACA

> *g c*

(37)

resolution and the mesh sensitivity analyses are indispensable to have.

2 The original and full version of the present investigation is found in (Veress et al., 2010)

Fig. 9. Lift coefficient in the function of angle of attack at the profile NACA 65-410

where 5 .1 is the solidity, 7.12 *cmc* is the chord and the *g* is the tangential spacing. The 1 and represents the angle between flow direction and the rotation of axis and the flow direction and chord (angle of attack) respectively in following. The boundary conditions are set to provide the same Reynolds number as in the experiment. The total inlet pressure is *ptot,in*=101750 [Pa], the total inlet temperature is *Ttot,in*=293.15 [K] and the static outlet pressure is *pstat,out*=101325 [Pa] over the H-type mesh (110×60). The pressure coefficients (38) along the profile are considered in the validation, where *ptot* is the upstream total pressure, 1 *q* and *stat*,1 *p* are the dynamic and static pressure respectively at wall surface position 1.

$$\mathbf{C}\_P = \frac{p\_{\text{tot}} - p\_{\text{stat},1}}{q\_1} \tag{38}$$

The measured parameters and the results of the simulation are compared with each other and the quantitative parameters of the pressure coefficients are shown in Fig. 10. The investigated variables of the calculation, at different angle of attack, are in a good correlation with the measurements (Emery et al., 1958). The difference between them is under the limit of the acceptance, the overall deviation is less then 8 %. The *Cp* values show higher dispersions near to the leading edge due to the geometrical inaccuracy (sharp edge).

Although the mathematical aspects of the applied methods as consistency, stability and convergence characteristics are strongly investigated and published in many articles (e.g. Barth et al., 2004), the validation of the Euler based CFD solver was completed in the present subchapter. The results of the analysis and measurements are compared with each other, for a 2D NACA 65-410 wing profile and its low speed cascade. The resembling shows acceptable agreement in engineering point of view. The average deviation between the real tests and the analyses is less then 8 % in both investigated cases, the accuracy of the numerical tool is reasonable, it can be used for further applications.

Application of Finite Volume Method in Fluid Dynamics and Inverse Design Based Optimization 23

Today, beside the developments of the central core of the fluid dynamics solvers, the different optimization techniques, related with CFD, are also under intensive research. In case of direct optimization techniques, an attempt has been made to find the optimal solution. They typically utilize some sort of search technique (e.g. gradient-based optimizer), stochastic based algorithms (e.g. evolutionary strategies, genetic algorithms), artificial neural networks or some other optimization methods. These procedures can be computationally expensive because several flow solutions must be calculated to specify the direction of deepest descent, fitness of individuals in the population, etc. in order to determine the shape changes (Lane, 2010). Furthermore, the required number of flow

In case of a specific set of the inverse design-type methods, the geometry modification is based on the prescribed set of the pre-defined variables at the wall by simple, fast and robust algorithms, which makes them especially attractive amongst other optimization techniques (Lane, 2010). The wall modification can be completed within much less flow solutions for inverse design techniques than for direct optimization methods. Hence, the inverse design methods typically being much more computationally efficient and they are very innovative to be used in practice. The main drawback of inverse design methods is that the designer should create target (optimum in a specific sense) pressure or velocity distributions that should correspond to the design goals and meet the required aerodynamic characteristics. However, it can be difficult to specify the expected pressure or velocity distribution that satisfies all design goals. Also, one cannot guarantee that an arbitrarily prescribed pressure/velocity distribution will provide mechanically correct surfaces or bodies (airfoils without trailing edge open or cross over for example). Hence, the one of the main goals of the following subchapters is to

The calculation process of the developed iterative type inverse design method is shown in Fig. 11. The procedure, first of all, requires an initial geometry and a required pressure distribution (*preq*) along the wall to be modified. The prescribed distribution can be the goal function of an optimization method or it can come from the industrial experiences and/or theory. The iterative cycle starts with the direct solution of the inviscid Euler solver. Completing the convergence criteria, if the target conditions are not reached, a new (opening) boundary condition is applied at the solid boundary to be redesigned or optimized. The required pressure distribution (*preq*) is imposed at the solid wall boundary, which is become locally opening as inlet or outlet, depends upon the evolved pressure differences between the boundary and computational domain. The outcome of this analysis is a velocity distribution along the wall, which is not necessarily parallel with it. The final step of the cycle is the wall modification. The wall becomes parallel with the local velocity vector corresponds to a new streamline of the flow field. The mentioned steps are repeated until the target distribution is reached by the direct analysis and so the new geometry is

All the contributions of the above presented procedure has been described in Chapter 2. except for the wall modification algorithm and the determination of the required pressure distribution (*preq* in Fig. 11.), which are the topic of the following paragraph and subchapter

**4. Finite volume method based optimization and inverse design for**

solutions increases dramatically with the number of design variables.

provide solutions for the above mention complaints.

available (Leonard & Van den Braembussche, 1990).

respectively.

**aeronautical applications**

Fig. 10. Pressure coefficient distribution over the NACA 65-410 compressor cascade profile at different (angle of attack) =30o, Re=2.45E5, and vinf, inlet=29 m/s The experimental data are from (Emery et al., 1958)

#### **4. Finite volume method based optimization and inverse design for aeronautical applications**

22 Finite Volume Method – Powerful Means of Engineering Design

Fig. 10. Pressure coefficient distribution over the NACA 65-410 compressor cascade profile

=30o, Re=2.45E5, and vinf, inlet=29 m/s The experimental

cp

cp

cp

0 20 40 60 80 100

Pressure coefficient distribution =90

0 20 40 60 80 100

percent chord

Pressure coefficient distribution =210

> numerical result experimental data

0 20 40 60 80 100

percent chord

percent chord

Pressure coefficient distribution = 10

> numerical results experimental data

 numerical results experimental data

0 20 40 60 80 100

Pressure coefficient distribution =70

> numerical results experimental data

0 20 40 60 80 100

Pressure coefficient distribution =150

> numerical results experimental data

0 20 40 60 80 100

percent chord

(angle of attack)

percent chord

 numerical results experimental data

percent chord

Pressure coefficient distribution = -50

> -0,5 0,0 0,5 1,0 1,5 2,0 2,5 3,0 3,5

at different

data are from (Emery et al., 1958)

cp

cp

cp

Today, beside the developments of the central core of the fluid dynamics solvers, the different optimization techniques, related with CFD, are also under intensive research. In case of direct optimization techniques, an attempt has been made to find the optimal solution. They typically utilize some sort of search technique (e.g. gradient-based optimizer), stochastic based algorithms (e.g. evolutionary strategies, genetic algorithms), artificial neural networks or some other optimization methods. These procedures can be computationally expensive because several flow solutions must be calculated to specify the direction of deepest descent, fitness of individuals in the population, etc. in order to determine the shape changes (Lane, 2010). Furthermore, the required number of flow solutions increases dramatically with the number of design variables.

In case of a specific set of the inverse design-type methods, the geometry modification is based on the prescribed set of the pre-defined variables at the wall by simple, fast and robust algorithms, which makes them especially attractive amongst other optimization techniques (Lane, 2010). The wall modification can be completed within much less flow solutions for inverse design techniques than for direct optimization methods. Hence, the inverse design methods typically being much more computationally efficient and they are very innovative to be used in practice. The main drawback of inverse design methods is that the designer should create target (optimum in a specific sense) pressure or velocity distributions that should correspond to the design goals and meet the required aerodynamic characteristics. However, it can be difficult to specify the expected pressure or velocity distribution that satisfies all design goals. Also, one cannot guarantee that an arbitrarily prescribed pressure/velocity distribution will provide mechanically correct surfaces or bodies (airfoils without trailing edge open or cross over for example). Hence, the one of the main goals of the following subchapters is to provide solutions for the above mention complaints.

The calculation process of the developed iterative type inverse design method is shown in Fig. 11. The procedure, first of all, requires an initial geometry and a required pressure distribution (*preq*) along the wall to be modified. The prescribed distribution can be the goal function of an optimization method or it can come from the industrial experiences and/or theory. The iterative cycle starts with the direct solution of the inviscid Euler solver. Completing the convergence criteria, if the target conditions are not reached, a new (opening) boundary condition is applied at the solid boundary to be redesigned or optimized. The required pressure distribution (*preq*) is imposed at the solid wall boundary, which is become locally opening as inlet or outlet, depends upon the evolved pressure differences between the boundary and computational domain. The outcome of this analysis is a velocity distribution along the wall, which is not necessarily parallel with it. The final step of the cycle is the wall modification. The wall becomes parallel with the local velocity vector corresponds to a new streamline of the flow field. The mentioned steps are repeated until the target distribution is reached by the direct analysis and so the new geometry is available (Leonard & Van den Braembussche, 1990).

All the contributions of the above presented procedure has been described in Chapter 2. except for the wall modification algorithm and the determination of the required pressure distribution (*preq* in Fig. 11.), which are the topic of the following paragraph and subchapter respectively.

Application of Finite Volume Method in Fluid Dynamics and Inverse Design Based Optimization 25

In the following two subchapters, two case studies of the application of the inverse design method have been presented. In the first one, the lift force is a goal function of an optimization procedure over the NACA 65-410 profile in external flow, meanwhile the increased static pressure ratio is the target condition of redesigning NACA 65-410 cascade geometry in the second case. The pressure distribution should be as low as possible over the solid surface of the suction side at given operational conditions for maximal profile loading. However, in order to reach the downstream conditions, the pressure must increase after the location of the maximum velocity. Stratford's limiting flow theory is used and coupled with the SQP (Sequential Quadratic Programming) nonlinear constraint optimization to provide the target pressure distribution represents the maximum lift force close but certain distance far from the separation in case of external flow. Stratford's limiting flow theory is implemented also in case of redesigning cascade geometry. The presented inverse design method is used to complete wall surface modification till the previously defined target pressure distributions are reached by means of the corresponding sequence of the inverse, wall modification and direct algorithms. The Euler equations are used for modelling basic physics for both cases. The standard cell centred finite volume method has been applied with Roe's approximated Riemann solver, MUSCL approach and Mulder limiter, which are

described in Chapter 2. The validation of the Euler solver is found in Subchapter 3.2.

**4.1 Airfoil optimization for maximal lift force by means of inviscid inverse design**

It has been pointed out in the introduction of the Chapter 4. that the inverse design methods require optimal pressure or velocity distributions to determine the adherent geometry. In order to maximise the lift force of the suction side of a profile at given and constant operational (boundary) conditions, the pressure distribution should be minimized. However, the adverse pressure gradient is appeared after the location of the maximum velocity (and minimum pressure) in order to recover downstream conditions. The adverse pressure gradient till the trailing edge should have limited in each discretized points to be just below the condition of causing separation. The maximum area bounded by the suction and the pressure side distributions in conjunction with the mentioned limited values of pressure gradients will provide the optimum solution as a target distribution to be specified

There are several existing methods for predicting separation as Goldschmied, Stratford, Head, and Cebeci-Smith for example (Smith, 1975). The accuracy these methods were examined several times. One of the output of these investigation shows that the operation of Goldschmied's method is unreliable. The other three are in reasonable agreement and Stratford's method tended to predict separation slightly early. The Cebeci-Smith method is appeared to be best and the Head method is a strong second one (Smith, 1975). Due to the good accuracy, simple expressions and conservative characteristics for predicting

Stratford has derived an empirical formula for predicting the point of separation in an

separation, Stratford's method has been used in followings (Veress et al., 2011).

arbitrary decelerating flow at the order of Re=10E6 (Stratford, 1959),

3 The original and full version of the subchapter is found in (Veress et al., 2011).

**method**3

for the inverse design method.

Fig. 11. Flowchart of the computational procedure of the iterative inverse design calculation

While the incoming and out coming velocity distribution (see *Vn* in Fig. 11.) is given at the solid wall, based on the inverse mode of the inviscid solver (see Subchapter 2.2 at opening wall boundary), the last step of the iterative design cycle is the modification of the geometry. The new position of the solid boundary coordinates is calculated by setting the wall to be parallel with the local velocity vector of the cell centre:

$$\Delta y\_i(\mathbf{x}\_i) = \sum\_{k=l\omega}^{l} \left(\frac{\upsilon\_k}{\mu\_k} \Delta \mathbf{x}\_k\right) \tag{39}$$

where *u* and *v* are the Cartesian component of the velocity vector at the wall. The geometry modification starts from the leading edge or inlet stagnation point till the trailing edge or the outlet stagnation point and completed in vertical directions (see Fig. 12.).

Fig. 12. Schematic view of the wall modification process based on the local velocity vector

Fig. 11. Flowchart of the computational procedure of the iterative inverse design calculation While the incoming and out coming velocity distribution (see *Vn* in Fig. 11.) is given at the solid wall, based on the inverse mode of the inviscid solver (see Subchapter 2.2 at opening wall boundary), the last step of the iterative design cycle is the modification of the geometry. The new position of the solid boundary coordinates is calculated by setting the wall to be

(*xy* )

*Lek*

where *u* and *v* are the Cartesian component of the velocity vector at the wall. The geometry modification starts from the leading edge or inlet stagnation point till the trailing

edge or the outlet stagnation point and completed in vertical directions (see Fig. 12.).

Fig. 12. Schematic view of the wall modification process based on the local velocity vector

 *i*

(39)

*k k k ii x u v*

parallel with the local velocity vector of the cell centre:

In the following two subchapters, two case studies of the application of the inverse design method have been presented. In the first one, the lift force is a goal function of an optimization procedure over the NACA 65-410 profile in external flow, meanwhile the increased static pressure ratio is the target condition of redesigning NACA 65-410 cascade geometry in the second case. The pressure distribution should be as low as possible over the solid surface of the suction side at given operational conditions for maximal profile loading. However, in order to reach the downstream conditions, the pressure must increase after the location of the maximum velocity. Stratford's limiting flow theory is used and coupled with the SQP (Sequential Quadratic Programming) nonlinear constraint optimization to provide the target pressure distribution represents the maximum lift force close but certain distance far from the separation in case of external flow. Stratford's limiting flow theory is implemented also in case of redesigning cascade geometry. The presented inverse design method is used to complete wall surface modification till the previously defined target pressure distributions are reached by means of the corresponding sequence of the inverse, wall modification and direct algorithms. The Euler equations are used for modelling basic physics for both cases. The standard cell centred finite volume method has been applied with Roe's approximated Riemann solver, MUSCL approach and Mulder limiter, which are described in Chapter 2. The validation of the Euler solver is found in Subchapter 3.2.

#### **4.1 Airfoil optimization for maximal lift force by means of inviscid inverse design method**3

It has been pointed out in the introduction of the Chapter 4. that the inverse design methods require optimal pressure or velocity distributions to determine the adherent geometry. In order to maximise the lift force of the suction side of a profile at given and constant operational (boundary) conditions, the pressure distribution should be minimized. However, the adverse pressure gradient is appeared after the location of the maximum velocity (and minimum pressure) in order to recover downstream conditions. The adverse pressure gradient till the trailing edge should have limited in each discretized points to be just below the condition of causing separation. The maximum area bounded by the suction and the pressure side distributions in conjunction with the mentioned limited values of pressure gradients will provide the optimum solution as a target distribution to be specified for the inverse design method.

There are several existing methods for predicting separation as Goldschmied, Stratford, Head, and Cebeci-Smith for example (Smith, 1975). The accuracy these methods were examined several times. One of the output of these investigation shows that the operation of Goldschmied's method is unreliable. The other three are in reasonable agreement and Stratford's method tended to predict separation slightly early. The Cebeci-Smith method is appeared to be best and the Head method is a strong second one (Smith, 1975). Due to the good accuracy, simple expressions and conservative characteristics for predicting separation, Stratford's method has been used in followings (Veress et al., 2011).

Stratford has derived an empirical formula for predicting the point of separation in an arbitrary decelerating flow at the order of Re=10E6 (Stratford, 1959),

 3 The original and full version of the subchapter is found in (Veress et al., 2011).

Application of Finite Volume Method in Fluid Dynamics and Inverse Design Based Optimization 27

measured from the very start of the flow, which begins as flat-plate, turbulent flow. The number *n* is a constant that Stratford finds to be about 6. The quantities *a* and *b* are arbitrary constants used in matching values and slopes in the two equations at the joining point, *p nnC* 12 *.* Of course, eq. (42) describes the beginning of the flow and eq. (43) the final part. The flow is an equilibrium flow that always has the same margin, if any, against separation. Two families of such flows have been computed; they are shown in Fig. 13. The features of the presented diagram, together with eq. (42) and (43) are found in (Smith, 1975). The method presented above is included in determining the pressure distribution at maximum lift force and at the limit of separation on the suction side for given far field

*p p p dxxCxCxC* 2 2 5.0

*dx*

where *p* is the static pressure at the given wall location and the other primitive variables correspond to free stream condition denoted by ∞ (downstream conditions are used in case of cascade design). The connection between *xC p* and *xC p* (pressure coefficient) is given by:

*u*

*pp*

2

*u*

*p xC*

*init*

0

*pp*

2 1

2

*u*

0

*pp*

2 1 1

subject to 0 *TE*

The reason of the constraint to be specified at the presented way is to fix trailing edge (*TE*) condition of Stratford's method and to minimize the disturbances of the optimal pressure distribution of the computational domain respect to the initial flow field. The posterior

The optimization procedure is divided by two sub steps. In the first sub step the physical connections between different parameters are described by Stratford's criteria to evaluate limiting pressure distribution. The pressure coefficient at the minimum pressure (*p*0) is given

0

0 50 70

*Mp.*

numerical test shows that the latter condition is not required, it can differ from zero.

2

*u*

0

*pp CC ,pp*

2

*TE*

2

*u*

*pp*

2 1

*pp C total*

0, 0,

*p pp*

*C CC*

2 00 0

*u*

minimize

2 1

*p*

The objective function is to

by:

1

2 1 2 00 0

*u*

*pp*

*Mp*

0 0

*pp pp*

2 1

1 (46)

*opt pp* (47)

2

(48)

0

*pp*

2

1

*Mp. pp*

*pp dx*

, *x* is the distance

(44)

(45)

In that two-part solution, *x*0 is the start of pressure rise, Re *xu* 000

conditions:

$$\frac{\overline{\mathbf{C}}\_{\nu} \left[ \mathbf{x} \left( d \overline{\mathbf{C}}\_{\nu} / d \mathbf{x} \right) \right]^{\mu 2}}{\left( 10^{-6} \operatorname{Re} \right)^{\mu 10}} = \mathbf{S} \tag{40}$$

where the canonical pressure distribution is

$$\overline{C}\_{\rho} = \frac{p - p\_o}{\frac{1}{2}\rho\_o u\_o^2} \tag{41}$$

and if 0 22 *d dxp* then *S* = 0.39 or if 0 22 *d dxp* then *S* =0.35. Additionally, *Cp* 7/4 . The flows under investigations consist first of a flat-plate flow. Hence, *x* is distance measured from the leading edge of the plate and Re 0*xu .* If the flows begin the pressure rise at a point *x*0 (it is the position of minimum pressure, *p*0 and maximum velocity, *u*0 )*,* left-hand side of eq. (40) starts from a zero value. The left-hand side then grows. When it reaches the limiting value of *S*, separation is said to occur. If *S* is held at its limiting value of 0.39 for 0 22 *d dxp* eq. (40) amounts to an ordinary differential equation for *xC* )(*p .* It is evident from eq. (40) that the equation describes a flow that is ready everywhere to separate. Stratford presents the following solutions (Stratford, 1959),

$$
\overline{\mathbf{C}}\_{\rho} = 0.645 \left[ 0.435 \operatorname{Re}\_{0}^{\prime \mu \tau} \left[ \left( \mathbf{x} / \mathbf{x}\_{0} \right)^{\mu \tau} - 1 \right] \right]^{2\mu} \quad \text{for} \quad \overline{\mathbf{C}}\_{\rho} \le (n-2) / (n+1) \tag{42}
$$

and

$$\overline{C}\_{\text{p}} = 1 - \frac{a}{[\overline{\{\text{x}/\text{x}\_{0}\} + b}]^{\text{p}2}} \quad \text{for} \quad \overline{C}\_{\text{p}} \ge (n-2)/(n+1) \tag{43}$$

Fig. 13. Stratford limiting flows at two values of unit Reynolds number (Smith, 1975)

In that two-part solution, *x*0 is the start of pressure rise, Re *xu* 000 , *x* is the distance measured from the very start of the flow, which begins as flat-plate, turbulent flow. The number *n* is a constant that Stratford finds to be about 6. The quantities *a* and *b* are arbitrary constants used in matching values and slopes in the two equations at the joining point, *p nnC* 12 *.* Of course, eq. (42) describes the beginning of the flow and eq. (43) the final part. The flow is an equilibrium flow that always has the same margin, if any, against separation. Two families of such flows have been computed; they are shown in Fig. 13. The features of the presented diagram, together with eq. (42) and (43) are found in (Smith, 1975).

The method presented above is included in determining the pressure distribution at maximum lift force and at the limit of separation on the suction side for given far field conditions:

$$\mathfrak{D}\left(\overline{\mathbb{C}}\_{\rho}(\mathbf{x})\right) = \mathfrak{D}\left(\mathbb{C}\_{\rho}(\mathbf{x})\right) = \oint \mathbb{C}\_{\rho}(\mathbf{x})d\mathbf{x} = \oint \frac{p - p\_{\circ \circ}}{1}d\mathbf{x} = \oint \frac{p - p\_{\circ \circ}}{0.5\eta p\_{\circ \circ}M\_{\circ \circ}^{2}}d\mathbf{x} \tag{44}$$

where *p* is the static pressure at the given wall location and the other primitive variables correspond to free stream condition denoted by ∞ (downstream conditions are used in case of cascade design). The connection between *xC p* and *xC p* (pressure coefficient) is given by:

$$\overline{C}\_{p} = \frac{p - p\_o}{\frac{1}{2}\rho\_o u\_o^2} = \frac{C\_p - C\_{p,o}}{1 - C\_{p,o}} = \frac{\frac{p - p\_o}{2} - \frac{p\_o - p\_o}{2}}{1 - \frac{p\_o - p\_o}{2}} = \frac{\frac{p}{2}\rho\_o u\_o^2}{p^{sat} - p\_o} = \frac{p - p\_o}{\frac{1}{2}\rho\_o u\_o^2} \tag{45}$$

The objective function is to

26 Finite Volume Method – Powerful Means of Engineering Design

21

2 00 0

Re435.0645.0 1 12 2 51

1 12 21

*nnCfor bxx*

*p* (42)

*a Cp p* (43)

eq. (43)

0 51 *C* 0 *xx p nnCfor n*

Fig. 13. Stratford limiting flows at two values of unit Reynolds number (Smith, 1975)

0

eq. (42)

Re10 (40)

(41)

*.* If the flows begin the pressure rise at a

 *S p p dxCdxC* 101 6

> 2 1 *u*

and if 0 22 *d dxp* then *S* = 0.39 or if 0 22 *d dxp* then *S* =0.35. Additionally, *Cp* 7/4 . The flows under investigations consist first of a flat-plate flow. Hence, *x* is distance measured

point *x*0 (it is the position of minimum pressure, *p*0 and maximum velocity, *u*0 )*,* left-hand side of eq. (40) starts from a zero value. The left-hand side then grows. When it reaches the limiting value of *S*, separation is said to occur. If *S* is held at its limiting value of 0.39 for 0 22 *d dxp* eq. (40) amounts to an ordinary differential equation for *xC* )(*p .* It is evident from eq. (40) that the equation describes a flow that is ready everywhere to separate.

*pp Cp*

where the canonical pressure distribution is

from the leading edge of the plate and Re 0*xu*

and

Stratford presents the following solutions (Stratford, 1959),

$$\text{minimize} \quad \overbrace{\text{N}[\overline{\mathbb{C}}\_r(\mathbf{x})]}^{1} \tag{46}$$

$$\mathbf{a}\mathbf{subject}\text{ to }p\_{\text{opt}}^{\text{TE}} - p\_{\text{init}}^{\text{TE}} = \mathbf{0} \tag{47}$$

The reason of the constraint to be specified at the presented way is to fix trailing edge (*TE*) condition of Stratford's method and to minimize the disturbances of the optimal pressure distribution of the computational domain respect to the initial flow field. The posterior numerical test shows that the latter condition is not required, it can differ from zero.

The optimization procedure is divided by two sub steps. In the first sub step the physical connections between different parameters are described by Stratford's criteria to evaluate limiting pressure distribution. The pressure coefficient at the minimum pressure (*p*0) is given by:

$$\mathbf{C}\_{p} = \mathbf{C}\_{p,o} = \frac{p\_o - p\_{\circ}}{\frac{1}{2}\rho\_{\circ}\mu\_{\circ}^2} = \frac{p\_o - p\_{\circ}}{0.5\eta p\_{\circ}M\_{\circ}^2} = \frac{p\_o - p\_{\circ}}{0.7p\_{\circ}M\_{\circ}^2} \tag{48}$$

Application of Finite Volume Method in Fluid Dynamics and Inverse Design Based Optimization 29

(optimum) one. The optimum pressure distribution belongs to the maximum area of the closed distribution of the pressure and suction side at the limit of separation in case of adverse pressure gradient flow conditions on the suction side. However, several points near to the leading edge of the suction side are modified to make the extremely high pressure gradient smoother. Moreover, an arbitrary (optimal) target pressure distribution often causes nonrealistic geometry as negative thickness, trailing edge opening or cross over. Based on several theoretical investigation and computational tests, it can be noticed, that the expected pressure distribution can not be arbitrary in case of subsonic flow due to the information propagation into the upstream (leading edge) direction along the streamline bounded by the wall. If the required pressure is differ from the initial one at the certain representative part of the near wall region, the flow can be retarded or sucked depends on the local conditions. This effect has an influence on the flow evolution starting from the leading edge and the pressure should be redistributed by considering higher or lower local kinetic energy along the stream line

Fig. 14. Pressure distribution of the initial (init), optimum (target) and result (of the inverse

The modified distributions have been imposed in the inverse design procedure to determine the geometry, which provides the optimal conditions. The inverse design method was converged after 10 iteration cycles of the inverse, wall modification and direct modes. The normal velocity distribution across the solid wall becomes near to zero at the last inverse subroutine, which represents that there is no need for any further steps, the pressure gradient is infinitesimally small (no flow) across the solid boundary. The corresponding results of the optimization procedure are found in Fig. 14. The target and optimized (result) pressure distribution are compared with each other and the deviation between them is negligible. The optimized geometry with Mach number and pressure distribution is shown in Fig. 15. The improvements are straightforward; the lower pressures at the suction side provide higher lift force in case of the optimum geometry. Further quantitative results are found in Fig. 16. The optimization was completed at zero angle of attack. However, the offdesign conditions show also the same order of improvements as the optimization at design

design based optimization procedure) cases (ps: pressure side and ss: suction side)

especially at the first couple mesh points of the leading edge.

where *p*0 and maximum velocity *u*0 is supposed to be constant starting from the leading edge of the suction side till the starting of the positive pressure gradient ( *x*0 ). The Mach numbers *M*0 at these points are calculated by:

$$\mathbf{C}\_{p} = \mathbf{C}\_{p,o} = \frac{p\_o - p\_{\circ \circ}}{0.7 p\_{\circ \circ} M\_{\circ}^2} = \frac{1}{0.7 M\_{\circ}^2} \left[ \left( \frac{1 + 0.2 M\_{\circ}^2}{1 + 0.2 M\_{\circ}^2} \right)^{\frac{r}{r-1}} - 1 \right] \tag{49}$$

The *T*0 , *u*0 and 0 are obtained by the energy equation of the isentropic flow and ideal gas law:

$$T\_o = T^{\text{total}} \left( 1 + \frac{\gamma - 1}{2} M\_o^2 \right)^{-1} \tag{50}$$

$$
\mu\_0 = \sqrt{\frac{\mathcal{Y}}{\mathcal{Y} - \mathbf{1}} \mathcal{R} \left( T^{\text{total}} - T\_0 \right) \mathbf{2}} \tag{51}
$$

$$
\rho\_0 = \frac{p\_0}{RT\_0} \tag{52}
$$

The *total* quantities correspond to the given operational (flight) or inlet boundary conditions.

A general way of determining pressure distribution starts with specifying a possible *p*0 . All parameter belongs to *p*0 can be calculated by eqs. (48)-(52). The next step is to find location *x*0 , which gives back the required trailing edge static pressure by using Stratford's equations (42) and (43) over *x*. Hence, the location of starting flow deceleration ( *x*0 ) and the Stratford's limiting pressure distribution till the required trailing edge pressure is the output of the first sub step of the optimization procedure. There are infinite possible pressure distribution existing of the presented method hence, the second sub step of the optimization procedure is the constraint optimization in order to determine the corresponding flow parameters and location belongs to the minimum pressure and maximum velocity point on the suction surface, which provides the maximum area bounded by the pressure distribution of the suction and pressure side of the profile. *p*0 , *T*0 , *u*0 , 0 , *x*0 and *xp* )( (by Stratford's criteria) parameters will be modified in the second sub step to satisfy (46) and (47).

The pressure side distribution is also modified by means of constraint optimization to maximize the area under the function restricted to less than or equal to the maximum pressure gradient or higher than or equal to the minimum pressure gradient respect to the original distribution.

NACA 65-410 profile has been used to provide initial geometry and flow field for the optimization. The boundary conditions are the followings: inlet total pressure: *ptot,in*=112800 [Pa]; inlet total temperature: *Ttot,in*=293.15 [K]; outlet static pressure: *pstat,out*=101325 [Pa] over the mesh size of 87×120. The pressure distributions of the given geometry are shown in Fig. 14. and they are noted as init (initial). The inverse design program modifies initial geometry till the result pressure distribution over the geometry gives back exactly the target

where *p*0 and maximum velocity *u*0 is supposed to be constant starting from the leading edge of the suction side till the starting of the positive pressure gradient ( *x*0 ). The Mach

2 2

70 1

0 0 2 1 1

1 *u*0 *TTR* 0 *total*

 *TT total M*

0 0 0 *RT p*

The *total* quantities correspond to the given operational (flight) or inlet boundary conditions.

A general way of determining pressure distribution starts with specifying a possible *p*0 . All parameter belongs to *p*0 can be calculated by eqs. (48)-(52). The next step is to find location *x*0 , which gives back the required trailing edge static pressure by using Stratford's equations (42) and (43) over *x*. Hence, the location of starting flow deceleration ( *x*0 ) and the Stratford's limiting pressure distribution till the required trailing edge pressure is the output of the first sub step of the optimization procedure. There are infinite possible pressure distribution existing of the presented method hence, the second sub step of the optimization procedure is the constraint optimization in order to determine the corresponding flow parameters and location belongs to the minimum pressure and maximum velocity point on the suction surface, which provides the maximum area bounded by the pressure distribution of the

The pressure side distribution is also modified by means of constraint optimization to maximize the area under the function restricted to less than or equal to the maximum pressure gradient or higher than or equal to the minimum pressure gradient respect to the

NACA 65-410 profile has been used to provide initial geometry and flow field for the optimization. The boundary conditions are the followings: inlet total pressure: *ptot,in*=112800 [Pa]; inlet total temperature: *Ttot,in*=293.15 [K]; outlet static pressure: *pstat,out*=101325 [Pa] over the mesh size of 87×120. The pressure distributions of the given geometry are shown in Fig. 14. and they are noted as init (initial). The inverse design program modifies initial geometry till the result pressure distribution over the geometry gives back exactly the target

2

*M.Mp.*

70

0

suction and pressure side of the profile. *p*0 , *T*0 , *u*0 ,

original distribution.

parameters will be modified in the second sub step to satisfy (46) and (47).

0

 1 201 201

0 are obtained by the energy equation of the isentropic flow and ideal gas

1 2

(52)

1

2

2 0

*M. M.*

*pp CC ,pp* (49)

(50)

(51)

0 , *x*0 and *xp* )( (by Stratford's criteria)

numbers *M*0 at these points are calculated by:

The *T*0 , *u*0 and

law:

(optimum) one. The optimum pressure distribution belongs to the maximum area of the closed distribution of the pressure and suction side at the limit of separation in case of adverse pressure gradient flow conditions on the suction side. However, several points near to the leading edge of the suction side are modified to make the extremely high pressure gradient smoother. Moreover, an arbitrary (optimal) target pressure distribution often causes nonrealistic geometry as negative thickness, trailing edge opening or cross over. Based on several theoretical investigation and computational tests, it can be noticed, that the expected pressure distribution can not be arbitrary in case of subsonic flow due to the information propagation into the upstream (leading edge) direction along the streamline bounded by the wall. If the required pressure is differ from the initial one at the certain representative part of the near wall region, the flow can be retarded or sucked depends on the local conditions. This effect has an influence on the flow evolution starting from the leading edge and the pressure should be redistributed by considering higher or lower local kinetic energy along the stream line especially at the first couple mesh points of the leading edge.

Fig. 14. Pressure distribution of the initial (init), optimum (target) and result (of the inverse design based optimization procedure) cases (ps: pressure side and ss: suction side)

The modified distributions have been imposed in the inverse design procedure to determine the geometry, which provides the optimal conditions. The inverse design method was converged after 10 iteration cycles of the inverse, wall modification and direct modes. The normal velocity distribution across the solid wall becomes near to zero at the last inverse subroutine, which represents that there is no need for any further steps, the pressure gradient is infinitesimally small (no flow) across the solid boundary. The corresponding results of the optimization procedure are found in Fig. 14. The target and optimized (result) pressure distribution are compared with each other and the deviation between them is negligible. The optimized geometry with Mach number and pressure distribution is shown in Fig. 15. The improvements are straightforward; the lower pressures at the suction side provide higher lift force in case of the optimum geometry. Further quantitative results are found in Fig. 16. The optimization was completed at zero angle of attack. However, the offdesign conditions show also the same order of improvements as the optimization at design

Application of Finite Volume Method in Fluid Dynamics and Inverse Design Based Optimization 31

The main goal of the first part of the subchapter is to show the outcome of the profile with maximum blade loading developed by the inverse design method coupling with SQP and Stratford's limiting flow theory (see Subchapter 4.1.). The result of the optimization was analyzed by the present Euler solver and ANSYS CFX. The same computational procedure has been used in both software except for the mesh at viscous mode of the CFX. y+ determines the first cell size next to the wall and has an influence on the cell number, which results 71X23 2D volumes. Concerning the boundary conditions, the following physical parameters were used: the inlet total pressure is ptot,in= 107853 [Pa], the inlet total temperature is Ttot,in= 298.42 [K], the inlet flow angle is 30° and the outlet static pressure is

The CFX was convergent after 110 iterations. The quantitative results of the analyses are shown in Fig. 17. The target and the result pressure distributions of the inviscid inverse design based optimization are presented over the final geometry beside the inviscid (Ansysss and Ansys-ps) and viscous (Ansys-ss-viscous and Ansys-ps-viscous) results of the CFX (ss: suction side and ps: pressure side). The average difference between the three approaches is less than 8 %, which is acceptable in engineering point of view. The main deviation is at the leading edge stagnation point; the static pressure in the Euler solver is higher compared with CFX. This unphysical feature is caused by the linear extrapolation of determining static

pressure at the ghost cell of the solid wall boundaries without averaging procedure.

2.2 2.7 3.2 3.7 4.2

due to the unexpected thick profile, which causes chocking.

**axial length [m]**

Fig. 17. Pressure distribution of the redesigned blade in case of recent Euler solver and the inviscid and viscous analysis of ANSYS CFX (ss: suction side and ps: pressure side)

It can be observed in Fig. 17 that the outlet static pressure is higher than the inlet one, so the cascade is working in a compressor mode, however, the static pressure ratio is negligible

**Pressure Distribution**

target\_ss

target\_ps

result\_ss

result\_ps

Ansys-ss

Ansys-ps

Ansys-ss-viscous

Ansys-ps-viscous

pstat,out=101325 [Pa].

88600

93600

98600

**pressure [Pa]**

103600

108600

point. The lift force coefficient is increased significantly around by 100 % in the investigated range of the angle of attack. The effect of drag force should be analyzed and considered with the aim of including it in the optimization process.

Fig. 15. Mach number and pressure distribution [Pa] of the result case (result of the inverse design based optimization procedure)

Fig. 16. Lift coefficient distribution in the function of angle of attack (the optimization has been completed at angle of attack 0)

#### **4.2 Application of the inverse design based semi-optimization in cascade flows**

Compressors and turbines are widely used in the vehicle engines as gas turbines for shaft power and jet engines. The axial compressors and turbines can be simplified to cascade geometry by extracting a cylindrical cut surface from them and laying out in 2D. The basic characteristics of the elementary flow field can be analysed and studied by this way.

NACA 65-410 profile has been used also in followings. The validation of the Euler solver for the cascade geometry is found in Subchapter 3.2.

point. The lift force coefficient is increased significantly around by 100 % in the investigated range of the angle of attack. The effect of drag force should be analyzed and considered with

 Fig. 15. Mach number and pressure distribution [Pa] of the result case (result of the inverse

Fig. 16. Lift coefficient distribution in the function of angle of attack (the optimization has

Compressors and turbines are widely used in the vehicle engines as gas turbines for shaft power and jet engines. The axial compressors and turbines can be simplified to cascade geometry by extracting a cylindrical cut surface from them and laying out in 2D. The basic

NACA 65-410 profile has been used also in followings. The validation of the Euler solver for

**4.2 Application of the inverse design based semi-optimization in cascade flows**

characteristics of the elementary flow field can be analysed and studied by this way.

the aim of including it in the optimization process.

design based optimization procedure)

been completed at angle of attack 0)

the cascade geometry is found in Subchapter 3.2.

The main goal of the first part of the subchapter is to show the outcome of the profile with maximum blade loading developed by the inverse design method coupling with SQP and Stratford's limiting flow theory (see Subchapter 4.1.). The result of the optimization was analyzed by the present Euler solver and ANSYS CFX. The same computational procedure has been used in both software except for the mesh at viscous mode of the CFX. y+ determines the first cell size next to the wall and has an influence on the cell number, which results 71X23 2D volumes. Concerning the boundary conditions, the following physical parameters were used: the inlet total pressure is ptot,in= 107853 [Pa], the inlet total temperature is Ttot,in= 298.42 [K], the inlet flow angle is 30° and the outlet static pressure is pstat,out=101325 [Pa].

The CFX was convergent after 110 iterations. The quantitative results of the analyses are shown in Fig. 17. The target and the result pressure distributions of the inviscid inverse design based optimization are presented over the final geometry beside the inviscid (Ansysss and Ansys-ps) and viscous (Ansys-ss-viscous and Ansys-ps-viscous) results of the CFX (ss: suction side and ps: pressure side). The average difference between the three approaches is less than 8 %, which is acceptable in engineering point of view. The main deviation is at the leading edge stagnation point; the static pressure in the Euler solver is higher compared with CFX. This unphysical feature is caused by the linear extrapolation of determining static pressure at the ghost cell of the solid wall boundaries without averaging procedure.

Fig. 17. Pressure distribution of the redesigned blade in case of recent Euler solver and the inviscid and viscous analysis of ANSYS CFX (ss: suction side and ps: pressure side)

It can be observed in Fig. 17 that the outlet static pressure is higher than the inlet one, so the cascade is working in a compressor mode, however, the static pressure ratio is negligible due to the unexpected thick profile, which causes chocking.

Application of Finite Volume Method in Fluid Dynamics and Inverse Design Based Optimization 33

Abbott H.; Doenhoff A. E. & Stivers, L. S. (1945). *NACA REPORT No. 824., Summary of Airfoil*

*ANSYS, Inc. (2010).* ANSYS Fluent User's guide*, ANSYS, Inc. Southpointe, 275 Technology Derive Canonsburg, PA 15317<EMAIL_ADDRESS>http://www.ansys.com, USA*  Barth, T. & Ohlberger, M. (2004). *Finite volume methods: foundation and analysis*, In: Stein, E.;

Blazek, J. (2005). *Computational Fluid Dynamics: Principle and Applications*, Elsevier, ISBN-13:

CFD Online Discussion Forum (2010). *http://www.cfd-online.com/Forums/main/75554-use-k-*

Emery, J. C.; Herrig, L. J; Erwin, J. R. & Felix, A. R. (1958). *Systematic two-dimensional cascade*

Gerolymos, G. A.; Sauret, E. & Vallet, I. (2003). *Oblique-Shock-Wave/Boundary-Layer Interaction*

Hirsch, C. (1990, 2007). *Numerical Computation of Internal and External Flows,* Wiley, ISBN 978-

Jones, W. P. & Launder B. E. (1973).*The Calculation of Low-Reynolds-Number-Phenomena with a*

Karki, K. & Patankar, S. V. (1989). *Pressure based calculation procedure for viscous flows at all*

Kuzmin, D. & Möller, M. (2004). *Algebraic Flux Correction II. Compressible Euler Equations*, http://www.mathematik.tu-dortmund.de/papers/KuzminMoeller2004a.pdf,

Lane, K. A. (2010). *Novel Inverse Airfoil Design Utilizing Parametric Equations*, MSC thesis, Faculty of California Polytechnic State University, San Luis Obispo. Lefebvre, M. & Arts, T. (1997). *Numerical aero-thermal prediction of laminar/turbulent flows in a*

Leonard, O. & Van den Braembussche, R. (1990). *Subsonic and Transonic Cascade Design*,

Manna, M. (1992). *A Three Dimensional High Resolution Compressible Flow Solver,* PhD thesis,

Menter, F. R. (1994). *Two-Equation Eddy-Viscosity Turbulence Models for Engineering*

McDonald, P. M. (1971). *The Computation of Transonic Flow through two Dimensional Gas*

and Turbomachinery Applications, May, 14-18, 1990, Belgium

*Applications,* AIAA Journal, vol. 32, no. 8, pp. 1598-1605, USA

Catholic University of Leuven, Belgium

*Turbine Cascades,* ASME paper 71-GT-89, USA

978-0-08-044506-9, ISBN-10: 0-08-044506-3, United Kingdom

Borst, R. & Hughes, T. (eds.) Encyclopedia of Computational Mechanics, John

*tests of NACA 65-series compressor blades at low speeds*, NACA Technical Report 1368,

*using Near-Wall Reynolds-Stress Models*, Université Pierre-et-Marie-Curie, AIAA 2003-3466, 33rd Fluid Dynamics Conference, 23-26 June 2003 Orlando, Florida, USA

*Two-Equation Model of Turbulence*, Int. J. Heat Mass Transl., Vol. 16, 1973, pp. 1119-

*speeds in arbitrary configurations*, AIAA Journal, Vol. 27, No 9, pp. 1167-1174, 1989,

*two-dimensional high pressure turbine linear cascade*, Second European Conference on Turbomachinery - Fluid Dynamics and Thermodynamics, Antwerp, Belgium, pp.

AGARD-VKI Special Course on Inverse Methods in Airfoil Design for Aeronautical

**6. References**

*Data,* USA

USA

1130

USA

Germany

401-409

Wiley & Sons, New York, USA

0-7506-6594-0, United Kingdom

*epsilon-k-omega-models.html*

Of course, the lower outlet static pressure – or higher mass flow – can improve design specification by means of increasing static pressure rise over the cascade due to the energy conversion. Moreover, the higher blade loading can also increase the static pressure ratio, which is shown in the second part of the present subchapter. The blade geometry variants in the function of the blade loading at the same mesh, solver settings, initial and boundary conditions are found in Fig. 18. with the corresponding pressure coefficients and distributions, which are based on the Stratford's limiting flow theory. The pressure side pressure distribution was the same at the three investigated cases. The physical boundary conditions are the followings: the inlet total pressure is ptot,in= 107853.4 [Pa], the inlet total temperature is Ttot,in= 298.4267 [K], the inlet flow angle is 30° and the outlet static pressure is pstat,out= 83325 [Pa]. Although the static pressure ratio is increased from 1.034 to 1.14 proportionally with blade loading and the absolute value of the pressure coefficient, further numerical and experimental investigations are indispensable to have for well established conclusions in the field of cascade optimization.

Fig. 18. Blade geometries (left side) at different Stratford's based pressure distributions (blade loadings) (right side) with the corresponding minimum pressure coefficients (ss: suction side and ps: pressure side)

#### **5. Acknowledgment**

This work has been supported by the Hungarian National Fund for Science and Research (OTKA) under the fund No. F 67555. The results reported in the chapter has been developed in the framework of the project "Talent care and cultivation in the scientific workshops of BME" project. This project is supported by the grant TÁMOP-4.2.2.B-10/1-2010-0009.

#### **6. References**

32 Finite Volume Method – Powerful Means of Engineering Design

Of course, the lower outlet static pressure – or higher mass flow – can improve design specification by means of increasing static pressure rise over the cascade due to the energy conversion. Moreover, the higher blade loading can also increase the static pressure ratio, which is shown in the second part of the present subchapter. The blade geometry variants in the function of the blade loading at the same mesh, solver settings, initial and boundary conditions are found in Fig. 18. with the corresponding pressure coefficients and distributions, which are based on the Stratford's limiting flow theory. The pressure side pressure distribution was the same at the three investigated cases. The physical boundary conditions are the followings: the inlet total pressure is ptot,in= 107853.4 [Pa], the inlet total temperature is Ttot,in= 298.4267 [K], the inlet flow angle is 30° and the outlet static pressure is pstat,out= 83325 [Pa]. Although the static pressure ratio is increased from 1.034 to 1.14 proportionally with blade loading and the absolute value of the pressure coefficient, further numerical and experimental investigations are indispensable to have for well established

Fig. 18. Blade geometries (left side) at different Stratford's based pressure distributions (blade loadings) (right side) with the corresponding minimum pressure coefficients (ss:

This work has been supported by the Hungarian National Fund for Science and Research (OTKA) under the fund No. F 67555. The results reported in the chapter has been developed in the framework of the project "Talent care and cultivation in the scientific workshops of

BME" project. This project is supported by the grant TÁMOP-4.2.2.B-10/1-2010-0009.

conclusions in the field of cascade optimization.

suction side and ps: pressure side)

**5. Acknowledgment**

**1. Introduction**

complete overview of currently used methods.

**2. Mathematical models**

**2.1 Shallow water equations**

general cross-section area described by

Balance laws arise from many areas of engineering practice specifically from the fluid mechanics. Many numerical methods for the solution of these balanced laws were developed in recent decades. The numerical methods are based on two views: solving hyperbolic PDE with a nonzero source term (the obvious description of the central and central-upwind schemes; (Kurganov & Levy, 2002; LeVeque, 2004)) or solving the augmented quasilinear nonconservative formulation (Gosse, 2001; Le Floch & Tzavaras, 1999; Pares, 2006). ` Furthermore, the methods can be interpreted using flux-difference splitting (or flux-vector splitting), or by selecting adaptive intervals and the transformation to the semidiscrete form (for example (Kurganov & Petrova, 2000)). We prefer the augmented quasilinear nonconservative formulation solved by the flux-difference splitting in our text. We try to formulate the methods in the most general form. The range of this text does not give the

**Numerical Schemes for Hyperbolic Balance**

**Laws – Applications to Fluid Flow Problems**

Marek Brandner, JiĢí Egermaier and Hana Kopincová

*The University of West Bohemia*

*Czech Republic*

In this section we describe the specific mathematical models based on hyperbolic balanced laws. There are many models that describe fluid flow phenomena but we are interested in the

We are interested in solving the problem related to the fluid flow through the channel with the

where *a* = *a*(*x*, *t*) is the unknown cross-section area, *q* = *q*(*x*, *t*) is the unknown discharge, *b* = *b*(*x*) is given function of elevation of the bottom, *g* is the gravitational constant and

= −*gabx* + *gI*2,

*at* + *qx* = 0, (1)

**2**

[*h*(*x*, *t*) − *η*]*σ*(*x*, *η*)*dη*, (2)

two type of them: models described open channel flow and urethra flow.

 *q*2 *a* + *gI*1 *x*

 *h*(*x*,*t*) 0

*qt* +

*I*1 =

### **Numerical Schemes for Hyperbolic Balance Laws – Applications to Fluid Flow Problems**

Marek Brandner, JiĢí Egermaier and Hana Kopincová *The University of West Bohemia Czech Republic*

#### **1. Introduction**

34 Finite Volume Method – Powerful Means of Engineering Design

Munz, C.-D.; Roller, S.; Klein, R. & Geratz, K. J. (2003). *The extension of incompressible flow*

Pope, S. B. (2000). *Turbulent Flows*, Cambridge University Press, ISBN<PHONE_NUMBER>866,

Roe, P. L. (1981). *Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes*,

Settles, G. (1975). *The shadowgram of Mach 3 airflow over an 18 degree compression corner,*

Dynamics Lab, Penn State University (PHD Thesis, Princeton, 1975), USA Smith, A. M. O. (1975), *High-Lift Aerodynamics*, Journal of Aircraft*,* Vol. 12 No. 6, pp. 501-530,

Stein, E.; Borst, R. & Hughes, T. (2004). *Finite volume methods: foundation and analysis,*  http://weberknecht.uni-muenster.de/num/publications/2004/BO04a/finvol\_

Stratford, B. S. (1959), *The Prediction of Separation of the Turbulent Boundary Layer*, Journal of

Veress, Á.; Gallina, T. & Rohács, J. (2010). *Fast and Robust Inverse Design Method for Internal*

Veress, Á.; Felföldi, A.; Gausz, T. & Palkovics, L. (2011). *Coupled Problem of the Inverse Design*

Wassgren, C. (2010). *Notes on Fluid Mechanics and Gas Dynamics,* lecture note, School of

Wilcox, D. C. (1998, (1988)). *Turbulence Modelling for CFD,* DCW Industries Inc. Second

Yee, H. C. (1989). *A class of high-resolution explicit and implicit shock-capturing methods*, VKI lecture series 1989-04, March 6-10, 1989; NASA TM-101088, Feb. 1989, Belgium Zucrow, M. J. & Hoffman, J. D. (1976). *Gas Dynamics Vol. 1,* Wiley, ISBN 0-471-98440-X, USA

*and Cascade Flows*, International Review of Aerospace Engineering (IREASE), ISSN

*and Constraint Optimization*, Applied Mathematics and Computation, DOI:

Journal of Computational Physics, Vol. 43 pp. 357-372

script.pdf, Edited by John Wiley & Sons, Ltd., USA

10.1016/j.amc.2011.08.110, Paper In Press, Corrected Proof.

edition, ISBN-10:<PHONE_NUMBER>, ISBN-13:<PHONE_NUMBER>153, USA

Mechanical Engineering, Purdue University, USA

Fluid Mechanics, Vol. 5. pp 1-16, USA

1973-7459 Vol. 3 N. 1. pp. 41-50.

196

USA

United Kingdom

*solvers to the weakly compressible regime*, Computers & Fluids, Vol. 32, No 2, pp. 173-

http://www.efluids.com/efluids/gallery/gallery\_pages/18degramp.htm, Gas

Balance laws arise from many areas of engineering practice specifically from the fluid mechanics. Many numerical methods for the solution of these balanced laws were developed in recent decades. The numerical methods are based on two views: solving hyperbolic PDE with a nonzero source term (the obvious description of the central and central-upwind schemes; (Kurganov & Levy, 2002; LeVeque, 2004)) or solving the augmented quasilinear nonconservative formulation (Gosse, 2001; Le Floch & Tzavaras, 1999; Pares, 2006). ` Furthermore, the methods can be interpreted using flux-difference splitting (or flux-vector splitting), or by selecting adaptive intervals and the transformation to the semidiscrete form (for example (Kurganov & Petrova, 2000)). We prefer the augmented quasilinear nonconservative formulation solved by the flux-difference splitting in our text. We try to formulate the methods in the most general form. The range of this text does not give the complete overview of currently used methods.

#### **2. Mathematical models**

In this section we describe the specific mathematical models based on hyperbolic balanced laws. There are many models that describe fluid flow phenomena but we are interested in the two type of them: models described open channel flow and urethra flow.

#### **2.1 Shallow water equations**

We are interested in solving the problem related to the fluid flow through the channel with the general cross-section area described by

$$a\_l + q\_\chi = 0,\tag{1}$$

$$q\_l + \left(\frac{q^2}{a} + gI\_1\right)\_\chi = -gab\_\chi + gI\_2.$$

where *a* = *a*(*x*, *t*) is the unknown cross-section area, *q* = *q*(*x*, *t*) is the unknown discharge, *b* = *b*(*x*) is given function of elevation of the bottom, *g* is the gravitational constant and

$$I\_1 = \int\_0^{h(\mathbf{x}, t)} [h(\mathbf{x}, t) - \eta] \sigma(\mathbf{x}, \eta) d\eta,\tag{2}$$

**2.2 Urethra flow**

We now briefly introduce a problem describing fluid flow through the elastic tube represented by hyperbolic partial differential equations with the source term. In the case of the male urethra, the system based on model in (Stergiopulos at al., 1993) has the following form

Numerical Schemes for Hyperbolic Balance Laws – Applications to Fluid Flow Problems 37

where *a* = *a*(*x*, *t*) is the unknown cross-section area, *q* = *q*(*x*, *t*) is the unknown flow rate

the cross-section of the tube under no pressure, *β* = *β*(*x*, *t*) is the coefficient describing tube compliance and *λ*(*Re*) is the Mooney-Darcy friction factor (*λ*(*Re*) = 64/*Re* for laminar flow).

where *μ* is fluid viscosity. This model contains constitutive relation between the pressure and

4*a*

*Re* = *ρq μa*

*p* = *a* − *a*0

**3. Conservative and nonconservative problems and numerical schemes**

**q**(*x*, 0) = **q**0(*x*), *x* ∈ **R**, The numerical scheme based on finite volume discretization in the conservation form can be

> *j* − Δ*t* Δ*x* (**F***n*

> > Δ*x*

from time level *tn* to *tn*+1

*xj*+1/2 *xj*−1/2

The relation (14) can be derived as the approximation of the integral conservation law at the

**q**(*x*, *tn*)*dx*−

**f**(**q**(*xj*+1/2, *t*))*dt* −

2*ρβ*2 *βx* − *q*2

4*a*2 *π*

**q***t* + [**f**(**q**)]*x* = **0**, *x* ∈ **R**, *t* ∈ (0, *T*), (13)

*j*+1/2 − **F***n*

*a λ*(*Re*), (10)

*a* ), *ρ* is the fluid density, *a*0 = *a*0(*x*) is

*π* , (11)

*j*−1/2). (14)

**f**(**q**(*xj*−1/2, *t*))*dt*

.

(16)

(**F***j*+1/2 − **F***j*−1/2). (15)

*t n*+1 *tn*

*β* + *pe*, (12)

*at* + *qx* = 0,

*x* = *a ρ a*0 *β x* + *a*2

*qt* + *q*2 *a* + *a*2 2*ρβ*

*Re* is the Reynolds number defined by

the cross section of the tube

written as follows,

interval

1 Δ*x*

where *pe* is surrounding pressure.

(we also denote *v* = *v*(*x*, *t*) as the fluid velocity, *v* = *q*

**3.1 Conservative problems and numerical schemes** We consider the conservation law in the conservative form

> **Q***n*+1 *j* = **Q***n*

> > *∂***Q***j ∂t* = − 1

> > > Δ*x*

− 1 Δ*t t n*+1 *tn*

We use also the semidiscrete version of (14)

**q**(*x*, *tn*+1)*dx* = 1

*xj*−1/2, *xj*+1/2

*xj*+1/2 *xj*−1/2

$$I\_2 = \int\_0^{h(\mathbf{x},t)} (h(\mathbf{x},t) - \eta) \left[ \frac{\partial \sigma(\mathbf{x}, \eta)}{\partial \mathbf{x}} \right] d\eta\_{\prime} \tag{3}$$

here *η* is the depth integration variable, *h*(*x*, *t*) is the water depth and *σ*(*x*, *η*) is the width of the cross-section at the depth *η*. The derivation can be found in e.g. (Cunge at al., 1980).

The first special case are the equations reflecting the fluid flow through the spatially varying rectangular channel

$$a\_l + q\_\chi = 0,\tag{4}$$

$$q\_l + \left(\frac{q^2}{a} + \frac{ga^2}{2l}\right)\_\chi = \frac{ga^2}{2l^2}l\_\chi - gab\_\chi$$

with *l* = *l*(*x*) being the function describing the width of the channel, and second one the system for the constant rectangular channel

$$h\_t + (hu)\_x = 0,\tag{5}$$

$$(hu)\_t + \left(hu^2 + \frac{1}{2}gh^2\right)\_x = -ghb\_x.$$

In the above equation, *h*(*x*, *t*) is the water depth and *u*(*x*, *t*) is the horizontal velocity. It also possible to add some friction term to the system described above. For example, we can write

$$h\_l + (hv)\_x = 0,$$

$$(hv)\_t + \left(hv^2 + \frac{1}{2}gh^2\right)\_x = -ghB\_x - gM^2 \frac{hv|hv|}{h^{7/3}},\tag{6}$$

where *M* is Manning's coefficient.

All of the presented systems can be written in the compact matrix form

$$[\mathbf{q}\_l + [\mathbf{f}(\mathbf{q}\_\nu \mathbf{x})]\_\mathcal{X} = \psi(\mathbf{q}\_\nu \mathbf{x}),\tag{7}$$

with **q**(*x*, *t*) being the vector of conserved quantities, **f**(**q**, *x*) the flux function and ψ(**q**, *x*) the source term. This relation represents the balance laws.

It is possible to use any augmented formulation, which is suitable for rewriting the system to the quasilinear homogeneous one. For example, we can obtain

$$\begin{aligned} h\_t + (hu)\_x &= 0, \\ (hu)\_t + \left(hu^2 + \frac{1}{2}gh^2\right)\_x &= -ghb\_{\mathbf{x}}, \\ b\_t &= 0, \end{aligned} \tag{8}$$

i.e.

$$
\begin{bmatrix} h \\ hv \\ b \end{bmatrix}\_t + \begin{bmatrix} 0 & 1 & 0 \\ -v^2 + gh \ 2v - gh \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} h \\ hv \\ b \end{bmatrix}\_x = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} . \tag{9}
$$

#### **2.2 Urethra flow**

2 Will-be-set-by-IN-TECH

here *η* is the depth integration variable, *h*(*x*, *t*) is the water depth and *σ*(*x*, *η*) is the width of the cross-section at the depth *η*. The derivation can be found in e.g. (Cunge at al., 1980).

The first special case are the equations reflecting the fluid flow through the spatially varying

with *l* = *l*(*x*) being the function describing the width of the channel, and second one the

In the above equation, *h*(*x*, *t*) is the water depth and *u*(*x*, *t*) is the horizontal velocity. It also possible to add some friction term to the system described above. For example, we can write

with **q**(*x*, *t*) being the vector of conserved quantities, **f**(**q**, *x*) the flux function and ψ(**q**, *x*) the

It is possible to use any augmented formulation, which is suitable for rewriting the system to

0 10 −*v*2 + *gh* 2*v* −*gh* 0 00

*ht* + (*hv*)*x* = 0,

�

= *ga*2

*∂σ*(*x*, *η*) *∂x*

2*l*2 *lx* − *gabx*,

= −*ghbx*.

= −*ghBx* − *gM*2 *hv*|*hv*|

= −*ghbx*,

⎤ ⎦ *x* = ⎡ ⎣ 0 0 0

*bt* = 0,

⎤ ⎦ ⎡ ⎣ *h hv b*

**q***t* + [**f**(**q**, *x*)]*x* = ψ(**q**, *x*), (7)

*ht* + (*hu*)*x* = 0, (8)

⎤

⎦ . (9)

�

*at* + *qx* = 0, (4)

*ht* + (*hu*)*x* = 0, (5)

*dη*, (3)

*h*7/3 , (6)

(*h*(*x*, *t*) − *η*)

*I*2 =

*qt* +

(*hu*)*t* +

� *hu*2 + 1 2 *gh*2 � *x*

system for the constant rectangular channel

(*hv*)*t* +

source term. This relation represents the balance laws.

⎡ ⎣ *h hv b*

⎤ ⎦ *t* + ⎡ ⎣

the quasilinear homogeneous one. For example, we can obtain

(*hu*)*t* +

� *hu*2 + 1 2 *gh*2 � *x*

where *M* is Manning's coefficient.

i.e.

� *hv*2 + 1 2 *gh*2 � *x*

All of the presented systems can be written in the compact matrix form

� *q*2 *a*

+ *ga*2 2*l* � *x*

rectangular channel

� *h*(*x*,*t*) 0

We now briefly introduce a problem describing fluid flow through the elastic tube represented by hyperbolic partial differential equations with the source term. In the case of the male urethra, the system based on model in (Stergiopulos at al., 1993) has the following form

$$\begin{aligned} a\_l + q\_\chi &= 0, \\ q\_l + \left(\frac{q^2}{a} + \frac{a^2}{2\rho\beta}\right)\_\chi &= \frac{a}{\rho} \left(\frac{a\_0}{\beta}\right)\_\chi + \frac{a^2}{2\rho\beta^2} \beta\_\chi - \frac{q^2}{4a^2} \sqrt{\frac{\pi}{a}} \lambda(\text{Re}), \end{aligned} \tag{10}$$

where *a* = *a*(*x*, *t*) is the unknown cross-section area, *q* = *q*(*x*, *t*) is the unknown flow rate (we also denote *v* = *v*(*x*, *t*) as the fluid velocity, *v* = *q a* ), *ρ* is the fluid density, *a*0 = *a*0(*x*) is the cross-section of the tube under no pressure, *β* = *β*(*x*, *t*) is the coefficient describing tube compliance and *λ*(*Re*) is the Mooney-Darcy friction factor (*λ*(*Re*) = 64/*Re* for laminar flow). *Re* is the Reynolds number defined by

$$Re = \frac{\rho q}{\mu a} \sqrt{\frac{4a}{\pi}}.\tag{11}$$

where *μ* is fluid viscosity. This model contains constitutive relation between the pressure and the cross section of the tube

$$p = \frac{a - a\_0}{\beta} + p\_{\varepsilon\prime} \tag{12}$$

where *pe* is surrounding pressure.

#### **3. Conservative and nonconservative problems and numerical schemes**

#### **3.1 Conservative problems and numerical schemes**

We consider the conservation law in the conservative form

$$\mathbf{q}\_t + [\mathbf{f}(\mathbf{q})]\_\mathbf{x} = \mathbf{0}, \; \mathbf{x} \in \mathbb{R}, \; t \in (0, T), \tag{13}$$

$$\mathbf{q}(\mathbf{x}, 0) = \mathbf{q}\_0(\mathbf{x}), \; \mathbf{x} \in \mathbb{R},$$

The numerical scheme based on finite volume discretization in the conservation form can be written as follows,

$$\mathbf{Q}\_{j}^{n+1} = \mathbf{Q}\_{j}^{n} - \frac{\Delta t}{\Delta \mathbf{x}} (\mathbf{F}\_{j+1/2}^{n} - \mathbf{F}\_{j-1/2}^{n}).\tag{14}$$

We use also the semidiscrete version of (14)

$$\frac{\partial \mathbf{Q}\_{j}}{\partial t} = -\frac{1}{\Delta \mathbf{x}} (\mathbf{F}\_{j+1/2} - \mathbf{F}\_{j-1/2}).\tag{15}$$

The relation (14) can be derived as the approximation of the integral conservation law at the interval *xj*−1/2, *xj*+1/2 from time level *tn* to *tn*+1

$$\begin{split} \frac{1}{\Delta x} \int\_{x\_{j-1/2}}^{x\_{j+1/2}} \mathbf{q}(\mathbf{x}, t\_{n+1}) d\mathbf{x} &= \frac{1}{\Delta x} \int\_{x\_{j-1/2}}^{x\_{j+1/2}} \mathbf{q}(\mathbf{x}, t\_n) d\mathbf{x} - \\ & \quad - \frac{1}{\Delta t} \left[ \int\_{t\_n}^{t\_{n+1}} \mathbf{f}(\mathbf{q}(\mathbf{x}\_{j+1/2}, t)) dt - \int\_{t\_n}^{t\_{n+1}} \mathbf{f}(\mathbf{q}(\mathbf{x}\_{j-1/2}, t)) dt \right]. \end{split} \tag{16}$$

variation of the solution is non-increasing, i.e.

**3.2 Nonconservative problems**

Δ*x*

*j*+1/2, **Q**+

[**A**−(**Q**−

*j*+1/2, **Q**+

these waves, which are related to the eigenvalues of matrix **A**(**q**).

Jacobi matrix of the system, fluctuations can be defined as follows

*j*+1/2, **Q**+

*j*+1/2, **Q**+

*j*−1/2, **Q**+

discontinuous initial condition. In the nonlinear case it has the form

is case of a continuous solution when the following equality holds

**q**˜�

**A**(**Q**−

**A**−(**Q**−

**A**+(**Q**−

**4.1 Riemann problem for conservative systems**

When the problem (22) is derived from the conservation form (13), i.e. **f**�

*∂***Qj** *∂t* = − 1

where **A**±(**Q**−

**4. Riemann problem**

where

∞ ∑ *j*=−∞

|**Q***n*+1

We consider the nonlinear hyperbolic problem in nonconservative form

*j*+1 − **Q***n*+1

*j* | ≤

Numerical Schemes for Hyperbolic Balance Laws – Applications to Fluid Flow Problems 39

for all time layers *tn*. This property is important for limitation of the oscillations in the solution.

**q**(*x*, 0) = **q**0(*x*), *x* ∈ **R**.

moving to the right or to the left. The directions are dependent on the signs of the speeds of

*j*−1/2) = **f**(**Q**−

*j*−1/2) = **f**(**Q**+

The Riemann problem is the special problem based on finite volume discretization with the

**q***t* + [**f**(**q**)]*x* = **0**, *x* ∈ **R**, *t* ∈ (0, *T*),

**Q***n*

We solve the transitions between two states, but this solution could not exists in general. These transitions can be rarefaction waves, shock waves or contact discontinuities. Rarefaction wave

**q**(*x*, 0) = **Q***n*

*j*+1/2) = **F**−

*j*+1/2, **Q**+

The numerical schemes for solving problems (22) can be written in fluctuation form

*j*+1/2) + **A**(**Q**−

∞ ∑ *j*=−∞

|**Q***n*

*j*+1 − **Q***n*

**q***t* + **A**(**q**)**q***x* = **0**, *x* ∈ **R**, *t* ∈ (0, *T*), (22)

*j*−1/2) + **A**+(**Q**−

*j*+1/2) are so called fluctuations. They can be defined by the sum of waves

*j*+1/2) − **f**(**Q**+

*j*+1/2 − **f**(**Q**−

*j*−1/2) − **F**+

*j* , *x* < *xj*+1/2,

*j*+1, *x* > *xj*+1/2.

*j*−1/2),

*j*+1/2),

*j*−1/2.

**q**(*x*, *t*) = **q**˜(*ξ*(*x*, *t*)), (26)

(*ξ*) = *α*(*ξ*)**r***p*(*ξ*) (27)

*j*−1/2, **Q**+

*j*−1/2)], (23)

(**q**) = **A**(**q**) is the

(24)

(25)

*j* | (21)

The previous relations lead to the following approximations of integral averages

$$\mathbf{Q}\_{j}^{n} \approx \frac{1}{\Delta \mathbf{x}} \int\_{\mathbf{x}\_{j-1/2}}^{\mathbf{x}\_{j+1/2}} \mathbf{q}(\mathbf{x}, t\_{\mathrm{n}}) d\mathbf{x}\_{\star}$$

$$\mathbf{F}\_{j+1/2}^{\mathrm{n}} \approx \frac{1}{\Delta t} \int\_{t\_{\mathrm{n}}}^{t\_{\mathrm{n}+1}} \mathbf{f}(\mathbf{q}(\mathbf{x}\_{j+1/2}, t)) dt. \tag{17}$$

Numerical fluxes **F***n j*+1/2 are usually defined by the approximate solution of the Riemann problem between states **Q***n j*+1 and **Q***n j* (this technique is called the flux difference splitting; it will be described in the following parts) or by the Boltzmann approach (flux vector splitting; it will be described later). In what follows, we use the notation **Q**+ *j*+1/2 and **Q**− *j*+1/2 for the reconstructed values of unknown function. Reconstructed values represent the approximations of limit values at the points *xj*+1/2. The most common reconstructions are based on the minmod function (see for example (Kurganov & Tadmor, 2000)) or ENO and WENO techniques (Cˇrnjaricˇ-Zicˇ at al., 2004).

If the exact solution of the problem has a compact support in the interval �0, *T*� then it is possible to show, that the scheme (15) is conservative, i.e.

$$\sum\_{j=-\infty}^{\infty} \mathbf{Q}\_j^{n+1} = \sum\_{j=-\infty}^{\infty} \mathbf{Q}\_j^n. \tag{18}$$

The uniqueness of discontinuous solutions to the conservation laws is not guaranteed. Therefore the additional conditions, based on physical considerations, are required to isolate the physically relevant solution. The most common condition is called entropy condition. The unique entropy satisfying weak solution *q* holds

$$[\eta(q)]\_{\mathfrak{f}} + [\mathfrak{g}(q)]\_{\mathfrak{x}} \le 0,\tag{19}$$

for the convex entropy functions *η*(*q*) and corresponding entropy fluxes *ϕ*(*q*) (for example see (LeVeque, 2004)).

If the conservation law (13) is solved by the consistent method in conservation form (15) the Lax-Wendroff theorem is valid (LeVeque, 2004): if the approximate function *Q*(*x*, *t*) converged to the function *q*(*x*, *t*)for the Δ*x*, Δ*t* → 0, the function *q*(*x*, *t*)is the weak solution of the problem (13).

Many theoretical results in the field of hyperbolic PDEs can be found in the literature. For example, for scalar problems with the convex flux function, the convergence of the some method to the entropy-satisfying weak solutions, is proven. It means that if the solution *q�* of the problem

$$\begin{array}{ll} q\_t + [f(q)]\_\mathbf{x} = \epsilon q\_{\mathbf{xx}}, & \mathbf{x} \in \mathbb{R}, 0 < t < T, \epsilon > 0, \\\ q(\mathbf{x}, 0) = q\_0(\mathbf{x}), & \mathbf{x} \in \mathbb{R}. \end{array} \tag{20}$$

exists and the limit *q*∗ = lim *�*→0+ *q�* exists too than *q*∗ is the entropy-satisfying weak solutions of the problem (20).

It is possible to define the appropriate properties of the methods. For example, in the case of scalar problem, TVD (Total Variation Diminishing) property, which ensures that the total variation of the solution is non-increasing, i.e.

$$\sum\_{j=-\infty}^{\infty} |\mathbf{Q}\_{j+1}^{n+1} - \mathbf{Q}\_{j}^{n+1}| \le \sum\_{j=-\infty}^{\infty} |\mathbf{Q}\_{j+1}^{n} - \mathbf{Q}\_{j}^{n}|\tag{21}$$

for all time layers *tn*. This property is important for limitation of the oscillations in the solution.

#### **3.2 Nonconservative problems**

4 Will-be-set-by-IN-TECH

 *xj*+1/2 *xj*−1/2

 *tn*+1 *tn*

it will be described in the following parts) or by the Boltzmann approach (flux vector

*j*+1/2 for the reconstructed values of unknown function. Reconstructed values represent the approximations of limit values at the points *xj*+1/2. The most common reconstructions are based on the minmod function (see for example (Kurganov & Tadmor, 2000)) or ENO and

If the exact solution of the problem has a compact support in the interval �0, *T*� then it is

The uniqueness of discontinuous solutions to the conservation laws is not guaranteed. Therefore the additional conditions, based on physical considerations, are required to isolate the physically relevant solution. The most common condition is called entropy condition. The

for the convex entropy functions *η*(*q*) and corresponding entropy fluxes *ϕ*(*q*) (for example see

If the conservation law (13) is solved by the consistent method in conservation form (15) the Lax-Wendroff theorem is valid (LeVeque, 2004): if the approximate function *Q*(*x*, *t*) converged to the function *q*(*x*, *t*)for the Δ*x*, Δ*t* → 0, the function *q*(*x*, *t*)is the weak solution of the problem

Many theoretical results in the field of hyperbolic PDEs can be found in the literature. For example, for scalar problems with the convex flux function, the convergence of the some method to the entropy-satisfying weak solutions, is proven. It means that if the solution *q�* of

*qt* + [ *f*(*q*)]*x* = *�qxx*, *x* ∈ **R**, 0 < *t* < *T*, *�* > 0,

It is possible to define the appropriate properties of the methods. For example, in the case of scalar problem, TVD (Total Variation Diminishing) property, which ensures that the total

∞ ∑ *j*=−∞

**Q***n*

**Q***n*+1 *j* =

splitting; it will be described later). In what follows, we use the notation **Q**+

∞ ∑ *j*=−∞ **q**(*x*, *tn*)*dx*,

*j*+1/2 are usually defined by the approximate solution of the Riemann

**f**(**q**(*xj*+1/2, *t*))*dt*. (17)

*j* . (18)

*j*+1/2 and

*j* (this technique is called the flux difference splitting;

[*η*(*q*)]*t* + [*ϕ*(*q*)]*x* ≤ 0, (19)

*q*(*x*, 0) = *q*0(*x*), *x* ∈ **R**. (20)

*q�* exists too than *q*∗ is the entropy-satisfying weak solutions of

The previous relations lead to the following approximations of integral averages

Δ*t*

**Q***n j* ≈ 1 Δ*x*

*j*+1 and **Q***n*

*j*+1/2 ≈ 1

**F***n*

Numerical fluxes **F***n*

(LeVeque, 2004)).

(13).

the problem

the problem (20).

exists and the limit *q*∗ = lim

*�*→0+

**Q**−

problem between states **Q***n*

WENO techniques (Cˇrnjaricˇ-Zicˇ at al., 2004).

unique entropy satisfying weak solution *q* holds

possible to show, that the scheme (15) is conservative, i.e.

We consider the nonlinear hyperbolic problem in nonconservative form

$$\mathbf{q}\_l + \mathbf{A}(\mathbf{q})\mathbf{q}\_\mathbf{x} = \mathbf{0}, \ x \in \mathbf{R}, \ t \in (0, T), \tag{22}$$

$$\mathbf{q}(\mathbf{x}, 0) = \mathbf{q}\_0(\mathbf{x}), \ x \in \mathbf{R}.$$

The numerical schemes for solving problems (22) can be written in fluctuation form

$$\frac{\partial \mathbf{Q}\_{\mathbf{j}}}{\partial t} = -\frac{1}{\Delta x} [\mathbf{A}^-(\mathbf{Q}\_{\mathbf{j}+1/2'}^- \mathbf{Q}\_{\mathbf{j}+1/2}^+) + \mathbf{A}(\mathbf{Q}\_{\mathbf{j}+1/2'}^- \mathbf{Q}\_{\mathbf{j}-1/2}^+) + \mathbf{A}^+(\mathbf{Q}\_{\mathbf{j}-1/2'}^- \mathbf{Q}\_{\mathbf{j}-1/2}^+)],\tag{23}$$

where **A**±(**Q**− *j*+1/2, **Q**+ *j*+1/2) are so called fluctuations. They can be defined by the sum of waves moving to the right or to the left. The directions are dependent on the signs of the speeds of these waves, which are related to the eigenvalues of matrix **A**(**q**).

When the problem (22) is derived from the conservation form (13), i.e. **f**� (**q**) = **A**(**q**) is the Jacobi matrix of the system, fluctuations can be defined as follows

$$\begin{array}{l} \mathbf{A}(\mathbf{Q}^{-}\_{j+1/2}, \mathbf{Q}^{+}\_{j-1/2}) = \mathbf{f}(\mathbf{Q}^{-}\_{j+1/2}) - \mathbf{f}(\mathbf{Q}^{+}\_{j-1/2}),\\ \mathbf{A}^{-}(\mathbf{Q}^{-}\_{j+1/2}, \mathbf{Q}^{+}\_{j+1/2}) = \mathbf{F}^{-}\_{j+1/2} - \mathbf{f}(\mathbf{Q}^{-}\_{j+1/2}),\\ \mathbf{A}^{+}(\mathbf{Q}^{-}\_{j-1/2}, \mathbf{Q}^{+}\_{j-1/2}) = \mathbf{f}(\mathbf{Q}^{+}\_{j-1/2}) - \mathbf{F}^{+}\_{j-1/2}. \end{array} \tag{24}$$

#### **4. Riemann problem**

#### **4.1 Riemann problem for conservative systems**

The Riemann problem is the special problem based on finite volume discretization with the discontinuous initial condition. In the nonlinear case it has the form

$$\begin{aligned} \mathbf{q}\_{l} + [\mathbf{f}(\mathbf{q})]\_{\times} &= \mathbf{0}, \; \mathbf{x} \in \mathbb{R}, \; t \in (0, T), \\ \mathbf{q}(\mathbf{x}, 0) &= \begin{cases} \mathbf{Q}\_{j}^{n} & \mathbf{x} < \mathbf{x}\_{j+1/2}, \\ \mathbf{Q}\_{j+1\nu}^{n} \; \mathbf{x} > \mathbf{x}\_{j+1/2}. \end{cases} \end{aligned} \tag{25}$$

We solve the transitions between two states, but this solution could not exists in general. These transitions can be rarefaction waves, shock waves or contact discontinuities. Rarefaction wave is case of a continuous solution when the following equality holds

$$\mathbf{q}(\mathbf{x},t) = \mathbf{\tilde{q}}(\mathbf{\tilde{q}}(\mathbf{x},t)),\tag{26}$$

where

$$\mathbf{\tilde{q}}'(\xi) = \mathfrak{a}(\xi)\mathbf{r}^p(\xi) \tag{27}$$

of continuous functions. The limiting behavior for *�* → 0 is strongly depend on "smoothing

Numerical Schemes for Hyperbolic Balance Laws – Applications to Fluid Flow Problems 41

Under a special assumption, the nonconservative product **A**(**q**)**q***x* can be understood as a Borel measure. In the following we introduce basic theorems and definitions, for details see

**Definition 4.1.** *A path* φ *in* Ω ∈ R*m is a family of smooth maps* �0, 1� × Ω × Ω → Ω *satisfying:*

≤ *k* |**q***r* − **q***l*|

*∂s* (*s*; **q***l*, **q***r*)

*∂s* (*s*; **q**2

*μ*(*o*) =

*o*

0 ), **q**(*x*+

**Remark 4.1.** *Borel measure μ is called nonconservative product and is usually written* [**A**(**q**)**q***x*]φ*.*

**Definition 4.2** (Weak solution)**.** *Let* φ *be a family of paths in the sense of definition 4.1. A function*

**A**(**q**) *∂***q** *∂x dx*,

0 )) *∂*φ

0 ) = lim *x*→*x*+ 0 **q**(*x*)*.*

*∂s* (*s*; **q**(*x*−

0 ), **q**(*x*+

(**q**)*then Borel measure,* [**A**(**q**)**q***x*]φ = **f**(**q**)*x is independent*

**q***t* + [**A**(**q**)**q***x*]φ = 0, (32)

0 )) *ds*,

*r* ∈ O *and almost all s* ∈ �0, 1�*.*

*l* , **q**2 *r* ) ≤ *K* **q**1 *l* − **q**2 *l* + **q**1 *r* − **q**2 *r*

**Theorem 4.1** (Dal Maso, Le Floch, Murat)**.** *Let* **q** : (*a*, *b*) → R*m be a function with bounded variation and* **A** : R*m* → R*m*×*m a locally bounded function. There exists a unique signed Borel*

(Gosse, 2001; Le Floch, 1989; Le Floch & Tzavaras, 1999; Pares, 2006). `

*•* φ(0; **q***l*, **q***r*) = **q***l and* φ(1; **q***l*, **q***r*) = **q***r*, ∀**q***l*, **q***r* ∈ Ω*,* ∀*s* ∈ �0, 1�, *• for each bounded set* O ∈ Ω*, there exists a constant k* > 0 *such that ∂*φ

*• for each bounded set* O ∈ Ω*, there exists a constant K* > 0 *such that*

*r* ) − *∂*φ

*l* , **q**1

*measure μ on* (*a*, *b*) *characterized by following properties:*

*1. if x* → **q**(*x*) *is continuous on an open set o* ∈ (*a*, *b*) *then*

*2. if x*0 ∈ (*a*, *b*) *is a discontinuity point of x* → **q**(*x*)*, then*

1

**A**

φ(*s*; **q**(*x*−

**q**(*x*) *and* **q**(*x*+

**q** ∈ [*L*∞(R × R+) ∩ BV*loc*(R × R+)]*m is a weak solution of system (31), if it satisfies*

0

0 ) = lim *x*→*x*− 0

*μ*(*x*0) =

*for any* **q***l*, **q***r* ∈ O *and almost all s* ∈ �0, 1�*,*

 *∂*φ *∂s* (*s*; **q**1

*for any* **q**1

*l* , **q**1 *r*, **q**2 *l* , **q**2

*where we denote* **q**(*x*−

*of the path* φ*.*

**Remark 4.2.** *In the case where* **A**(**q**) = **f**�

*as a bounded Borel measure on* R × R+*.*

out."

for any function *ξ*(*x*, *t*), where *α*(*ξ*) is a coefficient dependent on the function *ξ* and **R***p*(*ξ*) is the corresponding *p*-th eigenvector of the Jacobi matrix **f**� (**u**).

Shock waves and contact discontinuities are special cases of discontinuous solutions. The requirement that this solution should be (see (LeVeque, 2004)) a weak solution of the problem (25) leads to the following relation

$$s(\mathbf{q}^+ - \mathbf{q}^-) = \mathbf{f}(\mathbf{q}^+) - \mathbf{f}(\mathbf{q}^-),\tag{28}$$

where *s* is the speed of the propagation of the discontinuities. The relation (28) is known as the Rankine-Hugoniot jump condition. In some cases it is possible to construct the solution of Riemann problem as the sequence of the transitions between the discontinuous states (in the case of strong nonlinearity the discontinuous solution consists of the shock waves, in the case of linear degeneration solution consists of the contact discontinuities).

In the case of a linear problem **f**(**q**) = **Aq** it is known that the initial state of each characteristic variable *wp* is moving at a speed that corresponds to the eigenvalue *λp* of the matrix **A**. The solution is a system of constant states separated by discontinuities that move at speeds correspondent to eigenvalues. Therefore, the jump **Q***j*+1 − **Q***j* over *p*-th discontinuity can be expressed as,

$$(w\_{j+1}^p - w\_j^p)\mathbf{r}^p = \mathbf{a}^p \mathbf{r}^p,\tag{29}$$

where **r***p* is the *p*-th eigenvector of matrix **A**. The relation (29) represents the initial jump in the characteristic variable *wp* and at the same time **q** = **Rw**. Therefore, the solution of linear Riemann problem can be defined by the decomposition of the initial jump of the unknown function to the eigenvectors **r***p* of the Jacobi matrix **A** = **f**� (**q**)

$$\mathbf{Q}\_{j+1} - \mathbf{Q}\_j = \sum\_{p=1}^{m} a^p \mathbf{r}^p. \tag{30}$$

The discontinuities **W***p* = *αp***r***p* are called waves and they are propagated by the speeds *λp*. For details see (LeVeque, 2004).

#### **4.2 Riemann problem for nonconservative systems**

In this section we are interested in nonlinear systems in nonconservative form

$$\mathbf{q}\_t + \mathbf{A}(\mathbf{q})\mathbf{q}\_x = \mathbf{0}, \quad \mathbf{x} \in \mathbb{R}, \ t > 0,\tag{31}$$

$$\mathbf{q}\_0 \in [\mathbb{B}\mathbf{V}(\mathbb{R})]^m\_{\text{\textquotedblleft}\mathbf{R}}$$

where **q** ∈ R*m*, **q** → **A**(**q**) is smooth locally bounded matrix-valued map, matrix **A**(**q**) is strictly hyperbolic (diagonalizable, with real and different eigenvalues). We suppose that **A**(**q**)is not Jacobi matrix so it is not possible to rewrite nonconservative system in conservative form (13). Here, BV[(R)]*m* is function space contains functions with bounded total variation.

Above mentioned system makes sense only if **q** is differentiable. In the case when **q** admits discontinuities at a point, **A**(**q**) may admits discontinuities as well and **q***x* contains delta function with singularity at this point. Then **A**(**q**)**q***x* is product of Heaviside function with delta function and in general is not unique. There is possibility to smooth out" of this function over width *�*, for example by adding viscosity or diffusion. Then we get well defined product of continuous functions. The limiting behavior for *�* → 0 is strongly depend on "smoothing out."

Under a special assumption, the nonconservative product **A**(**q**)**q***x* can be understood as a Borel measure. In the following we introduce basic theorems and definitions, for details see (Gosse, 2001; Le Floch, 1989; Le Floch & Tzavaras, 1999; Pares, 2006). `

**Definition 4.1.** *A path* φ *in* Ω ∈ R*m is a family of smooth maps* �0, 1� × Ω × Ω → Ω *satisfying:*

$$\left|\frac{\partial\phi}{\partial s}(s; \mathbf{q}\_{I'} \mathbf{q}\_{I})\right| \leq k \left|\mathbf{q}\_{I} - \mathbf{q}\_{I}\right|$$

*for any* **q***l*, **q***r* ∈ O *and almost all s* ∈ �0, 1�*,*

6 Will-be-set-by-IN-TECH

for any function *ξ*(*x*, *t*), where *α*(*ξ*) is a coefficient dependent on the function *ξ* and **R***p*(*ξ*) is

Shock waves and contact discontinuities are special cases of discontinuous solutions. The requirement that this solution should be (see (LeVeque, 2004)) a weak solution of the problem

where *s* is the speed of the propagation of the discontinuities. The relation (28) is known as the Rankine-Hugoniot jump condition. In some cases it is possible to construct the solution of Riemann problem as the sequence of the transitions between the discontinuous states (in the case of strong nonlinearity the discontinuous solution consists of the shock waves, in the case

In the case of a linear problem **f**(**q**) = **Aq** it is known that the initial state of each characteristic variable *wp* is moving at a speed that corresponds to the eigenvalue *λp* of the matrix **A**. The solution is a system of constant states separated by discontinuities that move at speeds correspondent to eigenvalues. Therefore, the jump **Q***j*+1 − **Q***j* over *p*-th discontinuity can be

where **r***p* is the *p*-th eigenvector of matrix **A**. The relation (29) represents the initial jump in the characteristic variable *wp* and at the same time **q** = **Rw**. Therefore, the solution of linear Riemann problem can be defined by the decomposition of the initial jump of the unknown

The discontinuities **W***p* = *αp***r***p* are called waves and they are propagated by the speeds *λp*.

**q**0 ∈ [BV(R)]*m*,

where **q** ∈ R*m*, **q** → **A**(**q**) is smooth locally bounded matrix-valued map, matrix **A**(**q**) is strictly hyperbolic (diagonalizable, with real and different eigenvalues). We suppose that **A**(**q**)is not Jacobi matrix so it is not possible to rewrite nonconservative system in conservative form (13). Here, BV[(R)]*m* is function space contains functions with bounded total variation. Above mentioned system makes sense only if **q** is differentiable. In the case when **q** admits discontinuities at a point, **A**(**q**) may admits discontinuities as well and **q***x* contains delta function with singularity at this point. Then **A**(**q**)**q***x* is product of Heaviside function with delta function and in general is not unique. There is possibility to smooth out" of this function over width *�*, for example by adding viscosity or diffusion. Then we get well defined product

*m* ∑ *p*=1 (**u**).

*s*(**q**+ − **q**−) = **f**(**q**+) − **f**(**q**−), (28)

*j* )**r***p* = *αp***r***p*, (29)

*αp***r***p*. (30)

(**q**)

**q***t* + **A**(**q**)**q***x* = **0**, *x* ∈ R, *t* > 0, (31)

the corresponding *p*-th eigenvector of the Jacobi matrix **f**�

of linear degeneration solution consists of the contact discontinuities).

(*wp*

function to the eigenvectors **r***p* of the Jacobi matrix **A** = **f**�

**4.2 Riemann problem for nonconservative systems**

*j*+1 − *wp*

**Q***j*+1 − **Q***j* =

In this section we are interested in nonlinear systems in nonconservative form

(25) leads to the following relation

For details see (LeVeque, 2004).

expressed as,

*• for each bounded set* O ∈ Ω*, there exists a constant K* > 0 *such that*

$$\left| \frac{\partial \phi}{\partial s}(s; \mathbf{q}\_I^1, \mathbf{q}\_r^1) - \frac{\partial \phi}{\partial s}(s; \mathbf{q}\_I^2, \mathbf{q}\_r^2) \right| \le K \left( \left| \mathbf{q}\_I^1 - \mathbf{q}\_I^2 \right| + \left| \mathbf{q}\_r^1 - \mathbf{q}\_r^2 \right| \right)$$

*for any* **q**1 *l* , **q**1 *r*, **q**2 *l* , **q**2 *r* ∈ O *and almost all s* ∈ �0, 1�*.*

**Theorem 4.1** (Dal Maso, Le Floch, Murat)**.** *Let* **q** : (*a*, *b*) → R*m be a function with bounded variation and* **A** : R*m* → R*m*×*m a locally bounded function. There exists a unique signed Borel measure μ on* (*a*, *b*) *characterized by following properties:*

*1. if x* → **q**(*x*) *is continuous on an open set o* ∈ (*a*, *b*) *then*

$$
\mu(o) = \int\_{o} \mathbf{A}(\mathbf{q}) \frac{\partial \mathbf{q}}{\partial x} \, dx,
$$

*2. if x*0 ∈ (*a*, *b*) *is a discontinuity point of x* → **q**(*x*)*, then*

$$\mu(\mathbf{x}\_0) = \int\_0^1 \mathbf{A}\left(\phi(s; \mathbf{q}(\mathbf{x}\_0^-), \mathbf{q}(\mathbf{x}\_0^+))\right) \frac{\partial \phi}{\partial s}(s; \mathbf{q}(\mathbf{x}\_0^-), \mathbf{q}(\mathbf{x}\_0^+)) \, ds\_\prime$$

*where we denote* **q**(*x*− 0 ) = lim *x*→*x*− 0 **q**(*x*) *and* **q**(*x*+ 0 ) = lim *x*→*x*+ 0 **q**(*x*)*.*

**Remark 4.1.** *Borel measure μ is called nonconservative product and is usually written* [**A**(**q**)**q***x*]φ*.*

**Remark 4.2.** *In the case where* **A**(**q**) = **f**� (**q**)*then Borel measure,* [**A**(**q**)**q***x*]φ = **f**(**q**)*x is independent of the path* φ*.*

**Definition 4.2** (Weak solution)**.** *Let* φ *be a family of paths in the sense of definition 4.1. A function* **q** ∈ [*L*∞(R × R+) ∩ BV*loc*(R × R+)]*m is a weak solution of system (31), if it satisfies*

$$\mathbf{q}\_l + [\mathbf{A}(\mathbf{q})\mathbf{q}\_x]\_\phi = 0,\tag{32}$$

*as a bounded Borel measure on* R × R+*.*

and suppose family of path in the sense of definition 4.1. If **q**(*x*, *t*) is piecewice regular weak

Numerical Schemes for Hyperbolic Balance Laws – Applications to Fluid Flow Problems 43

where index *k* represents number of discontinuities in solution, *xk*(*t*) are discontinuity points

*j* − 1 Δ*x*

The amount of the quantity on the cell boundaries *xj*+1/2 can be splitted into contribution to cell *Ij*+1 and the contribution to the cell *Ij*. In other words, we split it into two terms **A**±

**Definition 4.4.** *Given a family of path* φ*, a numerical scheme is said to be path-conservation if it can*

*j*+1/2 = **A**±(**q**¯ *j*−*p*,..., **q**¯ *j*+*q*),

**A***j* = **A**(**q**¯ *j*−*p*,..., **q**¯ *j*+*q*),

**A**(φ(*s*; **q**−

*t* �*n*+1

*tn*

*a* (*Ij*)*dt*. By this way we get generalisation of classical conservative finite volume

*j*−1/2 + **A***j* + **A**−

1

0

∀**q***j* ∈ Ω, *j* = −*p*,..., *q*.

*xj*+1/2 �

*xj*−1/2

*is if the nonconservative system can be written as conservative one, path-conservative method become*

**A**±(**q**,..., **q**) = **0**, ∀**q** ∈ Ω, (43)

**A**(φ(*s*; **q***l*, **q***r*)) *∂*φ

(**q**) *then path-conservation scheme is consistent and conservative (that*

*k* , **q**+ *k* )) *∂*φ

*∂s* (*s*; **q**−

*k* , **q**+ *k* ) *ds*

*k* = lim *x*→*xk* (*t*)+ **q**(*x*, *t*), *δx*=*xk* (*t*) is Dirac measure at

*μφ*(*Ij*) *dt*. (41)

⎤

⎦*δx*=*xk* (*t*), (40)

*j*+1/2.

(44)

*t* � *n*+1 *tn μφ s dt* and

*j*+1/2]. (42)

*∂s* (*s*; **q***l*, **q***r*) *ds*

**A**(**q**)**q***x dx*. (45)

solution then for given time *t*, the Borel measure can be written in following form

⎡ ⎣ � 1

0

*k*

So the sum of this terms can be understood as a discrete representation of the

*j* − Δ*t* Δ*x* [**A**+

**A**(**q**)**q***x dx* + ∑

*k* = lim *x*→*xk* (*t*)− **q**(*x*, *t*), **q**+

**Q**¯ *n*+1 *j* = **Q**¯ *n*

**Q**¯ *n*+1 *j* = **Q**¯ *n*

**A**±, **A** : Ω*p*+*q*+1 → Ω *being continuous functions satisfying*

**A**±

**A**−(**q**−*p*,..., **q***q*) + **A**+(**q**−*p*,..., **q***q*) = �

**A**(**q**,..., **q**) =

*classical conservative and consistent). Notice here that* **A**± *are fluctuations, see 3.2.*

*μφ*(*o*) = *μφ*

**A***j* ≈

*1.*

*2.*

*3.*

*t* � *n*+1 *tn μφ*

methods, see (Pares, 2006), i.e. `

*be written in form (42), where*

**Remark 4.3.** *If* **A**(**q**) = **f**�

in the time *t* > 0 and **q**−

*a* (*o*) + *μφ*

the point *x* = *xk*(*t*). Then we can get

*s* = �

*o*

**Definition 4.3** (Entropy solution)**.** *Given an entropy pair* (*η*, *ϕ*)*(entropy, entropy flux) for (31), i.e. a pair of regular functions* Ω → R*, such that*

$$
\nabla \varphi(\mathbf{q}) = \nabla \eta(\mathbf{q}) \cdot \mathbf{A}(\mathbf{q}), \quad \forall \mathbf{q} \in \Omega. \tag{33}
$$

*A weak solution is said to be entropic if it satisfies the inequality*

$$\frac{\partial \eta(\mathbf{q})}{\partial t} + \frac{\partial \varphi(\mathbf{q})}{\partial \mathbf{x}} \le 0 \tag{34}$$

*in the sense of distribution.*

Function **q**(*x*, *t*) is weak solution if and only if, across a discontinuity with speed *ζ*, it satisfies generalized Rankine-Hugoniot condition,

$$-\mathcal{L}(\mathbf{q}\_{l'} - \mathbf{q}\_{l}) + \int\_{0}^{1} \mathbf{A}\left(\phi(s; \mathbf{q}\_{l'} \mathbf{q}\_{l'})\right) \frac{\partial \phi}{\partial s}(s; \mathbf{q}\_{l'} \mathbf{q}\_{l'}) \, ds = 0. \tag{35}$$

Now we define the Riemann problem for nonconservative strictly hyperbolic system:

$$\mathbf{q}\_l + \mathbf{A}(\mathbf{q})\mathbf{q}\_\mathbf{x} = \mathbf{0}, \quad \mathbf{x} \in \mathbb{R}, \; t > 0 \tag{36}$$

with an initial condition

$$\mathbf{q}(\mathbf{x},0) = \mathbf{q}\_0(\mathbf{x}) = \begin{cases} \mathbf{q}\_l \text{ pro } \mathbf{x} < 0, \\\\ \mathbf{q}\_l \text{ pro } \mathbf{x} > 0, \end{cases} \tag{37}$$

where **q***l*, **q***r* ∈ R*m* are vectors of constants.

**Theorem 4.2.** *Let* φ *be a family of path in the sense of definition 4.1. Assume that system (36) is strictly hyperbolic with genuinely nonlinear or linearly degenerate characteristic field and the family of path* φ *satisfies*

$$\frac{\partial \phi}{\partial \mathbf{q}\_1}(\mathbf{1}; \mathbf{q}\_0, \mathbf{q}\_0) - \frac{\partial \phi}{\partial \mathbf{q}\_1}(\mathbf{0}; \mathbf{q}\_0, \mathbf{q}\_0) = \mathbf{I}, \quad \forall \mathbf{q}\_0 \in \mathbb{R}^m. \tag{38}$$

*Then, for* |**q***r* − **q***l*| *small enough, the Riemann problem (36) and (37) has a solution* **q**(*x*, *t*) *with bounded variation, which depends only on x t and has Lax's structure. That is* **q**(*x*, *t*) *consist of m* + 1 *constant states separated by shock waves, rarefaction waves or contact discontinuities.*

The solution of the Riemann problem for nonconservative system is related to the solution of the Riemann problem for conservative system. The only difference is in the case of shock wave, precisely in Rankine-Hugoniot condition formulation.

Before we define a class of useful numerical methods, we introduce a brief motivation, which can be in details found in (Gosse, 2001). Suppose nonconservative system

$$\mathbf{q}\_l + \mathbf{A}(\mathbf{q})\mathbf{q}\_\mathbf{x} = \mathbf{0}, \quad \mathbf{x} \in \mathbb{R}, t > 0 \tag{39}$$

and suppose family of path in the sense of definition 4.1. If **q**(*x*, *t*) is piecewice regular weak solution then for given time *t*, the Borel measure can be written in following form

$$\mu^{\Phi}(o) = \mu\_{\mathfrak{a}}^{\Phi}(o) + \mu\_{\mathfrak{s}}^{\Phi} = \int\_{o} \mathbf{A}(\mathbf{q}) \mathbf{q}\_{\mathbf{x}} \, d\mathbf{x} + \sum\_{k} \left[ \int\_{0}^{1} \mathbf{A}(\phi(s; \mathbf{q}\_{k}^{-}, \mathbf{q}\_{k}^{+})) \frac{\partial \phi}{\partial s}(s; \mathbf{q}\_{k}^{-}, \mathbf{q}\_{k}^{+}) \, ds \right] \delta\_{\mathbf{x} = \mathbf{x}\_{k}(t) \prime} \tag{40}$$

where index *k* represents number of discontinuities in solution, *xk*(*t*) are discontinuity points in the time *t* > 0 and **q**− *k* = lim *x*→*xk* (*t*)− **q**(*x*, *t*), **q**+ *k* = lim *x*→*xk* (*t*)+ **q**(*x*, *t*), *δx*=*xk* (*t*) is Dirac measure at the point *x* = *xk*(*t*). Then we can get

$$\bar{\mathbf{Q}}\_{j}^{n+1} = \bar{\mathbf{Q}}\_{j}^{n} - \frac{1}{\Delta x} \int\_{t\_{n}}^{t\_{n+1}} \mu^{\Phi}(I\_{j}) \, dt. \tag{41}$$

The amount of the quantity on the cell boundaries *xj*+1/2 can be splitted into contribution to cell *Ij*+1 and the contribution to the cell *Ij*. In other words, we split it into two terms **A**± *j*+1/2. So the sum of this terms can be understood as a discrete representation of the *t* � *n*+1 *tn μφ s dt* and *t n*+1

**A***j* ≈ � *tn μφ a* (*Ij*)*dt*. By this way we get generalisation of classical conservative finite volume methods, see (Pares, 2006), i.e. `

$$
\bar{\mathbf{Q}}\_{j}^{n+1} = \bar{\mathbf{Q}}\_{j}^{n} - \frac{\Delta t}{\Delta \mathbf{x}} [\mathbf{A}\_{j-1/2}^{+} + \mathbf{A}\_{j} + \mathbf{A}\_{j+1/2}^{-}].\tag{42}
$$

**Definition 4.4.** *Given a family of path* φ*, a numerical scheme is said to be path-conservation if it can be written in form (42), where*

$$\mathbf{A}\_{j+1/2}^{\pm} = \mathbf{A}^{\pm} (\bar{\mathbf{q}}\_{j-p'} \dots \mathbf{\bar{q}}\_{j+q})\_{\prime}$$

$$\mathbf{A}\_{j} = \mathbf{A} (\bar{\mathbf{q}}\_{j-p'} \dots \mathbf{\bar{q}}\_{j+q})\_{\prime}$$

**A**±, **A** : Ω*p*+*q*+1 → Ω *being continuous functions satisfying 1.*

$$\mathbf{A}^{\pm}(\mathbf{q}\_{\prime}, \dots, \mathbf{q}) = \mathbf{0}, \quad \forall \mathbf{q} \in \Omega,\tag{43}$$

*2.*

8 Will-be-set-by-IN-TECH

**Definition 4.3** (Entropy solution)**.** *Given an entropy pair* (*η*, *ϕ*)*(entropy, entropy flux) for (31), i.e.*

*∂t* + *∂ϕ*(**q**) *∂x*

Function **q**(*x*, *t*) is weak solution if and only if, across a discontinuity with speed *ζ*, it satisfies

**A** (φ(*s*; **q***l*, **q***r*)) *∂*φ

**Theorem 4.2.** *Let* φ *be a family of path in the sense of definition 4.1. Assume that system (36) is strictly hyperbolic with genuinely nonlinear or linearly degenerate characteristic field and the family of*

*Then, for* |**q***r* − **q***l*| *small enough, the Riemann problem (36) and (37) has a solution* **q**(*x*, *t*) *with*

The solution of the Riemann problem for nonconservative system is related to the solution of the Riemann problem for conservative system. The only difference is in the case of shock

Before we define a class of useful numerical methods, we introduce a brief motivation, which

Now we define the Riemann problem for nonconservative strictly hyperbolic system:

*∂η*(**q**)

 1

0

**q**(*x*, 0) = **q**0(*x*) =

(1; **q**0, **q**0) − *∂*φ

wave, precisely in Rankine-Hugoniot condition formulation.

*∂***q**1

*constant states separated by shock waves, rarefaction waves or contact discontinuities.*

can be in details found in (Gosse, 2001). Suppose nonconservative system

∇*ϕ*(**q**) = ∇*η*(**q**) · **A**(**q**), ∀**q** ∈ Ω. (33)

**q***t* + **A**(**q**)**q***x* = **0**, *x* ∈ R, *t* > 0 (36)

(0; **q**0, **q**0) = **I**, ∀**q**0 ∈ R*m*. (38)

*t and has Lax's structure.  (35)

(37)

*a pair of regular functions* Ω → R*, such that*

generalized Rankine-Hugoniot condition,

where **q***l*, **q***r* ∈ R*m* are vectors of constants.

*bounded variation, which depends only on x*

*∂*φ *∂***q**1

*in the sense of distribution.*

with an initial condition

*path* φ *satisfies*

*A weak solution is said to be entropic if it satisfies the inequality*

− *ζ*(**q***r* − **q***l*) +

$$\begin{aligned} \mathbf{A}^{-}(\mathbf{q}\_{-p\prime}, \dots, \mathbf{q}\_{q}) + \mathbf{A}^{+}(\mathbf{q}\_{-p\prime}, \dots, \mathbf{q}\_{q}) &= \underset{\begin{subarray}{c} 0 \\ 0 \end{subarray}}{\text{A}} (\phi(s; \mathbf{q}\_{l\prime}, \mathbf{q}\_{l})) \frac{\partial \phi}{\partial s}(s; \mathbf{q}\_{l\prime}, \mathbf{q}\_{l}) \, ds \\ \forall \mathbf{q}\_{j} \in \Omega, j = -p\prime, \dots, q. \end{aligned} \tag{44}$$

*3.*

$$\mathbf{A}(\mathbf{q}, \dots, \mathbf{q}) = \int\_{\mathbf{x}\_{j-1/2}}^{\mathbf{x}\_{j+1/2}} \mathbf{A}(\mathbf{q}) \mathbf{q}\_{\mathbf{x}} \, d\mathbf{x}.\tag{45}$$

**Remark 4.3.** *If* **A**(**q**) = **f**� (**q**) *then path-conservation scheme is consistent and conservative (that is if the nonconservative system can be written as conservative one, path-conservative method become classical conservative and consistent). Notice here that* **A**± *are fluctuations, see 3.2.*

**5.1 Central schemes, local Lax-Friedrichs method**

where the approximation of local speed is defined by

*sj*+1/2 = max

The local waves can be based on the decomposition (49). It means

*∂***Q***j ∂t* = − 1

The numerical fluxes can be defined by the following relation as

*j*+1/2, **Q**+

*j*+1/2, **Q**+

*j*+1/2, **Q**+

*j*+1/2 − **Q**−

*j*+1/2) − **Ψ**±

*j*+1/2 = −*sj*+1/2(**Q**∗

**F***j*+1/2 = 1

2 [**f**(**Q**−

**A**−(**Q**−

**A**+(**Q**−

**A**(**Q**−

decomposition

*j*+1/2) − **f**(**Q**−

**Z**1

and the fluctuations as

**5.1.1 Steady states**

*j*+1/2 = **A**+

(*HV*)+

*j*+1/2 − (*HV*)−

(*HV*)+

 0

*j*+1/2

*j*+1/2 = (*HV*)−

Because there is no motion of the fluid, we have

problem.

**A**−

**f**(**Q**+

These schemes use estimate of upper bound of maximal local speed of the propagating discontinuities (for example see (Kurganov & Petrova, 2000)). They are based on the following

Numerical Schemes for Hyperbolic Balance Laws – Applications to Fluid Flow Problems 45

*j*+1/2 − **Q**∗

*j*+1/2)|, |*λp*(**Q**+

*j*+1/2) − *sj*+1/2(**Q**∗

*j*+1/2 = *sj*+1/2(**Q**+

(**F***j*+1/2 − **F***j*−1/2) + **Ψ***j* (52)

*j*+1/2 − **Q**−

*j*−1/2) − **Ψ**∓

*j* ,

*j*+1/2 = **0**, the **A**−

*j*+1/2)

*j*+1/2). (56)

*j*+1/2 − *H*−

 *A*−,1 *j*+1/2

*j*+1/2 + *H*−

*j*+1/2 − **Q**−

*j*+1/2)|}}. (50)

*j*+1/2 − **Q**∗

*j*+1/2), (49)

*j*+1/2). (51)

*j*+1/2) (53)

(54)

*j*+1/2 = **0**

. (55)

*j*+1/2 = *sj*+1/2(**Q**+

*p* {max{|*λp*(**Q**−

Δ*x*

*j*+1/2) = **Z**1

*j*+1/2) = **Z**2

*j*+1/2 = **A***j* = **0**. Here we present only the equality **A**+

can be derived by the similar way. The first component of decomposition (49) is

= *sj*+1/2(*H*+

*j*−1/2) = **f**(**Q**−

*j*+1/2 + **f**(**Q**+

*j*+1/2), **Z**2

We can use the scheme in conservation form based on numerical fluxes (15) or in nonconservative form based on the fluctuations (23). It can be used the scheme in the form

for the nonhomogeneous problems, where **Ψ***j* is a suitable approximation of the source term.

*j*+1/2)] − 1

*j*+1/2,

*j*+1/2.

The advantage of central schemes is the fact, that it is not necessary to solve a full Riemann

These schemes do not preserve general steady states for the problems with source terms, but only special one for shallow water equations. The steady state called "lake at rest", where *v* = 0, *h* + *b* = const. If we use the scheme in fluctuation form, we need the equalities

*j*+1/2 − *H*∗

 *A*+,1 *j*+1/2

*j*+1/2 = 0 ⇒ *H*∗

*j*+1/2)

*j*+1/2 = 1 2 (*H*+

−*sj*+1/2(*H*∗

2

*j*+1/2) − **f**(**Q**+

*sj*+1/2(**Q**+

#### **4.3 The properties of exact solution**

Consider a system in conservative form with special right hand side (this system agrees with shallow water system (5))

$$\mathbf{q}\_l + \mathbf{f}(\mathbf{q})\_\mathbf{x} = \mathbf{S}(\mathbf{q}) \sigma\_\mathbf{x}. \tag{46}$$

This system can be rewritten in a following homogeneous nonconservative form

$$\mathbf{u}\_l + \mathbf{A}(\mathbf{u})\mathbf{u}\_\mathbf{x} = \mathbf{0},\tag{47}$$

where

$$\mathbf{u} = \begin{bmatrix} \mathbf{q} \\ \sigma \end{bmatrix} \quad \text{and} \quad \mathbf{A}(\mathbf{u}) = \begin{bmatrix} \frac{\partial \mathbf{f}}{\partial \mathbf{u}}(\mathbf{u}) \ -\mathbf{S}(\mathbf{u}) \\ 0 \ 0 \end{bmatrix},$$

where *∂***f** *∂***u** (**w**) denotes Jacobian matrix of flux **f**. In order to define weak solutions of the system (47) we have to choose family of paths (see definition 4.1 and theorem 4.2). And following requirement become natural: if **u** = [**q**, *σ*] *T* is a weak solution of the nonconservative system (47) and *σ* is constant, then **q** must be weak solution of the homogeneous system of conservation laws, i.e. (46) with zero right hand side. By this requirement we impose addition condition of family of paths:

$$\text{if } \mathbf{u}\_l \text{ and } \mathbf{u}\_r \text{ are such that } \sigma\_l = \sigma\_r = \bar{\sigma}\_r \text{ then }$$

$$\Phi(s; \mathbf{u}\_l, \mathbf{u}\_r) = \bar{\sigma}\_r \quad s \in \left< 0, 1 \right>\. \tag{48}$$

For these systems, there are no difficulties with convergence. Moreover, the shock waves propagating in regions where *σ* is continuous are correctly captured independently of the choice of path. For details see (Le Floch at al., 2008).

#### **5. Approximate Riemann solvers**

There are many numerical schemes for solving (7) with different properties and possibilities of failing. The main types of the finite volume schemes are the central, upwind and central-upwind schemes. All these schemes relate together in variety of ways. We can interpret them as schemes with different deep knowledge about structure of the solution of Riemann problem. For example, the same scheme can be interpreted as the HLL solver (an approximate Riemann solver) or as the central–upwind scheme.

The main requirements on the numerical schemes are the consistency (in the finite volume sense, i.e. consistency with the flux function), the conservativity (if there is possibility to rewrite the problem to the conservative form it is required to have conservative numerical scheme), positive semidefiniteness (the scheme preserves nonnegativity of some quantities, which are essentially nonnegative from their physical fundamental) and the well-balancing (the schemes maintain some or all steady states which can occur). The next properties are the order of the schemes and stability. From the computational point of view there are other properties such as robustness, simplicity and computational efficiency. The second and third of these properties are typical for the central methods, but they are not common for scheme based on Roe's solver.

The steady states mean that the unknown quantities do not change in the time, i.e. **qt** = **0** in (7), and the flux function must balance the right hand side, i.e. [**f**(**q**)]*x* = ψ(**q**, *x*). Some schemes are constructed to preserve some special steady states for example called "lake at rest" in open channel flow problems, where there is no motion and the free surface is constant.

#### **5.1 Central schemes, local Lax-Friedrichs method**

These schemes use estimate of upper bound of maximal local speed of the propagating discontinuities (for example see (Kurganov & Petrova, 2000)). They are based on the following decomposition

$$\mathbf{f}(\mathbf{Q}^+\_{\mathbf{j}+1/2}) - \mathbf{f}(\mathbf{Q}^-\_{\mathbf{j}+1/2}) - \mathbf{P}^\pm\_{\mathbf{j}+1/2} = s\_{\mathbf{j}+1/2}(\mathbf{Q}^+\_{\mathbf{j}+1/2} - \mathbf{Q}^\*\_{\mathbf{j}+1/2}) - s\_{\mathbf{j}+1/2}(\mathbf{Q}^\*\_{\mathbf{j}+1/2} - \mathbf{Q}^-\_{\mathbf{j}+1/2}), \tag{49}$$

where the approximation of local speed is defined by

$$s\_{j+1/2} = \max\_{p} \{ \max \{ |\lambda^p(\mathbf{Q}\_{j+1/2}^-)|\_\prime |\lambda^p(\mathbf{Q}\_{j+1/2}^+)| \} \}. \tag{50}$$

The local waves can be based on the decomposition (49). It means

$$\mathbf{Z}\_{\mathbf{j}+1/2}^{1} = -s\_{\mathbf{j}+1/2}(\mathbf{Q}\_{\mathbf{j}+1/2}^{\*} - \mathbf{Q}\_{\mathbf{j}+1/2}^{-}), \qquad \mathbf{Z}\_{\mathbf{j}+1/2}^{2} = s\_{\mathbf{j}+1/2}(\mathbf{Q}\_{\mathbf{j}+1/2}^{+} - \mathbf{Q}\_{\mathbf{j}+1/2}^{\*}).\tag{51}$$

We can use the scheme in conservation form based on numerical fluxes (15) or in nonconservative form based on the fluctuations (23). It can be used the scheme in the form