Research on high-order flux 1 reconstruction method with adaptive 2 mesh for shock capturing

6 High order schemes have been developed for a quite long time, and many famous schemes arise like 7 Discontinuous Galerkin (DG), Spectral Difference (SD) and WENO schemes, etc. The Flux Reconstruction 8 (FR) scheme proposed by Huynh has attracted the attention of researchers for its simplicity and 9 efficiency. It’s written in differential form and bridges the DG and SD scheme s, which can be constructed 10 with a proper choice of parameter. In this paper, realize FR scheme based on the framework of the open 11 source Adaptive Mesh Refinement (AMR) library p4est. To capture shock sharply, the performance of 12 Localized Laplacian Artificial Viscosity (LLAV) and In-cell Piecewise Integrated Solution methods are 13 compared. As an important way of reducing computational cost, AMR technique integrated in p4est is 14 also combined with FR. The performance of the developed code is tested in both one dimensional and 15 two dimensional and get some quite attracting results.


21
There has been a great development in computational fluid dynamics in recent years. The Finite 22 Volume Method (FVM) has obtained a great success, and was applied in both commercial 23 application and academic research widely. However with the increasing requirement for efficiency 24 and precision, FVM shows a poor performance in the simulation of turbulence, capturing of shock, 25 combustion, or the vortex propagation. This motivates researchers to explore new methods to 26 overcome this difficulty. 27 There are basically two ways to minimize error. One is the development of high order schemes. 28 High order methods show a better performance than FVM by constructing high order distribution 29 function with local information. The well-known high order scheme Discontinuous Galerkin 30 method (DG) was first proposed by Reed [1] and was extended to CFD by Shu [2] later. It's easy to 31 implement and is stable, which is important in CFD research. It has attracted researches to modify 32 DG scheme. Later based on DG, many mutants like the SD [3], SV [4] [5] schemes of DG were 33 developed. 34 In 2007 Huynh [6] proposed Flux Reconstruction scheme. It is easy to implement and with a 35 proper choice of solution points and correction function, both DG scheme and SD scheme can be 36 constructed directly, and is simpler than both of them in some way. Because FR scheme only 37 involves local cells and the values on the boundary of adjacent cells, which makes the scheme 38 (1) 96 First consider partitioning the computational domain [a, b] into N non-overlapping cells, define 97 i-th cell with   To consider the interaction of neighboring cells, common value on boundary needs to be 113 considered. Use L and R to represent the value on left and right boundary, 114 Here the flux function is not continuous. Riemann flux id used as the common flux, in this paper 116 the well-known Roe scheme is chosen. 117 With correction function  () g , the continuous flux function can be written as 119 120 here the last two terms are the correction of left and right boundary. ( 1) 1, There are many types of correction function can be applied which can be found in Huynh's 124 research. In particular, Radau function makes the scheme of order 2K-1. In this paper, Radau 125 function is applied to get high accuracy. 126 To solve the conservation law, flux function needs to be differentiated. 127 By now the derivative of flux in computational domain is complete, with the transform equation 129 of x and  , we can get the derivative in physical domain (10) 131 And the conservation law becomes 132 Huynh pointed that in linear cases the choice of correction function has a great influence on the 134 stability and precision, the order of scheme is irrelevant with the type of solution points. In 135 Jameson's research [34], the type of solution point is important for it can affect the energy property 136 and aliasing error when conservation laws are non-linear. Gauss points can eliminate aliasing error. 137 In this paper Gauss points are chosen as solution points. Finally, explicit or implicit time marching 138 methods can be applied to solve the differential equations. 139 140

Two dimensional FR scheme 141
For quadrilateral mesh, it's easy to extend FR scheme to two dimensional without much effort. In 142 this paper the grid lines of applied mesh are all horizontal or vertical. And the calculation can be 143 conducted directly in one dimensional manner. 144 145 3

146
To overcome the problem of high order schemes in shock cases, shock capturing techniques need 147 to be applied. In this section the artificial viscosity and in-cell piecewise integrated solution 148 methods are introduced and compared. To improve computing efficiency the adaptive mesh 149 refinement technique is also presented. 150 151 5

Artificial Viscosity 152
Basically Artificial Viscosity technique is constructed based on local value, which is easy to 153 implement in high order schemes. But how to define the amount of viscosity is still an open 154 question for there are adjustable parameters that are problem dependent. 155 In this paper, the LLAV (Localized Laplacian Artificial viscosity) method proposed by Persson 156 [24] is applied, which has the property of limiting shock within a cell. The LLAV method adds a 157 Laplacian diffusivity term in control equations, which makes the equations become 158 On the right hand side is the Laplacian artificial viscosity term,  is the viscosity parameter, which 160 is relevant to the property of local flow. In Persson's research,  / hp , h is the cell size, P is the 161 order of solution polynomial. And  is expressed as follow 162 here e s is the smooth indicator which will be introduced shortly. In equation (13)  Adaptive mesh refinement method 203 In this paper, the AMR library p4est was applied to realize AMR efficiently. P4est [35] is an open 204 source AMR library provides many useful functions. Based on octree data structure, p4est has a 205 complete system of parallel generating mesh, refining mesh, coarsening mesh, load balancing, cell 206 searching, etc. It has a high parallel efficiency and is scalable. And the code is developed under the 207 framework of p4est. 208 209

Grid adaption criteria 210
In AMR, besides the complicate data structure and the algorithms, another important question is 211 the decision whether to refine or coarsen the grid during simulation. In fluid dynamics the common 212 choice is the gradient of physical variables, which regards that regions where flow changes rapidly 213 7 have important influence to the precision. And many physical variables could be chosen as 214 indicator like density, temperature, pressure, entropy, etc. 215 Basically, there are some differences in the choice of indicator for different flows. In flows that 216 contain discontinuities, the shock detector could be used as the indicator. When  − 0 e ss , which 217 means solution is smooth, the grid is marked "coarsen". Otherwise it indicates that discontinuities 218 exist, and the grid is marked "refine". In flows without discontinuity, the vortex flow is more 219 important. To accomplish a better vortex simulation, the following criteria presented in [36] is used 220

Nonconforming data projection 226
As mentioned above, FR scheme needs to calculate the common flux on each cell boundary at each 227 time step. When computational grid is uniform and stationary, the edges of neighboring cells have 228 the same length and flux points are collocated. And common numerical flux can be calculated 229 directly. When AMR is applied, the common edges may be not equal length, called non-conforming 230 edge. In this case, the flux points on each edge are not collocated, which is shown in Fig. 2

SOD shock tube 245
The SOD shock tube problem is a classic case of shock capturing, with the initial condition 246 In this paper, the grid size is dx=1/100. Fig. 3 shows the results with artificial viscosity and in-248 cell piecewise constant approximation methods. And they are compared with different order at 249 t=0.2s. From the results we can see that both the methods show satisfying results of shock 250 capturing, which solve shock precisely with little oscillation. On the other hand, the numerical 251 solution is more precise with the solution order increased, which is also the same as expected.

One dimensional Shu-Osher problem 257
The well-known Shu-Osher problem is the interaction of a density sine wave with a shock with 258 speed of Mach 3 moving to right. It contains small flow structures and shock. The initial condition 259 of Shu-Osher problem is 260 The computational domain is [-5, 5], the grid size is 1/20, the order of accuracy is 3, Fig. 4 shows 262 the density results at t=1.8s. From the results, we can see that both the methods are capable of 263 capturing small flow structure and there is no obvious oscillation near shock.

Shock-Vortex interaction 270
The shock-vortex interaction flow of the stationary shock of Mach 1.1 interacts with an isentropic 271 vortex. Here third order scheme is applied. The case has a computational domain of      0,2 0,1 , 272 the grid size is 1/60, the stationary shock is placed at x=0.5. Pre-shock initial condition is 273 ( 1) 8 Here  = 5 , the computational domain is     0 10,0 10 xy , the vortex center is located at (5, 5), 296 the boundary condition is periodic condition, and the grid is  20 20 . In this paper, the FR scheme 297 with AMR is tested with this case. Here give the results at t=10s. Fig. 6 shows the computational 298 grid at different time. Fig. 7 shows the density distribution with different order and with AMR. 299 From these results we can see that the grid could automatically refine and coarsen with the vortex 300 indicator without affecting precision, which greatly improves the simulation. 301

Oblique shock reflection 307
In this paper, mainly third order scheme is applied for 2D cases. First we tested the oblique 308 shock reflection problem. The initial condition is the oblique shock of Mach 3 reflected on an 309 inviscid wall with an inflow angle of 29 degree respect to x axis. The results are shown in Fig.   310 9, it's easy to find out that with AMR the shock is much sharper than the original one.

Double Mach reflection 319
The double Mach reflection problem is to simulate the flow of an oblique shock of Mach 10 is used 320 to test the scheme's ability of capturing strong shock. The double Mach reflection problem has a 321 computational domain of  [0,4] [0,1] , the lower boundary at  0 1/ 6 x and the left boundary are 322 set the post-shock initial condition, lower boundary at  1/ 6 4 x is set the reflect wall condition., 323 right boundary is set free outflow condition, upper boundary is set the same as the exact movement 324 of shock. The pre-shock and post-shock initial condition is 325 In this case, the initial grid size is 1/60, with third order simulation, and capture shock with 327 artificial viscosity. Fig. 10 shows the density contour and grid distribution at t=0.2s. Fig. 11 shows 328 the zoom in result of Double Mach Reflection. From these results, we can see that with artificial 329 viscosity, the shock could be captured precisely, and with the grid refinement in shock region, the 330 shock is sharp and flow structure is precise.

338
In this paper, FR scheme is realized based on the framework of AMR library p4est in two 339 dimensional. Compare the performance of shock capturing ability of the Artificial Viscosity and In-340 cell Piecewise Integrated solution methods. Both of them show pleasant results, the Artificial 341 Viscosity can have more precise results but is more complex and needs to adjust parameters which 342 are problem dependent. 343 With the help of p4est, the AMR is added to the code. From the simulation of both subsonic and 344 supersonic cases we can see that the AMR is successfully applied to FR scheme. 345 346

5.1
Further work 347 1. Only cases in Euler equations are simulated, the Navier-Stokes equations will be added to the 348 code; 349 2. Three dimensional simulation will be added to the system; 350 3. The curved boundary needs to be applied in real geometry; 351 4. The time marching method is Runge-Kutta method, which has a strict limitation to time step, 352 later the implicit time marching methods would be added. 353 354