Qin Wnglong,Lyu Hongqing,*,Wu Yizho,Zhou Shijie,Chen Zhengwu

aDepartment of Aerodynamics,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China

bChina Aerodynamics Research and Development Center,Mianyang 621000,China

Implicit discontinuous Galerkin method on agglomerated high-order grids for 3D simulations

Qin Wanglonga,Lyu Hongqianga,*,Wu Yizhaoa,Zhou Shijiea,Chen Zhengwub

aDepartment of Aerodynamics,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China

bChina Aerodynamics Research and Development Center,Mianyang 621000,China

High quality of geometry representation is regarded essential for high-order methods to maintain their high-order accuracy.An agglomerated high-order mesh generating method is investigated in combination with discontinuous Galerkin(DG)method for solving the 3D compressible Euler and Navier-Stokes equations.In this method,a fine linear mesh is first generated by standard commercial mesh generation tools.By taking advantage of an agglomeration method,a quadratic high-order mesh is quickly obtained,which is coarse but provides a high-quality geometry representation,thus very suitable for high-order computations.High-order discretizations are performed on the obtained grids with DG method and the discretized system is treated fully implicitly to obtain steady state solutions.Numerical experiments on several flow problems indicate that the agglomerated high-order mesh works well with DG method in dealing with flow problems of curved geometries.It is also found that with a fully implicit discretized system and a p-sequencing method,the DG method can achieve convergence state within several time steps which shows significant eff iciency improvements compared to its explicit counterparts.

1.Introduction

Discontinuous Galerkin(DG)methods have received growing interests in the past two decades for their various attractive features,such as high-order accuracy,great geometry flexibility,straightforward implementation of h/P adaptation and parallel computing.1–5The same as the classical finite element methods,DG methods can achieve high-order accuracy on grids by means of high-order polynomial approximation within elements while the physics of wave propagation is accounted for by means of solving Riemann problems at element interfaces as in upwind finite volume methods.

Despite these advantages,there are still some limitations of the DG methods.Particularly,it remains an open issue for generating meshes which can robustly produce accurate approximation of curved geometries.Bassi and Rebay showed that in the case of curved geometries,a high-order DG discretization with straight-sided two dimensional elements yields low order accurate and even physically wrong results,4and similar results are also proposed in publications.6–8In order to improve the quality of geometry representation,a high-order mesh is required.Landmann et al.9proposed a curved boundary representation method with local splining surface on triangles and quadrangles,which showed a good performance.The method was further extended to hybrid mesh and turbulent flows.10,11With the aid of CAD tools,such as non-uniform rational B-splines(NURBS)12and computational analysis programming interface(CAPRI)13,a more accurate boundary representation was achieved.Agglomeration algorithm isconsidered asan alternative efficient approach to provide high-order approximations of curved geometries.Combined with DG method,Bassi et al.14developed an agglomeration based discontinuous Galerkin algorithm for two-dimensional cases.Based on the multilevel algorithm library MGridGen,coarse polyhedral grids are obtained with the original fine grids,which is straightforward and robust.Compared with the results on grids with the same number of degrees offreedom,more accurate solutions have been obtained on coarse agglomerated grids due to its more accurate representation of the curved geometries.However this work is published only for two-dimensional cases and the library MGridGen is necessary as a pre-processing tool.

Efficient time integration schemes for discontinuous Galerkin method have been a research focus recently.Explicit timeintegration methods such as multistage Runge-Kutta schemes15–17have long been used in discontinuous Galerkin method for their advantages of easy implementation,low memory storage requirement and easy for parallelization.However for large-scale fluid problems,especially for solutions of the Navier-Stokes equations,the convergence rate slows down dramatically.In order to speed up convergence,a lot of implicit approaches have been developed.These range from Gauss-Seidel to Krylov subspace methods with a variety of preconditioners.18–21

Inspired by the agglomeration idea,an agglomerated highorder mesh generation method is investigated for generating three-dimensional coarse grids with fine boundary representations.This method has several attractive features:first,the surfaces are represented with the information of the original fine grids,which are accurately represented without any CAD aid;Secondly,this method can be implemented straightforwardly since no extra libraries are needed.Furthermore,this method can be directly extended for generating high-order highly stretched meshes without causing any cross-over problems,13which can be further used for turbulent flow computations.In order to speed up the convergence,a fully implicit temporal discretization technique is investigated,which demonstrates a good efficiency combined with the DG method.

2.Formulation of DG method

2.1.Governing equations

The Navier-Stokes equations governing unsteady compressible viscous flow can be expressed as

where the summation convention has been used.The conservative variable vector U,inviscid flux vector F and viscous flux vector G are defined by

where ρ,p and e denote the density,static pressure,and total energy per unit mass,respectively.uiis the velocity of the flow in the coordinate direction xiand δijis the Kronecker tensor.The components of the viscous stress tensor σijand the heat flux vector qjare given by

where T is the temperature.With the state equation for perfect gas,the system of equations can be completed,

In the above equations,μ denotes the molecular viscosity,γ the ratio of the specific heats and,Pr is the dimensionless Prandtl number,which takes 0.72 for air.Neglecting viscous effects,the left-hand-side of Eq.(1)represents the Euler equations governing unsteady compressible inviscid flows.22–26

2.2.DG discretization

With mixed formulation9,Eq.(1)can be reformulated as

where s is an introduced auxiliary variable,and fc（U）and fv（U，s）are the inviscid and viscous flux tensor.By introducing suitable piecewise polynomial test function φhand approximate solutions Uhat each cell K,we obtain the weak formulations:

where Ω is the domain, σ the boundary of Ω.The variables with superscript ‘–” denote the values inside the element K and n is the outward unit normal vector to the boundary.fc（Uh）·n and fv（Uh，sh）·n are the inviscid and viscous numerical fluxes,respectively.The inviscid numerical fluxes fc（Uh）·n can be handled in the same way as in the finite volume(FV)method,andhereweusethewellknownlocalLax-Friedrichs(LLF)27flux or Roe flux28for the approximate solution of the Riemann problem.As for the viscous numerical flux,many dealing approaches have been proposed.2,5,29,30The well-known BR231scheme is used in this paper.

2.3.Time integration

The weighted form of the DG discretized equations Eqs.(6)and(7)can be rewritten as an ordinary differential equation form:

where M,u and R（u）denote the mass matrix,the global solution vector and the residual vector.Using the backward Euler difference method,the following linear system is obtained:

where Δu=un+1-unand A= （M/Δt+ ∂R/∂u）.Since the matrix A is always a large sparse block matrix,it can be partitioned into the block diagonal D,upper U and lower L matrices:

Here we use the well-known block Gauss-Seidel method19,32to solve the linear systems.

where the superscript k+1 is the inner iteration index,and Δu0=0.We apply Eq.(12)iteratively until a given convergence tolerance or a maximum number of sweeps is reached.

3.Mesh agglomeration approach

For the traditional second-order accurate flow solvers,a good result can be achieved on a fine linear mesh and the surfaces of geometry are typically represented by a series offine linear face elements.However for DG methods,as has been shown in publication,4this seems not working well even with a lot of elements.The reason is that the linear face elements cannot provide accurate normal directions on the quadrature points which are very important in the computations using DG methods.Besides,there will be errors in the computations of the areas of the curved face elements.Instead,the mesh agglomeration method introduced below can be regarded as a local reconstruction of curved surfaces which can provide more accurate results(see Fig.1).

As illustrated in Fig.2,with eight hexahedral elements agglomerated,a new 27-node quadratic hexahedral element can be obtained.In this method,a structured block of elements with three direction index(i,j,k)is first outputted from commercial mesh generation tools.Since the numbers of nodes in each direction are known,the new elements can be defined easily with 27 control nodes which can be categorized into 4 groups:(A)corner vertices,when all the three direction indexes are even numbers;(B)edge vertices,when only one of the three direction indexes is odd number;(C)face vertices,when only two direction indexes are odd numbers;and(D)one volume vertex,when all the three direction indexes are odd numbers.

Fig.1 Outward unit normal vector n on different face elements.

Fig.2 Illustration of mesh agglomeration approach.

Cubic or higher order elements can also be achieved using this method.Quadratic elements are used in this paper because a curved surface can be represented well with 9 control points.For more complex geometries,more points and elements can be used.

Lagrange polynomials are used to transform between the numerical and physical space with the 27 control points within each element:

where xi,yi,ziare the coordinates of the control points,and Nijkis the Lagrange basis function which can be expressed with tensor product of 1D Lagrange basis functions φi（ξ）:

In order to test the capacity of the mesh agglomeration method in representing curved geometries,a case of compressible inviscid flow past a sphere is used(Ma∞=0.5).A fine linear mesh with 17280 hexahedral elements is first generated and is agglomerated into a quadratic mesh with 2160 hexahedral elements(Fig.3).Compressible Euler equations are then solved on the agglomerated quadratic mesh and the results are compared with that on 17280 linear elements.Fig.4 illustrates the computed absolute number of entropy production distributions(3rd order,P=2)and the entropy production is defined as

Fig.3 Mesh used for computing flow past a sphere.

Fig.4 Computed entropy production of simulating flow past a sphere.

where S=p/ργis the entropy.It can be seen that the number of entropy production is very large when using linear elements although it contains the same mesh information with the agglomerated coarse grids.Actually due to its inaccurate representation of the boundary surfaces,the residual oscillates which in comparison shows the superiority of the mesh agglomeration method.Besides,since the number of the computational agglomerated elements is 1/8 of the original linear elements,the computational time on the agglomerated grids is nearly 1/8 of that on the linear elements for each Newton iteration.Fig.5 gives a comparison of the pressure coefficient Cpon the symmetry plane.The result on a mesh of 13248 agglomerated quadratic elements is computed as a reference resolution.Consistent with the results of the entropy production,the pressure coefficient on the agglomerated grids matches well with the reference data and the data on the fine linear elements oscillates a lot,especially in the wake of the sphere.

4.Presentation of results

In this section,several numerical examples,including both inviscid and viscous test cases,are given to demonstrate the capacity of the presented high-order(up to P=3)implicit DG method on the agglomeratedquadratichexahedral elements.

Fig.5 Comparison of computed pressure coefficient on symmetry plane.

Fig.6 Computed results of simulating laminar flow past a flat plate.

4.1.Laminar flow past a flat plate

Fig.7 Illustration of boundary conditions and computed rates of convergence for subsonic flow through a channel with a bump.

Fig.8 Mach contours and pressure coefficient distribution for subsonic flow through a channel with a bump.

Fig.9 Logarithmic density residual vs time step(P=0–3)and CPU time(P=2)for subsonic flow through a channel with a bump.

A laminar flow over an adiabatic flat plate is first chosen to illustrate the accuracy of our DG code,as the classical Blasius solution can be used to measure the accuracy of the numerical solution.Agglomerated high-order mesh is used in this case although there are no curved geometries.The free-stream flow condition is given as Ma∞=0.5 and Re∞=100000.The computational space[-0.5，1]× [0，1]× [0，0.1]is discretized with（20+40）×16×1 agglomerated quadratic hexahedral elements,with 20×16×1 elements ahead of the flat plate and 40×16×1 elements for the flat plate.The computed profilesof the velocity components ux=u/u∞,in which x,y are coordinates and u,v are velocities in the cor responding directions.Rex=u∞x/ν is the Reynolds number and ν is the kinematic viscosity.It can be seen that an accurate prediction of velocity pro files is obtained,not only the uxprofile but also the vypro file,which is often regarded difficult to compute accurately.The computed skin friction coefficient Cfdistribution along the flat plate is shown in Fig.6(b)with logarithmic plot,which shows a good agreement with the Blasius solution,33indicating the high accuracy of the DG method.

4.2.An inviscid flow past a channel with a smooth bump

Fig.10 Mesh used for computing laminar flow past a sphere.

This test case is chosen to demonstrate the convergence accuracy and efficiency of the presented implicit DG method for solving internal flows on agglomerated high-order mesh.The problem underconsideration isa subsonicflow with Ma∞=0.5.The length,width and height of the channel are 3,0.5 and 0.8,respectively.The smooth bump is on the lower surface of the channel with a function de fined as 0.0625e-25x2from x=-1.5 to x=1.5.Four successively re fined quadratic grids are used to obtain the convergence rate,with point distribution of 6×3×2 for the coarse grids,11×5×3 for the medium grids,21×9×5 for the fine grids and 41×17×9 for the finest grids,respectively.Since this is also an isentropic flow problem,the L2-norm of the entropy production is used as the error measurement:

The definition of the entropy production is in Eq.(14).Fig.7 illustrates the used boundary conditions and the computed order of accuracy.As can be seen,the implicit DG method can achieve the designed order of O（hP+1）on the agglomerated quadratic grids.Fig.8 illustrates the computed mach contours on the finest mesh and the pressure coefficient distribution on xOy plane,which matches well with the symmetric property of the flow.Fig.9(a)shows the logarithmic density residual versus time step for the implicit method.With a P-sequencing method,33in which lower order solutions are used as the initial guess for the higher order solutions,a residual of 10-10can be obtained within several time steps.Fig.9(b)gives the logarithmic density residual versus the CPU time of P=2 solution.Both the explicit Runge-Kutta method and implicit Newton method are employed.It can be seen that the convergence history of the explicit method stalls somewhere,which in comparison demonstrate the efficiency of the fully implicit method.

Fig.11 Computed Mach number contours and streamlines on symmetry plane for laminar flow past a sphere.

4.3.Laminar flow past a sphere

A viscous flow past a sphere is considered in this case.The computation is performed at the Mach number of 0.5 and the Reynolds number of 118 based on the diameter of the sphere.This test case is chosen to demonstrate the capacity of the current studied method for the computation of viscous flows around a simple but highly-curved geometry.Since the model and flow are symmetric,we consider only half of the con figuration.Fig.10 shows the original linear fine grids with 33728 hexahedral elements and the resulting 4216 quadratic hexahedral elements after agglomeration.The non-slip boundary surface is represented with 48 quadratic face elements.The computed Mach number contours in the flow field and the velocity streamlines on the symmetry plane(P=3)are displayed in Fig.11,which agrees well with the results in Ref.19on 119390 tetrahedral elements using an reconstructed discontinuous Galerkin(RDG)method.19Fig.12 gives the logarithmic density residual versus time step and CPU time of the explicit and implicit method(P=2)for this viscous flow.Similar to the results of the inviscid cases,the implicit method used for laminar flow is also much more efficient than the explicit method.

4.4.An inviscid flow past an analytical 3D body of revolution

A subsonic flow past an analytical 3D body of revolution is considered.The computation is performed at a Mach number of 0.5 with an angle of attack 1°.This is a benchmark used in the first international workshop on high-order CFD method and it is aimed at testing high-order methods for the computation of external flow with a high-order curved boundary representation.As illustrated in Fig.13,three successively refined hexahedral grids are used in this case.The numbers of quadratic elements for the coarse,medium and fine grids are 6144,49152 and 393216,respectively.They are agglomerated from the linear meshes which contain 49152,393216 and 3145728elements.As the reference lift coefficient and drag coefficient are given as CLref=1.2201×10-4and CDref=1.298731×10-6,32the plot of CLand CDerror versus grid size h are computed to study the numerical order of accuracy(Fig.14).As expected,the designed orderof O（hP+1）are obtained onthe agglomerated grids.

4.5.Transonic flow past ONERA M6 wing

Fig.12 Logarithmic density residual vs time step and CPU time for laminar flow past a sphere.

Fig.13 Three successively refined quadratic meshes used for computing an inviscid flow past a 3D body of revolution.

Fig.14 Computed lift coefficient error and drag coefficient error vs grid size for inviscid flow past a 3D body of revolution.

Fig.15 Computing mesh and computed pressure contour offlow over an ONERA M6 wing.

A transonic flow over the ONERA M6 wing is considered to access the performance of the mesh agglomeration approach in representing the curved geometry of an aerospace configuration.The free-stream flow condition is given as Ma∞=0.699 and angle of attack α =3.06°.The mesh used in this computation consists of 51636 quadratic hexahedrons which are agglomerated from a fine linear mesh of 413088 hexahedrons(Fig.15(a)).Relatively more points are distributed around the leading edge for capturing the weak shock.The computed pressure contours(P=2)on the wing surface are shown in Fig.15(b),which are smooth except the weak discontinuities in the leading edge.Fig.16 shows the computed pressure coefficient distributions at six spanwise locations on the wing surface.The pressure coefficients are computed at the two nodes of each quadrilateral that intersect with the cutting plane and they are plotted by connecting the data of the two nodes with a straight line,which re flects the discontinuous nature of the DG solution truly.It can be seen from Fig.16 that the computed pressure coefficient distributions show a good agreement with the experimental data although there are some over-and under-shoots in the vicinity of the shock waves.

5.Conclusions

(1)An agglomerated mesh generation method is developed in this paper,which can produce a coarse mesh with relatively accurate representations of the curved geometries.

Fig.16 Computed pressure coefficient distributions at six spanwise locations for inviscid transonic flow over ONERA M6 wing.

(2)Convergence studies indicate that the DG solver has achieved the expected order of accuracy on the agglomerated quadratic grids for complex curved geometries.With a fully implicit discretized system and a psequencing method,DG method can quickly achieve the convergence state.

(3)Numerical experiments are conducted on flows past typical three-dimensional aerospace configurations and the computational results are in good agreement with the published data,which shows the clear superiority of the DG solver on the agglomerated high-order mesh.

Further work will be carried out including extending the agglomeration based high-order mesh method for more complex configurations and turbulent flow simulations.

Acknowledgments

This study was co-supported by the Aeronautical Science Foundation of China(No.20152752033),the National Natural Science Foundation of China(No.11272152)and the Open Project of Key Laboratory of Aerodynamic Noise Control.

1.Cockburn B,Shu CW.The Runge-Kutta discontinuous Galerkin method for conservation laws V:multidimensional systems.J Comput Phys 1997;141(2):199–224.

2.Cockburn B,Shu CW.The local discontinuous Galerkin method for time-dependent convection-diffusion systems.SIAM J Numer Anal 1998;35(6):2440–63.

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

4.Bassi F,Rebay S.High-order accurate discontinuous finite element solution of the 2d Euler equations.J Comput Phys 1997;138(2):251–85.

5.Arnold DN,Brezzi F,Cockburn B,Marini LD.Uni fied analysis of discontinuous Galerkin methods for elliptic problems.SIAM J Numer Anal 2002;39(5):1749–79.

6.Lubon C,Keßler M,Wagner S.A parallel CFD solver using the discontinuous Galerkin approachHigh performance computing in science and engineering.New York:Springer;1986.p.87–109.

7.Krivodonova L,Berger M.High-order accurate implementation of solid wall boundary conditions in curved geometries.J Comput Phys 2006;211(2):492–512.

8.Hindenlang FJ,Gassner GJ,Munz CD.Improving the accuracy of discontinuous Galerkin schemes at boundary layers.Int J Numer Meth Fluids 2014;75(6):385–402.

9.Landmann B,Kessler M,Wagner S,Kramer E.A parallel,highorder discontinuous Galerkin code for laminar and turbulent flows.Comput Fluids 2008;37(4):427–38.

10.Qin WL,Lu HQ,Wu YZ.High-order discontinuous Galerkin solution of N-S equations on hybrid mesh.Chin J Theor Appl Mech 2013;45(6):987–91[Chinese].

11.Qin WL,Lu HQ,Wu YZ.Discontinuous Galerkin solution of RANS equations on curved mesh.Acta Aerodyn Sin 2014;32(5):581–6[Chinese].

12.Hughes TJR,Cottrell JA,Bazilevs Y.Isogeometric analysis:CAD,finite elements,NURBS,exact geometry and mesh re finement.Comput Meth Appl Mech Eng 2005;194(39):4135–95.

13.Wang L,Anderson WK,Erwin T,Kapadia S.Solutions of highorder methods for three-dimensional compressible viscous flows.Report No:AIAA-2012-2836.Reston:AIAA;2012.

14.Bassi F,Botti L,Colombo A,Rebay S.Agglomeration based discontinuous Galerkin discretization of the Euler and Navier-Stokes equations.Comput Fluids 2012;61:77–85.

15.Jameson A.Solution of the Euler equations for two dimensional transonic flow by a multigrid method.Appl Math Comput 1983;13(3):327–55.

16.Spiteri RJ,Ruuth SJ.A new class of optimal high-order strongstability-preserving time discretization methods.SIAM J Numer Anal 2002;40(2):469–91.

17.Zhu HQ,Cheng Y,Qiu JX.A comparison of the performance of limiters for Runge-Kutta discontinuous Galerkin methods.Adv Appl Math Mech 2013;5(3):365–90.

18.Wang ZJ.High-order methods for the Euler and NavierStokes equations on unstructured grids.Prog Aerosp Sci 2007;43(1):1–41.

19.Jiang ZH,Yan C,Yu J.Implicit high-order discontinuous Galerkin method with HWENO type limiters for steady viscous flow simulations.Acta Mech Sinica 2013;29(4):526–33.

20.Jiang ZH,Yan C,Yu J.High-order implicit discontinuous Galerkin schemes for unsteady compressible Navier-Stokes equations.Chin J Aeronaut 2014;27(6):1384–9.

21.Zhang LP,Li M,Liu W,He X.An implicit algorithm for highorder DG/FV schemes for compressible flows on 2D arbitrary grids.Commun Comput Phys 2015;17(1):287–316.

22.Lu HQ,Sun Q,Qin WL.3D numerical solution of aeronoise with high-order discontinuous Galerkin method.Trans Nanjing Univ Aeronaut Astronaut 2013;30(3):227–31.

23.Luo H,Luo L,Nourgaliev R.A reconstructed discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids.J Comput Phys 2010;229(19):6961–78.

24.Xia Y,Luo H,Nourgaliev R.An implicit Hermite WENO reconstruction-based discontinuous Galerkin method on tetrahedral grids.Comput Fluids 2014;96(13):406–21.

25.Xia Y,Luo H,Nourgaliev R.An implicit reconstructed discontinuous Galerkin method based on automatic differentiation for the compressible flows on tetrahedral grids.Report No:AIAA-2013-0687.Reston:AIAA;2013.

26.Lu HQ,Xu YD,Gao YK,Sun Q,Qin WL.A CFD-based highorder discontinuous Galerkin solver for three dimensional electromagnetic scattering problems.Adv Eng Software 2015;83:1–10.

27.Toro EF,Spruce M,Speares W.Restoration of the contact surface in the HLL-Riemann solver.Shock Waves 1994;4(1):25–34.

28.Roe PL.Approximate Riemann solvers,parameter vectors,and difference schemes.J Comput Phys 1981;43:357–72.

29.Peraire J,Persson PO.The compact discontinuous Galerkin(CDG)method for elliptic problems.SIAM J Sci Comput 2007;30(4):1806–24.

30.Yu J,Yan C.Study on discontinuous Galerkin method for Navier-Stokes equations.Chin J Theor Appl Mech 2010;42(5):62–70[Chinese].

31.Bassi F,Crivellini A,Rebay S,Savini M.Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes and k-ω turbulence model equations.Comput Fluids 2005;34(4):507–40.

32.Lu HQ,Sun Q.A Straight forward hp-adaptivity strategy for shock-capturing with high-order discontinuous Galerkin methods.Adv Appl Math Mech 2014;6(1):135–44.

33.Boyd John P.The Blasius function in the complex plane.Exp Math 1999;8(4):381–94.

Qin Wanglong received the B.S.degree in aerospace engineering from Nanjing University of Aeronautics and Astronautics in 2011 and then became a Ph.D.candidate majoring in fluid mechanics there.His research interests are high-order numerical methods of the CFD.

Lyu Hongqiang received his Ph.D.degree from University of Leeds,U.K.,in 2006 and now he is a professor in Nanjing University of Aeronautics and Astronautics.His current research interests include high-order numerical methods of CFD,computational aero-acoustic and computational electromagnetic.

Wu Yizhao is a professor and Ph.D.supervisor at College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics.He received the Ph.D.degree from the same University in 1987.

Zhou Shijie received the B.S.degree in Aerospace Engineering from Nanjing University of Aeronautics and Astronautics in 2013 and then became a postgraduate majoring in fluid mechanics.

29 September 2015;revised 30 January 2016;accepted 6 May 2016

Available online 25 October 2016

Agglomeration;

Discontinuous Galerkin(DG);

High-order;

Implicit scheme;

Navier-Stokes equations

Ⓒ2016 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is anopenaccessarticleundertheCCBY-NC-NDlicense(http://creativecommons.org/licenses/by-nc-nd/4.0/).

*Corresponding author.

E-mail addresses:qinwanglong@126.com(W.Qin),hongqiang.lu@nuaa.edu.cn (H. Lyu), wyzao@nuaa.edu.cn (Y. Wu),523183133@qq.com(S.Zhou).

Peer review under responsibility of Editorial Committee of CJA.

Production and hosting by Elsevier

http://dx.doi.org/10.1016/j.cja.2016.10.004

1000-9361Ⓒ2016 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.

This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).