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

The Euler and Navier-Stokes (NS) equations are derived by applying Newton's second law to an infinitely small inviscid and viscous control volume respectively, and with the extension of the mass and energy conservation laws, they provide the highest level of approximations for the flow physics within approaching of continuum-mechanics based flow regime. The mentioned governing equations1 are the system of the nonlinear partial differential equations and they do not have general closed form solution as yet. However, in consideration of increasing expectations arisen from the industry and the high level evolution of the computer technology, the different numerical and optimization methods are developed and implemented in the complex framework of CFD (Computational Fluid Dynamics) programs. These software are widely spread in the engineering practice and they provide efficient solutions for different industrial problems. Beside the advantages of the virtual prototyping, the wide range of visualization tools and optimization, significant number of measurements can be replaced by the validated numerical analyses, which reduce time, capacity, cost and strongly contribute to the benefit of the company. Hence, followed by the short overview of different numerical methods generally used in CFD, a complete physical and mathematical interpretation is presented for a compressible NS and Euler solvers with validation and extension for coupling of inviscid flow solver with inverse design based optimization algorithm in a framework of the finite volume method.


Introduction
The Euler and Navier-Stokes (NS) equations are derived by applying Newton's second law to an infinitely small inviscid and viscous control volume respectively, and with the extension of the mass and energy conservation laws, they provide the highest level of approximations for the flow physics within approaching of continuum-mechanics based flow regime.The mentioned governing equations1 are the system of the nonlinear partial differential equations and they do not have general closed form solution as yet.However, in consideration of increasing expectations arisen from the industry and the high level evolution of the computer technology, the different numerical and optimization methods are developed and implemented in the complex framework of CFD (Computational Fluid Dynamics) programs.These software are widely spread in the engineering practice and they provide efficient solutions for different industrial problems.Beside the advantages of the virtual prototyping, the wide range of visualization tools and optimization, significant number of measurements can be replaced by the validated numerical analyses, which reduce time, capacity, cost and strongly contribute to the benefit of the company.Hence, followed by the short overview of different numerical methods generally used in CFD, a complete physical and mathematical interpretation is presented for a compressible NS and Euler solvers with validation and extension for coupling of inviscid flow solver with inverse design based optimization algorithm in a framework of the finite volume method.

Numerical methods for fluid flow
The continuum mechanics based NS equations describe flow physics in continuum and slip flow regime defined by Knudsen number as Kn<0.01 and 0.01<Kn<0.1 respectively (Zucrow & Hoffman, 1976).In case of Euler equations, in which the viscous (diffusion) effects are not considered, Re→∞ and Kn→0.The continuum approach does not count the individual molecules and instead, considers the substance of interest as being continuously distributed in space (Wassgren, 2010).The continuum based method requires that the smallest length Finite Volume Method -Powerful Means of Engineering Design 4 scale to be much larger than the microscopic length scale, typically the mean free path of a molecule (for gases).The Knudsen number for a generic engineering flow, in which the mean free path (λ) is around 1*10E-5 cm and the macroscopic length of interest (L) is 1 mm, is Kn=λ/L=0.0001(Wassgren, 2010).The other expectation is that the highest length scale of the spatial resolution must be small enough to accurately capture the parameters (e.g.density) to be approximately constant in the control volume and do not be affected by the physics or geometry.The continuum assumption is valid in the vast majority of the industrial applications and so the NS (and Euler) equations can be used in the most of the industrial applications.However, in general, they do not have closed form solution till now.Hence, different discretization methods were developed from the second part of the last century till now to have approximate results for the nonlinear partial differential system of the equations in each point of the temporal (in case of transient problem) and spatial discretization.Beside fulfilling consistency, stability and convergence characteristics as a measure of the mathematical correctness of the discretization methods, the final results of the numerical analyses are evoluted as a function of the applied governing equations, boundary conditions and the geometry.
The three most frequent discretization methods (in the percent of the available commercial CFD codes) are the finite difference (FDM) (~ 2 %), finite element (FEM) (~ 15 %) and finite volume (FVM) (~ 80 %) methods.The rest 3 % are consist mostly of Spectral, Boundary element, Vorticity type and Lattice gas or Lattice Boltzmann methods.
The most traditional and oldest methods applied for numerical solution of partial differential equations are the FDM.They are essentially based upon the properties of the Taylor series expansion and the method is only applicable to structured grids in practice (Manna, 1992).The accuracy of FDM method strongly depends upon the mesh size and its properties as stretch ratio, aspect ratio and skewness for example.Although the increasing number of mesh point improves the accuracy, it can lead to difficulties in solution procedures due to the matrix inversion at the algebraic system of the equations (Manna, 1992).
The FEM historically originated from structural mechanics.The physical domain is divided by cells or elements, they form the numerical mesh, which can be structured or unstructured providing higher flexibility for handling complex geometry than FDM (Manna, 1992).The field variables are approximated by linear combinations of known basis functions, which can be quite general with varying degrees of continuity at the inter-element boundaries (Hirsch, 2007).Its mathematical rigor and elegance makes the FEM algorithms very attractive and widely researched area in the CFD community (Manna, 1992).
The FVM is originally introduced by McDonald 1971 and they are based on the observation that the conservation laws have to be interpreted in integral form to preserve discontinuous solution as vortex sheets, contact discontinuities or shock waves.Similarly to the FEM, the physical domain is subdivided by cells.The flow field variables are evaluated in some discrete points on each cell and they are interpreted as average value over the finite volumes.The finite volumes are constructed as parts of one or more cells, which can, moreover, be either overlapped or non-overlapped.The conservation laws are then applied to the finite volumes to obtain the discrete equations.The possibility of modifying the shape and location of the control volumes associated to a given mesh allows large freedom in the choice of the function representation of the flow field.This property is not shared by either the FDM or FEM and it is mostly the reason of higher popularity of the FVM in the engineering applications (all paragraph from Manna, 1992).
Turn back to the physical interpretations of modelling, there are two main approaches can be distinguished.The Pressure-and the Density-based classes of methods have evolved distinct strategies for the discretization, non-linear relaxation and linear solution aspects underlying the computational schemes.Historically, the pressure-based approach was developed for lowspeed incompressible flows, while the density-based approach was mainly used for highspeed compressible flows.However, recently both methods have been extended and reformulated to solve and operate for a wide range of flow conditions beyond their traditional or original intent.In the Pressure-based methods the velocity is obtained from the momentum equations.The pressure field is extracted by solving a pressure or pressure correction equation, which is derived by manipulating continuity and momentum equations.The one of the most widespread methods is the SIMPLE (Semi-Implicit Pressure Linked Equations) family of schemes.Karki and Patankar developed the SIMPLER method for compressible flows, applicable for a wide range of speeds (Karki, Patankar 1989).Munz et al. extended the SIMPLE scheme for low Mach number flow employing multiple pressure variables, each being associated with different physical response (Munz et al. 2003).The time marching Densitybased methods represent a large class of schemes adopted for compressible flows and applied widely in computational fluid dynamics for modelling steady and transient, transonic, supersonic and hypersonic flows.The continuity equation is used to determine the density field while the pressure distribution is obtained from the equation of state.The velocity field is computed also from the momentum equations.Both approaches are now applicable to a broad range of flows (from incompressible to highly compressible), but the origins of the densitybased formulation may give it an accuracy (i.e.shock resolution) advantage over the pressurebased solver for high-speed compressible flows (ANSYS, Inc. 2010).Hence, this method is used in followings due to the high speed aeronautical applications.

Turbulence modelling
Concerning the physical level of modelling turbulent flows in CFD, the three most frequent approaches are the DNS (Direct Numerical Simulation), LES (Large Eddy Simulation) and RANS (Reynolds-averaged Navier-Stokes equations) based simulations, meanwhile for modelling laminar flows, there are no special treatment of the original NS equations are required.The closest model to the real phenomenon is the DNS, in which the Navier-Stokes equations are numerically solved without any turbulence model.It means that the whole range of spatial and temporal scales of the turbulence is resolved.The computational cost is extremely high, it is proportional to Re 3 (Pope, 2000).This method can not be used in the most part of the engineering practice due to the economical aspects.In case of LES, the large scales of turbulence are resolved directly providing more accurate results compared with RANS, meanwhile the geometry-independent, small scales and expensive structures are modelled.Although this approach is less computationally expensive compared with DNS, the industrial use is confined to low Reynolds number from purely computational consideration.The RANS based methods are the highest feasible level of approximation for the turbulent flows in general engineering applications.In this case, the parameters in the governing equations are averaged over a characteristic time interval in order to eliminate the influence of the turbulent fluctuations, meanwhile the unsteadiness of the other physical www.intechopen.comFinite Volume Method -Powerful Means of Engineering Design 6 phenomenon are preserved (Manna, 1992).The resulting equations are formally identical to the laminar Navier-Stokes equation, except for some extra terms, so called Reynolds stresses, which results from the non linear terms of the governing equations (Manna, 1992).The main goal of the different turbulence modelling is to provide functional relations, expressions and/or closure equations related with these new terms in order to couple and so make definite the governing equations (in discretized form).In the last forty years several turbulence models have been developed.Following the Boussinesq approximation (which expects the similarity of the mechanism of laminar and turbulent stresses) the Reynolds stresses can be expressed in terms of eddy viscosity, which is considered as the turbulent counterpart of the laminar molecular viscosity.Opposite to the molecular viscosity, which is a fluid property, the eddy viscosity is a function of the flow properties (Manna, 1992).The direct aim of the turbulence models is to identify the functional relations between the flow properties and the turbulent eddy viscosity.At the lowest level, they are based upon the mixing length concept, introduced by Prandtl in 1925, which effectively relates the turbulent shear stresses to the mean velocity gradients.These algebraic turbulence models are often called zero equation models.A higher degree of approximation is reached by solving additional equations written for the turbulence variables, as the transport equations for the turbulent kinetic energy and its rate of dissipation.This class of turbulence models is usually classified according to the number of additional equations applied for the turbulence variables, i.e. one equations models, two equation models or higher order closure like the Reynolds stress models (Manna, 1992).The two-equation models, especially the k-ε and k-ω based models, are the widest spread applications in the industry due to the best compromise of the physical accuracy and computational cost.However, beside the advantages, they have disadvantages also.
The first low Reynolds number k-ε model has been developed by Jones and Launder (1973) and suppose that the flow is fully turbulent.It is a computationally cheap and provides reasonable accuracy for a wide range of flows.However, the k-ε model performs poorly for complex flows involving severe pressure gradient, separation and strong streamline curvature.From the standpoint of aerodynamics, the most disturbing problem is the lack of sensitivity to adverse pressure-gradients.Under those conditions, the model predicts significantly too high shear-stress levels and thereby delays (or completely prevents) separation.Furthermore, it requires the application of damping mechanism for stabilization when the equations are integrated through the viscous sublayer (Menter, 1994).The standard k-ε model has been modified by many authors.
There are a significant number of alternative models that have been developed to overcome the shortcomings of the k-ε model.One of the most successful, with respect to the accuracy and the robustness, is the k-ω model of Wilcox (1988).It solves one equation for the turbulent kinetic energy k and a second equation for the specific turbulent dissipation rate ω.The most prominent advantages is that the equations can be integrated without additional terms through the viscous sublayer, which makes them y+ insensitive and provides straightforward application of the boundary conditions.This leads to significant advantages in numerical stability.The model performs significantly better under adverse pressure-gradient conditions and separation than the k-ε model and suitable for the complex boundary layer flows (e.g.external aerodynamics and turbomachinery).However, the k-ω model also has some shortcomings.The model depends strongly on the free stream values of ω that are specified outside the shear-layer.Another point of concern is that the model predicts spreading rates that are too low for free shear-layers, if the correct values are specified for ω (Menter, 1994).
In order to improve both the k-ε and the k-ω models Menter (1994) suggested to combine the two models called the SST (Shear Stress Transport) k-ω turbulence model.The use of a k-ω formulation in the inner parts of the boundary layer makes the model directly usable all the way down to the wall through the viscous sub-layer, hence the SST k-ω model can be used as a Low-Re turbulence model without any extra damping functions.The SST formulation also switches to the k-ε behaviour in the free-stream and thereby avoids the common k-ω problem that the model is too sensitive to the free-stream value of the turbulence variables (in particular ω).The further distinct of the SST turbulence model is the modified turbulence eddy-viscosity function.The purpose is to improve the accuracy of prediction of flows with strong adverse pressure gradients and pressure-induced boundary layer separation.The modification accounts for the transport of the turbulent shear stress, which is based on Bradshaw's assumption that the principal shear stress is proportional to the turbulent kinetic energy (Blazek, 2005).Due to the above mentioned characteristics, the Menter model has gained significant popularity in the aeronautical community and can be regarded as one of the standard approaches today.Despite of the improvements of the original SST model as SST-2003, SST-sust, SST-Vsust variants for example, there are also some complaints.The one of them is that the distance to the nearest wall has to be known explicitly.This requires special provisions on multiblock structured or on unstructured grids (Blazek, 2005).Also, the k-ε to k-ω switch can produce some unrealistic effective viscosity, which may not affect the results.Meanwhile the SST turbulence model provides similar benefits as standard k-ω (except for the free stream sensitivity), the dependency on wall distance can make it less suitable for the free shear flows compared to standard k-ω and it requires mesh resolution near to the wall (CFD Online Discussion Forum, 2010).The Wilcox's improved k-ω model (Wilcox, 1998) predicts free shear flow spreading rates more accurately than version 1988.The results of the benchmark simulations are in close agreement with flow measurements on far wakes, mixing layers, and plane, round, and radial jets, and it thus applicable to wall-bounded flows and free shear flows (ANSYS, 2010).Hence, this turbulence model has been implemented.

Finite volume method based compressible flow solver
Nowadays, in spite of disadvantages of turbulence closure models for RANS (Reynolds Averaged Navier-Stokes equations), they are at present the only tools available for the computation of complex turbulent flows of practical relevance.Their popularity comes from high efficiency in terms of accuracy and computational cost, which makes them widely used in commercial codes and related multidisciplinary applications.Hence, a detailed description of the physical and mathematical aspects of a RANS based compressible flow solver is presented in followings.
The governing equations in conservative form are derived by using density weighted averaging coupled with the time averaging of RANS.The code is based on structured, density-based cell centered finite volume method, in which the convective terms are discretized by Roe approximated Riemann method.The method of Roe is highly nondissipative and closely linked to the concept of characteristic transport.It is one of the most powerful linear Riemann solvers due to the excellent discontinuity-capturing property including shear waves.However, it is well-known that flux function mentioned above can produce non-physical expansion shocks that violate the entropy condition.This can be avoided by modifying the modulus of the eigenvalues for the non-linear fields.The method of Yee is used and discussed at the present case to cure the problem.Central discretization is applied for diffusive terms on a shifted mesh.MUSCL (Monotone Upstream Schemes for Conservation Laws) approach is implemented for higher order spatial reconstruction with Mulder limiter for monotonicity preserving.Wilcox k-ω two equations turbulence model is adopted and used (Wilcox, 1998).The explicit system of the equations is solved by the 4 th order Runge-Kutta method.The numerical boundary conditions are determined by the extrapolation technique for the NS solver and by the method of characteristics at the Euler solver.The interest is mostly in high speed industrial and aeronautical applications hence, the validation is completed for test cases are in the transonic, supersonic and subsonic flow regime as circular bump in the transonic channel and compression corner for the NS solver and flow over a wing profile and cascade for the Euler solver.The description of the benchmarks and the results are presented in Chapter 3.

Governing equations
In absence of external forces, heat addition, mass diffusion and finite rate chemical reaction, the unsteady two dimensional Navier-Stokes equations coupled with k-ω turbulence modelequations (Wilcox, 1998) in conservative, divergence and dimensional form are next:


of the Cartesian coordinate system, where x, y Є R and t Є R + .The conservative variables and convective fluxes are given by while the viscous fluxes and source terms are following: The bar over variables represents the Reynolds averaging over the characteristic time scale in order to separate and filter the small sized phenomena as turbulence fluctuation: For a supersonic or hypersonic compressible flow the local density is not constant and in case of turbulent flow it fluctuates also due to the pressure diffusion, dilatation, work and turbulent transport/molecular diffusion of turbulent energy.Hence, the instantaneous density can also be separated by averaging and fluctuating part, which requires the introduction of Favre averaging given by ( 5).
The F in superscripts and tides (~) over the parameters in ( 2) and ( 3) and in followings mean the Favre averaged parameters.Other relationships in ( 2) and ( 3) are given by ( 6)-( 9), in which p is the static pressure,  is the density, k is the turbulent kinetic energy, e is the internal energy u and v are the Cartesian components of velocity vector, v c is the specific heat at constant volume, T is the static temperature,  is the dynamic molecular are the constants in the Sutherland's formula to count the effect of temperature on dynamic viscosity and  is the specific turbulent dissipation rate.The terms in the expressions, which are related to values of j , i and k in indexes, range from 1 to 2. The closure expression of the NS equations is the ideal gas law (10).

T R p  
(10) The not mentioned parameters and expressions in the turbulence model equations are given by ( 11)-( 18) (Wilcox, 1998), where H is the Heaviside step function and c is the sound speed.The system of the nonlinear partial differential equations is already coupled.After discretization the system of algebraic equations can easily be solved.
Assuming a frictionless and isentropic flow, the NS equations -neglecting viscous and heat conducting terms -can be reduced to the Euler equations, which are the highest level approximation of the inviscid flow.

Boundary conditions
The numerical treatment of the boundary conditions strongly influences not only the convergence properties but the accuracy of the results in solving partial differential system of the equations.The physical boundary conditions secure the existence and uniqueness of the exact solution and numerical boundary conditions are supposed to ensure that various perturbations generated in the interior of the computational domain leave it without being reflected at the boundaries.Due to the convection dominated problem, the method of characteristic is used to determine the number and the exact values of the numerical boundary conditions in case of Euler equations, meanwhile extrapolation technique is applied for the NS equations.The direction of wave propagation (V n , V n , V n +c and V n -c) depends not only on the sign of the cell face normal velocity V n but also on the local speed of sound c.At the boundary, the number of physical boundary condition to be imposed equals the number of negative eigenvalues, which correspond to the incoming characteristics from the outside (boundary) to the computational domain.The need for numerical boundary conditions comes from the fact that the actual problem to be solved is formulated in terms of the conservative variables rather than Riemann invariants.Therefore, it is hard to impose the Dirichlet boundary conditions in the usual way.It is common practice to recover the boundary values by switching to the characteristic variables, evaluating the incoming Riemann invariants from the physical boundary conditions and extrapolating the outgoing ones from the interior of the computational domain (Kuzmin & Möller, 2004).
Concerning the inlet, it is examined whether the flow is supersonic or subsonic.First, consider the supersonic case, at which only incoming characteristics are available.Hence, total pressure, static pressure, total temperature and flow angle are imposed as physical boundary conditions and no numerical boundary conditions are required for the Euler equations.If the flow is subsonic, one outgoing characteristic is appeared (V n -c), so the two dimensional local Riemann problem belongs to that characteristic curve is solved by using physical and computed (existing) parameters and the total pressure, total temperature and flow angle are imposed as physical boundary conditions.The temperature and the components of the velocity vector are recovered by using ideal gas law and inlet flow angle, while the tangential velocity component is kept to be constant.Concerning the NS equations, additionally to the above mentioned specifications, the turbulent kinetic energy (k) and specific dissipation rate (ω) are imposed as physical boundary conditions and the static pressure is extrapolated from the computational domain in case of subsonic inlet.
If the outcoming flow is supersonic, there is no incoming characteristic hence, no physical boundary conditions are specified.If the flow is subsonic, there is one incoming characteristic hence, one parameter (static pressure) is imposed as physical boundary condition.The numerical boundary conditions are calculated by using characteristic variables (compatibility equations) in case of Euler equations, or they are extrapolated from the interior as the NS equations are implemented for viscous flow modelling.The static temperature is calculated by ideal gas law.
Concerning the Euler equations, the solid wall boundary conditions are considered as an outlet with the restriction of normal velocity is set to be zero across the wall.Hence, the numerical boundary conditions are calculated by using characteristic variables (compatibility equations) belongs to characteristic curves V n , V n and V n +c.The static temperature is calculated by ideal gas law also.In case of NS equations, the no-slip boundary condition is implemented at the solid walls, the velocity vectors are set to be zero.Assuming zero pressure gradients, the pressure is set equal to the one at the cell centre nearest to the wall.Adiabatic wall condition is used to determine temperature.The turbulent kinetic energy is zero at the wall and the specific dissipation rate is computed by suggestion of (Wilcox, 1998) assuming a rough wall.
If the flow field has any kind of periodicity, the calculation time can be reduced significantly by using periodic boundary condition, by which the rotationally or translationally shifted parameters are used in the cells at the boundaries.Periodic boundary condition is used before and after the profile to recover infinite blade number in cascades.
Meanwhile the expected pressure distribution is imposed at the solid wall boundary in the inverse mode of the inviscid solver, the opening boundary is used instead of solid wall to allocate the local flow direction -determined by the pressure difference between the boundary and computational domain -and its velocity V n (see Fig. 11.).The main outcome of the present mode is to have velocity profile over the geometry, which will be used for modifying mesh points in the wall modification module of the inverse design method (see Chapter 4.).
The detailed description of the presented boundary conditions for the Euler equations is found in (Veress et al., 2011).

Finite volume discretization
The finite volume method is a technique to handle the spatial derivatives that are appeared in the governing equations.The method is based on the integration of the equations over a finite volume.Then the integrals are transformed using the Gauss' divergence theorem where applicable.The physical meaning of the method is that fluxes flow through the faces of the finite volume while flux balance over the volume is satisfied.In the finite volume approach the first issue consists in evaluating the contour integral of the inviscid and viscous flux vectors in (1), hence the numerical flux functions are written in vector form given by ( 19).The x e  and y e  are the unit vectors in x and y directions.
It is convenient to define total fluxes normal to the boundary of elementary control volumes rather than making use of the individual Cartesians components, using the rotational invariance of the governing equations.Integrating system eq.( 1) over a control volume  , which is bounded by interface  and applying the Gauss' divergence theorem gives ( 20) and ( 21) (Manna, 1992), is the local outward pointing unit normal vector of the cell interface.The variables in conservative form, inviscid fluxes, source terms and viscous fluxes are the followings: where By means of finite volume discretization, in order to pass from a continuous to a discrete form, the unknown in a general finite volume of the partitioned computational domain is defined by ( 25), which corresponds to cell centre discretization.The vector j , i U has been interpreted as a mean value over the control volume.The fluxes are computed across quadrilateral cells in a structured grid, which can be seen in Fig. 1.The computational domain is divided by finite number of non overlapping finite surfaces or cells and the ( 21) is applied for each cell separately.It means that the second integral in ( 21) is replaced by summation over the all boundaries b N of the cell j , i and so eq.( 21) can be written in the general form of the semidiscrete expression over cell j , i as it is shown in (26).
x y e y e x , which takes into account the sign of the Jacobian matrices, or in other words the relevant propagation directions between the L and R states (Manna, 1992).The can be evaluated by linear wave decomposition where an unique average state (which is denoted by a hat) of the left and right states exist (Roe, 1981): For the ideal gas, Roe has shown that the matrix n D ˆ is equal to the Jacobian n D when it is expressed as a function of the variables  ˆ, u ˆ, v ˆ, and 0 h ˆ, which are weighted variables of the square root of density.0 h is the total enthalpy.Detailed information about the Roe's method of the approximate Riemann solver is found in (Roe, 1981).The method of Roe is highly non-dissipative and closely linked to the concept of characteristic transport.It is one of the most powerful linear Riemann solvers due to the excellent discontinuity-capturing property including shear waves.However, it is well-known that flux function mentioned above can produce non-physical expansion shocks that violate the entropy condition.This can be avoided, by modifying the modulus of the eigenvalues for the non-linear fields.The method of Yee (1989) is used at the present case.MUSCL (Monotone Upstream Schemes for Conservation Laws) approach is implemented for higher order spatial extension, by which the piece-wise constant distribution of the initial variables over the cell can be replaced by a piecewise linear or quadratic one.The mathematical deduction starts with the introduction of Taylor series expansion around point i.The results are found at (28) after discretization and integration.28) corresponds to a third order accurate space discretization in one dimensional problem (Manna, 1992).The spurious oscillations (wiggles) can occur with high order spatial discretization schemes due to shocks, discontinuities or sharp changes in the solution domain.Hence, in this case, Mulder limiter is implemented in the high resolution schemes for monotonicity preserving (Manna, 1992): where The same method was used for NS and for the turbulence model equations, however they were handled separately.
A simple central scheme is applied for the space discretization of the diffusive terms in (26) as follows: L U and R U are the conservative variables at the cell centres.The derivatives in the diffusive terms are determined at the centre of the cell interface.Hence, two new types of cells are formed, which are shifted by 2 1  i and 2 1  j directions respectively.The centres of the boundary of such cells are coincident with the centres and vertices of the original cells.In case of former situation, the parameters are known, because they are stored at the cell centre of the original cells.If the centres of the boundaries of the new cells are coincident with cell vertex of the original cells, the flow variables at the new cell boundary centre are the simple averages of the four neighbouring cell centre values of the original cells in case of quadrilateral mesh.Then, the derivatives into the x and y directions of the viscous flux function ( 23) can be obtained by using Green-Gauss theorem: where  is an arbitrary flow variable, ' j , i  is the area of the shifted cell i,j, k runs through the number of the boundaries of the shifted cell till the b N , which is the number of maximal boundary, kx ' n and ky ' n is the x and y components of the cell boundary normal unit vector of the interface k, k , ij  is the value of the flow variable at the given interface centre of the shifted cell i,j and k , ij '  is the length of the face k.At the boundaries of the computational domain, a series of the ghost cells are applied and filled with values, so no special treatment is necessary to determine the derivatives.The geometry of the ghost cells are extrapolated from the last two cells.
The integral of the source term in (21) are approximated as follows: is an average value over the cell j , i "  .Derivatives in the expressions are determined at the cell centres by using similar treatment as in case of diffusive term discretization.
A widely used class of non linear multi-stage time integration techniques is given by the Runge-Kutta (RK) schemes.They are usually designed to obtain higher order temporal accuracy with minimum computational storage and the large stability range with the specific coefficients, even though it has been often used for steady state calculations as herein.The 4 stages RK method (RK4) is used to solve the time derivatives of the conservative variables in (26) with for simplicity in each cell given by: for the Euler equations are the coefficients of RK4, n represents the parameters at previous time step and n+1 at the next time step over a cell.The RK4 index is denoted by k and it runs from 1 to m with its maximum value of 4 in step 2 (35).Due to the steady state assumption, the time accuracy is not required hence, the RK4 coefficients are applied to have high stability and smoothing properties of the upwind scheme with MUSCL reconstruction.In order to optimize the time step behind the stability criterion, the local time stepping has been used for every cells j , i as follows (Lefebvre, Arts, 1997): N is the number of cell boundaries and c is the sound speed.

Validation of the flow solver
The goal of the validation -in case of any calculation methods -is to provide information about the correct mathematical and physical operation of the simulation by means of comparing the results with real tests or other benchmarks especially referring to the application of flow physics under investigation.

Validation of the viscous flow solver
In the following sections, the numerical results are presented for transonic channel over circular bump and compression corner for validating frictional and heat conducted flow simulations.
The first test case is transonic channel over circular bump, in which the flow enters into the channel with Mach number 0.85 and a shockwave develops over the circular bump.The bump has 4.2 % maximum thickness.At the inlet, the total pressure, total temperature, and flow angles are specified as physical boundary conditions.The static pressure corresponds to the isentropic flow at Mach=0.85 is imposed at the outlet.Under these conditions the flow expands in the rear part of the bump up to a Mach number 1.2 and ends up into a week shock wave to allow the recovery of the free stream conditions.The results of FLUENT and own code are compared to each other at the same solver settings however, the k-ε turbulence model was used in the commercial program.The Mach number iso-lines show reasonable deflections (see Fig. 2.), the present method predicts the shockwave earlier.The shape of the geometry and the thickness of the boundary layer have a dominant effect on the shock wave evolution.Different numerical methods and turbulence models have different inherent mechanism to model boundary layer and shock wave-boundary layer interaction.The one of the criticisms against the k-ε turbulence is the lack of sensitivity to adverse pressure-gradients (see Subchapter 1.2.).The boundary layer seems to be thinner at the downstream of the circular bump in case of the commercial code compared to the own one.Hence, the shockwave triggered earlier in case of the present model.However, the differences between the two approaches in the entire computational domain are less than Fig. 3. Configuration and Schlieren photograph about compression corner at inlet Mach number 3 and with slope 18 ° (left side) (Settles, 1975) and Mach number distribution by the recent solver at inlet Mach number 2.85 and with slope 20 ° (right side) 5 % hence, based on the strongly validated commercial code, the accuracy of the in-house software can be accepted for this benchmark.
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, along which velocity distributions are measured is shown in Fig. 4.  -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 layer at upstream and the ω.

Validation of the inviscid flow solver
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 www.intechopen.comMost 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: p tot,in =112800 [Pa]; inlet total temperature: T tot,in =293.15[K]; outlet static pressure: p stat,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 by using NACA 65-410 profile also with: ,1 stat p are the dynamic and static pressure respectively at wall surface position 1. ,1 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 C p 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.

Finite volume method based optimization and inverse design for aeronautical applications
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 (p req ) 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 (p req ) 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 (p req in Fig. 11.), which are the topic of the following paragraph and subchapter respectively.
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.).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.

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.
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) 40) that the equation describes a flow that is ready everywhere to separate.Stratford presents the following solutions (Stratford, 1959),  (Smith, 1975) eq. ( 42) eq. ( 43) In that two-part solution, 0 x is the start of pressure rise, , 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, . 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: 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   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: p and maximum velocity 0 u is supposed to be constant starting from the leading edge of the suction side till the starting of the positive pressure gradient ( 0x ).The Mach numbers 0 M at these points are calculated by: The 0 T , 0 u and 0  are obtained by the energy equation of the isentropic flow and ideal gas law: 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 0 p .All parameter belongs to 0 p can be calculated by eqs.( 48)-( 52).The next step is to find location 0 x , 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 ( 0 x ) 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.0 p , 0 T , 0 u , 0  , 0 x and ) (x p (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: p tot,in =112800 [Pa]; inlet total temperature: T tot,in =293.15[K]; outlet static pressure: p stat,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 www.intechopen.com(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.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 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.

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.
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 p tot,in = 107853 [Pa], the inlet total temperature is T tot,in = 298.42[K], the inlet flow angle is 30° and the outlet static pressure is p stat,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.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.

Pressure Distribution
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

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.

Fig. 2 .
Fig. 2. Mach number distribution in transonic channel over circular bump test case (dotted line: FLUENT, continuous line: recent solver)

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 y EXP 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 layer at upstream and the ω.

Finite
Fig. 8. Mach number and pressure distribution over the profile NACA 65-410 at α=4 degrees angle of attack Fig. 9. Lift coefficient in the function of angle of attack at the profile NACA 65-410 where 5 . 1  

Fig. 11 .
Fig. 11.Flowchart of the computational procedure of the iterative inverse design calculation While the incoming and out coming velocity distribution (see V n 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:

Fig. 12 .
Fig. 12. Schematic view of the wall modification process based on the local velocity vector Fig. 13.Stratford limiting flows at two values of unit Reynolds number(Smith, 1975)

Fig. 14 .
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)

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

Fig. 17 .
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) 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 p tot,in = 107853.4[Pa], the inlet total temperature is T tot,in = 298.4267[K], the inlet flow angle is 30° and the outlet static pressure is p stat,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 .
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) ,