An Euler Solver Based on Locally Adaptive Discrete Velocities

17 Pages·0.18 MB·English

by  Balu Nadiga

An Euler Solver Based on Locally Adaptive Discrete Velocities
5
9
B. T. Nadiga
9
1 Theoretical Division and CNLS, Los Alamos National Lab.,

n Los Alamos, NM 87545
a
J

3
2

0
1
0
1
0
5
9
/
s
a
g B. T. Nadiga
-
p CNLS, MS B258, Los Alamos National Lab.
m
Los Alamos, NM 87545, USA
o
c e-mail: [email protected]
Tel: [phone]
Fax: [phone]
{ 2{
An Euler Solver Based on Locally Adaptive Discrete Velocities
B. T. Nadiga
Theoretical Division and CNLS, Los Alamos National Lab.,
Los Alamos, NM 87545
Abstract
A new discrete-velocity model is presented tosolve the three-dimensional Euler equations. The
velocities in the model are of an adaptive nature|both the origin of the discrete-velocity space
and the magnitudes of the discrete-velocities are dependent on the local (cid:13)ow| and are used in
a (cid:12)nite volume context. The numerical implementation of the model follows the near-equilibrium
(cid:13)ow method of Nadiga and Pullin [1]and results in a scheme which is second order in space (in the
smooth regions and between (cid:12)rst and second order at discontinuities) and second order in time.
(The three-dimensional code is included.) For one choice of the scaling between the magnitude
of the discrete-velocities and the local internal energy of the (cid:13)ow, the method reduces to a (cid:13)ux-
splittingschemebasedoncharacteristics. Asapreliminary exercise,theresultoftheSodshock-tube
simulation is compared to the exact solution.
1. Introduction
A discrete velocity gas is an ensemble of particles with each particle taking on one of a small
(cid:12)nite setofallowablevelocities [2,3]. Further, theinteractionbetween particles isde(cid:12)ned toachieve
the desired macro-behavior of the system, which is usually a set of partial di(cid:11)erential equations.
Such a discretization of the velocity space and de(cid:12)nition of the particle interactions (collisions
or relaxation or more generally redistribution) also form the basis for the lattice gas and lattice
Boltzmann techniques which have been developed over the last eight years [4,5,6 and references
therein]. In spite of the inherently compressible nature of discrete-velocity gases, lattice gases and
lattice-Boltzmann techniques, their applications in (cid:13)uid (cid:13)ow modeling have been restricted to the
incompressible orverylow Machnumber regimes [4,5,6andreferences therein,7,8,9]. In this article,
we present a discrete-velocity gas which for the (cid:12)rst time solves the compressible Euler equations
without introducing any artifacts ofthe velocity discretization overa wide range of Mach numbers.
This is achieved by letting the discrete-velocities of the model adapt to the local (cid:13)ow conditions
[10]. We show also how this scheme can be reduced to a characteristic-based (cid:13)ux-splitting scheme
for the Euler equations [11,12].
In the next section, the model is introduced and shown to reproduce the Euler equations.
Section 3 gives a step by step numerical evolution of the model on the lines of Nadiga and Pullin
[1], and ?? presents a the Sod shock tube simulation [13] with the model. In ?? we show the
reduction of the method to characteristic-based (cid:13)ux-splitting and end with some remarks in ??.
{ 3{
2. The Discrete-Velocity Model
The discrete-velocity model we consider has 27 velocities (qa;a= 0;1;:::;26):
qa = (qax;qay;qaz); (qax;qay;qaz)2 ((cid:0)q;0;q): (1)
p p
At any given point in space, the model thus has four di(cid:11)erent speeds 0;q; 2q, and, 3q. There
p
is one velocity with zero speed, six velocities with speed-q, twelve velocities with speed- 2q, and
p
eight velocities withspeed- 3q. Notethatthis model consistsofthe morefamiliar two-dimensional
nine-velocity model [14,15,16]at three di(cid:11)erent values of velocity in the z-direction ((cid:0)q;0;q),thus
comprising 27 velocities. Collisions between particles in the model are such that they individually
conserve mass, momentum, and energy. In particular, there are two types of collisions which
not only conserve mass, momentum and energy but which also change the speeds of the particles
involved. (The post-collision pair of speeds is not a mere permutation of the pre-collision speeds.)
p
Describing them in the plane of the collision, the (cid:12)rst type involves a speed- 2q particle colliding
with a stationary particle to result in two mutually perpendicular speed-q particles. The second
p
type involves a speed- 3q particle colliding with a stationary particle to result again in a pair of
particles moving mutually perpendicularly, but now one with speed q and the other with speed
p
2q.
2.1 The Stationary Equilibrium Distribution
Thermodynamic equilibrium of the model is de(cid:12)ned as the state of detailed balance of all
possible collisions. Consider the stationary equilibrium of such a gas: since there are no preferred
directions, the particles are identi(cid:12)ed by their speeds and there are thus four variables n0;n1;n2,
and n3, where n0 is the probability thata particle is stationary,n1 is a sixth of the probability that
a particle has a speed q, since there are six speed q velocities in the model, n2 is a twelth of the
p
probability that a particle has speed 2q, and n3 is an eigth of the probability that a particle has
p
speed 3q. The subscripts in n0;:::;n3 are the square of the speed divided by the unit of speed
q. Detailed balancing of collisions results in the equilibrium condition
2
n0n2 =n1;
(2)
n0n3 =n1n2:
The population densities (n0;n1;n2, and n3) are further constrained to satisfy the speci(cid:12)ed hydro-
dynamic quantities mass n and energy ne:
n= n0 +6n1+12n2+8n3;
(3)
2
ne = q (3n1+12n2+12n3):
{ 4{
The stationary (naqa = 0) equilibrium velocity distribution is obtained by the solution of above
four equations (2) and (3):
(cid:18) (cid:19)3
n e
n0 = 3(cid:0)2 2 ;
27 q
(cid:18) (cid:19)2
n e e
n1 = 3(cid:0)2 2 2;
27 q q
(cid:18) (cid:19)(cid:18) (cid:19)2 (4)
n e e
n2 = 3(cid:0)2 2 2 ;
27 q q
(cid:18) (cid:19)3
n e
n3 = 2 :
27 q
2.2 The Euler Equations
TheequationsweseektomodelaretheEulerequationswhichdescribethecompressibleinviscid
(cid:13)ow of a (cid:13)uid (here an ideal monatomic gas). Written in a conservative form [17], and with no
other body forces, they are
@(cid:26)
+r(cid:1)(cid:26)u= 0;
@t
@(cid:26)u
+r(cid:1)(pI+(cid:26)u(cid:10)u) = 0; (5)
@t
@(cid:26)et
+r(cid:1)f(p+(cid:26)et)ug= 0:
@t
2
where (cid:10) represents the binary outer product operator, I is the unit tensor, and et = e+u =2. The
pressure p is given by the ideal gas equation of state for a monatomic gas: p= 2=3(cid:26)e:
2.3 The Locally Adaptive Discrete Velocities
To represent the above equations as adiscrete-velocity gas, we consider the stationary(naqa =
0) 27-velocity gas discussed above under two simple transformations (see Fig.1):
(cid:15) The origin of the discrete-velocity space is translated to u(x;t), where u(x;t) is the tem-
porally and spatially varying velocity (cid:12)eld of the Euler equations.
(cid:15) The unit of discrete velocity, q is determined locally from the speci(cid:12)c internal energy (cid:12)eld
e(x;t) of the Euler equations:
p
q(x;t)= (cid:11)e(x;t) (6):
The scaling factor (cid:11) is a parameter in the model. To insure positivity of the distribution
(4), the restriction on (cid:11) is
2
< (cid:11) < 1: (7)
3
{ 5{
V
3(u,v+q,w)
(u-q,v+q,w)4 2(u+q,v+q,w)
(0,0,w)
U
(u-q,v,w)5 0 1(u+q,v,w)
(u,v,w)
(u-q,v-q,w)6 7(u,v-q,w) 8(u+q,v-q,w)
Fig.1 The origin of the discrete-velocity space is determined by the local (cid:13)ow velocity u=(u;v;w)
of the Euler system and the unit of discrete-velocity q by its speci(cid:12)c internal energy. The
original nine velocities (qax;qay;0)where(qax;qay)2((cid:0)q;0;q))areshownhere aftertheyhave
adapted themselves to the local macroscopic state. The schematic has been drawn on the
qaz = w plane, where w is the z-component of the macroscopic (cid:13)ow velocity used in the
Euler equations. Note that the allowable discrete velocities are now di(cid:11)erent at each point
in space and are di(cid:11)erent at the same point with time.
Figure1 shows the original nine velocities (qax;qay;0;(qax;qay)2((cid:0)q;0;q)) after they have adapted
themselves to the local macroscopic state. The schematic is a projection on to the qaz = w plane,
where w is the macroscopic velocity in z-direction used in the Euler equations (5). With this kind
of adaptation, the allowable velocities in the model are now ca(x;t) = qa(x;t)+u(x;t), i.e., the
allowable discrete-velocities are completely di(cid:11)erent at each point in space, and are di(cid:11)erent at the
same point in space at di(cid:11)erent times.
{ 6{
2.4 The Equivalence of the model to the Euler Equations
The density and energy of the discrete-velocity gas are set equal to that in (5) and again note
that naqa has been assumed 0. The equivalence of the discrete-velocity gas system to the system
of Euler equations might now be obvious. However, for the sake of completeness, we explain
the equivalence. Considering the evolution of the discrete-velocity gas, the (model) Boltzmann
equations [3,15]|a statement of the conservation of the number of particles with a particular
discrete-velocity|are
@na
+ca(cid:1)rna = Qa(n;n); a= 0;:::;26; (8)
@t
where Qa is the nonlinear collision operatorand the left hand side represents streaming of particles
with velocity ca. The zeroth, (cid:12)rst, and second order velocity moments of (8) give respectively,
noting that the moments of the collision terms on the right-hand side vanish owing to the mass,
momentum, and energy conserving nature of each collision,
@n
+r(cid:1)naca = 0;
@t
@naca
+r(cid:1)naca(cid:10)ca = 0; (9)
@t
@nac22a c2a
+r(cid:1)na ca = 0;
@t 2
where the overbar denotes averaging with respect to the discrete-velocities. From the translation
of the origin of the discrete-velocity gas,
naca = na(u+qa) = (cid:26)u; since naqa = 0:
naca (cid:10)ca =na(u+qa)(cid:10)(u+qa)
(10)
=nu(cid:10)u+naqa(cid:10)qa = pI+(cid:26)u(cid:10)u:
2 2 2
ca u +qa
na ca =na( +u(cid:1)qa)(u+qa)
2 2
2 2
qa u
=na u+(cid:26) u+na(qa(cid:1)u)qa
2 2
=(cid:26)etu+na(qa(cid:10)qa)u= (p+(cid:26)et)u:
The equivalence is thus complete. Note that in the above equation, since naqa = 0, the particles
are distributed symmetrically with respect to qa and so naqa(cid:10)qa is isotropic and reduces to the
pressure p.
For convenience, the Euler equations (5) and the moment equations (9) may be rewritten as
@f
+r(cid:1)G= 0
@t
!
(cid:26)
with f = (cid:26)u ; (11)
(cid:26)et
" #
2
ca
G = [(cid:26)u;pI+(cid:26)u(cid:10)u;(p+(cid:26)et)u]= naca;naca(cid:10)ca;na ca :
2
{ 7{
3. The Numerical Technique
Hyperbolic systems ofconservation laws like the Euler equations (5)in the present case, admit
weak solutions in the form of shocks. At these shocks, the gradients of the primary quantities
are in(cid:12)nite, and obviously this cannot be represented correctly in a shock-capturing numerical
scheme. Shock-capturing schemes are ones in which all grid points are treated in exactly the
same fashion irrespective of whether they are inside a shock or outside, and as opposed to shock-
tracking schemes which keep track of where the shocks are and treat grid points located in shocks
di(cid:11)erently from the others. While each have their advantages and disadvantages, used on present
daysupercomputers, shock-capturing methodsareclearly preferable because oftheirhomogenityof
computation. Numerically capturingshocksinnondissipativesystemsliketheEulerequationsposes
the problem of dispersive ripples: a phenomenon in which as a wavesteepens, energy (cid:13)ows into the
smaller scales and since there is no dissipation, accumulates in the smallest allowed wavelengths|
those of the grid-spacing. Thus with the formation of any shock, energy piles up in grid scale
oscillations and swamps out scales of interest. A way out of this problem is to dissipate energy
at the smallest length scales|the grid-scales|either explicitly or implicitly. Explicit arti(cid:12)cial
viscosity in general needs adjustments depending on the problem, while physics-based (as opposed
to the inviscid idealization represented by the Euler equations) implicit arti(cid:12)cial (from the point of
view of the Euler equations) viscosity is robust and results in numerical shock widths of the same
order as the actual viscous (Navier-Stokes in the present case) shock widths.
We use exactly such an implicit arti(cid:12)cial viscosity technique (also called a kinetic numerical
schemeowingtothephysical kinetic basisofthescheme)which wasdeveloped inNadigaand Pullin
[1] for the discrete-velocity framework. A brief description of its usage here follows: For simplicity,
consider the computational domain divided into uniform cubical cells, at the centroids of which
are stored the cell-averaged values of ((cid:26);(cid:26)u, and (cid:26)et). The evolution at each centroid proceeds as
follows:
Step 1: Calculate the local unit of discrete-velocity using (6) and (7).
Step 2: Using (4) and the macroscopic variables of density (cid:26) and speci(cid:12)c internal energy e,
calculate the population densities of the four speeds n0;n1;n2, and n3. Note that in light of
restriction (7)and owing to the fact that the stationarydistribution function was calculated based
on detailed balancing, the above four population densities are necessarily positive.
+ (cid:0)
Step 3: Calculate the split (cid:13)uxes G and G at the centroids using the de(cid:12)nitions
0 1
cax>0 cax>0 2 cax>0 cax>0 cax>0c2a
na cax na cax na caycax na cazcax na 2 cax
B C
G+ = B cay>0 cay>0 cay>0 2 cay>0 cay>0c2a C;
@na cay na caxcay na cay na cazcay na 2 cayA
caz>0 caz>0 caz>0 caz>0 2 caz>0c2a
na caz na caxcaz na caycaz na caz na 2 caz
0 1
cax<0 cax<0 2 cax<0 cax<0 cax<0c2a
na cax na cax na caycax na cazcax na 2 cax
B C
G(cid:0) = B cay<0 cay<0 cay<0 2 cay<0 cay<0c2a C; (12)
@na cay na caxcay na cay na cazcay na 2 cayA
caz<0 caz<0 caz<0 caz<0 2 caz<0c2a
na caz na caxcaz na caycaz na caz na 2 caz
{ 8{
cax>0
where na cax is the averagenacax takenonly overthediscrete velocities ca which haveapositive
x-component cax, etc..... For example, when u > 0, and w > 0, but v < 0,
0 1
2 2 cax>0c2a
uP +(u+q)Q u P +(u+q) Q uvP +(u+q)vQ uwP +(u+q)wQ na 2 cax
B C
G+ = B 2 cay>0c2a C;
@ (v+q)Q (v+q)uQ (v+q) Q (v+q)wQ na 2 cayA
2 2 caz>0c2a
wP +(w+q)Q wuP +(w+q)uQ wvP +(w+q)vQ w P +(w+q) Q na 2 caz
(13)
0 1
2 cax<0c2a
(u(cid:0)q)Q (u(cid:0)q) Q (u(cid:0)q)vQ (u(cid:0)q)wQ na 2 cax
B C
G(cid:0) = B 2 2 cay<0c2a C;
@vP +(v(cid:0)q)Q vuP +(v(cid:0)q)uQ v P +(v(cid:0)q) Q vwP +(v(cid:0)q)wQ na 2 cayA
2 caz<0c2a
(w(cid:0)q)Q (w(cid:0)q)uQ (w(cid:0)q)vQ (w(cid:0)q) Q na 2 caz
where P = (n0 +4n1 +4n2) is the density of particles in say the cax = u plane and Q = (n1 +
4n2+4n3) is the density of particles in say the cax = u+q or cax = u(cid:0)q plane. In the above two
cax>0c2a caz<0c2a
equations, na 2 cax, :::, na 2 caz have been left as such for compactness of notation.
Step 4: Assuming a linear distribution of the (cid:13)uxes within the cells in the direction under consid-
+ (cid:0)
eration, interpolate the split (cid:13)uxes G and G to the cell boundaries:
+ + + + + 1
(cid:15) interpolate G , G , G , G , and G to (i+ ;j;k),
11 12 13 14 15 2
(cid:0) (cid:0) (cid:0) (cid:0) (cid:0) 1
(cid:15) interpolate G , G , G , G , and G to(i;j(cid:0) ;k),and so on. (Note thatthe (cid:12)rstsub-
21 22 23 24 25 2
script of G corresponds to the coordinate direction and the second to one of the conserved
quantity.)
and apply the minmod limiter to the interpolated (cid:13)uxes:
+ 1 + 1 + +
G11(i+ ;j;k)= G11(i;j;k)+ minmod((cid:1)bckG11(i;j;k);(cid:1)fwdG11(i;j;k));
2 2
(cid:0) 1 (cid:0) 1 (cid:0) (cid:0)
G21(i;j(cid:0) ;k)= G21(i;j;k)(cid:0) minmod((cid:1)bckG21(i;j;k);(cid:1)fwdG21(i;j;k)); (14)
2 2
and so on, and where
+ + +
(cid:1)fwdG11(i;j;k)= G11(i+1;j;k)(cid:0)G11(i;j;k)
+
is the (cid:12)rst forward di(cid:11)erence of G in the x-direction at the centroid of cell (i;j;k).
11
(cid:0) (cid:0) (cid:0)
(cid:1)bckG21(i;j;k)= G21(i;j;k)(cid:0)G21(i;j(cid:0)1;k)
+
is the (cid:12)rst backward di(cid:11)erence of G in the y-direction at the centroid of cell (i;j;k) and so on.
21
And minmod is the one dimensional total-variation-diminishing operator as discussed in [18,19]:
(cid:26)
0 if sgn(p) 6= sgn(q)
minmod(p;q)= sgn(p) (15)
minfjpj;jqjg if sgn(p) = sgn(q),
{ 9{
with sgn(p) being the sign of p and jpj the absolute value of p.
Step 5: Calculate the (cid:13)uxes at the cell boundaries:
1 + 1 (cid:0) 1
G(i+ )= G (i+ )+G (i+ )
2 2 2
1 1 1 1
where (i+ )=(i+ ;j;k), (i;j+ ;k), and (i;j;k+ ) in turn.
2 2 2 2
Step 6: Advance the primary variables by one half of the time step to get the midpoint values:
ft0+(cid:1)t=2 =ft0 + (cid:1)t (cid:0)r(cid:1)Gt0(cid:1); (16)
2
!
G11 G11(i+ 12;j;k)(cid:0)G11(i(cid:0) 12;j;k)
where e.g. r(cid:1) G21 = +
(cid:1)x
G31
1 1 1 1
G21(i;j+ 2;k)(cid:0)G21(i;j(cid:0) 2;k) G31(i;j;k+ 2)(cid:0)G31(i;j;k(cid:0) 2)
+ + :
(cid:1)y (cid:1)z
Step 7: Repeat steps 1-4 using the midpoint values ft0+(cid:1)t=2 to obtain the (cid:13)uxes at the midpoint
Gt0+(cid:1)t=2, and take the full time step to obtain the new time level values ft1:
(cid:16) (cid:17)
ft1 = ft0 +(cid:1)t r(cid:1)Gt0+(cid:1)t=2 : (17)
Intheaboveprocedure the(cid:13)uxes wereinterpolated andlimited, instead, theprimaryquantities
f could be interpolated and limited. We have done both and the behavior of the twoare essentially
identical. The time step in the method is limited by the CFL stability criterion
(cid:18) (cid:19)
(cid:1)t
max jU (cid:6)qj (cid:20) 1; U 2 (u;v;w) and (cid:1)X 2 ((cid:1)x;(cid:1)y;(cid:1)z): (18)
(cid:1)X
The code (in Connection Machine Fortran) described by the above step-by-step procedure is in-
cluded in the Appendix. Notwithstandingthe somewhatcomplicated procedural description above,
the simplicity of the resulting code is evident.
4. A Numerical Example
The code presented in the Appendix, a Connection Machine Fortran implementation of the
second order (in space and time) scheme described in the previous section, was used to calculate
theSodshock-tubeproblem. TheSod testcasehastheinitial conditions ((cid:26)l = 1;ul = 0;pl = 1)and
((cid:26)r = 0:125;ur = 0;pr =0:1) corresponding to an initial pressure ratio of 10 and a density ratio of
8. Subscript ldenotesthelefthalfandsubscript rdenotestherighthalfattime0. Foramonatomic
gas,thisreduces totheinitial conditions ((cid:26)l = 1;ul = 0;el = 1:5)and((cid:26)r = 0:125;ur = 0;er = 1:2).
Only 128pointswereusedforthissimulation, atatimestepcorresponding toaCFLnumberof0.69.
(The (cid:11) in (6) was set at 10=9.) In Fig.2, the exact density, speci(cid:12)c internal energy, velocity and
pressure pro(cid:12)les (solid lines) are compared to the corresponding pro(cid:12)les (open diamonds) obtained
fromthediscrete velocity simulation. The agreementis good: The shock is typically 3-4cell-widths
and the edges of the rarefaction are not too badly rounded. The spreading of the contact surface
is however substantial as with other (cid:13)ux-splitting schemes.
{ 10{
1.2 2.0
1.0 Density
1.5
0.8 Rarefaction
0.6 1.0
Contact
0.4 Surface Specific Internal Energy
0.5
0.2
Shock
0.0 0.0
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6
X X
1.0 1.2
0.8 1.0 Pressure
0.6 0.8
0.4 0.6
0.2 Velocity 0.4
-0.0 0.2 Inital High Press. Initial Low Press.
and Density and Density
Region Region
-0.2 0.0
-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6
X X
Fig.2 The Sod shock-tube problem: At a dimensionless time of 0.288 after the diaphragm burst
(the diaphragm is initially at x=0),the inital discontinuity in pressure and density (pressure
ratioof10andadensity ratioof8)isresolved intoaright-going shockwavelocated atabout
x=0.5 , a left-going rarefaction centered about x=-0.3, and a right-going contact surface at
about x=0.2. The solid line is the exact solution and the diamonds are the result of the
simulation using the locally adaptive discrete velocities.
5. The recovery of Characteristic-based (cid:13)ux-splitting
From the point of view of (cid:12)nite di(cid:11)erencing, high-order shock-capturing techniques used for