Adaptive Runge_CKutta discontinuous Galerkin method for complex ...

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS Int. J. Numer. Meth. Fluids 2013; 73:847–868 Published online 9 July 2013 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/fld.3825

Adaptive Runge–Kutta discontinuous Galerkin method for complex geometry problems on Cartesian grid Jianming Liu1,2,3 , Jianxian Qiu1, * ,† , Ou Hu4 , Ning Zhao4 , Mikhail Goman3 and Xinkai Li3 1 School

of Mathematical Sciences, Xiamen University, Xiamen, China of Mathematics, Jiangsu Normal University, Xuzhou, China 3 Faculty of Technology, De Montfort University, Leicester, England 4 College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, China 2 Department

SUMMARY A Cartesian grid method using immersed boundary technique to simulate the impact of body in fluid has become an important research topic in computational fluid dynamics because of its simplification, automation of grid generation, and accuracy of results. In the frame of Cartesian grid, one often uses finite volume method with second order accuracy or finite difference method. In this paper, an h-adaptive Runge–Kutta discontinuous Galerkin (RKDG) method on Cartesian grid with ghost cell immersed boundary method for arbitrarily complex geometries is developed. A ghost cell immersed boundary treatment with the modification of normal velocity is presented. The method is validated versus well documented test problems involving both steady and unsteady compressible flows through complex bodies over a wide range of Mach numbers. The numerical results show that the present boundary treatment to some extent reduces the error of entropy and demonstrate the efficiency, robustness, and versatility of the proposed approach. Copyright © 2013 John Wiley & Sons, Ltd. Received 4 February 2013; Revised 29 May 2013; Accepted 8 June 2013 KEY WORDS:

RKDG; adaptive Cartesian grid; immersed boundary method; complex geometry; limiter; shock

1. INTRODUCTION In this paper, we investigate the h-adaptive Runge–Kutta discontinuous Galerkin (RKDG) method on Cartesian grid with ghost cell immersed boundary method and its applications for complex geometry. Recently, Cartesian grid methods have become popular in computational fluid dynamics (see [1–11] and their references), because such methods do not suffer from the complex grid generation and grid management requirements inherent to other methods and can be easily used in high order numerical schemes. Conceptually, the implementation of the Cartesian grid approach is much simpler than other grids. In this method, solid bodies are cut out of a single static background mesh, and their boundaries are represented by different types of cut cells, or solid bodies are equipped with ghost cells using the immersed boundary. When cut cells become very small, degenerate cells will be encountered. In this situation, numerical instability may occur when an explicit time step scheme is used in calculations. Some techniques have been employed to overcome those problems along with time step stability restrictions [1, 6–9]. Although there are many different techniques, ghost cell method (GCM) or immersed boundary method still obtains people’s favor because of its simplicity [2–5, 10, 11].

*Correspondence to: Jianxian Qiu, School of Mathematical Sciences, Xiamen University, Xiamen, China. † E-mail: [email protected] Copyright © 2013 John Wiley & Sons, Ltd.

848

J. LIU ET AL.

In addition to the simplification and automation of grid generation, Cartesian grid method also facilitates the adaptive mesh refinement (AMR). It is well known that for hyperbolic conservation law, most of AMR methods on Cartesian grid have only second order accuracy. Recently, there have been many research efforts on the high order accuracy Cartesian grid AMR method, in which essentially non-oscillatory or weighted essentially non-oscillatory (ENO/WENO) finite difference method is used [12–14]. However, there are limitations of the enhancement of the numerical accuracy when the interpolation schemes are implemented on the newly generated fine grids and coarsened child grids and the boundary conditions on the interface of coarse–fine grid [13]. According to the description of Shen et al. in [14], it is difficult to maintain high order accuracy across several levels of grids, and introducing a high order accurate scheme in the AMR setting makes it harder to maintain local mass conservation. Furthermore, high order scheme such as ENO/WENO will use broader grid stencil, which increases the difficulty on dealing with boundary conditions, in particular, demand more ghost points to retain the approximation accuracy on computational boundaries or the interface of coarse–fine grid. Discontinuous Galerkin (DG) method is a very effective high order numerical scheme for hyperbolic conservation law, which is developed by Cockburn et al. in a series of papers [15–19]. This method consists of many advantages, including flexible hp-adaptive ability, certain superconvergence property, and high parallelism of the resulting algorithm using explicit Runge–Kutta method. RKDG is not like the ENO/WENO method, the stencil used in the DG method is more compact. This is because the solution representation in each cell is kept independent of the solutions in other cells, in which communication occurs only with adjacent cells sharing a common edge. This makes the DG method extremely flexible for solving problems on complicated fluid boundaries and processing with adaptive grid. As an alternative of high order numerical algorithm during AMR, DG can also treat nonconforming element with hang node. During the process of the grid refinement or grid coarsening, for the prolongation and restriction of grid, DG method can use L2 projection to maintain local conservation. In addition, at the interface between the coarse and fine grid, the cell solution polynomial of the DG finite element space provides ready-made interpolation one. Because of these advantages, the RKDG method is a good alternative for high order accuracy numerical method on the adaptive Cartesian grid. However, the research on h-adaptive RKDG method is limited compared with the finite difference Cartesian grid AMR method. We refer, for instance, to Flaherty et al. [20–23] and Dedner et al. [24]. We also refer to Hartmann and Houston [25], where duality techniques were used for designing adaptive strategy. Recently, Zhu and Qiu [26, 27] systematically study the h-adaptive RKDG method with WENO limiter. To the authors’ knowledge, there are lacking practical applications on dealing with complex geometry with RKDG method on the Cartesian grids, such as flow past an airfoil. In fact, the high sensitivity to the accuracy of the boundary representation of DG and the appearance of cut cells near the wall limit its applications [28, 29]. In this paper, a discontinuous detector and a limiting technique are used to suppress oscillation, and a solution-based grid adaptive method with hybrid sensors is used during the grid refinement. Furthermore, we will focus on the implementation detail of the adaptive Cartesian RKDG method for a complex body. The paper is structured as follows. In Section 2, the governing equations and their numerical discretization are provided, including the review of the RKDG method, a detailed description of the h-adaptive RKDG method, and the limiting process near the discontinuous. The implementation of the ghost cell immersed boundary method coupling with the high-order accurate RKDG method is presented in Section 3. In Section 4, the computational results for two-dimensional (2D) benchmark problems with complex geometries are presented, illustrating the high order accuracy, robustness, and efficiency of the proposed h-adaptive RKDG scheme. Finally, a conclusion is given in Section 5.

2. GOVERNING EQUATIONS AND THE NUMERICAL METHOD The inviscid compressible Euler equations can be given in vector form explicitly expressing the conservation laws of mass, momentum, and energy, which are Copyright © 2013 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Fluids 2013; 73:847–868 DOI: 10.1002/fld

ADAPTIVE RKDG METHOD FOR COMPLEX GEOMETRY ON CARTESIAN GRID

849

@U C r F D 0, @t

(1)

where the conservative state vector U and the inviscid flux vectors F D .f , g/T are defined by 2 2 3 3 3 2 v u 6 6 u2 C p 7 7 uv 6 u 7 7,g D 6 7 ,f D 6 U D4 4 v 2 C p 5 . 4 5 v 5 uv E u.E C p/ v.E C p/ The variables p, , u, and v are the pressure, the density, and the two Cartesian components of the velocity vector, respectively, and E represents the total energy per unit mass. The pressure p is obtained using an equation of state for ideal gases: 1 2 2 p D . 1/ E u C v . (2) 2 2.1. Review of RKDG method To formulate the RKDG method, we first discretize (1) in space using the DG method. For simplicity, U is considered as a scalar function here only. If U is vector-valued, one simply proceeds with component by component. Assuming that Th is a grid coverage of where the domain is subdivided into a collection of nonoverlapping elements K, we seek the approximate solution Uh .t , x, y/ in the finite element space of discontinuous functions defined as ° ± Vh D vh 2 L1 ./ W vh jK 2 P k , 8K 2 Th , where P k being the space of polynomials with degree 6 k. Firstly, Equation (1) is multiplied by a test function v.x, y/ 2 Vh and integrated over cell K, then the following semi-discrete form of (1) is obtained by using integration by parts: Z Z XZ d U.x, y, t /vdxdy F .U / rvdxdy C F .U / nE e vds D 0, (3) dt K k e e2@K

R where nE e is the outward unit normal to the edge e. The volume integral term K F .U /rvdxdy can be computed either exactly or by a numerical quadrature with sufficiently high order of accuracy. The line integral in (3) is typically discretized by a Gaussian quadrature with sufficient accuracy: Z q X F .U / nE e vds jej !l F .U.Gl , t //E ne v.Gl /, e

lD1

where q is the number of Gaussian quadrature points. Because the numerical solution Uh is discontinuous between each element on the interfaces, ne must be replaced by a the flux F .U.Gl , t //E numerical flux function FO U .Gl /, U C .Gl /, nE e , which depends on both the inner and outer trace of Uh on @K, K 2 Th , and the unit outward normal vector nE e of edge e. It should be emphasized that the choice of the numerical flux function is independent on the finite element space employed. In this study, for the Euler equations, the HLLC approximate Riemann solvers [30] is used as the numerical flux. To complete the definition of the RKDG method, the following third-order total variation diminishing Runge–Kutta time discretization [31] is used to solve the semi-discrete form (3), this is u.1/ D un C tL.un , t n /, u.2/ D 34 un C 14 u.1/ C 14 tL u.1/ , t n C t , unC1 D 13 un C 23 u.2/ C 23 tL u.2/ , t n C 12 t , Copyright © 2013 John Wiley & Sons, Ltd.

(4)


850

J. LIU ET AL.

where L denotes the spatial discretization. In this study on h-adaptive RKDG method, for unsteady problems, global time steps are used in Runge–Kutta method, which are equal to the smallest local time step calculated by local maximum characteristic speed and the local cell size at each time level. For a steady problem, in order to accelerate the convergence to a steady solution, a local time step is used. 2.2. h-Adaptive RKDG method in 2D In fact, for h-adaptive Cartesian grid RKDG method, apart from the appearance of hanging cell, there is no main difference from general RKDG method. Cartesian grid in conjunction with tree data structure is a natural choice for solution with adaptive grid. In this work, we use a generalized binary tree (e.g., quad-tree in 2D, octree in 3D, etc.) data structure. We do not allow any cells with side length ratio > 2 or < 0.5 to be neighbors. This restriction results in what is known as a balanced tree. It should be noted here that such a tree is different from the one presented in [27], where they use a non-balance tree. For the RKDG method, if a local basis of P k is chosen and denoted as vl.K/ .x, y/ for l D 0, , Qk in cell element K, then numerical solution Uh .x, y, t / 2 Vh can be expressed as Uh .x, y, t /jK D

Qk X

UK.l/ .t /vl.K/ .x, y/,

(5)

lD0

where UK.l/ .t /.l D 0, , Qk / are the degrees of freedom. Commonly, the local orthogonal basis functions are adopted during the implementation and numerical calculation of RKDG method, which are 1 1, 1 , 2 , 12 , 1 2 , 22 3 inside the rectangular element K D x 1 , x 1 y K

1 D

2

KC

1 , , 3

1,y 1 K KC 2 2

2

, and

x xK y yK , 2 D , xK =2 yK =2

where .xK , yK / is the center of rectangle K, xK and yK are the lengths of K’s sides in the direction of x and y, respectively. Obviously, UK.0/ .t / is a cell average of Uh over K. In order to retain the approximation accuracy and the local conservation property, an L2 projection approach is used to obtain the new generated grid’s degrees of freedom [27]. Suppose Uh is already known on the mesh Th .tn /, and we need to determine the degrees of freedom UK.l/0 .tn /.l D 0, : : : , Qk / in the new cell K 0 2 Th .tnC1 /. Let Uh denote the L2 projection of Uh and satisfy the following equation: Z Z 0 .K 0 / 0 Uh jK 0 vl .x, y/dxdy D Uh vl.K / .x, y/dxdy. (6) K0

K0

Substituting (5) in (6), we have the degrees of freedom in the new generated grid K 0 . For simplicity, the implementation of the adaptive mesh for RKDG method is only based on P 2 polynomials on each element, that is, k D 2; however, if k > 2, the basic idea is the same as k D 2.

Figure 1. The sketches of coarsening (left) and refinement (right) in the adaptive mesh. Copyright © 2013 John Wiley & Sons, Ltd.



851

When four cells K1 , K2 , K3 , and K4 are flagged and coarsened to a new cell K 0 (see the left sketch in Figure 1), the new degrees of freedom are given by .0/ .0/ .0/ C U C U C U UK.0/0 D 14 UK.0/ , K K K 2 3 4 1 .1/ .1/ .1/ .1/ .1/ .0/ .0/ .0/ C U U C U UK 0 D 18 UK1 C UK2 C UK3 C UK4 C 38 UK.0/ K2 K3 K4 , 1 .2/ .2/ .2/ .0/ .0/ .0/ .0/ 3 UK.2/0 D 18 UK.2/ C U C U C U U C U C U C U K2 K3 K4 K K2 K3 K4 , 8 1 1 .3/ .3/ .3/ .2/ .2/ .2/ .2/ 1 3 (7) UK.3/0 D 16 C U C U C U C U U C U UK.3/ C U K2 K3 K4 K2 K3 K4 1 16 K1 3 9 C 16 UK.1/ C UK.1/ C UK.1/ UK.0/ UK.0/ C UK.0/ UK.1/ C 16 UK.0/ , 1 2 3 1 2 3 4 4 .4/ .4/ .4/ .4/ .4/ .1/ .1/ .1/ .1/ 1 15 UK 0 D 16 UK1 C UK2 C UK3 C UK4 C 32 UK1 C UK2 UK3 C UK4 , .5/ .5/ .5/ .2/ .2/ .2/ .2/ 1 15 UK.5/0 D 16 C U C U C U U C U C U UK.5/ C U K2 K3 K4 K1 K2 K3 K4 . 32 1 When a parent cell K is flagged and refined to four child cells K1 , K2 , K3 , K4 , (see the right sketch in Figure 1), the new degrees of freedom can be defined as UK.0/0 D UK.0/ C ai UK.1/ C bi UK.2/ C ai bi UK.3/ , i

UK.1/0 D 12 UK.1/ C .1/i ai bi UK.3/ C ai UK.4/ , i

UK.2/0 D 12 UK.1/ C .1/i jai bi jUK.3/ C bi UK.5/ , i

UK.3/0 D 14 UK.3/ ,

(8)

i

UK.4/0 i UK.5/0 i

D 14 UK.4/ , D 14 UK.5/ ,

and 1 1 1 ai D .1/i , b1 D b2 D , b3 D b4 D .i D 1, 2, 3, 4/. 2 2 2 It should be noted that in the original AMR finite difference method [32], one must impose correction of conservation law to make the fluxes into the fine grid across a coarse/fine cell boundary, as shown in Figure 2, equivalent to the flux out of the coarse cell after every time step advancing. However, the RKDG method is born with h-adaptive ability and can be applied in the presence of hanging nodes. In our study, a list is used to record the cell faces, and a loop over all cell faces

Figure 2. The treatment of hanging node in the adaptive Cartesian grid. Copyright © 2013 John Wiley & Sons, Ltd.


852

J. LIU ET AL.

is implemented to calculate the face flux integrations. At the coarse/fine cell face (as an example, see Figure 2), only the fine cell boundary ab and bc are added into cell face list. During time step advancing, numerical fluxes of face ab and bc of cell B and C are calculated, respectively. For the coarse cell A face, the numerical flux integral is divided into two partitions ab and bc, which are calculated by the same integral formula as the fine cell boundary face. This treatment will automatically maintain the flux conservation on both coarse and fine cell interface boundaries. 2.3. Discontinuous detector and limiting If there are strong discontinuities in the approximate solution, oscillations may occur or even the method may break down. To enhance the stability of the method and eliminate possible spurious oscillations in the approximate solution, some kinds of limiters must be added during the time evolution. Unfortunately, most of limiters frequently identify regions near smooth extrema as requiring limiting, and this typically results in a reduction of the optimal high-order convergence rate. In this study, the discontinuity detector introduced by Krivodonova et al. in [33] is used to distinguish regions where solutions are smooth or discontinuous. With such a method, limiting would only be used near discontinuities; thus, high-order accuracy would be preserved in smooth regions. The discontinuity detector works in the following way. First, partition the cell boundary @K into two portions @K and @K C , where the flow is into .V nE < 0, nE is the normal vector to @K/ and out of (V nE > 0/ cell K, respectively. Then, the discontinuity detector is defined as ˇ ˇR ˇ .Uk UnbK /ds ˇ @K , (9) IK D .kC1/=2 h j@K j k UK k and if IK > 1, UK is discontinuous; otherwise, UK is smooth. Note that here, h is the radius of the circumscribed circle in element K, nbK is the neighbor of K on the side of @K , and k UK k is a maximum norm based on local solution maximum at integration points. Considering it is effective both for the shock and the contact discontinuity, in this study, we use entropy as the discontinuity detection variable for Euler equations [33]. A local slope limiting introduced in [17, 18] is used near the discontinuous, which is defined as Ux

1 ,yK 2

KC

UxK ,y

1 KC 2

D UK.0/ C UQ Kx , UxC

K

1 ,yK 2

D UK.0/ C UQ Ky , UyCK ,y

1 K 2

D UK.0/ UQQ Kx , (10)

D UK.0/ UQQ Ky ,

where UQ Kx D UQ Ky D

Q Pk lD0 Q Pk lD0

x Q Pk .l/ .K/ UK.l/ vl.K/ xKC 1 , yK , UQQ K D UK vl xK 1 , yK , 2

lD0

2

y Q Pk .l/ .K/ UK.l/ vl.K/ xK , yKC 1 , UQQ K D UK vl xK , yK 1 . 2

lD0

x y We modify UQ Kx , UQQ K , UQ Ky and UQQ K in (10) as .0/ .0/ UQ Kx.mod / D m N UQ Kx , UeK UK.0/ , UK.0/ UwK , x x.mod / .0/ .0/ UQQ K Dm N UQQ K , UeK UK.0/ , UK.0/ UwK , .0/ .0/ UQ Ky.mod / D m N UQ Ky , UnK UK.0/ , UK.0/ UsK , y.mod / y .0/ .0/ UQQ K Dm N UQQ K , UnK UK.0/ , UK.0/ UsK ,

2

(11)

.0/ .0/ .0/ .0/ , UsK , UwK , and UnK denote the cell averages of the east, south, west, and north where UeK neighbor cells, respectively. m N is a modification of minmod function [16] defined by ² ifj˛1 j 6 Mx 2 ˛1 , m.˛ N 1 , , ˛m / D (12) m.˛1 , , ˛m /, otherwise,

Copyright © 2013 John Wiley & Sons, Ltd.



853

for the first two equations in (11), and with a change of x to y in the last two equations in (11). The minmod function m is defined by ´ sminj˛i j, ifs D sign.˛1 / D D sign.˛m /, i m.˛1 , , ˛m / D 0, otherwise, and M in (12) is the TVBM limiter constant, which can be chosen suitably for different problems. For the original TVBM limiter, M is taken to prohibit the degeneracy of accuracy at the smooth extrema and is made the resulting RKDG scheme retain its optimal accuracy. Unfortunately, to the date, there is not a very effective method to choose the coefficient M [17, 18]. In this study, we first use a discontinuity detector introduced by Krivodonova et al. to distinguish regions where solutions are smooth or discontinuous. Then, we use the TVBM limiter at those cells that have been detected as discontinuous cells. So, we suggest to choose small M . In the present study, except the special instructions, M D 10 is used in all examples. For uniform grids, there are no trouble in the limiting procedure; however, for AMR grid, the appearances of hanging nodes would make some neighbors not real cells with same refinement levels (see Figure 2, for example, the west neighbor F of cell A and the east neighbor A1 of cell B/. The cell averages of those imaginary ones, in this study, are calculated by using both first equation of (7) and (8). If there is a system, in order to improve the effectiveness on controlling oscillations, the limiting should be carried out in the local characteristic directions; the details for the system can be found in [17]. 2.4. Solution-based mesh refinement on Cartesian grids Sensors are employed to detect and localize physical flow phenomena. Because the divergence of velocity is direction independent and very effective in locating shock including strong shock and weak shock, the curl of velocity is direction independent and very effective in finding shear and vortex [34], and entropy finds shock and contact discontinuity well. In this work, for different flows, we will use different combination of those sensors, which includes a combination of the divergence of velocity and entropy as follows: 3

3

3

d i D jr V jdi2 , ei D jrS jdi2 jrp a2 rjdi2 ,

(13)

and a hybrid sensors of the divergence of velocity and the curl of velocity as follows: 3

3

d i D jr V jdi2 , ci D jr V jdi2 , (14) p for i D 1, 2, : : : , Nc , where Nc is the total number of cells and di D K (K is the cell volume). The standard deviation of the divergence of velocity, the entropy, and the curl of velocity are computed, respectively, as v v v uP uP uP u Nc 2 u Nc 2 u Nc 2 u u u d i d i t i D1 t i D1 t i D1 d i d D , e D , c D . Nc Nc Nc For different hybrid sensors of (13) and (14), a cell is flagged for refinement or coarsening if one of two possible conditions hold: If either d i > !1 d or ei > !2 e (ci > !2 c /, the cell is flagged for refinement; If both d i < !3 d and ei < !4 e (ci < !4 c /, the cell is flagged for coarsening;

where !l .l D 1, 2, 3, 4/ are adjustable coefficients based on different problems. Taking into account the computational efficiency, in this paper !1 and !2 are chosen between 1.0 and 1.5, and !3 and !4 are chosen between 0.1 and 0.4. In general, the choice of !l , l D 1, 2, 3, 4, only changes the total number of computational grid, but the computational results are not sensitive to the choice of !l . In this work, except the special instructions, we will always use hybrid sensors of (13) as the grid adaptive indicator. Copyright © 2013 John Wiley & Sons, Ltd.


854

J. LIU ET AL.

3. BOUNDARY TREATMENT It has been shown that the DG method for flow problems is highly sensitive to the accuracy of the boundary representation [28]. In order to cure the large errors arising near the boundary and reduce the pollution of the solution inside the domain, one often needs to use curved elements, which largely enhance the difficulties of geometry processing. Furthermore, in practical simulations, a curved geometry is often approximated by a straight-sided polygon. This produces larger numerical boundary error near the body boundary surface. This error may dominate the discretization error of the scheme and lead to a wrong solution. In this study, the body is immersed into the Cartesian volume grid, in other words, the quality of grid is not so good near the body boundary. Therefore, care should be taken when we deal with a complex geometry. In this paper, we adopt the idea of ghost cell immersed boundary method to treat nonconforming boundary for the RKDG method. Recently, Dadone et al. [4] use finite difference method to present some systemic results about a novel GCM for static body on Cartesian grid. In this paper, we extend it to treat the boundary in RKDG method on Cartesian grid. The basic idea is that the ghost cell values are developed on the basis of an assumption of flow field model in the vicinity of the wall consisting of a vortex flow with locally symmetric distribution of entropy S and the total enthalpy H per unit mass along a surface normal and also take into account the effect of curvature. If R is the local radius of curvature of the wall and VtE is the velocity component tangential to the body surface, then the flow model satisfies the normal momentum equation .V /2 @p D tE , @E n R

(15)

the nonpenetration boundary condition VnE D VE n D 0, and enforces antisymmetric normal derivative @S=@E n and @H=@E n along the surface normal to the body in the vicinity of the body. These entropy and total enthalpy distributions will produce zero normal derivatives when the flow is irrotational. Using finite difference scheme, the method has shown to produce superior accuracy compared with the classical surface boundary condition by Dadone et al. in [4]. In view of many good characters, the method has also been successfully used in second order accuracy finite volume method at body-fitted unstructured grid by Wang and Sun in [35]. In order to implement the ghost cell methodology for an immersed boundary, we need to proceed by first identifying the cells whose centers are inside the solid and determine the ghost cells. Solid cells can be identified by ray tracing. Because the stencil of RKDG is compact and only demand the nearest neighbor cells, hence, we only need to identify the closest row of solid cells near the body surface as the ghost cells (gc). Those cells are shown in Figure 3. Following this, we should derive a scheme that allows us to calculate the value of the variables at each of these ghost cell centers such

fluid cell

gc

gc

gc

gc

ip1 gc1

Solid body

ip gc

Figure 3. 2D schematic describing ghost cell method used in the current solver. Copyright © 2013 John Wiley & Sons, Ltd.



855

that the boundary conditions on the immersed boundary in the vicinity of the ghost cell are satisfied. Here, the GCM for two-dimensional inviscid flow on Cartesian grids use the following equations: pgc D pip ip gc VtE2gc VnE gc

VtE2ip

n, R 1 pgc D ip , pip 2 pip pgc 2 , D VtEip C 1 ip gc D VnE ip ,

(16)

where ip denotes the image point corresponding to his ghost cell (gc) and n indicates the distance between the ip and gc. Then, the nonpenetration boundary conditions are satisfied automatically. It is well known that the values at the image points ip can be achieved by some interpolation formulas with the surrounding fluid cells. Clearly, using any form of polynomial interpolation techniques across strong discontinuities will result in interpolation errors and may produce pseudo extremes. However, the inverse distance weighting interpolation normally neither generates any new extreme nor produces pseudo oscillations. In the DG frame, DG solutions are available in all the fluid cells, and one can use solution polynomial to interpolate the values at the reflected image points. However, it should be noted that in the presented boundary treatment, one situation will be encountered, which is the image point located in the ghost cell itself; see for instance the point ip1 shown in Figure 3. In this ghost cell, there is no DG solution polynomial. Of course, one can use DG solution in some neighboring fluid cells to extrapolate it. However, there are several fluid cells that are near this ghost cell. In this case, one must give a choice for using the neighboring cells. For convenience, it was chosen all the surrounding fluid cell averages to interpolate the value at the reflected point. It could also be the case that the surrounding interpolation stencil for image point of ghost cell contains other ghost cells. So in this work, for simplicity and robustness, the inverse distance weighting interpolation formula [5] at the image point is used. In order to evaluate the pressure pgc at a ghost cell in the first equation of (16), for an arbitrary curved boundary, it is also needed to estimate the local curvature. Here, we adopt the method suggested by Wang and Sun [35] for two-dimensional flow. Let us consider a curved boundary as shown in Figure 4. To estimate the curvature for body surface AB, first we use AB and the point to the left of the boundary face (point Lft) to make one estimate, and then we use A B and the right point Rgt to perform the other estimate. Then, the final radius can be obtained by a simple average of the two estimates. The local averaged radius is used as the final approximation of the local curvature. However, if either point A or B is a sharp corner, then the estimate would be only calculated by the three points to avoid the sharp corner. Recently, Krivodonova and Berger [29] presented a new boundary treatment on body-fitted unstructured grid to make the numerical velocity at every boundary integration point to coincide with the streamline direction near the body surface, that is, to be orthogonal to the ‘true’ approximated geometrical boundary rather than to the computational boundary. They impose the following condition at every integration point: V NE D 0,

(17)

B Rgt A

Lft

Figure 4. Estimation of local curvature for 2D body. Copyright © 2013 John Wiley & Sons, Ltd.


856

J. LIU ET AL.

approximate wall boundary computational boundary

fluid N n gc solid

Figure 5. Modification of the normal velocity.

in a small vicinity to the surface, where NE is the unit normal to the physical geometry. Using this technique, they did not use curved elements in body-fitted unstructured grid and obtained very accurate results. In this study, we develop the technique presented by Krivodonova to resolve the boundary velocity on adaptive Cartesian grid. In practice, an analytical description of the surface is usually not available. In our computations, we assumed that the only available information is the wall points. We approximated the physical geometry by the approach used to calculate the approximate curvature. In fact, the process is very simple; one only needs to make the normal direction of computational boundary nE in (16) be replaced by NE of approximated solid boundary (see Figure 5) as follows: pgc D pip ip gc D ip V E2

T gc

pgc pip

D V E2 C T ip

V E2

T ip

R

n,

1

,

2 1

pip pgc ip gc

(18)

,

VNE gc D VNE ip , where TE is the unit tangential direction corresponding to the new normal direction of NE . 4. NUMERICAL EXAMPLES In this section, numerical experiments will be carried out to demonstrate the accuracy and the effectiveness of the adaptive Cartesian grid RKDG method with the immersed boundary treatment proposed in this paper. In all numerical tests, an h-adaptive RKDG method with P 2 polynomial is used in each cell. 4.1. Shock vortex interaction In order to show the high order accuracy of the presented RKDG method on adaptive Cartesian grid, we first consider a shock and vortex interaction. This test case is adapted from Shu [36], who used it to show some advantages of high order methods. In order to gain a fine resolution, they often use a refined grid in the x-direction around the shock wave. Copyright © 2013 John Wiley & Sons, Ltd.


857


It contains an interaction between a stationary Mach 1.1 shock positioned at x D 0.5 and normal to the x-axis, and a moving isentropic vortex. The initial condition is set following the exact Rankine– Hugoniot condition, and a vortex is added to the main flow with its center at .x0 , y0 / D .0.25, 0.5/. p The left state of the shock is specified as .p, u, v, / D .1, 1.1 , 0, 1/. The vortex by is described

the perturbations to the velocities .u, v/, temperature .T D p=/, and entropy S D ln p mean flow and has the values

of the

2

u0 D " e ˛.1 / sin , 2 v 0 D " e ˛.1 / cos , 2 2˛.1 2 /

e T 0 D . 1/" 4˛ S 0 D 0,

,

p where D rrc , r D .x x0 /2 C .y y0 /2 , " D 0.3, rc D 0.05, and ˛ D 0.204. The computational domain is taken to be Œ0, 2 Œ0, 1 . We use a base grid of 100 50 with two levels of mesh refinement and set !1 D 1.2, !2 D 1.4, !3 D 0.2, and !4 D 0.4. The reflective boundary conditions are used at the top and bottom boundaries. For this special case, we use hybrid sensors (14) of the divergence and the curl of velocity as the solution adaptive sensors. The pressure contours and adaptive Cartesian grids at different time are shown in Figures 6 and 7. From Figure 6, when the vortex is going through the shock wave, even though its core is contaminated by the shock, the curved shock is still clear, and the restoration of the vortex is perfect. At time t D 0.35, the interactions between the shock wave and the vortex produce a Mach structure. At t D 0.8 (see Figure 7), one branch of the shock bifurcations has reached the top boundary and been reflected. By this time, the reflections near the upper boundaries are well captured. In Figure 8, we compare the results obtained by present h-adaptive RKDG method and the WENO5-LLF finite difference method [37]. Ninety contours are drawn for the pressure component in the range of .1.19, 1.37/. Here the result shown in the right of Figure 8 is calculated by an in-house code in a uniform grid with 250 100 points. This figure reveals that the present method is lower numerical dissipative than the WENO5-LLF finite 1

0.8

0.8

0.8

0.6

0.6

0.6

Y

Y

1

Y

1

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

0

1

0

0.6

X

(a) t = 0.05

(b) t = 0.2

(c) t = 0.35

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0.4

0.6

0.8

1

0.8

1

0.8

1

Y

0.8

Y

1

Y

1

0.2

0.4

X

1

0 0

0.2

X

0 0

0.2

0.4

0.6

0.8

1

0 0

0.2

0.4

0.6

X

X

X

(d) t = 0.05

(e) t = 0.2

(f) t = 0.35

Figure 6. Local views of 30 equally spaced pressure contours and adaptive Cartesian grids at time t = 0.05, t = 0.20 and t = 0.35, respectively, from left to right. Copyright © 2013 John Wiley & Sons, Ltd.


858

J. LIU ET AL. 1

0.8

0.8

0.6

0.6

Y

Y

1

0.4

0.4

0.2

0.2

0 0.4

0.6

0.8

1

1.2

0 0.4

1.4

0.6

0.8

1

X

X

(a) t = 0.6

(b) t = 0.8

0.8

0.8

0.6

0.6

1.4

1.2

1.4

Y

1

Y

1

1.2

0.4

0.4

0.2

0.2

0 0.4

0.6

0.8

1

1.2

0 0.4

1.4

0.6

0.8

1

X

X

(c) Adaptive grid at t = 0.6

(d) Adaptive grid at t = 0.8

Figure 7. Local view of 30 equally spaced pressure contours and adaptive Cartesian grids at time t D 0.6 and t D 0.8 .

1

0.8

0.8

0.6

0.6

Y

Y

1

0.4

0.4

0.2

0.2

0 0.4

0.6

0.8

1

1.2

1.4

0 0.4

0.6

0.8

1

X

X

(a) h-adaptive RKDG

(b) WENO5-LLF [37]

1.2

1.4

Figure 8. Ninety pressure contour levels from 1.19 to 1.37 at t D 0.6 with h-adaptive RKDG and WENO5-LLF [37] methods. Copyright © 2013 John Wiley & Sons, Ltd.


859


difference method [37] and the solution has some slight oscillations. Because of the utilization of AMR, the shock is more thin than WENO5-LLF finite difference method. 4.2. Subsonic flow around a circular cylinder

2

2

1

1

0

0

Y

Y

We consider a subsonic flow at Mach number M a1 D 0.38 through a circular cylinder with radius r D 0.5 and solve the problem using the presented RKDG method on the domain Œ4, 4 Œ4, 4 . Because the case is subsonic flow and the solution is smooth, we do not use any limiter here. In this test, the computational domain is divided by 50 50 base cells. In order to test the convergence of the numerical method, we consider same background grid with four different initial mesh refinement numbers (Nr D 0, 1, 2, 3, see the left of Figures 9 and 10) near the solid boundary. In order to show the ability of the presented method, a circular cylinder is approximated by a piecewise linear approximation (a polygon with 180 equal sides). The results were obtained by using presented RKDG method with ghost cell boundary treatment and integrating until a steady state was reached. The right of Figures 9 and 10 gives the Mach number contours. It should be noted that the plotting was carried out using cell average from inside of elements; no smoothing was applied. It is found that a solution obtained on the finest mesh is almost symmetric and does not have any visible wake. In our experiments, the boundary treatment does not stall the reduction of the residuals to reach a steady state. The number of time steps needed for convergence is increased; as expected, it is

-1

-1

-2

-2 -2

-1

0

1

2

-2

-1

X

0

1

2

1

2

X

2

2

1

1

0

0

Y

Y

(a) Nr = 0

-1

-1

-2

-2 -2

-1

0

1

2

-2

X

-1

0

X

(b) Nr = 1 Figure 9. Local views of computational grid and 16 equally spaced Mach number contours with Nr D 0 and Nr D 1. Copyright © 2013 John Wiley & Sons, Ltd.


860

2

2

1

1

Y

Y

J. LIU ET AL.

0

0

-1

-1

-2

-2

-2

-1

0

1

2

-2

-1

0

X

1

2

1

2

X

2

2

1

1

Y

Y

(a) Nr = 2

0

0

-1

-1

-2

-2

-2

-1

0

1

2

-2

-1

0

X

X

(b) Nr = 3 Figure 10. Local views of computational grid and 16 equally spaced Mach number contours with Nr D 2 and Nr D 3.

0

-5

Nr=0 Nr=1 Nr=2 Nr=3

log(Res.)

-10

-15

-20

-25

0

10000

20000

30000

40000

50000

Iteration

Figure 11. Convergence plot with Nr D 0, 1, 2, 3. Copyright © 2013 John Wiley & Sons, Ltd.



861

Table I. L1 =L2 errors in entropy for the circular cylinder. GCM (Equation (16)) Nr 0 1 2 3

GCM (Equation (18))

L1 ."ent /

L2 ."ent /

L1 ."ent /

L2 ."ent /

3.4276e-02 1.3535e-02 4.2376e-03 1.7150e-03

1.4776e-02 7.7288e-03 1.3669e-03 3.7130e-04

3.3800e-02 1.3397e-02 4.0243e-03 1.5812e-03

1.4307e-02 7.6334e-03 1.3506e-03 3.3752e-04

in inverse proportion to the reduced element size (h-refinement). The convergence histories with different grids are illustrated in Figure 11. To show the convergences of the proposed method, we also measure the L2 and L1 errors in entropy "ent defined as

p "ent D 1 (19) p1 1 where p1 and 1 are pressure and density of the free stream, respectively. Note that the entropy production serves as a good criterion to measure accuracy of the numerical solutions, because the flow under consideration is isentropic. The results with h-refinement are shown in Table I. Clearly, Table I shows that the modified version (18) of GCM produces slightly less errors than the original one (16). For instance, the modified version GCM method has about a 10% reduction in L2 error and an 8% reduction in L1 error (in the finest mesh) compared with the ones in the GCM method. 4.3. Shock diffraction over a circular cylinder This example concerns the shock/cylinder interaction problem, which is one of the well-studied test cases in the literature [10, 11, 38]. The flow consists of a planar shock moving at M a D 2.81 and impingingon a solid cylinder. The state variables at the right of the shock at time t D 0 are ., p, u, v/ D 1, 1 , 0, 0 , and quantities at the left of the shock are computed from the ones at the right side and the Mach number of shock wave using classical shock relations. The computational domain is Œ0, 1 Œ0.5, 0.5 , which is similar to the case presented in [10]. Wall boundary conditions are used for the top and bottom, and the Neumann condition is used for the outlet. At the inlet, a Dirichlet condition corresponding to a Mach 2.81 shock is enforced. The circular obstacle with a radius of 0.105 was positioned at x D 0.5 and y D 0. The initial discontinuity was located at x D 0.385. 1 The simulations were performed with a coarse base mesh of size x D y D 80 and four levels of mesh refinement with !1 D 1.1, !2 D 1.5, !3 D !4 D 0.3 was selected. Figure 12 shows the adaptive grids and the density contour images from the simulation of the shock and cylinder interaction at different times t D 0.2 and t D 0.5. The evolution of diffraction process includes regular reflection, transition to Mach reflection, complex wake/shock, and shock/shock interactions. The contact discontinuities, upper and lower triple points, complex wake structure, and the Mach systems are correctly reproduced by the present simulation. Furthermore, thanks to the high order accuracy RKDG method, the onset of instability was found behind the cylinder in the right of Figure 12. The computed trajectories of the upper triple point (denoted by AMRDG) are depicted in Figure 13. The present simulation agrees well with the numerical calculation presented in [10] and the experimental correction obtained in [39, 40]. 4.4. Supersonic flow around a triangle In this section, we consider the same example presented in [5, 38], where a supersonic flow past a solid body of triangular shape with height h D 0.5 and half angle D 20 deg, as shown on the left of Figure 14. The free stream Mach number is 2. In this geometry, a special procedure will be encountered, which is the problem of multivalued ghost points. It is often found near unresolved Copyright © 2013 John Wiley & Sons, Ltd.


862

J. LIU ET AL.

(a) t = 0.2

(b) t = 0.5 Figure 12. Adaptive grids and density pseudocolor images for a Mach 2.81 shock wave diffracting off a 1 stationary cylindrical obstacle, computed on a base mesh of size x D y D 80 with four levels of mesh refinement at time t D 0.2 and t D 0.5.

0.5 AMRDG Exp. Ref.

0.4

Y

0.3

0.2

0.1

0

0.4

0.5

0.6

0.7

0.8

0.9

1

X

Figure 13. Comparison of triple point trajectory. Copyright © 2013 John Wiley & Sons, Ltd.



863

1

0.5

0

Y

theta h

-0.5

-1 -0.5

0

0.5

1

1.5

X

Figure 14. Geometrical configuration for oblique shock wave analysis and the local view of the adaptive Cartesian grid.

thin surfaces and sharp corners, such as the cases at the sharp trailing edge of an airfoil or near the apex of triangle as shown in Figure 14. At sharp corners, one cell center inside the geometry may be the ghost cell center for one side of the corner surface, as well as for the other sides. Also, a ghost cell center pertaining to one side of a corner surface may be located inside the flow field on the other side of the corner. To handle this case, we first find the fluid cells, which are needed for multiple valued ghost points, and then set flags. Generally speaking, the ghost boundary conditions are given before time evolution. In this special situation, the multiple valued ghost points are evaluated during the time evolution, which accords to the fluid cell to find the corresponding ghost cell value. The rest of fluid cells are computed as usual. The advantage of this approach is that there is no need to allocate new memory to save the values on the multiple valued ghost points. The problem was solved by the present adaptive ghost cell Cartesian grid method on the domain Œ4, 5.5 Œ4.5, 4.5 . Initial grids (the domain is divided by 60 60 cells) have been refined six times near the body, then three levels of grid refinements are implemented during solution process with !1 D 1.1, !2 D 1.5, !3 D !4 D 0.3. The final mesh is shown on the right of Figure 14. The total number of grid points is 31,475. Figure 15 shows the Mach number contours and residual history obtained by the presented method. The positions of the attached shock waves are clearly indicated in the plots. 1 -2 MA

Y

0.5

0

-0.5

-1

-3 log(Res.)

3.8 3.6 3.4 3.2 3 2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2

-4

-5

-6

-7 -0.5

0

0.5

1

X

1.5

0

10000

20000

30000

40000

Iteration

Figure 15. Mach number contours and residual history. Copyright © 2013 John Wiley & Sons, Ltd.


864

J. LIU ET AL.

An oblique shock can be attached or detached depending on the values of the deflection angle and the upstream Mach number M1 . If the shock is attached to the triangle, its angle with the horizontal, ˇ, can be computed through the following relationship given by [38]: " # M12 sin2 ˇ 1 tan D 2 cot ˇ . (20) M12 . C cos 2ˇ/ C 2 For M1 D 2 and D 20 deg, the angle ˇ would be ˇ D 53.46 deg from Equation (20). Boiron et al. [38] presented numerical values of ˇ by using a high-resolution penalization method to solve viscous Navier–Stokes equations. Their numerical values of ˇ are 54.13 deg on 5122 grid points, 53.70 deg on 10242 grid points, and 53.56 deg on 15362 grid points, respectively. In our results for the Euler equations on the 31,475 grid points, which is much less than their grid points, the mean shock angle of ˇ is about ˇ 52.69 deg, which is comparable with theirs. 4.5. Transonic flow past a NACA0012 airfoil Now, let us consider an application of the present method on an airfoil. Because the airfoil has a sharp trailing edge, it is difficult to define the ghost cell boundary as the aforesaid triangle. In the transonic flow computation, the Mach number is 0.8 and the angle of attack 1.25 deg around the NACA0012 airfoil (chord = 1). The computational domain for this case is Œ6, 6.5 Œ6.5, 6.5 . The initial mesh .80 80/ was first refined four times near the body boundary. Four levels of solution-based refinement are carried out during the time evolution. The mesh is shown in the left of Figure 16. It should be noted that a sensor of the divergence of velocity d is only used in the process of solution-based refinement. The coefficients of sensors for adaptive refinement are !1 D 1.1 and !3 D 0.2. The density contours are shown on the right of Figure 16. The calculated pressure coefficient denoted by AMRDG is plotted on the left of Figure 17. The figure also gives a comparison with some references. In Figure 17, the result flagged by SU2 is calculated by Stanford University unstructured code SU2 [41], which uses a second order finite volume scheme on body-fitted unstructured grid to discretize the Euler equations. As shown in Figure 17, because of the usage of adaptive Cartesian grid, except that the calculated shock is more thin than SU2, the present result matches very well with the one calculated by SU2; both the strong shock and the weak shock are well resolved, and the positions of shock are also well located. Sj ö green and Petersson recently developed an embedded boundary method to solve Euler equations on Cartesian grid [2]. Using 1200 1200 grids, they obtained the pressure coefficients as shown in Figure 17. Their result slightly misses the position of the shock. We also compared our result with the one calculated by a body-fitted triangular grid RKDG method with WENO limiter by Zhu et al. in [42]. From the same figure, we find that heir result gives the right position of shock, but their pressure coefficients 1.5

1 0.5

Y

Y

0.5

0

0

-0.5 -0.5 -1

0

0.5

X

1

1.5

-1

-0.5

0

0.5

1

1.5

2

X

Figure 16. Adaptive Cartesian grid and 20 equally spaced density contours. Copyright © 2013 John Wiley & Sons, Ltd.



865

-5 -1 -6

log(Res.)

Cp

-0.5

0

-7

-8

0.5 AMRDG SU2 Sjogreen Zhu

1

-9

-10 1.5 0

0.2

0.4

0.6

0.8

1

0

20000

40000

60000

Iteration

X

Figure 17. Surface pressure coefficient distribution and residual history.

are slightly underestimated. Those comparisons show that the present method can give very good result, which is comparable with the general body-fitted finite volume method. At times, our result is even better than the one obtained by a high order RKDG method on body-fitted triangular grid. The relations between convergence histories and the number of iterations are shown on the right of Figure 17. 4.6. Supersonic flow through two staggerd NACA0012 airfoils Finally, we consider a supersonic flow through two staggered NACA0012 configuration with Mach number 2 and the attack of angle 0. They are shifted by a distance equal to half a chord in both the parallel and perpendicular directions. This case includes some complex fluid characters, such as shock reflection, Mach reflection, and contact discontinuity. The computational domain used in this case is also Œ6, 6.5 Œ6.5, 6.5 . The initial mesh .80 80/ was firstly refined four times near the body boundary. During the time evolution, we have restricted that the grid refinement does not exceed six levels, and the coefficients of sensors for solution adaptive refinement are !1 D 1.1, !2 D 1.4, !3 D !4 D 0.4. The adaptive Cartesian grid is displayed in Figure 18. Thanks to the combination of velocity and entropy as the adaptive sensors, the regions of contact discontinuity are also refined during the grid refinement. The computed Mach number contours and

1

Y

0.5

0

-0.5

0

0.5

1

1.5

2

X

Figure 18. Local view of the adaptive Cartesian grid. Copyright © 2013 John Wiley & Sons, Ltd.


866

J. LIU ET AL.

1.5 -0.5 1 0 0.5

Cp

Y

0.5 0

1 -0.5 AMRDG Li

1.5 -1 2 -1.5 -0.5

0

0.5

1

1.5

2

2.5

0

0.5

X

1

1.5

X

Figure 19. Local view of the Mach number contours and the surface pressure coefficient distribution.

the pressure coefficient distribution on the surface are shown in Figure 19. Clearly, the contact and shock wave are captured and resolved very well. It is also noted that the pressure coefficients match well with the one calculated by a nonet-Cartesian grid method [43]. 5. CONCLUSION The ghost cell Cartesian grid method provides an efficient and flexible alternative to the traditional body-fitted grid techniques. The RKDG method is an efficient high order scheme, which can easily deal with grid refinement and retain local conservation even in the case of hang node. In this paper, an h-adaptive Cartesian grid RKDG method is developed. A discontinuous detector and limiting technique is used to suppress oscillation, and a solution-based grid adaptive method with hybrid sensors is used during the grid refinement. In order to solve flow with complex geometry, in this paper, an h-adaptive Cartesian grid RKDG method with ghost cell immersed boundary method has been successfully developed and used to simulate steady and unsteady problems of the Euler equations. Although the boundary treatment is not conservative and at most in second order accuracy, numerical results show that the present method can achieve correct shock positions, high resolution, and some correct physical characters. It is demonstrated that the present method is very effective.

ACKNOWLEDGEMENTS

The research was partially supported by the National Science Foundation of China (No.11102179, No. 91230110), DMU PhD Studentship (De Montfort University, UK), and Science Foundation of Institutions of Higher Education of Jiangsu Province (No. 11KJB110016), China.

REFERENCES 1. Yang G, Causon DM, Ingram DM, Saunders R, Batten P. A Cartesian cut cell method for compressible flows part A: static body problems. Aeronautical Journal 1997; 101(1001):47–56. 2. Sj ö green B, Petersson N. A Cartesian embedded boundary method for hyperbolic conservation laws. Communications in Computational Physics 2007; 2(6):1199–219. 3. Forrer H, Jeltsch R. A higher-order boundary treatment for Cartesian-grid methods. Journal of Computational Physics 1998; 140(2):259–77. 4. Dadone A, Grossman B. Ghost-cell method for inviscid two-dimensional flows on Cartesian grids. AIAA Journal 2004; 42(12):2499–507. 5. Liu JM, Zhao N, Hu O. The ghost cell method and its applications for inviscid compressible flow on adaptive tree Cartesian grids. Advances in Applied Mathematics and Mechanics 2009; 1(5):664–82. Copyright © 2013 John Wiley & Sons, Ltd.



867

6. Ji H, Lien FS, Yee E. A robust and efficient hybrid cut-cell/ghost-cell method with adaptive mesh refinement for moving boundaries on irregular domains. Computer Methods in Applied Mechanics and Engineering 2008; 198:432–48. 7. Berger MJ, Helzel C, LeVeque RJ. H-box methods for the approximation of hyperbolic conservation laws on irregular grids. SIAM Journal on Numerical Analysis 2003; 41(3):893–918. 8. Helzel C, Berger MJ, LeVeque RJ. A high-resolution rotated grid method for conservation laws with embedded geometries. SIAM Journal of Scientific Computing 2005; 26(3):785–809. 9. Colella P, Graves DT, Keen BJ, Modian D. A Cartesian grid embedded boundary method for hyperbolic conservation laws. Journal of Computational Physics 2006; 211(1):347–66. 10. Sambasivan SK, UdayKumar HS. Ghost fluid method for strong shock interactions Part 2: immersed solid boundaries. AIAA Journal 2009; 47(12):2923–37. 11. Chaudhuri A, Hadjadj A, Chinnayya A. On the use of immersed boundary methods for shock/obstacle interactions. Journal of Computational Physics 2011; 230(5):1731–48. 12. Cecil T, Osher S, Qian J. Essentially non-oscillatory adaptive tree methods. Journal of Scientific Computing 2008; 35(1):25–41. 13. Li ST, Hyman JM. Adaptive mesh refinement for finite difference WENO schemes. Los Alamos Report LA-UR-0389272003, 2003. 14. Shen CP, Qiu JM, Christlieb A. Adaptive mesh refinement based on high order finite difference WENO scheme for multi-scale simulations. Journal of Computational Physics 2011; 230(10):3780–802. 15. Cockburn B, Shu CW. The Runge-Kutta local projection P 1 -discontinuous Galerkin finite element method for scalar conservation laws. Mathematical Modelling and Numerical Analysis 1991; 25(3):337–61. 16. Cockburn B. Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework. Mathematics of Computation 1989; 52(186):411–35. 17. Cockburn B, Shu CW. The Runge-Kutta discontinuous Galerkin finite element method for conservation laws V: multidimensional systems. Journal of Computational Physics 1998; 141(2):199–224. 18. Cockburn B, Lin SY, Shu CW. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: one dimensional systems. Journal of Computational Physics 1989; 84(1):90–113. 19. Cockburn B, Hou S, Shu CW. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case. Mathematics of Computation 1990; 54(190):545–81. 20. Flaherty JE, Loy RM, Shephard MS, Szymanski BK, Teresco JD, Ziantz LH. Adaptive local refinement with octree load-balancing for the parallel solution of three-dimensional conservation laws. Journal of Parallel and Distributed Computing 1997; 47:139–52. 21. Biswas R, Devine KD, Flaherty JE. Parallel, adaptive finite element methods for conservation laws. Applied Numerical Mathematics 1994; 14:255–83. 22. Remacle JF, Flaherty JE, Shephard MS. An adaptive discontinuous Galerkin technique with an orthogonal basis applied to compressible flow problems. SIAM Rivews 2003; 45(1):53–72. 23. Devine KD, Flaherty JE. Parallel adaptive hp-refinement techniques for conservation laws. Applied Numerical Mathematics 1996; 20:367–86. 24. Dedner A, Makridakis C, Ohlberger M. Error control for a class of Runge-Kutta discontinuous Galerkin methods for nonlinear conservation laws. SIAM Journal on Numerical Analysis 2007; 45(2):514–38. 25. Hartmann R, Houston P. Adaptive discontinuous Galerkin finite element methods for nonlinear hyperbolic conservation laws. SIAM Journal of Scientific Computing 2002; 24(3):979–1004. 26. Zhu HQ, Qiu JX. Adaptive Runge-Kutta discontinuous Galerkin methods using different indicators: one-dimensional case. Journal of Computational Physics 2009; 228(18):6957–76. 27. Zhu HQ, Qiu JX. An h-adaptive RKDG method with troubled-cell indicator for two-dimensional hyperbolic conservation laws. Advances in Computational Mathematics. DOI: 10.1007/s10444-012-9287-7. Published online 30 November 2012. 28. Bassi F, Rebay S. High-order accurate discontinuous finite element solution of the 2D Euler equations. Journal of Computational Physics 1997; 138(2):251–85. 29. Krivodonova L, Berger M. High-order accurate implementation of solid wall boundary conditions in curved geometries. Journal of Computational Physics 2006; 211(2):492–512. 30. Toro EF, Spruce M, Speares W. Restoration of the contact surface in the HLL Riemann solver. Shock Waves 1994; 4(1):25–34. 31. Shu CW, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics 1988; 77(2):439–71. 32. Berger M. Local adaptive mesh refinement for shock hydrodynamics. Journal of Computational Physics 1989; 82(1):64–84. 33. Krivodonova L, Xin J, Remacle JF, Chevaugeon N, Flaherty JE. Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws. Applied Numerical Mathematics 2004; 48:323–38. 34. De Zeeuw D. A quadtree-based adaptively-refined Cartesian-grid algorithm for solution of the Euler equations. PhD Thesis, University of Michigan, 1993. 35. Wang ZJ, Sun YZ. Curvature-based wall boundary condition for the Euler conditions on unstructured grids. AIAA Journal 2003; 41(1):27–33. Copyright © 2013 John Wiley & Sons, Ltd.


868

J. LIU ET AL.

36. Shu CW. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. NASA CR-97-206253, 1997. 37. Balsara DS, Shu CW. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. Journal of Computational Physics 2000; 160(2):405–52. 38. Boiron O, Chiavassa G, Donat R. A high-resolution penalization method for large Mach number flows in the presence of obstacles. Computers & Fluids 2009; 38(3):703–14. 39. Ripley R, Lien F, Yovanovich M. Numerical simulation of shock diffraction on unstructured meshes. Computers & Fluids 2006; 35(10):1420–31. 40. Kaca J. An Interferometric investigation of the diffraction of a planar shock wave over a semicircular cylinder. Inst for Aerospace Studies TR 269, Univ of Toronto, 1988. 41. Palacios F, Alonso JJ, Duraisamy K, Colonno MR, Aranake AC, Campos A, Copeland SR, Economon TD, Lonkar AK, Lukaczyk TW, Taylor TWR. Stanford University Unstructured (SU2): An open-source integrated computational environment for multi-physics simulation and design. AIAA Paper 2013-0287, 51st AIAA Aerospace Sciences Meeting and Exhibit, Grapevine, Texas, USA, January 7th–10th, 2013. 42. Zhu J, Qiu JX, Shu CW, Dumbser M. Runge-Kutta discontinuous Galerkin method using WENO limiters II: unstructured meshes. Journal of Computational Physics 2008; 227(9):4330–4353. 43. Li K, Wu ZN. Nonet-Cartesian grid method for shock flow computations. Journal of Scientific Computing 2004; 20(3):303–29.

Copyright © 2013 John Wiley & Sons, Ltd.
