Development of an improved spatial reconstruction technique for the HLL method and its applications

(1)

Development of an improved spatial reconstruction technique

for the HLL method and its applications

Matthew R. Smith

a

, K.-M. Lin

b

, C.-T. Hung

b

, Y.-S. Chen

c

, J.-S. Wu

b,⇑ a

National Centre for High-Performance Computing, National Applied Research Laboratories, No. 7 R&D Rd. VI, Hsinchu Science Park, Taiwan b

Department of Mechanical Engineering, National Chiao Tung University, 1001 Ta-Hsueh Road, Hsinchu 30050, Taiwan c

National Space Organization, National Applied Research Laboratories, 8F, 9 Zhan-Ye 1st Road, Hsinchu Science Park, Hsinchu, Taiwan

a r t i c l e

i n f o

Article history:

Received 27 February 2010

Received in revised form 24 August 2010 Accepted 17 September 2010

Available online 8 October 2010 Keywords:

Computational ﬂuid dynamics (CFD) Finite volume methods (FVM)

Total variable diminishing (TVD) schemes Plasma simulation

a b s t r a c t

The integral form of the conventional HLL fluxes are presented by taking integrals around the control volume centred on each cell interface. These integrals are demonstrated to reduce to the conventional HLL flux through simplification by assuming spatially constant conserved properties. The integral flux expressions are then modified by permitting the analytical inclusion of spatially linearly varying conserved quantities. The newly obtained fluxes (which are named HLLG fluxes for clarification, where G stands for gradient inclu-sion) demonstrate that conventional reconstructions at cell interfaces are invalid and can produce unstable results when applied to conventional HLL schemes. The HLLG method is then applied to the solution of the Euler Equations and Shallow Water Equations for var-ious common benchmark problems and finally applied to a 1D fluid modeling for an argon RF discharge at low pressure. Results show that the correct inclusion of flow gradients is shown to demonstrate superior transient behavior when compared to the existing HLL sol-ver and conventional spatial reconstruction without significantly increasing computational expense.

1. Introduction

The finite volume method for solution to partial differential equations forms the mainstay of modern computational fluid dynamics (CFD). The development of a conservative scheme for the solution of non-linear systems of hyperbolic conserva-tion laws by Godunov[1]introduced a revelation in the approach to the computation of fluid flow. The resulting family of solvers (known as Godunov type methods) discretizes the flow region into computational regions known as cells and use solutions to the Riemann problem at cell interfaces to calculate fluxes of conserved quantities in and out of control volumes. The original HLL method was developed by Harten et al.[2]as an approximate Riemann solver for use in a Godunov sol-ver. Rather than solving the Riemann problem analytically with knowledge of the behavior of the system (in many cases, an ideal gas), the HLL method solves for the flux in an intermediate region (or star region) between cells directly from the gov-erning partial differential equations. A control volume centred on a cell interface is integrated over space and time assuming spatially and temporally constant properties. Through the introduction of the star region bounded by two propagating waves, the flux across the interface can be mathematically determined. This allows for the flexible solution of a various number of conservative hyperbolic systems, such as the Shallow Water Equations and the Euler Equations.

The keystone for the HLL method is the assumption of the presence of a single intermediate state bounded to two prop-agating waves. This is a valid assumption for two variable hyperbolic systems[3]and thus the HLL method is still used for

⇑Corresponding author. Tel.: +886 3 573 1693; fax: +886 3 611 0023. E-mail address:[email protected](J.-S. Wu).

Contents lists available atScienceDirect

Journal of Computational Physics

(2)

the solution of the Shallow Water Equations by various authors[3–5]. However, for solutions to other sets of governing equa-tions, such as the Euler Equaequa-tions, this assumption is invalid[3]due to the presence of two or more intermediate states, each separated by a distinct propagating wave. One solution to this problem was proposed by Toro et al.[6]by the introduction of a contact discontinuity in the HLL expressions. The HLLC (where C stands for Contact surface) provides results with less numerical dissipation.

A conventional approach for the extension of the method to higher order spatial accuracy is performed through the inclu-sion of spatially varying conserved quantities[3,4]. While this does not solve the fundamental problems associated with HLL, it does aid in the removal of numerical dissipation. In such an approach, gradients are used to reconstruct conserved quan-tities at cell interface locations, following which the ordinary ﬂux expressions are employed for the solution to the Riemann problem.

Presented here is a modification to the HLL method where allowance is made for the inclusion of gradient terms within the flux expressions themselves. First, the original HLL fluxes are presented in integral form and shown to reduce to the ori-ginal HLL expressions under the assumption of spatially constant properties. The integral form is then re-evaluated when conserved quantities are permitted to vary linearly in space. The resulting fluxes are then combined with the traditional spa-tial reconstruction concept to produce the HLLG scheme, where the G represents the inclusion of gradient terms. This meth-od is then applied to the solution of the Shallow Water Equations, then the Euler Equations and finally the fluid mmeth-odeling equations for RF discharge at low pressure. The results show an improvement in the numerical diffusion associated with the traditional HLL higher order extensions without significantly increasing computational expense.

2. Numerical method

2.1. Integral form of the original HLL ﬂux expressions

The ﬂuxes developed by Harten et al.[2]are presented here in their complete integral form.Fig. 1shows a control volume in x t space covering the region between cells i and i + 1 centred on the interface separating the cells at x = 0. The region is temporally bound by the limits t = 0 and t = T. At t = 0, waves moving at velocities SL(<0) and SR(>0) move away from the

dis-continuity between the cell interface. If the governing differential equation is given by

@U @t þ

@FðUÞ

@x ¼ 0; ð1Þ

the conditions inside the control volume in the region between [0, T] and [xL, xR] can be described by the integral: Z0 xL U0Ldx þ Z xR 0 U0Rdx þ ZT 0 FLdt Z T 0 FRdt ¼ Z SLT xL UTLdx þ Z SRT SLT UTdx þ Z xR SRT UTRdx; ð2Þ

where the subscripts L and R represent conditions in the left and right cells, respectively, while the superscript 0 and T rep-resents conditions at time t = 0 and t = T, respectively. The subscript * reprep-resents conditions in the star region between the propagating waves. Rearranging to solve in terms of this star region, we obtain:

ZSRT SLT UT dx ¼ Z 0 xL U0 Ldx þ Z xR 0 U0 Rdx þ Z T 0 FLdt Z T 0 FRdt Z SLT xL UT Ldx Z xR SRT UT Rdx: ð3Þ

This equation assumes nothing regarding variation of ﬂuxes F or conserved quantities U within space and time. The only assumptions made are in the presence of the two propagating waves surrounding a single intermediate region at x = 0.Fig. 2

(3)

shows the revised x t diagram focusing on the left cell only. Using the same method of integrating conserved quantities over space and ﬂuxed quantities over time obtains:

Z 0 xL U0 Ldx þ Z T 0 FLdt Z T 0 F¼ Z SLT xL UT Ldx þ Z 0 SLT UT dx: ð4Þ

We further assume that the average state over the region between x = SLT and x = SRT is the same as the average between

the region x = SLT and x = 0, or: 1 SLT Z 0 SLT UTdx ¼ 1 SRT SLT Z SRT SLT UTdx; ð5Þ or: Z 0 SLT UTdx ¼ SL SR SL Z SRT SLT UTdx: ð6Þ

This assumption provides us the knowledge required to solve for our desired quantity, the ﬂux F*. Combining Eqs.3,4,6

provide the complete integral form of the total ﬂux F* over time step T:

Z T 0 F¼ Z 0 xL U0Ldx ZSLT xL UTLdx þ Z T 0 FLdt Z 0 SLT UTdx ¼ Z 0 xL U0 Ldx ZSLT xL UT Ldx þ Z T 0 FLdt þ SL SR SL Z 0 xL U0Ldx Z SLT xL UTLdx þ Z xR 0 U0Rdx Z xR SRT UTRdx þ Z T 0 FLdt Z T 0 FRdt : ð7Þ

2.2. Simpliﬁcation of the integral form

The integral form of the HLL ﬂuxes can be evaluated by making assumptions about the integral expressions. Assuming the regions to the left of the left moving wave and to the right of the right moving wave remain untouched over the period of a time step T, the following simpliﬁcations can be made:

Z 0 xL U0Ldx Z SLT xL UTLdx Z 0 SLT ULdx; Z xR 0 U0Rdx Z xR SRT UTRdx ZSRT 0 URdx: ð8Þ

Note these assumptions still permit spatially varying quantities of U and only assume that regions outside propagating waves SLand SR. From a (more mathematical) perspective, this step represents a simple resizing of the control volume

em-ployed for the flux computation. Using these assumptions, the modified equation simplifies to:

Z T 0 F¼ Z 0 SLT ULdx þ Z T 0 FLdt þ SL SR SL Z 0 SLT ULdx þ Z SRT 0 URdx þ Z T 0 FLdt Z T 0 FRdt : ð9Þ

(4)

Notice now that the equation for the total flux across the interface is now only governed by integrals bounded by the left and right moving waves. This control volume is shown together with the remaining cell topography inFig. 3. By assuming that the quantities U remain spatially constant (i.e. monotonic in nature) and the fluxes are temporally constant we find (for example): Z 0 SLT ULdx ¼ SLTUL; Z T 0 FLdt ¼ FLT: ð10Þ

The consequence of the latter assumption is that the ﬂux FLis not permitted to change over time (or alternatively, FL

rep-resents the average value of F at that location) and must be calculated accordingly. After applying the said simpliﬁcations and evaluating the complete integral expression in Eq.(7), the original HLL ﬂux expressions[2]are obtained.

2.3. Inclusion of spatial variation of conserved quantities

The assumption of constant conserved quantities U made in Eq.(10)results in spatially ﬁrst order accuracy. We can in-stead allow U to vary linearly in space such that:

ULðxÞ ¼ Uiþ 1 2 L þ dU dx L x; ð11Þ

where superscript i + 1/2 indicates conditions at the cell interface where x = 0 and subscript L refers to conditions on the left side of the interface. This is demonstrated inFig. 4. The gradient term is calculated using ﬁnite differences together with a selected slope limiter to maintain positivity. A commonly used example is the MINMOD limiter[7]:

dU dx eff ¼ MINMOD dU dx F ; dU dx B ; ð12Þ

Fig. 3. Revised control volume between propagating waves moving with speeds SLand SR. The ﬂuxes through the left and right surfaces of the control volume are calculated using values of U at locations SRTand SLTthus allowing the assumption of constant FLand FRover time T.

Fig. 4. Spatial reconstruction of conserved quantity U. The coordinate system must be consistent with that used in the HLL ﬂux derivation. The revised control volume (shown) is bounded by SLDt < x < SRDt.

(5)

where ‘‘eff” indicates an effective gradient, subscripts F and B are the forward and backward differences, respectively. The MINMOD (or Minimum Modulus) function is given by:

MINMOD½a; b ¼

0; IF SIGNðabÞ < 0;

a; IF ðSIGNðabÞ > 0Þ AND ðjaj < jbjÞ; b; IF ðSIGNðabÞ > 0Þ AND ðjbj < jajÞ: 8

> < >

: ð13Þ

A general improvement to the MINMOD limiter is the Monotonized Central Difference (MC) limiter[8]. This has been used in the previous applications of HLL applied to solutions of the Euler and Shallow Water equations[3,4]in addition to the kinetic theory based Euler solver QDS (Quiet Direct Simulation)[9]. The MC limiter is deﬁned as:

dU dx eff ¼ MINMOD dU dx C ;

a

MINMOD dU dx F ; dU dx B : ð14Þ

The variable alpha is used to control dissipation in the scheme and is set to

a

= 1 or 2 to limit dissipation to a minimum (setting

a

= 0 recreates a spatially ﬁrst order scheme). Assuming a uniform cell widthDx, the interface values can be found using the gradient information such that:

Uiþ12 L ¼ ULþ dU dx L

D

x 2 ! U iþ1 2 R ¼ UR dU dx R

D

x 2 : ð15Þ

Substituting these obtains the simple expressions:

ULðxÞ ¼ ULþ dU dx L x þ

D

x 2 ; URðxÞ ¼ URþ dU dx R x

D

x 2 : ð16Þ

This allows an alternate evaluation of the integral form of the HLL flux expressions. Substituting and evaluating the inte-gral obtains the final flux:

F 1 SR SL F ULþ dU dx L

D

x 2 þ SLT SR F UR dU dx R

D

x 2 SRT SL SRSL UR dU dx R

D

x 2 SRT 2 ULþ dU dx L

D

x 2 þ SLT 2 : ð17Þ

The average ﬂux represented in Eq.(17)can be rewritten in terms of the original HLL ﬂux expression:

F F U1L SR F U1R SL SRSL U2R U 2 L SR SL ; ð18Þ

where superscripts 1 and 2 indicate reconstruction at locations at distances ST and ST/2 away from the cell interface respec-tively. Care must be taken with the reconstruction when SL> 0 or SR< 0 to avoid over-reconstruction in regions each

respec-tive cell when upwind ﬂuxes are used.

2.4. Comparison against conventional higher order extension

The original HLL scheme assumes monotonic conditions across each cell and is thus ﬁrst order accurate in space. In the decades which have followed since HLL’s developments, a large family of methods for extension of spatial accuracy have been developed. Many of these methods use cell-averaged states to compute gradients for use in reconstruction at cell inter-faces[10–16]. For example, recent developments such as General Riemann Problem (GRP) schemes[14–16]employ spatial reconstruction to determine conditions at the interface:

Un_jþ1=2;þ¼ Unj þ dU dx n j

D

x 2 ; U n jþ1=2;¼ U n jþ1 dU dx n jþ1

D

x 2 ; ð19Þ

where superscript n indicates conditions at step n and +/ represents the left/right side of the interface at j + 1/2, respec-tively. Following this reconstruction, GRP schemes differ only from conventional second order methods through employing analytical expressions for dU

dt

n

jþ1=2to advance the solution[15]: Unþ1=2_jþ1=2 ¼ Un_jþ1=2þ dU dt n jþ1=2

D

t 2; ð20Þ where Un jþ1=2¼ R A _Un jþ1=2;þ;U n jþ1=2;

is the solution to the Riemann problem at step n. These conditions are then used to com-pute the ﬂuxes across the interface. The inclusion of this step provides the extension to higher order temporal accuracy with-out relying on more approximate temporal extensions as often employed in conventional MUSCL schemes [17]. In the

(6)

absence of this additional step, such schemes are identical to (temporally first order accurate) MUSCL type extensions[15]. The integral (balance) approach demonstrated here is not applicable to the additional step employed by GRP type schemes since investigations presented here are limited to first order temporal accuracy (refer to integrals given in Eq.(10)). Follow-ing conventional reconstruction at cell interfaces and application of the original HLL fluxes upon reconstructed states, the flux is: F F ULþ dUdx L Dx 2 SR F UR dUdx R Dx 2 SL SRSL UR dUdx R Dx 2 ULþ dUdx L Dx 2 SR SL : ð21Þ

Although a trivial exercise in mathematics, comparison against Eqs.(17) and (21)reveals that conventional spatial recon-struction (as employed by MUSCL and/or GRP schemes) is inconsistent when HLL is used for flux computation. The formal inclusion of gradients in the HLL flux derivation reveals that two reconstructions are required, at distances ST and 0.5ST on each side from the cell interface respectively, for the flux function evaluation and the U term typically associated with the anti-diffusion component. This might seem trivial given the demonstrated performance of these schemes[3,4]. Indeed, when using most popular limiting functions (e.g. MINMOD) the flaws associated with this oversight are very difficult to perceive for some problems. However, demonstration of flaws resulting from this oversight is possible.

2.5. Implementation of HLLG ﬂuxes

The calculation of the modiﬁed HLL (or HLLG) ﬂuxes in one dimension is shown by the following steps and supplemented by the diagram shown inFig. 5:

1. Estimate the wave speeds SLand SRusing the cell-centred values. The value of these wave speeds differ depending on the

system of equations solved. Review into the computation of wave speed estimates lies outside the scope of this investi-gation – details can be found in[3–5].

2. Reconstruct U at the locations within the control volume shown inFig. 5:

f1 L ¼ min

D

x 2 þ SLT;

D

x 2 ; f1 R¼ max

D

x 2 þ SRT;

D

x 2 ; f2 L ¼ min

D

x 2 þ SLT 2 ;

D

x 2 ; f2 R¼ max

D

x 2 þ SRT 2 ;

D

x 2 ; U1L¼ ULþ fL1 dUL dx eff ; U1R¼ UR fR1 dUR dx eff ; U2 L¼ ULþ fL2 dUL dx eff ; U2 R¼ UR fR2 dUR dx eff : ð22Þ

3. Calculate the ﬂux across the interface. If either SL> 0 or SR < 0, the respective upwind ﬂuxes are used

SL>0 : F¼ F U1L

;SR<0 : F¼ F U1R

. Otherwise, the revised HLL (or HLLG) ﬂux is used:

Fig. 5. Reconstruction locations for the revised (HLLG) expressions showing the control volume employed during the analysis. Fluxes F are computed using values of U reconstructed at the edge of the control volume while the value of conserved quantities U employed in the HLLG ﬂux expression are reconstructed in the centre of each half of the control volume.

(7)

F F U1 L SR SR F U1R SL SRSL U2R U 2 L SR SL : ð23Þ

3. Results and discussion

3.1. 2D shallow water equations – dam break

An ideal two dimensional circular dam break problem is considered. The two dimensional dam break simulation is per-formed to test the shock capturing properties of HLLG. The governing equations employed are the Shallow Water Equations, deﬁned as: @U @tþ @F @xþ @G @y¼ 0; @ @t / /u /

v

2 6 4 3 7 5 þ_@x@ /u /u2_þ/2 2 /u

v

2 6 4 3 7 5 þ_@y@ /

v

/u

v

/

v

2_þ/2 2 2 6 4 3 7 5 ¼ 0; ð24Þ

where / = gh, g is gravitation acceleration (9.81) and h is the water level. The values u and

v

are the velocity in the x and y directions, respectively. A simple method for estimating the wave speeds is SL= MIN[VL aL, VR aR] and SR= MAX[VL+ aL, V-R+ aR][3–5]where aL¼ ffiffiffiffiffi /L p and aR¼ ffiffiffiffiffiffi /R p

are the characteristic speeds of propagation and VL, VRare the velocities of ﬂow

normal to the interface. A circular region of water of elevation hINis elevated and separated from the surrounding water of

elevation hOUTby an inﬁnitely thin separating dam. The initial water levels are: h0ðxÞ ¼ hIN¼ 2:5 m; ðx xcÞ2þ ðy ycÞ 2 6_R; hOUT¼ 0:5 m; ðx xcÞ2þ ðy ycÞ 2 >R; ( ð25Þ

where R is the radius of the circular region R = 2.5 m and subscript c indicates the location of the centre of the elevated region (xc, yc). The water is assumed to be at rest (u =

v

= 0). The computational domain extends over L = H = 40 m and results are

presented for a regular Cartesian computational grid of 200 by 200 cells after a ﬂow time of 5 s. The time step is adaptive to ensure a constant maximum CFL of 0.5. The results are expected to be radially symmetric in nature due to the radially symmetric nature of the initial conditions.

Various figures of the resulting water level after 5 s are shown inFig. 6. The oscillations shown in the HLL results strongly affect the flow field. The oscillations are weaker along the horizontal (h = 0°) and vertical (h = 90°) lines. This may mislead one to interpret the errors as due to poor mesh alignment. Indeed, employing a circular grid will remove most of these oscilla-tions. However, these errors are not characteristic of errors resulting from direction decoupling[18]and are not consistent with the HLLG results, as also demonstrated by the water levels shown inFig. 6. There are major differences between the diagonal and horizontal results shown in the HLLG solutions, thus ruling out direction coupling and poor grid alignment as the cause of the oscillations present in the HLL results. The spatial reconstruction of conserved properties is performed using finite differences with the MINMOD limiter[7], which is known to be diffusive. Thus, the same simulations are re-peated using the Monotonized Central Difference (MC) limiter[8]with

a

= 2. The influence of the instabilities presented in the HLL solution is dominant, as demonstrated inFig. 7. However, the influence of the new limiting function has only min-imal influence on the HLLG results, which still remain radially symmetric with only minor oscillations in regions where flows are unaligned with the computational grid. The above shows that the use of HLLG scheme greatly improves the water shal-low equation simulation and reduces sensitivity to the selected limiter when compared to the use of conventional HLL reconstruction.

3.2. 1D Euler equations – shock/acoustic wave interaction

The one dimensional Euler Equations are applied to Shu and Osher’s simulation of a shock/acoustic wave interaction

[19,20]. The one dimensional Euler Equations, without source terms in conservation form, are:

@U @tþ @F @x¼ 0; @ @t

q

u E 2 6 4 3 7 5 þ_@x@

q

u

q

u2_{þ p} uðE þ pÞ 2 6 4 3 7 5 ¼ 0; ð26Þ

where

q

is the gas density, u is the velocity, p is the pressure computed from the ideal gas law (p =

q

RT) and the energy E is given by E =

q

(1/2u2_{+ C}

VT) where Cv is the speciﬁc heat constant at constant volume and T is the temperature. The initial conditions

(8)

4 6 x < 5. The gas is inviscid and ideal (

c

= 1.4) and the flow is allowed to develop until flow time t = 1.8. Wave speed com-putation is similar to that employed above except that the characteristic propagation speed is now the local speed of sound Fig. 6. Values of geopotential (U) taken from solutions to the 2D dam break problem as computed by HLL (left) and HLLG (right). Each simulation employed 200 200 cells and a constant maximum CFL of 0.5. Gradients employed for reconstruction are calculated by finite differences and limited by the MINMOD limiter. Contours are displayed at levels of 0.5.

(9)

a = (

c

RT)1/2_{. Resulting density contours from both the proposed HLLG solver and HLL solver are presented in}_{Fig. 8}_{. Each}

sim-ulation uses 400 cells with a constant maximum CFL of 0.1 enforced. The MC limiter (

a

= 2) is employed by both solutions to Fig. 7. Values of geopotential (U) taken from solutions to the 2D dam break problem as computed by HLL (left) and HLLG (right). Each simulation employed 200 200 cells and a constant maximum CFL of 0.5. Gradients employed for reconstruction are calculated by ﬁnite differences and limited by the MC limiter. Contours are displayed at levels of 0.5.

(10)

limit the spatial reconstruction to ﬁrst order in spurious regions. It can be seen that the result obtained by HLLG provides a clo-ser comparison to the benchmark solution shown than the conventionally reconstructed second order HLL results.

3.3. 2D Euler equations – blast wave

The two dimensional Euler Equations are applied to the simulation of a temperature driven radially symmetric blast wave

[18]. The two dimensional Euler Equations are:

@U @t þ @F @xþ @G @y¼ 0; @ @t

q

u

q

v

E 2 6 6 6 4 3 7 7 7 5þ @ @x

q

u

q

u2_{þ p}

q

u

v

uðE þ pÞ 2 6 6 6 4 3 7 7 7 5þ @ @y

q

v

q

v

u

q

v

2_{þ p}

v

ðE þ pÞ 2 6 6 6 4 3 7 7 7 5¼ 0; ð27Þ

where

q

is the gas density, u and

v

are velocities in the x and y directions, respectively, and p is the pressure computed from the ideal gas law (p =

q

RT). The internal energy E is given by E =

q

(1/2(u2₊

_v

2_{) + C}

VT). A region of high temperature gas of

radius r = 0.1L is separated by an imaginary diaphragm and surrounded by a region of low temperature gas THT1L ¼ 100

. At t = 0, the diaphragm is removed and the resulting blast wave propagates outwards in a radially symmetric fashion. The resulting contours of density and velocity are shown inFig. 9. Both solvers are shown to possess characteristic direction decoupling errors with higher regions of density along the diagonal[18]. The conventional extension of HLL to sec-Fig. 8. One dimensional shock/acoustic wave interaction problem. Results show proﬁles of density for the conventionally extended HLL solver and the proposed HLLG extension. All simulations employ 400 cells with an enforced maximum CFL = 0.1. HLL and HLLG simulations using 2000 cells are employed as benchmark results.

(11)

ond order spatial accuracy using the MC limiter results in a large number of asymmetric, spurious oscillations which are not present in the HLLG results.

Fig. 9. Two dimensional temperature driven blast wave problem contours of density [top], x velocity [center] and y velocity [bottom] from proposed HLLG reconstruction [right] and conventional HLL reconstruction [left]. Gradients employed for reconstruction are calculated by ﬁnite differences and limited by the MC limiter (a= 0). Both simulations employ 100 100 cells with an enforced maximum CFL number of 0.5. Results are shown after ﬂow time t = 0.25.

(12)

3.4. 2D Euler equations – supersonic ﬂow over a forward facing step

The proposed HLLG method is applied to the simulation of supersonic ﬂow over a forward facing step[21]. The problem geometry and boundary conditions are shown inFig. 10. An ideal gas (

c

= 1.4) with an initial Mach number of M = 3 every-where. This problem is commonly used to verify gas flow solvers. The density contours computed by both methods are pre-sented inFig. 11after the flow has been advanced to flow time t ¼ 4:0. In this instance, the differences between the HLLG and the conventionally extended HLL solution are negligible.

Fig. 10. Supersonic flow over a forward facing step problem. All surfaces are reflective with exception of the supersonic outflow located at x = 3, 0.2 6 y 6 1.0 and the inflow located at x = 0.0.

Fig. 11. Supersonic ﬂow over a forward facing step problem after t = 4 showing contours of density computed by conventionally extended second order HLL [top] and the proposed HLLG second order extension [bottom]. Contours of density are shown every 0.25. The simulation employs 300 100 cells with CFL = 0.25.

(13)

3.5. 2D Euler equations – shock bubble interaction

We apply the proposed reconstruction method to the simulation of a shock wave moving over a low density bubble[19]. The shock-bubble problem solved is shown inFig. 12. The gas is assumed ideal (

c

= 1.4). Symmetry can be exploited along y = 0 and only half the problem need be considered. The flow is advanced until flow time t ¼ 0:2. All simulations employed 200 200 cells with varying CFL numbers (0.5–0.75). Presented inFig. 13are numerical Schlierens (indicative of gradients of density) demonstrating similar oscillatory instabilities in regions where flows are misaligned with the computational grid. While the HLLG results demonstrate these errors (which are dependent on the CFL number employed), the numerical ‘‘wrig-gles” present in the HLL solution are much more severe. Thus, for any given CFL number, the result provided by the HLLG reconstruction will prove less-prone to numerical instability than the original extension of HLL to higher order accuracy. 3.6. 1D argon RF plasma simulation

Finally, the simulation of parallel-plate argon RF (f = 13.56 MHz) coupled capacitively discharge (CPC) is investigated using the HLLG solver. The motion of charged and neutral species is governed by a series of general transport equations:

@ne @t þ @

C

e @x ¼ Se; @ni @t þ @

C

i @x ¼ Si; @nn @t þ @

C

n @x ¼ Sn; @ðneh

e

iÞ @t þ @

C

E @x ¼ e

C

E E þ SE; ð28Þ

where subscripts e, i and n represent the electron, ion and neutral species number density, respectively. By applying drift-diffusion approximation[15], the ﬂux expressions (C) for each equation are:

C

e¼

l

eEne De @ne @x ;

C

i¼

l

iEni Di @ni @x;

C

e¼ 5 3

l

eEneh

e

i 5 3Di @ðneh

e

iÞ @x ; ð29Þ

where the coefﬁcients of mobility (

l

) and diffusion coefﬁcients (D) control the general motion of each species. The electron mobility and diffusion coefﬁcients are functions of the electron temperature, computed from the electron energy density h

e

i ¼3

2kBTe, with values taken from BOLSIG [22]. Charged species are inﬂuenced by the presence of an electric ﬁeld

E = rU, where the potential ﬁeld as predicted by the Poisson Equation:

@ @x

e

@

U

@x ¼ e Xni X ne : ð30Þ

The potential ﬁeld is solved implicitly using a conventional ﬁnite difference discretization and either a conjugate gradient or GMRES solver. A semi-implicit formulation[23]is required for computation of the source terms in Eq.(30)– otherwise the

Fig. 12. Shock-Bubble interaction problem with a shock wave located at x = 0 propagating into a low-density bubble of equal pressure. All simulations exploit symmetry around y = 0.

(14)

solution (regardless of approach) quickly becomes unstable. The diffusivity of ions is given by the Einstein relation with coef-ﬁcients of mobility computed from[24]:

l

i¼ 0:574T

Ppffiffiffiffiffiffiffiffim

a

; ð31Þ

where T and P are the background gas temperature and pressure, respectively, m is the reduced mass (given in AMU) and

a

is the polarizability of the background gas (more details can be found in[24]). Neutral species use diffusion coefﬁcients com-puted from the Chapman–Enskog theory of gases[25]. The source terms present in the governing ﬂuid modeling equations (Eq.(28)) (for example, Se) describe the production of species resulting from plasma and chemical reactions. These reactions

considered here are[26–29]:

e_{þ Ar ! e}_{þ Ar;} e_{þ Ar ! 2e}_{þ Ar}þ_; e_{þ Ar ! e}_{þ Ar} m; e_{þ Ar ! e}_{þ Ar} r; e_{þ Ar} m! 2eþ Arþ; e_{þ Ar} m! eþ Arr; 2Arm! Ar þ Arþþ e; Arr! Ar; Ar þ Arm! 2Ar: ð32Þ

Fig. 13. Numerical Schlierens taken from simulation results of HLL (left) and HLLG (right). Both simulations show results after t = 0.2 with CFL numbers of 0.5 (top) and 0.75 (bottom) using 200 200 cells. Gradients employed for reconstruction are calculated by ﬁnite differences and limited by the MINMOD limiter.

(15)

The investigated problem is described inFig. 14with a small region of space L = 2 cm separating two inﬁnitely large plates. One plate is grounded while another plate has an applied potential of 100 V oscillating according to a sinewave with frequency of 13.56 MHz. The gas separating the plates is of pressure P = 10 (Torr) with temperature T = 300 (K) with initial number densities (electrons and argon ions) of ne= np= 1E16 m3. The applied potential varies according to a sin wave with

frequency 13.56 MHz. Boundary conditions depend on the species: neutral species are treated using Neumann across each surface while charged species employ the following ﬂux:

C

i¼ a

l

iEniþ 1 4ni

v

i;th;

C

e¼ a

l

eEneþ 1 4ne

v

e;th;

v

i;e;th¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi 8kBTe

p

mi;e s ; ð33Þ

where a in this case is an integer dependent on the direction of the flux under the influence of the electric field (a = 1 if

l

E points toward the boundary, a = 0 otherwise). The initial condition is advanced in an unsteady time-stepping procedure until the pseudo-steady state is reached. This result is shown inFig. 15compared to a conventional second order TVD scheme using the same solution parameters. The electrons, being very light and easily inﬂuenced by the electric ﬁeld, tend to move back and forth across the region between the plates. This motion can be demonstrated by examining the distribution of elec-trons at various stages during the RF cycle.

Due to the very large electron convective velocities present (and resulting large Peclet numbers), the time step size must be restricted to quite a low number to ensure electrons do not transverse more than one cell in any given time step. In addi-tion, in the sheath regions where large electric fields are present, the upwind feature of HLLG is prevalent: usually only ever the left or right flux is used directly since the convective speeds are so large. These two features of this flow reveal that the pseudo-steady state for both the HLLG and the conventional second order TVD scheme are almost indistinguishable. How-ever, the unsteady progression to steady state is not: the HLLG scheme requires more RF cycles to reach the pseudo-steady state due to the reduced numerical dissipation in the scheme. In regions just outside the sheath region (where electric fields are very low), the solution is diffusion-bound. In this instance, the accuracy of HLLG provides an excellent solution; however, the increased number of RF cycles is required to reach the pseudo-steady state, despite the increase in accuracy. However, HLLG may become important in several transient simulations of discharge problems.

4. Conclusion

Presented is the derivation of the integral form of the HLL flux expressions which are then applied with spatially recon-structed values of conserved properties over cell widths. The obtained flux expressions, referred to as HLLG fluxes, proves that conventional extension to higher order through reconstruction at cell interfaces is mathematically inconsistent with the assumptions employed during the flux derivation and effectively results in an over-reconstruction of solution properties. The resulting flux expressions demonstrate that several reconstruction locations are required within the cell width as op-posed to the single reconstruction location at the cell interface. Through solutions of the Euler and shallow water equations, conventional extension of HLL to second order using various limiting functions is demonstrated to produce oscillations while the mathematically complete extension employed by HLLG is shown to be less prone to such oscillations. The method is fi-nally applied to the simulation of Argon RF plasma and compared to a conventional second order TVD method. In the course of the simulation, very large Peclet numbers are encountered in the sheath regions, causing the HLLG scheme effectively re-duces to a conventional second order upwind scheme. Thus, the comparison between the HLLG and conventional scheme shows virtually identical results, demonstrating that the additional accuracy present in HLLG flux expressions may not be Fig. 14. Schematic of simple 1D argon plasma layout. A sin wave potential is applied to the left plate (100 V, f = 13.56 MHz) while the right plate is grounded. A sheath region is expected to exist near the electrodes where the net charge is non-zero.

(16)

essential for the pseudo-steady plasma simulations. However, this additional desired accuracy may be required for plasma simulations of differing conﬁguration outside the scope of this investigation.

Acknowledgments

We would like to thank the ﬁnancial support of National Science Council and Ministry of Economic Affairs of Taiwan through Grant Nos. 96-2628-E-009-136-MY3 and 98-EC-17-A-07-S2-0043, respectively.

References

[1] S.K. Godunov, Finite difference method for the numerical computation of discontinuous solutions of the equations of ﬂuid dynamics, Math. Sbornik 47 (1959) 271–306.

[2] A. Harten, P.D. Lax, B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev. 25 (1) (1983) 35–61. [3] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer-Verlag, New York, 1997.

[4] D.M. Causon, D.M. Ingram, C.G. Mingham, G. Yang, R.V. Pearson, Calculation of shallow water ﬂows using a Cartesian cut cell approach, Adv. Water Res. 23 (5) (2000) 545–562.

[5] D.M. Ingram, D.M. Causon, C.G. Mingham, Developments in Cartesian cut cell methods, Math. Comput. Simul. 61 (3) (2003) 561–572. [6] E.F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock Waves 4 (1) (1994) 25–34. [7] P.L. Roe, Characteristic-based schemes for the Euler equations, Ann. Rev. Fluid Mech. 18 (1986) 337–365.

[8] B. van Leer, Towards the ultimate conservative difference scheme. III: Upstream-centered ﬁnite-difference schemes for ideal compressible ﬂow, J. Comput. Phys. 23 (1977) 263–275.

[9] M.R. Smith, H.M. Cave, J.-S. Wu, M.C. Jermy, Y.-S. Chen, An improved quiet direct simulation method for Eulerian ﬂuids using a second order scheme, J. Comput. Phys. 228 (6) (2009) 2213–2224.

[10] B. van Leer, Towards the ultimate conservative difference scheme, V. A second order sequel to Godunov’s method, J. Comput. Phys. 32 (1979) 101–136. [11] N.E. Kolgan, Application of the minimum-derivative principle in the construction of ﬁnite-difference schemes for numerical analysis of discontinuous

solutions in gas dynamics, Uchenye Zapiski TsaGI [Sci. Notes of Central Inst. of Aerodynamics] 3 (6) (1979) 68–77 (in Russian).

[12] N.E. Kolgan, Finite-difference schemes for computation of three dimensional solutions of gas dynamics and calculation of a ﬂow over a body under an angle of attack, Uchenye Zapiski TsaGI [Sci. Notes of Central Inst. of Aerodynamics] 6 (2) (1975) 1–6 (in Russian).

Fig. 15. Comparison of solutions from proposed HLLG scheme and a conventional second order TVD scheme applied to the solution of the ﬂuid modeling equations. Results are shown for the pseudo-steady state at various points during the RF cycle – 0 (top, left),p/2 (top, right),p(bottom, left) and 3p/2 (bottom, right). Results show number densities of Ar, Ar+

, Arm

and e_{as a function of location between the plates, electron temperature and potential field.} The gradient of the potential field (i.e. the electric field) is very low in regions outside the sheath regions.

(17)

[13] V.A. Titarev, E.F. Toro, ADER schemes for three-dimensional non-linear hyperbolic systems, J. Comput. Phys. 204 (2) (2005) 715–736. [14] M. Ben-Artzi, J.-Q. Li, G. Warnecke, A direct Eulerian GRP scheme for compressible fluid flows, J. Comput. Phys. 218 (1) (2006) 19–43. [15] E. Han, J.-Q. Li, H.-Z. Tang, An adaptive GRP scheme for compressible fluid flows, J. Comput. Phys. 229 (5) (2010) 1448–1466.

[16] J.-Q. Li, G.-Z. Chen, The generalized Riemann problem method for the shallow water equations with bottom topography, Int. J. Numer. Methods Eng. 65 (6) (2005) 834–862.

[17] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.

[18] M.R. Smith, M.N. Macrossan, M.M. Abdel-Jawad, Effects of direction decoupling in ﬂux calculation in ﬁnite volume solvers, J. Comput. Phys. 227 (8) (2008) 4142–4161.

[19] M. Cada, M. Torrilhon, Compact third-order limiter functions for finite volume methods, J. Comput. Phys. 228 (11) (2009) 4118–4145. [20] C.-W. Shu, S. Osher, Efficient implementation and essentially non-oscillatory shock capturing schemes II, J. Comput. Phys. 126 (1989) 32–87. [21] P. Woodward, P. Collela, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1984) 115–174. [22] SIGLO, BOLSIG: user-friendly Boltzmann solver from the SIGLO series, Webpage last visited on 5th Feb, 2010. Available from:

<http://www.siglo-kinema.com/bolsig.htm>.

[23] H.C. Kim, F. Iza, S.S. Yang, M. Radmilovic´-Radjenovic´, J.K. Lee, Particle and ﬂuid simulations of low-temperature plasma discharge: benchmarks and kinetic effects, J. Phys. D: Appl. Phys. 38 (2005) 283–301.

[24] D. Herrebout, A. Bogaerts, M. Yan, R. Gijbels, W. Goedheer, E. Dekempeneer, One-dimensional ﬂuid model for an rf methane plasma of interest in deposition of diamond-like carbon layers, J. Appl. Phys. 90 (2) (2001) 570–579.

[25] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, second ed., John Wiley & Sons, 2007.

[26] M. Hayashi, Bibliography of Electron and Photon Cross Sections with Atoms and Molecules Published in the 20th Century: Argon, Technical Report NIFS-DATA-72, National Institute for Fusion Science, Japan, 2003.

[27] E.A. Bogdanov, A.A. Kudryavtsev, L.D. Tsendin, R.R. Arslanbekov, V.I. Kolobov, V.V. Kudryavtsev, The inﬂuence of metastable atoms and the effect of the nonlocal character of the electron distribution on the characteristics of the positive column in an argon discharge, Tech. Phys. 49 (6) (2004) 698–706. [28] T.V. Rakhimova, O.V. Braginsky, V.V. Ivanov, T.K. Kim, J.T. Kong, A.S. Kovalev, D.V. Lopaev, Y.A. Mankelevich, O.V. Proshina, A.N. Vasilieva, Experimental

and theoretical study of RF plasma at low and high frequency, IEEE Trans. Plasma Sci. 34 (3) (2006) 867–877.

[29] M.W. Kiehlbauch, D.B. Graves, Modeling argon inductively coupled plasmas: the electron energy distribution function and metastable kinetics, J. Appl. Phys. 91 (6) (2002) 3539.