1 Introduction

In many computational fluid dynamics (CFD) applications, the physical domain and flow have axisymmetric features. The process of solving such a three-dimensional problem is often quite expensive, and the two-dimensional cylindrical coordinate system is able to reduce the amount of work greatly by taking the advantage of symmetry. Therefore, the study of cylindrical symmetric compressible fluid flows is significant for the simulations of some implosion or explosion problems.

A large body of literatures involving the cylindrical coordinate calculation have been published, e.g., the Lagrangian method^[1-5], the Eulerian method^[6], and the arbitrary Lagrangian Eulerican method (ALE)^[7-9]. Different from the researches focusing on the grid moving strategies^[10-20], this paper pays more attention to the numerical formulations of the ALE method.

When extending the two-dimensional Godunov method in Cartesian coordinates to the cylindrical geometry, there are always some extra issues associated with the cylindrical algorithm. The first issue is about the solution singularity when a shock wave hits the symmetrical axis. Due to the lack of a strict mathematical theory in this aspect, the singular behavior stemming from the geometric source has not be solved completely in theory. However, in the numerical calculation, some scholars have tried to do some studies. Matsuo et al.^[21] simulated the cylindrically converging shock and demonstrated the behavior when the wave approached the axis. Noh^[22] took note of the strong “wall heating” phenomenon when a shock reflects from the symmetry axis and introduced the heat viscosity for this problem. Li et al.^[23] designed a new boundary condition in the symmetric center for a one-dimensional Lagrangian method and hoped to remove the singularity in their solver of the generalized Riemann problem (GRP). The second issue involves the relation between the spherical symmetry preservation and the conservation of a numerical scheme. Here, the spherical symmetry preservation means that a two-dimensional cylindrical scheme should keep the strict symmetry for all physical quantities when a one-dimensional spherical problem is computed on an equal-angle-zoned grid. It is also the key point of this paper and needs to be discussed furthermore.

It is well-known that there are some different discretized manners when cylindrical coordinate problems are dealt with. One is to take the cylindrical problem as a purely plane problem with the addition of a geometric source term (see Refs. [24]-[25]), and the corresponding numerical scheme is constructed for the plane problem. This approach has an advantage that the existing numerical solver for the planar problem can be used directly without modification, and the source term can be incorporated into an ordinary differential equation (ODE) system directly. However, some authors regard that this method has a severe drawback, i.e., the original multidimensional conservation property has lost. Another way to derive the discrete algorithm is from a three-dimensional control volume cell and finally construct a two-dimensional numerical scheme based on cylindrical elements. The method preserves the strict conservation in the three-dimensional sense but loses the one-dimensional spherical symmetry. Usually, it is denoted as a ‘volume weighted scheme’ or a classical ‘control volume scheme’ in literatures (see Refs. [3]-[6]). In order to restore the symmetry property, some different methods have been proposed and developed. Margolin and Shashkov^[26] used the ‘curvilinear grid’ to connect the nodes so that the spherical symmetry could be exactly preserved even on unequal-angle-zoned grids. Cheng and Shu^[3] proposed a modification to the classical control volume method in which the source term was replaced by an average of two contact pressures on two specific edges. Zabrodin et al.^[9] changed the initial values in a one-dimensional Riemann solver to let the force acted on the azimuthal direction be zero, and thus preserved the spherical symmetry eventually. There are still other methods^{[5, 24]}. Among all these approaches, the most widely used one is the so-called area weighted difference^{[5, 27-28]}, which usually occurs in the context of a staggered Lagrangian frame and preserves the spherical symmetry strictly.

In this paper, the area weighted scheme is extended from a Lagrangian formulation to an ALE one. The method is suitable for the arbitrary grid, and its spherical symmetry property can be proved strictly. In addition, the aforementioned three kinds of ALE formulations are compared in numerical experiments to illustrate which one is better when the strong shock wave is simulated.

The outline of this paper is as follows. In Section 2, the governing equations written on the cylindrical geometry are given. In Section 3, the traditional ALE methods including the control volume scheme and the plane scheme are described. The spherical symmetry property is discussed in Section 4. An area weighted scheme for the ALE method is proposed in Section 5, and the numerical results are illustrated in Section 6. Finally, the conclusions are summarized in Section 7.

2 Governing equations and discretization in space

The governing equations of the inviscid flow in two-dimension are as follows:

(1)

where the vectors of the conserved variables and flux vectors are

where ρ, p, and E are the fluid density, the pressure, and the total energy, respectively, and u=(u, v) is the fluid velocity. The equation of state is in the form of

where γ is the specific heat ratio.

In many literatures, the system of Eq. (1) is also written as

(2)

which has exactly the same form as a purely plane problem, but with the addition of a geometric source term on the right-hand side,

(3)

It is convenient, for the subsequent discretizations, to recast Eqs. (1) and (2) in the more fundamental control volume formulation^[29], which holds for an arbitrary moving control volume,

(4)

(5)

where w is the moving velocity of the control volume Ω(t) and is defined at the boundary ∂Ω(t). N is the unit outward normal direction on the boundary of Ω. If w=u, the system can be reduced to a Lagrangian formulation, and if w=0, it has an Eulerian form.

If we introduce the material derivative , then we have the differential equation with the Lagrangian form, i.e.,

(6)

where τ=1/ρ is the specific volume. Notice that the momentum equation in Eq.(6) poses a Cartesian form if the factor r on the both sides is eliminated simultaneously.

Remark 1 Although there exist different forms for the inviscid axisymmetry problem, Eqs.(1), (2), and (6) have the same weak solution, i.e., the solutions in the smooth region are the same, and the jump conditions across a discontinuity are the same too. The RankineHugoniot condition is

(7)

where D is the propagation speed of a discontinuity, and L and R represent the left and the right sides of the discontinuous interface, respectively.

Proposition 1  For any unit vector N=(n_x, n_r)^T, two-dimensional fluxes satisfy the following rotational invariance property:

(8)

where

(9)

Please see Ref. [25] for the proof.

3 Discrete schemes 3.1 Notations on generic polygonal grid

On the xr-plane, each cell A_c of the mesh is assigned a unique index c. We use the index f to denote a generic neighboring cell A_f which has a common edge denoted as k, and N_k=(n_x, n_r)k is the unit normal vector of the cell A_c on the edge k (see Fig. 1). Each vertex of the mesh is assigned a unique index q, and and are defined to be the sets of vertices and edges of A_c, respectively, i.e.,

The physical quantities such as the density ρ_c, the pressure p_c, the velocity u_c, and the energy E_c are defined in the cell center of A_c. The moving velocity w_q of the grid is defined on the node q.

Some geometry quantities such as the volume and the area can be given as follows:

(10)

(11)

where (x_q, r_q) is the coordinate of the node x_q.

3.2 Control volume schemes

We discrete Eq. (4) by utilizing the idea of the Godunov’s scheme on the unstructured mesh,

(12)

where N_k is the unit outward norm of the edge k between A_c and A_f. U_k^w and F_k^w are the rotated conservative quantity and flux on the edge k, respectively. L_k is the length of the edge, and is the weighting function at the cell face k, i.e.,

where w_k is the moving velocity of the cell edge k, which is determined by the nodal moving velocities w_q and w_q⁺, and w_k is its normal projection, i.e.,

(13)

In this semi-discrete scheme, we assume that the conservative quantities and fluxes U_k^w and F_k^w are constants along the edge k. Furthermore, we assume that all quantities appearing in conservative and flux functions are set at the time tⁿ for simplicity. It means that the above geometric quantities L_k and are also defined at the time tⁿ.

The source term S_vdxdr has a variety of discrete forms, but the simplest one is

(14)

where denotes the area of the cell A_c, and is the cell centered source term at the time t_n. Then, the fully discrete scheme is

(15)

where H_k^w=F_k^w −U_k^w w_k is the flux cross the moving edge k. The update of the mesh depends on the nodal moving velocity,

(16)

Remark 2  The discrete scheme (15) is called ‘the control volume’ or ‘the volume weighted’ scheme. It can be understood as a limit case of a finite volume scheme on a specific threedimensional control volume. Thus, the scheme is conservative strictly in the three-dimensional meaning.

Similarly, a purely plane algorithm with a geometric source term can be expressed as

(17)

where the geometric source term is the cell center value at the cell A_c, and r_cⁿ is the r-direction coordinate at the cell geometric center.

3.3 Approximate Riemann solver

The numerical conservative quantities and fluxes U_k^w and F_k^w between the cells A_c and A_f can be obtained from many methods. The most popular way is to solve a one-dimensional Riemann problem along the normal direction of the edge k, i.e.,

(18)

the initial values are

(19)

and the numerical flux in Eq. (15) can be expressed by

(20)

where U_k^w=U ^w(U_L, U_R) and F_k^w=F^w(U_L, U_R) are the Riemann solutions of the state and the flux function along the grid velocity x/t=w_k, respectively. In this paper, we use a robust approximate Riemann solver Harten-Lax-van Leer-contact (HLLC)-2D^[30], which has a better behavior than the classical HLLC solver^[31].

4 Spherical symmetry preservation on equal-angle polar grid

The spherical symmetry preservation is a very important property when one-dimensional spherical problems are computed on an equal-angle-zoned grid. It is because the small deviation from the spherical symmetry may be exaggerated by some instabilities and induce very large errors. Unfortunately, the control volume scheme (15) cannot preserve the spherical symmetry. Without loss of generalization, we use the Eulerian method with a fixed grid to demonstrate this conclusion. Now, the momentum equation in the numerical scheme (15) can be reduced to

(21)

where ρ_k^w, u_k^w, and v_k^w denote the states on the cell interface k.

With an equal-angle-zoned grid shown in Fig. 2, the normal vectors of the cell boundaries are

where the subscripts r, l, t, and b are neighboring cell indices of the cell A_c, and correspond to the right, left, top, and bottom, respectively. ø is the angle between the radial direction passing through the cell center and the x-axis, and the θ is the angle between two adjacent radial lines. Other geometric parameters are as follows:

where and are the radii of the cell nodes 1 and 2 in Fig. 2.

The volume and the area of the cell A_c are

For a one-dimensional spherical flow, the initial density, the pressure, and the total energy are only functions of the radius, and the velocity is in the radial direction, i.e.,

where T_r is a unit vector perpendicular to N_r.

On the four edges of the cell A_c, from a classical approximate Riemann solver, we have

(22)

Multiply Eq. (21) with the vector T_r. Then, after some calculations, there is

(23)

Obviously, the symmetry property of the control volume scheme (21) is preserved if and only if the following condition holds:

(24)

In a classical approximate Riemann solver, usually, . It is because of the presence of the normal velocity between the edge 23 or 34. Under the circumstance, the contact pressure p is always larger or smaller than p_cⁿ. Thus, the control volume scheme cannot maintain the one-dimensional spherical symmetry.

In order to preserve the spherical symmetry strictly, a simple modification^[3] is to replace the pressure p_cⁿ in the source term by (p_t^w + p_b^w)/2. At this time, the symmetric condition (24) holds naturally when the calculated problem is a spherical one-dimensional one. Note that the discrete momentum equation (21) now can be rewritten as

However, the drawback of the method is also obvious owing to its dependence on the polar grid. A general grid cannot distinguish the azimuthal direction, for example, the method does not work on any unstructured grid.

If we hope to find a numerical method which works on any unstructured grid and keeps the spherical symmetry on a polar grid, then the area weighted scheme is a good choice.

5 Area weighted approach

In the Lagrangian framework, the most widely used method in the cylindrical problem is the area weighted method. In this approach, the momentum equation in Eq. (6) is integrated on the area element dxdr rather than the volume element rdxdr due to its Cartesian form. The discrete formulation of the momentum equation is

(25)

However, when the numerical method is extended to the Eulerian and ALE method, the way to perform an integration on the area element is not easy to implement since the momentum equation no longer has the above Cartesian form. But if we view the area weighted method as a discrete technique to the source term, just as pointed out in Ref. [5], then we may readily generalize the method to an ALE formulation.

Rewriting Eq. (12), we have

(26)

where the geometric source is represented as

(27)

Notice that , and two terms on the left side can be approximated as

where p_k^w is the pressure flux on the moving interface k.

Now, the source can be expressed as

(28)

where is an argument normal vector. Thus, the resulting area weighted scheme is

(29)

Remark 3  When Eq. (29) is degenerated to the Lagrangian formulation, the momentum equation is discretized as

(30)

Compared with Eq. (25), the two area weighted schemes are not completely the same.

Remark 4  The discretization (28) of the source S_{v, 2} is completely different from that of the source term S _v. One might question whether the split form of Eq. (27) breaks the RankineHugoniot jump condition of the original problem, since it creates non-divergence product terms of physical quantities. However, this kind of worry is needless. In fact, the split formula p=(rp)_r − rp_r is a linear representation, especially the coefficient r is continuous. Thus, the split does not destroy any discontinuity relation. Further, it does not influence the weak solution of the original problem.

For a one-dimensional spherical flow, we may prove that the area weighted scheme (29) maintains the spherical symmetry strictly if the equal-angle polar grid keeps a radial motion. In fact, the discrete momentum equation is

(31)

Notice that the symmetrical polar grid, at this time, keeps a radial movement, which means w_t=w_b=0. After some algebraic calculations, we have

Obviously, the obtained velocity is along the radial direction. Similarly, we can prove that ρ_cⁿ⁺¹ and E_cⁿ⁺¹ are independent of ø. In summary, we have the following proposition.

Proposition 2  The area weighted scheme (29) is spherically symmetrical.

6 Numerical experiments

In this section, we present three test cases to assess the performances to the control volume scheme, the area weighted scheme, and the plane plus geometrical source scheme.

6.1 Spherical Noh test case

The spherical Noh problem^[22] is a challenging test case. The initial dimensionless state is

where ø=arctan(r/x) is the polar angle between the radial direction passing through the point and x-axis. The computational domain is a quarter circle with a radius of 1. The initial grid is an equi-angular polar grid, and the grid moving manner is determined by a least square approach^{[29, 32]}, which can be regarded as an approximate Lagrangian method. We use the first-order scheme to calculate this problem and obtain the final grids at t=0.6 (see Fig. 3). Obviously, the symmetry preservation of the control volume scheme is broken, but the symmetry from other two discrete formulations is kept rather well.

6.2 Spherical Sod problem

The spherical Sod problem is a classical example. The computational domain is defined in the polar coordinates by , where , and ø=arctan(r/x) (see Fig. 4). The initial dimensionless values are

We compare the three schemes in an ALE framework. The adaptive moving grid strategy comes from Ref. [11]. The initial grid and the finally generated grid from the area weighted scheme are illustrated in Fig. 4. The other two resulting grids obtained from the planar scheme and the control volume scheme are almost the same. The densities are illustrated in Fig. 5, in which each cell value is plotted. Note that we have used a second-order MUSCL reconstruction technique to improve accuracy, where MUSCL stands for a monotone upstream centered scheme for conservation laws^[33]. The notation ‘2nd plane ALE’ means that the result is obtained from the second-order plane scheme with a geometrical source on the adaptive moving grid, and ‘2or area ALE’ and ‘2nd volume ALE’ correspond to the second-order area weighted scheme and the volume weighted scheme, respectively. The results show that the differences among the three formulations are little for such a problem with the moderate strength shock.

6.3 Sedov problem in cylindrical coordinate system

We consider the Sedov problem for a point blast in a uniform medium with the spherical symmetry. An analytic solution was obtained by Sedov^[34] using self-similarity arguments. We model this problem on a domain [0, 1.2] × [0, 1.2] in a dimensionless form. The initial state is

A finite source internal energy at the origin is e₀=0.425536. It will produce an infinite strength shock that is located at the radius 1 at the time t=1.

In this test example, three formulations are firstly compared by use of the Eulerian method. Figure 6 shows the density contours using a second-order MUSCL reconstruction on a triangular mesh. As seen here, although the triangular grid is no longer symmetric, the density contours obtained from the three schemes still keep the spherical symmetry characteristics rather well. The shock fronts from the area weighted scheme and the control volume scheme are correct, while the shock position from the ‘plane plus source term’ scheme is wrong. In order to further observe the behaviors of three schemes, let us see what would happen in the Lagrangian method. Figure 7 compares three different grids at t=1. Once more, both the area weighted and control volume schemes have the correct shock wave speed whereas the ‘plane plus source term’ scheme propagates faster than the exact solution. The reason to produce such phenomena is not yet clear. However, based on the observation to the ‘plane plus source term’ scheme in Fig. 7, we can find that the error (including the grid distortion) near the symmetry axis is larger than that in the vicinity of the r-axis. The result shows that the error might come from the influence of the singular source term. The further discussion is beyond the scope of this paper.

7 Conclusions

In this paper, we pay attention to the spherical symmetry property, extend the concept of the area weighted scheme from the classical Lagrangian schemes to the ALE method, and remain the area weighted scheme in the framework of control volume algorithms. Based on the work, three different discrete formulations are described. They are the control volume method, the area weighted scheme, and a plane approach with the addition of source term. Numerical examples are devised to compare and contrast the performances of the three schemes.

All schemes are applicable for problems with the moderate strength shock. However, when there exists the strong shock wave in the vicinity of the symmetry axis, the plane plus source term method may suffer from the influence of the singular source, whereas the control volume scheme and the area weighted scheme are free of this flaw.

When one prefers first-order accurate numerical methods, the area weighted scheme is a good choice because it has the better symmetry than the control volume one. For secondorder accurate algorithms, the control volume and area weighted schemes usually have similar numerical behaviors.