Interface flux reconstruction method based on optimized weight essentially non-oscillatory scheme

Peixun YU,Junqing BAI,Hi YANG,Song CHEN,Ki PAN

aNorthwestern Polytechnical University,Xi’an 710072,China

bAircraft Strength Research Institute of China,Laboratory of Aeronautical Acoustics and Dynamics,Xi’an 710065,China

cHubei Space Aircraft Research Institute,Wuhan 430040,China

1.Introduction

In the past few decades,high order finite difference methods have been developed in ComputationalAeroAcoustics(CAA).Most of these methods evolve from the Dispersion-Relation-Preserving(DRP)scheme proposed by Tam and Webb,1or the compact finite difference scheme by Lele.2Though having gained great success in CAA,the high order finite difference methods cannot be applied to practical CAA multi-scale problems with complex geometry difficultly,since it requires that the flux across the grid interface has to be continuous.When facing aeroacoustics problems with complex geometries,multi-size structured grids are usually employed,which greatly improves the capability of application in some practical engineering problems such as aircraft high lift systems and landing gears.However,how to construct the interface flux efficiently and stably is still a hot issue.

Aimed at multi-scale CAA problems,Tam and Kurbatskii3proposed low-dissipation&low-dispersion DRP scheme by optimizing the interpolation parameters with comparing multiple grid scales ratio.Besides,Tam improved the stencil DRP scheme and artificial selective damping terms which are used in the interface flux construction,in order to compute the flux accurately.

WENO scheme was initially proposed by Liu and Osher4based on Essentially Non-Oscillatory(ENO)scheme.5The weights depend on the local smoothness of the data.Smoothness measurements cause stencils that span large flow field gradients to have relative small weights; any candidate stencil containing a shock receives a nearly zero weight.In completely smooth regions,weights revert to optimal values,where optimal value is defined by maximum order of accuracy or maximum bandwidth.In the following years,Jiang and Shu6,7cast WENO scheme in to finite difference form.This scheme,which is referred to as WENO-JS hereinafter,can be capable of resolving shocks with high resolution.However,WENO-JS is too dissipative for the detailed simulation of turbulent flow.Martin et al.8,9developed two new formulations of a symmetric WENO method for the direct numerical simulation of compressible turbulence.The schemes are designed to maximize order of accuracy and bandwidth,while minimizing dissipation.The formulations and the corresponding coefficients are introduced.Numerical solutions to canonical flow problems are used to determine the dissipation and bandwidth properties of the numerical schemes.In addition,the suitability and accuracy of the bandwidth-optimized schemes for direct numerical simulations of turbulent flows are assessed in decaying isotropic turbulence and supersonic turbulent boundary layers. Wu et al.10presenteda maximum order preserving optimized WENO scheme which is a weighted average of the maximum order scheme and an optimized scheme.Hou et al.11modified the weights of WENO scheme to be smoother,which can eliminate the fluctuation within the grid variation.Lin and Hu12presented two groups of WENO schemes based on the dissipation and dissipation,which are investigated for computational aeroacoustics.

In this paper,the modified WENO scheme is deduced firstly.Then,we apply modified WENO scheme to establish the model of interface flux.On this basis,the accuracy of the model in multi-scale grid problems has been verified by several standard verification cases.The numerical simulation results with interface flux method are presented in Section 5,and the brief conclusions are given in Section 6.

2.Optimized WENO scheme

In the design of a traditional WENO scheme,the practice is to maximize the order of accuracy of the WENO scheme given the size of the difference stencil.However,high order schemes may not be the best for CAA problems.Aimed at short waves,the finite difference scheme needs to maintain its lowdissipation and low-dispersion property,as shown by Tam and Webb.To fix the idea,the initial value problem associated with the scalar wave equation is considered as follows:

where x is coordinate,t is time,u is function of x and t,a is constant number.

Given a uniform grid xi=iΔx with the same grid spacing Δx,the semi-discretized form of Eq.(1)is

where i is integer number,aui+1/2and aui-1/2are numerical fluxes which depend on k(r+s+1=k,r≥0,s≥0;r and s are different integer numbers respectively)grid points including xiitself,i.e.,

Here crjcould be obtained by achieving k-th order accuracy in Taylor series truncation expansion.Considering five-point stencil WENO-JS scheme,according to Ref.6,it would be divided to 3 candidate stencils{S1，S2，S3}.Each of candidate stencils has three nodes,as shown in Fig.1.The second order polynomial approximation uk(x)=akx2+bkx+ckcould be obtained by using the function value of u(x)in each candidate stencil,where k=1，2，3.

In each candidate stencil,the formula for uk(x)could be obtained by Taylor series expansion at xjas follows:

When x=Δx/2,Eq.(4)can be approximated by

Then the numerical flux in five-point-stencil WENO scheme could be denoted with weights as

where ω is weight coefficient.

As for CAA problems,it is very important to simulate the short waves with a limited stencil scheme as much as possible.We may equate Eq.(6)and Eq.(1)to yield

Table 1 Different free parameter ω0and λ.

where

Consider the Fourier transform of the left-hand and right-hand sides of Eq.(7).Then using a Taylor expansion,we can find

When function u(x)is discontinuous,smoothness factor has to be defined to change the redistribution of weights,which could guarantee that the node value approximation is essentially non-oscillatory around the discontinuity.Due to the design method proposed by Jiang and Shu6,in this fivepoint-stencil WENO scheme,the smoothness factor βmcould be denoted as

Then the non-linear weights are

In order to avoid the denominator being zero,ε=10-6,m=0，1，2.

By calculating the non-linear weighting parameters,the flux term of Eq.(6)could be denoted as

For numerical stability,modification has been made to the flux term of WENO scheme by the Lax-Friedrichs flux-vector splitting method,and the flux term of Eq.(1)could be denoted as

3.Modified WENO scheme verification

In this section,the comparison of dispersion and dissipation between the modified WENO scheme and the DRP scheme has been conducted with 1D convection equation,in which multiple initial conditions and different weights have been taken into account.

Consider 1D linear convection equation

The initial condition is

values,the dissipation of the scheme would be different.The scheme acts as central when β=0,and acts as upwind when β=1.According to the β definition in Ref.12,the definition of β within[s1，s2]has been modified as shown below:

The grid scale has been set as Δx=0.004,with time step Δt=0.0001 and the number of iterations being 10000.The time marching method used is six-level four-order HALERK6 scheme.13,14Fig.4 shows the comparison of different spatial discretion schemes with analytical solution with grid scale as Δx=0.004.It is also shown in these two figures that:(A)the DRP scheme could induce serious oscillations.Though the oscillation could be attenuated by adding artificial dissipation,the oscillation near the discontinuity is still very severe.(B)All the WENO scheme solutions match the analytical solutions very well,with the accuracy getting better when spatial grid scale is smaller.(C)As for the same WENO scheme,the computational result varies very little while using different flux splitting methods.

4.Interface flux reconstruction

DRP spatial discretion scheme coupled with artificial dissipation is usually applied to the computation of flux at the interface.With the function test results in Part 3,when the magnitude of discontinuity at the interface is unknown,the amount of the artificial dissipation could not be determined with this method.The flux dissipation would be too much when the artificial dissipation is excessive,while it could result in oscillations when the dissipation is inadequate.In this paper,modi fied WENO scheme is applied at the interface,with hybrid Lax-Friedirchs flux-vector splitting method15,16used for flux splitting,and high order DRP scheme has been used inside the grid elements,as shown in Fig.5.F is flux term.The flux at grid node A of the interface of Block 1 could be discreetly solved.The grid node variable of Block 2 could be interpolated from the neighbor grid nodes.

WENO schemes are applied at the interface of blocks,which well keep the numerical stability without damaging the overall computational accuracy.Besides,in order to further improve the computational stability,WENO_opt6 is applied at node C,while WENO_opt4 is applied at node B and WENO_opt1 at node A.

5.Numerical simulation

In order to verify the feasibility of the interface flux reconstruction in multi-block patch/non-patch interface,two study cases have been conducted.

5.1.Gauss impulse propagation

One Gauss impulse sound source has been set at the center of the computational domain,as shown in Fig.6,with the background Mach number as zero.At the initial time t=0,the initial value of the Gauss sound source is set as

where u′and v′are component of fluctuation velocity,ρ′is fluctuation density,p′is fluctuation pressure.

In order to avoid the numerical contamination from the sound wave reflection,Perfect Match Layer(PML)nonreflection boundary condition17has been applied at the farfield.The unified unitless time step dt=0.5 has been used in sound field computation.Three different grid topologies have been used,which are named as Case 1(Δd1=1，Δd2=1),Case 2 (Δd1=1，Δd2=1/2), and Case 3 (Δd1=1/2，Δd2=1).As shown in the right part of Fig.6,it is the grid topology of Case 2 con figuration.

Figs.7 and 8 includes the sound pressure distribution contour of the Gauss impulse at t=100Δt and the sound pressure curve at the location of y=0 station.It is shown in the sound propagation contour that sound wave could propagate uniformly to the space with constant speed in different kinds of patched grids.From the comparison of the computational results and the analytical solutions,it could be concluded that the interface flux reconstruction method could be effectively applied in sound wave propagation.Fig.9 is comparison of sound pressure by using three grid topologies.Compared to Case 1,the calculation results of Case 2 and Case 3 are more closer to analytical solution.

5.2.A Harmonic energy source scattering from multiple cylinders

In order to further prove the simulation capability of the interface flux reconstruction method applied to complex con figuration problems,the periodical point sound source propagation across triple 2D cylinder test case from the NASA 4th CAA workshop18–20has been selected for verification.

As shown in Fig.10,three cylinders of various diameters have been laid out in the computational domain.The left cylinder is at(x=-4，y=0)with diameter being 1,the upper right and lower right cylinder are at (x=3，y=4) and(x=3，y=-4)respectively,with both of the diameters being 0.75.At the center of the computational domain,(x=0，y=0)is a periodical point sound source,the formula of which is as follows:

where S is sound source term.Slide wall boundary has been applied in order to simulate the interference,diffraction and other complex phenomenon of the sound waves,while the PML non-reflection boundary condition17has been used at the far- field.Seven-point four-order DRP scheme is used for spatial discretion,seven-point stencil dissipation scheme is used inside the grid block,the flux reconstruction scheme with modified WENO scheme is applied at the interface,and time advancing method used is the 6-level 4-order Runge-Kutta scheme.21The size of the complete computational domain is(x ∈ [-9，9]，y ∈ [-9，9]),the topology of the grid domain is shown in Fig.10.Two different grid distribution methods have been used in this case.At the interface near the three cylinders,1–1 grid node matching and 1–2 patch are both used,with the grid details at the cylinder as shown in Fig.11.Grid scale Δx=0.025 has been used for non-patch grid sound field simulation,with unified unitless time step Δt=0.002 and the number of iterations as 105.

Figs.12–14 are the sound pressure distribution contour of the periodical sound source propagation procedure at different moments.From the comparison,it could be shown that(A)interference always occur between the sound waves regardless of any grid subdivision methodology,so does the diffraction;(B)the computational results by 1–2 patch grid could obtain higher definition of sound wave compared to the 1–1 node matching results.

In order to quantitatively analyze the capability of interface flux reconstruction in the program,the RMS of the sound pressure distribution on the centroid line y=0 together with that on the cylinder surface is presented,as shown in Figs.15 and 16.From the comparison,it is shown that the peak value of the sound wave computed in this paper is slightly less than that of exact solution.In Fig.16,the sound wave peak value on the right cylinder surface is also slightly less than that of exact solution.The possible reason that causes the discrepancy might be that interface of block had more large dissipation.However,the total trend has shown that the current results well match that of the references,and the patch grid could even obtain better results.

6.Conclusions

For multi-scale grid problems,interface flux reconstruction method based on modi fied WENO scheme has been introduced.Based on the methodology of DRP scheme,the modified WENO scheme has been deduced.

(1)For different weights of the phase errors and amplitude errors,the influence of this scheme on the grid definition and dissipation has been analyzed.Then a hybrid flux vector splitting method with treatment to grid discontinuity has been developed and verified using sine wave function and hybrid wave function.The results have shown that modified WENO scheme could effectively simulate non-discontinuous wave and continuous wave problems with enough grid definition.

(2)According to the characteristic of the dissipation of the modified WENO scheme, flux reconstruction method has been developed using amplitude error accumulation at the discrete grid node on the interface.The feasibility and accuracy of the developed flux reconstruction method applied in multi-scale grid problem have been verified by comparing the numerical result of Gauss impulse sound source radiation with the analytical solutions.

(3)The accuracy of the developed method based on multiscale grid applied to complex con figuration problems has been verified by analysis of triple cylinder interference case.

References

1.Tam CKW,Webb JC.Dispersion-relation-preserving finite difference schemes for computational acoustics.J Comput Phys 1993;107(2):262–81.

2.Lele SK.Compact finite difference schemes with spectral-like resolution.J Comput Phys 1992;103(1):16–42.

3.Tam CKW,Kurbatskii KA.Multi-size-mesh multi-time-step dispersion-relation-preserving scheme for multiple-scales aeroacoustics problems.Int J Comput Fluid Dynam 2014;17(2):119–32.

4.Liu XD,Osher S,Chan T.Weighted essentially non-oscillatory schemes.J Comput Phys 1994;115(1):200–12.

5.Shu CW,Osher S.efficient implementation of essentially nonoscillatory shock capturing schemes.J Comput Phys 1988;77(2):439–71.

6.Jiang GS,Shu CW.efficient implementation of weighted ENO schemes.J Comput Phys 1996;126(1):202–28.

7.Zhang SH,Jiang GS,Shu CW.Improvement of convergence to steady state solutions of Euler equations of Euler equations with the WENO schemes.J Sci Comput 2011;47(2):216–38.

8.Martin MP,Taylor EM,Wu M.A bandwidth-optimized WENO scheme for the effective direct numerical simulation of compressible turbulence.J Comput Phys 2006;220(1):270–89.

9.Martin MP.DNS of hypersonic turbulent boundary layers.Part I:Initialization and comparison with experiments.J Fluids 2007;570(1):347–64.

10.Wu CH,Zhao N,Tian LL.An improved hybrid compact WENO scheme.Acta Aerodynamica Sinica 2013;31(4):477–81[Chinese].

11.Hou ZX,Yi SH,Li H.Analysis and improvement of high precision,high resolution WENO schemes.J Nat Univ Defense Technol 2003;25(1):17–20[Chinese].

12.Lin SY,Hu JJ.Parametric study of weighted essentially nonoscillatory schemes for computational aeroacoustics.AIAA J 2001;39(3):1002–14.

13.Allampalli V,Hixon R,Nallasamy M.High-accuracy large-step explicit Runge-Kutta(HALE-RK)schemes for computational aeroacoustics.J Comput Phys 2009;228(18):3837–50.

14.Calvo M,Franco JM.A new minimum storage Runge-Kutta scheme for computational acoustics.J Comput Phys 2004;201(2):1–12.

15.Kang L,Lee CH.A non- flux-splitting WENO scheme with low numerical dissipation.Sci China Tech Sci 2011;41(4):460–73.

16.Zhu HJ,Yan ZG,Liu HY.Properties of Osher flux with entropy fix with entropy fix in high-order WCNS.Acta Aeronautica et Astronautica Sinica 2017;38(5):1–10[Chinese].

17.Hu FQ,Li XD,Lin DK.Absorbing boundary conditions for nonlinear Euler and Navier-Stokes equations based on the perfectly matched layer technique.J Comput Phys 2008;227(1):4398–424.

18.Dahl MD.Fourth computational aeroacoustics(CAA)workshop on benchmark problems.Washington D.C.:NASA Glenn Research Center;2004.Report No.:NASA/CP-2004-212954.

19.Gao JH.A block interface flux reconstruction method for numerical simulation with high order finite difference scheme.J Comput Phys 2014;241(3):1–17.

20.Gao JH,Yang ZG,Li XD.An optimized spectral difference method for CAA problems.J Comput Phys 2013;231(1):4848–66.

21.Liu L,Li XD,Hu FQ.Non-uniform time-step explicit Runge-Kutta discontinuous Galerkin method for computational aeroacoustics.J Comput Phys 2010;229(19):6874–97.