c
TU Delft, The Netherlands, 2006
DG AND MUSCL APPROACHES : COMPARISON OF
SOME RECENT VARIANTS FOR AEROACOUSTICS AND
VORTEX COMPUTATIONS
Michel Borrel∗
O.N.E.R.A.
BP 72, 92322 Chˆatillon cedex, France
∗e-mail: michel.borrel@onera.fr
Key words: MUSCL, Discontinuous Galerkin Method, Planar Wave, Vortex Flow Abstract. Several recent MUSCL and DG developments are presented and tested in this study on very simple test cases concerning plane wave propagation , a vortex and the interaction of a sinusoidal wave with a shock. These two approaches are currently in constant progress and now can be applied to new fields such as non linear aeroacoustics or vortex aerodynamics, fields formerly reserved to centered approaches.
1 INTRODUCTION
In this preliminary study, we compare the behavior of two classical numerical ap-proaches in compressible fluid dynamics , the MUSCL and the Discontinuous Galerkin (DG) method, applied to aeroacoustics and vortex computation, these two fields being related as noise is often produced from vortex fluctuations.
The calculation of the unsteady effects with a minimum of dissipation and numerical dispersion is then of primary importance and it is necessary to re-examine, accordingly, these methods and their multiple variants which often were established and validated for steady flows with strong discontinuities (cf. [1] or [2] for example); in addition, these methods were recently the object of important improvements (cf. for example [3], [4], [9], [10]) which must be integrated and tested. As we are interested here only in propagation problems, we will neglect the viscous effects.
2 MUSCL and DG approaches
Following the works of many authors, MUSCL and DG are now standard for solv-ing conservation law equations. In our code, both MUSCL and DG-P 1 approaches are implemented.
∂tW +
− →
∇ · F(W) = 0 (1)
where W is the conservation variable vector :
W = (ρ, ρ−→U , ρE)T (2)
and F the flux vector : f = ( ρ u , ρ u2+ p , ρ uv , ρ uw , u (ρE + p ) )T g = ( ρ v , ρ uv , ρ v2+ p , ρ vw , v (ρE + p ) )T h = ( ρ w , ρ uw , ρ vw , ρ w2+ p , w (ρE + p ) )T (3)
ρ is the density,−→U the velocity , E the total energy and p the pressure given by the state law for perfect gas (the specific heat ratio γ is supposed to be constant).
2.2 Weak formulation
At each time t ∈ [0, T ], an approximate solution Wh(x, y, z, t) is defined from an initial
solution W0(x, y, z) and boundary conditions. The computational 3D domain Ω is parti-tioned using an i,j,k structured mesh. The centre of each cell Ωijkis noted xOijk, yOijk, zOijk.
In each cell, the solution is expanded on a local basis of degree one polynomials: ∀(x, y, z) ∈ Ω, ∀t ∈ [0, T ], Wh(x,y,z,t) = X i,j,k X `=0,3 p`ijk(x, y, z) ¯W`ijk(t) (4) where {p`ijk, ` = 0, 3} is the local basis of Legendre polynomials :
p0
ijk(x, y, z) = 1 if (x, y, z) ∈ Ωijk, 0 else
p1
ijk(x, y, z) = x − xOijk if (x, y, z) ∈ Ωijk, 0 else
p2ijk(x, y, z) = y − yOijk if (x, y, z) ∈ Ωijk, 0 else
p3ijk(x, y, z) = z − zOijk if (x, y, z) ∈ Ωijk, 0 else
(5)
¯
W`ijk are the degrees of freedom on Ωijk which approximate the mean and gradient cell
values. In both approaches, the discretisation is built from the weak formulation of (1) obtained by multiplying by a test function ϕh in a space Vh and by integrating by parts
on Ω : Z Ωijk ϕijk ∂ ∂tWhdΩ+ I ∂Ωijk ϕijkF(W˜ h) · −→n dσ − Z Ωijk − → ∇ϕijk· F(Wh)dΩ = 0 (6) with :
Vhn = {ϕh ∈ L1(Ω) / ∀(i, j, k), ϕh|Ωijk = ϕijk ∈ Pn(Ωijk)} (7)
where Pn(Ω
ij) represents the polynomials of degree at least n on Ωijk. For DG, we
polynomials as a basis leads to the ODE system : d dt ¯ W0ijk Z Ωijk dΩ = − I ∂Ωijk ˜ F · −→n dσ d dt X l=1,3 ¯ Wlijk Z Ωijk plijk pmijkdΩ = − I ∂Ωijk pmijkF · −˜ →n dσ + Z Ωijk F ·−→∇pm ijkdΩ (m = 1, 2, 3) (8)
For MUSCL, only the first equation is relevant and for DG, if a Cartesian grid is used, the mass matrix of the second equation is diagonal with :
M = ∆x ∆y ∆z ∆x2 12 , ∆y2 12 , ∆z2 12 I (9)
2.3 Reconstruction techniques and slope limiters
Introducing nonlinear slope limiters that ensure stability but doesn’t degrade the ac-curacy of the method is a problem that was largely studied in the last thirty years either for structured or unstructured grids. With structured grids, the problem is reduced easily to the one dimensional case considering only the values interpolated on both sides of each interface. The triad introduced by G. Billet in [5] and [6] was used successfully here:
WRijk−1/2 = Wijk − 14[(1 − κ`)∆+` + (1 + κ`)∆−` ] WLijk+1/2 = Wijk + 14[(1 − κ`)∆−` + (1 + κ`)∆+` ]
∆+` = M ax(0, M in( A , ω`B) sgn(Wijk+1− Wijk)
∆−` = M ax(0, M in( B , ω`A) sgn(Wijk− Wijk−1)
A = (Wijk+1− Wijk) sgn(Wijk− Wijk−1)
B = (Wijk− Wijk−1) sgn(Wijk+1− Wijk)
(10)
In this formula, the indices and their increments are taken relatively to each grid direction and we refer to the papers quoted for the precise definition of parameters κ` and ω`
which uses the five values Wijk−2, ..., Wijk+2. Let us notice that this technique does not
reconstruct a solution in V1
h since to each interface is attached a different polynomial in
each cell.
According to our experience with DG, it is necessary to resort to slope limiters only for particular problems, (cf [7] [8]), for which the use of a cut-off technique, which consists, near discontinuities for example, in locally going down to a P0approximation, has revealed
insufficient.
Slope limiters for DG-P1 with the triad consist in taking the Minmod value between the corresponding degrees of freedom and gradients calculated with the limiter in each grid direction :
¯
2.4 P2 MUSCL variant for DG
A similar reconstruction technique can be applied to DG thus making it possible to obtain an evaluation of the second derivatives. It is important to use an orthogonal or quasi-orthogonal polynomial basis such as the Legendre polynomials so that the weak formulation (6) does not add additional terms in order to remain consistent. With the polynomial basis (5), it is enough to add in Cartesian grid :
pxx ijk(x, y, z) = 1 2(x − xOijk) 2 − ∆x2 24 if (x, y, z) ∈ Ωijk, 0 else pyyijk(x, y, z) = 12(y − yOijk) 2 −∆y2 24 if (x, y, z) ∈ Ωijk, 0 else pzzijk(x, y, z) = 12(z − zOijk) 2 − ∆z2 24 if (x, y, z) ∈ Ωijk, 0 else
pxyijk(x, y, z) = (x − xOijk)(y − yOijk) if (x, y, z) ∈ Ωijk, 0 else
pyzijk(x, y, z) = (y − yOijk)(z − zOijk) if (x, y, z) ∈ Ωijk, 0 else
pxzijk(x, y, z) = (x − xOijk)(z − zOijk) if (x, y, z) ∈ Ωijk, 0 else
(12)
The polynomial development (4) becomes : ∀(x, y, z) ∈ Ωijk, Wh(x, y, z) =
X
`=0,3
p`ijk(x, y, z) ¯W`ijk + pxxijk(x, y, z) ¯Wxxijk+ pxyijk(x, y, z) ¯Wxyijk + ...
(13) As for the classical MUSCL approach, to spare memory, the reconstructed quantities ¯
Wxxijk, can be calculated direction by direction. Moreover, by using the Gauss-Taylor or Taylor variants presented below, the cross terms ¯Wxyijk, do not occur. Several formulas were tested in this study, among which:
¯
Wxxijk = Wijk+1 − 2 Wijk + Wijk−1
!
/ ∆ξ2 or ¯
Wxxijk = ave( ¯W`ijk+1 − ¯Wijk` , ¯W`ijk − ¯W`ijk−1) / ∆ξ , ` = 1, 2, 3
(14)
where the ave function is equal to either the arithmetic mean or to a classical limiter Minmod or Superbee. First results in linear propagation of an acoustic wave on the Shu-Osher test case show that one cannot use the same formula for all the problems if one wants to optimize results .
2.5 Flux integrals Computation
Three types of flux integrals are to be computed with DG and only one with MUSCL : I ∂Ωijk ˜ F · −→n dσ , I ∂Ωijk pmijkF · −˜ →n dσ and Z Ωijk F ·−→∇pm ijkdΩ (15)
Roe or Lax-Friedrich with DG. Flux integrals are evaluated by numerical quadrature formulas. With MUSCL, a one point quadrature formula is used :
I ∂Ωijk F (x, y, z) dσ M U SCL≈ X k=1,6 |∂Ωk ijk|F (xck, yck, zck) (16) where {∂Ωk
ijk, k = 1, 6} represents the control volume faces; (xck, yck, zck) the
coor-dinates of their center. With DG, three variants were tested for the evaluation of the surface integrals. The first one consists classically in using a four point Gauss quadrature formula : I ∂Ωijk F (x, y, z) dσ DG−Gauss≈ 1 4 X k=1,6 |∂Ωk ijk| X i=1,4 F (xgi, ygi, zgi) (17)
where {(xgi, ygi, zgi) , i = 1, 4} represents the coordinates of the Gauss points. This
vari-ant, called DG-Gauss, implies carrying out four evaluations of numerical fluxes per face; it is more expensive in time computing but it can be regarded as the reference version since it can be extended to a polynomial basis of unspecified degree if one wishes to go up in order. A second variant, called DG-Gauss-Taylor, exploits the idea of Van de Vegt of a Taylor development to the first order of the numerical flux to avoid the evaluation at each inte-gration point : starting from the same formula of inteinte-gration as (17), the values at the integration points are obtained by means of the Taylor development at the center value :
F (xgi, ygi, zgi) ∂Ωk ijk DG−Gauss−T aylor ≈ F (xck, yck, zck) + ξ1 ∂F ∂ξ1 + ξ2 ∂F ∂ξ2 (18) where (ξ1, ξ2) represents the surface local coordinates. This method requires the
calcu-lation of the jacobian matrices ∂F ∂ξi = ∂F ∂WR ∂WR ∂ξi + ∂F ∂WL ∂WL ∂ξi ( i = 1, 2) (19)
which can be evaluated in an approximated way, for example with Roe fluxes : ∂ ˜FRoe(WR, WL)
∂WR/L · −
→n ≈ 1
2R (Λ ± |Λ|) R
−1 (20)
where R, Λ and |Λ| respectively represent the right eigenvectors matrix , the diagonal eigenvalues matrix and the eigenvalue modules associated to the Jacobian. Use of formula (4) allows then an explicit evaluation.
The third variant, called DG-Taylor, introduced by van der Vegt and van der Ven [9], consists in integrating the Taylor development (18) in the form of a polynomial in ξ1 and
I ∂Ωijk F (x, y, z) dσ DG−T aylor≈ X k=1,6 I ∂Ωk ijk F (xck, yck, zck) + ξ1 ∂F ∂ξ1 + ξ2 ∂F ∂ξ2 ∂ξ1∂ξ2 (21)
In the three DG variants, the volume integrals attached to fluxes are evaluated by means of an eight point quadrature formula :
I Ωijk F (x, y, z) dΩ DG≈ 1 8|Ωijk| X m=1,8 F (xgm, ygm, zgm) (22) 2.6 Time discretization
Time integration of system (8) written in a compact form dW
dt = L(W) (23)
is realized by means of an explicit p step Runge-Kutta method (RKp):
W0 = Wn W1 = W0 + C1∆t L(W0) W2 = W0 + C2∆t L(W1) · · · · · Wp = W0 + Cp∆t L(Wp−1) Wn+1 = W p (24)
where ∆t is the time step calculated by means of a CFL type criterion and {C1, C2, ...Cp}
are constants ranging between 0 and 1 classified by ascending order (Cp = 1 for
con-sistency). The advantage of this method of integration, in addition to being easily im-plementable, is that it requires identical memory storage for any value of p. Moreover, RK3 with the optimized parameters {.17, .5, 1} is of order three with good properties of dissipation (TVD, cf Van Leer). For CPU reasons attached to large 3D calculations, we restricted ourselve to RK1 (Euler scheme), RK2 (C1 = .5, Heunn type scheme) and RK3,
schemes which appeared sufficient for the application presented.
2.7 RK1(L-W) variant for DG
development in time (of the order of the selected polynomial i e. first order in our case), then to replace the time derivative by a space derivative by means of the quasilinear sys-tem (1) with the flux Jacobian matrix (and its powers for the higher order developments) : Wijk(x, y, z, t) = Wijk(x, y, z, 0) + t ∂W ∂t = Wijk(x, y, z, 0) − t dF dW · − → ∇W (25)
After time integration :
1 ∆t Rtn+1 tn W¯ 0 ijk(t) dt = ¯W 0 ijk(0) − ∆t2 ( dfdWW¯ 1 ijk(0) + dg dWW¯ 2 ijk(0) + dhdWW¯ 3 ijk(0)) 1 ∆t Rtn+1 tn W¯ ` ijk(t) dt = ¯W ` ijk(0) , ` = 1, 2, 3 (26)
the Jacobian matrices are :
df dW = 0 1 0 0 0 bQ22 − u2 (2 − b)u −bv −bw b −uv v u 0 0 −uw w 0 u 0 u[(b − 1)Q22 − a2 b ] − a2 b + Q2 2 − bu
2 −buv −buw γu
, dg dW = ... with b = γ − 1, a2 = γp/ρ and Q2 =−→U ·−→U .
As many terms are identical or null in these Jacobian matrices , this correction costs only little CPU time compared to a complete Runge-Kutta time step. The weak formu-lation then naturally provides a RK1 scheme in time, flux integrals being evaluated as previously with the new degrees of freedom (26) and the volume integrals evaluated at Gauss space/time points :
Z Ωijk ϕijk(Wn+1h − W n h) dΩ + ∆t I ∂Ωijk ϕijkF(W˜ n+1/2h ) · −→n dσ − Rtn+1 tn dt Z Ωijk − → ∇ϕijk· F(Wh(t))dΩ = 0 (27)
Finally, let us remark that a similar variant can be devised for MUSCL.
3 COMPUTATIONAL RESULTS
3.1 Plane wave, 1D case
One considers the problem of an unsteady plane wave aligned with a Cartesian grid (cf [11]). It is a 1D problem but with a 3D treatment: the grid used is uniform in the three directions with 80 cells along the x axis, 40 along the y axis the and 2 along the z axis which is degenerated ). The calculation field is 200mm long and the size of the cells are 2, 5mm. The simulation consists, starting from an initial state defined by the uniform flow at Mach number 0, 4 along the x axis in imposing a (static) pressure fluctuation downstream
p(t) = ¯p + ∆p sin(2π t/T ) (28)
where the amplitude ∆p is 89 P a for a total pressure of 101 kP a and the period defined by T = N ∆x / c with c the speed of sound and N the wave number which was set equal to 20, which corresponds to a signal of approximately 6700Hz, characteristic of turboshaft engines problems. The wave goes upstream, where a non reflection condition based on characteristics is imposed, while at top and bottom of the field of calculation a condition of periodicity is applied and on the sides a condition of symmetry (see figure 1a). With this data file, the length of the field of calculation corresponds to six and a half wave length , each wavelength being discretized on approximately twelve cells , which can seem relatively low but which answers a reasonable criterion of memory storage for large 3D calculations 3D.
Figure 1: Plane wave propagation : a)2D mesh view b)pressure field at steady state
To refine our comparison and to remove the effects of the upstream non reflection condition , all following calculations were stopped at time t = 200. Figure 3 shows the effects of limiters and time discretization for MUSCL : the dependence in ∆t with RK1 is obvious as well as the dissipation of the M inmod limiters.
x P/ Pre f 0 50 100 150 200 1.21 1.211 1.212
MUSCL (Roe) - triad - RK1 (CFL .1) MUSCL (ausm+) - triad - RK1 (CFL .1) MUSCL (Roe) - triad - RK2 DG (Roe) - RK2
Figure 2: Plane wave propagation: DG and MUSCL comparison at steady state
x P/ Pre f 0 50 100 150 200 1.21 1.211 1.212 DG-RK2 triad-RK2 triad-RK1 (CLF .1) triad-RK1 (CFL .4) Minmod-RK1
DG / MUSCL , flux de Roe , t=200
Figure 3: Plane wave propagation: limiter influ-ence on MUSCL
x P/ Pre f 50 100 150 1.21 1.211 1.212 DG-Roe MUSCL-Roe MUSCL-van Leer MUSCL-ausm+ DG / MUSCL (triad) , RK2 , t=200
Figure 4: Plane wave propagation: numerical flux formulation effect on MUSCL
x P/ Pre f 0 50 100 150 200 1.21 1.211 1.212 RK2 RK1 (CFL .4) RK1 (CFL .2) RK1 (CFL .1) RK1 (CFL .05)
MUSCL - triad , flux de Roe , t=200
Figure 5: Plane wave propagation:time discreti-sation effect on MUSCL
Figure 6 compares various DG variants: identical results were obtained with DG-Gauss (reference variant) and DG-DG-Gauss-Taylor with a gain of 15% in CPU time in favor of the second variant; very close results are obtained for DG-Gauss with Roe fluxes and Lax-Friedrich fluxes , with this time a gain of 10 % in CPU time for the latter; if DG-Gauss-RK1 is unstable (what is a known result and already analyzed), it appears that DG-Taylor behaves in a similar way to the MUSCL-triad: with RK2 one obtains rather diffusive results and with RK1, for an optimal CFL of 0, 1 one obtains a minimum dissi-pation with now a CPU time saving of 22% compared to a DG-Gauss calculation.
Figure 7 shows that with an identical number of degrees of freedom (in 1D), therefore on a grid twice as coarse, DG-RK2 behaves like MUSCL-RK2, but with a gain of 20% in CPU time.
The CPU times obtained on a scalar computer SGI-Origin-2000 for these calculations is presented in the following table:
MUSCL MUSCL DG-Gauss-T DG-Gauss DG-Gauss DG-Gauss-T DG-Taylor
RK2 RK1 RK2 RK2 RK2 RK2 RK2
CFL=.4 CFL=.1 CFL=.2 CFL=.2 CFL=.2 CFL=.2 CFL=.2
Roe Roe 2dX Roe Roe Lax-Friedrich Roe Roe
298 594 242 1700 1546 1450 1330
x P/ Pre f 0 50 100 150 200 1.21 1.211 1.212 Gauss-Roe-RK2 Gauss-LaxFr-RK2 Gauss-RK1 (CFL .1) Taylor-RK2 Taylor-RK1 (CFL .1) DG , t=200
Figure 6: Plane wave propagation: DG variants comparison x P/ Pre f 0 50 100 150 200 1.21 1.211 1.212 DG, dX (CFL .2) DG, 2dX (CFL .2) MUSCL-triad, dX (CFL .4)
DG - MUSCL , RK2 , flux de Roe , t=200
Figure 7: Plane wave propagation: DG/MUSCL comparison for 2 meshes
On figures 8, 9, 10, 11, we increased the size of the calculation field following the x axis in order to carry out simulations over longer times. Figure 8 shows a comparison DG/MUSCL with RK2: if the quality of DG results, it should be noted however that the DG overcost is of a factor 5 in CPU time.
x P/ Pre f 1000 1500 2000 2500 1.2095 1.21 1.2105 1.211 1.2115 1.212 DG-GT-RK2 MUSCL-RK2 DG / MUSCL , Nc=12 , t=2500
Figure 8: Plane wave propagation: DG/MUSCL with RK2 comparison x P/ Pre f 0 50 100 150 200 1.21 1.211 1.212 DG, dX (CFL .2) DG, 2dX (CFL .2) MUSCL-triad, dX (CFL .4)
DG - MUSCL , RK2 , flux de Roe , t=200
Figure 9: One meter length plane wave propaga-tion: MUSCL results with RK1 and 2 CFL values
RK1 (LW) for DG and MUSCL. The excellent results obtained with DG-MUSCL (P2) were carried out without limiters on second derivatives, as introducing limiters would be going back to the basic solution. The results with RK1 (LW) are better than the basic solutions for each approach, the improvement for MUSCL being spectacular while for DG, it is more a saving of CPU time of 25 % in spite of numerical fluxes at the interfaces which were evaluated at the Gauss space/time points, i.e with twice more evaluations. The CPU time ratio between the DG calculation and the MUSCL calculation, carried out with the same CFL, is of four.
x P/ Pre f 1000 1500 2000 2500 1.2095 1.21 1.2105 1.211 1.2115 1.212 DG-P1/MUSCL2 DG-P1 DG , Nc=12 , t=2500
Figure 10: One meter length plane wave prop-agation:DG(P1) comparison with and without MUSCL(P2) correction x P/ Pre f 1000 1500 2000 2500 1.2095 1.21 1.2105 1.211 1.2115 1.212 DG MUSCL DG / MUSCL , RK1(L-W) , Nc=12 , t=2500 CFL = .2
Figure 11: One meter length Plane wave propa-gation: DG-MUSCL(P2) and RK1(L-W) variants evaluation
3.2 Oblique plane wave, 2D case
The problem corresponds to a test case suggested in the European project TurboNoise CFD to evaluate computer codes on noise propagation: the flow is identical to the pre-ceding case but the acoustic wave propagation is carried out with a 45 degrees angle with the grid. The downstream pressure fluctuation is related now to t and y:
p(t) = ¯p + ∆p sin(2π (t/T + y/L)) (29)
where L = N ∆y , with N the wave number and ∆y the space step . The boundary conditions are identical to the preceding case: it is a question of testing the no-reflection upstream boundary conditions. The technique uses characteristic relations which consists in defining these as:
where −→n is the boundary normal, either the value in the border cell or the value of a state of reference, according to the sign of the associated eigenvalue {−→U · −→n ± c , −→U · −→n } and to then deduce the value of the flux.
We tested the idea suggested in [11] which consists in taking for −→n the unit vector aligned with−→∇ p whose value is evaluated starting from the computed field at the border. We present on figure 12 the improvement made by this treatment which makes it possible to reduce the stray reflections, without however removing them.
Figure 12: Oblique wave propagation: no reflection boundary condition influence
3.3 Advection of an isolated vortex
In order to test MUSCL and DG approaches for turbulent flows, we chose the problem of the advection of a two-dimensional isolated and isentropic vortex proposed per H. Yee [12], an exact solution the Euler equations . This vortex is initialized by adding to the principal flow a disturbance on the speed and temperature fields :
(δu, δv) = 2πβ e1−r22 (−y + y0, x − x0) δT = −(γ−1)β8γπ22e1−r 2 (31)
where r = q(x − x0)2+ (y − y0)2 is the distance to the center of the vortex which moves
with a velocity−→U∞ and β = 5 the intensity of the vortex . The initial solution is written
In our calculations, we took ρ∞ = 1, T∞ = 1 and
−→
U∞ = (1, 0, 0). At every moment tn,
the solution is a simple translation of the initial field: Wn = W0((x, y, z) −−→U∞tn).
A strategy of grid refinement of the AMR type was used with a Cartesian and uniform coarse grid ( DeltaX = 0, 5) for which the diameter of the vortex is represented on approximately eight coarse cells. Calculations were carried out until dimensionless time t = 50 which corresponds to an advection on 20 diameters of the vortex whereas its core has carried out approximately six rotations. 2800 iterations in time are necessary with a time step corresponding to a CFL of 0,2. Calculations with a refinement by two and four were carried out (figure 13). With the same level of refinement, DG calculation is four times more expensive than a MUSCL calculation , on the other hand the diffusion and especially the dispersion starts to be very important, even with the best variant ( Roe flux and triad reconstruction). To obtain results comparable with those obtained by DG, a refinement by 4 is necessary, the overcost in computing times of DG on MUSCL being now only of 20 %.
3.4 Shu-Osher test case
It is a 1D problem of a shock at Mach 3 meeting a density fluctuating flow . This interaction leads to high frequencies harmonics behind the shock, difficult to simulate numerically and only on very refined grids. The initial conditions are:
ρ = 3.857143 u = 2.629369 p = 10.33333 si −5 ≤ x < −4
ρ = 1 + 0.2 sin5x u = 0. p = 1. si 5 ≥ x ≥ −4
-5 0 5 0 1 2 3 4 5 Sol. Ref. DG ρ a) ρ
DG, RK3, Roe triad, 300 mailles ρ
Figure 14: Shu and Osher test case: DG compu-tation on a 300 cell mesh
-1 0 1 2 3 3 4 Sol. Ref. MUSCL DG ρ
MUSCL/DG, RK3, Roe triad, 300 mailles ρ
b)
Figure 15: Shu and Osher test case: DG/MUSCL comparison
The results are not however entirely satisfactory since, because of the introduction of limiters, the results with DG on 200 cells are definitely worse than results obtained with MUSCL on 400 cells (i.e. with the same number of degrees of freedom) (figure 15). On the other hand, one can see on figure 16 that variant DG-RK1 (LW) - MUSCL (P2) provides as good results on a 300 cell mesh as MUSCL-(RK2) on a 400 cell mesh with an identical CPU time cost .
-1 0 1 2 3 3 4 Sol. Ref. DG-RK3 DG-RK1(LW)-MUSCL(P2) ρ
DG, flux de Roe + triad, 300 mailles ρ a) -1 0 1 2 3 3 4 Sol. Ref. MUSCL ρ
MUSCL, RK2, flux de Roe + triad, 400 mailles ρ
b)
Figure 16: Shu and Osher test case: a)DG variants comparison (gauss and RK1(LW)- MUSCL(P2)); b)MUSCL computation on a 400 cell mesh
4 CONCLUSION
of a sinusoidal wave with a shock. It appears that these two approaches are currently in constant progress and become now applicable to new fields such as non linear aeroacoustics or vortex aerodynamics, fields formerly reserved to centered approaches. In aerodynamics, traditional approaches (cf [2] - [8]), obtained , on the same grid, a factor of about 12 in CPU time between a MUSCL calculation and a DG calculation , resulting on the one hand from a factor three due to computation cost per point and iteration and on the other hand from a more severe criterion of stability for DG (factor two for the CFL plus two for the time scheme). Considering in addition the problems of memory storage , in particular for the implicit scheme, all this explains why the DG approach has only few followers and especially in structured grids. For aeroacoustics, the traditional criterion of stability founded on the Fourier analysis does not seem sufficient any more and a more restrictive criterion founded on a maximum displacement of the wave during a time step seems more suitable, thus a limitation of the time step of the same order of magnitude for MUSCL and DG.
For calculations presented, with equal precision, CPU time cost of a DG calculation (initial version) compared to a MUSCL calculation is only now a factor of 2 to 3 according to the selected variant.
For an equal precision and similar memory storage , if one eliminates RK1 calculations which present too strong a dependence in time, MUSCL is even 20% times more costly than DG. Because of larger potentialities for the DG approach such as , for example, the rise in order while keeping the compactness or the use of non structured grids while keeping the precision, the results presented confirm the interest of DG, for aeroacoustics as well as unsteady aerodynamics.
In addition, first results with new developments MUSCL (P2) and RK1 (LW) seem very promising concerning precision and the speed of calculations; their evaluation must be continued. Lastly, let us finish by noting that effective absorbing boundary conditions must be established for the two MUSCL and DG approaches .
REFERENCES
[1] B. Berde and M. Borrel Comparison off High-order Godunov-type design for the Euler equations one irregular meshes, ECCOMAS, Paris, (1994).
[2] B. Berde. Study of a discontinuous finit cells method of Gakerkine type: resolution of the equations of Euler and Navier-Stokes on irregular grids, Thesis of the University Pierre and Marie Curie, (1995).
Comput. Sc. and Eng., vol. 9, Springer Ed. High-Order Methods for Computational Physics, (1999).
[4] B. Cockburn. Discontinuous Galerkin Method for Problems Convection-Dominated, Reading Notes in Comput. Sc. and Eng., vol. 9, Springer Ed., High-Order Methods for Computational Physics, (1999).
[5] G. Billet. A simple algorithm to improve the accuracy off TVD MUSCL approach, Proceedings of the 7th International Conference one Hyperbolic Problems, Theory,
Numerics and Applications, ETH Zurich, (1998).
[6] G. Billet and O. Louedin. Adaptive limiters for improving the accuracy of the MUSCL approach for unsteady flow, J. Comput. Phy., 170, 161–183, (2001).
[7] C. Drozo, M. Borrel and A. Lerat Discontinuous Galerkin schemes for the com-pressible Navier-Stokes equations, Proceedings of the 16th ICNMFD, Springer Ed.,
266–271, (1998).
[8] C. Drozo. Method of Discontinuous the Galerkine type for the resolution of the Navier-Stokes equations into compressible, Thesis of ENSAM, (1998).
[9] J.J.W. van der Vegt and H. van der Ven. Space time discontinuous Galerkin finite compressible element method with dynamic grid motion for inviscid flows. II. Efficient flow squaring, Comput. Meth. Appl. Mech. Eng., 191, 4747–4780, (2002).
[10] M. Dumbser and C.-D. Munz. ADER Discontinuous Galerkin Designs for Aeroa-coustics, Preprint submitted to Elsevier Science.
[11] V. Couaillier. Effective multidimensional non reflective boundary condition for CFD calculations applied to turboengine aeroacoustics prediction, ISABE Paper 2005-1185, 17th Int. Symp. on Airbreathing Eng., (2005).
