QJMAM

Download as pdf or txt
Download as pdf or txt
You are on page 1of 18

Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS


GAS OF VARIABLE DENSITY
by ANTIM CHAUHAN† and RAJAN ARORA‡
(Department of Applied Science and Engineering, Indian Institute of Technology Roorkee,
Roorkee 247667, India)
and
AMIT TOMAR§

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


(Department of Mathematics, Amity University, Noida 201303, India)
[Received 29 August 2019. Revise 21 November 2019. Accepted 21 December 2019]

Summary

The converging problem of cylindrically or spherically symmetric strong shock wave collapsing
at the axis/centre of symmetry, is studied in a non-ideal inhomogeneous gaseous medium. Here,
we assume that the undisturbed medium is spatially variable and the density of a gas is decreasing
towards the axis/centre according to a power law. In the present work, we have used the perturbation
technique to the implosion problem and obtained a global solution that also admits Guderley’s
asymptotic solution in a very good agreement which holds only in the vicinity of the axis/centre of
implosion. The similarity exponents together with their corresponding amplitudes are determined
by expanding the flow parameters in powers of time. We also refined the leading similarity
exponents near the axis/centre of convergence. We compared our calculated results with the
already existing results and found them in good agreements up to two decimal places. Shock
position and flow parameters are analysed graphically with respect to the variation of values
of different parameters. It is observed that an increase in the density variation index, adiabatic
exponent and Van der Waals excluded volume, causes the time of shock collapse to decrease due
to which the shock acceleration gets increased and shock reaches the axis/centre much faster.

1. Introduction
The problem of imploding shock waves has been studied and investigated by many researchers all
over the world for various real thermodynamically imperfect media and become a field of continuing
research interest as it generates high temperature and pressure plasmas at the centre of convergence.
High temperatures and pressures at the centre of convergence of cylindrical or spherical shock waves
give rise to the various interesting applications such as detonation, fusion initiation and destruction of
kidney and ureteral stones (1). In the last few decades, due to the theoretical and practical importance
of converging shock waves in various astrophysical situations, plasma physics, nuclear science and
engineering physics, their study is receiving much attention from researchers and scientists in the
fields ranging from gas dynamics to condensed matter (2–4). The shock waves are also employed
in laboratories to maintain high temperature to study various processes which take place in gases.

† <antimchauhan1@gmail.com>
‡ <rajanfpt@iitr.ernet.in>
§ Corresponding author <amitmath.14@gmail.com>

Q. Jl Mech. Appl. Math, Vol. 73. No. 2 © The Author(s), 2020. Published by Oxford University Press;
all rights reserved. For Permissions, please email: journals.permissions@oup.com.
Advance Access publication 15 February 2020. doi:10.1093/qjmam/hbaa003

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 101 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

102 A. CHAUHAN ET AL.

All of these processes involve a cylindrical or spherical shock wave travelling through a region of
inhomogeneities. The shock waves generated by the cylindrical and spherical piston are of great
interest from both mathematical and physical point of views.
The natural scaling invariance of a problem can be determined from the dimensional analysis that
presents a way to find a similarity solution of a governing system of hyperbolic partial differential
equations (PDEs), a problem involving symmetry. The two most common examples of it are point
implosion and explosion problems which were studied by Taylor (5) and Sedov (6). The solutions
obtained from dimensional analysis are known as self-similar solution of first kind. Usually, these

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


types of solutions do not satisfy the imposed boundary conditions globally. However, by using
asymptotical analysis in fixed domains, it is possible to determine the symmetric solutions. Zeldovich
and Raizer (7) have given simple scaling arguments to obtain similarity solutions with an illustration
of self-similar or invariant nature of the scaled solutions. Guderley (8) was the first who theoretically
analysed the cylindrical or spherical symmetric converging shock wave near the axis/centre of
convergence, which is a remarkable example of a similarity solution of the second kind (also see the
references (9, 10)). Lazarus and Richtmyer (11) discussed the known results for similarity solutions
for the flow problems of a strong converging shock wave and extended that work in different ways. By
using Lie group of transformations method, Sekhar and Sharma (12) presented similarity solutions
to the nonlinear three-dimensional unsteady Euler equations of gas dynamics, in some cases, they
also obtained exact solutions. Ponchaut et al. (13) have studied the implosion problem in order to
obtain the most accurate results and approaches. Towards the study of shock phenomena, the works
of Jeffrey (14), Sari et al. (15), Pandey et al. (16), Boyd et al. (17), Whitham (18), Ramsey et al.
(19), Rogers and Castell (20), Hirschler and Gretler (21) and Haque and Broadbridge (22) provide
the major contributions.
As the strength of the converging shock increases, the effect of non-ideal gas becomes significant
and must be involved for theoretical study (23). Evidently, non-ideal gases are not exactly described
by the ideal gas model. There are always deviations from the ideal gas model in the behaviour of a
non-ideal gases, in particular with regard to shock wave phenomena. Thus, shock waves in non-ideal
gases exhibit richer behaviour than that predicted by the ideal gas model. In particular, the non-ideal
gas is exemplified by the Van der Waals equations of state. In various practical applications such as
chemical processes, aerospace engineering, nuclear reactions, mechanical engineering systems, the
study of shock-related phenomena through a non-ideal gas is found to be important. Wu and Roberts
(24) analysed the flow field behind the converging shock waves through a Van der Waals gas. Zhao
et al. (25) presented the complete classification of shock waves in a Van der Waals fluid. A particular
exact solution for the quasilinear hyperbolic system of PDEs governing planar or radially symmetric
unsteady motion in perfectly conducting non-ideal magnetogasdynamics has been obtained by Bira
et al. (26) by using Lie group symmetry analysis. Towards the study of converging shock in a
non-ideal gas, we mentioned the work by (10, 19, 27–29).
The situation in which a strong converging waves can propagate through a non-uniform medium
of decreasing density and finally reaching to the destination, where the density is almost zero (30),
is of great interest in astrophysics as it is highly relevant to the problem of the origin of cosmic rays
(7, 30, 31). The propagation of shock waves through a medium of variable-density have been treated
by Sakurai (32) and Rogers (33). Results obtained by them are more applicable to the shock formed
in the deep interior of stars. By considering the initial density of gas is either constant or decreasing
towards the centre according to a power law, a problem of converging strong shock wave in an
perfect quiescent gas with zero pressure have been investigated by Hafner (31) by mean of power
series representation of the non-dimensional flow velocity. Madhumita and Sharma (34) examined

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 102 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 103

the strong converging shock waves problem in an ideal gas of varying density which varies according
to the power law while the problem of strong converging shock waves in a non-ideal gas of constant
density was investigated by Arora and Sharma (27) and Tomar et al. (35) by using an alternative
approach proposed by Vandyke and Guttmann (36). By expanding the flow parameter in power of
time, a problem of cylindrical converging strong shock wave in non-ideal gas in presence of magnetic
field has been solved by Chauhan et al. (37). An arbitrary number of terms for the series solution
of realistic model of a spherical blast wave, in different regions have been calculated by Haque et
al. (38). In the present work, we extended the problem treated by Arora and Sharma (27) by taking

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


the medium as inhomogeneous non-ideal gas of varying density instead of constant density, with
possible applications in the study of explosions in rotating stars. We have used the method suggested
by the Van Dyke and Guttmann (36) and obtained the global solution to the shock implosion problem
in a non-ideal gas of varying density. From the obtained global solution which is valid throughout
the flow field, we extracted the Guderley’s local solution near the axis/centre of symmetry. Also,
we made the comparison of the leading similarity exponents obtained by us with already existing
similarity exponents in references (27, 31) as shown in Table 4, and found that the results are in very
good agreement. The effects of the variation of the initial density index, adiabatic exponent of the
gas and Van der Waals excluded volume on the flow field behind the shock wave are investigated.
It is observed that the gas pressure remains bounded in the region behind the shock wave even it is
stationary in most of the region except in the vicinity of the front where it attains a maximum while
the density decreases ahead of shock front at the same rate at which the square of the front velocity
increases.
The article is summarised as follows: in section 2, we have introduced the system of basic equations
governing one dimension unsteady, adiabatic, planar, cylindrically or spherically symmetric motion
of a non-ideal gas together with initial conditions and Rankine–Hugoniot (R–H) conditions. Also, in
section 2, the flow parameters and shock path are analysed by expanding the solution in Taylor’s series
in powers of the time and the coefficients in Taylor’s series expansion are calculated by writing the
programs in software package ‘Mathematica7’. In section 3, the coefficients in the series expansion
have been analysed together with the radius of shock waves as shock reaches the axis/centre for
different values of parameters m, γ , δ and b. In section 4, the similarity exponents together with
the corresponding amplitudes are calculated and are listed in Table 4. In section 5, we compared
the results computed by us with the already existing results by Arora and Sharma (27) and Hafner
(31). We analysed the flow parameters and shock location graphically, in section 5. Finally, section 6
concludes the article.

2. Basic equations and solutions in power of time


We consider a cylindrical or spherical piston of an initial radius R0 that bounds a quiescent non-ideal
gas, and the density of gas initially varies according to the power law. At time t = 0, the piston begins
contracting with very high constant velocity V which is larger than the acoustic speed of the medium
and producing a cylindrical or spherical symmetric strong shock wave whose position R(t) has to be
found out (see Fig. 1). We take the initial conditions for the undisturbed medium as follows:

v = 0, ρ = ρ0 = ρc (r/R0 )δ , p = p0 , (2.1)

here v, ρ, p denote the outward radial velocity, density, pressure, respectively, while ρc , p0 , δ are
the suitable positive constants.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 103 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

104 A. CHAUHAN ET AL.

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


Fig. 1 The constant velocity converging cylindrical/spherical piston problem (4).

The system of equations governing one dimension unsteady, adiabatic, planar, cylindrically or
spherically symmetric motion of a non-ideal gas are given by (18, 27):

∂ρ ∂ρ ∂v mρv
+v +ρ + = 0, (2.2)
∂t ∂r ∂r r
∂v ∂v 1 ∂p
+v + = 0, (2.3)
∂t ∂r ρ ∂r
 
∂p ∂p ∂ρ ∂ρ
+v − a2 +v = 0, (2.4)
∂t ∂r ∂t ∂r

where t denotes the time and r denotes the spatial coordinate which is radial in cylindrically or
spherically symmetric flow, b is a non-ideal parameter of gas, γ is the adiabatic coefficient and m is
0 for planar flow, 1 for cylindrically and 2 for spherically symmetric flow. The speed of sound a is
defined as a = (γ p/ρ(1 − bρ))1/2 .
In order to analyse, how the solution can be affected by the deviation of an ideal gas, the equation
of state of the medium is assumed to be in the form of the Van der Waals type (Wu and Roberts
(24)):
p(1 − bρ) = ρRT , (2.5)

where R and T denote the gas constant and the absolute temperature, respectively.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 104 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 105

Behind the shock front r = R(t), the R–H jump conditions are given by:
2(1 − b)
v= Ṙ,
γ +1
γ +1
ρ= ρ0 , (2.6)
γ − 1 + 2b
2(1 − b)
p= ρ0 Ṙ2 ,
γ +1

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


where Ṙ(t) denotes the shock speed, and ρ0 is the initial density of the gas.
Let us assume that there is no flow through the piston whose condition is given by

v = −V at r = R0 − Vt, (2.7)

where, v is the gas velocity relative to the fixed coordinates in which the gas is at rest in the piston
with radius R0 . For simplicity, we calculate the distance x = R0 − r and the corresponding velocity
u = −v directed inward. By considering the canonical nature of basic flow, we take a new variable
of the following form
2 x 
y= −k , (2.8)
(γ − 1) Vt
where the constant k has to be determined. The variable y varies from−2(k − 1)/(γ − 1) at the piston
to unity at the basic position of shock wave. The value of this constant k = 1 for an ideal gas.
We introduce the dimensionless variables by taking lengths to R0 , density to ρc , speed to V ,
pressure to ρc V 2 , time to R0 /V and non-ideal parameter b to 1/ρc . The system (2.2)–(2.4) and
boundary conditions defined in (2.6) can be rendered dimensionless by using the dimensionless
variable. The obtained system of PDEs together with boundary conditions in terms of dimensionless
variables are transformed into new coordinate system (y, t) as given below:
      
1 ∂u 1 ∂ρ 1 ∂ρ
1 − k + (γ − 1)y t ρ + (γ − 1)t + u − k − (γ − 1)y
2 ∂y 2 ∂t 2 ∂y
1
(γ − 1)mtρu, = (2.9)
2
 
1 ∂u 1 ∂u ∂p
ρ u − k − (γ − 1)y + (γ − 1)tρ + = 0, (2.10)
2 ∂y 2 ∂t ∂y
  
1 ∂p ∂ρ
u − k − (γ − 1)y ρ(1 − bρ) − γp
2 ∂y ∂y
 
1 ∂p ∂ρ
+ (γ − 1)t ρ(1 − bρ) − γp = 0. (2.11)
2 ∂t ∂t
Let the shock front be at distance χ which is measured inward from the initial position of piston
(i.e. χ = R0 − R). Then, the boundary conditions at shock (2.6) and piston (2.7), become
2(1 − b) (γ + 1) 2(1 − b)
u= χ̇ , ρ= (1 − χ )δ , p= (1 − χ )δ χ̇ 2 , (2.12)
(γ + 1) (γ − 1 + 2b) (γ + 1)

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 105 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

106 A. CHAUHAN ET AL.

at y = (2/(γ − 1)) (χ /t − k) and


−2(k − 1)
u=1 at y= . (2.13)
(γ − 1)
By considering that the solution is analytic in time, we expand the unknown position χ of the
shock wave and all the flow variables in Taylor as follows



Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


χ (t) = χn t n , (2.14)
n=1

and similarly, we expand the flow variables as



 ∞
 ∞

u= Un t n−1 , ρ= Rn t n−1 , p= Pn t n−1 . (2.15)
n=1 n=1 n=1

On substituting (2.15) into (2.9)–(2.11) and using (2.12), (2.13), (2.14), we found the coefficients
for the first-order approximation U1 , R1 and P1 by equating the coefficients of like powers of t

(γ + 1) (γ + 1) (γ + 1) 2 + b(γ − 1)
U1 = 1, R1 = , P1 = , χ1 = , k= .
(γ − 1 + 2b) 2(1 − b) 2(1 − b) 2(1 − b)
The coefficients for the second-order approximation U2 , R2 and P2 satisfy the following first-order
linear ODEs:
 
 1  1
R1 U2 − k − 1 + (γ − 1)y R2 + (γ − 1)R2
2 2
1
= m(γ − 1)R1 U1 , (2.16)
2
 
1 1  
(γ − 1)R1 U2 − k − 1 + (γ − 1)y R1 U2 + P2 = 0, (2.17)
2 2
(γ − 1)
((1 − bR1 )R1 P2 − γ P1 R2 )
2
  
1  
− k − 1 + (γ − 1)y (1 − bR1 )R1 P2 − γ P1 R2 = 0. (2.18)
2

Here, first-order derivative with respect to y is denoted by the prime. The boundary condition
(2.13) on piston becomes

−2(k − 1)
Un (y) = 0 at y= , for n = 2, 3, 4, ... (2.19)
γ −1
and boundary condition (2.12) at shock becomes
 
4(1 − b) δ(γ + 1)2 γ +1 2
U2 (1) = χ2 , R2 (1) = − , P2 (1) = 4χ2 − δ. (2.20)
(γ + 1) 2(γ − 1 + 2b)(1 − b) 2(1 − b)

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 106 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 107

From the first-order linear PDEs (2.16)–(2.18), we found the following conditions
  
U2 (y) = R2 (y) = P2 (y) = 0. (2.21)

By using the boundary conditions (2.19) and (2.20) together with the condition (2.21), we obtained
the coefficients for the second-order approximation U2 , R2 and P2 as follows:
 
δ(γ − 1)(γ + 1) + 2mγ (γ − 1 + 2b)  b(γ + 1) 

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


U2 = + (γ − 1)y , (2.22)
4(2γ − 1)(γ − 1 + 2b) (1 − b)
(1 − b)(γ − 1)2 (γ + 1)   
R2 = 2m(γ − 1 + 2b) − δ(γ + 1) 1 − y ,
2(γ − 1 + 2b) (2γ − 1)
2

δ(γ + 1)2  
− b + (γ − 1)(1 − b)y , (2.23)
2(γ − 1 + 2b)2 (1 − b)
γ (γ + 1)  
P2 = 2m(γ − 1 + 2b) − δ(γ + 1) , (2.24)
4(1 − b)2 (2γ − 1)
(γ + 1)  
2
χ2 = 2mγ (γ − 1 + 2b) + δ(γ − 1) . (2.25)
16(1 − b)2 (2γ − 1)

In a similar manner as in (2.21), we can show that the coefficients Un , Rn and Pn for nth order
approximation are (n − 1)th degree polynomials in y of the following form for n > 2:


n 
n 
n
Un (y) = Unj yj−1 , Rn (y) = Rnj yj−1 , Pn (y) = Pnj yj−1 . (2.26)
j=1 j=1 j=1

The nth-order approximation can be calculated by substituting the values from (2.26) into the
system (2.9)–(2.11) together with the shock conditions (2.12). We obtained a system of 3n + 1 linear
algebraic equations in 3n + 1 coefficients Unj , Rnj , Pnj (j = 1, 2, …, n) and χn by equating the like
powers of y and t; thereafter on solving the obtained system of linear algebraic equation, we found
the coefficients of nth-order approximation. Up to the second-order approximation, we determine
the position of shock wave as follows:

(γ + 1) (γ + 1)  
χ (t) = t+ 2mγ (γ − 1 + 2b) + δ(γ − 1)(γ + 1) t2 + · · · .
2(1 − b) 16(1 − b)2 (2γ − 1)
(2.27)

From the above expressions, we observe that in the absence of density stratification (δ = 0), the
above results reduce to the results obtained by Arora and Sharma (27) while for the ideal gas (b = 0)
of varying density (δ  = 0), the above results reduce to the results obtained by Madhumita and Sharma
(34). Also, we clearly see that the results obtained by Van dyke and Guttmann (36) can be easily
recovered from the above obtained results. As we can see that after the second-order approximation,
the hand calculations become very lengthy, so we compute the succeeding terms of the series (2.14)
by writing the program in software package ‘Mathematica7’. We write the program in two parts for

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 107 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

108 A. CHAUHAN ET AL.

the computation of coefficient of χn , Un , Rn and Pn , where n ≥ 3. The first part is devoted to the
generation of the system of algebraic equations by using (2.15) and (2.26) while in second part we
evaluate the coefficients χnk , Unk , Rnk and Pnk by solving the algebraic equations obtained in the
first part. All the calculations have been done for planar, cylindrical and spherical geometry with
adiabatic coefficient γ = 7/5, 5/3, 6/5, non-ideal parameter b = 0, 0.01, 0.05, 0.1 and density exponent
(see (30, 31)) δ = 0, 0.5, 1, 2. In particular, we calculated first 31 coefficients in the series expansions
of the shock location χ (t) for the following sets of parameters (a) m = 0, γ = 7/5, b = 0, δ = 0.5,
(b) m = 1, γ = 5/3, b = 0.1, δ = 1 and (c) m = 2, γ = 7/5, b = 0.01, δ = 2, which are listed in

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


Table 1. We rounded off all the coefficients up to 12 significant digits for the reduction of round-off
and truncation errors.

3. Analysis of coefficients of series for shock position


It is clear from Table 1, that the nearest singularity of series function χ (t) lies on a positive real axis
because all the coefficients of series (2.14) which are shown in Table 1 are positive. Also, we can
see from Table 1 that the coefficients of series (2.14) have steady increment which indicates that
the radius of convergence of the series (2.14) for shock position is less than unity. The power series
solution break down at the time t = 1 because at this time the piston itself reaches to the axis/centre.
Now, we will show that the singularity which occurs at time t = 1 is the Guderley’s singularity. Now,
we assume that near the collapse of shock wave, it admits the Guderley’s similarity solution whose
position is given by


R(t) = 1 − χ (t) = 1 − χn t n
 
t α1
∼ A1 1 − as t → tc , (3.1)
tc

here, the leading similarity exponent is denoted by α1 and the corresponding amplitude is denoted
by A1 .
Therefore,
 
χn 1 1 + α1
∼ 1− as n → ∞. (3.2)
χn−1 tc n

From the above expression (3.2), we see that the value of 1/tc can be approximated by the ratio
χn /χn−1 , for the large value of n. By using this result, we have constructed a sequence of ratios
χn /χn−1 for the values of n = 22, 23,…, 31 and refined this estimate of 1/tc by forming a Appendix
Neville Table A1 for the radius of convergence tc and Appendix Neville Table A2 for estimating
1/tc by following the similar procedure by Gaunt and Guttmann (39) and Tomar et al. (35). We also
estimated the leading similarity exponent α1 as shown in Table 2 for the following set of parameters
m = 1, γ = 5/3, δ = 1 and b = 0.1.
From Appendix Table A1, we can observe that the radius of convergence tc of the series (2.16) is
less than unity which we have already observed earlier from Table 1. The time of shock collapse are
calculated for different set of values of parameters m,γ , b and δ and are listed in Appendix Table A3.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 108 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 109

Table 1 Coefficients in the Taylor series expansion (2.14)

Planar (m = 0) Cylindrical (m = 1) Spherical (m = 2)


n γ = 7/5, b = 0, δ = 0.5 γ = 5/3, b = 0.1, δ = 1 γ = 7/5, b = 0.01, δ = 2
1 1.2 1.481481481481 1.212121212121
2 0.0400000000000 0.411522633745 0.363228242016
3 0.0284444444444 0.489130496142 0.542169022751

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


4 0.0218105263158 0.596495982775 0.665160178909
5 0.0180139127305 0.836182073392 0.867032293613
6 0.0159410419371 1.28917776615 1.33342400080
7 0.0149514600526 2.09732370385 2.17175979057
8 0.0146820336046 3.54272897946 3.55189188851
9 0.0149358773168 6.16121714890 5.97315357536
10 0.0156152896211 10.9621424012 10.3697032462
11 0.0166825492595 19.8617603418 18.2736545950
12 0.0181382047637 36.5262627324 32.5126361422
13 0.0200099378234 68.0134460431 58.5718314231
14 0.0223476702518 127.985965446 106.722778655
15 0.0252223927629 243.028824912 195.977257159
16 0.0287273602632 465.115467765 362.227292356
17 0.0329810079862 896.289337758 673.962340538
18 0.0381313606656 1737.68552569 1261.41359752
19 0.0443619558961 3387.18186754 2372.27364161
20 0.0518994647392 6634.44878271 4480.68974608
21 0.0610233178318 13051.5298772 8497.90474002
22 0.0720777634132 25776.8099427 16176.4732706
23 0.0854869125699 51092.0488428 30893.4591956
24 0.101773478758 101601.284029 59176.6825452
25 0.121582103537 202651.064933 113671.661417
26 0.145708389026 405320.506092 218909.818198
27 0.175135042252 812748.844308 422565.240705
28 0.211076892180 1.63358041786 × 106 817465.212349
29 0.255036985616 3.29062238533 × 106 1.58463557484 × 106
30 0.308876526517 6.64207029745 × 106 3.07758233314 × 106
31 0.374902123903 1.34325300902 × 107 5.98761389238 × 106

4. Structure of local singular solution


By using the theory of Baker and Hunter (40), Van Dyke and Guttmann (36) and Chauhan et al. (37)
determined all the real exponents and their corresponding amplitudes for an ideal and non-ideal gas
of constant density in the absence and presence of magnetic field, respectively. In this section, we
calculated all the real exponent and their corresponding amplitudes, listed in Table 3, for an non-ideal
gas of varying density by following the procedure mentioned in references (36, 37, 41).

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 109 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

110 A. CHAUHAN ET AL.

Table 2 Neville Table for similarity exponent α1 for m = 1, γ = 5/3, δ = 1 and b = 0.1

n e0n Linear Quadratic Cubic Quartic


27 0.67589175 0.67402659 0.67326890 0.67263939 0.67198451
28 0.67582300 0.67396679 0.67318937 0.67252667 0.67185036
29 0.67575695 0.67390767 0.67310961 0.67241838 0.67174155
30 0.67569336 0.67384914 0.67302969 0.67231039 0.67160846

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


31 0.67563200 0.67379111 0.67294960 0.67220208 0.67147098

5. Results and discussion


A problem of converging shock wave has been discussed for a non-ideal gas of varying density which
obeys a power law. All the flow parameters are extended from the piston to the axis/centre of collapse
and are analysed by writing them in a Taylor series expansions. By using the perturbation series
approach (see Van Dyke and Guttmann (36)), we found global solution to the implosion problem
which also recovers the Guderley’s local self-similar solution near the axis/centre of symmetry. From
the coefficients in the perturbation series, we extracted the all possible similarity exponents together
with their corresponding amplitudes in Guderley’s local expansion and listed them in Table 3. For
different values of parameters m, γ , b and δ, the values of first dominant similarity exponent α1 are
shown in Table 4, and compared with the results obtained by Arora and Sharma (27) and Hafner
(31).
We have shown the Neville Table 2 in section 3 and Appendix Neville Tables A1 and A2 for
estimating the similarity exponent α1 , the time tc of shock wave to collapse at axis/centre of symmetry
and the reciprocal radius of convergence 1/tc , respectively, for the following set of parameters
m = 1, γ = 5/3, δ = 1 and b = 0.1. We observe that the values of leading similarity exponent (α1 )
shown in the Neville Table 2 and obtained by Pade approximation as shown in Table 3 are in good
agreement up to two decimal places. From the Guderley’s local self-similar solution, we can obtain
only the first dominant similarity exponent but the present method can provide us the other less
dominant similarity exponents and their corresponding amplitudes which are shown in the Table 3.
We noticed that the value of α1 is always less then unity, the reason behind this fact is that the shock
accelerated continuously and become unbounded as t → tc but less rapidly than (t − tc )−1 . From
Table 4, we also observe that increment in the values of any of the parameters m, γ , b and δ cause
the leading similarity exponent α1 to decrease due to which the shock acceleration gets increased as
it reaches the axis/centre of symmetry. The leading similarity exponent is diminished by non-ideal
perturbations causing the perturbed acceleration to increase.
The computations for large gamma (say gamma = 30), which is of great interest to certain parts of
the explosives community, where solids under high pressure must be modelled, have been done. We
observe that for this large value of γ , the time of shock collapse reduces to a very small value, and
shock collapses much faster. We found that as the value of b is increasing from 0 to 0.5, the value
of leading similarity exponent α1 is decreasing but if the value of b exceeds this range, the value of
leading similarity exponent again starts to increase. We found that the value of α1 becomes greater
than unity for b ≥ 0.7.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 110 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 111

Table 3 Similarity exponents and corresponding amplitudes for different values of m, γ , δ and b

m γ δ b Exponents Amplitudes
1 7/5 1 0.1 α1 = 0.68512567 A1 = 0.99486456
α2 = 2.31774664 A2 = 0.01854242
1 7/5 2 0.1 α1 = 0.60182799 A1 = 0.99803348
2 7/5 2 0.1 α1 = 0.52965070 A1 = 0.99854792
α2 = 0.53042091 = 0.00267264

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


A2
1 7/5 1 0.01 α1 = 0.71505122 A1 = 0.98704064
1 7/5 2 0.01 α1 = 0.62898650 A1 = 0.98856979
2 7/5 1 0.01 α1 = 0.62491570 A1 = 1.02759150
α2 = 2.47754911 A2 = 0.02148560
α3 = 6.71762715 A3 = 0.00002955
2 7/5 2 0.01 α1 = 0.55795501 A1 = 0.95497668
α2 = 2.55753804 A2 = 0.01501358
1 7/5 2 0.001 α1 = 0.63204375 A1 = 0.98840240
α2 = 2.83261748 A2 = 0.04094717
α3 = 3.13143796 A3 = 0.02098456
2 7/5 0 0.05 α1 = 0.67164484 A1 = 0.98401360
α2 = 2.08579887 A2 = 0.02297566
1 5/3 1 0 α1 = 0.689476826 A1 = 0.99007483
α2 = 2.15818921 A2 = 0.01855308
1 5/3 1 0.1 α1 = 0.67151726 A1 = 0.99775471
α2 = 2.30359117 A2 = 0.01291417
1 5/3 2 0.1 α1 = 0.58237740 A1 = 0.99443810
α2 = 2.32790498 A2 = 0.01572973
2 5/3 2 0.1 α1 = 0.50762511 A1 = 1.00328856
α2 = 2.37342469 A2 = 0.01051652
α3 = 4.20624277 A3 = 0.00048313
1 5/3 1 0.01 α1 = 0.68744588 A1 = 0.99076530
α3 = 2.16868018 A3 = 0.01802677
1 5/3 1 0.001 α1 = 0.68927047 A1 = 0.99014287
α2 = 2.15918986 A2 = 0.01850143
1 5/3 2 0.001 α1 = 0.60003668 A1 = 0.99182734
α2 = 2.24200957 A2 = 0.01897527

We have drawn graphs for the flow variables to explain the effects of variation of the adiabatic
coefficient (γ ), the wave front curvature (m), the density stratification (δ) and the non-ideal parameter
(b) of gas.
Figures 2(a,c)–7(a,c) show that as we move towards the piston, the velocity decreases
monotonically behind the shock wave while density shows increasing trend in the region behind
the shock. This is due to the area contraction or due to the geometrical convergence, which causes
the velocity to decrease and density to increase. In the region behind the shock wave, the gas which is
highly compressed by the shock, gets cooled down and cause the density to decrease ahead of shock

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 111 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

112 A. CHAUHAN ET AL.

Table 4 Comparison of computed values of the leading similarity exponents with Arora and
Sharma (27) and Hafner (31)

m γ δ b Computed α1 Arora and Sharma Hafner


(27) (31)
0 7/5 0.5 0 0.90339043 0.90619443
0 7/5 1 0 0.82946229 0.83185109
1 7/5 0 0.1 0.77817425 0.77411046

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


1 7/5 1 0.1 0.68512567
1 7/5 2 0.1 0.60182799
1 7/5 2 0 0.63240405 0.62834171
2 7/5 2 0 0.56201448 0.55880816
2 7/5 0 0 0.71993541 0.71717450
2 7/5 0 0.01 0.70966368
2 7/5 1 0.01 0.62491570
2 7/5 2 0.01 0.55795501
1 7/5 0 0.05 0.80566231 0.80292649
2 7/5 0 0.05 0.67164484 0.66154487
1 5/3 0 0.1 0.77553472 0.77939911
1 5/3 1 0.1 0.67151726
1 5/3 2 0.1 0.58237740
2 5/3 0 0 0.69075929 0.68837682
2 5/3 2 0 0.53088525 0.52638596
2 5/3 0 0.1 0.63260656 0.63276963
2 6/5 0 0.01 0.74170083 0.73449201
1 6/5 2 0 0.67045234 0.67176320
2 6/5 0 0 0.74646017 0.75714181
2 6/5 2 0 0.60660533 0.60734072

front at the same rate at which the square of the front velocity increases. Figures 2(b)–7(b) show that
the pressure remains bounded in the region behind the shock wave even it is stationary in most of
the region except in the vicinity of the front where it attains a maximum. From the figures, it is clear
that the amplification mechanism of flow convergence is stronger for the spherical geometry than
the cylindrical geometry. Further, Figures 2(a,c)–6(a,c), demonstrate that as the values of Van der
Waals excluded b and ambient density exponent δ increase, the velocity shows an increasing trend
while density decreases. From the Figures 2(d)–7(d), we see that the cylindrical and spherical shock
trajectories coincide in the neighbourhood of the piston indicating that the geometrical convergence
effects are small initially at large radii. In fact, in the intermediate region, shock accelerates gently
while near the axis/centre of symmetry it accelerates rapidly due to adiabatic compression of the
shocked state because of flow area convergence.
Appendix Table A3 shows that an increase in either of the parameters b, m, γ and δ, causes the
time of shock collapse tc to decrease i.e with the increase of these parameters, the shock reaches
the axis/centre much faster. The Figures 2(d)–7(d) show that shock path is an increasing function of

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 112 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 113

(a) (b)

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


(c) (d)

Fig. 2 Flow variables and shock path profiles for cylindrically symmetric flow fields for density stratification
δ = 1 and γ = 7/5.

time t. Also, this growth is further enhanced by an increase in the values of the parameters b, m, γ
and δ.

6. Conclusion
In this article, an imploding strong shock wave problem for a cylindrically or spherically symmetric
flow collapsing at the axis/centre of symmetry, is studied in a non-ideal inhomogeneous gas medium
i.e the undisturbed medium is spatially varying in which the density of gas is decreasing towards
the axis/centre according to a power law. We applied the perturbation series technique to the
implosion problem in which all the flow parameters are extended from the piston to the axis/centre
of collapse, and are analysed. Global solutions to the implosion problem recovering Guderley’s local
self-similar solution in the vicinity of the axis/centre of collapse are found. From the coefficients
in the perturbation series, we extracted all the possible similarity exponents together with their
corresponding amplitudes in Guderley’s local expansion. The values of first dominant similarity
exponents α1 calculated by us are compared with the results obtained by Arora and Sharma (27) and
Hafner (31). All the flow parameters and shock trajectory in the region extending from the piston
to the axis/centre of collapse, are analysed graphically with respect to the variation of different
parameters.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 113 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

114 A. CHAUHAN ET AL.

(a) (b)

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


(c) (d)

Fig. 3 Flow variables and shock path profiles for cylindrically symmetric flow fields for density stratification
δ = 2 and γ = 7/5.

(a) (b)

(c) (d)

Fig. 4 Flow variables and shock path profiles for spherically symmetric flow fields for density stratification
δ = 1 and γ = 5/3.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 114 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 115

(a) (b)

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


(c) (d)

Fig. 5 Flow variables and shock path profiles for spherically symmetric flow fields density stratification δ = 2
and γ = 5/3.

(a) (b)

(c) (d)

Fig. 6 Flow variables and shock path profiles for cylindrically symmetric flow fields for γ = 5/3 and b = 0.1.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 115 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

116 A. CHAUHAN ET AL.

(a) (b)

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


(c) (d)

Fig. 7 Flow variables and shock path profiles for cylindrically and spherically symmetric flow fields for
γ = 7/5, 5/3, b = 0.01 and density stratification δ = 2.

Acknowledgements
The first author Antim Chauhan acknowledges the research support from ‘University Grant
Commission (Govt of India)’ (Sr. No. 2121541039 with Ref No. 20/12/2015 (ii)EU-V).

References
1. K. Takayama and T. Saito, Shock wave/geophysical and medical applications, Annu. Rev. Fluid
Mech. 36 (2014) 347–379.
2. E. Bertschinger, Cosmological detonation waves, Astrophysical Journal, 295 (1985) 1–13.
3. M. Y. Jaffrin and R. F. Probstein, Structure of a plasma shock wave, Phys. Fluids 7 (1964).
4. S. D. Ramsey and R. S. Baty, Piston driven converging shock waves in a stiffened gas, Phys.
Fluids 31 (2019) 086106.
5. G. I. Taylor, The formation of a blast wave by a very intense explosion, theoretical discussion.
Proc. R. Soc. Lond. A, 201 (1950) 159–174.
6. L. I. Sedov, Similarity and Dimensional Methods in Mechanics (Academic Press, New York
1959).
7. Y. B. Zeldovich and Y. P. Raizer, Physics of Shock Waves and High Temperature Hydrodynamic
Phenomena, Vol. II (Academic Press, New York, NY 1967).
8. G. Guderley, Starke kugelige und zylindrische Verdichtungsstosse in der Nahe des Kugelmit-
telpunktes bzw der Zylinderachse, Luftfahrtforschung 19 (1942) 302–312.
9. J. D. Logan and J. D. J. Perez, Similarity solutions for reactive shock hydrodynamics, SIAM J.
Appl. Math. 39 (1980) 512–527.
10. Z. Boyd, S. D. Ramsey and R. S. Baty, On the existence of self-similar converging shocks in
non-ideal materials, Quart. J. Mech. Appl. Math. 70 (2017) 401–417.
11. R. B. Lazarus and R. D. Richtmyer, Similarity solutions for converging shocks, Los Alamos
Scientific Laboratory Report, LA-6823-MS, Los Alamos, NM (1977).

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 116 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

CONVERGING SHOCK WAVES IN A VAN DER WAALS GAS 117

12. T. R. Sekhar and V. D. Sharma, Similarity solutions for three dimensional Euler equations using
Lie group analysis, Appl. Math. Comput. 196 (2008) 147–157.
13. N. F. Ponchaut, H. G. Hornung, D. I. Pullin and C. A. Mouton, On imploding cylindrical and
spherical shock waves in a perfect gas, J. Fluid Mech. 560 (2006) 103–122.
14. A. Jeffrey, The formation of magnetoacoustic shocks, J. Math. Anal. Appl. 11 (1965) 139–150.
15. R. Sari, N. Bode, A. Yalinewich and A. MacFadyen, Slightly two- or three-dimensional self-
similar solutions, Phys. Fluids 24 (2012) 087102.
16. M. Pandey, R. Radha and V. D. Sharma, Similarity analysis and exact solutions of
magnetogasdynamic equations, Quart. J. Mech. Appl. Math. 61 (2008) 291–310.

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


17. Z. M. Boyd, S. D. Ramsey and R. S. Baty, On the existence of self-similar converging shocks
for arbitrary equation of state, Quart. J. Mech. Appl. Math. 70 (2017) 401–417.
18. G. B. Whitham, Linear and Nonlinear Waves (Wiley-Interscience, New York 1974).
19. S. D. Ramsey, E. M. Schmidt, Z. M. Boyd, J. F. Lilieholm and R. S. Baty, Converging shock
flows for a Mie-Grüneisen equation of state, Phys. Fluids 30 (2018) 046101.
20. C. Rogers and S. P. Castell, A parametrisation and non uniform shock waves in quasi-one-
dimensional non-steady oblique field magneto-gasdynamics, Acta Mech. 13 (1972) 255–263.
21. T. Hirschler and W. Gretler, Similarity analysis of strong converging spherical shock waves in
radiating gas, Acta Mech. 154 (2002) 159–177.
22. E. Haquea and P. Broadbridge, Exact solution of a boundary value problem describing the
uniform cylindrical or spherical piston motion, Appl. Math. Model. 35 (2011) 3434–3442.
23. M. Liverts and N. Apazidis, Limiting temperatures of spherical shock wave implosion, Phys.
Rev. Lett. 116 (2016) 014501.
24. C. C. Wu and P. H. Roberts, Structure and stability of a spherical shock wave in a Van der Waals
gas, Quart. J. Mech. Appl. Math. 49 (1996) 501–543.
25. N. Zhao, A. Mentrelli, T. Ruggeri and M. Sugiyama, Admissible shock waves and shock induced
phase transitions in a Van der Waals fluid, Phys. Fluids 23 (2011) 086101.
26. B. Bira, T. Raja Sekhar and G. P. Raja Sekhar, Collision of characteristic shock with weak
discontinuity in non-ideal magnetogasdynamics, Comput. Math. Appl. 75 (2018) 3873–3883.
27. R. Arora and V. D. Sharma, Convergence of strong shock in a Van Der Waals gas, SIAM J. Appl.
Math. 66 (2006) 1825–1837.
28. R. Arora, A. Tomar and V. P. Singh, Similarity solutions for strong shocks in a non-ideal gas,
Math. Model. Anal., 17 (2012) 351–365.
29. M. Pandey, B. D. Pandey and V. D. Sharma, Symmetry groups and similarity solutions for
the system of equations for a viscous compressible fluid, Appl. Math. Comput. 215 (2009)
681–685.
30. A. Sakurai, On the problem of a shock wave arriving at edge of a gas, Commun. Pure Appl.
Math. 13 (1960) 353–370.
31. P. Hafner, Strong convergent shock waves near the center of convergence: a power series solution,
SIAM J. Appl. Math. 48 (1988) 1244–1261.
32. A. Sakurai, Propagation of spherical shock waves in stars, J. Fluid Mech. 1 (1956) 436–453.
33. M. H. Rogers, Analytic solutions for blast wave problem with an atmosphere of varying density,
Astrophys. J. 125 (1957) 478–493.
34. G. Madhumita and V. D. Sharma, Propagation of strong converging shock waves in a gas of
variable density, J. Eng. Math. 46 (2003) 55–68.
35. A. Tomar, R. Arora and A. Chauhan, Propagation of strong shock waves in a non-ideal gas, Acta
Astronaut. 159 (2019) 96–104.

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 117 101–118
Copyedited by: ES MANUSCRIPT CATEGORY: Research article

118 A. CHAUHAN ET AL.

36. M. Van Dyke and A. J. Guttmann, The converging shock wave from a spherical or cylindrical
piston, J. Fluid Mech. 120 (1982) 451–462.
37. A. Chauhan, R. Arora and A. Tomar, Convergence of strong shock waves in a non-ideal
magnetogasdynamics, Phys. Fluids 30 (2018) 116105.
38. E. Haque, P. Broadbridge and P. L. Sachdev, Expansion of high pressure gas into air—a more
realistic blast wave model, Math. Comput. Model. 50 (2009) 1606.
39. D. S. Gaunt and A. J. Guttmann, Asymptotic analysis of coefficients, Phase Transitions and
Critical Phenomena, Vol. 3 (eds. C. Domb and M. S. Green; Academic, New York 1974)
181–243.

Downloaded from https://academic.oup.com/qjmam/article/73/2/101/5736526 by Amity University user on 20 April 2021


40. G. A. Baker and D. L. Hunter, Methods of series analysis II. Generalized and extended methods
with applications to the Ising model, Phys. Rev. B7 (1973) 3377–3392.
41. G. A. Baker, The theory and application of the Pade approximant method, Advances in
Theoretical Physics, Vol. 1 (ed. K. A. Brueckner; Academic Press, New York 1965) 1–58.

Appendix

Table A1 Neville Table for estimating the radius of convergence tc for m = 1, γ = 5/3, δ = 1
and b = 0.1

n e0n Linear Quadratic Cubic Quartic


27 0.46776519 0.46771172 0.46773619 0.46773647 0.46773532
28 0.46776334 0.46771347 0.46773620 0.46773633 0.46773548
29 0.46776168 0.46771504 0.46773621 0.46773629 0.46773604
30 0.46776017 0.46771645 0.46773622 0.46773626 0.46773608
31 0.46775880 0.46771772 0.46773622 0.46773623 0.46773604

Table A2 Neville Table for estimating 1/tc for m = 1, γ = 5/3, δ = 1 and b = 0.1

n e0n Linear Quadratic Cubic Quartic


27 2.005200408 2.13805336 2.13796479 2.13795607 2.13796087
28 2.009944928 2.13804697 2.13796392 2.137956671 2.13796024
29 2.014362041 2.13804119 2.13796319 2.13795683 2.13795787
30 2.018484505 2.13803595 2.13796256 2.13795695 2.13795771
31 2.022340849 2.13803118 2.13796203 2.13795707 2.13795788

Table A3 Values of tc , the time taken by the shock to collapse

γ b tc (m = 1, δ = 1) tc (m = 2, δ = 2)
7/5 0.01 0.62394129425 0.49943750532
7/5 0.1 0.54536888812 0.42805032338
5/3 0 0.63331233665 0.50825547344
5/3 0.01 0.53709321107 0.42033031820
5/3 0.1 0.47461503536 0.36655780676

[19:41 4/5/2020 OP-QJMA200003.tex] QJMAM: The Quarterly Journal of Mechanics & Applied Mathematics Page: 118 101–118

You might also like