=Paper=
{{Paper
|id=Vol-1839/MIT2016-p32
|storemode=property
|title= About verification of calculation methods of the shock waves
|pdfUrl=https://ceur-ws.org/Vol-1839/MIT2016-p32.pdf
|volume=Vol-1839
|authors=Valentin Kuropatenko,Elena Shestakovskaya
}}
== About verification of calculation methods of the shock waves==
https://ceur-ws.org/Vol-1839/MIT2016-p32.pdf
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
About Verification of Calculation Methods of the Shock
Waves
Valentin Kuropatenko1,2 and Elena Shestakovskaya2
1
Russian Federal Nuclear Center-Zababakhin All–Russia Research Institute of Technical
Physics,
456770, Snezhinsk, Russian Federation
2
South Ural State University(National Research University),
454080, Chelyabinsk, Russian Federation
v.f.kuropatenko@yandex.ru,leshest@list.ru
Abstract. Mathematical modeling is now a key tool of research into dynamic
processes in continuum mechanics. Each particular problem is solved with al-
ready existing or newly developed models and methods, whose properties are
determined from a priori study into stability, approximation, monotonicity etc
within linear approaches. The accuracy of difference schemes is mainly evalu-
ated through comparison between calculated results and reference solutions. The
paper discusses some problems which have analytical solutions. These are shock
convergence, the dynamic compression of a gas sphere, and some problems with
stationary shocks.
Keywords: shock, analytical solution, ideal gas, spherical symmetry, stationary
shock
1 Introduction
The properties of difference schemes that approximate the conservation laws are often
evaluated with a priori methods which involve studies into stability, approximation,
monotonicity, conservatism, distraction etc. It should however be noted that most of
these methods are developed for acoustic approximations and simple equations of state.
In continuum mechanics, the properties of a mathematical model may notably differ
from what linear theory predicts due to nonlinearities induced by real-world equations
of state, shocks, plasticity and other material properties.
Linear theory looses its rigor when applied to nonlinear equations. The importance
of the convergence theorem [1] is strongly exaggerated because it is still proved for
linear equations and not for nonlinear ones, and real calculations are done not for van-
ishing but finite ∆x and ∆t. A very vivid discussion of stability and convergence criteria
and their rigidity can be found in [2].
The calculation of shock waves strong discontinuities in all material properties re-
quires special attention. On the shock surface the conservation laws take the form of
nonlinear algebraic equations which relate the values of quantities across the shock.
Entropy jumps as all the other functions do. This is the fundamental difference between
a shock and a wave where the quantities vary continuously. Flows with shocks are often
357
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
simulated with homogeneous methods which treat the strong shock as a layer of a finite
width comparable with the cell size. This ability of difference schemes is called distrac-
tion [3]. Since the states across the shock are related by the Hugoniot, there must be a
mechanism which allows entropy to grow in the shock distraction region. Physical vis-
cosity and heat conduction in continuum mechanics equations cannot give a distraction
width of several cell sizes. The proposal by Neumann and Richtmyer to use a math-
ematical ’viscosity’ [4] seems to resolve the problem. Their method has gained wide
acceptance. Pseudo-viscosity is taken in different forms linear, quadratic, or linear-
quadratic [4-7]. But the method does not ensure convergence to the exact solution if
the form of pseudo-viscosity changes. So, the author of [8] gives an example where
different schemes with different viscosities converge to different solutions in the limit.
In [9], there is an example of spherical convergence where energy dissipation defined
by pseudo-viscosity is shown to be several times higher than energy dissipation due to
plasticity.
Advantages and disadvantages of a mathematical model can be seen from compar-
ison between its predictions and analytical solutions. The paper discusses some prob-
lems which have exact analytical solutions. It gives their statements (initial and bound-
ary conditions, equations of state, and physical parameters) and solutions in the form of
formulas or tables. In order to verify performance of a difference scheme, one needs to
solve the problem numerically and compare the result with the exact solution.
The problems are broken into two groups for stationary and non-stationary shocks.
In problems with stationary shocks, derivatives in the exact solution are zero every-
where beyond the distraction zone and hence approximation errors are also zero. In
the distraction zone, the derivatives and approximation errors reach high values. Here
the strong shock is smeared over several cells where entropy differs. These problems
help verify real shock distraction, monotonicity, entropy variation, and the dependence
of calculated results on the relation between steps in space and time, and on cell size
(the number of points in the mesh). All these properties reveal themselves differently
for strong and weak shocks. Shock strength is characterized be the difference between
pressures behind and before the shock.
In problems with non-stationary shocks, the derivatives and derivative-dependent
approximation errors are high beyond the distraction zone. This group includes shock
convergence and spherical shell convergence problems. In the last problem, the bound-
ary conditions and released energy are adjusted so as to keep material density constant
despite large pressure and velocity gradients.
Some analytical solutions are used for comparison with results obtained with the
difference schemes which are based on the energy dissipation method described in [10-
13].
2 Stationary shock
Consider a material with parameters P0 , V0 , E0 , U0 which do not change with time. At
a time t0 its left boundary instantaneously starts moving at a constant positive velocity,
producing a shock wave which propagates into the material. The equations
358
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
WV+ + U+ = WV− + U− , (1)
WU+ − P+ = WU− − P− , (2)
Wε+ − P+ U+ = Wε− − P− U− (3)
relate the material states P− = P0 , V− = V0 , E− = E0 , U− = U0 with the state after
the discontinuity P = P+ , V = V+ , E = E+ , U = U+ before and behind the shock, and
the shock velocity W. The number of quantities is larger than the number of equations
(1) - (3) + equation of state. To solve this system of equations requires that one of the
quantities be taken as parameter. Let it be U. Take the equation of state (EOS) for ideal
gas in the form
PV = (γ − 1) E (4)
and transform equations (1) - (3) to the dependence of P on U and other zero-subscripted
quantities
√︃
)︃2
γ + 1 ∆U 2 γ + 1 ∆U 2 γP0
(︃
P = P0 + + + ∆U 2 , (5)
4 V0 4 V0 V0
where ∆U = U − U0 . From equations (1) - (3) and (5) we find P, W, V and E:
W = (P − P0 ) / (U − U0 ), V = V0 −(U − U0 ) /W, E = E0 +0, 5 (P + P0 ) (V − V0 ) . (6)
For condense matter, a simple equation of state has the form
2
P = (n − 1) ρE + C0k (ρ − ρ0k ) , (7)
where ρ = 1/V - is material density, ρ0k - is its density at a point with coordinates
T = 0, P = 0 and C 0k - is sound velocity at this point. For EOS (7) equations (1) - (3)
transform to the Hugoniot equation
√︃
)︃2
n+1 γ+1
(︃
ρ0 ∆U 2 + 2
(︁ )︁
P = P0 + ρ0 ∆U 2 + ρ0 ∆U 2 nP0 + ρ0k C0k . (8)
4 4
Equations (6) for W, V, and E remain the same.
For convenience, we treat all quantities in equation (5) - (8) as dimensionless. That
is why both in gas and in condense matter, the initial dimensionless density is unity.
Conversion to density in g/cm3 is done through multiplying by the constant used for
conversion to dimensionless density. All the other quantities are treated similarly.
Problem 1. Strong shock in monatomic gas. At t = 0, a region 0 ≤ x0 ≤ 1 is
occupied with gas described by EOS (4) with parameters γ = 5/3, P0 = 0, ρ0 = 1,
E 0 = 0, U 0 = 0. Here x0 – is the Eulerian coordinate at t = 0. In Lagrangian difference
schemes x0 is a Lagrangian coordinate. At t > 0, U = 1 is specified on the left boundary
(x0 = 0) and U = 0 is on the right one (x0 = 1).
The quantities behind the shock front and front velocity W are determined from
equations (1) - (3): ρ = 4, E = 0.5, P = 4/3, W = 4/3. At t = 0.375, the shock is at a
point x0 = 0.5 and the analytical solution is determined by
359
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
P = 1.33333, ρ = 4.0, E = 0.5, U = 1.0, for x0 ≤ 0.5, and
P = 0, ρ = 1.0, E = 0, U = 0, for x0 > 0.5
Figures 1 and 2 depict P(x0 ) and U(x0 ) at t = 0.575. The solid lines show analytical
solutions and the marked ones show calculations with the difference scheme from [12].
The calculations were done with a uniform mesh of N = 100 points in x0 and Courant
number 0.5.
P U
1.2
0.8
0.8
0.4
0.4
0 0.0
0.4 0.45 0.5 0.55 x0 0.4 0.45 0.5 0.55
Fig. 1. Problem 1. P(x0 ) at t = 0.375 Fig. 2. Problem 1. U(x0 ) at time t = 0.375
Problem 2. Strong shock in monatomic gas described by EOS (4) and γ = 1.25
(ethylene). All the other parameters are the same as in Problem 1: ρ0 = 1, P0 = 0,
E0 = 0, U0 = 0. The boundary conditions are also the same: U(x0 = 0, t) = 1, U(x0 =
1, t) = 0.
The quantities behind the shock front and front velocity are determined from equations
(1) - (3): ρ = 9, E = 0.5, P = 1.125, W = 1.125. At t = 0.44444 the shock is at a
point x0 = 0.5, and the analytical solution is determined by
P = 1.125, ρ = 9.0, E = 0.5, U = 1 for x0 ≤ 0.5, and
P = 0, ρ = 1.0, E = 0, U = 0, for x0 > 0.5.
Figures 3 and 4 depict P(x0 ) and U(x0 ) at t = 0.44444. The solid lines show analytical
solutions and the marked ones show calculations with the difference scheme from [12].
The calculations were done with a uniform mesh of N = 100 points in x0 and Courant
number 0.5.
Problem 3. The weak shock wave in monatomic gas. At t = 0, a region 0 ≤ x0 ≤ 1
is occupied with monatomic ideal gas described by EOS (4) with parameters γ = 5/3,
ρ0 = 1, P0 = 1, ρ0 = 1, E0 = 1.5, U0 = 0. At t > 0, U = 0.5 is specified on the
left boundary (x0 = 0) and U = 0 is on the right one (x0 = 1). Behind the shock,
ρ = 1.428573, E = 1.925, P = 1.833333, and W = 1.666666. The analytical solution at
t = 0.3 is determined by
360
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
P U
0.9
0.8
0.6
0.4
0.3
0.0 0
0.4 0.45 0.5 0.55 0.4 0.45 0.5 0.55
Fig. 3. Problem 2. P(x0 ) at t = 0.44444 Fig. 4. Problem 2. U(x0 ) at t = 0.44444
P = 1.83333, ρ = 1.42857, E = 1.925, U = 0.5 for x0 ≤ 0.5, and
P = 1.0, ρ = 1.0, E = 1.50, U = 0 for x0 > 0.5.
Figures 5 and 6 depict P(x0 ) and U(x0 ) at t = 0.3. The solid lines show analytical
solutions and the marked ones show calculations with the difference scheme from [13].
The calculations were done with a uniform mesh of N = 100 points in x0 and Courant
number 0.5.
P U
1.8
0.4
1.6
1.4
0.2
1.2
1.0 0
0.3 0.4 0.5 0.3 0.4 0.5
Fig. 5. Problem 3. P(x0 ) at t = 0.3 Fig. 6. Problem 3. U(x0 ) at t = 0.3
361
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Problem 4. Strong shock in condense matter. At t = 0, a region 0 ≤ x0 ≤ 1 is
occupied with condense matter described by EOS (7) which at ρ0k = 1 and C0k = 1
takes the form
P = (n − 1) ρE + ρ − 1, (9)
At t = 0, the parameters are ρ0 = 1, E0 = 0, P0 = 0, U0 = 0, n = 3. For t > 0, U = 2
is specified on the left boundary (x0 = 0) and U = 0 is on the right one (x0 = 1). The
equation for P is obtained from (8):
√︃
)︃2
n+1 n+1
(︃
2
P= ρ0 U + ρ0 U + ρ0 U 2 .
2
4 4
P U
2.0
6
1.5
1.0
3
0.5
0 0.0
0.4 0.45 0.5 0.55 0.4 0.45 0.5 0.55
Fig. 7. Problem 4. P(x0 ) at t = 0.118034 Fig. 8. Problem 4. U(x0 ) at t = 0.118034
The analytical solution at t = 0.118034 is determined by
P = 8.47214, ρ = 1.89443, E = 2.0, U = 2.0 for x ≤ 0.5, and
P = 0, ρ = 1.0, E = 0, U=0 for x > 0.5.
Figures 7 and 8 depict P(x0 ) and U(x0 ) at t = 0.118034. The solid lines show analytical
solutions and the marked ones show calculations with the difference scheme from [13].
The calculations were done with a uniform mesh of N = 100 points in x0 and Courant
number 0.5.
Problem 5. Weak shock in condense matter. At t = 0, a region 0 ≤ x0 ≤ 1 is
occupied with a material described by EOS (9) with parameters: ρ0 = 1, E0 = 2,
P0 = 4, U0 = 0, n = 3. For t > 0, U = 1 on the left boundary and U = 0 on the right
one.
The analytical solution at t = 0.105448 is determined by
P = 8.74166, ρ = 1.26726, E = 3.34359, U = 1.0 for x0 ≤ 0.5, and
362
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
P U
8
0.8
7
6
0.4
5
4 0
0.4 0.45 0.5 0.55 0.4 0.45 0.5 0.55
Fig. 9. Problem 5. P(x0 ) at t = 0.105448 Fig. 10. Problem 5. U(x0 ) at t = 0.105448
P = 4.0, ρ = 1.0, E = 2.0, U=0 for x0 > 0.5.
Figures 9 and 10 depict P(x0 ) and U(x0 ) at t = 0.105448. The solid lines show analytical
solutions and the marked ones show calculations with the difference scheme from [13].
The calculations were done with a uniform mesh of N = 100 points in x0 and Courant
number 0.5.
3 The motion of a spherical layer of compressible ideal fluid
The problem of bubble collapse in fluid, or spherical shell convergence, arises in con-
nection with cavitation corrosion of propellers. Solutions to the problem can be found
in [14-16]. The full continuum mechanics model for 1D spherically symmetric flow
of ideal compressible continua includes mass conservation, motion and internal energy
equations:
∂V ∂r2 U
− 4π = 0, (10)
∂t ∂M
∂U ∂P
+ 4πr2 = 0, (11)
∂t ∂M
∂E ∂V
+P = 0. (12)
∂t ∂t
Equations (10)–(12) are written in Lagrangian coordinates. Here the partial deriva-
tives with respect to time are substantial derivatives.
There no incompressible matter in nature. Mechanics simply considers a wide
class of flows where density remains constant in time. Density constancy is often un-
derstood as incompressibility, i.e., β s = 0 and C 2 = ∞. It is not true. The property of
flow is not the property of matter.
363
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
As a rule [16], the models of ’incompressible’ fluid do not include the energy con-
servation law and the equation of state. This makes them internally contradictive. As
follows from (12), in ideal compressible fluid, E = const at V = const. But in this case
from the equation of state P = P(V, E) we obtain that P = const, too. So, on the one
hand, P varies, and on the other hand, P is constant. This contradiction can be removed
if assume that fluid is not adiabatic, i.e., there is a source of energy in it. Then equation
(12) is written as
∂E ∂V ∂q
= −P + . (13)
∂t ∂t ∂t
Equations (10), (11) and (13) allow solutions where density is constant. To keep the
density of fluid constant requires energy. As follows from the theory of equations of
state [17], in fluid, thermal pressure and energy, PT and ET , are related by the equation
PT V = Γ(V)ET . Hereafter for V = const the equation is taken in the form
PV0 = ΓE, (14)
where Γ = const, P = PT , E = E T . Since P(t, M) is a solution to equations (11)
and (12), then the dependence E(t, M) which follows from (13) is quite specific in each
flow. It is defined by the necessity of meeting the condition V = const.
Here we limit ourselves to flows where specific volume is independent of either M,
or t, i.e., V = const. Also, we assume that fluid is compressible, i.e., its compressibility
βS is nonzero.
For V = V0 , the(︁ equation
⧸︁ )︁ that relates the Eulearian coordinate r and the Lagrangian
coordinate dM = 4πr2 V0 dr can be integrated from M = 0 at r = rB to an arbitrary
finite M
)︁1/3
r = r3B + 3V0 M/4π .
(︁
(15)
Here rB is the time dependent coordinate of the bubble boundary. At V = V0 , equation
(10) has the solution
r2 U = f (t) . (16)
Since f is independent of M, equations (15) and (16) are valid for arbitrary M. On the
bubble boundary where M = 0, equation (16) takes the form
r2
r2B U B = f (t) , U = U B B2 . (17)
r
where UB is bubble boundary velocity.
Find U B from (17) and substitute in the bubble boundary motion equation
drB
(︃ )︃
= UB. (18)
dt M
Integrate (18) together with (15), to obtain the dependence of rB on t
364
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
(︃ ∫︁ t )︃1/3
rB = r3B0 + 3 f (t) dt . (19)
t0
It is seen from (17) and (19) that the motion of the bubble boundary is completely
defined by f (t). If f (t) < 0, then U B < 0 too, i.e., the bubble collapses. The boundary
convergence time tf is found from (19) at r B = 0
∫︁ t f
r3B0 + 3 f (t) dt = 0. (20)
t0
Following [16], consider shell convergence with zero pressure on the inner bound-
ary and f (t) = U B0 r2B0 = const. At time t0, the radius and velocity of the inner bound-
ary, rB0 and U B0 < 0, are defined. In accord with (17) at time t0 , velocity depends on
radius
U = U B0 (rB0 /r)2 . (21)
Let all quantities on the outer boundary be subscripted ”a”. Assume that the shell mass
is equal to Ma . The coordinate of the outer boundary, ra, relate to that of the inner
)︁1/3
boundary rB as ra = r3B + b , where b = 3V 3 3
(︁
4π Ma = ra0 − r B0 . Pressure and velocity on
0
the outer boundary are
⎛(︃ )︃3 ⎞−2/3
⎜⎜⎜ ra0 t − t0 ⎟⎟⎟⎟
Ua = U B0 ⎜⎝ ⎜ − , (22)
rB0 t f − t0 ⎠
⎟
)︃⎞−4/3 ⎞⎟
⎜⎜⎜ t f − t −4/3 ⎜⎜⎜ ra0 3
2 ⎜(︃
⎛
U B0
⎛(︃ )︃ (︃
t t
)︃
− 0
Pa = ⎟⎟ ⎟⎟⎟ .
⎟⎟⎟ ⎟⎟
− ⎜⎜⎝ − (23)
2V0 t f − t0 rB0 t f − t0 ⎠ ⎠
⎜⎝⎜
The values Ua0 and Pa0 are found from (22) and (23) at t = t0 . In Lagrangian coordinates
the pressure, velocity and released energy are defined by
2 ⎜(︃ ⎞−4/3 ⎞
U B0
)︃−4/3 ⎛
⎜⎜⎜ t f − t ⎜⎜⎜ t f − t
⎛
3V0 M ⎟⎟⎟ ⎟⎟⎟⎟
P= −⎝ + ⎟⎠ ⎟⎠ ,
2V0 ⎝ t f − t0 t f − t0 4πr3B0
⎜ ⎜
⎞−2/3
⎜⎜⎜ t f − t 3V0 M ⎟⎟⎟
⎛
U = U B0 ⎝⎜ + , (24)
t f − t0 4πr3B0
⎟⎠
3 ⎜⎛ ⎞−7/3 (︃
∂q 2U B0 ⎜⎜⎜⎜⎜⎜ t f − t t f − t −7/3 ⎟⎟⎟
⎛ )︃ ⎞
3V0 M ⎟⎟⎟
= + − ⎟⎟⎠ .
∂t ΓrB0 ⎝ t f − t0 4πr3B0 t f − t0
⎜⎜⎝ ⎟⎠
As the reference problem we consider the motion of 10% shell.
Problem 6. The motion of a 10%-shell. At t0 = 0, rB0 = 1, ra0 = 1.1, V0 = 1,
and Ma = 1.38649. The velocity of the inner boundary is U B0 = 1. The EOS of shell
material with parameters ρ0k = 1, C0k = 1 and Γ = 2 has the form P = 2ρE + ρ − 1. At
the initial time, pressure, specific internal energy and velocity in the shell are defined
by
365
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
)︃−4/3 ⎛ )︃−4/3 ⎞
3M 1 ⎜⎜⎜⎜ 3
(︃ (︃
P (M) = 1 − 1 + , E (M) = ⎜⎝1 − 1 + ⎟⎟⎠ ,
⎟⎟⎟
4π 4 4π
⎞2/3
⎜⎜ 4πr3B0 ρ0 ⎟⎟⎟
⎛
U = U B0 ⎜⎜⎝ ⎟⎠ .
3M
)︁−2/3
Boundary conditions: at t ≥ 0, PB = 0, MB = 0, Ua = − 1, 13 − 3t
(︁
at t ≥ 0,
M = M a . From (20) we found t f = 1/3. At t ≥ 0, energy release as a function of t
and M is defined by
)︃−7/3
dq 3M
(︃
= 1 − 3t + − (1 − 3t)−7/3 .
dt 4π
For t ≥ 0 and M a ≥ M ≥ 0 the solution has the form
⎛ )︃−4/3 ⎞
1 ⎜⎜⎜⎜ 3M
(︃
E (t, M) = + ⎟⎟⎠ , P = 2E, ρ = 1,
⎟⎟⎟
(1 − 3t) −4/3
− 1 − 3t
4 4π
⎜⎝
)︃−2/3
3M
(︃
U (t, M) = − 1 − 3t + .
4π
Figures 11 and 12 show pressure and velocity as functions of m = M/M a at t1 = 0, 30,
t2 = 0, 32.
P U
-2
30
-4
20
-6
10
-8
0 -10
0.0 0.2 0.4 0.6 0.8 m 0 0.2 0.4 0.6 0.8 m
Fig. 11. Problem 6. Fig. 12. Problem 6.
P(m) at t1 = 0.3 (line 1) and t2 = 0.32 (line2) U(m) at t1 = 0.3 (line 1) and t2 = 0.32 (line2)
366
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
4 Shock in a gas sphere
In different years, there were published a number of papers [18-22] with self-similar
solutions to shock convergence in infinite ideal gas. Shock convergence in a gas sphere
of finite radius is considered in [23,24]. At t = t0 , pressure in gas is P0 = 0, density
ρ0 = const, velocity U0 = 0, and specific internal energy E0 = 0. The boundary of the
sphere is at a point (r0 , t0 ). Velocity on the boundary is Ug0 < 0. In other words, velocity
jumps on the boundary, producing a shock wave which moves into the sphere. At the
time when the shock converges, t f , its coordinate rw is zero. The equation of motion
which satisfies all these conditions is
tf − t n
(︃ )︃
rw = r0 (25)
t f − t0
for n > 0 Its differentiation gives shock velocity
t f − t n−1
(︃ )︃
D = D0 , (26)
t f − t0
where
⧸︁ (︁ )︁
D0 = −r0 n t f − t0 . (27)
Flow parameters are defined by
∂ρ ∂ρ ∂U 2ρU
+U +ρ + = 0,
∂t ∂r ∂r r
∂U ∂U 1 ∂P
+U + = 0, (28)
∂t ∂r ρ ∂r
∂P ∂P ∂U 2U
(︃ )︃
+U + γP + = 0.
∂t ∂r ∂r r
For solving the problem we change from the variables t and r to variables t and
ξ(t, r). The function ξ(t, r) is taken such as to remain constant on the shock. Its simplest
form reads as
)︃n
r t f − t0
(︃
ξ= . (29)
r0 t f − t
Now express P, ρ and U as functions of time multiplied by functions of ξ:
P = α p (t) Π (ξ) , ρ = αρ (t) δ (ξ) , U = αu (t) M (ξ) . (30)
Choose Π (ξ) , δ (ξ) , M (ξ) such that to allow them at ξ = 1 take the values
γ+1 2 2
δw = , Mw = , Πw = . (31)
γ−1 γ+1 γ+1
With these δw , Π w , M w the function αρ , αu and α p take the forms
367
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
t f − t0 1−n t f − t0 2(1−n)
(︃ )︃ (︃ )︃
αρ = ρ0 , αu = D0 , α p = ρ0 D20 . (32)
tf − t tf − t
After appropriate manipulation for conversion to the functions Π, δ, M and variables t,
ξ, we obtain equations for functions which only depend on ξ:
2Mδ n−1
(M − ξ) δ′ + δM′ + = 0, δM + M′ δ (M − ξ) + Π ′ = 0, (33)
ξ n
2 (n − 1) 2γMΠ
Π + Π ′ (M − ξ) + γΠM′ + = 0. (34)
n ξ
For M ′ , δ′ , Π ′ , these equations give a system of linear homogeneous equations. If its
determinant
Z = (M − ξ) γΠ − δ (M − ξ)2
(︁ )︁
(35)
is nonzero, the system has a unique solution. At the point ξ* where Z = 0, the matrix
of coefficients and the augmented matrix of coefficients should be considered. At this
point their ranks are identical and equal to 2, and all third-order minors are zero, hence
the system of equations (33) and (34) has a unique solution. It is easy to show that with
the zero third-order minors we come to
(n − 1) ξ* (2 (M* − ξ* ) − γM* ) + 2γnM* (M* − ξ* ) = 0. (36)
The value of n is found from the condition that equation (36) holds simultaneously with
Z* = (M* − ξ* ) γΠ* − δ* (M* − ξ* )2 = 0.
(︁ )︁
(37)
From equations (36) and (37) we find the appropriate values of n for each γ. This
solution was used to evaluate the accuracy of some shock calculation methods.
Problem 7. A cold gas sphere of radius rg0 = 1 with parameters P0 = 0, ρ0 = 1,
U0 = 0, Ug0 = 1, γ = 5/3. The boundary condition is defined through reverse transition
from t, ξ and Π, δ, M to t, M and P, ρ, U. Pressure and boundary velocity as functions of
time are presented in Table 1. Pressure, density and velocity profiles at t = 0.4, 0.45, 0.5
(marked 1, 2, 3) are shown in Figs. 13-15. The solid lines show the analytical solution
derived in this work, the lines with circles show calculations by the VOLNA code [25]
with no shock smearing, and the dashed lines show VOLNA calculations with shock
smearing. The calculations were done on a uniform mesh of 100 points in r at t = t0 .
In Fig.14, entropy traces are seen in the dashed density profiles, which are a result of
shock smearing on the boundary. Figure 16 depicts M(ξ), (ξ) and δ(ξ) for 1 ≤ ξ ≤ 5.
Table 1.
Problem 7. The boundary condition in the case of γ = 5/3
No. t U P No. t U P
1 0 -1.00000 1.333333 17 0,30 -1.045146 2.445445
368
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
Table 1. Continuation.
No. t U P No. t U P
2 0.04 -1.008791 1.415094 18 0,31 -1.044921 2.517083
3 0.07 -1.015161 1.484529 19 0,32 -1.044473 2.592849
4 0.10 -1.021256 1.562232 20 0,33 -1.043787 2.673060
5 0.13 -1.026982 1.649664 21 0,34 -1.042842 2.758063
6 0.16 -1.032224 1.748622 22 0,35 -1.041619 2.848236
7 0.18 -1.035380 1.822074 23 0,36 -1.040097 2.943996
8 0.20 -1.038211 1.902413 24 0,37 -1.038251 3.045798
9 0.22 -1.040657 1.990554 25 0,38 -1.036058 3.154141
10 0.23 -1.041717 2.037877 26 0,39 -1.033488 3.269573
11 0.24 -1.042655 2.087565 27 0,40 -1.030514 3.392698
12 0.25 -1.043463 2.139778 28 0,42 -1.023220 3.664748
13 0.26 -1.044129 2.194694 29 0,45 -1.008354 4.149469
14 0.27 -1.044643 2.252502 30 0,50 -0.969947 5.233492
15 0.28 -1.044992 2.313410 31 0,55 -0.899013 6.144230
16 0.29 -1.045165 2.377642
P ρ
8
10
6
5 4
2
3 2 1
3 2 1
0
0 0.2 0.4 R 0 0.2 0.4 R
Fig. 13. Problem 7. pressure profiles at times Fig. 14. Problem 7. density profiles at times
0.4 (1), 0.45 (2), 0.5 (3) 0.4 (1), 0.45 (2), 0.5 (3)
369
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
U M δ
3 2 1
-1
5
-2
M
δ
-3
0 0.2 0.4 R 1 2 3 4 ξ
Fig. 15. Problem 7. velocity profiles at times Fig. 16. Problem 7. profiles of dimensionless
0.4 (1), 0.45 (2), 0.5 (3) M(ξ), Π(ξ), δ(ξ)
References
1. P.D. Lax and R.D. Richtmyer, Survey of stability of linear finite difference equations. Com-
mun. Pure and Appl. Math. – 1956 – V. 9 – P. 267–293.
2. P.J. Roach, Computational Fluid Dynamics. Albuquerque: Hermosa Publishers, 1976.
3. V.F. Kuropatenko, and I.R. Makeyeva, Investigation of discontinuity distraction in shock cal-
culation methods. J. Mathematical Modeling – 2006 – . 18 – No. 3 – P. 120–128.
4. J. von Neumann, and R.D. Richtmyer, A method for the numerical calculation of hydrody-
namic shocks. Appl. Phys. – 1950 – V. 21 – No. 3 – P. 232–237.
5. M.L. Wilkins, Simulation of elastic-plastic flows. In Computational Methods in Hydrodynam-
ics, Moscow, MIR Publishers, 1967, P. 212–263.
6. A.A. Samarsky, and V.Y. Arsenin, Numerical solution of hydrodynamic equations with differ-
ent viscosities. J. Comp. Math. and Mathem. Phys. – 1961 – V. 1 – No. 2 – P. 357–380.
7. V.F. Kuropatenko, A pseudo-viscosity form. Transactions of the Siberian Branch of USSR
Academy of Sciences – 1967 – No. 3 – Is. 3 – P. 81–82.
8. V.F. Dyachenko, Numerical analysis of discontinuous solutions to quasi-linear systems. J.
Comp. Math. and Mathem. Phys. – 1961 – V. 1 – No. 6 – P. 1127–1129.
9. V.F. Kuropatenko, and Y.N. Andreyev, Simulation of dynamic processes in spherical and
cylindrical shells. J. Computational Continuum Mechanics – 2010 – V. 3 – No. 4 – P. 53–
67.
10. V.F. Kuropatenko, A shock capture method. SUSU Bulletin, Mathematical Modeling and
Programming Series – 2014 – V. 7 – No. 1 – P. 62–75.
11. V.F. Kuropatenko, A shock calculation method. Proceedings of USSR Academy of Sciences
– 1960 – V. 3 – No. 4 – P. 771–772.
12. V.F. Kuropatenko, I.A. Dorovskikh, and I.R. Makeyeva, The influence of difference scheme
properties on mathematical modeling of dynamic processes. J. Computing Technologies –
2006 – V. 11 – Part 2 – P. 9–11.
13. V.F. Kuropatenko, and M.N. Yakimova, A Method of Shock Calculation. J. Computational
and Engineering Mathematics – 2015 – V. 2 – No. 2 – P. 60–70.
370
Mathematical and Information Technologies, MIT-2016 — Mathematical modeling
14. E.I. Zababakhin, Energy cumulation and its bounds. J. Achievements in Physics – 1965 –
V. 85 – Is. 4 – P. 721–726.
15. K.V. Brushlinsky, and Y.M. Kazhdan, On self-similar solutions to some hydrodynamic prob-
lems. J. Achievements in Mathematics – 1963 – V. 18 – No. 2 – P. 3–23.
16. V.F. Kuropatenko, Collapse of spherical cavities and energy cumulation in compressible ideal
fluid. J. Combustion and Detonation Physics – 2015 – V. 51 – No. 1 – P. 57–65.
17. V.F. Kuropatenko, Equations of state for low-temperature dense plasma components. En-
cyclopedia of Low-temperature Plasma – B Series – V. VIII – Part 2 – Moscow: JANUS-K
Publishers – 2008 – P. 436–450.
18. G. Hunter, On the collapse of an empty cavity in water. J. Fluid Mechan. – 1960 – V. 8 –
No. 2 – P. 241–263.
19. L.I. Sedov, Unsteady flows of compressible fluid. Proceedings of USSR Academy of Sci-
ences – 1945 – V. 47 – No. 2 – P. 94–96.
20. K.P. Stanyukovich, Self-similar solutions to hydrodynamic equations with central symmetry.
Proceedings of USSR Academy of Sciences – 1945 – V. 48 – No. 5 – P. 331–333.
21. A.N. Krayko, Fast cylindrically and spherically symmetric strong compression of ideal gas.
J. Applied Mathematics and Mechanics – 2007 – V. 71 – No. 5 – P. 774–760.
22. S.P. Bautin, Mathematical modeling of strong gas compression. Novosibirsk, Nauka Pub-
lishers, 2007.
23. V.F. Kuropatenko, E.S. Shestakovskaya, and M.N. Yakimova, Dynamic compression of a
cold gas sphere. Proceedings of Academy of Sciences – 2015 – V. 461 – No. 5 – P. 530–532.
24. V.F. Kuropatenko, E.S. Shestakovskaya, and M.N. Yakimova, Shock in a gas sphere. SUSU
Bulletin Mathematical Modeling and Programming Series – 2015 – V. 8 – No. 4 – P. 14–29.
25. V.F. Kuropatenko, V.I. Kuznetsova, G.V. Kovalenko, G.I. Mikhaylova, and G.N. Sapozh-
nikova, VOLNA code and an inhomogeneous difference technique for unsteady flows of com-
pressible continua. J. Problems in Nuclear Science and Technology – Mathematical Modeling
of Physical Processes Series – 1989 – Is. 2 – P. 9–25.
371