Moscow state university by name M.V. Lomonosov, Research Institute for Mechanics, Moscow, Russia
Received Date: May 25, 2017; Accepted Date: July 15, 2017; Published Date: July 23, 2017
Citation: Tunik YV (2017) Instability of Contact Surface in Cylindrical Explosive Waves. Fluid Mech Open Acc 4: 168. doi: 10.4172/2476-2296.1000168
Copyright: © 2017 Tunik YV. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Fluid Mechanics: Open Access
In this paper are developed modifications of the Godunov scheme, based on Kolgan's scheme of the second order of accuracy in the spatial variables for smooth solutions. It is constructed schemes of the first and the variable order of approximation, which exceed the Godunov scheme in accuracy. Referencing to the system of differential equations for propagation of flat sound waves in a gas at rest, the Kolgan scheme and the first-order schemes obtained are investigated onto the ability to ensure the non decrease of entropy, that is, to product of physically justified numerical solutions. The test problems of nonlinear gas dynamics on the decay of a discontinuity in a pipe and the transformation of a non uniformity in a plane-parallel flow are solved. Cylindrical explosion task is considered as the main one. The stability of a contact discontinuity behind a blast wave is investigated numerically in the Cartesian and polar coordinate systems. Analysis of obtained and published solutions does not confirm the instability of the contact discontinuity which initially has the circular shape. Change of the shape of initially perturbed break is largely caused by the instability of Taylor, not Richtmyer-Meshkov. Calculations are partially fulfilled using supercomputer “Lomonosov” of Moscow state university.
Euler equations; Perfect gas; Blast; Non decreasing of entropy; Accuracy; Contact discontinuity; Richtmyer-Meshkov; Taylor instability
The explosion task was first posed and solved by Sedov  and independently by von Neumann  for a point initiation and shock waves of high intensity. In subsequent works various aspects of blasting processes are studied. The problem of propagation of blast waves are solved subject to backpressure and the finite size of the region of initiation. The history of the explosion theory is in detail expound [3-4]. Development of finite-difference numerical methods [5-21] have allowed to simulate discontinuous solutions of the equations of inviscid gas dynamics. Currently, the problem of the explosion is considered as a test for modern numerical schemes of a high order of accuracy. An overview of these methods can be found, for example [22-24]. Performed calculations taking into account the backpressure and the finite region of the blast initiation indicate the development of instability of the contact discontinuity . It has been suggested that this is a development of the Richtmyer-Meshkov instability .
In the present work for the solution of the problem of cylindrical blast taking into account the backpressure and the finite region of the explosion initiation there are used modifications of the scheme of S.K. Godunov, which increase the computation accuracy with respect to spatial variables. We numerically investigate the question of stability of the contact discontinuity.
The original Godunov's scheme  was presented and analyzed on the example system of equations describing the propagation of plane sound waves in resting gaseous medium:
Here - the density and speed of sound in the unperturbed environment; and - the perturbation of velocity and pressure.
The Godunov scheme for the system of eqn. (1) on the uniform computational grid with step h on the spatial coordinates is equivalent ratios ()
Here and below upper indexes refer to variables at the upper time layer, the subscripts relate to lower layer; CFL =τc0 / h - Courant number, τ - time step, Un+1/2 and Pn+1/2 values of u and p at the boundary between cells n and n+1:
The stability condition: CFL ≤ 1, determines the maximum time step at which perturbations generated at the boundary of the cell come up to the centre and back after the interaction with the oncoming perturbation.
In the case of non-linear equations of hydrodynamics, the magnitudes like Un+1/2 and Pn+1/2 represent flows of mass, momentum and energy, calculated using flow parameters at the boundary between neighboring cells. The required values of pressure, density and velocity are determined by solving the self-similar problem of break decay, and the Courant number is calculated by the maximum speed of propagation of perturbations .
The scheme of Godunov – Kolgan (SGK)
A first modification of the Godunov scheme is proposed by Kolgan , according to which the value of any parameter is attributed only to the center of the cell with number n. On the borders n+1/2 and n-1/2 parameter values are adjusted with account of linear distribution within the nth cell:
Here - the distance between centre and bound of the cell. The derivative in direction l is being calculated numerically and therefore ambiguously. Studies  proposed to use its minimum value for determining the quantities and , that is for any n on a uniform grid
Here min mod(A,B)=A, if |A| ≤ |B| , or else B, additionally min mod(A,B)=0, if A⋅B<0.
The Godunov scheme (2-3) under Kolgan modification is reduced to relations
In the case of non-linear equations of hydrodynamics the flows through cell borders are calculated by the way solution of the problem of break decay for the parameters corrected according eqn. (5). It should be noted that in this case use of self-similar solutions of the problem of break decay cannot be considered correct, because the parameters in cells are not constant. To take into account the change of parameters within the computational cell, strictly speaking, necessary to use the solution of the generalized Riemann problem . However, we can consider that in the Kolgan scheme cell is divided into two parts, in each of which a parameter has its constant value, distinct, in general case, from values in the cell center with the difference dφ|n. Then the self-similar solution of the problem of break decay is correct on half of the cell. Basing on the physical interpretation for the stability condition for the Godunov scheme, the Kolgan scheme must be stable when CFL ≤ 0.5.
Kolgan's and Godunov's schemes are attractive because clarity and simplicity, which gives them an advantage in modeling gas dynamic flows with nonequilibrium physicochemical processes. To generalize the Kolgan scheme we take the parameter values not at the border, but at intermediate point L, between the center and the boundary of the calculated cell. Then the increment of the unknown function will be specified by the magnitude , where α∈[0,0.5), and the relations (5) are written as:
When α=0 it is the Godunov scheme of the first order, and when α=0.5 it is the Kolgan scheme of the second order in the spatial variables. By analogy with the scheme of Godunov, the stability condition for the obtained scheme of Godunov-Kolgan (SGK) in the case of system (1) can be represented as CFL ≤ 1- α.
Entropy of the numerical solution in Godunov-Kolgan scheme
The non-decreasing entropy provides physically justified numerical solutions. In contrast to the Godunov scheme the corresponding problem exists in the case of schemes of a high accuracy and is discussed, for example, [16,26-28]. Here an ability of the Godunov- Kolgan scheme to satisfy the entropy condition is considered on the example of the system (1). The difference analogue of SGK for this system can be represented in the form
The residual term in the equation for the speed
From (10) it is clear that while α∈[0,0.5) the SGK has a first order of approximation, both in time and spatial variables. When α=0.5 this scheme has second order in spatial variables. If take in account equality
whose are consequence of the system (1), relations (9) can be rewritten
The solution of eqn. (11) can be considered as an approximate solution of the propagation equations of flat sound waves in the resting viscous gas with viscosity coefficient
The smaller the viscosity coefficient, the numerical solution is more close to the exact solution of system (1) for inviscid gas. If =0, the scheme has the second order of approximation with respect to both time and space for the linear system (1). Since the SGK have to guarantee the non-decreasing of entropy, the viscosity coefficient cannot be negative. Then from (12) it follows that the Courant number should satisfy the condition CFL≤(1−2α). When α=0 (Godunov scheme) this this inequality coincides with the condition of stability. If the intermediate point L approaches to the boundary of cell then allowable time step is reduced, because α→0.5. So condition of the non-decreasing entropy is an additional restriction on the time step for SGK. If α=0.5 the Courant number must be zero. Therefore, the Kolgan scheme does not guarantee the non-decreasing of entropy in numerical modeling of non-stationary processes, that may lead to physically unreasonable results.
Scheme with a variable approximation (SVA)
It is possible to consider SGK with a variable parameter α: .
Here any flow parameter on the lower time layer in the nth cell, β - the additional unit order parameter, which defines the reaction rate of the scheme to high spatial gradients χ. If spatial gradient , the scheme has first order (α≈0), when is small the scheme has the second order in the spatial variables (α≈0.5). So the obtained scheme is the scheme of a variable approximation (SVA).
Test tasks is being solved on the base of two-dimensional unsteady Euler equations for perfect nonviscous gas with constant ratio of specific heats γ=1.4.
The task of break decay in a pipe
It is being solved in rectangle -5≤x≤5 and 0 ≤ y ≤ 1 in the Cartesian coordinate system. At the initial moment pressure p and density ρ to the left of the line x=0 are equal unit: p1=1 and ρ1=1, on right p=p2=0.1, ρ=ρ2=0.125. The gas is at rest everywhere [29,30].
Computational grid has 100 × 10 cells. The calculation continues for the time during which a perturbation does not reach left and right boundaries of the computational domain, therefore parameters on the right and left border keep the initial values. On horizontal boundaries the normal derivatives of all parameters are set to zero.
The task can be solved in one-dimensional approximation; however, a two-dimensional formulation is used for checking the correctness of calculations in cells at the top and bottom borders.
As a result of a decay of high pressure region in low pressure region is propagated shock wave D, which is followed by density discontinuity C, onto left is propagated rarefaction wave RW (Figure 1a and 1b). At the maximum possible values of the Courant number for the stable accounts the Kolgan scheme (α=0.5), SGK (α=0.4) and SVA with β=1 give similar results in throughout computational domain. Loses in accuracy of the Godunov scheme is particularly noticeable in the regions of the contact discontinuity and of the rarefaction wave (Figure 1c and 1d).
On the other hand the solution obtained by the Godunov scheme is monotonic throughout the computational domain, in particular, in the vicinity of the point of initial discontinuity x=0 (Figure 2a). The Kolgan scheme leads to a non-monotonic parameter distributions in this area. The monotony of the scheme Kolgan rises with decreasing Courant number (Figure 2b), that corresponds to the findings in the first section.
Development of non-uniformity
Development of non-uniformity in a plane-parallel gas flow is considered in rectangle, which in the Cartesian coordinate system (x,y) is limited by direct lines y=0, y=1, x=0 and x=40. At the initial moment the pressure is equal to unit in whole region: p=p0=1, the projection of the velocity on the y-axis is zero, longitudinal component of the velocity u=1. The flow has uniform density: ρ=ρ0=1 everywhere except the limited region a<x<b, where ρ=ρ1 ≠ ρ0. Parameters on inlet do not change over time. On other boundaries the normal derivatives are set to zero.
It is obviously, that the exact solution can be obtained by parallel shift of the initial distribution along abscissa axis. In this task interest is represented by diffusion properties of computation schemes.
Two options are considered: ρ1=0.2 and ρ1=2. In both cases, the maximal smearing of the non-uniformity occurs when using a Godunov scheme (Figure 3a). SGK (with α=0.4) simulates the exact solution much more accurately than the scheme of Godunov. The Kolgan scheme gives the solution the most close to exact, which is providing almost constant density profile and preserving the extreme values almost at the initial level (Figure 3b). Similar properties has SVA at β=1 (Figure 3c).
Solution of the problem in Cartesian coordinate system
From the beginning the task is being solved in formulation, similar to ref. : in the Cartesian system of co-ordinates (x, y) are given two circles with centre at the origin. First of them with radius R1 limits the region of high pressure. Another circumference in the general case has the radius Rс ≥ R1 and sets the position of the density break in a rest gas at initial moment.
In ref.  there are presented the results of solving this problem at Rс=R1 with use six modern numerical schemes of high order approximation. The first four schemes have the second order and the last two the third order of accuracy. All indicate the development of instability of the initial contact discontinuity given by the circle. The distortion of the contact surface in ref.  is explained by the instability of Richtmyer-Meshkov . On the other hand, the development of the flow in a blast wave at breaking up of a high pressure cylindrical area is in detail described in ref.  just as in a lot of early papers on theory of an explosion without any mention about instability of contact discontinuity.
In this work it is considered that the pressure p=p1=1 in region r ≤ R1, and p=p0=0.1 if r > R1. The gas density is ρ=ρ1=1 for r ≤ RC and ρ=ρ0 when r > RC. The computational domain is a square in the first quadrant of the plane (x, y) with side a=2, which, as a rule, is divided into 400 equal parts. Calculations are made on the base of two-dimensional unsteady Euler equations. On the axes are set conditions normal flow lack. Parameter values at the other two boundaries are taken from the computational domain.
The first series of calculations is carried out at RC=0.25, R1=0.25 and ρ0=0.2. As a result of disintegration of high-pressure zone, in the domain with low pressure propagates a circular shock front, followed by a contact discontinuity. Towards the center goes a rarefaction wave, which after reflection continues to reduce the gas density behind the contact discontinuity. In the central part gas becomes over rarefied. That leads to the formation of a convergent shock wave. The reflection of this wave generates a secondary shock wave coming from the center. (Figure 4a-4c) corresponds to the moment when the leading front of the blast wave is at a distance of 1.8 from the center, and the secondary shock wave is already formed.
By Godunov scheme isochores are concentric circles (Figure 4a). Such development of the explosion agrees with description in ref. .
The results of calculations according to the scheme Kolgan substantially depend on the Courant number CFL. When CFL=0.5 (Figure 4b) isochores produce the chaotic and implausible picture. Monotonous flow is formed when CFL=0.1 (Figure 4c). The same result are obtained in the case of SGK (α=0.4) and SVA at CFL=0.25.
In this task the Kolgan scheme improves the flow picture not essentially. On the grid 800 × 800 the scheme of Godunov leads to very close results (Figure 5a-5d).
The gas expansion in the central zone ensures a rapid alignment of the density near the contact discontinuity, the initial break becomes invisible (Figure 5a). High density gradient on the (Figure 5b) corresponds to the secondary shock wave.
In the second series of calculations, the radius of the contact discontinuity is twice the radius region of high pressure: RC=0.5, R1=0.25 and ρ0=0.2. As a result of the decay of the high pressure region at the distance R1 from center is formed a second contact discontinuity 4 (Figure 6a and 6b). In contrast to the previous case rarefaction waves have not time to eliminate the primary break of density before the formation of the converging shock wave. In flow behind the leading shock front 1 is being formed a ring area of high density 3 (Figure 6a and 6b). Over time, the calculations according to the Kolgan scheme (CFL=0.1) the surface of the initial contact 2 becomes distorted (Figure 6a), like in ref. . Similar results give calculations with use SGK (α=0.4) and SVA (β=1).
However, it should be noted that when using a rectangular computational grid in a Cartesian coordinate system, any circle, in particular, that which defines the initial form of the contact discontinuity, is approximated non-uniformly along the contour. So each step of numerical integration automatically perturbs the contact surface. In fact, we have another formulation of the cylindrical explosion task, solution of which does not answer the question about the stability of the contact break. In the calculations according to the Godunov scheme contact discontinuity remains circular form (Figure 6b). Distortions of the contact surface are localized near the boundaries of the computation domain. In this case that may be the consequence of a bad choice of coordinate system and computational grid.
The solution of the problem in polar coordinate system
To test the hypothesis  about the instability of the contact discontinuity and to eliminate errors caused by the bad choice of the computational grid, the cylindrical blast problem is being solved in polar coordinate system: x=rcos(θ), y=rsin(θ). In the variables (r, θ) the region of computation is the circular sector with central angle π/2 and a radius equals 2. In the basic variant the computational domain radius is divided into 400 parts and the angle by 180 cells. On the axes θ=0 and θ=π/2 are set periodic boundary conditions. Such a construction of mesh eliminates the disturbance of the circular surface of initial discontinuity and ensures the symmetry of the flow with cylindrical waves.
The considered numerical schemes give almost identical results without any signs of the contact discontinuity instability, which initially is set as the circle with radius RC=0.5, while the radius of the high pressure region R1=0.25, ρ0=0.2 (Figure 7a).
In this regard, the interest is the case of the circular break with an initial sinusoidal perturbation (Figure 7b). In the calculations, the amplitude of the initial perturbation AMP=0.02. Impulsive interaction of the blast wave with the contact discontinuity does not alter the ratio of the amplitude and of perturbation wavelength of the contact line KS (Figure 8a). In the flow are formed conditions for rise of the perturbation amplitude: the pressure difference on the contact in the area of the concavity is smaller than near of the convexity (Figure 8b).
Figure 8:and b: Isochores (a) and isobars (b) on the pressure background after the passage of the shock front of a blast wave through initially perturbed contact discontinuity. c and d: Isochores (c) and isobars (d) on the background of pressure after changing of the sign of the contact line curvature.
As a result of interaction of incident shock front and the break the convexities and concavities of the contact line are reversed (Figure 8c) that is typical for the Richtmyer-Meshkov instability. However, in the rarefaction wave which follows the shock front, the mechanism of the Richtmyer-Meshkov instability does not actuate, because the pressure difference along the discontinuity line is quickly equated (Figure 8d).
However, all the schemes, including the Godunov scheme, indicate the growth of the amplitude of the perturbation. The weakening of the blast wave leads to a deceleration of the flow behind the leading front. Gas with high density from the inner side of the contact discontinuity locally "is injected" into the region of low density in front of the contact discontinuity line (Figure 9). It is result of actuation the mechanism of Taylor instability : the inertia of the dense gas replaces gravity.
Thereafter, the shape of the contact line is altered by damped waves of compression and rarefaction, which are periodically reflected from the contact surface and the axis of symmetry. On the discontinuity surface appear signs of propagation of longitudinal disturbances (Figure 10). But the gap retains a periodic structure. Over time, the change in the shape and amplitude of the perturbation on the contact surface practically stops due to the fall of velocity in the central part of the blast wave.
Thus, the cylindrical blast wave increases the initial perturbation of the primary contact surface, but to a certain limit. The reason for this gain is the Taylor instability.
Analysis of obtained and published solutions does not confirm the instability of the contact discontinuity which initially has the circular shape. Variation of the shape of initially perturbed break is largely caused by the instability of Taylor, not Richtmyer-Meshkov.
Schemes of first order developed in this work have more high accuracy than the Godunov scheme, and they produce results close to the Kolgan scheme second order in the spatial variables. Promising is the scheme, order of accuracy to which changes from the first to the second in the spatial variables depending on the gradient of flow parameters.
Despite the relatively high diffusion properties, the Godunov scheme successfully solves the problems of gas dynamics, in particular, simulates the development of instability of the contact discontinuity in the blast problem. This scheme guarantees the non decrease of entropy, and therefore remains an actual tool of numerical modeling along with modern schemes of high order accuracy.