CC BY
25
Hybrid numerical method for modelling of 2D viscous compressible flows
Bondarev A.E.
Keldysh Institute of Applied Mathematics Russian Academy of Sciences
bond@keldysh. ru
Abstract. This paper describes a numerical method for modeling of compressible viscous time-dependent 2D flows. Considered numerical method is based on hybrid implicit finite-difference scheme (WW scheme). The scheme can be referred to class of two-parametrical finite-difference schemes. Having second order accuracy for time and space and unconditional stability the scheme has also internal artificial viscosity regulated by the choice of weight parameters. The feature of controlled artificial viscosity allows one to avoid undesirable oscillations in solution. Being simple and effective the method is applied to some practical CFD problems such as: jets interaction, separation problems, optimizing analysis.
Key words: Navier-Stokes Equations, Viscous Flows, Hybrid Scheme
1 Introduction
Numerical modeling of viscous compressible flows based on solving the Navier-Stokes equations. There are many well-known numerical methods intended for 2D viscous compressible flows such as Maccormack method [MacCormack, 1970;MacCormack, 1982], Beam-Warming method [Beam&Warming, 1978] and many others. The methods for viscous compressible flows are thoroughly described in survey [Shang, 1985].
These papers form some general requirements for numerical methods in practice. To compute viscous compressible time-dependent flows one should apply numerical methods with following properties:
- At least second order accuracy in space and time,
- Unconditional stability,
- The absence of limitations for time-step,
- Convenience for programming and computing.
To obtain the most exact results of mathematical modelling for different physical processes the researcher should use the best qualities of different finite-difference schemes combined by different ways.
For instance, if solution is smooth enough in some parts of computational field, it is natural to apply high-order accuracy finite-difference schemes. To prevent the appearance of nonphysical oscillations in the vicinity of solution discontinuities one should apply schemes with artificial viscosity or monotone schemes having first order accuracy. Hybrid finite-difference schemes are applied with purpose to combine the most useful properties of different numerical schemes in one computational field.
For the simplest case hybrid scheme can be written as combination
S*G1+(1-S)*G2 (1)
where S - hybrid coefficient ( or weight parameter), G1 and G2 - numerical schemes, having different useful properties. For example, G, - first-order accuracy scheme, G2 - high-order accuracy scheme.
Most of applying to practical CFD problems finite-difference schemes is hybrid. According to review [Kulikovsky et al., 2001], such well-known algorithms as FCT (flux corrected transport), different types of TVD (total variation diminishing) schemes, ENO (essentially non-oscillatory) and WENO (weighted essentially non-oscillatory) schemes and many other popular schemes can be referred to the class of hybrid finite-difference schemes. Detailed description and classification for different types of hybrid finite-difference schemes is considered in review [Kulikovsky et al., 2001].
So, hybrid schemes are very useful because the researcher can choose the best from the properties of different schemes. At the same time some difficulties appear. The researcher needs detail analysis of properties and possible limitations for hybrid coefficients (or weight parameters) to keep the correspondence between physical process and numerical model.
The present paper describes numerical method, which agrees with mentioned above general requirements. It's ADI-method modification using hybrid implicit finite difference scheme. The scheme has second order accuracy in space and time. Also the scheme (we'll call this scheme as WW-scheme) is unconditionally stable and simple for programming. Except these properties WW-scheme has one interesting and useful feature. When non-linear problem with strong shocks is solved, one has to reduce undesirable solution oscillations. There are two ways for this. The first way is concerned with procedure of smoothing between time-steps. The second way consists in adding some terms with artificial viscosity to basic equations. Both ways require more calculations and complicate algorhytm. The present numerical method doesn't need these ways. Needed for stabilization of solution artificial viscosity is an internal property of WW-scheme. One can regulate the artificial viscosity by the choice of weight parameters. This property is quite suitable for practical applications.
2 Application to Burgers equation
Before discussing the solution of the complete Navier-Stokes equations, it's worthwhile to consider the basic elements of WW-scheme applied to Burgers equation.
a/ a/ a2/ (2)
— + c— — v—— dt dx dx
This equation is considered as one-dimensional analogy of the Navier-Stokes equations. We'll consider the simplest case with constant viscosity v and convection velocity c.
Hybrid numerical method for modelling of 2D viscous compressible flows Denote finite difference operators for approximation of space derivatives:
W=fL - &, sjr = /;.;, - 2/;+& (3)
With the help of finite difference operators (3) one can write the class of schemes with properties depending of the choice two parameters S1 and S2. We'll call these parameters as «weights». These schemes can be written in a form
fr1=sj? +(i-s2)<un]+
2 2 Ax
+-^[s2s2frl+v-s2)s2fin]
(Ax)
It is easy to note that for S2 = 0.5 finite-difference scheme (4) is a linear combination of Crank-Nicolson scheme type, having second order accuracy for time and space variables, and Lax scheme, having significant artificial viscosity.
1—S
Choosing sk =-1 the additional term containing artificial viscosity can be
T
presented as
h2 d2u (5)
2 dx2
Thus, WW-scheme can be considered as hybrid implicit unconditionally stable finite-difference scheme having second order accuracy for time and space variables. This scheme includes terms with artificial viscosity preventing non-physical oscillations near the shocks. The artificial viscosity is regulated by the choice of weight parameters. The important advantage of WW-scheme is the fact that one can directly control the artificial viscosity by regulating the meaning of weight parameter sk. This hybrid scheme can be quite effective being applied to computations of 2D viscous flows.
3 Application to Navier-Stokes equations
The unsteady compressible 2D form of Navier-Stokes equations is used for modeling viscous flows. These equations can be written in non-dimensional form as
3/ df
at ox
1 - H
--ja 2-
Re p y
3/ = 1 3 a S\+f
3y Re p 1 dx dx dy 2J" dy J
(6)
where
7 = (p,u,v,T)
(J)
Новые информационные технологии в автоматизированных системах 2014
4 у 3 Рг 4 у
F = (Fp,Fu,Fv,FT) The components for vector F = (Fp,Fu,Fv,FT) can be written as
Эм ЭУ — + —
дх Эу
\
■2-1, У
F _1ЭР+ A •
/? Эх Эу Эх 3 Эх Эу
"2 Э M [Л Эу i
3 Эх [y J У Эх JRc.P'
/г _ 1 dP | J Э ди_ .2 v_
р Эу )Эх Эу 3 Эу Эх 3 у^Эу
\
—+ 2—
У)
R
Rе.р
Эм ЭУ .У —+— + у —
Эх Эу у у
'a«Y ГэИ2
ЧЭху
+ 7
/ \2 V
У;
ЭУ ЭМ
—+—
Эх Эу
2/
л2
Эи Эу .У —+— + ] —
Эх Эу у.
(8)
in terms of density p, velocity components u and v for directions x and y, viscosity coefficient /i, temperature T, heat conductivity Â, Mach number , Reynolds number ReM, Prandtl number Pr, specific heat ratio y. Index j defines the type of problem under consideration. If j = 0, one considers flat 2D problem. For j = 1 we consider axis-symmetric problem in cylindrical coordinates.
The pressure P can be obtained by non-dimensional equation of state
P = pT (9)
We assume that viscosity coefficient ¡i and heat conductivity coefficient A are connected as
¿=m=M{t) (10)
The system of equations (6) is solved in assumption that each sought-for gasdynamic function p,u,v,Tcan be defined from corresponding equation of system (6). Each equation of system (6) can be presented as
df +hdf rd
di+aYx+bYy = cYx
dx
d
+ H
(ID
where a,b,c,d,e,g are corresponding coefficients. Part //contains the terms with mixed derivatives and the terms with functions (and derivatives) corresponding to other defining equations from (6).
The modified ADI-method is used for numerical solving the equations (6). For instance one can write equation (6) for x coordinate direction in a form
dt dx dx I dx J
a . df d
A = —b— + e— dy dy
(12)
+ H
Such equation can be solved by means of WW scheme (4) with corresponding choice of weight parameters (5). Using ADI-method one should apply the scheme to each coordinate direction by turns. The linearization of the non-linear coefficients is obtained by the procedure of linear extrapolation for time-steps. Using WW-scheme the solution of equation (12) amounts to solving corresponding tridiagonal matrix equation. This procedure should be implemented for each direction by turns.
4 Applications
Being quite simple and suitable for programming described above method was applied to many standard CFD problems for viscous compressible flows such as flows near obstacles, wake flows, boundary layers near the surface, jet interaction with obstacle etc. The method can be considered as quite effective for such 2D problems. It allows to define shock waves and vortex zones in the flows having good agreement with corresponding experimental results.
These properties allow using of numerical method for some more complicated problems. Also the method is used as a basic method for some program systems. The examples are described below.
Described numerical method is used for testing of program tool intended for hybrid finite-difference schemes optimization [Bondarev et al., 2013]. This paper contains the description of developed program tool Burgers2. The program tool is intended for tuning and optimization of computational properties for hybrid finite-difference schemes applied to Burgers equation. One-dimensional model Burgers equation describes propagation of disturbances for dissipative medium. The equation has exact solution, so it is widely used for tuning-up of computational tools. Described program tool is based on combining of optimization problem solution and visual data presentation. Visual presentations of maximal error surface and error function are implemented as program tool features. User is able to visualize error function distribution for
any chosen moment of time. These visual presentations allow analyzing and control computational properties of hybrid finite-difference schemes under consideration. Users have possibility of creating hybrid finite-difference schemes and analyzing computational properties for chosen grid template provided by program tool. Visual presentation of optimization problem solution allows finding of suitable weight coefficients for hybrid finite-difference scheme under consideration. WW scheme was used as test for program. The optimization solution allows to find minimal value for artificial viscosity (5) preventing oscillations. The choice of minimal value sk allows decreasing of error up to 0,5% in comparison with exact solution [Bondarev et al., 2013].
Another example of numerical method usage is described in paper [Bondarev, 2008a]. This paper presents a program set for numerical simulation and visual presentation of viscous heat-conductive flows in channels. Described above method is used as a basic method for this program set. The program set allows user simulating and visualizing flowfield in channels for different types of boundary conditions and initial conditions. Also the program set allows to solve different types of inverse problems for channels [Bondarev, 2008b].
Another one example of numerical method application is devoted to practical analysis of unsteady circulating zones transformation. The results are presented in paper [Bondarev, 2010]. The problem of unsteady interaction of the supersonic viscous flow with jet obstacle is considered.
This obstacle appears due to co-current underexpanded jet exhausting from the nozzle. The nozzle is placed to external supersonic viscous flow. Expanding jet propagates on the external surface of the nozzle and creates obstacle in external flowfield. The obstacle disturbs external flow and circulating zone appears ahead the obstacle.
We consider a problem containing time-dependent boundary condition for underexpanded jet. Jet pressure ratio was set at the nozzle edge as time-dependent function n = n(t) = Pa/P„ (where Pa- jet pressure, P„ - external flow pressure). The full system of time-dependent Navier-Stokes equations for viscous compressible heat-conductive flow is used as mathematical model. The dependence n = n(t) is chosen as linear function. It allows one to set different rates of pressure ratio growth up to n = 100.
As a result of calculations of direct problem a new space-time structure is formed. Increasing the rate of pressure ratio growth one obtains new space-time structure in the vicinity of circulating zone ahead the jet obstacle. This new structure is described thoroughly in paper [Bondarev, 2010].
This problem is considered also from the point of view of parametric optimizing analysis in paper [Bondarev&Galaktionov, 2013]. The formation of new space-time structure in the flow is considered as unsteady event in question. The rate of pressure ratio growth is chosen as control parameter. The case of four determining parameters is considered. These four parameters are Mach
number - M„, Reynolds number - Re„, Prandtl number -Pr^ and Shx - Strouhal number for the problem under consideration. For each fixed set of these numbers (Af^Re^Pr^S/O the inverse problem is solved by varying pressure ratio growth rate until the onset of the new structure formation in the flow. Characteristic parameters vary in definite ranges [Bondarev&Galaktionov, 2013].
So for each set of characteristic (M^Re^Pr^SAJ one defines the exact values for crucial velocity of pressure ratio growth V* when the new flow structure appears. Parallel algorithm is implemented for computations. For the space of determining parameters two types of grids are chosen: 5 and 10 points for each determining parameter. It requires computing 625 and 10000 inverse problems. The computations are performed by parallel cluster K100 (Keldysh Institute of Applied Mathematics RAS, Moscow, Russia). MPI technology is applied to control parallel computations.
As a result of approach application four-dimensional data array is obtained. This array contains numerical presentations of crucial velocity V* dependence on four determining parameters (Af^Re^Pr^S/O.
The analysis of variances for each characteristic parameter and construction of different 3D data projections for various triplets of determining parameters allow to decrease the number of dimensions up to three.
So we are able to consider 3D array V* = v\M. Analyzing visual presentation of the array one can approximate the isosurfaces by planes. For the purpose of rough estimation the sought-for dependence can be approximated by plane. It allows one to get average estimation of V* and its dependence on determining parameters as
V'=V\M_,Pi_,ShJ = -0.1M_+Q.U5Pi_ + 0.2ASh_ (13)
This example illustrates that described above numerical method can be used for inverse problems and optimization problems of 2D viscous compressible time-dependent flows.
5 Conclusions
The paper describes a numerical method for simulating of 2D viscous compressible time-dependent flows. The method is based on a hybrid finite difference scheme having artificial viscosity regulating by the choice of scheme weight parameters. Being quite simple and suitable for programming the method is implemented in many applications for modeling of 2D viscous compressible flows.
6 Acknowledgements
This work is supported by RFBR grants (project 13-01-00367A and project 14-01-00769A).
7 REFERENCES
[Beam&Warming, 1978] Beam R., Warming R.F. "An implicit factor scheme for the compressible Navier-Stokes equations", AIAA Journal, 16(4), 385-402 (1978). [Bondarev, 2008a] Bondarev A.E. Razrabotka programmnogo kompleksa dlya chislennogo modelirovaniya I vizual'nogo predstavleniya vyazkih teploprovodnyh techeniy v kanalah [The elaboration of program set for numerical simulation and visual presentation of viscous heat-conductive flows in channels], Preprints of Keldysh Inst. Of Applied Mathematics, 73, 16 p. (2008).
[Bondarev, 2008b] Bondarev A.E. Reshenie obratnyh zadach v programmnom komplekse dlya modelirovaniya vyazkih techeniy v kanalah [The inverse problems in program set for simulation of viscous flows in channels], Preprints of Keldysh Inst. Of Applied Mathematics, 74, 12 p. (2008).
[Bondarev, 2010] Bondarev A.E. "Analiz nestacionarnogo vzaimodeystviya vyazkogo sverhzvukovogo potoka so struynoy pregradoy s pomoschyu resheniya zadachi optimizacii opredelyauschih parametrov [Analysis of viscous supersonic flow interaction with jet obstacle by determining parameters optimization]", Proc. of the 13-th Workshop "New information technologies", Moscow, 212-217 (2010).
[Bondarev et al., 2013] Bondarev A.E., Bondarenko A.V, Galaktionov V.A., Mihailova T.N., .N., Ryzhova I.G., "Design of program tool Burgers2 for hybrid finite difference schemes optimization and visualization", Scientific Visualization, 5(1), 26-37 (2013). [Bondarev&Galaktionov, 2013] Bondarev A.E , Galaktionov V.A. "Parametric Optimizing Analysis of Unsteady Structures and Visualization of Multidimensional Data", Int. J. of Modeling, Simulation and Scientific Computing, 4(supp01), 13 p., DOI 10.1142/S1793962313410043, (2013).
[Kulikovsky et al., 2001] Kulikovsky A.G., Pogorelov N.V., Semenov A.Yu. Matematicheskie voprosy chislennogo resheniya giperbolicheskih sistem uravneniy [Mathematical questions of hyperbolic systems numerical solution], Moscow, 'PhysMathLit" Ed., 608 p. (2001).
[MacCormack, 1970] MacCormack R.W. "Numerical solution of the interaction of a shock wave with a laminar boundary layer", Proc. of the Second Int. Conf. on numerical methods in fluid dynamics, Univ.of California, 151-163 (1970).
[MacCormack, 1982] MacCormack R.W. "A numerical method for solving the equations of compressible viscous flow", AIAA Journal, 20(9), 1275-1281 (1982).^ [Shang, 1985] Shang J.S. "An assessment of numerical solutions of the compressible Navier-Stokes equations", Journ. of Aircraft, 22(5), 353-370 (1985).
