Shock capturing schemes with local mesh adaptation for high speed compressible flows on three dimensional unstructured grids

Computers & Fluids 70 (2012) 126–135

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Shock capturing schemes with local mesh adaptation for high speed compressible flows on three dimensional unstructured grids Vinh-Tan Nguyen ⇑, Hoang Huy Nguyen, Matthew Antony Price, Jaik Kwang Tan Institute of High Performance Computing, 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore

a r t i c l e

i n f o

Article history: Received 16 February 2012 Received in revised form 7 August 2012 Accepted 9 September 2012 Available online 2 October 2012 Keywords: Computational fluid dynamics Shock capturing Mesh adaptation Unstructured grids Solution adaptivity

a b s t r a c t This paper contributes towards a more complete approach to capture very strong shocks in various applications of high speed compressible Navier–Stokes flows including blasts and explosions using second order finite volume method on unstructured grids. The HLLC Riemann solver is employed to solve for fluxes at cell interfaces with second order approximation of local Riemann states, thus obtaining second order accuracy. In order to stabilize solutions due to high order approximation of solutions in the presence of discontinuities, several strategies are presented in this work. Slope limiters are first explored on unstructured grid to maintain monotonicity of the solution reconstruction following local extremum diminishing (LED) or total variation diminishing (TVD) criteria. The hybrid HLLC/HLLE scheme is appended to eliminate shock instabilities in very strong shock cases. To improve resolution of shocks, a local mesh adaptation scheme is used to increase mesh resolution in areas of high gradients. The scheme only regenerates mesh locally and is proven to be robust and efficient for capturing of unsteady shock propagation applications. Comparisons on the accuracy and performance of different methods on various applications are drawn to suggest a more robust and efficient method for capturing shocks on unstructured grids. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction High speed compressible flows are fundamental to many engineering problems involving hypersonic flows, blasts and explosions, shock wave propagations. A wide range of numerical methods have been developed for simulations of compressible flows in the presence of very strong shocks. In particular, second order finite volume methods on unstructured grids have become desirable computations in many industrial applications. In these high speed compressible flow applications, strong shocks in the solutions present a formidable challenge in resolving discontinuities. It is then essential to derive an effective approach to accurately capture flow features including strong shocks. In capturing shocks, Godunov-type schemes [4] has been popular for its conservative treatment of resolving discontinuities. Under a basic assumption of piece-wise constant solution across elements in the Godunov-type methods, intercell fluxes need to be constructed from the current solution in order to compute the solution in the next stage. The intercell fluxes are solutions to the so-called local Riemann problem defined at cell interfaces. The local Riemann problems may be solved analytically, if desired. However, existence of exact solutions to the Riemann problems is ⇑ Corresponding author. Tel.: +65 6419 1591; fax: +65 6467 4350. E-mail address: [email protected] (V.-T. Nguyen). 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.09.011

not always trivial; in many cases, iterative methods are employed to solve for intercell fluxes. Alternatively, solutions to local Riemann problems can be solved approximately with the original motivation that it can provide the solutions more cheaply than exact solvers. Since the pioneering work by [10,11], the research on approximate Riemann solvers has received considerable progress over the last few decades and formed a foundation for extensive studies on shock capturing. The combination of Godunov-type schemes with approximate Riemman solvers are widely adopted and become a standard for benchmarking of any new schemes. Among those Roe’s scheme [10], Harten–Lax–van Leer (HLL) schemes [3] and its variants are popular choices in many computations. Despite huge success over a large class of problems, there are fundamental deficiencies of the Godunov-type framework, especially for very strong shock applications. In [1] it reported and classified failures of various approximate Riemann solvers on multidimensional problems. The limitations of the Riemann solvers to shock-capturing properties was catalogued according to their failing modes including expansion shock, negative internal energy, slowly-moving shock, the carbuncle phenomenon, kinked Mach stem, and odd–even decoupling. Details of these problem can be found in the reference cited therein [1] in which the author proposed a strategy of using combined fluxes with more dissipative flux such as HLLE in the shock area. In [2], the so-called shock instability in multidimensional

127

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

flows and the nonexistence of a solution for strong receding flows are major problems in such failings of approximate Riemann solvers. The failures were attributed to the dissipative term in mass fluxes defined in [2]; thus it was suggested that making the pressure contribution in dissipative terms vanish by controlling the magnitude of wave speed forms a shock-stable scheme. Having mentioned the failures of approximate Riemann solvers in the framework of Godunov-type schemes, it is worth noting that low resolution is also considered one of the biggest shortcomings of such schemes. For first order methods, it requires very fine grids for better resolutions, particularly for shock wave interactions. Higher order accurate methods are thus desirable for numerical simulations of high speed compressible flows. Extension of first order finite volume methods to higher order is a favourable option since it only requires one to employ higher order reconstructions of solutions before feeding them to local Riemann solvers. While second order finite volume schemes have been widely adopted in many industrial codes, higher-than-second-order methods remains a formidable challenge for large scale problems, especially the construction of such schemes on unstructured grids. There has been much progress in this research direction over the last few decades including the development of essentially non-oscillatory (ENO) methods, discontinuous Galerkin methods [13] and their variants. However, it is necessary to stabilize the solution to eliminate oscillations due to high order approximations as discontinuities develop in the solutions. In many cases, the inherent diffusion from numerical discretizations is not sufficient and a special treatment is required to prevent oscillations from growing in the solutions. Among various treatments, artificial viscosity and limiters are frequently employed to stabilize shocks and discontinuities. In using artificial viscosity [12], numerical dissipation is directly added into systems to help stabilize the discretization. The approach has been widely used with successful results for shock capturing in various discretization methods; however, adding viscosity certainly reduces stability of the system especially for very strong shock cases. Alternatively, in using limiters, the solution is stabilized by reducing the interpolating order across the shocks; therefore it may affect the order of accuracy in the vicinity of the shock front. Essentially, employment of stabilizing approaches will introduce deficiencies to high order discretizations. Despite considerable progress in deriving high order methods for high speed compressible flows, one has yet found robust and efficient remedies fully addressing all the above issues related to this type of applications. In this work, we first review our numerical method of finite volume discretization for compressible Navier–Stokes flows on unstructured grids. In the next section, a brief description of the method is given. The edge-based vertex-centered finite volume approach requires the construction of dual mesh by connecting edge midpoints, element centroids, face centroids in such a way that only one grid node is contained in each control volume. Edge-based data structure together with median dual control volume efficiently facilitate discretization of inviscid and viscous fluxes of compressible Navier–Stokes equations. In particular, HLLC flux scheme is employed as an approximate Riemann solver for solving flux functions at volume interfaces. The HLLC scheme is appended by a second order approximation of local Riemann states; thus obtaining second order accuracy. The second order reconstruction in turn makes the scheme unstable for strong shocks and discontinuities. In the next section, various treatments to stabilize the solutions are proposed and implemented including various slope limiters and a hybrid HLLC approach to overcome instability of the second order method for strong shocks. These remedies essentially reduce the approximation to first order accuracy in the vicinity of discontinuities. Therefore, we proposed an robust and effective mesh adaptive approach to increase the resolution in

areas of steep gradients. In our proposed approach, a solution sensor is first designed to detect regions of high gradient; subsequently a robust adaptivity procedure was applied to regenerate three dimensional grid in those regions. Employing this local mesh adaptivity strategy, it is able to enhance the resolution of solutions in the presence of strong shocks with affordable computational cost. 2. Unstructured grid finite volume discretization 2.1. Governing equations The unsteady inviscid compressible flows are governed by the time-dependent, Euler equations on a three-dimensional Cartesian domain X  R3 , with surface @ X, can be expressed in integral form as

Z X

@U dx þ @t

Z

F j nj dx ¼ 0;

ð1Þ

@X

where the conventional summation is employed and nj is the outward unit normal vector to @ X. The unknown vector of the conservative variables and inviscid flux tensor are given by

0

1

0

1

quj q C B qu u þ pd C B B 1 j B qu1 C 1j C C C B B C; F j ¼ B qu2 uj þ pd2j C: U¼B q u 2 C C B B C C B B @ qu3 uj þ pd3j A @ qu3 A uj ðq þ pÞ q

ð2Þ

Here q denotes the fluid density, ui the i’th component of the velocity vector and  the specific total energy. The system is closed by assuming the gas to be calorically perfect, thus setting

p ¼ qRT;

ð3Þ

and

 ¼ cv T þ

1 uk uk ; 2

ð4Þ

where R is the real gas constant and cv = cp  R is the specific heat at constant volume. In this expression, cp is the specific heat at constant pressure. Throughout this work, the ratio of the specific heats,

c¼

cp ; cv

ð5Þ

is set to c = 1.4 which is the value for air at standard conditions. 2.2. Edge-based vertex-centred unstructured FV method The computational domain X is subdivided into a set of non-overlapping tetrahedral elements using a Delaunay mesh generation process with automatic point creation. To enable the implementation of a cell vertex finite volume solution approach, a median dual mesh is constructed by connecting edge midpoints, element centroids and face centroids such that only one node is present in each control volume. Each edge of the grid is associated with a segment of the dual mesh interface between the nodes connected to the edge. The dual mesh interface inside the computational domain surrounding node I is denoted CI, while the parts of the dual situated on the computational boundary are termed CBI , so that CI \ CBI ¼ ;. The triangular facets which define the control volume interface surrounding node I are denoted by CKI . In three dimensions, the segment of the dual mesh associated with an edge is a surface. This surface is defined using triangular facets, where each facet is connected to the midpoint of the edge, a neighbouring element centroid and the centroid of an element face connected to the edge. With this dual mesh definition, the control

128

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

preservation as well as resolving isolated shocks exactly. However, the original HLL scheme cannot resolve contact discontinuities exactly. Building on the HLL flux, HLLC scheme restores the missing contacts and shear waves by assuming a wave configuration consisting of three waves (contact, shock or rarefaction waves) separating four constant states. The approximated flux function is constructed by averaging intermediate states in the exact solution such that isolated shocks and contact discontinuities are exactly resolved. The HLLC flux across the facets CIJ is computed as follows:

F IJHLLC

Fig. 1. Illustration of the dual mesh surrounding an internal node I. The dual is n o constructed by a closed set of planar triangular facets CKI : Each facet only touches a single edge and the set of facets touching the edge spanning between nodes I and J is termed CIJ.

volume can be thought of as constructed by a set of tetrahedra with base on the dual mesh. Fig. 1 shows a sample control volume at point I and its notations of faces and neighbouring cells. The FV discretization transforms surface and volume integrals into a sum of face and control volume integrals and approximates them to second order accuracy. To perform the numerical integration of the inviscid fluxes, a set of coefficients is calculated for each edge using the dual mesh segment associated with the edge. The values of these coefficients for an internal edge are given by the formula

C IJj  nIJj ¼

X

CK

ACK nj I ;

K2CIJ

ð6Þ

I

DIJj ¼

X K2CBIJ

CKI

ACK nj ; I

ð7Þ

where CBIJ is the set of dual mesh facets on the computational boundary touching the edge between nodes I and J. The contribution of the inviscid flux over the control volume surface for node I is then computed as

Z

F j nj dx 

@ XI

X F IJ ;

ð8Þ

J2KI

if SLIJ P 0 if SLIJ < 0 6 SIJ if SIJ < 0 6 SRIJ

;

ð10Þ

if SRIJ < 0

R R where U LIJ  U I , U L IJ , U IJ and U IJ  U J are the four constant state data, respectively. The acoustic wave speeds SLIJ ; SRIJ and SIJ are given as in [4],



SIJ ¼







qJ qRIJ SRIJ  qRIJ  qI qLIJ SLIJ  qLIJ þ pI  pJ 





qJ SRIJ  qRIJ  qI SLIJ  qLIJ



:

ð11Þ

One can refer to [4,15] for more details on the choices of wave speed for HLLC solver. In the above expressions,

qLIJ ¼ uIj nIJj ;

qRIJ ¼ uJj nIJj ;

ð12Þ

R and U L IJ ; U IJ are the conserved variables in the star region separated by the contact,

1 1   C B j B pK pK nIJ C IJ C !B K   K B uK þ C SIJ  qIJ B j C qK SKIJ qKIJ C; B K  B C SIJ  SIJ B C S p qKIJ C B EK pK A @ q þ IJ  IJ K  0

0

1

qK

B C K C B U K IJ  @ ðquj Þ A ¼ qK K

CKI

where ACK is the area of facet CKI and nj is the outward unit normal I vector of the facet from the viewpoint of node I, Fig. 1. In a similar way, the contribution to boundary integrals from the computational boundary can be expressed using a coefficient

8   L > > > F U IJ > >     > > L > < F U LIJ þ SL U L IJ  U IJ     ¼ R > > F U RIJ þ SR U R > IJ  U IJ > > >   > > : F U RIJ

ðqEÞ

K

ð13Þ

qK SKIJ qKIJ

where K = L, R denotes either left or right states. The pressure in the star region is defined as

   K K pK qKIJ  SIJ þ pK : IJ ¼ qK qIJ  SIJ

ð14Þ

In the above Eq. (13), the state solutions in the star region are computed with an assumption of piece-wise constant solutions across elements, thus rendering the scheme to first order accurate. For second order HLLC scheme, it requires second order approximation of left and right conserved variables for local Riemann problems at the facets CIJ. The approximation is expressed as follows:

e L ¼ U L þ dU IJ U IJ IJ

ð15Þ

where KI denotes the set of nodes connected to node I by an edge and F IJ is the inviscid numerical flux in IJ direction obtained from solving local Riemann problems at the facets of node I’s control volume,

where dUIJ and dUJI are the gradient approximations of solution from node I and J, respectively

F IJ  F  nIJ :

dU IJ  $U I  hIJ :

ð9Þ

In this current work, the numerical flux is computed using HLLC flux scheme. 2.3. Second order HLLC flux scheme The HLLC flux scheme [4] is a modification of the original HLL scheme [3], a very efficient and robust Riemann solver which has several good properties including entropy satisfaction, positivity

e R ¼ U R þ dU JI ; U IJ IJ

ð16Þ

Reconstruction of solution gradients can be done by various ways in the framework of finite volume methods. The most straight forward way of computing gradients of solutions is employing weighted average using midpoint rule, h j UI

r

2 3 IJ X IJ 1 4X C j ¼ Dj U I 5: ðU I þ U J Þ þ V I J2C 2 B I

J2CI

ð17Þ

129

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

This results in a second order approximation on simplex meshes satisfying minimal criteria and smooth structured meshes; however, it requires a wide stencil of five points in each direction. Alternatively, one can use an edge-based compact scheme as in [9] where derivatives at midpoint are computed by splitting into two components, normal and tangential direction to the edge. Following [8,14], a corrected least square method can be used for gradient reconstruction on unstructured grids. The basic idea is to cast the gradient reconstruction into a least-square formulation and find the approximation of gradient rUI at given vertex I such that

hIJ UJ  UI  rU I ¼ ; jhIJ j jhIJ j

8J 2 KI :

ð18Þ

This problem can in turn be written as a minimisation problem:

 W I ¼ arg

 min kA I W  DU I k2 ;

ð19Þ

W2R3 ;J2KI h

where ðA I ÞIJ ¼ jhIJIJ j and ðDU I ÞJ ¼

U J U I . jhIJ j

Gradient at node I is then set

rUI  WI. Solution to this weighted least square problem can be found by using the iterative process as proposed in [8]: U 0II

1. Compute initial guess of the gradient, r

using (18),

 1 rU 0I ¼ A TI A I A TI DU I :

ð20Þ

2. Iteratively compute least square estimation of gradient rU nI until convergence



rU nI ¼ A TI A I

1

  1 : A TI DU I þ DrU n1 I 2

ð21Þ

 1 where A TI A I A TI is known as pseudo-inverse of A I . In this work, both reconstruction schemes, weighted average and iterative least square methods, were implemented and tested out for strong shock applications. It is known that second order reconstructions make the numerical scheme unstable in the presence of strong shocks thus requiring stabilization of solutions. Several treatments were proposed and implemented in the present work including various slope limiters and a hybrid HLLC approach to overcome instability of the second order method for strong shocks. As these treatments essentially reduce the approximation to first order accuracy in the vicinity of discontinuities, we proposed a robust and effective mesh adaptive approach to increase the resolution in the areas of steep gradients.

be consistent, dUIJ = dUJI. This condition can be achieved by using minmod function,

8 > < minj ; aj > 0; 8j minmodða1 ; a2 Þ ¼ maxj ; aj < 0; 8j : > : 0; otherwise

ð24Þ

Monotonicity preservation is ensured by either the total variation diminishing (TVD) [17] or local extremum diminishing (LED) [16] conditions. While it is more favourable to express those conditions on Cartesian structured grids, see [4], it is though more complex for the case of unstructured grids. In the current work, several limiters were implemented for various tests. 3.1.1. LED limiter For unstructured grids, [16] derived a LED based monotonicity preserving condition as

 U I hIK 6 maxðU J  U I Þ; minðU J  U I Þ 6 r J2KI

J2KI

8 K 2 KI :

ð25Þ

The Barth and Jespersen (BJ) limiter is modified as in [8]

8    maxJ ðU J U I Þ IK > sg n maxrJUðUI hJ U ; if rU IK hIK > 0 > r U h Þ > I IK I <    minJ ðU J U I Þ UIK ¼ IK sg n minrJ UðUI hJ U ; if rU IK hIK < 0 ; > rU I hIK IÞ > > : 1; otherwise with the sigmoid function defined as sg n ðtÞ ¼

t ð1þt n ÞðnÞ 1

ð26Þ

. This modified

limiter is directional and obtains gradual gradient normalisation as compared to the original BJ limiter which might introduce excessive dissipation to the discretization on unstructured grids. 3.1.2. TVD limiters In the framework of unstructured grids, implementation of TVD limiters are challenging due to difficulty in identification of directional neighbour information. In [21], several TVD limiters were proposed on unstructured grids where the information on upwind and downwind direction were identified locally for each element’s interfaces with its neighbours. In our current implementation, it is assumed without loss of generality that upwind direction is from node I to J at their control volume intersection. The reconstructed gradient is thus computed as

dU IJ  /I ðRIJ ÞrU I  hIJ ;

ð27Þ

where RIJ is defined as

2rU I  hIJ  1: UJ  UI

3. Stabilization methods for strong shocks

RIJ ¼

3.1. Slope limiters

This is used with various TVD slope limiters including common ones such as:

As stated above, the second order reconstruction will produce spurious oscillations in the vicinity of discontinuities and high gradients. Therefore, slope limiters are required to ensure non-oscillatory solutions. The slope limiter is generally employed to limit the gradient as follows:

 U I ¼ / rU I : r I

ð22Þ

van Leer :

/VL ðRÞ ¼

ð28Þ

4R

ðR þ 1Þ2 2R van Albada : /VA ðRÞ ¼ 2 R þ1   2 2R ; minmod : /mm ðRÞ ¼ min Rþ1 Rþ1

In this expression, the solution gradient reconstructed in the previ U I using slope ous section rUI is replaced by the limited gradient r limiting operator /I, thus

Superbee :

dU IJ  /I rU I hIJ :

3.2. Hybrid HLLC shock capturing

ð23Þ

In literature, construction of slope limiters is subjected to the constraint of preserving monotonicity and conservative property. In the context of finite volume method, conservative property implies that the reconstructed gradient at the cell interface has to

ð29Þ

/SB ðRÞ ¼ maxð0; minð1; 2RÞ; minð2; RÞ:

It was concluded in [1] that any Godunov-type scheme built upon a single Riemann solver has drawbacks and occasionally fails due to shock instabilities, especially for very high resolution schemes. Liou [2] analyzed causes of shock instabilities using mass

130

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

fluxes and attributed failures to the dissipative pressure term in the mass fluxes. In order to control the instabilities, artificial dissipations can be appended to the discretization leading to a reduction in the resolution of the scheme. Alternatively, it was suggested that a combination of complementary solvers would be able to resolve the issues. For the HLLC scheme, it was shown in [6,7] that the pressure dissipation term in mass flux is non-zero thus it cannot avoid shock instability. In contrast, the HLLE scheme [5] was shown in [2] that it has the property of zero pressure dissipation term thus HLLE is a shock-stable scheme. The HLLE scheme assumes three constant states separated by two acoustic waves with the numerical flux function computed as,

F IJHLLE

8   > if SLIJ P 0 F U LIJ > > > < L R R L ¼ SR F ðU IJ ÞSL F ðU IJ ÞþSL SR ðU IJ U IJ Þ if SLIJ < 0 6 SRIJ : SR SL > >   > > : F UR if SRIJ < 0 IJ

ð30Þ

The wave speeds were estimated from the local data and the Roe average eigenvalues,

  ~IJ  a ~IJ ; 0 SL ¼ min qLIJ  aLIJ ; q   ~IJ þ a ~IJ ; 0 ; SR ¼ max qRIJ þ aRIJ ; q

ð31Þ

where quantities denoted by tildes are the Roe averages. In this work, we followed suggestions by [1] to employ the HLLE scheme in the vicinity of strong shock areas while the HLLC scheme was used in the rest of the domain. That is,

F IJ ¼ bs F IJHLLC þ ð1  bs ÞF IJHLLE :

ð32Þ

where bs is a switching function to identify cell interfaces (facets) in the vicinity of a strong shock. A simple switching function based on local pressure ratio [1] was used

( bs  bps ¼

1:0; if

jpR pL j minðpL ;pR Þ

< p ;

0:0; otherwise:

ð33Þ

In this expression, p is some threshold parameter indicating strength of shocks which will be resolved using the HLLE scheme. Thought the threshold parameter is problem dependent, in our experience, p = 1.5–2.0 results in robust performance of the hybrid scheme.

3.3.1. Solution sensor Given a solution on the current mesh, it is essential to detect regions in the domain containing high gradients of the solution. The solution gradient was computed by either weighted average or least square methods as described in the earlier section. The underlying idea was that the mesh must be clustered in the area of high gradients and coarsened in the least changed solution. The mesh adaptation is released by determining optimal nodal mesh spacing from a scalar solution r chosen from conserved variables Un. The desired mesh spacing in the domain can be obtained from the principle of equal distribution of errors as suggested in [19]:





d2    d2P  2   dx  e

r

¼ C  constant;

ð34Þ

where dP denotes the local spacing at node P. In this expression, the second order derivative can be obtained by projection or variational recovery. The directional nodal spacing can then be determined as

d2ðiÞP jkiP j ¼ C ¼ d2min jkjmax ;

ð35Þ

where d(i)P is the nodal spacing in direction xi, kiP is the eigenvalues of second order derivatives and dmin is the user-specified minimum spacing. This process requires the recovery of second order derivative and estimation of eigenvalues, thus introduces more overhead in computation. Alternatively, we proposed a simple nodal spacing criteria based on first order derivative,

dP ¼ dmin þ ðdmax  dmin Þe



~ r2 r P 2c2

;

ð36Þ

~ rP measures gradient of where dmax is the maximum spacing and r solution at P as

~ rP ¼ r

jrrP j  jrrjmin : jrrjmax

In the above expression c is a controlling parameter in the range of 0.1–0.5 specifying the width of deleting area. This nodal spacing criteria only requires the evaluation of first order derivatives of chosen solution variables and is much less expensive than the earlier expression (35) involving second order derivatives and eigenvalue computations. 3.3.2. Local mesh adaptivity As the optimal nodal spacing is determined from given solutions, the mesh will be adapted in the areas where there is a significant difference between the desired spacing and current mesh spacing. A node will be marked as deleted if

3.3. Solution-based adaptivity

jdP  dP j 6 d : dP

Applying the above stabilization methods induces deficiencies to the numerical schemes as the limiters and hybrid HLLC approach reduce the accuracy to first order around strong shocks and discontinuities. Especially for unsteady problems, it propagates to the whole domain and pollutes the resolutions. In order to obtain more accurate results, mesh must be clustered in the areas of high gradients. This can be achieved by using either uniformly fine grids or solution adaptive meshes. The former can be very expensive for computation and requires a-priori knowledge of the solutions which makes it unsuitable for unsteady problems. The latter poses a challenge in adaptive mesh generation but more cost effective. In this work, we proposed a local solution adaptivity approach for better resolutions of shock waves in unsteady compressible flows. The approach consists of solution detection scheme to identify the areas required to remesh and a robust meshing procedure to regenerate three dimensional grid locally in those marked areas.

In this expression,  dP is the current mesh spacing at node P and d is the user-specified threshold. In all of the current applications d = 0.3–0.7 provides a reasonable cut-off threshold. Those deleted nodes will then be added to form local deleted regions or holes for remeshing while the deleted mesh is kept as local background meshes of the holes. The empty holes will then be filled by performing volume mesh generation locally. For each of the holes, the boundary faces are identified and clustered as the hole’s surface discretization and the isotropic Delaunay triangulation is applied to generate tetrahedral mesh based on the above local surface mesh and local background mesh. In many cases, when the holes intersect with geometrical surface of the computational domain, it is necessary to regenerate surface mesh on the domain surface. From geometrical definition of the domain surface and the intersecting edge boundary, the advancing front algorithm is employed to generate the local surface mesh. From the newly generated surface mesh, the same procedure of isotropic Delaunay triangulation is

ð37Þ

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

131

Fig. 2. Local mesh adaptivity procedure.

applied to create volume mesh for the hole. Local meshes from the holes is then assembled with the original mesh to form the new volume mesh of the domain. Current solutions then need to be updated on the newly generated mesh. One can refer to [18] for details of the procedure which can be summarised as in Fig. 2. It should be noted that local mesh can be refined or coarsened depending on the local optimal spacings. The solution from the previous mesh is interpolated to the adapted mesh by a weighted linear method. Higher order interpolations can also be implemented to improve solution reconstruction.

4. Numerical examples 4.1. Sod’s shock tube problem We first tested various implemented shock capturing schemes on the one dimensional Sod’s shock tube problem [4]. The exact solution to this popular Riemann problem can be found and consists of a shock wave, right travelling contact wave and a left rarefaction wave. In our computation, a one dimensional set-up was cast in a three dimensional box of 1 unit length and 0.1 unit square cross section. The diaphragm position is located at x = 0.3 and symmetry boundary conditions are applied at the side faces to resemble the one-dimensional set-up. The domain was discretized into an unstructured mesh of uniform size with 200 points along the length of the channel as shown in Fig. 5. The flow was initialized as a jump in conserved quantities across the diaphragm. Fig. 3 show density and pressure profiles along the centerline of the domain at time T = 0.2 using different shock capturing schemes. It can be seen that the second order HLLC flux (HLLC2) significantly improves the accuracy over the first order scheme. Among the second order schemes, the hybrid scheme and the LED limiter gives similar resolution across the shock and

contact discontinuities, even though the shock profile shows a small kink. The LED limiter with a clipping factor helps to smoothen out shock profiles while maintaining the resolution across rarefaction waves. In Fig. 3, it also shows the results of different types of limiters for this problem with better shock capturing using hybrid scheme and LED with auto clipping. Next we employed mesh adaptation with second order HLLC (HLLC2) and LED limiter for this problem where mesh features were locally changed according to solution propagation. In this case, the finest mesh size was taken as half of the size of the uniform mesh, and the largest mesh size is double of the size in uniform mesh. This resulted in the adapted mesh having at most the same number of elements as many as the uniform one. It can be clearly seen from Fig. 4 that mesh adaptation is able to improve resolution across the rarefaction and contact discontinuity while it performed better in capturing shocks than those methods operating on the fixed uniform mesh. Fig. 5 shows the mesh clustering around the shock, contact and rarefaction wave while the mesh is much coarser outside of the area where the solution is least changing. 4.2. The forward-facing step problem We consider the problem of Mach 3 flow over a forward facing step channel studied extensively in [20] and later by many others in literature for different shock capturing schemes. In this work, we studied the problem in three dimensions with the channel of 1 unit width, 3 units length and 1 unit depth containing a step of 0.2 unit length and located at 0.6 unit length. Reflecting boundary conditions are applied along top, bottom and side boundaries; while inflow and outflow boundary condition are set at inlet and outlet, respectively. The corner of the step is a singular points thus affecting the solution around the corner. It is well-known that this may

132

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

Fig. 3. Comparison of (a) density and (b) pressure profile for shock tube problem at T = 0.2 using different shock capturing schemes. Second order HLLC gives better resolution while the LED with clipping and superbee limiter show improved shock capturing.

Fig. 4. Comparison of (a) density and (b) pressure profiles for shock tube problem at T = 0.2 using local adaptation with different limiters on uniform mesh.

Fig. 5. Results of pressure, density contours and mesh used for Sod’s shock tube problem. Comparison of contours for (a) fixed mesh with HLLC/HLLE hybrid approach and (b) adapted mesh at T = 0.2 with second order HLLC and LED limiter.

cause an entropy layer downstream of the corner in the bottom wall. This layer, when interacting with reflected shock waves potentially creates a spurious Mach stem at the bottom wall.

In this example, the second order HLLC scheme with LED limiter is employed for solving the flow over the step. Fig. 6 shows meshes used for the current computation including a unstructured uniform

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

133

Fig. 6. (a) Uniform mesh and (b) adapted mesh at t = 3.0 with second order HLLC and LED limiter.

Fig. 7. Mach number profile of flow over a forward facing step with (a) uniform mesh and (b) adaptivity with second order HLLC and LED limiter.

mesh of size d0 = 0.02 and an adapted mesh of the smallest size equal to half of the uniform size dmin = 0.5d0. However, in the adapted mesh, the mesh is only clustered at the area of high gradient while it is much coarser elsewhere, thus resulting in about the same total number of elements. In Fig. 7, Mach number profiles at time t = 3.0 are shown using both shock capturing strategies. It can be seen that the adaptivity is able to capture much better shock fronts compared to the solution on the uniform mesh. In our adaptivity scheme, the mesh is adapted following the propagation of the shocks, shown in Fig. 8 to better resolve unsteady features of the flows. Fig. 9 shows contours of density and entropy at t = 0.3

using different shock capturing strategies. Finer and well resolved shock fronts were again observed for the adaptivity scheme using HLLC with LED limiter. The entropy layer at the bottom of the domain is less severe in the adapted mesh than uniform one. 4.3. The double Mach reflection problem Next, we consider a problem of Mach 10 normal shock wave travelling over a ramp. This is another typical example taken from [20] and commonly used for testing of shock capturing schemes. In this work, a 2D plane is extruded in the z direction to make a full

Fig. 8. Adapted mesh follows shock propagation in the flows.

Fig. 9. Density and entropy (A = p/qc) contours of flow over a forward facing step using (a) HLLC/HLLE hybrid approach and (b) adaptivity with second order HLLC and LED limiter.

134

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

Fig. 10. Meshes used for computation of flow over a ramp: (a) uniform mesh approach and (b) adapted mesh at t = 0.2 with second order HLLC and LED limiter.

3D channel with a ramp of 30°. Reflected boundary conditions are applied at the top and bottom wall, while the z-plane is forced with symmetric boundary conditions. The exact plane shock of Mach 10 is prescribed at the inlet. The Mach 10 shock wave moving in x direction interacts with a ramp to form a double Mach reflection including two Mach stems and two contact discontinuities. The challenge for this problem lies in the ability to capture the second weaker contact and shock accurately as well as a jet formed at the

lower wall due to the interaction of contact discontinuity and the wall. Results from using different shock capturing schemes for this problem are shown in Figs. 10 and 11 for computational grids and density profiles. It can be seen again that the use of adaptivity with second order HLLC and LED limiter is able to resolve the shock better with apparent flow structures compared to results on fine uniform mesh. The numerical shock instabilities are also avoided by using HLLC2 with LED limiter in those cases.

Fig. 11. Density contour of flow over a ramp using (a) hybrid approach on uniform mesh and (b) adaptivity with second order HLLC and LED limiter.

Fig. 12. (a) Pressure history at three stations d = 1, 2, 3 m away from the charge and (b) comparison of peak pressure using different schemes. The adaptivity with HLLC2 and hybrid improve peak pressure results.

V.-T. Nguyen et al. / Computers & Fluids 70 (2012) 126–135

135

Fig. 13. Comparison of impulse and arrival time using different schemes against empirical data.

4.4. Explosion in a box

References

The air blast resulting from a 1.0 kg charge of TNT is simulated to test the ability of the proposed approach to capture strong shocks associated with explosive detonations. The explosion is modelled using a bursting sphere approach where the volume of the explosive is initialized with an ideal gas having a high density and energy based on the detonation energy of the explosive material. To compare with standard Ref. [22], the charge size is assumed to be 1.0 kg TNT which has a radius of 5.3 cm. The computational domain has dimensions of 5.5  2.5  2.5 m and utilises two planes of symmetry. Monitoring points inside the domain record the shock pressures over a distance of 5 m along the x-axis (symmetry). The baseline mesh is generated with a constant element size of 2 cm along the x-axis (symmetry) to a radius of 0.5 m, and a background mesh size of 10 cm. A cell size of 1 cm was used for the explosive charge to improve the initial solution gradients. Fig. 12a shows the comparison of peak pressure at various stations of d = 1, 2, 3 m away from the charge. Second order HLLC scheme (HLLC2) improves the peak pressure as well as the pressure profile in all stations. It is also apparent that local mesh adaptation helps to improve the solutions. The peak pressure was highest using local mesh adaptation with HLLC2 hybrid scheme. In Fig. 12b, the peak pressure was compared against the empirical data obtained from [22]. It can be seen that the current result are very close to the empirical ones especially for the results from HLLC2 hybrid scheme with mesh adaptivity. The impulse and arrival time are shown in Fig. 13 where the results from current work shows good agreement with the empirical data.

[1] Quirk JJ. A contribution to the great Riemann solver debate. Int J Numer Meth Fluid 1994;18:555–74. [2] Liou MS. Mass flux schemes and connection to shock instability. J Comput Phys 2000;160:623–48. [3] Harten A, Lax PD, van Leer B. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev 1983;25:35–61. [4] Toro EF. Riemann solvers and numerical methods for fluid dynamics. 2nd ed. Springer; 1999. [5] Einfeldt B, Munz CD, Roe PL, Sjogreen B. On Godunov-type methods near low density. J Comput Phys 1991;92:273–95. [6] Kim Sung Don, Lee Bok Jik, Lee Hyoung Jin, Jeung In-Seuck. Robust HLLC Riemann solver with weighted average flux scheme for strong shock. J Comput Phys 2009;228:7634–42. [7] Kim Sung Don, Lee Bok Jik, Lee Hyoung Jin, Jeung In-Seuck, Choi Jeong-Yeol. Realization of contact resolving approximate Riemann solvers for strong shock and expansion flows. Int J Num Meth Fluids 2010;62:1107–33. [8] Remaki L, Hassan O, Morgan K. New limiter and gradient reconstruction method for HLLC finite volume scheme to solver Navier–Stokes equations. In: V European conference on computational fluid dynamics ECCOMAS CFD; 2010 [9] Hassan O, Sorensen KA, Morgan K, Weatherill NP. A method for time accurate turbulent compressible fluid flow simulation with moving boundary components employing local remeshing. Int J Numer Meth Fluids 2007;53:1243–66. [10] Roe PL. Approximate Riemann solvers, parameter vectors, and difference schemes. J Comput Phys 1981;43:357. [11] Osher S, Solomon F. Upwind difference schemes for hyperbolic systems of conservation laws. Math Comp 1982;38:339. [12] Von Neumann J, Richmyer RD. A method for numerical calculation of hydrodynamic shocks. J Appl Phys 1950;21:232–7. [13] Cockburn B, Shu C-W. Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, review article. J Sci Comp 2001;16:173–261. [14] Mavriplis DJ. Revisiting the least squares procedure for gradient reconstruction on unstructured meshes. AIAA paper; 2003. p. 2003–3986. [15] Batten P, Clarke N, Lambert C, Causon DM. On the choice of wavespeeds for the HLLC Riemann solver. SIAM J Sci Comput 1997;18:15531570. [16] Barth TJ, Jespersen DC. The design and application of upwind schemes on unstructured meshes. In: 27th Aerospace sciences meeting. AIAA paper 890366; January 1989. [17] Berger M, Aftosmis MJ. Analysis of slope limiters on irregular grids. In: 43rd AIAA aerospace sciences meeting. AIAA paper 2005-0490. January 2005. [18] Nguyen HH, Nguyen V-T, Hassan O. Local solutioned-based mesh adaptation for unstructured grids in complex three-dimensional domains. In: The 21st international meshing roundtable: mesh modelling for simulations and visualization, San Jose, California, USA; October 2012. [19] Peraire J, Vahdati M, Morgan K, Zienkiewicz OC. Adaptive remeshing for compressible flow computations. J Comput Phys 1987;17:449–66. [20] Colella P, Woodward P. The piecewise parabolic method (PPM) for gasdynamical simulations. J Comput Phys 1984;54:174–201. [21] Darwish MS, Moukalled F. TVD schemes for unstructured grids. Int J Heat Mass Transfer 2003;46:599611. [22] Kingery CN, Bulmash G. Airblast parameters from TNT spherical air burst and hemispherical surface bursts. ARBRL-TR-02555; April 1984.

5. Conclusions In this work, various shock capturing schemes were implemented and validated in the framework of second order finite volume methods on unstructured grid in three dimensions. The implemented schemes were tested for different problems including standard shock tube case as well as very strong shock examples, particularly explosion cases. The local mesh adaptation significantly improved resolution of solutions in the presence of shocks. The numerical examples showed good performance in shock capturing as compared against available literature and empirical data.