CC BY
4
DOI: 10.17516/1997-1397-2022-15-4-500-509 УДК 519.688
Compact MacCormac-type Schemes Applied for Atmospheric Escape Problem
Kseniya D. Gorbunova*
Siberian Federal University Krasnoyarsk, Russian Federation Institute of Computational Modelling SB RAS Krasnoyarsk, Russian Federation
Nikolai V. ErkaeV
Institute of Computational Modelling SB RAS Krasnoyarsk, Russian Federation
Received 27.02.2022, received in revised form 04.06.2022, accepted 11.06.2022
Abstract. This article discusses the advantages of the compact MacCormack-type scheme over its original version when modeling the hydrodynamic outflow of planetary atmospheres caused by the absorption of ultraviolet radiation from the Sun in the upper atmosphere. The results of calculations with different parameters of the problem are presented.
Keywords: MacCormack-type method, compact difference scheme , planetary atmospheres, intensity of atmospheric loss.
Citation: K. D. Gorbunova, N. V. Erkaev, Compact MacCormac-type Schemes Applied for Atmospheric Escape Problem, J. Sib. Fed. Univ. Math. Phys., 2022, 15(4), 500-509. DOI: 10.17516/1997-1397-2022-15-4-500-509.
Introduction
In computational hydrodynamics, the MacCormack method [1,2] is a widely used discretization scheme for numerical solution of hyperbolic partial differential equations.This 2-nd order finite difference approximation method was first introduced by Robert W. MacCormack in 1969. D.Gottlieb and E.Turkel proposed a modification [3] of the MacCormack scheme, known as scheme 2-4. It has a 2-nd order of accuracy in time and 4-th order in space. This scheme has been successfully used for solving a wide range of gas-dynamic and acoustic problems. Further improvements of the MacCormack method in literature [4-17] were aimed at increasing the approximation order on the base of explicit or implicit stencils for derivative approximations. The implicit stencils are used in the so called compact schemes, which have a particular advantage that do not require a significant extension of the stencil size with an increase in the order of approximation. A new series of compact MacCormack-type schemes with one-sided implicit templates is presented in [14], where explicit boundary stencils are given, boundary conditions are described and the results of test calculations are shown. The authors conclude that the compact methodology provides a significant performance advantage over previous explicit MacCormack-type schemes.
Previously, the classical MacCormack scheme was used in [19, 20] to simulate the hydrodynamic outflow of planetary atmospheres caused by the absorption of extreme ultraviolet (EUV) radiation from the Sun in the upper atmosphere. This process plays a very important role in
* [email protected] [email protected] © Siberian Federal University. All rights reserved
the evolution of planets. However, application of the original MacCormack scheme for this task was limited to a rather narrow range of parameters for which the scheme worked steadily. This is because of large gradients of density> pressure and especially temperature in vicinity of planets, where the original MacCormack scheme produces growing non-physical oscillations. In this regard, our paper considers the advantages of using a modification of the MacCormack scheme of the 4th order with a compact approximation.
1. Description of the compact schemes
We consider two modified MacCormack-type schemes (4/2 and 4/4) described by R. Hixon and E. Turkel [14]. The peculiarity of this method is that the derivative operator splits into two finite difference operators
d = DF+DB, (1)
2 V '
where DF, DB are the forward and backward derivative operators, respectively. In 4/2 compact scheme, the spatial derivatives of a function f (x) are determined by the following equations
ADB-i + (1 - A)DB = ---(fi - fi-i),
-X (2) ADF+i + (1 - A)DF = Axx(fi+1 - fi),
were A = 2 ± ^^ and the sign selection depends on the choice of the boundary stencils.
We choose sign "minus" according to [14].
Taylor series expansion of the forward and backward operators gives:
DF = f< + f(II) - f(IV) + O(Ax)4,
DB = f> - ^2f(II) + (MV2f(IV) + a(Ax)4.
Hence this scheme provides the 4th order accuracy for derivative approximation (2). The 4/4 compact scheme implies larger stencil size as follows:
1 d—+2 db=^fi+i+3f - 5 /.-a
3D+ + 3DF = ÂX(6fi+1 - 3fi - 6fi)-The Taylor series expansion of these operators also shows a 4th order of accuracy:
DB = f> + ^ f(IV ) + O(Ax)4,
(3)
(4)
DF = IL - f(IV) + O(Ax)4.
(5)
We consider unsteady conservative hyperbolic system of equations written in the vector form as follows
Ut + {G(U )}a = Q(U), (6)
where U is the vector of unknown functions, G(U) is the given vector function of U, and Q is the source function depending generally on U, x. To solve this equation, we apply the four-stage
Runge-Kutta time marching method (RK4) using the compact MacCormack schemes
h(1) = -AtDF [G(Un)] + AtQ(Un)
h(2) = -AtDB [G(Un + a2h(1))] + AtQ(Un + a2
h(3) = -AtDF [G(Un + a3h(2))] + AtQ(Un + a3
h(2)),
h(4) = -AtDB [G(Un + a4 h(3))] + AtQ(Un + a4 h(3)).
(7)
Here, At is the time step, a2 = ^= ^,a4 = 1 are the coefficients for the RK4. DF and DB
are the forward and backward spatial derivative approximations determined by equations (2) or (4). The function on the next time step is calculated as follows
Un+1 = Un + 5ih(1) + Ô2h(2) + S3h(3) + Ô4h(4),
(8)
were ¿i = 1, S2 = 1, S3 = 1, S4 = 1 are also the coefficients for the RK4 method. 6 3 3 6 R. Hixon and E. Turkel [14] noted that the 4/2 compact scheme using the RK4 time marching
scheme has 4-th order of accuracy in space, and the time accuracy would depend on nonlinear
accuracy of the time advancement algorithm. Earlier R. Hixon [12] mentioned that the RK4
method would have a 4-th order of accuracy in linear problems, and in nonlinear problems an
order less than the 4-th.
Each compact scheme requires an appropriate choice of the boundary stencils. In particular,
explicit five-point boundary stencils for the 4/2 compact scheme derived in [14] are given by
DF
43 -
DB .
i mtn
' 25 17
—12+ Т2Т3У
fi + ( 4 -5
25
6л/з) fi+3 ( 4 12%/3
fi+4
6V3/ ]1
fi+1 - 3 -
(з - з23)
Ax'
-I — + J7 , ' 12 12 V3J
)
+ (3 + 6^
fi + ( 4 + 5
13 ^ fi+3 -( i + ^73
df
i max
-(3+
17
25
12 + 12V3J
fi
4 +
fi+4
25
25 N
673) ~ 1
Ax'
fi+1 - 3 +
И)
13
DB
i max
6V3
6V2J fi-3 + (i + йЫ fi-4
V 1 2 12^7 V 6^
fi-1 + 3 +
1
Ax'
(3 + 3^)
fi-1 + 3 -
(4 - ^ 44 -
1
5 f4
4 12V3J
1
Ax'
(3 - 7)
fi-2-
fi-2
(9)
Also the boundary stencils for the 4/4 compact scheme are given by
D
F
i min
DB
DF
i max
D
B
i max
- Щ f-+
37 ) fi+1
- ( "6" j fi+2 + ( -9- ) fi+3 - ( 18 ) fi+4
2 ä)
- S )fi+( ? ) fi+1 - ( ? ) fi+2+( y ) fi+3 - ( ) ) fi+4
S)fi - 2 35 )fi-1+2 17 )fi-2 - 2 ï)fi-3+2 9) fi-4.
5
1
Ax' 1
Ax'
f)fi - (?)fi-1 + (f)fi-2 - (f)fi-3 + (18)fi4
1
Ax' 1
Ax'
R.Hixon noted [21] that the numerical realization of the boundary stencils for compact scheme might have a great effect on stability and accuracy of the scheme, since the error can propagate from the boundaries into the computational domain.
2. Statement the problem of atmospheric loss
We consider a 1-D unsteady hydrodynamic equations in spherical coordinates which were applied in [19] for simulation of upper atmosphere of planets exposed by EUV solar radiation. Numerical solution of this problem was obtained on a base of the original two-stage MacCormack scheme, which is of 2-nd order accuracy in both space and time,
d[pr2] + d[pvr2 ]
dt
dr
0,
d[pvr2] + d[r2(pv2 + P )]
dt
dr
+ M + d\v( ^ + Y—1 )
dt dr
-pß + 2Pr
pv2 yP
Y + Y—1,
(11)
—pv ß + qr2
Here, t and r are the dimensionless time and radial distance, respectively; 3 is the constant parameter equal to the ratio of the gravitational and thermal energies of a gas particle at the lower boundary of the computational domain which is near the planetary surface. In the astrophysical literature, 3 is called as "Jeans parameter", characterizing ability of particles to overcome the planetary gravitation. The gas parameters p,v,P,j, q are the dimensionless mass density, radial velocity, pressure, adiabatic index and normalized EUV heating rate, respectively. Here the radial distance, density and pressure are normalized to their values at the lower boundary R0, p0 and P0, respectively; the velocity is normalized to \JP0/p0
According to [19], the heating rate can be approximated by the analytical function depending on the density distribution in the atmosphere
q(r) = 0.767Bp(r) exp - 0.6as(s3 + 1)p(s)ds
0-
(12)
Here a = ap0R0/m, where a is the cross section of EUV absorbtion by the atmospheric atom, and m is the mass of atom; B = 0.096 IEuv, where Ieuv is the intensity of the EUV radiation at the upper boundary far away from the planetary surface.
Initial conditions with isothermal density distribution are given by
pinit = exp
ß( - - 1 r
(13)
Vinit = 0.1(r - 1),
Pinit/pinit 1,
where pinit,vinit, Pinit are the mass density, radial velocity and pressure at the initial time. At the lower boundary (r=1) we set conditions for the density and pressure,
P =1, P =1. (14)
The velocity at this boundary is determined by the solution.
The upper boundary (r = rm) has to be chosen at sufficiently large distance rm, where the flow becomes supersonic. In such case, physical quantities at this boundary are determined by the
2
d
r
solution in the interior region. In our numerical scheme, we calculate the upper boundary velocity, density and pressure at each time step, using relations along three outgoing characteristics:
dv dP dP n dp —pc ± — = 0, ---c2^- =0,
dt± dt± dto dto /-, r\
dr dr (15)
—— = v ± c, —— = v,
dt± dto
where d/dt+— = d/dt + (v ± c)d/dr, and d/dt0 = d/dt + vd/dr. Intensity of atmospheric loss can be written as
L = p(rm)v(rm)r2m, (16)
where p(rm) and v(rm) are the mass density and radial velocity at the upper boundary of the computational domain, and rm is the radial distance to the boundary. When solving the problem, we are interested in evolution of the loss rate L, as well as radial profiles of the velocity (v), pressure (P), density (p) and temperature (T = P/p) for various parameters ¡3. During time relaxation, radial distributions of the atmospheric parameters evolve to a stationary flow structure with supersonic velocity at the upper boundary. An existence of such supersonic outflow regime is directly related to the EUV heating source.
3. Calculation algorithm
For computational convenience, we redefine the variables and introduce new functions:
p' = pr2, P' = Pr2 = p'T,
Sv P' S = p'v, G = Sv + P', E =— +
W _ Sv2 + YvP'
2 ^ y - 1' (17)
2 y - 1
Using new notations (17), we transform hydrodynamic equations (11) to those with respect to three unknown functions p',E,W:
dp' dS _ 0 dt dr '
dS + dG(E,S,p') _ + 2(y - 1) (E - S2A
dt dr r2 r V 2p' J '
dE + dW(E,Sp') _ R + + H -Trr +--«-_ -S^2 + q +
dt dr r2
G(E, S, p') _ (y - 1)E S2,
2p
(18)
S2
yE - (Y - 1)S1 , 2p'
rm 1 o' (s)
q'(r) _ 0.767Bp'(r)exp(- 0.6a(s3 + ds),
S
W(E, S, p') _ -p
where H is the additional stabilization term specified below. As mentioned by [14], stability of the scheme is provided by the Courant-Friedrichs-Lewy (CFL) condition, which yields the upper limit for time step At
Ar
At _ K-^--, c _JYp, K< 1. (19)
max(|v| + c) y p
To solve equations (18), we apply two compact MacCormack-type schemes RK4-4/2 and RK4-4/4. The first one is described above by formulas (2), (7) with boundary stencils (9), and the second one is related to (4), (7) with boundary stencil (10).
In cases of large parameters 3, the heating source term has typically a very steep enhancement for increasing radial distance just prior its maximum point. This leads to the corresponding very large gradient of the temperature before its maximum. This temperature behaviour might cause growing oscillations of solution in the region of enhanced gradient. These oscillations can be suppressed by a stabilizing additive inserted in the right hand side of the energy equation. We use the stabilizing term defined as
-=< (p'dT) (20)
where A is a constant coefficient, and h is the grid size.
4. Results
In our simulations, we apply the basic parameters similar to those used in [19]
5
Ieuv = 10, r e [1, 20], y = 3, a = 1.593e4.
The criterion for choosing the upper boundary is the supersonic velocity condition. This condition will be satisfied if the upper boundary is set at the distance rm = 20. The calculations were performed on five uniform grids (3000, 4000, 8000, 16000, 20000 points) for different parameters 3 (30, 40, 50, 60) and the CFL number K = 0.5.
For the stabilizing term we use the following grid approximation
nn = Xh
T+l _ Tn pn + pn+l _ Tin _ T— pn + PÏ-!
_Ti+1 _ Ti 2 Ti _ Ti-1 2
(21)
where A = 2, n is the number of time step, and i is the number of spatial grid point.
In each case, calculations were performed during sufficiently long time to obtain a stationary solution with constant value of the loss rate L. A stationary state of the solution is characterized also by the integral energy conservation
qp'dr = L3 + L (-(rm)2 --(ro)2 + ^(T(rm) - T(r0))) (22)
We have to note about restrictions of usage initial conditions (13), which can provide a stable calculations only for moderate 3 parameters (3 < 50). Taking this into account, we used results for smaller 3 as initial conditions when calculated for larger ¡3. For example, the final results for 3 = 50 were used as initial conditions to calculate case 3 = 60, and so on.
Fig. 1 shows variations of the loss rate L as function of time for different values of the parameter 3 in the range from 30 to 60. This parameter would be very large in cases of giant planets, heavy gas particles or low temperature of the atmosphere. One can see that the stationary loss rate is larger, when the beta parameter is smaller.
The computational complexity of the problem of the hydrodynamic atmospheric outflow is associated with a certain radial dependence of the heating rate, the gradient of which increases rather sharply as we approach the planet. This yields a consequent sharp increase in temperature, the maximum of which is substantially higher for larger values of the main parameter 3, as can be seen in Fig. 2. Radial gradients of the density, pressure and temperature become as stronger as larger the parameter 3.
Time
Fig. 1. Time dependence of the atmospheric loss rate for different parameters 3
Fig. 2. The atmospheric temperature as function of the radial distance
Fig. 3 shows the radial velocity as the increasing function of the radial distance. Contrary to temperature, the velocity maximum, reached at the upper boundary, somewhat decreases with an increase in the parameter 3.
Fig. 4 at the top panel (a) demonstrates the stationary velocity maximum reached at the upper boundary in dependence on the grid step. In addition, at the bottom panel (b) it shows
2.4
>
o
0
2
4
6
8
10
12
14
16
18
20
Fig. 3. The radial velocity profiles
the approximation errors 6 (for the grid step h = 0.004) corresponding to the original MacCormac scheme and the modified compact scheme described above.,draft
For the exact stationary solution, expression (23) must vanish. Fig. 4 indicates clearly that the compact MacCormac scheme provides much better accuracy compared to the original scheme. Moreover, applying the original MacCormac scheme for the same problem we get numerical instability for 3 > 40, while the compact MacCormac-type scheme used in our work allows us to obtain stable solutions in wide range of the 3 parameter.
Conclusions
We analyzed an application of the compact MacCormack-type scheme for solving the hy-drodynamic problem of atmospheric escape due to absorption of the extreme solar radiation, which is important for planetary evolution models. The gradients of the physical quantities (temperature, density, pressure) increase crucially with increasing Jeans parameter ¡3, which makes computational difficulties. In particular, a stabilizing additive was introduced to suppress possible oscillations appeared in the region of the very steep rise in temperature in front of its peak. The results were compared with those presented in [19,20], which were obtained using the original MacCormack scheme. Our study reveals a significant advantage of the compact MacCormack scheme RK4-4/2, which allows one to solve the problem considered above for large values of the Jeans parameter (3 > 40), while the scheme used in [19,20] becomes unstable in this range of the parameter.
In this work we used the compact scheme on a uniform spatial grid, and an extension of the method on non-uniform grids would make further substantial improvement.
(23)
a)
2.257
13
. £ 2.256 2.255
1 1.5 2 2.5 3
b) h x 10-3
0.2 0.1 ^ 0 -0.1 -0.2
6 8 10 12 14 16 18
r
Fig. 4. a) The stationary velocity maximum as function of the grid step; b) Approximation error 5 defined by formula (23) in cases of the original MacCormac and the modified compact MacCormac-type schemes, ¡3 =30
This work is supported by the Krasnoyarsk Mathematical Center and financed by the Ministry of Science and Higher Education of the Russian Federation in the framework of the establishment and development of regional Centers for Mathematics Research and Education (Agreement no. 075-02-2022-873).
References
[1] R.W.MacCormack, The effect of viscosity in hypervelocity impact cratering, AIAA Paper 69-354, 1969.
[2] R.W.MacCormack, Lecture notes in physics, Proceedings of the Second International Conference on Numerical Methods in Fluid Dynamics, 1971.
[3] D.Gottlieb, E.Turkel, Dissipative Two-Four Method for Time Dependent Problems, Mathematics of Computation, 30(136)(1976).
[4] S.K.Lele, Compact Finite Difference Schemes with Spectral-like Resolution, Journal of Computational Physics, 103(1992), no. 1.
[5] Z.Haras, S.Ta'asan, Finite difference schemes for long-time integration, Journal of Computational Physics, 114(1994), no. 2.
[6] J.W.Kim, D.J.Lee, Optimized compact finite difference schemes with maximum resolution, AIAA journal, 34(1996), no. 5, 887-893. DOI: 10.2514/3.13164
[7] S.T.Yu, K.C.Hsieh, Y.L.P.Tsai, Simulating waves in flows by Runge-Kutta and compact difference schemes, AIAA journal, 33(1995), no. 3. DOI: 10.2514/3.12470
[8] R.Agarwal, K.Huh, A dispersion-relation preserving fourth-order compact time-domain/frequency-domain finite-volume method for computational acoustics, AIAA Paper 96-0277, 1996. DOI: 10.2514/6.1996-277
[9] C.K.W.Tam, J.C.Webb, Dispersion-relation-preserving finite difference schemes for computational acoustics, Journal of computational physics, 107(1993), no. 2.
- ---Original scheme -Compact scheme 1 ». 0 / / ij! ^
" " »' u ' Ai iji^-i 1 ' 1
10] D.Zingg, H.Lomax, H.Jurgens, An optimized finite-difference scheme for wave propagation problems, AIAA Paper 93-0459, 1993.
11] C.A.Kennedy, M.H.Carpenter, Several new numerical methods for compressible shear-layer simulations, Applied Numerical Mathematics, 14(1994), no. 4.
12] R.Hixon, On Increasing the Accuracy of MacCormack Schemes for Aeroacoustic Applications, AIAA Paper 97-1586, 1997.
13] R.Hixon, E.Turkel High-accuracy compact MacCormack-type schemes for computational aeroacoustics, AIAA Paper 98-0365, 1998.
14] R.Hixon, E.Turkel, Compact Implicit MacCormack-Type Schemes with High Accuracy, Journal of Computational Physics, 158(2000).
15] S.A.Jordan, Optimization, resolution and application of composite compact finite difference templates, Applied numerical mathematics 61(2011), no. 1. DOI: 10.1016/J.APNUM.2010.08.008
16] S.Yazdani, D.R.Hixon, N.Mansouri, Hybrid MacCormack-type Schemes for Computational Aeroacoustics, AIAA Paper 2015-3134, 2015. DOI: 10.2514/6.2015-3134
17] A.Rona, I.Spisso, et al., Optimised prefactored compact schemes for linear wave propagation phenomena, Journal of Computational Physics, 138(2017). DOI: 10.1016/j.jcp.2016.10.014
18] R.Ciegis, O.Suboc, High order compact finite difference schemes on nonuniform grids, Applied Numerical Mathematics, 132(2018). DOI: 10.1016/j.apnum.2018.06.003
19] N.V.Erkaev, H.Lammer, et al,. XUV exposed non-hydrostatic hydrogen-rich upper atmospheres of terrestrial planets. Part I: Atmospheric expansion and thermal escape, Astrobiol-ogy Journal,, 13(2013), no. 11, 1011-1029. DOI: 10.1089/ast.2012.0957
20] N.V.Erkaev, H.Lammer, et al., EUV-driven mass-loss of protoplanetary cores with hydrogen-dominated atmospheres: the influences of ionization and orbital distance, Monthly Notices of the Royal Astronomical Society, 460(2016), no. 2. DOI: 10.1093/mnras/stw935
21] R.Hixon, A New Class of Compact Schemes, 36th AIAA Aerospace Sciences Meeting and Exhibit, 1998.
Применение компактных разностных схем для задач гидродинамического истечения атмосфер планет
Ксения Д.Горбунова
Сибирский федеральный университет Красноярск, Российская Федерация Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Николай В. Еркаев
Институт вычислительного моделирования СО РАН Красноярск, Российская Федерация
Аннотация. В этой статье рассматриваются преимущества компактной схемы типа МакКормака над ее оригинальной версией при моделировании гидродинамического истечения атмосфер планет, обусловленного поглощением ультрафиолетового излучения Солнца в верхних слоях атмосферы. Приведены результаты расчетов с разными параметрами задачи.
Ключевые слова: схема типа МакКормака, компактное дифференцирование, атмосфера планеты, интенсивность потери атмосферы.
