h-adaptive Discontinuous Galerkin Method for Laminar Compressible Navier-Stokes Equations on Curved Mesh

Sun Qiang,Lyu Hongqiang,Wu Yizhao

College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,P.R.China

h-adaptive Discontinuous Galerkin Method for Laminar Compressible Navier-Stokes Equations on Curved Mesh

Sun Qiang,Lyu Hongqiang*,Wu Yizhao

College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,P.R.China

An h-adaptive method is developed for high-order discontinuous Galerkin methods(DGM)to solve the laminar compressible Navier-Stokes(N-S)equations on unstructured mesh.The vorticity is regarded as the indicator of adaptivity.The elements where the vorticity is larger than a pre-defined upper limit are refined,and those where the vorticity is smaller than a pre-defined lower limit are coarsened if they have been refined.A high-order geometric approximation of curved boundaries is adopted to ensure the accuracy.Numerical results indicate that highly accurate numerical results can be obtained with the adaptive method at relatively low expense.

h-adaptivity;high-order discontinuous Galerkin methods(DGM);N-S equations;high-order boundary approximation

0 Introduction

Discontinuous Galerkin(DG)methods[1-13]have received increasing attention in computational fluid dynamics in recent years due to various attractive features.Bassi and Rebay[3-5]developed a high-order discontinuous finite element method to solve the Euler and Navier-Stokes equations. Cockburn and Shu[6,7]devised a high-order accurate total variation bounded(TVB)Runge-Kutta discontinuous Galerkin(RKDG)method to simulate the nonlinear systems of conservation laws. More recently,high-order DG methods have been applied to solve various engineering problems[8-18]on unstructured grid.In fact,DG methods are similar to finite element methods which can achieve higher-order accuracy via using high-order polynomial approximation inside elements.Moreover,upwind scheme can be easily implemented through using appropriate numerical fluxes over element interfaces.In addition,DG methods lead to compact space discretization formulae for both the Euler and the Navier-Stokes(N-S)equations. The compactness of the methods has advantages for parallel implementation.

Despite these advantages,DG methods still need to be improved in many respects,such as the shock-capturing and the huge computational expense caused by the high-order polynomial approximation[8-15].Usually,high discontinuity only exits locally,e.g.the boundary layer in the flow field.It will cost high expense to capture them by enhancing the order of the polynomials or generating more dense mesh globally.Adaptive DG methods are helpful to solve such problems. Thanks to the simple communication at element interfaces,elements with″hanging nodes″can be treated as elements without hanging nodes,which simplifies local mesh h-refinement.In addition, the communication at element interfaces allows different orders between neighboring elements. Several adaptive DG methods[19-23]have been developed to improve the accuracy and reduce the computational expense.

It has been proved that high-order DG methods are inaccurate at curved solid walls if a piece-wise linear approximation of the geometry of the boundary is employed[3],and a higher-order boundary representation is necessary to ensure the accuracy of the solution.In the paper,an h-adaptive strategy is developed for DG methods to simulate compressible laminar N-S equations on highly accurate boundary.Because of the high intensity of the vorticity in boundary layer and shedding vortex regions,vorticity is used as the sensor of the h-adaptivity.Eor the steady case,a Newton method[24]is employed to solve the nonlinear discrete systems and the Block-Gauss Seidel[10,11]method is used to solve the resulting sparse linear system at each nonlinear iteration.The time integration of the unsteady case presented below can be accomplished by means of an explicit method. The four-stage Runge-Kutta scheme is used in the paper.Since DG methods are relatively sensitive to the initial guess,a hierarchical solution procedure is suggested[10,12].

*Corresponding author,E-mail address:hongqiang.lu@nuaa.edu.cn.

How to cite this article:Sun Qiang,Lyu Hongqiang,Wu Yizhao.h-adaptive discontinuous Galerkin method for laminar compressible Navier-Stokes equations on curved mesh[J].Trans.Nanjing Univ.Aero.Astro.,2016,33(5):566-575.

http://dx.doi.org/10.16356/j.1005-1120.2016.05.566

1 Governing Equations

The two-dimensional laminar N-S equations can be written as follows

where the conservative variables u and the cartesian components fe(u)and ge(u)of the inviscid (Euler)flux function Fe(u)are given by

whereρ,P and e denote the density,pressure and the total internal energy per unit mass respectively.u and v are the velocity components.The total

The cartesian components fv(u,?u)and gv(u,?u)of the viscous flux function Fv(u,?u) are given by

2 DG Discretization

The weak formulation of Eq.(1)can be obtained by multiplying a″test function″v,integrating over the domainΩand performing integration by parts

where F(u,?u)=Fe(u)-Fv(u,?u),?Ωis the boundary ofΩ.

The integrals over the domainΩcan be expanded into the sum of integrals over a collection of non-overlapping triangle elements{E}.The semi-discrete equations for element E can be written as

where?E is the boundary of E.In each element, the functions uhand vh,which are the approximations to u and v,are given by

where the expansion coefficients Uiand Videnotethe degrees of freedom of the numerical solution and the test function in element E,φithe n (shape)basis functions of degree p.Since Eq.(6) must be satisfied for any element E and function vhand vhare a linear combination of n shape functionsφi,Eq.(6)is equivalent to the following system

Note that,no global continuity is enforced on u andφi,discontinuities are allowed over element interfaces.Elux terms are not unique at element interfaces due to the discontinuous function approximation.The physical normal flux F(uh,?uh)·n in Eq.(8)is replaced by a numerical flux H u-h,?u-h,u+h,?u+h,

( )n,which is calculated using the internal u-h,external interface state u+hand the normal vector n pointing outward from E.The numerical flux for the inviscid part of the equations can be analogous to those employed in upwind finite volume methods.In the paper,the LLE scheme is used[10,24].

In the context of the DG method an auxiliary variableθ=?u is introduced for the treatment of the viscous flux.Since DG methods obey to the hyperbolic systems of conservation laws,the following system of two first-order equations is obtained

Similar to the treatment of Eq.(1),the weak formulation of the first equation can be obtained by applying the DG discretization

θcan be written as the following formulation

The numerical flux function Hhincludes the inviscid numerical fluxand the viscous numerical flux function Hv.In the paper,Hvis given by the average value of the viscous fluxes on the interface.

3 Relaxation

By assembling together all the elemental contributions,the semidiscrete equations can be written as

where M is the mass matrix,U the global vector of the degrees of freedom,and() R U the residual vector.Due to the block diagonal structure of M, the time integration of the unsteady problem presented below can be accomplished by means of an explicit method.In the paper,the TVB Runge-Kutta schemes is used[6].

Eor steady problems,the Newton method[23]is used to solve the nonlinear system in Eq.(15)

where w is the under-relaxation factor andΔUnis obtained by solving the following linear system

In order to improve the conditioning of the linear system(17),a pseudo-time derivative is introduced to original discrete system[25]

4 Adaptive Strategy

The entire computation procedure starts from solving p=0(p is the order of the basic functions)solution on a very coarse initial grid.As the order p increases,the mesh needs to be refined in the region where the solution is not smooth enough.After a number of iteration steps,the solution in the region mentioned may become smooth enough.Then,the elements which have been refined should be coarsened to reduce the computational expense when the solution becomes smooth enough.

The vorticity v exists everywhere in the viscous flow.Moreover,due to the great velocity gradient in the shear layer and the vortices,it can be huge in these regions.In the paper,the vorticity v is used as the adaptivity sensor.

4.1 Mesh refinement

During the computational process,elements should be refined when v is larger than the pre-defined upper limit.The″father″e(cuò)lement which needs to be refined is divided into four″child″e(cuò)lements(See Fig.1).

It has been proved that the high-order DG method is inaccurate at curved solid walls if only a piecewise linear approximation of the boundary is employed and a higher-order boundary representa-tion can improve the accuracy of the solution[3]. In the paper,the edges on the solid wall of the boundary elements are represented by a high-order polynomial.The designed high-order(Sixthorder)curve can represent the real solid wall precisely(See the dash lines in Fig.1).

Fig.1 Element refinement

If the boundary elements need to be refined, the mid-point of the boundary edge is found according to the designed curve and the new″child″e(cuò)lements are generated by connecting it with the mid-points of the other two edges as shown in Fig.1.Two of the new″child″e(cuò)lements are on the solid wall and their boundary edge will also be represented by the high-order polynomial(See the dash lines in Fig.1).

In order to avoid significant gradient in mesh size between neighboring elements,a smoothing strategy is employed.If element e in Fig.2 needs to be refined,the neighboring element f must also be refined to avoid extreme difference in local mesh size.In another word,the maximum difference between refinement times of neighboring elements is 1.At the same time,a minimum mesh size hminis pre-defined and the refinement will stop when the element′s genomic size reaches hminto avoid the unlimited refinement of the element.

Fig.2 Smoothing strategy

4.2 Mesh coarsening

During the computation,solution in the elements which have been refined may become smooth enough.In order to reduce the expense, these elements can be coarsened.In this case,the four″child″e(cuò)lements will merge into one″father″e(cuò)lement where the vorticity v is smaller than the pre-defined lower limit.

In Fig.3,the solid lines indicate the elements which are not on the solid wall.If the″father″e(cuò)lement is on the solid wall,its boundary edge will also be represented by a high-order polynomial to approach the real wall(See the dash lines in Fig. 3).Like in the refinement,in order to avoid extreme difference in local mesh size,″child″e(cuò)lements fi(i=0,1,2,3)cannot be merged(See Fig.2)if e is coarsened.In another word,mesh coarsening is the inverse operation of the mesh refinement and the max difference between refinement times of neighboring elements is 1.In addition,the initial element cannot be coarsened.

Fig.3 Element coarsening

4.3 Data storage structure of grid

To ensure the high program transportability, the mesh adaptivity works as an independent module.It only changes the mesh and flow field solver module is not impacted.In the paper,all the information of the points,edges and elements is stored in a list structure.Fig.4 demonstrates the refinement of element E.

Fig.4 Element increasing

In Fig.4,E denotes the″father″e(cuò)lement which needs to be refined.The center″child″e(cuò)lement remains the same index as E,and the other three around the center″child″e(cuò)lement(See Fig.2)range at the end of the list.The same method is applied to the points and edges.Data transmission between grid module and flow field solver module will work well without any influence from mesh adaptivity.

Similarly,Fig.5 shows that four″child″e(cuò)le-ments merge into a″father″e(cuò)lement.

Fig.5 Element decreasing

The three″child″e(cuò)lements fiwill be merged to the center″child″e(cuò)lement E,and accordingly the index of the elements after f3should be subtracted 3.

Two flags are attached to each element to indicate the initial index and refinement times during the mesh refinement.The mesh coarsening will work according to these two flags and the relationship between neighboring elements.

5 Numerical Results

5.1 Transonic flow around NACA0012 airfoil

Eirstly,the subsonic viscous flow around the NACA0012 airfoil(Ma=0.5,α=0°,Re=5 000) is simulated.The initial mesh contains 478 elements,260 grid points,and there are 32 grid points on the solid boundary(See Fig.6).

Fig.7 demonstrates the Mach contours and the vorticity distribution obtained when p=4 on the initial grid.It is obvious that the solution is not smooth enough because of the coarse grid in the boundary layer,which suggests smaller mesh size in this region to improve the accuracy of the solution.

In order to improve the accuracy of the solution and reduce the computational expense,the local adaptive method introduced above is applied. Fig.8(a)shows the final mesh and the vorticity v distribution,where only the local mesh near the boundary and in wake region is refined.Fig.8(b) shows the final local mesh after adaption compared with the initial mesh.

Fig.9 depicts the Mach contours obtained with the adaptive method,where the solution is much smoother than that obtained on the initial mesh in Fig.7.In addition,Fig.9(b)shows the streamlines and the two symmetrical vortices[3]in the wake region are captured well.

Fig.6 Initial grid

Fig.8 Vorticity distribution and local grid after adaption

Fig.7 Mach contours and vorticity v distribution when p=4

Fig.9 Mach contours and streamlines after adaption

Fig.10 Cpand Cfdistribution

5.2 Von Karman vortex street

The well-known Von Karman vortex street is simulated,where Ma=0.1,α=0,Re=150.The initial mesh is shown in Fig.11,which contains 1 114 elements and there are only 12 points on the solid boundary(See Fig.11).

Fig.12 demonstrates the Von Karman vortex street obtained when p=4 on the initial grid.The solution is not smooth and the resolution of the vortices is low because the initial mesh is not fine enough.Unfortunately,shedding vortices are not captured precisely on the initial coarse mesh even if the high-order scheme is applied.On the other hand,the intensity of the vortices is low due to the big numerical dissipation which is mainly caused by the large mesh size.It is suggested that the mesh in these regions should be refined.

Adaptive method introduced above is used in this case and the final local grid and Von Karman vortex street are shown in Fig.13.Due to the smaller mesh size in the boundary layer region andwake region after adaption,the solution is much smoother compared with that in Fig.12 and the intensity of the vortices is enhanced because of the low dissipation.

Fig.11 Initial grid

Fig.12 Von Karman vortex street on initial mesh

Fig.13 Von Karman vortex street with adaptive method

Fig.14 Von Karman vortex street with mesh adaptivity

In the paper,refinement and coarsening always work simultaneously.The elements where the vorticity is greater than the pre-defined upper limit are refined,and those where the vorticity is smaller than the pre-defined lower limit are coarsened to reduce the computation and storage expense if they have been refined.Fig.14 shows the Von Karman vortex street with the adaptive method introduced above on several time of one period.

The Von Karman vortex street is a quasisteady case because of its periodicity.The number of elements in the entire domain should vary in a small scale for which it is a good case to test the behavior of the introduced method.Fig.15 shows the history of the number of the elements over the entire domain,where the element number varies around 2 500.

The evolution of the drag and lift coefficient in time is shown in Fig.16 while the period of vortex shedding(Strouhal number)is found to be St=0.185.In Table 1,the variations of lift coefficient Cland drag coefficient Cdare documented with amplitudes and mean values and the results of sixth-order finite difference scheme are also included as a reference[26]for comparison.It is clearly that the accuracy of lift and drag coefficient is drastically improved by using the adaptive method.

Fig.15 Number of element during iteration

Fig.16 Variation of lift and drag coefficient

Tab.1 Comparison between lift and drag coefficient

6 Conclusions

An h-adaptive DG method is developed to solve the two-dimensional compressible laminar N-S equations.The vorticity v works well as the sensor of the mesh adaptivity in the subsonic viscous flow cases.In order to ensure the accuracy of the solution,a high-order approximation boundary is designed to approach the real solid wall during the h-adaptivity.Numerical results show that grid refinement and coarsening only work in the local region and the accuracy of the solution is improved at low expense.

Acknowledgement

This work was supported by the National Natural Science Eoundation of China(11272152).

[1] REED W H,HILL T R.Triangular mesh methods for the neutron transport equation:Los Alamos Report LA-UR-73-479[R].1973.

[2] LESAINT P,RAVIART P A.On a finite element method for solving the neutron transport equation[J]. Mathematical Aspects of Einite Elements in Partial Differential Equations,1974(33):89-123.

[3] BASSI E,REBAY S.A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[J]. Journal of Computational Physics,1997,131(2): 267-279.

[4] BASSI E,REBAY S.High-order accurate discontinuous finite element solution of the 2D Euler equations [J].Journal of Computational Physics,1997,138 (2):251-285.

[5] BASSI E,REBAY S.A high order discontinuous Galerkin method for compressible turbulent flows [M].[S.l.]:Springer Berlin Heidelberg,2000:77-88.

[6] COCKBURN B,SHU C W.The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems[J].Journal of Computational Physics,1998,141(2):199-224.

[7] COCKBURN B,SHU C W.The local discontinuous Galerkin method for time-dependent convection-diffusion systems[J].SIAM Journal on Numerical Analysis,1998,35(6):2440-2463.

[8] LU H,BERZINS M,GOODYER C E,et al.Adaptive high-order discontinuous Galerkin solution of elastohydrodynamic lubrication point contact problems[J].Advances in Engineering Software,2012, 45(1):313-324.

[9] LU H,XU Y,GAO Y,et al.A CED-based high-order discontinuous Galerkin solver for three dimensional electromagnetic scattering problems[J].Advances in Engineering Software,2015,83:1-10.

[10]LU H,SUN Q.A straight forward hp-adaptivity strategy for shock-capturing with high-order discontinuous Galerkin methods[J].Advances in Applied Mathematics and Mechanics,2014,6(1):135-144.

[11]LU H,SUN Q,QIN W L.3D numerical solution of aero-noise with high-order discontinuous Galerkin method[J].Transactions of Nanjing University of Aeronautics&Astronautics,2013,30(3):227-231.

[12]LV Hongqiang,ZHU Guoxiang,SONG Jiangyong, et al.High-order discontinuous Galerkin solution of linearized Euler equations[J].Chinese Journal of Theoretical and Applied Mechanics,2011,43(3): 621-624.(in Chinese)

[13]LV H Q.High-order discontinuous Galerkin solution of low-Re viscous flows[J].Modern Physics Letters B,2009,23(3):309-312.

[14]LV H Q,XU Y,GAO Y,et al.A high-order discontinuous Galerkin method for the two-dimensional time-domain Maxwell′s equations on curved mesh [J].Advances in Applied Mathematics&Mechanics,2016,8(1):104-116.

[15]LV H Q,CAO K,WU L B Y,et al.High-order mesh generation for discontinuous Galerkin methodsbased on elastic deformation[J].Advances in Applied Mathematics&Mechanics,2016,8(4):693-702.

[16]Zhang Laiping,Liu Wei,He Lixin,et al.A class of discontinuous Galerkin/finite volume hybrid schemes based on the″static re-construction″and″dynamic reconstruction″[J].Chinese Journal of Theoretical and Applied Mechanics,2010,42(6):1013-1022.(in Chinese)

[17]YU J,YAN C.An artificial diffusivity discontinuous Galerkin scheme for discontinuous flows[J].Computers&Eluids,2013,75:56-71.

[18]YU Jian,YAN Chao.Study on discontinuous Galerkin method for Navier-Stokes equations[J].Chinese Journal of Theoretical and Applied Mechanics,2010, 42(5):962-969.(in Chinese)

[19]YANG X,JAMES A J,LOWENGRUB J,et al.An adaptive coupled level-set/volume-of-fluid interface capturing method for unstructured triangular grids [J].Journal of Computational Physics,2006,217 (2):364-394.

[20]REMACLE J E,LI X,SHEPHARD M S,et al.Anisotropic adaptive simulation of transient flows using discontinuous Galerkin methods[J].International Journal for Numerical Methods in Engineering,2005, 62(7):899-923.

[21]XU Yun,YU Xijun.Adaptive discontinuous Galerkin methods for hyperbolic conservation laws[J].Chinese Journal of Computational Physics,2009,26(2):159-168.(in Chinese)

[22]WU Di,YU Xijun.Adaptive discontinuous Galerkin method for Euler equations[J].Chinese Journal of Computational Physics,2010,27(4):492-500.(in Chinese)

[23]HARTMANN R,HOUSTON P.Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations[J].Journal of Computational Physics,2002,183:508-532.

[24]XIA Yidong,WU Yizhao,LV Hongqiang,et al. Parallel computation of a high-order discontinuous Galerkin method on unstructured grids[J].Acta Aerodynamica Sinica,2011,29(5):537-541.(in Chinese)

[25]HILLEWAERT K,CHEVAUGEON N,GEUZAINE P,et al.Hierarchic multigrid iteration strategy for the discontinuous Galerkin solution of the steady Euler equations[J].International Journal for Numerical Methods in Eluids,2006,51(9/10):1157-1176.

[26]INOUE O,HATAKEYAMA N.Sound generation by a two-dimensional circular cylinder in a uniform flow[J].Journal of Eluid Mechanics,2002,471: 285-314.

Mr.Sun Qiang received B.S.degree from College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics in 2011.In September 2011,he became a post-graduate student at the Aerodynamics Department. His research is focused on the discontinuous Galerkin method and mesh adaptivity.

Prof.Lv Hongqiang received B.S.degree from Aerodynamics Department of Nanjing University of Aeronautics and Astronautics in 1999 and M.S.degree in 2002 from the same university.He received Ph.D.and Doctor degrees from University of Leeds in 2006.In the same year,he joined in College of Aerospace Engineering of Nanjing University of Aeronautics and Astronautics.His current research interest includes the application of the discontinuous Galerkin methods and numerical simulation of the Maxwell equations.

Prof.Wu Yizhao received B.S.degree from University of Science and Technology of China in 1968,and received M. S.degree in 1981 and Ph.D.degree in 1987 from Nanjing Aviation Institute.He is currently a professor and doctoral supervisor at College of Aerospace Engineering of Nanjing University of Aeronautics and Astronautics,and his research interests are computational fluid dynamics and multigrid.

(Executive Editor:Xu Chengting)

V211.3 Document code:A Article ID:1005-1120(2016)05-0566-10

(Received 14 April 2015;revised 18 August 2015;acceped 5 September 2015)