CC BY
34
UDC 519.63+532.5-1/-9
Modification of the Numerical Code for Gas-Dynamical Flows
in Cylindrical Coordinates
E. A. Filistov
Department of Physics Moscow State University of Civil Engineering 26, Jaroslavskoje shosse, Moscow, Russia, 129337
The goal of this article is to develop a robust and accurate numerical method for solving hyperbolic conservation laws in three dimensions.
The basic equations are the three-dimensional Euler equations describing the motion of an inviscid gas.
The mathematical description of the model is represented by the system of equations of continuity, motion and energy (three dimensional nonstationary partial differential equations). We used the equation for adiabatic motion in this article.
The numerical method for solution of the gas-dynamical equations in strict divergent form has been used in this work. The three-dimensional numerical code for perfect non-stationary gas-dynamical flows simulation in cylindrical coordinates is constructed. This code is based on the explicit quasimonotonic, first-order TVD scheme. This scheme admit introduction of the limits on the anti-diffusion flows, which enhances the approximation order (to third order in the spatial coordinates) with minimal numerical dissipation and preservation of the monotonicity of the scheme.
In order to ensure numerical stability, the time step is restricted by a well-known Courant-Friedrich-Lewy stability condition.
The proposed scheme is comparable to the high order over the classical TVD schemes. Our scheme has the added advantage of simplicity and computational efficiency. The numerical tests which were fulfiled by the author in additional researches, validated the robustness and effectiveness of the proposed scheme.
Key words and phrases: gas dynamics, numerical simulation.
1. Introduction
The theoretical foundations of high-resolution TVD schemes for homogeneous scalar conservation laws and linear systems of conservation laws have been firmly established through the work of Harten [1], Sweby [2], and Roe [3-5]. These TVD schemes seek to prevent an increase in the total variation of the numerical solution, and are successfully implemented in the form of flux-limiters or slope limiters for scalar conservation laws and systems. However, their application to conservation laws in strict divergent form written in cylindrical coordinates is still not fully developed and can be improved. In this work we construct the three-dimensional numerical code for simulation of a perfect non-stationary gas-dynamical flows. This code is based on the explicit quasimonotonic TVD-scheme. We use strictly divergent form laws of conservation in cylindrical coordinates. This work is continuation and development of the author work [6] where study is limited to the calculation of two-dimensional free flows.
2. Basic Equations of Gas Dynamics
The basic equations of gas dynamics govern the motion of a perfect non-stationary, inviscid gas-dynamical flows by conserving mass, momentum and energy (three dimensional nonstationary partial differential equations). They are simplifications of the more general Navier-Stokes equations which include the effects of viscosity on the flow (Landau and Lifshitz, 1979 [7]).
Received 25th February, 2013.
2.1. Basic Equations
Using the divergence form, we have the following macroscopic equations:
dtp = -dj (pvJ)
(1)
for the conservation of mass of flow. This expression is often known at the continuity equation;
dt(pv *) = -d0 (pv V + P g") (2)
for the conservation of momentum;
t. 2
9t (^y + pE^J = -dj (v> (py + p £ + P
))
(3)
for the conservation of total energy.
Here p (r, t) denotes the fluid density, v (r, t) — the bulk flow velocity, P(r, t) — the thermal pressure, s = e(P, p) — the internal energy per unit mass (specific internal energy), gij — the metric tensor, r — the radius vector, and t — the time. In order to close this system of equations (1)-(3) and fully describe the fluid with this model, we need a fourth equation, an equation of state.
In condition of local equilibrium, the scalar pressure P can be expressed as a function of two other thermodynamic variables through an equation P = P(p, s). An equation of state is necessary, and we assume the perfect gas equation, P = (7 — 1)ep, where 7 is the usual ratio of heat capacities (the index of adiabatic).
2.2. Equations on a Cartesian Coordinate System
In three space dimensions the equations of gas dynamics describing the motion of an inviscid gas may be written in the following compact dimensionless differential vector form of the conservation law in an inertial Cartesian coordinate system:
dtq + + dy g + dz h = 0,
(4)
where q is the vector of conserved variables (the state vector), and F, G, H are the flux vectors in the corresponding coordinate directions given by
q
p pu pv pw
pu F2 puv puw
Pv , F = puv , G = G3 , H = pvw
pw puw pvw H4
\pe) puh pvh \pwhj
(5)
and v = {u, v, w}, u, v, w, are the components of velocity along the x, y and z directions,
?2 „-¡2 p p
£ =T + £,
v2 P h = - + £ + -, 2 p
(7 - 1)P '
F2 = pu2 + P , G3 = pv2 + P , H4 = pw2 + P .
Equations (4) are nonlinear hyperbolic partial differential equations for the unknown q and must be solved with suitable boundary conditions. Even though in (5) F, G, H are expressed as functions of the primitive variables p, u, v, w, P they can also be expressed as functions of the conserved variables q.
Hyperbolicity of system equations (4) requires that the matrixes dF/dq, dG/dq, dH/dq have real eigenvalues and be diagonalizable.
2.3. Governing Equations in Cylindrical Coordinates
The governing equations can be re-expressed in any coordinate system. For many engineering applications involving non-rectangular general geometries, a non-orthogonal body-fitted coordinate system (£, r/, s, t) is desirable, which is formally related to the physical coordinate system (x, y, z, t) as follows:
T = T (t), £ = £(x,y,z,t), 'q = 7](x,y,z,t), ç = ç (x,y,z,t).
(6)
In equations (6) t is the time, and £, r/, ? are the three curvilinear coordinate directions.
Making use of certain mathematical manipulations, the Cartesian derivatives can be are replaced by the curvilinear counterparts, and the gas-dynamical equations can be recast in the so called "strong conservation form". The details of the manipulations can be found in [8], for example. The transformed equations written in strong conservation form in the generalized curvilinear coordinate system read:
dT q + dz F + Ôv G + dç H = 0,
where
q, F
j,
- & q + &F + Cy G + & H),
G = 1(vt q + VxF + % G + Vz H), H = 1(st q + Ç*F + Çy G + H),
(7)
J is the Jacobian of the transformation given by:
j =
D(x,y,z)
£,x £,z
Vx Vy Vz
Çx Çy Çz
(8)
In an inertial frame in cylindrical coordinates (r, ft, z) the Jacobian of the transformation (8) is J = r-1, and equations of gas dynamics (4) also can be written in the strict divergent form of hyperbolic conservation laws. The system is:
dtq + drF + d#G + dzH = 0 . (9)
The expressions (7) in this case take the following form:
f P\ / pvr \
pu pvr u + Pßi2
r ■ Pv , F = r ■ pVr V + Pßi3
pw pvr w
\pE) \ pvr h
( \ / pw \
pv$u + P&12 pwu
PV#V + P&13 , H = r ■ pwv
pv$ w pw2 + P
pv# h \ pwh /
where we used the the local density p = p(r, ft, z, t), the local velocity v = v (r, ft, z, t) = (vr,v#,w), a12 = cos ft, a13 = sin ft, b12 = — sin ft, b13 = cos ft, u = vra12 + v#b12, v = vra13+Vfib13, w are the Cartesian components of velocity, the pressure P = P(r, ft, z,t),
and
r, v2 ^ , ^ , V2 P
E = + £ , P = (7 - 1) ep, h = —+ - + £, 2 2 p
as before.
3. Numerical Code
In this section, we construct the three-dimensional numerical code for perfect non-stationary gas-dynamical flows simulation on the Eulerian cylindrical grid. This code is based on the explicit quasimonotonic TVD-scheme. Our approach applies to the Roe method [3-5]. This an explicit finite-difference scheme is first order of approximation.
3.1. Numerical Approach
In cylindrical coordinates this scheme is represented as follows. Consider the system of equations (9), which here we rewrite in finite-difference form
A<i,j,k + Fi+1/2,j,k — Fi — 1 /2,j,k + Gqi,j + 1/2,k - Gi,j-i/2,k + Hi,j,k+1 /2 — HCi,j,k—1/2 = „
T + Ar + At + A =
where Aqi,J-,fc = q^+l - q^^ and q^+j1 = q(ri, tj, zk, tn+i) is the solution at the time step n + 1 (n is the time level); t, Ar, At?, Az are the steps of the grid. The numerical fluxes after using the definition of Jacobian Matrices,
= w • 8 = OG • c = f •
oq oq oq
are calculated as
q , — Fi,3,k + ~Fi+1,j,k 1 I\m(±.*\|A ,m < i+1/2,j,k = -2--(<C )|A ^+1/2,^(<C ) •
A^ = diag{vr — c, vr , vr , vr , vr + c} , A S\+1/2,j,k = ((ri+1pi+1,j,k — r-ipi,j,k) T q*C* (Vri+1,j,k — Vri,j,k^ •
A S2+1 /2,j,k = 2"*'?* (V#i+1,j,k — V#i,j,k) ,
2c
A S i+1/2,i,fc = 2^*0*(wi+1,j,k — Wi,j,k)
A Si+1/2,j,fc = 21*2 (C*2((ii+1,3,k — %j,k) — (Ti+1Pi+1,j,k — ripi,j,k^ ; the matrix of the corresponding right eigenvectors takes the form:
nA =
( 1 0 0 2 1
u* — C*d12 2c*612 0 2u* u* + c*a12
v* — c*a13 2 c* bV3 0 2v* v* + c*a13
* w* 0 2 2 w* w*
h* — c* v* 2 c* v* 2 c*w* v*2 h* + c*v* /
where
c=\/7i p/p
is the adiabatic speed of sound, q = rp, and in which q*, v *, Roe average.
They can be obtained from
Q* = VQi+i,j,k Qi,j,k j
^ * \/Qi+l,j,k Vi+\,j,k + JQi,j,k Vi,j,k v = _v--- -
y/Qi+i,j,k + a/ Qi,j,k h* = hl+l,J,k + ^^
^<li+\,j,k + y/qi,j,k
f
* I VQi,3,k Ci,j,k + //Qi+l,j,k ci+l,j,k , P y/Qi,j,k Qi+\,j,k _
C = v-m-+ -+ 2( m-+ m ^ (vi+i,i,k - viJ,k)2 .
Similar expressions can be written down for Gi,j+i/2,k and for IIi,j,k+i/2.
Note that the first order scheme is highly diffusive. This diffusive property is not desirable since it spreads out the original discontinuity and flattens out the peaks in the solution. Therefore usually an explicit finite-difference scheme with flux correction in Chakravarthy-Osher [9-12] form is used. This scheme is first order of approximation in time and third order of approximation in space and is oscillation-free near discontinuities.
3.2. Stability and Boundary Conditions
A characteristic feature of explicit difference schemes is the limits on the time step t, which is governed by the stability criteria. In order to ensure numerical stability, the time step is restricted by a well-known Courant-Friedrich-Lewy stability condition [13], which can be written as
f A^ 4
t = a min< ——- > •-;-— . (10)
г,™\ |A™|J 5 - <p + (1 + y '
where 0 < a ^ 1 is the Courant number, and Ar is the minimum grid spacing; tp and P is free parametres, and last should satisfy to the condition 1 <¡3 < /9max, here
B =
^max — .
1 - y
If tp = 1/3 the scheme there is the third order of approximation [9,12]. With known grid spacing and flow conditions, the time step is evaluated using Eq. (10).
The boundary conditions used usually, are divided into two different types: the land boundary and the open boundary. For the land boundary, the velocity normal to the land is set to zero to represent no flux through the boundary. At the open boundary, it is necessary to solve a boundary Riemann problem [14].
4. Conclusions
The numerical method for solution of the gas-dynamical equations in strict divergent form has been used, in this paper. The three-dimensional numerical code for perfect non-stationary gas-dynamical flows simulation in cylindrical coordinates is constructed. This code is based on the explicit quasimonotonic TVD-scheme. The theory of numerical schemes for homogeneous scalar conservation laws is well established. Total Variation Diminishing (TVD) schemes have proved to be particularly successful at capturing shock waves and discontinuous solutions.
The investigated TVD scheme produce satisfactory results for the selected relevant test cases (see, for example, [15,16]). An adapted one-step third order scheme gives a very good accuracy of the solution in smooth regions and in the proximity of the shock. This motivates the use of TVD-like schemes for inhomogeneous problems, however, although care needs to be taken in the inclusion of the source terms.
References
1. Harten A. High Resolution Schemes for Hyperbolic Conservation Laws //J. Comp. Phys. — 1983. — Vol. 49. — Pp. 357-393.
2. Sweby P. K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws // SIAM J. Numer. Anal. — 1984. — Vol. 21, No 5. — Pp. 995-1011.
3. Roe P. L. The Use of the Riemann Problem in Finite-Difference Scheme // Lect. Notes Phys. — 1980. — Vol. 141. — Pp. 354-359.
4. Roe P. L. Approximate Riemann Solvers, Parameter Vectors and Difference Schemes // J. Comp. Phys. — 1981. — Vol. 43, No 2. — Pp. 357-372.
5. Roe P. L. Characteristic-Based Schemes for the Euler Equations // Ann. Rev. Fluid Mech. — 1986. — Vol. 18. — Pp. 337-365.
6. Филистов Е. А. Модификация 2D численного кода для газодинамических течений в полярных координатах // Вестник РУДН, Серия «Математика. Информатика. Физика». — 2013. — № 2. — С. 150-158. [Filistov E. A. Modification of the Two-Dimensional Numerical Code for Gas-Dynamical Flows in Polar Coordinates. Bulletin of Peoples' Friendship University of Russia, Moscow. — 2013. — No 2. — P. 150-158. — (in russian). ]
7. Landau L. D., Lifshitz E. M. Fluid Mechanics. — Oxford: Pergamon Press, 1979. — P. 433.
8. Tang L. Improved Euler Simulation of Helicopter Vortical Flows: Ph.D. thesis. — University of Maryland, 1998. — P. 227.
9. Chakravarthy S., Osher S. A. A New Class of High Accuracy TVD Schemes for Hyperbolic Conservation Laws // AIAA Pap. — 1985. — No 85-0363. — Pp. 1-11.
10. Einfeldt B. On Godunov-Type Methods for Gas Dynamics // SIAM J. Numer. Anal. — 1988. — Vol. 25, No 2. — Pp. 294-318.
11. Вязников К. В., Тишкин В. Ф, Фаворский А. П. Построение монотонных разностных схем повышенного порядка аппроксимации для систем уравнений гиперболического типа // Математическое моделирование. — 1989. — Т. 1, № 5. —
C. 95-120. [Vyaznikov K.V., Tishkin V.F., Favorskii A. P. Construction of Monotone Difference Schemes of an Increased Order of Approximation for Systems of Equations of Hyperbolic Type // Mat. Model. — 1989. — Vol. 1, No 5. — P. 95120. — (in russian). ]
12. Кузнецов О. А. Численное исследование схемы Роу с модификацией Эйн-фельдта для уравнений газовой динамики // Препринт ИПМ им. М. В. Келдыша РАН. — 1998. — № 43. — С. 44. [Kuznetsov O. A. Numerical Analysis of the Roe's Scheme with Einfeldt's Modification for the Equations of Gasdynamics. — http://library.keldysh.ru/preprint.asp?id=1998-43. — (in russian). ]
13. Courant R., Friedrichs K. O., Lewy H. Uber die partiellen Differen-zengleichungen der mathematischen Physik // Math. Ann. — 1928. — Vol. 100. — P. 32.
14. Approximate Riemann Solvers in FVM for 2D Hydraulic Shock Wave Modeling /
D. H. Zhao, H. W. Shen, J. S. Lai, G. Q. Tabios // J. of Hydraulic Engineering. — 1996. — Vol. 122, No 12. — Pp. 692-702.
15. Filistov E. A. Polygonal Structure of Spiral Galaxies // Astronomy Reports. — 2012. — Vol. 56, No 1. — Pp. 9-15.
16. Filistov E. A. Correlation of a Bar and Spirals at Formation of Galaxies // Astronomical and Astrophysical Transactions. — 2012. — Vol. 27, No 2. — Pp. 215-220.
УДК 519.63+532.5-1/-9 Модификация 3В численного кода для газодинамических течений в цилиндрических координатах
Е. А. Филистов
Кафедра физики Московский государственный строительный университет Ярославское шоссе, д. 26, Москва, Россия, 129337
Цель этой статьи состоит в том, чтобы построить надёжный и точный численный код для решения трёхмерных газодинамических уравнений.
Математическое описание модели представлено системой уравнений непрерывности, движения и энергии. В работе использовано уравнение адиабатического потока невязкого газа.
Для расчёта нестационарных течений идеального газа применён эффективный экономичный метод с использованием полностью консервативной разностной схемы строго дивергентных газодинамических уравнений в эйлеровых переменных в цилиндрических координатах. На основе явной квазимонотонной ТУБ-схемы первого порядка аппроксимации построен ЗБ-численный код для моделирования газового потока. Схема допускает введение ограничителей антидиффузионных потоков, повышающих порядок аппроксимации (до 3-го порядка по пространственным координатам), с минимальной численной диссипацией, и сохраняющих свойство монотонности.
Числовая устойчивость обеспечивается ограничением временного шага известным условием Куранта—Фридрихса—Леви.
Представленная схема отвечает высокому порядку классических схем ТУБ и обладает дополнительнным преимуществом простоты и вычислительной эффективности. Числовые тесты, выполненные автором, показали надежность и эффективность предложенной схемы.
Ключевые слова: газовая динамика, численное моделирование.
