A symmetry preserving scheme for three-dimensional Lagrangian radiation hydrodynamic simulations of ICF capsule implosion

Computers and Fluids 195 (2019) 104317

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

A symmetry preserving scheme for three-dimensional Lagrangian radiation hydrodynamic simulations of ICF capsule implosion Shaodong Guo, Mingyu Zhang∗, Haibing Zhou, Jun Xiong, Shudao Zhang Institute of Applied Physics and Computational Mathematics, Beijing 100094, China

a r t i c l e

i n f o

Article history: Received 2 July 2019 Revised 24 September 2019 Accepted 9 October 2019 Available online 10 October 2019 Keyword: Symmetry preserving Radiation hydrodynamics 3D Cartesian coordinate system Mimetic finite difference method Lagrangian hydrodynamics

a b s t r a c t A symmetry preserving scheme for three-dimensional (3D) Lagrangian radiation hydrodynamic simulations of an inertial confinement fusion (ICF) capsule implosion is presented. Violation of spherical symmetry exists in conventional scheme in 3D Cartesian coordinate system. Following the idea of Caramana et al. (E.J. Caramana, et al., A Compatible, Energy and Symmetry Preserving Lagrangian Hydrodynamics Algorithm in Three-Dimensional Cartesian Geometry, Journal of Computational Physics, 157 (20 0 0) 89–119), a symmetry preserving scheme in 3D Cartesian coordinate system is extended to Lagrangian radiation hydrodynamics. This scheme is specially designed for spherical (radial) grids. Appropriate calculations of auxiliary points in the discretization of Lagrangian hydrodynamic equations are given. Small modification is applied to the flux-related Jacobian matrix in the diffusion scheme. The modified diffusion scheme can guarantee the symmetrical flux. The symmetry preserving properties of the developed scheme for Lagrangian radiation hydrodynamics are proved in theory. The developed symmetry preserving scheme is applied to the 3D simulations of an ICF capsule implosion. Numerical results show that the symmetry preserving scheme for Lagrangian radiation hydrodynamics helps to remove unphysical instability in the 3D simulations of ICF capsule implosion. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The Support Operator Method (SOM), also named Mimetic Finite Difference method (MFD), is widely used to solve diffusion equations and Lagrangian gas dynamics equations [1]. The SOM has discrete approximations that preserve important properties of continuum equations. For the diffusion equation, the discrete flux operator and divergence operator are carefully designed to satisfy duality relations using the SOM. This diffusion scheme [2,3] is second-order accurate on both smooth and non-smooth meshes either with or without material discontinuities, although it has a non-local stencil for the diffusion operator. A significant development is the introduction of the local support operator method [4], where both cellcentered and face-centered unknowns are used to discretize the flux operator. The new local SOM has a local stencil and yields a sparse matrix representation for the diffusion operator. Many extensions to this approach have been developed recently. It can be found that the approach is used to construct mimetic discretization on general polygonal meshes [5], arbitrary polyhedral meshes

∗

Corresponding author. E-mail addresses: [email protected] [email protected] (M. Zhang). https://doi.org/10.1016/j.compfluid.2019.104317 0045-7930/© 2019 Elsevier Ltd. All rights reserved.

(S.

Guo),

[6], meshes with local refinement [7], and non-conformal meshes [8]. For the Lagrangian gas dynamics equations, using the SOM, the discrete divergence of the velocity is derived to be consistent with the change of volume of a fluid element to satisfy the geometric conservation law, and the staggered compatible discretization of the momentum and internal equations is designed to satisfy the conservation of the total energy [9]. This compatible or mimetic discrete scheme is also applied to solve Lagrangian gas dynamics equations on general polygonal meshes [10], unstructured polyhedral meshes [11], including sub-zonal pressure forces [12], and artificial viscosity forces [13,14]. For more details of the SOM or MFD, a complete history and main ideas of the SOM are comprehensively reviewed in the reference [1]. Since there is no symmetrical grid for a sphere in 3D Cartesian coordinate system, small departure from 1D spherical symmetry always exists. The small departure will be amplified in the 3D simulations of an ICF capsule implosion. Numerical results will be polluted by the unphysical numerical instability. A symmetry preserving scheme helps to remove numerical errors due to the amplification of small departure from spherical symmetry. And it improves our predictive capability of numerical methods in the 3D simulations of an ICF capsule implosion. It is extremely important to keep one-dimensional (1D) spherical symmetry in two dimensional (2D) or three dimensional (3D)

2

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

Cartesian coordinate systems during numerical simulations of inertial confinement fusion (ICF). Research on symmetry preserving scheme has focused on Lagrangian hydrodynamics in 2D and 3D and Arbitrary Lagrangian-Eulerian (ALE) in 2D. Caramana [15] et al. presented symmetry preserving scheme which preserve cylindrical and spherical symmetry in Cartesian and cylindrical coordinate system. Symmetrical results can be obtained on the grid with unequal angular zoning. Further Caramana [16] et al. extended the 2D symmetry preserving scheme to 3D. The developed scheme exactly preserves 1D spherical symmetry in 3D Cartesian coordinate system. The scheme is suitable for both structured and unstructured grids. Margolin [17] et al. developed symmetry preserving discretizations for Lagrangian hydrodynamics by using a curvilinear grid. A symmetry preserving version of artificial edge viscosity was given. Cheng [18,19] et al. proposed a cell-centered Lagrangian scheme which preserves 1D spherical symmetry in 2D cylindrical coordinate system and conserves mass, momentum, and total energy. Velechovsky [20,21] et al. presented symmetry- and boundspreserving momentum remap for staggered ALE methods. Kenamond [22] et al. proposed a compatible, total energy conserving and symmetry preserving algorithm for 2D ALE hydrodynamics on general polygonal meshes. All these previous researches concentrate on Lagrangian hydrodynamics. The ICF capsule implosion is an extremely complex multi-scale and multi-physics process which can be described by the radiation hydrodynamics equations. In order to reduce the effects of small departure from 1D spherical symmetry due to discrete errors, it is necessary to study a symmetry preserving scheme for 3D radiation hydrodynamic equations for the ICF capsule implosion. In current work a symmetry preserving scheme for radiation hydrodynamics in 3D Cartesian coordinate system is developed following the idea of Caramana [15,16] et al. The developed scheme exactly preserves 1D spherical symmetry in 3D Cartesian coordinate system. The ability of the developed scheme to preserve symmetry is proved in theory. And 3D numerical applications in ICF capsule implosion are given. It should be pointed out that this scheme is specially designed for spherical or radial grid, not for general grids. This paper presents the newly developed symmetry preserving scheme for 3D simulations of an ICF capsule implosion in details. In Section 2, the governing equations of Lagrangian radiation hydrodynamics in ICF capsule implosion are given. The SOM to discretize the radiation hydrodynamics equations is presented. The main idea of the developed symmetry preserving scheme for radiation hydrodynamics is shown. The details of the developed scheme are proposed. In Section 3, 3D numerical tests for an ICF capsule implosion are given. The ability of the developed scheme to preserve the symmetry is shown. Numerical tests include hydrodynamically equivalent implosion, one-temperature radiation hydrodynamics with and without small perturbations. Finally, the concluding remarks with the future work are given. 2. Numerical algorithm 2.1. The equation of radiation hydrodynamics: one-temperature diffusion approximation In this paper, the diffusion approximation to the radiation transport is used in the radiation hydrodynamics equations. The diffusion approximation, called one-temperature approximation, is considered in our numerical models. In the one-temperature approximation model, it assumes that the electrons and ions are in equilibrium at the same temperature (the material temperature). The conservation equations for mass, momentum and total energy are

dρ = −ρ∇ · v dt

(1)

ρ

dv = −∇ (Pm + Pr ) dt

(2)

ρ

d ( em + er ) = −∇ · Fr − (Pm + Pr )∇ · v + WL . dt

(3)

In Eqs. (1)–(3), ρ is the density, v is the fluid velocity, (Pm , em ) and (Pr , er ) are the pressure and energy for the material and radiation respectively, WL is the energy source and Fr is the radiation heat flux

Fr = −κr ∇ Tr

(4)

where κ r is the conductivity. The thermodynamic quantities are derived either from the ideal gas model or from the data of realistic equation of state. 2.2. MFDM for the above radiation hydrodynamics The operator splitting method is used to solve the onetemperature equations. First a hydrodynamics step is done to solve Eqs. (1)–(3) while ignoring the diffusion term, the energy exchange term and the energy source term. In this step, the staggered compatible hydrodynamics algorithm is used to discretize the momentum and energy equations.

Mp

 p dv p = fz dt

(5)

 d ez =− f pz · v p dt

(6)

z∈S( p)

Mz

p∈S(z )

where vp is the velocity defined at the grid vertex, and ez is the energy defined in the cell center. fp z or fz p is the corner force deriving from the pressure, the artificial viscosity and the sub-zonal pressure. The summations in the right-hand terms are performed over the stencils of the zone or point. The stencil S(p) for a point is each zone that contains the point, and the stencil S(z) for a zone is each point that is a vertex of the zone. Mp is the point mass and Mz is the zone mass or cell mass, both of which are defined using the corner mass mp z and mz p as follows:



Mz =

mzp

(7)

mzp .

(8)

p∈S(z )

Mp =

 z∈S( p)

For details of the calculations of the corner mass and corner force, one can be found in the references [9,12,14,16]. After the hydrodynamics step, an intermediate energy or temperature from EOS is obtained. Next, with the intermediate temperature and coefficients such as conductivity, a diffusion step is started to simulate the heat diffusion and energy deposition due to the diffusion term, the energy exchange term and the energy source term in Eqs. (1)–(3). In the diffusion step, the equations for the temperature can be written as the following generic form

ρCv α

∂ Tα = −∇ · Fα + Wα ∂t

(9)

where α can be e, i or r, denoting the electrons, ions or radiation field. Cvα , Tα , and Wα are the corresponding heat capacity, the temperature, and the source term, respectively. Fα is the heat flux, defined as

Fα = −κα ∇ Tα where κ α is the corresponding conductivity.

(10)

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

3

In this step, the local SOM is used to discretize the diffusion Eq. (9). Following the idea of the local SOM [4,6], the divergence operator is discretized as the prime operator

∇ ·F =

s 1 f k Ak V

(11)

k=1

where V is the volume of the cell, Ak is the area of the kth face of the cell, and fk is the face-normal component of the flux F, which is located at the centers of the kth face. fk can be calculated as

f k = F · nk .

(12)

Then according to the adjoint relationship of the discrete operators, the discretized flux vector F˜ is derived as follows

F˜ = G−1 AT˜

(13)

where F˜ is expressed in terms of the face-normal components of the flux as F˜ =[ f1 , f2 , · · · , fk , · · · , fs ]T . In the right hand side, T˜ =[Tc − T1 , Tc − T2 , · · · , Tc − Tk , · · · , Tc − Ts ]T , where Tc is the temperature defined at the cell center, T1 , T2 , , Tk , , Ts are the temperatures defined at the face center, A =[A1 , A2 , · · · , Ak , · · · , As ]T , and G is given by

G=



PnT Jn−1 K −1 Jn−T PnVn

Fig. 1. Initial grid construction for 1D spherically symmetrical problem.

(14)

n

In the Eq. (14), the summation is performed over all nodes of the cell. Jn is the Jacobian matrix given by Jn = [nf 1, nf 2, nf 3 ], where f1 , f2 and f3 are the faces adjacent to the node n, and nf 1, nf 2, nf 3 are the outward unit normal vectors of these faces. Pn is the permutation matrix for the node n, which is used to convert the short vectors involving the faces adjacent to the node n into the long vectors involving all of the faces of the cell. Vn is the volume of the corner cell containing the node n. Substituting the discrete divergence operator (11) and flux operator (13) to the Eq. (3), we can finally obtain the spatially discrete diffusion equation for cell centers. Next, the continuity conditions of the temperature and the normal component of the flux are exerted across the cell interfaces. The continuity conditions require that c2 Tkc1 1 = Tk2

(15)

fkc11 = − fkc22

(16)

where the cell c1 shares its k1 th face with the k2 th face of the cell c2 . With the discrete diffusion equations and the continuity of flux equations, we obtain the closed system of equations for cell-center and face-center temperature which can be solved using many efficient iterative methods. In the numerical simulation, we use the BoomerAMG (a parallel implementation of the algebraic multigrid method) preconditioner and conjugate gradient solver by the HYPRE library [23]. After the diffusion step, the final temperature can be obtained. Then the final energy and pressure can be computed from the density, temperature and the EOS at the end of the time step. 2.3. Symmetry preserving scheme 2.3.1. Statement of symmetry: mathematics and numerics The problem of symmetry in mathematics requires that the initial conditions (ICs) and the boundary conditions (BCs) are symmetrical. Therefore, the corresponding numerical system is required to have symmetrical ICs and BCs. Besides these restrictions, the initial mesh has to represent the symmetry, such as 1D spherical symmetry. As mentioned in [16], the initial mesh construction is not completely free and has to be specified according to the

problem of symmetry. Fig. 1. shows the initial radial mesh construction for the problem of 1D spherical symmetry in 3D Cartesian coordinate system. Any nodes can be represented by (i, j, k). The index i means that the distance from the origin O (the sphere center) is ir. The index j means that the polar angle is jθ . And the index k means that the azimuthal angle is k. The mesh in Fig. 1 has the following properties: a) the edges represented by the same j and k are along the radial direction; b) the nodes represented by the same i are on the same sphere surface; c) the cells having the same direction and solid angles are equal. In the current method, the density, pressure, and specific internal energy are defined in cell centers, and the coordinate and velocity are defined at grid nodes. The initial mesh in Fig. 1 has the property of 1D spherical symmetry. To keep the spherical symmetry, the velocity must be along the radial direction and the variation of the cell variables must be uniform between two i-surface. Usually it is required that the pressure gradients are along the radial direction and have the uniform magnitude between two isurface. The radiative fluxes have the same requirement. In mathematics, the increments of the variables (the velocity, density, pressure, and specific internal energy) and discrete radiative fluxes are required to only depend on the radial coordinate. 2.3.2. Symmetry preserving scheme for Lagrangian hydrodynamics In Lagrangian hydrodynamics, the gradient operator is discretized as below

∇p =

 1 Al (−pi + pi−1 ) Mp

(17)

where Mp is the nodal mass, pi is the cell pressure, and Ai is the surface area vector. Usually the sum of the surface area vectors Ai is not guaranteed to be along the radial direction. Therefore, the gradient operators cannot preserve 1D spherical symmetry. According to [16], the sum of the surface area vectors is modified to keep the component along the radial direction. The modified gradient operator can preserve 1D spherical symmetry. The calculations of auxiliary points which include the centers of edge, face, and cell are very important. From the Eq. (17), the ratio of the surface area to the corresponding volume is required to be uniform between i-surfaces. So the auxiliary points are required to be similar between i-surfaces. In Section 3.2, numerical simulation

4

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

Fig. 2. The grid for the flux-related Jacobian Matrix.

shows that inappropriate calculations of auxiliary points result in destroy of 1D spherical symmetry. 2.3.3. Symmetry preserving scheme for the diffusion equation The symmetry preserving scheme for the diffusion equation has the same requirement for the mesh construction as that of Lagrangian hydrodynamics. And also ICs and BCs have to be 1D spherically symmetrical. The radiative fluxes are required to be along the radial direction and have uniform magnitude between i-surfaces. Since the diffusion scheme is implicit, the linear system for the diffusion problem should have the following property. The equations along i-surface are equivalent. Usually the linear system for the diffusion problem with symmetrical ICs and BCs cannot guarantee the 1D spherical solution, even though the initial mesh construction is 1D spherically symmetrical. Here we show that this problem can be solved by simple modification of the flux-related Jacobian matrix in the Eq. (14), which will be given in the following part. The modified linear system will preserve 1D spherical symmetry on the spherical mesh. The Jacobian matrix for the node 1 in Fig. 2 is given as below



J1 = nle f t ,n f ront ,nlower and



J1T J1 =

1 n f ·nle n f ·nlo



n f ·nle 1 nle ·nlo

(18) n f ·nlo nle ·nlo 1



Fig. 3. The error as a function of the cell width via (a) the symmetry preserving scheme and (b) the conventional support operator scheme.

(19)

where n is the face normal vector, and f, le, lo represent front, left, lower, respectively. In order to obtain the radiative fluxes, which are along the radial direction, the Jacobian matrix is modified in the following way. Only the face normal vector nlo is modified as

nlom = (nlo · n51 )n51

(20)

where n51 is the vector from node 5 to node 1. Then nlo in the Eq. (19) is replaced with nlo m , and the Eq. (19) is modified as



J1T J1

m



=

1 n f ·nle 0

n f ·nle 1 0

0 0 |nlo · n51 |2



(21)

The same manipulation is applied to the calculations of all the radiative fluxes.

Finally, the variation of the cell temperature can be obtained as

T =



C1 R4up K

ρ cv R3up − R3lower +

2 (T − Tup )t

C2 R4lower K



ρ cv R3up − R3lower

2 (T − Tlower )t

(22)

where H is the height as shown in Fig. 2, C1 , C2 are constants, ρ is the density, cv is the capacity, and K is the conductivity. All of these depend only on the radial coordinate r with the assumption of initial spherical symmetric distributions of these physical quantities. The details of the derivation of the Eq. (22) can be found in the Appendix A. Then the variation of the temperature in (22) can

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

Fig. 4. The temperature distribution of the outer surface via (a) the symmetry preserving scheme and (b) the conventional scheme.

5

6

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

be rewritten as

T = f ( r )

(23)

which implies that the computed temperature will be spherically symmetric if the initial conditions (ICs) are symmetrical. Therefore, the developed diffusion scheme can preserve spherical symmetry. From the Eq. (21), we find that the construction of the Jacobian matrix requires that each vertex of the cell is shared by three faces of the cell. This condition can be satisfied in the hexahedral, tetrahedral, and prismatic grids, except pyramidal grid. In order to implement the symmetry preserving scheme, it is necessary to split pyramidal grid into two tetrahedral grids. The details can be found in the Appendix A. 3. Numerical discussion In this section, numerical tests in 3D Cartesian coordinate system are studied to prove the symmetry preserving property of the developed scheme. First, we test the accuracy and convergence of the scheme using a two-material spherical diffusion problem which has an analytical solution. Then, we perform several numerical tests of ICF target implosion include the hydrodynamically equivalent implosion, one-temperature approximation implosion of an ICF target. The performance of the developed scheme is shown in the following part.

Fig. 5. The schematic of an ICF target in the hydrodynamically equivalent implosion.

3.1. Accuracy test of the developed diffusion scheme First, we validate our new scheme using a two-material spherical diffusion problem which has an analytical solution. In order to make a comparison, the traditional support operator scheme without any symmetrical modification is also used to solve this problem. The problem can be described as follows:

ρC



 ∂T 1 ∂ ∂T = 2 r2 K + 1 + r 2 r = x2 + y2 + z 2 ∂t ∂r r ∂r

(24)

for r ∈ [0, 1], where



K=

K1 0 ≤ r ≤ 0.5 K2 0.5 < r ≤ 1.0

(25)

The problem domain consists of a two-region sphere. The inner region is defined by 0 < r < 0.5, and the outer region is defined by 0.5 < r < 1.0. The diffusion coefficient is discontinuous at r = 0.5. The temperature of the outer surface is specified as follows:

T |r=1 =

16 15

(26)

The exact steady solution to this two-material problem is

T =

⎧ 2 4 ⎪ ⎨ 16 + 43 + 11 − r − r 0 ≤ r ≤ 0.5 15

960K1

15

60K2

64K2

6K1

20K1

2 4 ⎪ ⎩ 16 + 13 − r − r 0.5 < r ≤ 1.0

6K2

(27)

20K2

Thus this problem is spherically symmetrical, and the solution is only dependent on the radial coordinate r. We solve this problem in 3D Cartesian coordinate system on spherical meshes with K1 = 1.0, K2 = 2.0, ρ = 1.0, C = 1.0. We perform several calculations with seven meshes of various sizes. The cell width along the radius direction of each mesh is 1/10, 1/20, 1/30, 1/40, 1/60, 1/80 and 1/100 respectively. The maximum error for each calculation is computed. This error is defined as follows:





Emax = max Tiexact − Tinumerical  i

(28)

We make a linear fit to the logarithm of the error as a function of the logarithm of the cell width. Fig. 3a shows the error as a function of the cell width using the symmetry preserving scheme. The

Fig. 6. The time history of the pressure at the boundary for the hydrodynamically equivalent implosion of the ICF target. Table 1 Comparison of the error and convergence rate of the symmetry preserving scheme and support operator scheme. Symmetry preserving scheme

Support operator scheme

h

Error

Rate

Error

Rate

1/10 1/20 1/30 1/40 1/60 1/80 1/100

1.22 × 10−2 3.51 × 10−3 1.62 × 10−3 9.25 × 10−4 4.17 × 10−4 2.36 × 10−4 1.52 × 10−4

1.80 1.91 1.94 1.96 1.98 1.98

1.27 × 10−2 3.57 × 10−3 1.64 × 10−3 9.33 × 10−4 4.20 × 10−4 2.37 × 10−4 1.52 × 10−4

1.83 1.93 1.95 1.97 1.98 1.99

slope of this linear function is 1.91. Perfect second-order convergence corresponds to a slope of 2.0. Thus the convergence rate of our scheme is near to second order. Fig. 3b shows the error and the convergence rate of the conventional scheme. Its convergence rate is 1.93, which is a little faster than that of our scheme. The errors and convergence rates for these calculations are given in the Table 1. Both schemes give comparable errors and convergence rates on all meshes. Next we study the symmetry properties of two schemes. Fig. 4 shows temperature distributions of the outer surface of the sphere computed by the traditional support operator scheme and the symmetry preserving scheme respectively. From the figure, the temperature of our scheme can preserve 1D spherical symmetry very well while the traditional scheme results in small departure from the symmetry. Although this departure has only the magnitude of 10−5 , it will be amplified obviously when the strong compression exists. This will be shown in next sections. From this test we find the convergence rate of our diffu-

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

7

Fig. 7. The pressure distribution and grid via (a) the conventional and (b) the symmetry preserving scheme for the ICF target in the hydrodynamically equivalent implosion at t = 5.0 ns.

sion scheme is very close to that of the traditional support operator scheme, although some symmetry preserving modifications are made to it. 3.2. Hydrodynamically equivalent implosion in an ICF target In this section, a hydrodynamically equivalent implosion of an ICF target is studied. The ICF target consists of deuterium-tritium (DT) ice and gas layers, as shown in Fig. 5. The outer and inner radii of DT ice layer are 950 μm and 870 μm, respectively. The initial density of DT ice and gas are 0.25 g/cm3 and 0.3 mg/cm3 , respectively. Ideal gas equation of state (EOS) is used in the simulation. The hydrodynamically equivalent pressure acts on the outer

Table 2 Pressure vs time for the boundary condition.

t(0.1 ns) p(Mbar)

t0

t1

t2

t3

t4

26.130 1

31.485 4

33.462 16

50 64

60 100

interface of DT ice layer. The time history of pressure at boundary is shown in Fig. 6. The corresponding values of pressure vs time in Fig. 6 can be found in Table 2. The BCs with uniform pressure and small perturbation (1 + 0.01cos2θ ) are simulated respectively.

8

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

Fig. 8. The shape of (a) the DT ice and (b) the DT gas of the ICF target. (c)The pressure distribution and grid via the symmetry preserving scheme for the ICF target in the hydrodynamically equivalent implosion at t = 5.9 ns.

Hydrodynamically equivalent implosion in ICF target with uniform pressure BC is simulated via the conventional scheme and the symmetry preserving scheme. Fig. 7(a) shows the pressure distribution and the grid of the cross section at 5.0 ns by using of the conventional scheme. Numerical errors due to discretization are amplified in the numerical integration. The result departs from 1D spherical symmetry. The pressure force calculated via the conventional scheme does not point in the radial direction. Fig. 7(b) shows the corresponding result by using of the symmetry preserving scheme. The pressure distribution and the grid of the cross section at 5.0 ns preserve 1D spherical symmetry very well. It justifies that the symmetry preserving scheme can remove numerical errors in 3D Cartesian coordinate system effectively. Further Fig. 8 shows the shape of DT ice and gas in the ICF target and corresponding pressure distribution and the grid at 5.9 ns. The implosion is close to the retarding stage at this time. At the stage of retarding nu-

merical errors due to small departure from 1D spherical symmetry will be amplified obviously via the conventional scheme. However, as shown in Fig. 8, the shape of DT ice and gas keep spherical perfectly and the pressure at the cross section still preserve 1D spherical symmetry. The symmetry preserving scheme keeps working well even at the state of strong compression. As mentioned in Section 2.3.2, the appropriate calculation of auxiliary points in the developed symmetry preserving scheme is very important. Inappropriate calculation of auxiliary points, for example the arithmetic averages of the surrounding points are used as the centers of face and cell for tetrahedral and prismatic grids, will result in small departure from 1D spherical symmetry. Numerical error due to small departure will be amplified during the simulation of the ICF capsule implosion. Fig. 9 shows the pressure distribution and the grid in the center region of ICF target with inappropriate and appropriate geometrical calculation for

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

9

Fig. 9. The pressure distribution and grid in the center region with (a) inappropriate and (b) appropriate geometrical calculation for auxiliary points via the symmetry preserving scheme for the ICF target in the hydrodynamically equivalent implosion at t = 5.9 ns.

auxiliary points via the symmetry preserving scheme at 5.9 ns. As shown in Fig. 9(a), the grid in the center region of ICF target with inappropriate calculation for auxiliary points seems to preserving 1D spherical symmetry. However, the pressure in the center region already departs from 1D spherical symmetry obviously. Corresponding results with appropriate calculation for auxiliary points are shown in Fig. 9(b). The pressure distribution in the center region preserves 1D spherical symmetry perfectly. The pressure of cells with the same I index has uniform magnitude. Fig. 10 gives the grid in the center region at 6.3 ns. The grid with inappropriate calculation of auxiliary points has very large distortion, as shown in Fig. 10(a). The grid with appropriate calculation of auxiliary points preserves 1D spherical symmetry very well. It justifies the importance of appropriate calculation of auxiliary points in the simulation of the 1D spherically symmetrical problem. Then we compare our hydrodynamic results to 1D simulation with traditional support operator scheme for 1D. Because this case is a spherically symmetrical problem, the 1D simulation can be used as a reference solution. As shown in Fig. 11, the time history of the position of the outer and inner interfaces of DT ice in 3D compares with that in 1D. From the figure, both interfaces of 3D simulation agree well with the reference result (1D), except that the inner interface is a little slower during a short period when it begins to accelerate.

Next, the hydrodynamically equivalent implosion in the ICF target with a small perturbation (1 + 0.01cos2θ ) pressure BC is simulated via the symmetry preserving scheme. Fig. 12 shows the pressure distribution at the cross section of the ICF target and the shape of DT gas at t = 5.9 ns. The pressure is higher at the region near the polar axis around the outer interface of DT gas. At this time DT gas looks like a ring. Numerical instability has been removed effectively by using of the symmetry preserving scheme. Since the initial small perturbation is independent on azimuthal ϕ , the boundary condition has cylindrical symmetry. Therefore, the distribution of the pressure should preserve cylindrical symmetry during the evolution of the capsule implosion process. As shown in Fig. 12, the pressure distribution is cylindrically symmetric. The pressure distribution ranges from 50 to 700Gpa, which is comparable with that of the previous case of uniform pressure BC. Therefore, the results in Fig. 12 are qualitatively correct in physics. 3.3. One-temperature approximation implosion in a spherical target In this section, we simulate the one-temperature approximation implosion of a spherical target which is composed of three different material zones. As shown in Fig. 13, the outer radius of the sphere is 1110 μm. The radii of material interfaces are 950 μm and 870 μm, respectively. The initial densities of three materials are 1.1 g/cm3 , 0.25 g/cm3 and 0.0 0 03 g/cm3 , respectively. In the simu-

10

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

Fig. 11. The position of the outer and inner interfaces of DT ice moves with time. Table 3 Radiation temperature vs time for the boundary condition. Time(0.1 ns)

Tr (MK)

Time(0.1 ns)

Tr (MK)

t0 = 0 t1 = 1.0 t2 = 71.47 t3 = 72.47 t4 = 95.89 t5 = 96.89

Tr0 = 0.0003 Tr1 = 1.0 Tr2 = 1.0 Tr3 = 1.40 Tr4 = 1.40 Tr5 = 1.96

t6 = 106.7 t7 = 117 t8 = 135.5 t9 = 146.9 t10 = 147.4 t11 = 150.0

Tr0 = 1.96 Tr1 = 2.87 Tr2 = 3.25 Tr3 = 3.25 Tr4 = 0.1 Tr5 = 0.1

at the boundary is shown in Fig. 14. Corresponding values of the temperature vs the time in Fig. 14 can be found in Table 3. This section is divided into three parts. Section 3.3.1 simulates only the pure diffusion problem without any hydrodynamic motion. The uniform radiation temperature on the boundary is exerted to maintain the spherically symmetry of the problem, which is devised to test the effectiveness of the symmetry preserving scheme for the diffusion equation. Section 3.3.2 simulates the onetemperature implosion problem where both the hydrodynamics and the diffusion are present. Similarly, the uniform boundary temperature condition is given to test the ability to preserve spherical symmetry of the scheme for the equations coupling diffusion to hydrodynamics. In Section 3.3.3, we use the small perturbation (1 + 0.005cos2θ ) temperature boundary condition instead of the uniform temperature. In this case, we test the accuracy and robustness of the scheme when the symmetry is absent.

Fig. 10. The grid in the center region with (a) inappropriate and (b) appropriate geometrical calculation for auxiliary points via the symmetry preserving scheme for the ICF target in the hydrodynamically equivalent implosion at t = 6.3 ns.

lation, the ideal gas equation of state with γ = 5/3 is used.

P = ( γ − 1 )ρ e

(29)

The heat conductivity is taken in the form κ = cρ m Tn where c, m and n are constants. c = 100, m = −2, n = 3 are used for all the materials. Irradiation on the outer surface of the target is applied via setting the radiation temperature of the outer surface according to the time history curve. The time history of the radiation temperature

3.3.1. Pure diffusion in the spherical target In this section, pure diffusion in the target is simulated. There is not any motion and the grids do not move. Fig. 15 shows the radiation temperature distribution of the outer surface and cross section of the target via the conventional diffusion scheme and the symmetry preserving diffusion scheme at t = 15.0 ns. Results via the conventional diffusion scheme are shown in Fig. 15(a). The temperature distribution departs from 1D spherical symmetry obviously. However, the results via the symmetry preserving diffusion scheme, as shown in Fig. 15(b), preserve 1D spherical symmetry. The symmetry preserving diffusion scheme removes numerical errors effectively. From the Fig. 15(a), the errors for conventional scheme are more obvious along the vertical axis than those along the horizontal axis. It results from the difference of the topologies of grids along the vertical axis and the horizontal axis. In the

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

11

Fig. 12. The pressure distribution of the ICF target in the hydrodynamically equivalent implosion with a small perturbation (1 + 0.01cos2θ ) at t = 5.9ns: (a) cross section, (b) DT gas.

polar region around the vertical axis, there are prismatic grids. In other regions, there are hexahedral grids. Therefore, the topologies of grids along the vertical axis are different from those in other regions. Different topologies may lead to the asymmetry of the geometrical objects, such as face centers, cell centers and normal vectors, if they are not calculated appropriately. Then the departure of geometrical symmetry will lead to differences of the discrete flux around the vertical axis and other regions. It is one of the key issues to calculate the auxiliary points appropriately in the symmetry preserving scheme. The symmetry error is introduced to evaluate the spherical symmetry quantitatively. The average temperature of a spherical layer R is defined as follows:

T¯ (R ) =

N (R ) 1  Ti (R ) N (R )

(30)

i=1

Fig. 13. The schematic of an ICF target for the one-temperature approximation.

4 (t8,Tr8)

Tr(MK)

3

(t9,Tr9)

(t7,Tr7)

where Ti (R) denotes the temperature of the grid belonging to the spherical layer R, and N(R) is the total number of the grids belonging to the layer. The average temperature is treated as the base value for the symmetry evaluation of the spherical layer R. The relative symmetry error of the layer R can be then estimated using both the maximum norm



i

2 1 0

(t1,Tr1)

(t0,Tr0)

0

(t3,Tr3)

Er L2 ( R ) =

(t6,Tr6)





Ti (R ) − T¯ (R ) T¯ (R )

2 .

(32)

For the whole spherical computational mesh, the overall symmetry error can be estimated using the maximum value of all the spherical layers in the mesh.

(t4,Tr4)

100 t(0.1ns)

N (R ) i=1

(t2,Tr2)

50

(31)

and the mean-square norm



(t5,Tr5)





Er max (R ) = max  Ti (R ) − T¯ (R ) /T¯ (R )

(t10,Tr10)

150

Fig. 14. The time history of the radiation temperature at the boundary for the target implosion.

Er max = max (Er max (R ) )

(33)

ErL2 = max (ErL2 (R ) )

(34)

R

R

With the above definitions, quantitative overall symmetry errors with the conventional and symmetry preserving scheme are calculated and given in Table 4. As shown in Table 4, symmetry errors with the symmetry preserving diffusion scheme are 5–6 order of

12

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

Temperature (MK)

t = 15.0 ns

t = 15.0ns R = 483 m

Temperature (MK)

0.800

1.000

0.794 0.789 0.783

0.922 0.844 0.767

0.778 0.772

0.689 0.611

0.767 0.761 0.756

0.533 0.456 0.378

0.750

0 .300

(a) Temperature (MK)

t = 15.0 ns

t = 15.0ns R = 483 m

Temperature (MK)

0.800

1.000

0.794 0.789 0.783

0.922 0.844 0.767

0.778 0.772

0.689 0.611

0.767 0.761 0.756

0.533 0.456 0.378

0.750

0.300

(b) Fig. 15. The radiation temperature distribution via (a) the conventional and (b) the symmetry preserving scheme at t = 15.0 ns during the diffusion in the target.

Table 4 The symmetry errors for diffusion in the ICF target.

Conventional scheme Symmetry preserving scheme

Ermax

ErL2

1.96 × 10−1 3.48 × 10−6

1.51 × 10−2 2.78 × 10−7

magnitude lower than those with the conventional scheme. Radial distributions of symmetry errors of both schemes are also shown in Fig. 16. Symmetry errors in the whole target are obviously reduced by using of the symmetry preserving diffusion scheme. 3.3.2. One-temperature approximation implosion in the target with uniform BC In this section, One-temperature approximation Eqs. (1)–(3) with uniform temperature boundary conditions are solved to simulate the implosion in the spherical target via the conventional and the symmetry preserving scheme. Fig. 17 shows the radiation temperature distribution of the outer surface and cross section of the target. Fig. 17(a) shows that results with the conventional scheme depart from 1D spherical symmetry obviously. Fig. 17(b) shows that results with the symmetry preserving scheme

Fig. 16. The radial distributions of symmetry errors via the conventional and the symmetry preserving diffusion scheme at t = 10.93 ns.

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

13

Fig. 17. The radiation temperature distribution via (a) the conventional and (b) the symmetry preserving scheme at t = 10.7 ns during the one-temperature approximation implosion in the target.

are spherically symmetrical. Quantitative comparisons of symmetry errors are shown in Fig. 18. Comparing with results via the conventional scheme, it can be found that the maximum error is significantly reduced by 9 or 10 orders of magnitude via the symmetry preserving scheme.

3.3.3. One-temperature approximation implosion in the target with small perturbation Instead of a uniform temperature in Section 3.3.2, a perturbed temperature which varies as T = Tr (1 + 0.005cos2θ ) in θ direction is exerted on the outer surface of the target. The temporal variation of the temperature Tr is the same as those of the Section 3.3.1 and 3.3.2. In order to reduce numerical errors, the symmetry preserving scheme is used in the simulation. The grid with the temperature contour on the outer surface and cross section of the inner gas in the target at the time t = 13.0 ns are shown in Fig. 19. The time history of the shape and the temperature distribution of the inner gas region from 12.9 ns to 15 ns are given in Fig. 20. The inner gas is compressed and reach the retarding stage at around t = 13.3 ns. The shape of the inner gas is close to that of hydrodynamically equivalent implosion at the stage of retarding. Then the inner gas expands and stretches along the polar direction. The result is qualitatively correct in physics. Numerical results show that the sym-

Fig. 18. The radial distributions of symmetry errors via the conventional and the symmetry preserving diffusion scheme in the one-temperature approximation implosion in the target at t = 10.41 ns.

14

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

tained by the small modification of the flux-related Jacobian matrix. It is proved in theory that the developed scheme guarantees 1D spherical symmetry. It is emphasized that the appreciate calculation of auxiliary points are important for the symmetry preserving scheme. The developed symmetry preserving scheme is applied to the 3D simulation of an ICF target implosion. Hydrodynamics equivalent and one-temperature approximation radiation implosions in the ICF target are simulated. Numerical results show that the developed symmetry preserving scheme can guarantee 1D spherical symmetry and remove the numerical instability effectively in the 3D simulation of the ICF target implosion. The developed symmetry preserving scheme will be applied further to 3D simulations of an ICF target with different approximations and boundary conditions. The idea of the symmetry preserving scheme can be extended to ALE method or other numerical methods when the symmetrical property should be emphasized in numerical simulations. Declaration of Competing Interest Fig. 19. The radiation temperature distribution of DT gas in the one-temperature approximation implosion with a small perturbation (1 + 0.005cos2θ ) at t = 13.0 ns.

None. Acknowledgements

metry preserving scheme helps to remove the numerical error in the 3D simulation of a target implosion effectively. 4. Conclusion A newly developed symmetry preserving scheme for Lagrangian radiation hydrodynamics is present. The developed scheme is ob-

The authors acknowledge the support from National Natural Science Foundation of China (Grant Nos. 11205016, 11372050, 11371066, and 11371068) and projects by CAEP (Grant Nos. 2014B0202031, and 2013A0101004). The authors also thank Dr. Yingyan Wu for her constructive work for the developed code on the symmetry preserving scheme for Lagrangian hydrodynamics.

Fig. 20. The time history of the shape and the temperature distribution in the DT gas region of the target with a small perturbation (1 + 0.005cos2θ ) from t = 12.9 ns to 15 ns. (a) t = 12.9 ns (b) t = 13.1 ns (c) t = 13.3 ns, (d) t = 13.5 ns (e) t = 14.0 ns (f) t = 14.5 ns (g) t = 15.0n.

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

15

Appendix A. The proof of the symmetry preserving property of the developed diffusion scheme on the spherical mesh In this appendix, it is demonstrated that new developed scheme preserves the spherical symmetry on a spherical mesh. First, it is assumed that the temperature of cell center Tc and face center Tf at the (n)th time level are symmetrical and independent on the polar angle θ and the azimuthal angle ϕ (as shown in Fig. 1), and Tc and Tf can be expressed as the function of the radial coordinate r

Tc n = Tc n (r )T f n = T f n (r ).

(A.1)

Next it will be proved that the temperature at the (n + 1)th time level Tn +1 is also independent on θ and ϕ . According to the developed scheme of the Section 2.2 in the main text, the change of the temperature can be obtained by

ρ cv

s Tc n+1 − Tc n 1 = f i Ai t V

(A.2)

i=1

where fi is the normal component of the flux across the ith face of the cell, Ai is the area of the ith face of the cell, ρ is the density, and cv is the heat capacity. It is obvious that ρ cv /t only depends on the radial coordinate r. The term on the right hand side of Eq. (A.2) is related to the geometries of the mesh as well as the discrete flux. In current simulation, the mesh is shown in Fig. 21, where exists four kinds of grids with different topologies. Fig. 21(a) shows the grids in the central layer where the tetrahedral grids locate in polar zones and pyramidal grids locate in other zones. Fig. 21(b) shows the outer layer where prismatic grids locate in polar zones and hexahedral grids locate in other zones. In the following, we will present the discrete forms of Eq. (A.2) on these four kinds of grids. Case 1: tetrahedral grid A tetrahedral cell OABC is shown in Fig. 22. Vertex O in the sphere center and vertex A, B, C locate on a sphere surface. According to Eq. (19), the flux-related Jacobian matrix for vertex A is given by



JA T JA =

1  OAB  OAC · n n  ABC · n  OAB n

 OAC  OAB · n n 1  ABC · n  OAC n



 ABC  OAB · n n  ABC .  OAC · n n 1

(A.3)

The normal vector of the face ABC related to vertex A is modified as follows



 AABC ∗ n

=

 OAB × n  OAC n  ·n |n OAB × n OAC | ABC



 OAB × n  OAC n . |n OAB × n OAC |

(A.4)

Using Eq. (A.4), the Jacobian matrix for vertex A is rewritten as



JA T JA

∗



=

1  OAB  OAC · n n 0

 OAB · n  OAC n 1 0



0 0 . ∗ ∗  AABC  AABC · n n

(A.5)

Fig. 22. A tetrahedral cell.

In a similar way, the Jacobian matrix for vertex B and C can be given by





T

∗

T

∗

JB JB

JC JC



=

 =

where



B ∗  ABC n

=

∗

=

  CABC n

1  OBC  OAB · n n 0

 OBC · n  OAB n 1 0

0 0 ∗ ∗  BABC  BABC · n n

1  OAC  OBC · n n 0

 OAC · n  OBC n 1 0

0 0 ∗ ∗  CABC  CABC · n n

 OAB  OBC × n n  ·n |n OBC × n OAB | ABC  OBC  OAC × n n  ·n |n OAC × n OBC | ABC





 OAB  OBC × n n

 OBC  OAC × n n

(A.7)

(A.8)

|n OBC × n OAB | |n OAC × n OBC |

(A.6)

.

(A.9)

According to Eq. (13), we define the face-normal component of the flux as T F˜ = ( fOAB , fOBC , fOAC , fABC ) .

(A.10)

By substituting Eqs. (A.5)–(A.7) into Eq. (13) and using the corner cell volume Vn = V/4, the flux matrix can be obtained as follows



F˜ =



−1

PnT Jn−1 K −1 Jn−T PnVn

n

⎛

G11 4K ⎜G21 = ⎝G V 31 0 where H2 A

Fig. 21. The grids of a spherical mesh: (a) the central layer (b) the outer layer.





G12 G22 G32 0

G44 = AABC /(

G13 G23 G33 0

AT˜

⎞⎛

0 T 0 ⎟⎜T 0 ⎠⎝T G44 T

1 ∗ · ∗ A n nA ABC ABC

+

⎞

− TOAB − TOBC ⎟ − TOAC ⎠ − TABC

1 ∗ · B ∗ B n nABC ABC

(A.11)

+

1 ∗ · C n nCABC ∗ ABC

)=

ABC , and TOAB , TOBC , TOAC , TABC are the temperature of the face 3R2 OAB, OBC, OAC and ABC, respectively, and K is the conductivity. Because the radial coordinate r of the center of OAB, OBC, OAC are equal to that of cell center, which is guaranteed by appropriate

16

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

Similarly, the normal vector of the face ABCD related to vertex A is modified as



 ADAB n

∗

=

 BAE  DAE × n n  ·n |n DAE × n BAE | DAB



 BAE  DAE × n n

|n DAE × n BAE |

.

(A.17)

The other normal vectors such as B ∗, n D ∗, n E ∗, n H ∗  ABC  CBCD ∗ , n  CDA  HEF  FEF G ∗ , n  FGGH ∗ , n  GHE n can also be computed in the same way. With these normal vectors, we can obtain the Jacobian matrixes which have similar forms to Eqs. (A.5)–(A.7). The face-normal components of fluxes has the form T F˜ = ( fAEHD , fAEF B , fBF GC , fCGHD , fABCD , fEF GH ) ,

(A.18)

which can be obtained using Eqs. (13) and (14). Thus



−1



F˜ =

PnT Jn−1 K −1 Jn−T PnVn

n

⎛

AT˜

⎞⎛

G11 G12 G13 G14 00 T ⎜G21 G22 G23 G24 00⎟⎜T ⎜G31 G32 G33 G34 00⎟⎜T ⎟⎜ = K⎜ ⎜G41 G42 G43 G44 00⎟⎜T ⎝ ⎠⎝ 0 0 0 0G55 0 T 0 0 0 0 0G66 T where G55 = AABCD /( AABCD HABCD

1 ∗ H ∗ H n · nGHE GHE

calculations in Section 2.3.2. The temperature of these face centers are equal to the temperature of cell center

(A.12)

Vn,ABCD =

)=

AEF GH HEF GH 2 4Vn,EF GH REF GH 2

V 8

(A.19)

1 1 1 ∗ B ∗ + C ∗ C ∗ + D ∗ D ∗) = B   n · nABC n · nCDA n · nBCD BCD ABC CDA 1 1 AEF GH /( E ∗ E ∗ + F ∗ F ∗ + G ∗1 G ∗ +    n · nHEF n · nEF G n · nF GH HEF EF G F GH

1 ∗ A ∗ A n · nDAB DAB

G66 =

4Vn,ABCD RABCD 2

Fig. 23. A hexahedral cell.

T = TOAB = TOBC = TOAC .

2

⎞

− TAEHD − TAEF B ⎟ − TBF GC ⎟ ⎟ − TCGHD ⎟ ⎠ − TABCD − TEF GH

Vn,EF GH =

+

and corner cell volume is computed by

V . 8

(A.20)

By substituting Eq. (A.12) into Eqs. (A.10) and (A.11), we obtain the all face-normal components of fluxes of the cell as follows:

Considering the temperature of the face AEHD, AEFB, BFGC, CGHD are equal to the temperature of cell center, the remaining non-zero fluxes in Eq. (A.12) are given by

fOAB = fOBC = fOAC = 0

fABCD =

fABC =

4H 2 K AABC (T − TABC ) 3R2 V

(A.13)

(A.14)

where H is the height from the sphere origin O to the face ABC. fABC represents the normal component of the flux in radial direction across the face ABC from the neighbor cell with the same θ and ϕ , and fOAB , fOBC , fOAC represent the normal components of the fluxes across the face OAB, OBC, OAC from the neighbor cell with the same radius. Thus Eq. (A.13) shows that there is no flux between cells with the same radius coordinate r. By using of Eqs. (A.13) and (A.14), Eq. (A.2) becomes

T 4H 2 K AABC 2 ρ cv = (T − TABC ). t 3R2 V 2

(A.15)

Using V = AABC H/3, we obtain the variation of the temperature

T =

12K

ρ cv R2

(T − TABC )t.

(A.16)

According to Eq. (A.1), the cell center temperature T and face center temperature TABC are only dependent on the radial coordinate r. The density ρ , the heat capacity cv , and the conductivity K are also dependent on the radial coordinate r. Therefore, the variation of the temperature at the (n + 1)th time level Tn +1 can preserve the spherical symmetry. Case 2: Hexahedral grid A hexahedral cell is illustrated in Fig. 23, where the vertex A, B, C, D locate on a sphere surface with the same radial coordinate r and the vertex E, F, G, H locate on another sphere surface.

=

2K AABCD HABCD 2 V RABCD 2

(T − TABCD ) fEF GH

2K AEF GH HEF GH 2 V REF GH 2

(T − TEF GH ).

(A.21)

Considering the similarity of the face ABCD and EFGH, with these fluxes the variation of the temperature can be deduced as follows

T =



18R4ABCD K

ρ cv R3ABCD − R3EF GH +



2 (T − TABCD )t

18R4EF GH K

ρ cv R3ABCD − R3EF GH

2 (T − TEF GH )t

(A.22)

where the right hand side is only dependent on the radial coordinate r. Therefore, the sphere symmetry can be preserved in a hexahedral cell. Case 3: prismatic grid For a prismatic cell illustrated in Fig. 24, the normal vectors of the face ABC and DEF are also modified to be along the radial direction. Similar to the case of hexahedral cell, the variation of the temperature at the (n + 1)th time level can be deduced and given by

T =



18R4ABC K

ρ cv R3ABC − R3DEF +



2 (T − TABC )t

18R4DEF K

ρ cv R3ABC − R3DEF

2 (T − TDEF )t

(A.23)

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

17

Fig. 24. A prismatic cell.

Fig. 26. Two tetrahedral cells by splitting a pyramid cell.

According to Eqs. (A.5)–(A.7), the construction of the Jacobian matrix requires that each vertex is shared by three faces of the cell. In the case of the tetrahedral cell, the vertex A is shared by the face OAB, OAC and ABC. This condition is satisfied in tetrahedral, hexahedral and prismatic cells. However, this condition is not satisfied in pyramidal. As shown in Fig. 25, the vertex O is shared by four faces. Here we split the pyramid cell into two tetrahedral cells illustrated in Fig. 26. Then the normal vector modifications can be made and the flux-related Jacobian matrixes can be obtained in each tetrahedral cell in the same way of the Case 2. For the front tetrahedral cell OABC of Fig. 26, the normal vector of the face ADC related to the vertex A is given by



 AABC ∗ n

=

 OAC  OAB × n n  ·n |n OAB × n OAC | ABC



 OAC  OAB × n n

|n OAB × n OAC |

.

(A.24)

For the back tetrahedral cell OADC, the normal vector of the face ADC related to the vertex A is



 AADC n

Fig. 25. A pyramidal cell.

which is identical to that of the hexahedral cell. The spherical symmetry can be preserved in a prismatic cell. Case 4: pyramidal grid In a pyramidal cell illustrated in Fig. 25, vertex O is the sphere center and vertex A, B, C, D locate on a sphere surface.

∗

=

 OAC  OAD × n n  ·n |n OAD × n OAC | ADC



 OAC  OAD × n n

|n OAD × n OAC |

.

(A.25)

After eliminating the additional unknowns of the internal face OAC which is introduced by the splitting process, the flux matrix for the pyramidal cells can be obtained as follows

⎛G G G G 0 ⎞⎛T − T 11 12 13 14 OAB ⎜G21 G22 G23 G24 0⎟⎜T − TOBC F˜ = K ⎜G31 G32 G33 G34 0⎟⎜T − TOCD ⎝ ⎠⎝ G41 G42 G43 G44 0 0 0 0 0G55

T − TCDA T − TABCD

⎞

⎟ ⎟. ⎠

(A.26)

18

S. Guo, M. Zhang and H. Zhou et al. / Computers and Fluids 195 (2019) 104317

The only non-zero flux is given as follows

fABCD =

4H 2 K AABCD (T − TABCD ). 3R2 V

(A.27)

Substituting Eq. (A.27) into the conservation Eq. (A.2), the variation of the temperature with time can be obtained as follows

T =

12K

ρ cv R2

(T − TABCD )t

(A.28)

which is similar to the Eq. (A.16) of the tetrahedral grid. Therefore, we prove that the developed scheme can preserve the spherical symmetry in tetrahedral, hexahedral, prismatic, and pyramidal grids. Reference [1] Lipnikov K, Manzini G, Shashkov M. Mimetic finite difference method. J Comput Phys 2014;257:1163–227. [2] Shashkov M, Steinberg S. Support-operator finite-difference algorithms for general elliptic problems. J Comput Phys 1995;118:131–51. [3] Shashkov M, Steinberg S. Solving diffusion equations with rough coefficients in rough grids. J Comput Phys 1996;129:383–405. [4] Morel J, Roberts R, Shashkov M. A local support-operator diffusion discretization scheme for quadrilateral r-z meshes. J Comput Phys 1998;144:17–51. [5] Kuznetsov Y, Lipnikov K, Shashkov M. The mimetic finite difference method on polygonal meshes for diffusion-type problems. Comput Geosci 2004;8:301–24. [6] Lipnikov K, Shashkov M, Svyatskiy D. The mimetic finite difference discretization of diffusion problem on unstructured polyhedral meshes. J Comput Phys 2006;211:473–91. [7] Lipnikov K, Morel J, Shashkov M. Mimetic finite difference methods for diffusion equations on non-orthogonal non-conformal meshes. J Comput Phys 2004;199:589–97. [8] Berndt M, Lipnikov K, Shashkov M, Wheeler MF, Yotov I. A mortar mimetic finite difference method on non-matching grids. Numer Math 2005;102(2):203–30.

[9] Caramana EJ, Burton DE, Shashkov MJ, Whalen PP. The construction of compatible hydrodynamics algorithms utilizing conservation of total energy. J Comput Phys 1998;146(1):227–62. [10] Campbell JC, Shashkov MJ. A compatible Lagrangian hydrodynamics algorithm for unstructured grids. Selcuk J Appl Math 2003;4(2):53–70. [11] Jia ZP, Liu J, Zhang SD. An effective integration of methods for second-order three-dimensional multi-material ALE method on unstructured hexahedral meshes using MOF interface reconstruction. J Comput Phys 2013;236:516–62. [12] Caramana EJ, Shashkov MJ. Elimination of artificial grid distortion and hourglass-type motions by means of Lagrangian subzonal masses and pressures. J Comput Phys 1998;142(2):521–61. [13] Campbell JC, Shashkov MJ. A tensor artificial viscosity using a mimetic finite difference algorithm. J Comput Phys 2001;172:739–65. [14] Caramana EJ, Shashkov MJ, Whalen PP. Formulations of artificial viscosity for multi-dimensional shock wave computations. J Comput Phys 1998;144(1):70–97. [15] Caramana EJ, Whalen PP. Numerical preservation of symmetry properties of continuum problems. J Comput Phys 1998;141(1):174–98. [16] Caramana EJ, Rousculp CL, Burton DE. A compatiable, energy and symmetry preserving Lagrangian hydrodynamics algorithm in three-dimensional Cartesian geometry. J Comput Phys 20 0 0;157:89–119. [17] Margolin L, Shashkov M. Using a curvilinear grid to contruct symmetry-preserving disretizations for Lagrangian gas dynamics. J Comput Phys 1999;149:389–417. [18] Cheng J, Shu CW. A cell-centered Lagrangian scheme with the preservation of symmetry and conservation properties for compressible fluid flows in two-dimensional cylindrical geometry. J Comput Phys 2010;229:7191–206. [19] Cheng J, Shu C. A Lagrangian scheme with the preservation of symmetry and conservation in cylindrical geometry: preliminary study. Proc Comput Sci 2010;1(1):1903–11. [20] Velechovsky J, Kucharik M, Liska R, Shashkov M. Symmetry-preserving momentum remap for ALE hydrodynamics. J Phys 2013;454(1):012003. [21] Velechovsky J, Kucharik M, Liska R, Shashkov M, Vachal P. Symmetry- and essentially-bound-preserving flux-corrected remapping of momentum in staggered ALE hydrodynamics. J Comput Phys 2013;255:590–611. [22] Kenamond M, Bement M, Shashkov M. Compatible, total energy conserving and symmetry preserving arbitrary Lagrangian-Eulerian hydrodynamics in 2D rz-Cylindrical coordinates. J Comput Phys 2014;268:154–85. [23] Hypre high performance preconditioners. Center for applied scientific computing. USA: LLNL; 2006.