，，*

1.College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，P.R.China；2.Key Laboratory of Aerodynamic Noise Control，China Aerodynamics Research and Development Center，Mianyang 621000，P.R.China；3.State Key Laboratory of Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，P.R.China（Received 3 June 2019；revised 11 November 2019；accepted 19 November 2019）

Abstract:For the numerical simulation of compressible flows，normally different mesh sizes are expected in different regions.For example，smaller mesh sizes are required to improve the local numerical resolution in the regions where the physical variables vary violently（for example，near the shock waves or in the boundary layers）and larger elements are expected for the regions where the solution is smooth.h-adaptive mesh has been widely used for complex flows.However，there are two difficulties when employing h-adaptivity for high-order discontinuous Galerkin（DG）methods.First，locally curved elements are required to precisely match the solid boundary，which significantly increases the difficulty to conduct the“refining”and“coarsening”operations since the curved information has to be maintained.Second，h-adaptivity could break the partition balancing，which would significantly affect the efficiency of parallel computing.In this paper，a robust and automatic h-adaptive method is developed for high-order DG methods on locally curved tetrahedral mesh，for which the curved geometries are maintained during the h-adaptivity.Furthermore，the reallocating and rebalancing of the computational loads on parallel clusters are conducted to maintain the parallel efficiency.Numerical results indicate that the introduced h-adaptive method is able to generate more reasonable mesh according to the structure of flow-fields.

Key words：h-adaptivity；discontinuous Galerkin（DG）method；curved mesh；tetrahedral mesh；compressible flows

0 Introduction

Many flow fields have multiscale structures which are difficult to numerically capture since it is impossible to give a proper mesh before the computation.Adaptive meshing techniques are introduced to increase the solution accuracy，as well as to reduce the computational time and computer memory.There are three typical adaptive meshing techniques：r-adaptivity，h-adaptivity and p-adaptivity.The r-adaptivity makes the mesh distribution match the flow structures by changing the positions of the grid points［1-6］.The h-adaptivity refines the region including small structures by dividing the original elements into smaller elements，and coarsens the region where the solution is smooth enough by agglomerating smaller elements into larger elements［7-11］.The p-adaptiviy methods are only suited to high-order methods and they employ different orders in different regions instead of changing the mesh［12-17］.Besides，the three adaptive meshing techniques can be used simultaneously.For example，the h-adaptivity changes grids and the p-adaptivity uses different orders for different elements.

Due to the limit of the number of elements and the highest order，the r-adaptivity and the p-adaptivity methods may not be able to provide sufficiently accurate solutions for complex flows.Therefore，the h-adaptivity is adopted in this paper.Babuska introduced h-adaptivity to finite element method to improve the accuracy in 1975［7］.Lentini and Pereyra combined h-adaptivity and finite difference method in the same year［18］.Then，h-adaptivity is widely used for the finite volume method［19-21］.In recent years，some researchers have applied h-adaptivity and discontinuous Galerkin（DG）methods to solve 2D flows［22-29］.Wu and Yu［30］employed h-adaptivity and DG method to solve the steady-state 3D Euler equations.

Compared with the traditional finite volume method，there are extra difficulties to implement the h-adaptivity for high-order DG methods on tetrahedral mesh.The curved geometry has to be maintained during the refining and coarsening operations.Furthermore，re-balancing operation is necessary for ensuring the efficiency of the parallel computation.In this paper，a 3D h-adaptive high-order DG method with a re-balancing parallel strategy is introduced to solve compressible flow problems on locally curved tetrahedral mesh.The details of the computation scheme are given，which includes the numerical method，h-adaptivity indicator，3D h-adaptivity and the parallel computing.Then，the numerical results are displayed.

1 Computation Scheme

1.1 Numerical method

1.1.1 Governing equations

The Navier-Stokes equations in conversation form are given by

whereUare the conversation variables.Fc(U)andFv(U，?U)are the inviscid and viscous flux terms，respectively.

whereρis the density，pthe pressure，Ethe total internal energy per unit mass，hthe total enthalpy per unit mass，andτthe stress tensor.u，v，andware the velocity components.

1.1.2 DG discretization

Employing the well-known BR2 scheme［31］，Eq.（1）can be discretized as

whereφiare the basis functions，，

where

The final discrete system can be written as

The above system can be solved using the following implicit relaxation method.u(0)=u(n)

1.2 h-adaptivity indicator

In order to capture the shocks，an artificial viscosity approach［29］is employed for transonic and supersonic cases，for which an artificial viscosity term is introduced to the original governing equations

whereεis defined as the artificial viscosity coefficient in elemente

wherese=log10Se，ε0=De/p，s0=log10(1/p4)，κis an constant andDethe size of elemente.Seis given as

Sinceεtends to 0 in smooth region and becomes larger in the region with great gradient such as the region near shocks，it can be used as an hadaptivity indicator.For the viscous flows without shocks，the magnitude of the vorticity can be employed as the h-adaptivity indicator.

1.3 h-adaptivity

Tetrahedron elements are employed in the entire computational domain here，and the h-adaptivity operations consist of refining and coarsening.Fig.1（a）displays the refining procedure，where the original“father”element on the left is divided into eight“child”elements on the right.To ensure the accuracy of the solid boundaries，the elements near the solid boundaries are usually curved.Fig.1（b）demonstrates the way to divide a curved element into smaller ones，where“child”elements are also curved to accurately match the geometry of“father”element.The coarsening procedure is exactly the inverse operation of those shown in Fig.1.Note that the way to describe curved elements is not unique.For example，we can use both high-order polynomials and governing points located on the edges and faces of an element to define the element geometry.In this paper，high-order polynomials are used，from which all required curved information can be computed，thereby maintaining the curved geometries during the refining and coarsening operations.

Fig.1 Element refining

A generation label is created for each element over the computational domain.All of the original elements are labelled 0 and the“child”elements generated after refiningntimes will be labelledn.To avoid the significant difference of the neighboring elements in size，the difference in the generation label of the neighboring elements are limited to be 1.For example，once the element“e”in Fig.2 is refined，the element“f”will also be refined to keep the local mesh size relatively smooth.

Fig.2 Refining the neighboring elements

The elements over the computational domain are numbered as shown in Fig.3.Once the elementEin Fig.3（a）is refined into eight“child”elements，one inherits the numberEand the others are numbered fromN+1 toN+7.Fig.3（b）indicates the inverse operation，where the elementsEandf1—f7are agglomerated back to their“father”elementE.This dynamic ordering rule maintains the chain structure of the original mesh，for which the computational module does not need to be modified during the h-adaptivity.

Fig.3 Data structure of elements during refining and coarsening

1.4 Parallel computing

It is well known that DG methods are well suited for parallel computing.The computation in this paper starts from partitioning the initial mesh intoNpartitions which have nearly equal number of elements with each other.However，after conducting h-adaptivity，the partition load balance would be broken，which means that the number of elements of each partition may be quite different.The unbalanced partition load would result in poor efficiency of the parallel computing.

This paper re-balances the mesh partitions after h-adaptivity.First，the elements to be refined or coarsened are marked.Second，all the processes send the marked h-adaptivity information to the first process.Third，the first process conducts the refining and coarsening over the entire domain.Finally，the first process redivides the entire computational domain intoNpartitions and distributes them to the corresponding processors.Then the computation continues on the new mesh.

2 Numerical Results

2.1 Inviscid flow around ONERA M6 wing

To validate the h-adaptivity method for the shock capturing，the typical inviscid transonic flow around the ONERA M6 wing is numerically simulated，for whichMa∞=0.84，α=3.06.The 3D Euler equations are solved using the DG method（p=2），where the artificial viscosity term is added into the Euler equations to capture the shocks，as shown in Eq.（3）.

Fig.4 presents the results computed on the initial mesh.Fig.4（a）shows the mesh distribution and there are 31 044 tetrahedral elements in total over the entire domain.The contours of the pressure coefficient are displayed in Fig.4（b）.The distribution ofεis given in Fig.4（c），which matches the shocks well.In the region whereεis relatively large，the elements are going to be refined to improve the resolution of the shocks.

Fig.4 Numerical results obtained on the initial mesh of ONERA M6 wing

Using the h-adaptivity method developed previously，a more accurate numerical result is obtained after only one time h-adaptivity.Figs.5（a），（b）display the final mesh distribution on the surface of the wing and the refined elements in the space，respectively.Fig.5（c）presents the new contours of the pressure distribution，where the shocks are much more accurately captured than that in Fig.4（b）.The distributions of the pressure coefficients on different spanwise profiles are given in Fig.6，where the shocks obtained after the h-adaptivity are sharper than those captured on the initial mesh and match better with the experimental results.

Fig.5 Numerical results obtained on the final mesh of ONERA M6 wing

Fig.6 Distributions of pressure coefficients at different spanwise profiles of ONERA M6 wing

2.2 Inviscid flow around pitching NACA0012 wing

To verify the applicability of the introduced hadaptivity method to unsteady flows，the inviscid flow around a periodically pitching NACA0012 wing is numerically simulated，for which the infinite Mach number is 0.755.The wing pitches aroundZ=0.25 periodically and the attack angle varies asα(τ)=α0+αmsin(2Kτ)，whereα0=0.016°，αm=2.51°，K=0.081 4.The DG method of 3rd-order accuracy（p=2）is employed and the explicit Runge-Kutta time-stepping is used.

Fig.7 Initial mesh of oscillating NACA0012 wing

Fig.7 displays the global view and the local view of the initial mesh，for which there are 20 312 elements in total.Fig.8 gives the contours of Mach number obtained on the initial mesh at different time steps during one periodic cycle.Note that significant numerical oscillation can be observed near the shocks，which means the numerical solution is not accurate enough.The mesh distributions with hadaptivity during the pitching in one period cycle are presented in Fig.9，where the maximum time of refinement is limited as two.In Figs.9（a），（b），the elements on the upper surface of the wing near the shocks are refined.Figs.9（c），（d）display the mesh distributions when the shocks appear on the lower surface，for which the elements near the shocks are refined and the previously refined regions in Figs.9（a），（b）are coarsened.Fig.10 displays the Mach contours at different time steps.During the h-adaptivity，the total number of the elements over the entire computational domain remains less than 42 000 since the refining and the coarsening operations are conducted simultaneously.

Fig.9 Refined mesh of oscillating NACA0012 wing

Fig.11 displays the distributions of the pressure coefficient on theZ=0 profile at different attack angles，where the results computed on the adaptive mesh match better with the experiments than those obtained on the initial mesh，and the numerical oscillations near the shocks are significantly reduced by the h-adaptivity.The comparison of the hysteresis loops is given in Fig.12，where the one computed on adaptive mesh is closer to the experimental one than the result obtained on the initial mesh.

Fig.10 Mach contours computed on the final mesh of oscillating NACA0012 wing

Fig.11 Comparisons of pressure coefficient profiles at different time steps

Fig.12 Comparison of hysteresis loops

2.3 Viscous flow around sphere

For the viscous flows without shocks，we use the magnitude of vorticity as an h-adaptivity indicator.Here we use the 3rd-order DG（p=2）to solve the viscous flow around sphere，for whichMa∞=0.5 andRe=118.The governing equations are the N-S equations without the artificial viscosity term in Eq.（3）.

The initial mesh is given in Fig.13 which contains 2 412 elements.The steady-state flow field is solved using the widely used Newton method.The pressure contours and the streamlines obtained on the initial mesh are displayed in Fig.14（a），where the pressure distribution is not smooth enough.The distribution of the magnitude of the vorticity is given in Fig.14（b），where the magnitude of the vorticity in the region near the solid surface and part of the wake zone is relatively large.

The initial mesh is refined one time according to the distribution of the vorticity.Figs.15（a），（b）display the final mesh in the symmetric profile and in the spaces，respectively.The Mach contours in Fig.16（a）are smooth and the accuracy of the pressure distribution is improved compared to that in Fig.14（a）.

Fig.13 Initial mesh for sphere

Fig.14 Numerical results on the initial mesh

Fig.15 The final mesh

Fig.16 Numerical results on the final mesh

3 Conclusions

（1）An h-adaptive DG method on locally curved tetrahedral mesh is developed，which contains both refining and coarsening operations.The curved geometry is well maintained during both refining and coarsening，which ensures the applicability of the obtained mesh to high-order DG methods.

（2）To ensure the parallel efficiency，a re-balance strategy for parallel computing is conducted after the h-adaptivity，which can provide a nearly equivalent computational load on different processors.However，one possible disadvantage of the rebalance method is that it may be time-consuming when the number of elements is large since only one processor is used to conduct the h-adaptivity and the re-balancing.

（3）Both steady-state and unsteady-state cases are numerically tested.Numerical results indicate that the accuracy can be significantly improved with the introduced h-adaptive method at a relatively low expense since only local regions are refined according to the given h-adaptivity indicators.Furthermore，the choice of the h-adaptivity indicators is flexible.

Typical test cases are numerically computed to validate the introduced h-adaptive DG method.In the future，we should consider further improving the efficiency of the h-adaptivity and the re-balancing operation for more complex computational fluid dynamics（CFD）problems，such as the large eddy simulation（LES）for complex geometries.

AcknowledgementsThis work was supported by the funding of the Key Laboratory of Aerodynamic Noise Control（No.ANCL20190103），the State Key Laboratory of Aerodynamics（No.SKLA20180102），the Aeronautical Science Foundation of China（Nos.2018ZA52002，2019ZA052011），the National Natural Science Foundation of China（Nos.61672281，61732006）.

AuthorsDr.AN Wei received the B.S.degree in aircraft design and engineering from Nanjing University of Aeronautics and Astronautics（NUAA）in 2017.He is currently an doctoral candidate majoring in fluid mechanics of NUAA.His research domain is the discontinuous Galekin method.

Prof.LYU Hongqiang is a Ph.D.tutor in College of Aerospace Engineering of NUAA.His research interests are computational fluid dynamics，artificial intelligence application，computational electromagnetism，and aerodynamic noise.He has published more than 30 academic papers on wellknown domestic and foreign journals.

Author contributionsProf.LYU Hongqiang designed the whole study.Dr.AN Wei performed the simulation，and wrote the manuscript.Dr.AN Wei and Mr.HUANG Zenghui performed the result investigation and figure generation.

Competing interestsThe authors declare no competing interests.