- Email: [email protected]
A numerical study of 2D detonation waves with adaptive finite volume methods on unstructured grids
Accepted Manuscript A numerical study of 2D detonation waves with adaptive finite volume methods on unstructured grids
Guanghui Hu
PII: DOI: Reference:
S0021-9991(16)30633-7 http://dx.doi.org/10.1016/j.jcp.2016.11.041 YJCPH 6991
To appear in:
Journal of Computational Physics
Received date: Revised date: Accepted date:
4 December 2015 6 November 2016 28 November 2016
Please cite this article in press as: G. Hu, A numerical study of 2D detonation waves with adaptive finite volume methods on unstructured grids, J. Comput. Phys. (2016), http://dx.doi.org/10.1016/j.jcp.2016.11.041
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A numerical study of 2D detonation waves with adaptive finite volume methods on unstructured grids Guanghui Hua,b a Department b UM
of Mathematics, University of Macau, Macao S.A.R., China Zhuhai Research Institute, Zhuhai, Guangdong Province, China
Abstract In this paper, a framework of adaptive finite volume solutions for the reactive Euler equations on unstructured grids is proposed. The main ingredients of the algorithm include a second order total variation diminishing Runge-Kutta method for temporal discretization, and the finite volume method with piecewise linear solution reconstruction of the conservative variables for the spatial discretization in which the least square method is employed for the reconstruction, and weighted essentially nonoscillatory strategy is used to restrain the potential numerical oscillation. To resolve the high demanding on the computational resources due to the stiffness of the system caused by the reaction term and the shock structure in the solutions, the h-adaptive method is introduced. OpenMP parallelization of the algorithm is also adopted to further improve the efficiency of the implementation. Several one and two dimensional benchmark tests on the ZND model are studied in detail, and numerical results successfully show the effectiveness of the proposed method. Keywords: reactive Euler equations; finite volume methods; unstructured grids; mesh adaptive methods;
1. Introduction The development of the numerical methods for the chemical reactive flows, which is an important application for multimaterial flows, has attracted considerable attention in the last few years[12, 25, 29]. The reactive Euler equations are one of the most important models for describing the physical phenomenon such as detonation. To date, there have been a lot of numerical methods for solving reactive Euler equations. For example, the finite difference methods[10, 14, 15], the gas-kinetic schemes[29], the finite volume methods[18, 36], the discontinuous Galerkin schemes[39], etc. In all these works, excellent numerical results can be observed. However, the structured mesh grids are employed in all these works, which restrains their application in the practical problems such as the numerical simulations for the rotating detonation engine. In the practical problems, the computational domain always has complicated geometry, hence the unstructured mesh grids become a popular choice to handle the complicated domain. For some early numerical methods for reactive Euler equations on unstructured meshes, people may refer to [31, 37] for the finite element methods, etc.
Preprint submitted to Journal of Computational Physics
November 30, 2016
Besides the ability on handling complicated geometry of the domain, the other feature the numerical methods desired for practical applications is the efficiency. In the reactive Euler equations, the existence of the shock structure in the computational domain makes that the mesh size around the shock region should be small enough to resolve the shock well, which means that the size of the time step in the numerical methods should also be small enough to guarantee the simulation stability. Moreover, the stiffness of the reactive Euler system makes the choice of the size of the time step more severe. On the other hand, since the computational domain is bounded, it is necessary to remove or reduce the unphysical wave reflected by the outflow boundary, from the reliability point of view. To do so, several techniques have been developed. For example, in [30], an extrapolation technique is employed, where the conserved quantities of the system are relaxed to the condition at infinity with a given formula. In this case, the subsonic outflow boundary condition can be imposed trivially. In [15], the perfectly matched layer (PML) method is introduced in the simulation, and excellent numerical results are obtained. A more straightforward way is to use a sufficiently large domain. The artificial effect introduced by the wave absorbing strategies mentioned above can be effectively avoided with this strategy. However, the requirement on the size of the domain, together with the requirement on the sizes of the mesh grids and time step, will result in the highly demanding on the computational resource. This issue motivates the introduction of the adaptive mesh methods. There are three different kinds of mesh adaptive methods in the market. The h-adaptive methods refine and/or coarsen the mesh grids locally, and have been widely used in the computational fluid dynamics[9, 20, 28]. The r-adaptive methods redistribute the mesh grids in the domain to better resolve the singularity in the solutions, people may refer to [5, 17, 27, 35] for the detail of the methods and their applications in different areas. The p-adaptive methods locally upgrade the order of the approximate polynomial in the element. In the practice, p-adaptive methods always work together with the h-adaptive methods, people may refer to [33] for their applications in CFD. For the simulations of the detonation wave with the mesh adaptive methods, there have been some pioneer works[2, 9, 31]. In our previous work[7, 4, 6, 20], a general framework for h-adaptive mesh methods has been developed and applied for the computational fluid dynamics and density functional theory in 2D/3D cases. In this paper, such h-adaptive mesh framework will be introduced to solve the reactive Euler equations. To further improve the implementation efficiency, the OpenMP parallelization of the algorithm will also be discussed. The numerical results shown in the last section will demonstrate the effectiveness of the proposed method. In this paper, we develop a general framework of using adaptive finite volume methods for solving the reactive Euler equations on unstructured meshes. The numerical discretization of the Euler equations includes a second order total variation diminishing Runge-Kutta method for the temporal discretization, and finite volume method with piecewise linear solution reconstruction for the spatial discretization. To effectively remove the potential numerical oscillation introduced by the solution reconstruction, we follow [21] to use
2
WENO reconstruction. The reconstruction patch for the element which has an edge located on the boundary of the domain is designed specially to avoid the poor gradient reconstruction. To reduce the unphysical effect generated by the reflection of the wave from the outflow boundary, we simply use a sufficiently large domain [29]. The h-adaptive methods are employed in this paper to enhance the implementation efficiency. To further improve the efficiency, the algorithm is also parallelized with OpenMP, and detail on the implementation is discussed. A variety of numerical simulations successfully show the numerical convergence and effectiveness of our method on solving the reactive Euler equations. The rest of the paper is organized as follows. In the next section, the reactive Euler equations are introduced in detail. In Section 3, the numerical discretization of reactive Euler equations are derived. Then the solution reconstruction, the mesh adaptive methods, and some issues on OpenMP parallelization of our algorithm will also be introduced in depth. In Section 4, the initial and boundary conditions used in the simulations are introduced first, then numerical results from a variety of tests are demonstrated to show the numerical convergence and effectiveness of our method. Finally, we give the conclusion of this work, and also some future work on this topic.
2. Model Statement In this paper, we focus on the two dimensional reactive Euler equations with two reaction species, i.e., reactant and production species. There is no intermediate species created during the reaction, and such single one-step and irreversible reaction can be described by k
→ P roduct, Reactant − where k denotes certain reaction rate. By ignoring the viscosity of the fluid, the convection and reaction phenomenon of the system during the detonation can be modeled by the following reactive Euler equations, ∂ U + ∇ · F (U ) = S(U ), ∂t
(1)
where U , F (U ), and S(U ) denote the conservative variable, flux, and source term, respectively. Their expressions are given as follows, ⎡ ⎤ ⎡ ρ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ρr ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ U = ⎢ ρu ⎥ , F (U ) = ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ρv ⎥ ⎢ ⎣ ⎦ ⎣ E
⎤ ρu,
ρv
ρr u,
ρr v
⎡
⎤ 0
⎥ ⎢ ⎥ ⎢ ⎥ ⎢ −ρr k ⎥ ⎢ ⎥ ⎢ ⎥ , and S(U ) = ⎢ 0 ρuv ⎥ ⎢ ⎥ ⎢ ⎢ 0 ρv 2 + P ⎥ ⎦ ⎣ (E + P )v 0
ρu2 + P, ρuv, (E + P )u,
3
⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎦
(2)
where ρ, ρr , u = (u, v)T , E, and P , stand for the density of the fluid, the density of the reactant, velocity, total energy, and pressure, respectively. The density of the fluid ρ, the density of the reactant ρr , and the density of the product ρp satisfy the relation ρ = ρr + ρp . In the source term S(U ), the non-zero term −ρr k represents the change of the reactant in the reaction. In this paper, the Arrhenium reaction rate is employed for the reaction rate. Hence, k has the following form k = Ke−Ea /T , where K is the pre-exponential factor, Ea is the activation energy, and T = P/ρ is the temperature (with gas constant normalized to 1 in this paper). To close the above system, the following equation of state is used E = P/(γ − 1) + 0.5ρ(u2 + v 2 ) + ρr q0 ,
(3)
where γ is the ratio of the specific heat, and q0 is the amount of heat released per unit mass by the reaction.
3. Numerical Discretization and Some Numerical Issues In this section, the numerical discretization for the reactive Euler equations is described in detail first. Then some numerical issues including the initial and boundary conditions, numerical stability, etc., as well as acceleration techniques such as h-adaptive methods and OpenMP parallelization are introduced. For the spatial discretization, the finite volume method is employed, not only because of its numerical conservative property which is very important for preserving the desired physical behavior of the numerical solution, but also because of its flexibility on handling the domain with complex geometry. To deliver a quality solution reconstruction, the linear reconstruction based on the least square method is employed, with a WENO strategy to restrain the possible nonphysical oscillation of the numerical solution around the shock. In the temporal discretization, the second order TVD Runge-Kutta scheme is adopted, which is also positive on removing the nonphysical oscillation. To describe the discretization specifically, let Ω denotes the computational domain of the problem. In tri this paper, we consider the unstructured meshes, and use T = {Ki }N i=1 to denote a triangulation of the
domain Ω, and Ntri is the total number of the triangle element in T . We suppose that there is no nontrivial overlap for every two triangle elements, i.e., the overlap part of two elements could be a vertex, an edge, or nothing, and that the union of the all elements is the domain Ω. The temporal domain we used in the NT −1 simulation is [0, Tend ], and its partition is denoted by {[tn , tn+1 ]}n=0 , where NT means the total number
of the segments in the interval [0, Tend ].
4
3.1. Runge-Kutta Finite Volume Discretization To derive the semi-discretized numerical scheme first, one takes the integral of governing equations (1) in a control volume Ki . Then the divergence theorem gives ∂ U dxdy + F (U ) · nds = S(U )dxdy, Ki ∂t ∂Ki Ki
(4)
where ∂Ki denotes the boundary of the i-th triangle element Ki in T , and n denotes the unit out normal vector w.r.t. ∂Ki . By introducing the cell average of conservative variable U and the source S(U ) in Ki , i.e., 1 ¯= 1 U U dxdy, and S¯ = S(U )dxdy, |Ki | Ki |Ki | Ki the equations (4) become 1 ∂ ¯ U+ ∂t |Ki |
¯ F (U ) · nds = S,
(5)
∂Ki
where |Ki | is the area of the element Ki . In the simulation, the flux F (U ) is replaced by the numerical flux F¯ , which is determined by a local Riemann problem with the initial conditions given by the solutions in the element Ki and its neighbor Kj . In our work, the Lax-Friedrichs flux is adopted for the numerical flux. With the numerical flux F¯ , the equations become 1 ∂ ¯ U = S¯ − ∂t |Ki |
F¯ · nds := RHS,
(6)
∂Ki
which is an ODE system now, and an ODE solver is needed for a full discretization. For the numerical simulation of gas dynamics with shock structures, non-oscillatory property of the numerical method is crucial not only for the numerical stability, but also for the reliability of the numerical solution. Consequently, we use the second order TVD Runge-Kutta method[16] in the temporal discretization, which results in the following full discretization of the reactive Euler equations (1) in the domain Ki × [tn , tn+1 ]. ⎧ ⎪ ¯ (1) = U ¯ n + δtRHS n , ⎪ U ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ U ¯ n + 0.5U ¯ (1) + 0.5δtRHS (1) , ¯ n+1 = 0.5U
(7)
¯ (1) . where RHS (1) means the evaluation of the RHS at U Besides the above Runge-Kutta finite volume method for reactive Euler equations, an alternative and competitive way is the splitting strategy[38] in which the two procedure, i.e., the advection of the reactant and the reaction, are supposed to happen in order. Since both the pure advection process and the stiff ordinary differential equations can be handled efficiently by classical methods, the splitting strategy may bring significant improvement on the efficiency. However, it has been shown in [38] that the splitting method 5
is only exact for some model problem. Further, it is nontrivial to design high order splitting methods. Hence, in this paper we use the above coupled process to describe the advection and reaction. To get a stable simulation, the selection of the size of the time step δt should follow the CourantFriedrichs-Lewy (CFL) condition. The idea behind the CFL condition is that for a given element, the incident wave is not allowed to move across the element in one single time step. To apply the CFL condition in the simulation, the wave speed and the cell size need to be determined. In terms of the characteristic analysis of the quasi-linear system derived from the reactive Euler equations (1), five eigenvalues obtained are λ1 = u − c, λ2 = λ3 = λ4 = u, and λ5 = u + c, where c denotes the sound speed. For the element size, the following approximation is used for it, si = 2 |Ki |/π,
(8)
where |Ki | means the area of the i-th element Ki . The idea of such formula is to use the diameter of a circle which has the same area to Ki as the size of the element, which is simple to calculate and reasonable for the mesh with regular element. With these parameters, the suggested time step δti for the element Ki is given by δti = σcf l si /(|u| + c),
(9)
where σcf l is the CFL number which is taken as 0.8 in this paper. Then the global time step δt is chosen as δt = min{δti }. i
(10)
3.2. Solution Reconstruction In the framework of the Godunov finite volume schemes, the available quantity from each governing cell is the cell average of each conservative variable. To deliver a high resolution method, solution variation in each element is needed, which can be obtained from the solution reconstruction. In this paper, we focus on the linear reconstruction of the conserved quantities in each element. Two components are needed for the solution reconstruction, i.e., the selection of the reconstruction patch for each element, and the reconstruction method. Below we describe the selection of the reconstruction patch first, then the reconstruction method follows. In the previous papers [20, 21], the reconstruction patch for the linear solution reconstruction has been discussed. Especially in [20], the influence of the different reconstruction patch for the h-adaptive methods is discussed and tested. In this work, we use two different reconstruction patches shown in Fig. 1 to reconstruct the solutions. For the first kind of the reconstruction patch for the current element K0 , the element which has a common edge with K0 will be selected as one element in the patch. With this principle, the patches P1 and P1B will be used as the reconstruction patches for the internal element and the element located on the domain boundary, 6
Figure 1: Left: the patch P1 = {Ki , i = 0, 1, 2, 3}, and the patch P2 = {Ki , i = 0, 1, . . . , 13} for an internal element K0 . Right: the patch P1B = {Ki , i = 0, 1, 2}, and the patch P2B = {Ki , i = 0, 1, . . . , 7} for an element K0 located on the boundary of the domain.
respectively. The drawback for P1B is obvious. First, there are only three elements (or even two in some extreme situation) in the patch, which would result in a determined (or under-determined in the extreme situation) system. Even worse, the distribution of all elements in the patch is almost along the domain boundary, which would cause the determined system ill-posed. The negative influence by using the patches P1 and P1B has been observed in [20] for steady Euler equations. Similar unsatisfied simulations are also obtained in this paper for reactive Euler equations. Consequently, in the simulation with uniform mesh, the patches P1 and P2B are employed. It is noted that in P2B , the patch is effectively enlarged, and the element along the normal direction w.r.t. the domain boundary is also introduced which effectively improves the condition number of the reconstruction system. Hence, numerical convergence of the simulation is observed successfully in the numerical experiments in the last section. When h-adaptive method is used, the mesh becomes nonuniform. Then the sizes of the neighbor elements for K0 could be very different. In [20], it has been shown that the use of P1 and P2B would bring negative effect on solving steady Euler equations, and to improve the robustness of the algorithm, the patches P2 and P2B are used for h-adaptive simulations. With this reconstruction patch, the smooth convergence of the steady state simulations is obtained. In this paper, it is found that the use of the patches P1 and P2B can handle the simulations with h-adaptive simulations well. Hence in the last section, we still use the patches P1 and P2B for the solution reconstruction. For simplicity of the description, we use Pi to denote the reconstruction patch for the element Ki in the following description. Once the patch is available, it is ready to do the solution reconstruction in each element. In this paper, the linear reconstruction is considered, and the least square method is used. More
7
specifically, in the element Ki , a linear variation of the conservative variable is supposed, i.e., ¯i + a1 (x − xi ) + a2 (y − yi ), P (x, y) = U
(11)
¯i is the cell average of the conservative quantity on the element Ki , and (xi , yi ) is the barycenter of where U Ki . It is noted that for a linear function, its evaluation on the barycenter of a given element is equal to its average on the element. Then with a given patch Pi = {Kj }, the gradient (a1,i , , a2,i )T in (11) is calculated by solving arg min a1,i ,a2,i
¯j − U ¯i − a1,i (xj − xi ) − a2,i (yj − yi )||22 . ||U
(12)
Kj ∈Pi ,j=i
Limiting process is necessary for the second order or higher order solution reconstruction to restrain the possible nonphysical oscillation around the shock. Due to its ability on delivering nonoscillatory and high order numerical solutions, WENO method has been widely used in a variety of practical simulations. In this paper, we follow [19, 21] to design the WENO process as follows. First, it is assumed that there are M (m)
(m)
candidate gradients {a1,i , a2,i , m = 1, 2, . . . , M } in the element Ki after the solution reconstruction. For example for the element K0 in Fig. 1 (left one), four candidate gradients can be obtained by solving (12) for K0 on the patches {K0 , K1 , K2 , K3 }, {K0 , K1 , K4 , K5 },{K0 , K2 , K8 , K9 }, and {K0 , K3 , K11 , K12 }, if the first kind of patch is used. Then for each candidate gradient, a smoothness indicator is defined as follows, (m) (m) (m) (a1,i + a2,i )dxdy, m = 1, 2, . . . , M. Si = Ki
(m) ω ˜i
(m) Si ) β ,
(m)
= 1/( + the weight assigned to each candidate gradient is designed as ωi = By introducing (m) (k) ˜ i . Here β and are two constants, and we use β = 2 and = 1.0e − 04 in the simulations. The ω ˜i / k ω final reconstructed gradient on the element Ki is given by T (m) (m) (m) (m) T ωi a1,i , ωi a2,i . (a1,i , a2,i ) = m
m
For the WENO reconstruction for the second kind of patches, i.e., P2 and P2B , it follows the above calculations. We refer to [22, 23] for the implementation detail. 3.3. h-Adaptive Methods In the simulation of the detonation wave, resolving the detonation front effectively and dynamically is crucial for a quality simulation. Since the variation of the conservative quantities in this region is relatively much larger than that in the other region, a uniform dense mesh in the whole domain is not appropriate from the implementation efficiency point of view. Consequently, mesh adaptive methods become the possible choice to resolve this issue. In this paper, we consider to use the local refinement technique to handle the adaptive process. This technique has been already adopted for the steady Euler equations, see [20, 22, 23, 28] and references 8
therein. For the h-adaptive technique, an efficient management strategy on handling the local refinement and coarsening of the mesh is crucial for the efficiency of the implementation. In the previous work, we use the hierarchy geometry tree (HGT) data structure, which is based on the fork-tree data structure in 2D case, to handle the mesh refinement and coarsening. Briefly speaking, the original mesh consists of a set of root elements. During the simulation, a series of the subtrees for those root elements would be generated according to the error indicators produced by the numerical solutions. By collecting the leaf nodes of all subtrees, a nonuniform mesh will be generated for the computational domain. With the HGT, the CPU time consumed on the local refinement and/or coarsening of the mesh is nothing compared with the whole CPU time. Please refer to [26] for detail of HGT, and also [4, 6] for 3D case. In the following, we will try to answer the other two questions on the quality realization of the mesh adaptive methods for the detonation wave simulation, i.e., where to refine the mesh, and how to update the conservative solutions. The first question mentioned above is related to the generation of the error indicator. During the adaptive refinement process, an error indicator will be assigned to each element. Compared with a given tolerance, it would be determined that whether an element will be refined into four subelements, and whether a set of subelements will be cast into a big element. In the detonation simulation, there is a shock for the pressure around the detonation front, and the peak pressure around the detonation front is also a very important quantity. Hence, it seems like that the gradient of the pressure would be a natural choice for the indicator generation. However, since the variation of the pressure around the detonation front is pretty dramatic, the gradient of the pressure would be very sensitive to any small numerical perturbation. Although the perturbation can be reduced systematically with reducing the mesh size, the quality of the numerical solution can still be affected since the existence of the transverse detonation instability for a two dimensional simulation. In the implementation of the algorithm, if the gradient of the pressure is employed to deliver an error indicator, some clustered locally refined mesh grids would be generated along the transverse direction of the detonation front, which might trigger the transverse instability numerically. Hence, to provide a uniform simulation environment is important, and a possible way to resolve the issue is to resort to some methods such as level set method [3] and phase field method [34] for the interface tracking, then the indicator for each element is generated with the knowledge of the interface. In this paper, we focus on the semi-2D and 2D simulations for the Zel‘dovich-Neumann-D¨oring (ZND) models. Hence, the structure of the shock is almost one dimensional with perturbations caused by the transverse instability of the simulations. Since the shock position can be predicted in advance [29, 30, 40], it is a reasonable choice to use skew normal distribuction function to generate a distribution of the error indicator along the x-direction. The skew normal distribution is given as f (x|α, μ, σ 2 ) = 2Φ(x|α)φ(x|μ, σ 2 ),
9
(13)
where φ(x|μ, σ 2 ) is the normal distribution φ(x|μ, σ 2 ) = √
1 2σ 2 π
exp (−
(x − μ)2 ), 2σ 2
and Φ(x|α) is cumulative distribution function given by 2 1 1 + erf ( √ ) . Φ(x|α) = 2 αx
(14)
(15)
In the above formulas, σ 2 and μ stand for the variance and the mean of the distribution, respectively, and α is a parameter to control the skew of the distribution. In the simulations, the parameter μ is chosen as the position of the leading shock which is given by μ = max x
Px , ρ
(16)
and the parameters σ 2 = 64, and α = 10. After the skew normal distribution is determined, the component of the error indicator for the element Ki from the normal distribution is given by 2 i indK nd = |Ki |f (xi |α, μ, σ ),
(17)
where |Ki | is the area of the element Ki , and the xi is the x-component of the centroid of the element Ki . There are two advantages on using the above refinement indicator. First, by choosing the position of the leading shock as μ, the indicator would reach the maximum value around the leading shock, which means that the mesh grids in this area would be much denser than that in other area. Second, the distribution function would reach zero slowly and smoothly, which means that the mesh grids would become coarse gently. This is positive for the simulation since solutions in the reaction zone can also be resolved well. The excellent results can be observed in the last section by using this kind of indicators. For the solution update from the old mesh to the new mesh, we need to take care of two things, i.e., conservation of the cell average and the high resolution of the numerical solution. With the HGT data structure, this process could be realized easily. Let To and Tn denote the old mesh and new mesh, respectively. From the property of the HGT data structure, there are following four relations between each two elements from two meshes, i.e., for Ko ∈ To and Kn ∈ Tn , we have the following relations, K o ≺ Kn ,
Ko K n ,
Ko = Kn , and Ko has nothing to do with Kn .
Here Ko ≺ Kn means Ko is a child of the element Kn , and Ko Kn means Ko is the parent element of Kn , and Ko = Kn means that two elements are the same node in the HGT. Suppose that the relation between two elements Ko and Kn is known, the solution update from Ko to Kn will be implemented as follows in this paper.
10
¯n + = U ¯o ∗ |Ko |/|Kn |; • Ko ≺ Kn : U ¯ n = Uo (X K ); • Ko K n : U n ¯n = U ¯o , • Ko = K n : U K means the barycenter of the element K. Uo means the where |K| means the area of the element K, and X reconstructed solution. It is noted that with HGT data structure, the relation between each two elements from two meshes can be obtained very efficiently, which results in an efficient implementation on the solution update. 3.4. OpenMP Parallelization The h-adaptive methods introduced in last subsection could effectively enhance the algorithm efficiency by controlling the amount of mesh girds in the simulation. To further improve the algorithm efficiency, the OpenMP parallelization technique will be considered in this subsection. OpenMP[1] is an API that supports shared memory multiprocessing programming in C/C++. Since the library for the numerical simulations in this paper is developed by C++, and the hardware is a DELL Precision T5610 workstation with 12 CPU Cores and 64 Gigabytes memory, it is a very natural choice to accelerate the implementation by parallelizing the algorithm with OpenMP. In the programming of the algorithm with OpenMP, an important issue is to resolve the data dependence among those threads. A straightforward way is to use the the clauses provided by OpenMP such as critical and atomic to avoid the data race. However, the efficiency of the implementation would be affected negatively, since in some sense the serial operations are forced with those clauses. Hence, a better way is to logically separate the write operations, then the calculation task can be assigned to every thread directly without data race. It is suggested that the edge/face-based data structure is preferred for the finite volume discretization from the efficiency point of view, since each edge in the mesh would be visited exactly once, instead of twice for an element-based data structure, in the numerical flux calculation and the cell average update. However, it is different to parallelize the above two operations. In the numerical flux calculation, there are two read operations to extract the left and right state of the conservative quantities from two adjacent elements, and only one write operation to save the calculated flux value on the current edge. Consequently, the for loop for the calculations of the numerical flux on the edges can be parallelized directly. The situation is totally different for the cell average update operation. In the update process, there is one read operation to extract the numerical flux from the current edge, and two write operations on two adjacent elements to accrue the change of the cell average. They are the two write operations which result in the data dependence among those threads, so the direct parallelization of the for loop for the cell average update is no longer acceptable. To resolve this issue, one may either generate a series of sets of independent edges first, then parallelize the 11
edge-based cell average update in each set, or implement an element-based cell average update instead. By element-based operation here, it means that we loop the element in the domain, then loop the edges in the element to calculate the new cell average. With such operation, except for the edge located on the domain boundary, each edge will be visited twice. However, the data dependence issue would be resolved effectively. In this paper, we use the element-based cell average update in the parallel simulation. For the other main components in the finite volume methods such as the solution reconstruction and the implementation of the Runge-Kutta module, either there is no data dependence, or the data dependence can be resolved trivially.
4. Numerical Tests The purpose of the numerical tests in this paper are two-fold. First, the numerical convergence of our method will be examined by a variety of benchmark problems, and the Zel‘dovich-Neumann-D¨ oring (ZND) model will serve the examination. Then the effectiveness of the proposed h-adaptive method will be demonstrated. 4.1. Initial Condition In this paper, we follow [15] to use the solutions of 1D ZND model as the initial conditions of our 2D simulations. The following is a brief description on how to generate the initial solutions for the 2D simulations with specific setup. Suppose that the distribution of the mass fraction f1 = ρr /ρ in the domain is known, and the x-component of the initial position of the detonation wave front is denoted by x0df . Then in the flame zone, the x-component velocity is given by γ u=− γ+1 where
ξ(f1 ) =
1 γ (D + ) γ+1 D
1 D+ D
2 −2
γ−1 γ+1
+
ξ(f1 ),
q0 (1 − f1 ) +
(18)
D2 γ + γ−1 2
.
Once u is obtained, the density and the pressure are calculated by the following formulas. ρ=−
D , u
P = Du + D2 + 1.
(19)
In the above two formulas (18) and (19), D is the detonation wave velocity, and is derived from the Chapman2 Jouguet velocity DCJ [13] and an overdriven factor f by 2 D2 = f DCJ ,
2 where DCJ = γ + (γ 2 − 1)q0 +
(20)
2
(γ + (γ 2 − 1)q0 ) − γ 2 .
In summary, once the distribution of the mass fraction f1 is known, the other quantities such as the x-component velocity u, density ρ and the pressure P can be deduced from (18), (19). Together with (24) 12
and the equation of state (3), all quantities in the reactive Euler equations (1) are obtained. Hence, the only thing left is how to find the distribution of f1 . In [41], the mass fraction f1 is generated by solving the following ODE, Kx0df f1 −Ea /T df1 =− e , dx u
(21)
ρ = 1, P = 1, T = 1, and u = −D,
(22)
with the boundary condition,
which denote the quantities in the unburned zone. Instead, in the simulations in this paper, the following function is adopted to simply approximate the distribution of the mass fraction f1 , ⎧ 1 ⎪ ⎪ tanh(π(x − x0df )) + 1 , x < x0df , ⎪ ⎨ 2 f1 = ⎪ ⎪ ⎪ ⎩ 1, elsewhere.
(23)
Although (23) is only a rough approximation for the solution of ODE (21), it can successfully give a jumpstart of the detonation wave simulation. We refer to [14] for a different f1 . For the simulations which need an initial perturbation along the y-direction, the y-component of the velocity is given by the following transversely sinusoidal perturbed planar ZND wave, i.e., ⎧ ⎪ ⎪ A sin (2k1 πy/L), x0df − Ww ≤ x ≤ x0df , ⎪ ⎨ v= ⎪ ⎪ ⎪ ⎩ 0, elsewhere,
(24)
where A, k1 , and L denote the amplitude, wave number and the wavelength of the perturbation, and Ww is the width of the perturbed zone along the x-axis. 4.2. Boundary conditions The computational domain for the following numerical simulations is a rectangle channel. The setup of the flow of the fluid is as follows. The reactant flows into the domain from the right boundary, and the production leaves the domain from the left boundary. The slip wall is employed for the lower and upper boundaries of the domain. Under the assumption of the ZND model, the velocity of the fluid in the domain is shifted by the negative Chapman-Jouguet velocity. Hence, the subsonic inflow boundary condition and subsonic outflow boundary condition are adopted on the right and left boundaries of the domain, respectively. The reflective boundary condition is adopted for both lower and upper boundaries of the domain to simulate the solid wall. To impose the subsonic inflow boundary, the classical way of using Riemann invariant can be used directly. However, since the configuration of the fluid in the far field outside the outflow boundary is not 13
clear, it brings difficulty to set up the boundary condition when the physical boundary condition is needed. To resolve this issue, several strategies are developed to make the numerical simulations reliable. In [30], an extrapolation technique is employed, where the conserved quantities of the system are relaxed to the condition at infinity with a given formula. In this case, the subsonic outflow boundary condition can be imposed trivially. In [15], the perfectly matched layer (PML) method is introduced in the simulation, and excellent numerical results are obtained. Besides the above two strategies, a straightforward way to handle the above issue is to use a sufficiently long computational domain along the x-axis, then the boundary values for those conservative quantities are obtained by the extrapolation. A formula for the relation between the length of the domain and the time can be found in [24, 29]. In the numerical simulations in this paper, we will try to use sufficiently long domain to restrain the nonphysical boundary effect. 4.3. Semi-2D ZND simulations The purpose of the numerical examples presented in this subsection is to verify the validity of our algorithm. The parameters used in the simulations are given in Table 1. To unify the length of the halfreaction zone L1/2 , the pre-exponential factor K is chosen as 145.69 and 230.75 for stable and unstable cases, respectively. γ
E0
q0
f
Stable
1.2
50.0
50.0
1.8
Unstable
1.2
50.0
50.0
1.6
Table 1: Parameters used in the semi-2D ZND simulations.
The parameters in Table 1 are used in previous works [13, 14, 29] to test the one dimensional ZND model. In this paper, we focus on numerical methods for two dimensional reactive Euler equations. Hence, a semi2D configuration for the numerical experiments are adopted in the simulations, i.e., the velocity component v in the equations is initialized as 0 in the whole domain, while the other quantities are initialized by the analytical solutions for the one dimensional ZND model. With this configuration, theoretically, the numerical solutions obtained from the semi-2D problems should be consistent with the solutions from one dimensional problems. However, due to its intrinsic instability for the multi-dimensional detonation wave, it is nontrivial to reproduce 1D results. From the stability analysis presented in [36], it is known that only a discrete set of the transverse disturbance whose wavelengths are compatible with the channel width can exist when the solid wall boundary condition is adopted, and the wavelengths have the following relations with the width of the channel, wi = 2W/i, i = 1, 2, 3, . . . ,
14
where wi and W stand for the wavelengths of the transverse disturbance and the channel width, respectively. From the above formula, it can be seen that the largest wavelength is w1 = 2W . Hence, if the number of the mesh grids along the y-direction is not big enough, the disturbance introduced by the discretization error along the y-direction will not be resolved well, and our semi-2D simulations will reproduce the numerical results given by 1D ZND model. Based on the above discussion, the computational domain used in the following two numerical experiments is [0, 130] × [−0.0433, 0.0433]. The initial mesh for this domain is given in Fig. 2. The mesh size is around 0.1. To study the numerical convergence, the simulations need to be implemented on a successively refined meshes. To remove the transverse instability when the mesh grids along y-direction become dense, the uniform refinement of the mesh grids is abandoned. Instead, the total length of the domain along x-direction and the regularity of the triangle element is kept, only the mesh size becomes 0.05. In this case, the mesh is similar to the coarse one, and the domain now becomes [0, 130] × [−0.0217, 0.0217]. The successively refined meshes in this subsection are generated by this strategy. In the following simulations, the initial position of the detonation front is xdf = 80.
Figure 2: Three fractions of the initial mesh used in Subsection 4.3.
For the simulations with the stable parameters given in Table 1, the results obtained on three successively refined meshes are given in Fig. 3. Since the parameter K = 145.69 makes L1/2 the length of the half-reaction zone the unit, there are around 10 mesh grids in the zone along the x-direction for the initial mesh, and 20 and 40 for the two refined meshes. In all three simulations, the expected stable state of the detonation is obtained smoothly. The numerical convergence can also be observed from the results successfully, i.e., the front pressure approaches the reference value 75.79 [8] by refining the mesh. 85
80 75 70
85
Ref. mesh size 0.05
Front Pressure
mesh size 0.1 Ref.
Front Pressure
Front Pressure
85
80 75 70
0
20
40
60 Time
80
100
Ref. mesh size 0.025
80 75 70
0
20
40
60 Time
80
100
0
20
40
60
80
100
Time
Figure 3: Mesh refinement study of the pressure history at the shock front for the stable detonation wave simulations. Reference value is 75.79[8].
Two remarks are given below for the results shown in Fig. 3. First, as we can observe from the results, there is a large oscillation of the pressure at the beginning of the simulations. This can be explained by the 15
use of the initial distribution of f1 given by (23), which is actually only an approximation for the 1D ZND model. This inconsistency introduces the large oscillation. Second, the oscillation in the numerical solution disappears slower in the dense mesh case than that in the coarse mesh case. This can be partially explained as follows. With denser mesh, the inconsistency between the exact ZND solution and the initial function can be resolved better. Since there is no instability mode for the simulation with stable parameters, all perturbations introduced by the initial condition go to zero eventually. In Fig. 4, the numerical results obtained from four successively refined meshes with the unstable parameters given in Table 1 are shown. The initial mesh is given by Fig. 2 with the mesh size 0.1, and the mesh sizes for three successively refined meshes are 0.05, 0.025, and 0.0125, respectively. With the pre-exponential factor K = 230.75, the length of the half-reaction zone is unit, which means that there would be 10, 20, 40, and 80 mesh grids in the region for the four meshes, respectively. Corresponding to these parameters, it is a 1D unstable detonation, and the reference top pressure is 101.1 ± 0.2[11, 29]. From Fig. 4 (top left), it can be seen that the simulation fails to present pulsation structure of the front pressure. Actually, the solution becomes stable with the time evolution. This can be explained that the mesh is too coarse to preserve the instability mode of the detonation wave. From the following three simulations implemented on three successively refined meshes, the periodic unstable structures are successfully observed, and the peak pressure approaches the reference value with the refinement of the mesh grids, which shows the numerical convergence of the method for this case.
Figure 4: Mesh refinement study of the pressure history at the shock front for the stable detonation wave with the parameters in the second row of Table 1.
4.4. Cell structure simulation The numerical test in this subsection is a two dimensional case, and the simulation will show that the regular cell structure would be generated after the detonation front, with the parameters given in Table 2. In the simulations, xdf = 62, and K = 1134363.64. To initially perturb the velocity v along the y-axis, the formula (24) with the parameters A = 0.1, L = 10, k1 = 4, and Ww = 1 are used. 16
γ
E0
q0
f
1.2
20.0
2.0
1.1
Table 2: Parameters used in the cell structure simulation.
First, the following three computational domains are tested to check the boundary effect on the simulations, i.e., [0, 34.64] × [−5, 5], [0, 69.28] × [−5, 5], and [0, 103.92] × [−5, 5]. Fig. 5 shows the initial mesh on the domain [0, 69.28] × [−5, 5].
Figure 5: Computational domain [0, 69.28] × [−5, 5].
In all three initial meshes, the side length of the equilateral triangle are all equal to 11.55. In the first test, each initial mesh are uniformly refined 7 times, i.e., the side length becomes 0.09. Then the simulation on each mesh runs till t ≈ 200. The record of front pressure of each simulation is presented in Fig. 6. From the figure, it can be observed that for the domain with length 34.64, the peak pressure keeps decaying from around t = 120 (upper one), while the periodic pulsation structure is preserved well for the domains with length 69.28 and 103.92, respectively (lower one). Furthermore, it is found that there is almost no difference between the records of the front pressure from two domains (Fig. 6, lower one). Hence, in the following simulations, the domain [0, 69.28] × [−5, 5] and the initial mesh shown in Fig. 5 are used. Remark 1. Since its intrinsic instability for the multi-dimensional detonation wave, small perturbation along the transverse direction of the leading shock may cause the different behavior of the numerical solution. To remove the numerical perturbation introduced from the solution reconstruction due to the nonuniform distribution of the mesh grids, the equilateral triangle meshes are employed in the simulations. To demonstrate the influence of the mesh grids on the numerical solutions, the simulations with the parameters given in Table 2 are reimplemented in the domain [0, 34.64] × [−5, 5] with two difference initial meshes shown in Fig. 7. The mesh in Fig. 7 (left one) is uniformly refined 7 times, while it is 5 times for the right mesh in the figure so that the mesh sizes for both meshes are almost the same. The results on the front pressure are shown in Fig. 8, in which it can be observed that formation of the periodic pulsation structure with the mesh generated by EasyMesh[32] is much earlier than that generated by the equilateral triangle mesh. It can also be observed that with the refinement of the mesh, the formation of the pulsation structure is delayed, but still earlier than that for equilateral triangle mesh case. This obviously shows the sensitivity of the numerical solution to the perturbation. In the following simulations, the mesh shown in Fig. 5 is used as the initial mesh. In the first simulation, the mesh is generated from the initial mesh by uniformly refining 6 times, the side length of the triangle is around 0.18. Then the simulations are reimplemented on three successively refined meshes. The records of 17
Font Pressure Font Pressure
6.5 5.5 4.5 3.5 2.5
Domain [-5,5] X [0, 34.64]
6.5 5.5 4.5 3.5
Domain [-5,5] X [0, 69.28] Domain [-5,5] X [0, 103.92]
2.5 0
20
40
60
80
100
120
140
160
180
200
Time Figure 6: The records of the front pressure obtained from simulations on three meshes with different length.
Figure 7: Two different meshes in the domain [0, 34.64] × [−5, 5]. Left one is an equilateral triangle mesh, and right one is generated by EasyMesh[32].
.
Front Pressure
7 Equilateral with 57921 points EasyMesh with 41121 points EasyMesh with 197569 points
6 5 4 3 0
20
40
60
80
100
Time Figure 8: The records of the front pressure obtained from equilateral triangle mesh and meshes generated from EasyMesh.
18
the front pressure obtained from first three simulations are shown in Fig. 9, and the pressure distribution and the vorticity for each simulation at t ≈ 30 are shown in Fig. 10. From Fig. 9, the intrinsic instability of the multi-dimensional detonation is shown obviously. Although with the equilateral triangle mesh, the stable state of the detonation can last a long time, and with the refinement of the mesh, it can last longer, the periodic pulsation structure of the front pressure eventually appears. Fig. 10 shows the pressure and the vorticity at the time t ≈ 30 before the formation of the pulsation structure. It can be seen that the regular cell structure is generated successfully, and with the denser mesh, the resolution of the structure becomes better. In Fig. 11, the pressure and vorticity at t ≈ 130 are shown. In this time, the periodic pulsation structure has been formed already.
Front Pressure
7
27105 mesh grids 107457 mesh grids 427905 mesh grids
6 5 4 3 0
20
40
60 Time
80
100
120
Figure 9: The records of the front pressure obtained from simulations on three meshes with different sizes.
4.5. Simulations with h-adaptivity In this subsection, the adaptive module of the proposed method is tested. For the convenience on the comparison, the parameters in Table 2 are still used in this subsection, and the mesh shown in Fig. 5 is used as the initial mesh. The other parameters used in the simulation are same to that in the last subsection. The results obtained from the simulation are shown in Figs. 12 and 13. In Fig. 12, the record of the front pressure from adaptive method is shown. As a comparison, the records obtained from two successively refined meshes are also shown in the figure. It can be seen that in the adaptive method, there are 85254 mesh grids in the domain, and the record of the front pressure from the adaptive method almost duplicate that from the uniform mesh which has 427950 mesh grids (subfigure in Fig. 12). This is reverified by the pressure and vorticity in Fig. 13 (top left and top middle, receptively) that they coincide with the ones shown in Fig. 10 (third row) very well. In this case, over 80% mesh grids are saved with the help of the adaptive method to reach almost the same accuracy. Fig. 13 (bottom one) shows the distribution of the mesh grids in the whole domain, while the Fig. 13 (top right one) shows the distribution of the mesh grids around the leading shock region. It can be observed that with the indicators generated by (17), most mesh grids are 19
Figure 10: The pressure (left column) and the vorticity (right column) obtained from simulations on four successively refined meshes at t ≈ 30.
20
Figure 11: The pressure (left column) and the vorticity (right column) obtained from simulations on the mesh with 427905 grid points at t ≈ 130.
Front Pressure
7
Adaptive mesh with 85254 points Equilateral with 107457 points Equilateral with 427905 points
3.8 3.7 3.6
6 5
15 17 19 21 23 25
4 3 0
10
20
30
40 Time
50
60
70
80
Figure 12: The record of the front pressure obtained from the adaptive method, and the comparison to the records obtained from two successively refined equilateral triangle meshes. The subfigure shows the detail of the comparison.
21
Figure 13: The pressure (top left), the vorticity (top middle), the local mesh (top right) around the leading shock, and the mesh (bottom) in the whole domain, obtained at t ≈ 30.
located around the leading shock region and the reaction region, which guarantees the numerical accuracy of the simulations. There are only a small amount of the mesh grids in other regions, which guarantees the efficiency on the simulation effectively.
Conclusion In this paper, an adaptive finite volume method is proposed for the reactive Euler equations. The numerical discretization includes a second order TVD Runge-Kutta method for the temporal discretization, and a finite volume method with linear solution reconstruction for the spatial discretization. To restrain the potential numerical oscillation introduced by the high order solution reconstruction, the WENO methodology is employed as the limiting strategy. To enhance the efficiency of the algorithm, the h-adaptive methods are introduced in this work. The indicator and the solution update method are developed from the consideration of the efficiency and the numerical accuracy. The OpenMP technique is also introduced in the algorithm to further improve the efficiency. A variety of the numerical tests are implemented to show the numerical convergence of the proposed method. The effectiveness of the h-adaptive method is also successfully demonstrated. The unstructured meshes have advantage for the practical simulations when the geometry of the domain is complicated. However, the numerical results in this paper show that the numerical behavior of the detonation simulation is sensitive to the perturbation introduced by the solution reconstruction on the unstructured meshes. Hence, a quality solution reconstruction method for the unstructured mesh grids is under considering. The other issue is related to the error indicator for the mesh adaptation. Although
22
the strategy used in this paper works well, a fine error indicator is needed to handle more complicated simulations, and to further enhance the efficiency. The related results will be reported in the forthcoming paper. Acknowledgements The work of G. H. Hu was partially supported by 050/2014/A1 from FDCT of Macao S.A.R., MYRG201400109-FST from University of Macau, and National Natural Science Foundation of China (Grant No. 11401608). The author would like to thank Prof. W. S. Don, Prof. Z. Gao, Mr. P. Li, and Prof. R. Li for their helpful suggestions and discussion on this work. [1] http://openmp.org. [2] Marco Arienti. A numerical and analytical study of detonation diffraction. PhD thesis, California Institute of Technology, 2003. [3] Tariq D. Aslam, John B. Bdzil, and D.Scott Stewart. Level set methods applied to modeling detonation shock dynamics. Journal of Computational Physics, 126(2):390 – 409, 1996. [4] Gang Bao, Guanghui Hu, and Di Liu. An h-adaptive finite element solver for the calculations of the electronic structures. Journal of Computational Physics, 231(14):4967 – 4979, 2012. [5] Gang Bao, Guanghui Hu, and Di Liu. Numerical solution of the Kohn-Sham equation by finite element methods with an adaptive mesh redistribution technique. Journal of Scientific Computing, 55(2):372– 391, 2013. [6] Gang Bao, Guanghui Hu, and Di Liu. Real-time adaptive finite element solution of time-dependent Kohn-Sham equation. Journal of Computational Physics, 281:743–758, 2015. [7] Gang Bao, Guanghui Hu, and Di Liu. Towards translational invariance of total energy with finite element methods for Kohn-Sham equation. Communications in Computational Physics, 19:1–23, 2016. [8] Ralf Deiterding. Parallel adaptive simulation of multi-dimensional detonation structures. September 2003. [9] Ralf Deiterding. A parallel adaptive method for simulating shock-induced combustion with detailed chemical kinetics in complex domains. Computers & Structures, 87(11–12):769 – 783, 2009. Fifth MIT Conference on Computational Fluid and Solid Mechanics. [10] Hua-Shu Dou, Her Mann Tsai, Boo Cheong Khoo, and Jianxian Qiu. Three-dimensional numerical simulation for detonation waves using WENO schemes. American Institute of Aeronautics and Astronautics, 45th AIAA Aerospace Sciences Meeting and Exhibit, 2007. 23
[11] Jerome J. Erpenbeck. Stability of idealized one-reaction detonations. Physics of Fluids, 7(5):684–696, 1964. [12] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). Journal of Computational Physics, 152(2):457 – 492, 1999. [13] Wildon Fickett and Willam C. Davis. Detonation: Theory and Experiment. Dover Publications, 2011. [14] Zhen Gao, Wai Sun Don, and Zhiqiu Li. High order weighted essentially non-oscillation schemes for onedimensional detonation wave simulations. Journal of Computational Mathematics, 29:623–638, 2011. [15] Zhen Gao, WaiSun Don, and Zhiqiu Li. High order weighted essentially non-oscillation schemes for two-dimensional detonation wave simulations. Journal of Scientific Computing, 53(1):80–101, 2012. [16] Sigal Gottlieb and Chi-Wang Shu. Total variation diminishing Runge-Kutta schemes. Mathematics of Computation, 67:73–85, 1998. [17] P. He and H.Z. Tang. An adaptive moving mesh method for two-dimensional relativisitic hydrodynamics. Communications in Computational Physics, 11:114–146, 2012. [18] William D Henshaw and Donald W Schwendeman. An adaptive numerical scheme for high-speed reactive flow on overlapping grids. Journal of Computational Physics, 191(2):420 – 447, 2003. [19] Changqing Hu and Chi-Wang Shu. Weighted essentially non-oscillatory schemes on triangular meshes. Journal of Computational Physics, 150:97–127, 1999. [20] Guanghui Hu. An adaptive finite volume method for 2D steady Euler equations with WENO reconstruction. Journal of Computational Physics, 252(0):591 – 605, 2013. [21] Guanghui Hu, Ruo Li, and Tao Tang. A robust WENO type finite volume solver for steady Euler equations on unstructured grids. Communications in Computational Physics, 9:627–648, 2011. [22] Guanghui Hu, Xucheng Meng, and Nianyu Yi. Adjoint-based an adaptive finite volume method for steady Euler equations with non-oscillatory k-exact reconstruction. Computers & Fluids, 139:174 – 183, 2016. [23] Guanghui Hu and Nianyu Yi. An adaptive finite volume solver for steady Euler equations with nonoscillatory k-exact reconstruction. Journal of Computational Physics, 312:235 – 251, 2016. [24] P Hwang, R P Fedkiw, B Merriman, T D Aslam, A R Karagozian, and S J Osher. Numerical resolution of pulsating detonation waves. Combustion Theory and Modelling, 4(3):217–240, 2000. 24
[25] Smadar Karni. Hybrid multifluid algorithms. SIAM Journal on Scientific Computing, 17(5):1019–1039, 1996. [26] Ruo Li. On multi-mesh h-adaptive methods. Journal of Scientific Computing, 24(3):321–341, 2005. [27] Ruo Li and Tao Tang. Moving mesh discontinuous Galerkin method for hyperbolic conservation laws. Journal of Scientific Computing, 27(1-3):347–363, 2006. [28] Ruo Li, Xin Wang, and Weibo Zhao. A multigrid block LU-SGS algorithm for Euler equations on unstructured grids. Numerical Mathematics: Theory, Methods and Applications, 1:92–112, 2008. [29] Y.S. Lian and K. Xu. A gas-kinetic scheme for multimaterial flows and its application in chemical reactions. Journal of Computational Physics, 163(2):349 – 375, 2000. [30] Z. Liang and L. Bauwens. Cell structure and stability of detonations with a pressure-dependent chainbranching reaction rate model. Combustion Theory and Modelling, 9(1):93–112, 2005. [31] E. Loth, S. Sivier, and J. Baum. Adaptive unstructured finite element method for two-dimensional detonation simulations. Shock Waves, 8(1):47–53, 1998. [32] Bojan Niceno. EasyMesh. http://web.mit.edu/easymesh_v1.4/www/easymesh.html. [33] J.T. Oden and L. Demkowicz. h-p adaptive finite element methods in computational fluid dynamics. Computer Methods in Applied Mechanics and Engineering, 89(1–3):11 – 40, 1991. Second World Congress on Computational Mechanics. [34] Y. Sun and C. Beckermann. Sharp interface tracking using the phase-field equation. Journal of Computational Physics, 220:626–653, 2007. [35] Tao Tang. Moving mesh methods for comutational fluid dynamics. Contemporary Mathematics, American Mathematical Society, 383:141–173, 2005. [36] B. D. Taylor, A. R. Kasimov, and D. S. Stewart. Mode selection in weakly unstable two-dimensional detonations. Combustion Theory and Modelling, 13(6):973–992, 2009. [37] Fumiya Togashi, Rainald L¨ ohner, and Nobuyuki Tsuboi. Numerical simulation of h2/air detonation using unstructured mesh. Shock Waves, 19(2):151–162, 2009. [38] Eleuterio F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer Berlin Heidelberg, 2009.
25
[39] Cheng Wang, Xiangxiong Zhang, Chi-Wang Shu, and Jianguo Ning. Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations. J. Comput. Phys., 231(2):653–665, January 2012. [40] Daniel N. Williams, Luc Bauwens, and Elaine S. Oran. Detailed structure and propagation of threedimensional detonations. Symposium (International) on Combustion, 26(2):2991–2998, 1996. [41] Zeng-Chan Zhang, S.T. John Yu, Hao He, and Chang Sin-Chung. Direct calculations of two- and threedimensional detonations by an extended CE/SE method. 39th Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, 2001.
26
Guanghui Hu
PII: DOI: Reference:
S0021-9991(16)30633-7 http://dx.doi.org/10.1016/j.jcp.2016.11.041 YJCPH 6991
To appear in:
Journal of Computational Physics
Received date: Revised date: Accepted date:
4 December 2015 6 November 2016 28 November 2016
Please cite this article in press as: G. Hu, A numerical study of 2D detonation waves with adaptive finite volume methods on unstructured grids, J. Comput. Phys. (2016), http://dx.doi.org/10.1016/j.jcp.2016.11.041
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A numerical study of 2D detonation waves with adaptive finite volume methods on unstructured grids Guanghui Hua,b a Department b UM
of Mathematics, University of Macau, Macao S.A.R., China Zhuhai Research Institute, Zhuhai, Guangdong Province, China
Abstract In this paper, a framework of adaptive finite volume solutions for the reactive Euler equations on unstructured grids is proposed. The main ingredients of the algorithm include a second order total variation diminishing Runge-Kutta method for temporal discretization, and the finite volume method with piecewise linear solution reconstruction of the conservative variables for the spatial discretization in which the least square method is employed for the reconstruction, and weighted essentially nonoscillatory strategy is used to restrain the potential numerical oscillation. To resolve the high demanding on the computational resources due to the stiffness of the system caused by the reaction term and the shock structure in the solutions, the h-adaptive method is introduced. OpenMP parallelization of the algorithm is also adopted to further improve the efficiency of the implementation. Several one and two dimensional benchmark tests on the ZND model are studied in detail, and numerical results successfully show the effectiveness of the proposed method. Keywords: reactive Euler equations; finite volume methods; unstructured grids; mesh adaptive methods;
1. Introduction The development of the numerical methods for the chemical reactive flows, which is an important application for multimaterial flows, has attracted considerable attention in the last few years[12, 25, 29]. The reactive Euler equations are one of the most important models for describing the physical phenomenon such as detonation. To date, there have been a lot of numerical methods for solving reactive Euler equations. For example, the finite difference methods[10, 14, 15], the gas-kinetic schemes[29], the finite volume methods[18, 36], the discontinuous Galerkin schemes[39], etc. In all these works, excellent numerical results can be observed. However, the structured mesh grids are employed in all these works, which restrains their application in the practical problems such as the numerical simulations for the rotating detonation engine. In the practical problems, the computational domain always has complicated geometry, hence the unstructured mesh grids become a popular choice to handle the complicated domain. For some early numerical methods for reactive Euler equations on unstructured meshes, people may refer to [31, 37] for the finite element methods, etc.
Preprint submitted to Journal of Computational Physics
November 30, 2016
Besides the ability on handling complicated geometry of the domain, the other feature the numerical methods desired for practical applications is the efficiency. In the reactive Euler equations, the existence of the shock structure in the computational domain makes that the mesh size around the shock region should be small enough to resolve the shock well, which means that the size of the time step in the numerical methods should also be small enough to guarantee the simulation stability. Moreover, the stiffness of the reactive Euler system makes the choice of the size of the time step more severe. On the other hand, since the computational domain is bounded, it is necessary to remove or reduce the unphysical wave reflected by the outflow boundary, from the reliability point of view. To do so, several techniques have been developed. For example, in [30], an extrapolation technique is employed, where the conserved quantities of the system are relaxed to the condition at infinity with a given formula. In this case, the subsonic outflow boundary condition can be imposed trivially. In [15], the perfectly matched layer (PML) method is introduced in the simulation, and excellent numerical results are obtained. A more straightforward way is to use a sufficiently large domain. The artificial effect introduced by the wave absorbing strategies mentioned above can be effectively avoided with this strategy. However, the requirement on the size of the domain, together with the requirement on the sizes of the mesh grids and time step, will result in the highly demanding on the computational resource. This issue motivates the introduction of the adaptive mesh methods. There are three different kinds of mesh adaptive methods in the market. The h-adaptive methods refine and/or coarsen the mesh grids locally, and have been widely used in the computational fluid dynamics[9, 20, 28]. The r-adaptive methods redistribute the mesh grids in the domain to better resolve the singularity in the solutions, people may refer to [5, 17, 27, 35] for the detail of the methods and their applications in different areas. The p-adaptive methods locally upgrade the order of the approximate polynomial in the element. In the practice, p-adaptive methods always work together with the h-adaptive methods, people may refer to [33] for their applications in CFD. For the simulations of the detonation wave with the mesh adaptive methods, there have been some pioneer works[2, 9, 31]. In our previous work[7, 4, 6, 20], a general framework for h-adaptive mesh methods has been developed and applied for the computational fluid dynamics and density functional theory in 2D/3D cases. In this paper, such h-adaptive mesh framework will be introduced to solve the reactive Euler equations. To further improve the implementation efficiency, the OpenMP parallelization of the algorithm will also be discussed. The numerical results shown in the last section will demonstrate the effectiveness of the proposed method. In this paper, we develop a general framework of using adaptive finite volume methods for solving the reactive Euler equations on unstructured meshes. The numerical discretization of the Euler equations includes a second order total variation diminishing Runge-Kutta method for the temporal discretization, and finite volume method with piecewise linear solution reconstruction for the spatial discretization. To effectively remove the potential numerical oscillation introduced by the solution reconstruction, we follow [21] to use
2
WENO reconstruction. The reconstruction patch for the element which has an edge located on the boundary of the domain is designed specially to avoid the poor gradient reconstruction. To reduce the unphysical effect generated by the reflection of the wave from the outflow boundary, we simply use a sufficiently large domain [29]. The h-adaptive methods are employed in this paper to enhance the implementation efficiency. To further improve the efficiency, the algorithm is also parallelized with OpenMP, and detail on the implementation is discussed. A variety of numerical simulations successfully show the numerical convergence and effectiveness of our method on solving the reactive Euler equations. The rest of the paper is organized as follows. In the next section, the reactive Euler equations are introduced in detail. In Section 3, the numerical discretization of reactive Euler equations are derived. Then the solution reconstruction, the mesh adaptive methods, and some issues on OpenMP parallelization of our algorithm will also be introduced in depth. In Section 4, the initial and boundary conditions used in the simulations are introduced first, then numerical results from a variety of tests are demonstrated to show the numerical convergence and effectiveness of our method. Finally, we give the conclusion of this work, and also some future work on this topic.
2. Model Statement In this paper, we focus on the two dimensional reactive Euler equations with two reaction species, i.e., reactant and production species. There is no intermediate species created during the reaction, and such single one-step and irreversible reaction can be described by k
→ P roduct, Reactant − where k denotes certain reaction rate. By ignoring the viscosity of the fluid, the convection and reaction phenomenon of the system during the detonation can be modeled by the following reactive Euler equations, ∂ U + ∇ · F (U ) = S(U ), ∂t
(1)
where U , F (U ), and S(U ) denote the conservative variable, flux, and source term, respectively. Their expressions are given as follows, ⎡ ⎤ ⎡ ρ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ρr ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ U = ⎢ ρu ⎥ , F (U ) = ⎢ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎢ ρv ⎥ ⎢ ⎣ ⎦ ⎣ E
⎤ ρu,
ρv
ρr u,
ρr v
⎡
⎤ 0
⎥ ⎢ ⎥ ⎢ ⎥ ⎢ −ρr k ⎥ ⎢ ⎥ ⎢ ⎥ , and S(U ) = ⎢ 0 ρuv ⎥ ⎢ ⎥ ⎢ ⎢ 0 ρv 2 + P ⎥ ⎦ ⎣ (E + P )v 0
ρu2 + P, ρuv, (E + P )u,
3
⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎦
(2)
where ρ, ρr , u = (u, v)T , E, and P , stand for the density of the fluid, the density of the reactant, velocity, total energy, and pressure, respectively. The density of the fluid ρ, the density of the reactant ρr , and the density of the product ρp satisfy the relation ρ = ρr + ρp . In the source term S(U ), the non-zero term −ρr k represents the change of the reactant in the reaction. In this paper, the Arrhenium reaction rate is employed for the reaction rate. Hence, k has the following form k = Ke−Ea /T , where K is the pre-exponential factor, Ea is the activation energy, and T = P/ρ is the temperature (with gas constant normalized to 1 in this paper). To close the above system, the following equation of state is used E = P/(γ − 1) + 0.5ρ(u2 + v 2 ) + ρr q0 ,
(3)
where γ is the ratio of the specific heat, and q0 is the amount of heat released per unit mass by the reaction.
3. Numerical Discretization and Some Numerical Issues In this section, the numerical discretization for the reactive Euler equations is described in detail first. Then some numerical issues including the initial and boundary conditions, numerical stability, etc., as well as acceleration techniques such as h-adaptive methods and OpenMP parallelization are introduced. For the spatial discretization, the finite volume method is employed, not only because of its numerical conservative property which is very important for preserving the desired physical behavior of the numerical solution, but also because of its flexibility on handling the domain with complex geometry. To deliver a quality solution reconstruction, the linear reconstruction based on the least square method is employed, with a WENO strategy to restrain the possible nonphysical oscillation of the numerical solution around the shock. In the temporal discretization, the second order TVD Runge-Kutta scheme is adopted, which is also positive on removing the nonphysical oscillation. To describe the discretization specifically, let Ω denotes the computational domain of the problem. In tri this paper, we consider the unstructured meshes, and use T = {Ki }N i=1 to denote a triangulation of the
domain Ω, and Ntri is the total number of the triangle element in T . We suppose that there is no nontrivial overlap for every two triangle elements, i.e., the overlap part of two elements could be a vertex, an edge, or nothing, and that the union of the all elements is the domain Ω. The temporal domain we used in the NT −1 simulation is [0, Tend ], and its partition is denoted by {[tn , tn+1 ]}n=0 , where NT means the total number
of the segments in the interval [0, Tend ].
4
3.1. Runge-Kutta Finite Volume Discretization To derive the semi-discretized numerical scheme first, one takes the integral of governing equations (1) in a control volume Ki . Then the divergence theorem gives ∂ U dxdy + F (U ) · nds = S(U )dxdy, Ki ∂t ∂Ki Ki
(4)
where ∂Ki denotes the boundary of the i-th triangle element Ki in T , and n denotes the unit out normal vector w.r.t. ∂Ki . By introducing the cell average of conservative variable U and the source S(U ) in Ki , i.e., 1 ¯= 1 U U dxdy, and S¯ = S(U )dxdy, |Ki | Ki |Ki | Ki the equations (4) become 1 ∂ ¯ U+ ∂t |Ki |
¯ F (U ) · nds = S,
(5)
∂Ki
where |Ki | is the area of the element Ki . In the simulation, the flux F (U ) is replaced by the numerical flux F¯ , which is determined by a local Riemann problem with the initial conditions given by the solutions in the element Ki and its neighbor Kj . In our work, the Lax-Friedrichs flux is adopted for the numerical flux. With the numerical flux F¯ , the equations become 1 ∂ ¯ U = S¯ − ∂t |Ki |
F¯ · nds := RHS,
(6)
∂Ki
which is an ODE system now, and an ODE solver is needed for a full discretization. For the numerical simulation of gas dynamics with shock structures, non-oscillatory property of the numerical method is crucial not only for the numerical stability, but also for the reliability of the numerical solution. Consequently, we use the second order TVD Runge-Kutta method[16] in the temporal discretization, which results in the following full discretization of the reactive Euler equations (1) in the domain Ki × [tn , tn+1 ]. ⎧ ⎪ ¯ (1) = U ¯ n + δtRHS n , ⎪ U ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ U ¯ n + 0.5U ¯ (1) + 0.5δtRHS (1) , ¯ n+1 = 0.5U
(7)
¯ (1) . where RHS (1) means the evaluation of the RHS at U Besides the above Runge-Kutta finite volume method for reactive Euler equations, an alternative and competitive way is the splitting strategy[38] in which the two procedure, i.e., the advection of the reactant and the reaction, are supposed to happen in order. Since both the pure advection process and the stiff ordinary differential equations can be handled efficiently by classical methods, the splitting strategy may bring significant improvement on the efficiency. However, it has been shown in [38] that the splitting method 5
is only exact for some model problem. Further, it is nontrivial to design high order splitting methods. Hence, in this paper we use the above coupled process to describe the advection and reaction. To get a stable simulation, the selection of the size of the time step δt should follow the CourantFriedrichs-Lewy (CFL) condition. The idea behind the CFL condition is that for a given element, the incident wave is not allowed to move across the element in one single time step. To apply the CFL condition in the simulation, the wave speed and the cell size need to be determined. In terms of the characteristic analysis of the quasi-linear system derived from the reactive Euler equations (1), five eigenvalues obtained are λ1 = u − c, λ2 = λ3 = λ4 = u, and λ5 = u + c, where c denotes the sound speed. For the element size, the following approximation is used for it, si = 2 |Ki |/π,
(8)
where |Ki | means the area of the i-th element Ki . The idea of such formula is to use the diameter of a circle which has the same area to Ki as the size of the element, which is simple to calculate and reasonable for the mesh with regular element. With these parameters, the suggested time step δti for the element Ki is given by δti = σcf l si /(|u| + c),
(9)
where σcf l is the CFL number which is taken as 0.8 in this paper. Then the global time step δt is chosen as δt = min{δti }. i
(10)
3.2. Solution Reconstruction In the framework of the Godunov finite volume schemes, the available quantity from each governing cell is the cell average of each conservative variable. To deliver a high resolution method, solution variation in each element is needed, which can be obtained from the solution reconstruction. In this paper, we focus on the linear reconstruction of the conserved quantities in each element. Two components are needed for the solution reconstruction, i.e., the selection of the reconstruction patch for each element, and the reconstruction method. Below we describe the selection of the reconstruction patch first, then the reconstruction method follows. In the previous papers [20, 21], the reconstruction patch for the linear solution reconstruction has been discussed. Especially in [20], the influence of the different reconstruction patch for the h-adaptive methods is discussed and tested. In this work, we use two different reconstruction patches shown in Fig. 1 to reconstruct the solutions. For the first kind of the reconstruction patch for the current element K0 , the element which has a common edge with K0 will be selected as one element in the patch. With this principle, the patches P1 and P1B will be used as the reconstruction patches for the internal element and the element located on the domain boundary, 6
Figure 1: Left: the patch P1 = {Ki , i = 0, 1, 2, 3}, and the patch P2 = {Ki , i = 0, 1, . . . , 13} for an internal element K0 . Right: the patch P1B = {Ki , i = 0, 1, 2}, and the patch P2B = {Ki , i = 0, 1, . . . , 7} for an element K0 located on the boundary of the domain.
respectively. The drawback for P1B is obvious. First, there are only three elements (or even two in some extreme situation) in the patch, which would result in a determined (or under-determined in the extreme situation) system. Even worse, the distribution of all elements in the patch is almost along the domain boundary, which would cause the determined system ill-posed. The negative influence by using the patches P1 and P1B has been observed in [20] for steady Euler equations. Similar unsatisfied simulations are also obtained in this paper for reactive Euler equations. Consequently, in the simulation with uniform mesh, the patches P1 and P2B are employed. It is noted that in P2B , the patch is effectively enlarged, and the element along the normal direction w.r.t. the domain boundary is also introduced which effectively improves the condition number of the reconstruction system. Hence, numerical convergence of the simulation is observed successfully in the numerical experiments in the last section. When h-adaptive method is used, the mesh becomes nonuniform. Then the sizes of the neighbor elements for K0 could be very different. In [20], it has been shown that the use of P1 and P2B would bring negative effect on solving steady Euler equations, and to improve the robustness of the algorithm, the patches P2 and P2B are used for h-adaptive simulations. With this reconstruction patch, the smooth convergence of the steady state simulations is obtained. In this paper, it is found that the use of the patches P1 and P2B can handle the simulations with h-adaptive simulations well. Hence in the last section, we still use the patches P1 and P2B for the solution reconstruction. For simplicity of the description, we use Pi to denote the reconstruction patch for the element Ki in the following description. Once the patch is available, it is ready to do the solution reconstruction in each element. In this paper, the linear reconstruction is considered, and the least square method is used. More
7
specifically, in the element Ki , a linear variation of the conservative variable is supposed, i.e., ¯i + a1 (x − xi ) + a2 (y − yi ), P (x, y) = U
(11)
¯i is the cell average of the conservative quantity on the element Ki , and (xi , yi ) is the barycenter of where U Ki . It is noted that for a linear function, its evaluation on the barycenter of a given element is equal to its average on the element. Then with a given patch Pi = {Kj }, the gradient (a1,i , , a2,i )T in (11) is calculated by solving arg min a1,i ,a2,i
¯j − U ¯i − a1,i (xj − xi ) − a2,i (yj − yi )||22 . ||U
(12)
Kj ∈Pi ,j=i
Limiting process is necessary for the second order or higher order solution reconstruction to restrain the possible nonphysical oscillation around the shock. Due to its ability on delivering nonoscillatory and high order numerical solutions, WENO method has been widely used in a variety of practical simulations. In this paper, we follow [19, 21] to design the WENO process as follows. First, it is assumed that there are M (m)
(m)
candidate gradients {a1,i , a2,i , m = 1, 2, . . . , M } in the element Ki after the solution reconstruction. For example for the element K0 in Fig. 1 (left one), four candidate gradients can be obtained by solving (12) for K0 on the patches {K0 , K1 , K2 , K3 }, {K0 , K1 , K4 , K5 },{K0 , K2 , K8 , K9 }, and {K0 , K3 , K11 , K12 }, if the first kind of patch is used. Then for each candidate gradient, a smoothness indicator is defined as follows, (m) (m) (m) (a1,i + a2,i )dxdy, m = 1, 2, . . . , M. Si = Ki
(m) ω ˜i
(m) Si ) β ,
(m)
= 1/( + the weight assigned to each candidate gradient is designed as ωi = By introducing (m) (k) ˜ i . Here β and are two constants, and we use β = 2 and = 1.0e − 04 in the simulations. The ω ˜i / k ω final reconstructed gradient on the element Ki is given by T (m) (m) (m) (m) T ωi a1,i , ωi a2,i . (a1,i , a2,i ) = m
m
For the WENO reconstruction for the second kind of patches, i.e., P2 and P2B , it follows the above calculations. We refer to [22, 23] for the implementation detail. 3.3. h-Adaptive Methods In the simulation of the detonation wave, resolving the detonation front effectively and dynamically is crucial for a quality simulation. Since the variation of the conservative quantities in this region is relatively much larger than that in the other region, a uniform dense mesh in the whole domain is not appropriate from the implementation efficiency point of view. Consequently, mesh adaptive methods become the possible choice to resolve this issue. In this paper, we consider to use the local refinement technique to handle the adaptive process. This technique has been already adopted for the steady Euler equations, see [20, 22, 23, 28] and references 8
therein. For the h-adaptive technique, an efficient management strategy on handling the local refinement and coarsening of the mesh is crucial for the efficiency of the implementation. In the previous work, we use the hierarchy geometry tree (HGT) data structure, which is based on the fork-tree data structure in 2D case, to handle the mesh refinement and coarsening. Briefly speaking, the original mesh consists of a set of root elements. During the simulation, a series of the subtrees for those root elements would be generated according to the error indicators produced by the numerical solutions. By collecting the leaf nodes of all subtrees, a nonuniform mesh will be generated for the computational domain. With the HGT, the CPU time consumed on the local refinement and/or coarsening of the mesh is nothing compared with the whole CPU time. Please refer to [26] for detail of HGT, and also [4, 6] for 3D case. In the following, we will try to answer the other two questions on the quality realization of the mesh adaptive methods for the detonation wave simulation, i.e., where to refine the mesh, and how to update the conservative solutions. The first question mentioned above is related to the generation of the error indicator. During the adaptive refinement process, an error indicator will be assigned to each element. Compared with a given tolerance, it would be determined that whether an element will be refined into four subelements, and whether a set of subelements will be cast into a big element. In the detonation simulation, there is a shock for the pressure around the detonation front, and the peak pressure around the detonation front is also a very important quantity. Hence, it seems like that the gradient of the pressure would be a natural choice for the indicator generation. However, since the variation of the pressure around the detonation front is pretty dramatic, the gradient of the pressure would be very sensitive to any small numerical perturbation. Although the perturbation can be reduced systematically with reducing the mesh size, the quality of the numerical solution can still be affected since the existence of the transverse detonation instability for a two dimensional simulation. In the implementation of the algorithm, if the gradient of the pressure is employed to deliver an error indicator, some clustered locally refined mesh grids would be generated along the transverse direction of the detonation front, which might trigger the transverse instability numerically. Hence, to provide a uniform simulation environment is important, and a possible way to resolve the issue is to resort to some methods such as level set method [3] and phase field method [34] for the interface tracking, then the indicator for each element is generated with the knowledge of the interface. In this paper, we focus on the semi-2D and 2D simulations for the Zel‘dovich-Neumann-D¨oring (ZND) models. Hence, the structure of the shock is almost one dimensional with perturbations caused by the transverse instability of the simulations. Since the shock position can be predicted in advance [29, 30, 40], it is a reasonable choice to use skew normal distribuction function to generate a distribution of the error indicator along the x-direction. The skew normal distribution is given as f (x|α, μ, σ 2 ) = 2Φ(x|α)φ(x|μ, σ 2 ),
9
(13)
where φ(x|μ, σ 2 ) is the normal distribution φ(x|μ, σ 2 ) = √
1 2σ 2 π
exp (−
(x − μ)2 ), 2σ 2
and Φ(x|α) is cumulative distribution function given by 2 1 1 + erf ( √ ) . Φ(x|α) = 2 αx
(14)
(15)
In the above formulas, σ 2 and μ stand for the variance and the mean of the distribution, respectively, and α is a parameter to control the skew of the distribution. In the simulations, the parameter μ is chosen as the position of the leading shock which is given by μ = max x
Px , ρ
(16)
and the parameters σ 2 = 64, and α = 10. After the skew normal distribution is determined, the component of the error indicator for the element Ki from the normal distribution is given by 2 i indK nd = |Ki |f (xi |α, μ, σ ),
(17)
where |Ki | is the area of the element Ki , and the xi is the x-component of the centroid of the element Ki . There are two advantages on using the above refinement indicator. First, by choosing the position of the leading shock as μ, the indicator would reach the maximum value around the leading shock, which means that the mesh grids in this area would be much denser than that in other area. Second, the distribution function would reach zero slowly and smoothly, which means that the mesh grids would become coarse gently. This is positive for the simulation since solutions in the reaction zone can also be resolved well. The excellent results can be observed in the last section by using this kind of indicators. For the solution update from the old mesh to the new mesh, we need to take care of two things, i.e., conservation of the cell average and the high resolution of the numerical solution. With the HGT data structure, this process could be realized easily. Let To and Tn denote the old mesh and new mesh, respectively. From the property of the HGT data structure, there are following four relations between each two elements from two meshes, i.e., for Ko ∈ To and Kn ∈ Tn , we have the following relations, K o ≺ Kn ,
Ko K n ,
Ko = Kn , and Ko has nothing to do with Kn .
Here Ko ≺ Kn means Ko is a child of the element Kn , and Ko Kn means Ko is the parent element of Kn , and Ko = Kn means that two elements are the same node in the HGT. Suppose that the relation between two elements Ko and Kn is known, the solution update from Ko to Kn will be implemented as follows in this paper.
10
¯n + = U ¯o ∗ |Ko |/|Kn |; • Ko ≺ Kn : U ¯ n = Uo (X K ); • Ko K n : U n ¯n = U ¯o , • Ko = K n : U K means the barycenter of the element K. Uo means the where |K| means the area of the element K, and X reconstructed solution. It is noted that with HGT data structure, the relation between each two elements from two meshes can be obtained very efficiently, which results in an efficient implementation on the solution update. 3.4. OpenMP Parallelization The h-adaptive methods introduced in last subsection could effectively enhance the algorithm efficiency by controlling the amount of mesh girds in the simulation. To further improve the algorithm efficiency, the OpenMP parallelization technique will be considered in this subsection. OpenMP[1] is an API that supports shared memory multiprocessing programming in C/C++. Since the library for the numerical simulations in this paper is developed by C++, and the hardware is a DELL Precision T5610 workstation with 12 CPU Cores and 64 Gigabytes memory, it is a very natural choice to accelerate the implementation by parallelizing the algorithm with OpenMP. In the programming of the algorithm with OpenMP, an important issue is to resolve the data dependence among those threads. A straightforward way is to use the the clauses provided by OpenMP such as critical and atomic to avoid the data race. However, the efficiency of the implementation would be affected negatively, since in some sense the serial operations are forced with those clauses. Hence, a better way is to logically separate the write operations, then the calculation task can be assigned to every thread directly without data race. It is suggested that the edge/face-based data structure is preferred for the finite volume discretization from the efficiency point of view, since each edge in the mesh would be visited exactly once, instead of twice for an element-based data structure, in the numerical flux calculation and the cell average update. However, it is different to parallelize the above two operations. In the numerical flux calculation, there are two read operations to extract the left and right state of the conservative quantities from two adjacent elements, and only one write operation to save the calculated flux value on the current edge. Consequently, the for loop for the calculations of the numerical flux on the edges can be parallelized directly. The situation is totally different for the cell average update operation. In the update process, there is one read operation to extract the numerical flux from the current edge, and two write operations on two adjacent elements to accrue the change of the cell average. They are the two write operations which result in the data dependence among those threads, so the direct parallelization of the for loop for the cell average update is no longer acceptable. To resolve this issue, one may either generate a series of sets of independent edges first, then parallelize the 11
edge-based cell average update in each set, or implement an element-based cell average update instead. By element-based operation here, it means that we loop the element in the domain, then loop the edges in the element to calculate the new cell average. With such operation, except for the edge located on the domain boundary, each edge will be visited twice. However, the data dependence issue would be resolved effectively. In this paper, we use the element-based cell average update in the parallel simulation. For the other main components in the finite volume methods such as the solution reconstruction and the implementation of the Runge-Kutta module, either there is no data dependence, or the data dependence can be resolved trivially.
4. Numerical Tests The purpose of the numerical tests in this paper are two-fold. First, the numerical convergence of our method will be examined by a variety of benchmark problems, and the Zel‘dovich-Neumann-D¨ oring (ZND) model will serve the examination. Then the effectiveness of the proposed h-adaptive method will be demonstrated. 4.1. Initial Condition In this paper, we follow [15] to use the solutions of 1D ZND model as the initial conditions of our 2D simulations. The following is a brief description on how to generate the initial solutions for the 2D simulations with specific setup. Suppose that the distribution of the mass fraction f1 = ρr /ρ in the domain is known, and the x-component of the initial position of the detonation wave front is denoted by x0df . Then in the flame zone, the x-component velocity is given by γ u=− γ+1 where
ξ(f1 ) =
1 γ (D + ) γ+1 D
1 D+ D
2 −2
γ−1 γ+1
+
ξ(f1 ),
q0 (1 − f1 ) +
(18)
D2 γ + γ−1 2
.
Once u is obtained, the density and the pressure are calculated by the following formulas. ρ=−
D , u
P = Du + D2 + 1.
(19)
In the above two formulas (18) and (19), D is the detonation wave velocity, and is derived from the Chapman2 Jouguet velocity DCJ [13] and an overdriven factor f by 2 D2 = f DCJ ,
2 where DCJ = γ + (γ 2 − 1)q0 +
(20)
2
(γ + (γ 2 − 1)q0 ) − γ 2 .
In summary, once the distribution of the mass fraction f1 is known, the other quantities such as the x-component velocity u, density ρ and the pressure P can be deduced from (18), (19). Together with (24) 12
and the equation of state (3), all quantities in the reactive Euler equations (1) are obtained. Hence, the only thing left is how to find the distribution of f1 . In [41], the mass fraction f1 is generated by solving the following ODE, Kx0df f1 −Ea /T df1 =− e , dx u
(21)
ρ = 1, P = 1, T = 1, and u = −D,
(22)
with the boundary condition,
which denote the quantities in the unburned zone. Instead, in the simulations in this paper, the following function is adopted to simply approximate the distribution of the mass fraction f1 , ⎧ 1 ⎪ ⎪ tanh(π(x − x0df )) + 1 , x < x0df , ⎪ ⎨ 2 f1 = ⎪ ⎪ ⎪ ⎩ 1, elsewhere.
(23)
Although (23) is only a rough approximation for the solution of ODE (21), it can successfully give a jumpstart of the detonation wave simulation. We refer to [14] for a different f1 . For the simulations which need an initial perturbation along the y-direction, the y-component of the velocity is given by the following transversely sinusoidal perturbed planar ZND wave, i.e., ⎧ ⎪ ⎪ A sin (2k1 πy/L), x0df − Ww ≤ x ≤ x0df , ⎪ ⎨ v= ⎪ ⎪ ⎪ ⎩ 0, elsewhere,
(24)
where A, k1 , and L denote the amplitude, wave number and the wavelength of the perturbation, and Ww is the width of the perturbed zone along the x-axis. 4.2. Boundary conditions The computational domain for the following numerical simulations is a rectangle channel. The setup of the flow of the fluid is as follows. The reactant flows into the domain from the right boundary, and the production leaves the domain from the left boundary. The slip wall is employed for the lower and upper boundaries of the domain. Under the assumption of the ZND model, the velocity of the fluid in the domain is shifted by the negative Chapman-Jouguet velocity. Hence, the subsonic inflow boundary condition and subsonic outflow boundary condition are adopted on the right and left boundaries of the domain, respectively. The reflective boundary condition is adopted for both lower and upper boundaries of the domain to simulate the solid wall. To impose the subsonic inflow boundary, the classical way of using Riemann invariant can be used directly. However, since the configuration of the fluid in the far field outside the outflow boundary is not 13
clear, it brings difficulty to set up the boundary condition when the physical boundary condition is needed. To resolve this issue, several strategies are developed to make the numerical simulations reliable. In [30], an extrapolation technique is employed, where the conserved quantities of the system are relaxed to the condition at infinity with a given formula. In this case, the subsonic outflow boundary condition can be imposed trivially. In [15], the perfectly matched layer (PML) method is introduced in the simulation, and excellent numerical results are obtained. Besides the above two strategies, a straightforward way to handle the above issue is to use a sufficiently long computational domain along the x-axis, then the boundary values for those conservative quantities are obtained by the extrapolation. A formula for the relation between the length of the domain and the time can be found in [24, 29]. In the numerical simulations in this paper, we will try to use sufficiently long domain to restrain the nonphysical boundary effect. 4.3. Semi-2D ZND simulations The purpose of the numerical examples presented in this subsection is to verify the validity of our algorithm. The parameters used in the simulations are given in Table 1. To unify the length of the halfreaction zone L1/2 , the pre-exponential factor K is chosen as 145.69 and 230.75 for stable and unstable cases, respectively. γ
E0
q0
f
Stable
1.2
50.0
50.0
1.8
Unstable
1.2
50.0
50.0
1.6
Table 1: Parameters used in the semi-2D ZND simulations.
The parameters in Table 1 are used in previous works [13, 14, 29] to test the one dimensional ZND model. In this paper, we focus on numerical methods for two dimensional reactive Euler equations. Hence, a semi2D configuration for the numerical experiments are adopted in the simulations, i.e., the velocity component v in the equations is initialized as 0 in the whole domain, while the other quantities are initialized by the analytical solutions for the one dimensional ZND model. With this configuration, theoretically, the numerical solutions obtained from the semi-2D problems should be consistent with the solutions from one dimensional problems. However, due to its intrinsic instability for the multi-dimensional detonation wave, it is nontrivial to reproduce 1D results. From the stability analysis presented in [36], it is known that only a discrete set of the transverse disturbance whose wavelengths are compatible with the channel width can exist when the solid wall boundary condition is adopted, and the wavelengths have the following relations with the width of the channel, wi = 2W/i, i = 1, 2, 3, . . . ,
14
where wi and W stand for the wavelengths of the transverse disturbance and the channel width, respectively. From the above formula, it can be seen that the largest wavelength is w1 = 2W . Hence, if the number of the mesh grids along the y-direction is not big enough, the disturbance introduced by the discretization error along the y-direction will not be resolved well, and our semi-2D simulations will reproduce the numerical results given by 1D ZND model. Based on the above discussion, the computational domain used in the following two numerical experiments is [0, 130] × [−0.0433, 0.0433]. The initial mesh for this domain is given in Fig. 2. The mesh size is around 0.1. To study the numerical convergence, the simulations need to be implemented on a successively refined meshes. To remove the transverse instability when the mesh grids along y-direction become dense, the uniform refinement of the mesh grids is abandoned. Instead, the total length of the domain along x-direction and the regularity of the triangle element is kept, only the mesh size becomes 0.05. In this case, the mesh is similar to the coarse one, and the domain now becomes [0, 130] × [−0.0217, 0.0217]. The successively refined meshes in this subsection are generated by this strategy. In the following simulations, the initial position of the detonation front is xdf = 80.
Figure 2: Three fractions of the initial mesh used in Subsection 4.3.
For the simulations with the stable parameters given in Table 1, the results obtained on three successively refined meshes are given in Fig. 3. Since the parameter K = 145.69 makes L1/2 the length of the half-reaction zone the unit, there are around 10 mesh grids in the zone along the x-direction for the initial mesh, and 20 and 40 for the two refined meshes. In all three simulations, the expected stable state of the detonation is obtained smoothly. The numerical convergence can also be observed from the results successfully, i.e., the front pressure approaches the reference value 75.79 [8] by refining the mesh. 85
80 75 70
85
Ref. mesh size 0.05
Front Pressure
mesh size 0.1 Ref.
Front Pressure
Front Pressure
85
80 75 70
0
20
40
60 Time
80
100
Ref. mesh size 0.025
80 75 70
0
20
40
60 Time
80
100
0
20
40
60
80
100
Time
Figure 3: Mesh refinement study of the pressure history at the shock front for the stable detonation wave simulations. Reference value is 75.79[8].
Two remarks are given below for the results shown in Fig. 3. First, as we can observe from the results, there is a large oscillation of the pressure at the beginning of the simulations. This can be explained by the 15
use of the initial distribution of f1 given by (23), which is actually only an approximation for the 1D ZND model. This inconsistency introduces the large oscillation. Second, the oscillation in the numerical solution disappears slower in the dense mesh case than that in the coarse mesh case. This can be partially explained as follows. With denser mesh, the inconsistency between the exact ZND solution and the initial function can be resolved better. Since there is no instability mode for the simulation with stable parameters, all perturbations introduced by the initial condition go to zero eventually. In Fig. 4, the numerical results obtained from four successively refined meshes with the unstable parameters given in Table 1 are shown. The initial mesh is given by Fig. 2 with the mesh size 0.1, and the mesh sizes for three successively refined meshes are 0.05, 0.025, and 0.0125, respectively. With the pre-exponential factor K = 230.75, the length of the half-reaction zone is unit, which means that there would be 10, 20, 40, and 80 mesh grids in the region for the four meshes, respectively. Corresponding to these parameters, it is a 1D unstable detonation, and the reference top pressure is 101.1 ± 0.2[11, 29]. From Fig. 4 (top left), it can be seen that the simulation fails to present pulsation structure of the front pressure. Actually, the solution becomes stable with the time evolution. This can be explained that the mesh is too coarse to preserve the instability mode of the detonation wave. From the following three simulations implemented on three successively refined meshes, the periodic unstable structures are successfully observed, and the peak pressure approaches the reference value with the refinement of the mesh grids, which shows the numerical convergence of the method for this case.
Figure 4: Mesh refinement study of the pressure history at the shock front for the stable detonation wave with the parameters in the second row of Table 1.
4.4. Cell structure simulation The numerical test in this subsection is a two dimensional case, and the simulation will show that the regular cell structure would be generated after the detonation front, with the parameters given in Table 2. In the simulations, xdf = 62, and K = 1134363.64. To initially perturb the velocity v along the y-axis, the formula (24) with the parameters A = 0.1, L = 10, k1 = 4, and Ww = 1 are used. 16
γ
E0
q0
f
1.2
20.0
2.0
1.1
Table 2: Parameters used in the cell structure simulation.
First, the following three computational domains are tested to check the boundary effect on the simulations, i.e., [0, 34.64] × [−5, 5], [0, 69.28] × [−5, 5], and [0, 103.92] × [−5, 5]. Fig. 5 shows the initial mesh on the domain [0, 69.28] × [−5, 5].
Figure 5: Computational domain [0, 69.28] × [−5, 5].
In all three initial meshes, the side length of the equilateral triangle are all equal to 11.55. In the first test, each initial mesh are uniformly refined 7 times, i.e., the side length becomes 0.09. Then the simulation on each mesh runs till t ≈ 200. The record of front pressure of each simulation is presented in Fig. 6. From the figure, it can be observed that for the domain with length 34.64, the peak pressure keeps decaying from around t = 120 (upper one), while the periodic pulsation structure is preserved well for the domains with length 69.28 and 103.92, respectively (lower one). Furthermore, it is found that there is almost no difference between the records of the front pressure from two domains (Fig. 6, lower one). Hence, in the following simulations, the domain [0, 69.28] × [−5, 5] and the initial mesh shown in Fig. 5 are used. Remark 1. Since its intrinsic instability for the multi-dimensional detonation wave, small perturbation along the transverse direction of the leading shock may cause the different behavior of the numerical solution. To remove the numerical perturbation introduced from the solution reconstruction due to the nonuniform distribution of the mesh grids, the equilateral triangle meshes are employed in the simulations. To demonstrate the influence of the mesh grids on the numerical solutions, the simulations with the parameters given in Table 2 are reimplemented in the domain [0, 34.64] × [−5, 5] with two difference initial meshes shown in Fig. 7. The mesh in Fig. 7 (left one) is uniformly refined 7 times, while it is 5 times for the right mesh in the figure so that the mesh sizes for both meshes are almost the same. The results on the front pressure are shown in Fig. 8, in which it can be observed that formation of the periodic pulsation structure with the mesh generated by EasyMesh[32] is much earlier than that generated by the equilateral triangle mesh. It can also be observed that with the refinement of the mesh, the formation of the pulsation structure is delayed, but still earlier than that for equilateral triangle mesh case. This obviously shows the sensitivity of the numerical solution to the perturbation. In the following simulations, the mesh shown in Fig. 5 is used as the initial mesh. In the first simulation, the mesh is generated from the initial mesh by uniformly refining 6 times, the side length of the triangle is around 0.18. Then the simulations are reimplemented on three successively refined meshes. The records of 17
Font Pressure Font Pressure
6.5 5.5 4.5 3.5 2.5
Domain [-5,5] X [0, 34.64]
6.5 5.5 4.5 3.5
Domain [-5,5] X [0, 69.28] Domain [-5,5] X [0, 103.92]
2.5 0
20
40
60
80
100
120
140
160
180
200
Time Figure 6: The records of the front pressure obtained from simulations on three meshes with different length.
Figure 7: Two different meshes in the domain [0, 34.64] × [−5, 5]. Left one is an equilateral triangle mesh, and right one is generated by EasyMesh[32].
.
Front Pressure
7 Equilateral with 57921 points EasyMesh with 41121 points EasyMesh with 197569 points
6 5 4 3 0
20
40
60
80
100
Time Figure 8: The records of the front pressure obtained from equilateral triangle mesh and meshes generated from EasyMesh.
18
the front pressure obtained from first three simulations are shown in Fig. 9, and the pressure distribution and the vorticity for each simulation at t ≈ 30 are shown in Fig. 10. From Fig. 9, the intrinsic instability of the multi-dimensional detonation is shown obviously. Although with the equilateral triangle mesh, the stable state of the detonation can last a long time, and with the refinement of the mesh, it can last longer, the periodic pulsation structure of the front pressure eventually appears. Fig. 10 shows the pressure and the vorticity at the time t ≈ 30 before the formation of the pulsation structure. It can be seen that the regular cell structure is generated successfully, and with the denser mesh, the resolution of the structure becomes better. In Fig. 11, the pressure and vorticity at t ≈ 130 are shown. In this time, the periodic pulsation structure has been formed already.
Front Pressure
7
27105 mesh grids 107457 mesh grids 427905 mesh grids
6 5 4 3 0
20
40
60 Time
80
100
120
Figure 9: The records of the front pressure obtained from simulations on three meshes with different sizes.
4.5. Simulations with h-adaptivity In this subsection, the adaptive module of the proposed method is tested. For the convenience on the comparison, the parameters in Table 2 are still used in this subsection, and the mesh shown in Fig. 5 is used as the initial mesh. The other parameters used in the simulation are same to that in the last subsection. The results obtained from the simulation are shown in Figs. 12 and 13. In Fig. 12, the record of the front pressure from adaptive method is shown. As a comparison, the records obtained from two successively refined meshes are also shown in the figure. It can be seen that in the adaptive method, there are 85254 mesh grids in the domain, and the record of the front pressure from the adaptive method almost duplicate that from the uniform mesh which has 427950 mesh grids (subfigure in Fig. 12). This is reverified by the pressure and vorticity in Fig. 13 (top left and top middle, receptively) that they coincide with the ones shown in Fig. 10 (third row) very well. In this case, over 80% mesh grids are saved with the help of the adaptive method to reach almost the same accuracy. Fig. 13 (bottom one) shows the distribution of the mesh grids in the whole domain, while the Fig. 13 (top right one) shows the distribution of the mesh grids around the leading shock region. It can be observed that with the indicators generated by (17), most mesh grids are 19
Figure 10: The pressure (left column) and the vorticity (right column) obtained from simulations on four successively refined meshes at t ≈ 30.
20
Figure 11: The pressure (left column) and the vorticity (right column) obtained from simulations on the mesh with 427905 grid points at t ≈ 130.
Front Pressure
7
Adaptive mesh with 85254 points Equilateral with 107457 points Equilateral with 427905 points
3.8 3.7 3.6
6 5
15 17 19 21 23 25
4 3 0
10
20
30
40 Time
50
60
70
80
Figure 12: The record of the front pressure obtained from the adaptive method, and the comparison to the records obtained from two successively refined equilateral triangle meshes. The subfigure shows the detail of the comparison.
21
Figure 13: The pressure (top left), the vorticity (top middle), the local mesh (top right) around the leading shock, and the mesh (bottom) in the whole domain, obtained at t ≈ 30.
located around the leading shock region and the reaction region, which guarantees the numerical accuracy of the simulations. There are only a small amount of the mesh grids in other regions, which guarantees the efficiency on the simulation effectively.
Conclusion In this paper, an adaptive finite volume method is proposed for the reactive Euler equations. The numerical discretization includes a second order TVD Runge-Kutta method for the temporal discretization, and a finite volume method with linear solution reconstruction for the spatial discretization. To restrain the potential numerical oscillation introduced by the high order solution reconstruction, the WENO methodology is employed as the limiting strategy. To enhance the efficiency of the algorithm, the h-adaptive methods are introduced in this work. The indicator and the solution update method are developed from the consideration of the efficiency and the numerical accuracy. The OpenMP technique is also introduced in the algorithm to further improve the efficiency. A variety of the numerical tests are implemented to show the numerical convergence of the proposed method. The effectiveness of the h-adaptive method is also successfully demonstrated. The unstructured meshes have advantage for the practical simulations when the geometry of the domain is complicated. However, the numerical results in this paper show that the numerical behavior of the detonation simulation is sensitive to the perturbation introduced by the solution reconstruction on the unstructured meshes. Hence, a quality solution reconstruction method for the unstructured mesh grids is under considering. The other issue is related to the error indicator for the mesh adaptation. Although
22
the strategy used in this paper works well, a fine error indicator is needed to handle more complicated simulations, and to further enhance the efficiency. The related results will be reported in the forthcoming paper. Acknowledgements The work of G. H. Hu was partially supported by 050/2014/A1 from FDCT of Macao S.A.R., MYRG201400109-FST from University of Macau, and National Natural Science Foundation of China (Grant No. 11401608). The author would like to thank Prof. W. S. Don, Prof. Z. Gao, Mr. P. Li, and Prof. R. Li for their helpful suggestions and discussion on this work. [1] http://openmp.org. [2] Marco Arienti. A numerical and analytical study of detonation diffraction. PhD thesis, California Institute of Technology, 2003. [3] Tariq D. Aslam, John B. Bdzil, and D.Scott Stewart. Level set methods applied to modeling detonation shock dynamics. Journal of Computational Physics, 126(2):390 – 409, 1996. [4] Gang Bao, Guanghui Hu, and Di Liu. An h-adaptive finite element solver for the calculations of the electronic structures. Journal of Computational Physics, 231(14):4967 – 4979, 2012. [5] Gang Bao, Guanghui Hu, and Di Liu. Numerical solution of the Kohn-Sham equation by finite element methods with an adaptive mesh redistribution technique. Journal of Scientific Computing, 55(2):372– 391, 2013. [6] Gang Bao, Guanghui Hu, and Di Liu. Real-time adaptive finite element solution of time-dependent Kohn-Sham equation. Journal of Computational Physics, 281:743–758, 2015. [7] Gang Bao, Guanghui Hu, and Di Liu. Towards translational invariance of total energy with finite element methods for Kohn-Sham equation. Communications in Computational Physics, 19:1–23, 2016. [8] Ralf Deiterding. Parallel adaptive simulation of multi-dimensional detonation structures. September 2003. [9] Ralf Deiterding. A parallel adaptive method for simulating shock-induced combustion with detailed chemical kinetics in complex domains. Computers & Structures, 87(11–12):769 – 783, 2009. Fifth MIT Conference on Computational Fluid and Solid Mechanics. [10] Hua-Shu Dou, Her Mann Tsai, Boo Cheong Khoo, and Jianxian Qiu. Three-dimensional numerical simulation for detonation waves using WENO schemes. American Institute of Aeronautics and Astronautics, 45th AIAA Aerospace Sciences Meeting and Exhibit, 2007. 23
[11] Jerome J. Erpenbeck. Stability of idealized one-reaction detonations. Physics of Fluids, 7(5):684–696, 1964. [12] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). Journal of Computational Physics, 152(2):457 – 492, 1999. [13] Wildon Fickett and Willam C. Davis. Detonation: Theory and Experiment. Dover Publications, 2011. [14] Zhen Gao, Wai Sun Don, and Zhiqiu Li. High order weighted essentially non-oscillation schemes for onedimensional detonation wave simulations. Journal of Computational Mathematics, 29:623–638, 2011. [15] Zhen Gao, WaiSun Don, and Zhiqiu Li. High order weighted essentially non-oscillation schemes for two-dimensional detonation wave simulations. Journal of Scientific Computing, 53(1):80–101, 2012. [16] Sigal Gottlieb and Chi-Wang Shu. Total variation diminishing Runge-Kutta schemes. Mathematics of Computation, 67:73–85, 1998. [17] P. He and H.Z. Tang. An adaptive moving mesh method for two-dimensional relativisitic hydrodynamics. Communications in Computational Physics, 11:114–146, 2012. [18] William D Henshaw and Donald W Schwendeman. An adaptive numerical scheme for high-speed reactive flow on overlapping grids. Journal of Computational Physics, 191(2):420 – 447, 2003. [19] Changqing Hu and Chi-Wang Shu. Weighted essentially non-oscillatory schemes on triangular meshes. Journal of Computational Physics, 150:97–127, 1999. [20] Guanghui Hu. An adaptive finite volume method for 2D steady Euler equations with WENO reconstruction. Journal of Computational Physics, 252(0):591 – 605, 2013. [21] Guanghui Hu, Ruo Li, and Tao Tang. A robust WENO type finite volume solver for steady Euler equations on unstructured grids. Communications in Computational Physics, 9:627–648, 2011. [22] Guanghui Hu, Xucheng Meng, and Nianyu Yi. Adjoint-based an adaptive finite volume method for steady Euler equations with non-oscillatory k-exact reconstruction. Computers & Fluids, 139:174 – 183, 2016. [23] Guanghui Hu and Nianyu Yi. An adaptive finite volume solver for steady Euler equations with nonoscillatory k-exact reconstruction. Journal of Computational Physics, 312:235 – 251, 2016. [24] P Hwang, R P Fedkiw, B Merriman, T D Aslam, A R Karagozian, and S J Osher. Numerical resolution of pulsating detonation waves. Combustion Theory and Modelling, 4(3):217–240, 2000. 24
[25] Smadar Karni. Hybrid multifluid algorithms. SIAM Journal on Scientific Computing, 17(5):1019–1039, 1996. [26] Ruo Li. On multi-mesh h-adaptive methods. Journal of Scientific Computing, 24(3):321–341, 2005. [27] Ruo Li and Tao Tang. Moving mesh discontinuous Galerkin method for hyperbolic conservation laws. Journal of Scientific Computing, 27(1-3):347–363, 2006. [28] Ruo Li, Xin Wang, and Weibo Zhao. A multigrid block LU-SGS algorithm for Euler equations on unstructured grids. Numerical Mathematics: Theory, Methods and Applications, 1:92–112, 2008. [29] Y.S. Lian and K. Xu. A gas-kinetic scheme for multimaterial flows and its application in chemical reactions. Journal of Computational Physics, 163(2):349 – 375, 2000. [30] Z. Liang and L. Bauwens. Cell structure and stability of detonations with a pressure-dependent chainbranching reaction rate model. Combustion Theory and Modelling, 9(1):93–112, 2005. [31] E. Loth, S. Sivier, and J. Baum. Adaptive unstructured finite element method for two-dimensional detonation simulations. Shock Waves, 8(1):47–53, 1998. [32] Bojan Niceno. EasyMesh. http://web.mit.edu/easymesh_v1.4/www/easymesh.html. [33] J.T. Oden and L. Demkowicz. h-p adaptive finite element methods in computational fluid dynamics. Computer Methods in Applied Mechanics and Engineering, 89(1–3):11 – 40, 1991. Second World Congress on Computational Mechanics. [34] Y. Sun and C. Beckermann. Sharp interface tracking using the phase-field equation. Journal of Computational Physics, 220:626–653, 2007. [35] Tao Tang. Moving mesh methods for comutational fluid dynamics. Contemporary Mathematics, American Mathematical Society, 383:141–173, 2005. [36] B. D. Taylor, A. R. Kasimov, and D. S. Stewart. Mode selection in weakly unstable two-dimensional detonations. Combustion Theory and Modelling, 13(6):973–992, 2009. [37] Fumiya Togashi, Rainald L¨ ohner, and Nobuyuki Tsuboi. Numerical simulation of h2/air detonation using unstructured mesh. Shock Waves, 19(2):151–162, 2009. [38] Eleuterio F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer Berlin Heidelberg, 2009.
25
[39] Cheng Wang, Xiangxiong Zhang, Chi-Wang Shu, and Jianguo Ning. Robust high order discontinuous Galerkin schemes for two-dimensional gaseous detonations. J. Comput. Phys., 231(2):653–665, January 2012. [40] Daniel N. Williams, Luc Bauwens, and Elaine S. Oran. Detailed structure and propagation of threedimensional detonations. Symposium (International) on Combustion, 26(2):2991–2998, 1996. [41] Zeng-Chan Zhang, S.T. John Yu, Hao He, and Chang Sin-Chung. Direct calculations of two- and threedimensional detonations by an extended CE/SE method. 39th Aerospace Sciences Meeting and Exhibit, American Institute of Aeronautics and Astronautics, 2001.
26