Adaptive Mesh Refinement with High-Order Scheme ...

AIAA 2014-0077 AIAA SciTech Forum 13-17 January 2014, National Harbor, Maryland 52nd Aerospace Sciences Meeting

Adaptive Mesh Refinement with High-Order Scheme for an Unstructured Pressure-Based Solver H. Q. Yang1, Z. J. Chen2, and Andrzej Przekwas3 CFD Research Corp. Huntsville, AL 35805 and Jonathan Dudley4 USAF Research Laboratory (AFRL/RWWC), Eglin AFB, Florida 32542-6810

Downloaded by Hong Yang on January 19, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-0077

Abstract In this study, high-order schemes of up to 4th-order spatial and temporal accuracy are combined with dynamic adaptive mesh refinement to study two types of fluid flow: incompressible vortical flow and highly compressible shock flow. The adaptive mesh refinement strategy starts with an initial body-fitted unstructured grid and then identifies the logical elements or zones where additional resolution is needed at run time. These zones are then "block refined" into grid patches that are comprised of even finer zones. The unstructured high-order adaptive method is demonstrated for a channel flow with shock reflections and transitional incompressible flow over a SD7003 airfoil. It is shown that in order to resolve turbulence to the Kolmogorov scale, grid refinement is crucial in regions of finite Reynolds stress (more than just in regions of high pressure or velocity gradient), as this is where the transition occurs. The adaptive scheme is evaluated using different metrics related to turbulence features including vorticity magnitude, mean Reynolds-stress, and Q-criterion. The study shows that the mean Reynolds stress-based grid adaptation is the best approach for transitional flow as it takes the turbulent flow physics into consideration and focuses the grid refinement on the most critical regions of the flow field.

I. Introduction The accuracy of a CFD computation depends on two kinds of error sources: modeling errors and numerical or discretization errors. It has been shown that, for a CFD simulation of a standard airplane at cruise conditions, the numerical accuracy mainly depends on discretization errors [1]. It is well known that, in the asymptotic range, discretization errors are proportional to the local size of the mesh as O(hp), where h represents the local size of the mesh elements and p is the formal discretization order of the numerical method. The reduction of the numerical error may be accomplished by a twofold strategy [2]: 1) to increase the order of the discretization scheme and/or 2) to increase the number of mesh elements. These two techniques are known as p and h adaptation, respectively. The p or even hp adaptation [3] techniques are usually considered within the finite-element community, where the order p can be prescribed in a flexible way. In the finite-volume community, however, h refinement is usually preferred because of the very nature of this discretization approach. In an unsteady fluid dynamics simulation, a continuous domain of interest is discretized into a grid of many individual elements. This grid may be static and established at the beginning of the computation, or it may be dynamic, tracking the features of the result as the computation progresses. To resolve the flow features that are much smaller than the overall scale of the problem, and that evolve in time, one can either include many more static grids to cover the region of interest, or apply a dynamically adapting scheme. A dynamic gridding scheme has the advantages of increased computational and storage savings over a static grid approach and complete control of grid resolution based on the flow physics.

1

Chief Scientist, CFD Research Corp., 701 McMillian Way, Huntsville, AL 35806, AIAA Senior Member. Technical Fellow, CFD Research Corp., 701 McMillian Way, Huntsville, AL 35806, AIAA Member 3 Senior VP and CTO, CFD Research Corp, 701 McMillian Way, Huntsville, AL 35805, and AIAAA Associate Fellow 4 Research Engineer, AFRL/RWWC, 101W Eglin Blvd., AIAA Senior Member. 2

1 American Institute of Aeronautics and Astronautics 092407

Copyright © 2014 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved.

Many dynamic mesh approaches rely on a combination of mesh modification operators such as node insertion, element reconnection, and node movement. For example, some existing Adaptive Mesh Refinement (AMR) techniques employ a Cartesian grid where the mesh is square in 2D and cu cubic bic in 3D [4,5], as shown in Figure 1a. 1 This type of grid system has limitationss in resolving curved boundaries, where a “stair case” often has to be used to model a surface that is actually curved.


In this study, we combine high-order order schemes up to 4th order in space and 4th order in time with dynamic adaptive mesh refinement to study vortical flow at the incompressible limit and ccompressible ompressible flow with shock waves. Our present adaptive mesh refinement uses a body body-fitted, unstructured grid id system, as shown in Figure 1b. 1 During the simulation run,, the adaptive mesh refinement scheme will identify the logical elements, or zones, where additional resolution is needed. Those zones needing more detail are then "block refined" into grid patches that have finer zones. In other words, to show more detail, one original zone might become four zones for quad and three zones for a triangular grid in 2D, as illustrated in Figure 11b.

Figure 1. Conventional Cartesian grid and the present body body-fitted, fitted, unstructured adaptive grid. In the following, the standard pressured pressured-based based CFD solution methodology will be first presented. The implementation of the proposed high-order order spatial schemes using the successive differentiation approach will then be derived. Extension of the 2nd order backward temporal scheme to 3rd and 4th order accuracy will also be described. Several benchmark problems will be used to verify the achievement of the designed order-of-accuracy order both for spatial and temporal schemes schemes. The problems of shock wave in a channel with multiple reflections and transitional flow over the SD7003 airfoil will be used for validation and demonstration of the adaptive mesh refinement.

II. Numerical Method Computational Fluid Dynamics technology for solving Navier Navier-Stokes Stokes equations can be classified into two categories: density-based methods and pressure pressure-based methods. The density-based method,, originally developed for external compressible flow problems, retain retains density as a dependent variable in the solution of continuity equation. The density-based method is formulated on a mathematically rigorous evaluation of the inviscid flux at control volume cell faces. By incorporating the nonlinear wave properties into the flux functions in the form of Riemann problems and characteristic equations equations, the density-based solvers are both robust and accurate for high-speed compressible flows. Linearization of the flux function results in a coupled system of equations that are linked via a flux Jacobian matrix. At each control volume volume, a set of momentum, continuity uity and energy equations are solved simultaneously. Some examples of commercial and produc production level density-based based codes are: are OVERFLOW [6], USM3D [7], Loci/CHEM [8],, FUN3D [9], CFD-FASTRAN [10], and CFD++ [11]. Although the density-based density method has been successfully applied to high high-speed compressible flows, in the incompressible limit numerical 2 American Institute of Aeronautics and Astronautics 092407

difficulties arise in the calculation of pressure from the equation of state. These include 1) round off error due to using density as a primary variable, and 2) a time step (or Courant-Freiderichs-Lewy number) constraint due to the infinite acoustic speed. Various techniques have been proposed to extend the applicability of these schemes to incompressible flows including preconditioning, dual time stepping, or through the addition of a 'pseudo' (artificial) compressibility term into the continuity inuity equation equation. The convergence rate is, however, highly dependent on the choice of problem dependent parameter values values.


On the other hand, in the pressure--based method, the continuity equation is formulated such that pressure is a primary variable and density is found from the equation of state. There ar aree no limitations, due to the type of variables, to using a pressure-based based method for compressible flow calculations. The pressure pressure-based methods have been used to develop the calculation procedure procedures that are valid for or the entire spectrum of Mach numbers. As the pressure-based solver can formulate the Navier-Stokes equations in a sequential solution procedure, it is efficient in terms of memory requirement and simulation time. Pressure-based solvers can also handle multiphase flow problems in industrial and aerospace rospace applications. Some examples of commercial and production-level production pressurebased solvers are: CFD-ACE+ [12],, ANSYS ANSYS-Fluent [13], Star-CD [14], Loci/Stream [15],, and CoBi in this study. study Table 1 lists some feature comparisons between pressure pressure-based and density-based CFD solvers.

Table 1. Feature ture Comparison between Density Density-Based and Pressure-Based ased CFD Solvers

This paper represents our continuous ontinuous effort in developing a multi-physics, pressure-based based fluid dynamics solver for structured and unstructured grids [16 [16-18], and in developing a fully-coupled fluid-structure structure dynamics solver [19]. In the following, the governing equation ations will first be presented, and a description of the discretization procedure will follow. Then, details of the high-order order spatial and temporal schemes will be presented. A. Governing Equations solved are the conservation of momentum, mass, and energy: The governing equationss to be solve

∂ρu i ∂p ∂τ ij ∂ + ( ρu j u i ) = − + ∂t ∂x j ∂x i ∂x j

i = 1,2,3 .


(1)

∂ρ ∂ρu i + = 0. ∂t ∂x i

(2)

∂ρht ∂ ∂ ∂T ∂p ∂ ( ρu j ht ) = k (τ ij u j ) . + + + ∂t ∂x j ∂x j ∂x j ∂t ∂x j

(3)

The shear stress is written as:


 ∂u ∂u j  2 ∂u k − µ τ ij = µ  i + δ .  3 ∂x ij ∂ x ∂ x i  k  j

(4)

The momentum and energy equations can be cast into a standard transport equation for a general dependent variable f as:

∂ρf ∂ ( ρu j f ) ∂ ∂f + = Df +Sf . ∂t ∂x j ∂x j ∂x j

(5)

B. Discretization By integrating Eq. 5 over a control volume P as shown in Fig. 2 (a), one has:

∫

Vol

∂ρf ∂f + ∫ ρu j n j fdΓ = ∫ D f n j dΓ + ∫ S f dVol , Γ Γ Vol ∂t ∂x j

(6)

where Vol is control volume, Γ is face area of the control volume, Df is diffusivity, and Sf is a source term. Function f is 1 for the continuity equation, ui for the momentum equation and ht for energy equation, respectively.

a)

b)

Figure 2. General control volume in 2D. P is the main control volume and Bi is a neighbor control volume. Γ is the boundary of the control volume and n is the face normal vector. a) relation of P to its neighbor; b) relation of face e to the surrounding cell center and corner and centers. The discretization of the governing differential equations is obtained using a cell-centered finite volume approach, where cell-averaged flow variables are placed in the center of the control volume. The control volume discretization of the governing equations is expressed in flux conservation form and involves numerical evaluation of convection and diffusion fluxes at all cell boundaries and a transient term at the control volume center. For 4 American Institute of Aeronautics and Astronautics 092407

convenience, we will present the details of the flux evaluation at a selected face, denoted by a subscript e (see Fig. 2b), but the same evaluation applies to all faces of the control volume. The discretization will be presented for the convective and diffusive fluxes, and the transient term. C. Convective Flux The convective flux can be evaluated as:

Fe = ∫ ρu j n j fdΓ = ( ρVn A) e f e ,

(7)

Γ

where fe is the value at the cell face center. The example interpolations for the 1st order upwind and 2nd order central schemes are: 1st order upwind scheme:


(a)

fe = fE ; 2nd order central scheme:

(b)

f e = f P , ; if Vn ≥ 0 .

(8)

if Vn ≤ 0 .

(9)

f e = ε P f P + (1 − ε P ) f E .

(10)

Implementation of higher-order schemes for unstructured grids is very challenging. It involves reconstruction of nodal values, cell-center local values and gradients based on the cell averages. Based on that information, cell face fluxes are evaluated from the left and right states. For the conventional 2nd order schemes to avoid numerical oscillations in regions of steep gradients, the reconstruction step is expressed as [17]:

f f = f P + ψ∇f • (r f − rP ) ,

(11)

where ψ is a limiter function. Example limiter functions are the Barth limiter [20] and the Venkatakrishnan limiter [21]. The second order upwind scheme can be expressed as:

f e = f P + (∇f ) P • (re − rP )

if Vn ≥ 0 .

(12)

f e = f E + (∇f ) E • (re − rE )

if Vn < 0 .

(13)

D. High-Order Schemes The technical approach to obtain high-order accuracy has been described in detail in references [22,23]. In general, one can also write the face value of a variable in a compact form as:

f face = f center + ∆xi where:

∂f 1 ∂2 f 1 ∂3 f + ∆x i ∆x j + ∆x i ∆x j ∆x k ∂xi 2! ∂xi ∂x j 3! ∂xi ∂x j ∂x k

∆xi = xi , face − xi ,center

The Gauss-Green theorem will be applied to f,

∂f ∂2 f ∂3 f , , and ∂xi ∂xi ∂x j ∂xi ∂x j ∂x k

For example, we have:


(14)

∂f 1 = ∑ fΓ AΓ ni ∂xi Vol Once the field for

(15)

∂f ∂2 f is built, can be computed as: ∂xi ∂xi ∂x j  ∂f  ∂2 f 1   AΓ n j = ∑ ∂xi ∂x j Vol  ∂xi Γ

(16)


With the second derivative computed, the third derivative can also be computed as:

 ∂2 f ∂3 f 1 = ∑ ∂x i ∂x j ∂x k Vol Γ  ∂x i ∂x j

  AΓ n k  Γ

(17)

The key to extending the current high-order scheme to multiple dimensions is in recovering the cell-centered solution values from the cell-averages. In multiple dimensions, Eq. (14) is written as:

f face = f center + ( ∆xi − xi , center )

∂f 1 ∂2 f 1 ∂3 f (18) + (∆xi ∆x j − xij ,center ) + (∆xi ∆x j ∆xk − xijk ,center ) ∂xi 2 ∂xi ∂x j 6 ∂xi ∂x j ∂xk

The verifications of the high-order accuracy for many problems can be found in references [22] and [23]. E. Diffusive Flux The diffusion flux is divided into two parts [17]: one is main diffusion flux, Gem, and another is the cross diffusion flux Gec as following:

Ge = Ge + Ge = ∫ D f ∇f • n1 dΓ ≅ ( D f ∇f ) e • n1e Ae . m

c

(19)

Γ

By discretizing the gradient f at the face center e using a central difference scheme, the main and cross diffusion flux can be written as: m

Ge = ( D f ) e

c

Ge = ( D f ) e

Ae ( fE − fP ) , (rE − rP ) • n 1e

(20)

Ae (n1 • n 2 ) e ( f c 2 − f c1 ) , rE − rP

(21)

where fc1 and fc2 are values of f at corner nodes. In the above expression, the main diffusion flux will be included in the left-hand side matrix, and the cross diffusion flux will be treated explicitly as source term. F. Transient Term To account for transient effects, the governing equations must be discretized in time for the transient term. The popular temporal algorithms currently in use for aerodynamics are the first- and second-order accurate Backward Differentiation Formulae (BDF1 and BDF2). The BDF is a linear multistep method in that for a given function and


time, it approximates the derivative of that function using information from already computed times, thereby increasing the accuracy of the approximation. The BDF is expressed as follows: (a). BDF1: the 1st order backward Euler: n

1  ∂ρf  ( ρf ) n − ( ρf ) n −1 .   = ∆t  ∂t 

(22)

1 3 1   ∂ρf  ( ρf ) n − 2( ρf ) n −1 + ( ρf ) n − 2  .   =  ∆t  2 2   ∂t 

(23)

[

]

(b). BDF2: the 2nd order scheme:


n

To achieve high-order temporal accuracy, the 3rd order and the 4th order implicit backward differentiation formulae (BDF) are implemented in CoBi. The recipes for the 3rd and the 4th order BDF temporal schemes are: (c). BDF3: the 3rd order scheme:

n

1 11 3 1  ∂ρf   ( ρf ) n − 3( ρf ) n −1 + ( ρf ) n − 2 − ( ρf ) n −3    =  ∆t  6 2 3  ∂t  

(24)

(d). BDF4: the 4th order scheme:

1  25 4 1  ∂ρf   ( ρf ) n − 4( ρf ) n −1 + 3( ρf ) n −2 − ( ρf ) n −3 + ( ρf ) n− 4   n =  ∆t  12 3 4  ∂t  

(25)

G. General Linear Equation Form After summing up all cell face fluxes, the discretized equation can be expressed in the following algebraic form:

a P f P = ∑ a nb f nb + S fc ,

(26)

where anb are the link coefficients which involve contributions of convection and diffusion fluxes at neighbor points. As a result, the solution of a linear system of equations is necessary. The left-hand side matrix of the linear equation has a constant bandwidth with four non-zero entries for triangular cells, five non-zero entries for quadrilateral and tetrahedral cells, and nine non-zero entries for hexahedral cells. For the arbitrary polyhedra, the bandwidth includes all direct neighbors of the polyhedra. H. Pressure Correction Equation To formulate the pressure equation from the continuity equation, let us define face velocity as:

Ve = (1 − ε )VP + εVE + (1 − ε ) H PT (∇p) P + εH ET (∇p) E − H eT (∇p) e .

(27)

Define velocity as a function of pressure from momentum equation:

Ve ' = − H e (∇p' ) e . 7 American Institute of Aeronautics and Astronautics 092407

(28)

By substituting into the continuity equation:

∆ ( ρV ) P + ∑ (ρV f • A f ) * +( ρV f • A f )' = 0 . ∆t f

[

]

(29)

One can get the samee algebraic form of equation ((26).


I. Solution Procedure 1. Solve momentum equation first, and get estimate of velocity using existing pressure field; 2. Compute mass flux at interface; 3. Solve pressure correction equation using estimated veloci velocity field; 4. Update pressure field and velocity field, so that velocity will satisfy continuity equation; 5. The new velocity field may not satisfy continuity equation, so back to step one, and iterate until solution is converged.

III. Adaptive Mesh Refinement Adaptive daptive mesh refinement (AMR) has been used for over 20 years in density density-based based CFD solvers on Cartesian grids for accurate ate simulations of shock waves by Berger in 1984, 1989 [24,25].. However, AMR has not been demonstrated in pressure ressure based codes using mixed (pol (polyhedral) grids to resolve complex flow phenomena. A parallel version of solution-based based AMR in CoBi with 3D mesh refinement (near steep gradients) and coarsening (in smooth solution zones) has been developed at CFDRC [26] and it can be used to simulate complex fluid flow and shock waves such as nonlinear acoustic fields generated by blast shock waves or high speed je jets ts e.g. on the aircraft carrier and calculate pressure loads on un/protected human ears.

Figure 3. 2D cell ell and 3D surface refinement by mid mid-point insertion.


Approach Dynamic mesh esh adaptation methods can provide significant efficiency in CFD simulations, especially especia in simulations of vortical flow and shock capturing simulations. A mesh adaptation algorithm has as been implemented into the CoBi solver framework. When a cell is flagged for refinement, tthe basic strategy is to insert a mid-point mid in a plane surface and find the mid-point point along all the edges of that surface. By connecting the mid-point mid inside the surface and mid-point point on the edges, one can break one cell into multiple cells as illustrated in Figure 3. During this process, the conformal nature of a mixed element grid is maintained. For example in 2D, a triangular cell is divided into 3 cells, and a quadrilateral drilateral cell is divided into 4 cells. A polyhedral cell containing n sides is divided into n cells. In 3D the volume is divided first from the cell faces, as in 2D, and then the volume is divided. The second level of adaptation is to break new gener generated sub-cells cells into 4 cells again, as illustrated in Figure 4.


P2 P1

P1

Figure 4. Sub-cell break up to 4 new cells. The general algorithm can be designed into a tree tree-based structure shown in Figure 5. This tree data structure allows refinement and coarsening.

Figure 5. Tree based cell adaptation algorithm. The general flow chart for the adaptation is given in Figure 6. The selection of a reliable adaptation parameter is a key aspect to reduce the errors ors in the computation [27]. One classical approach for adaptation indicators in the literature is based on gradient and undivided difference-based feature detectors.. This approach, known as a featured-based indicator, can detect ect flow phenomena, when taking for instance the primitive variables into account. In the current implementation, the criterion to adapt a cell is based on the calculation of the gradient of selected flow quantities, such as the gradient of pressure, density, density and temperature, or shear stress, vorticity, and so on. Thee general strategy is defined as

λ=

∇ϕ ∇ϕ max

 > λo , refining = < λ1 , coarseing

,

(30)

where φ is a flow quantity, such as pressure, density, temperature or velocity magnitude. While ∇ ϕ can be vorticity, shear stress.

λo , λ1

are user input constants.



Q-Criterion Based Adaption In a typical vortical flow problem, the mesh adaptation is done based on the local vorticity. As shown in Figure 1b, although the grid refinement is adapted to the highly separated flow region, there are significant amount of grid cells being placed near the laminar boundary layer zone on both suction aand nd pressure sides (see Figure 1b). From a flow physics point of view, it is the transi ansition zone on the upper surface that contributes to the drag increase. The refinement is crucial in the non-zero zero Reynolds stress zone instead of the boundary layer. Actually, the Reynolds stress is nearly zero near the wall.

Figure 6. Mesh adaptation flowchart. In this study, different criteria eria for grid adaptation have been investigated. One of the choices is the Q-criterion, Q which is defined as a combination of vorticity and strain rate.

Q=

1 1 (Ω ij Ω ij − S ij S ij ) = (ω 2 − 2S ij S ij ) , 2 4

(31)

where Ωij and Sij are the anti-symmetric symmetric and symmetric part of the velocity gradient gradient, respectively, that is:

Ω ij =

∂u j 1 ∂u i ∂u j 1 ∂u ( − ) and S ij = ( i + ). 2 ∂x j ∂xi 2 ∂x j ∂xi

(32)

The Q-criterion thus us represents the balance between the rate of vorticity and the rate of strain. In the core of a vortex, Q > 0 since vorticity increases as the center of the vortex is approached. Thus regions of positive Q-criterion Q correspond to 3D vortical structures. A Although lthough this definition measures the difference between the rotation and strain rate magnitudes, this value is still dependent upon the characteristic length and velocity scales inherent to the specific problem. To yield a suitable non non-dimensional form, if normalized by ||S||2, a threshold function can be obtained that is of the form: 2

f threshold

1 Ω = ( 2 − 1) . 2 S

(33)

This threshold function will be evaluated in each cell which is marked for refinement if the th resulting function value is greater than a pre-specified specified value. Irrotational flow occurs when fthreshold →−1/2, −1/2, solid body rotation happens when fthreshold →∞,, and all positive values represent regions of swirl, i.e., regions where the rotational strengt strength is larger than the strain rate strength. Mean Reynolds-Stress Based Adaptat ation Another choice is to use the mean Reynolds Reynolds-stress by which only the turbulence region on the upper surface of the wing will be refined.. To build the refinement criterion, the time averaged values for the u velocity component, for example, can be determined as: 10 American Institute of Aeronautics and Astronautics 092407

n

∑ u(t )∆t i

u=

i

i =1

,

n

∑ ∆t

(34)

i

i =1

and its perturbation can be defined as:

u(ti )' = u(ti ) − u .

(35)

For the Reynolds stresses component of “uv”, then we have: n

∑ u (t )' v(t )' ∆t i


ρ u ' v' = ρ

i

i

i =1

(36)

n

∑ ∆t

i

i =1

Non-dimensionalization can be made using the local mean velocity.

IV. Results and Discussions This section will discuss the results using the above new high-order schemes on 1D scalar transport equation. Application of high order scheme for compressible flow will also be demonstrated. Finally adaptation examples using Q-criterion and mean Reynolds-stress will be given. A. One Dimensional Scalar Transport Problem To demonstrate the current high-order scheme, we first consider a one dimensional scalar transport equation of the form:

∂f ∂f + = 0. ∂t ∂x

(37)

The analytical solution to the above equation is:

f ( x, t ) = f 0[( x − t ),0] .

(38)

The solution depends on the initial condition. First a sinusoidal wave is used to evaluate the accuracy of the current high-order schemes. The computational domain has a size of 2 in x direction and a periodic boundary condition is applied at both ends. The initial condition is:

f 0 = sin(πx);

0 < x < 2.

(39)

After the wave has traveled half of the domain at t=1 using a very small CFL of 1.0e-3 with 2nd order BDF2, the solutions using different spatial schemes are compared to analytical solution. The L2 errors are plotted in Figure 7. One can see that: a) The conventional 2nd order scheme, either central or upwind, has the highest L2 error compared to the other schemes. b) The 3rd order upwind scheme exhibits its formal order-of-accuracy of 3rd order; c) The 3rd order central scheme has an observed order-of-accuracy of 4th order; d) The 4th order upwind scheme exhibits 4th order accuracy as expected; and e) The 4th order central scheme exhibits 5th order accuracy.



To ensure a 1.0e-3 L2 error, one only requires 10 points when using the 4th order central scheme, scheme while one requires 40 points when using the 3rd order upwind scheme. This is in contrast to the 80 points that t would be required when using the 2nd order schemes schemes. Thus the savings in required cell count for the 4th order central scheme over the 2nd order scheme is a factor of 8 for this 1D problem.

Figure 7. Solution error as a function of number of cells for the present high order schemes B. Verification of High-Order Order Temporal Schemes To demonstrate the current high-order order temporal scheme, we again take the previous one dimensional scalar transport equation (37) with initial nitial condition of equation ((39). The grid resolution is kept at 80 for a domain from x=0 to x=2. Periodic boundary conditions are applied at the two ends. Figure 8 shows the time-step step size dependence of different spatial schemes when using the traditional BDF2 temporal scheme for solution of the above scalar transport equation. It is obs observed erved that for the 2nd order spatial schemes, when the CFL number is less than 1.0, the solution is independent of time time-step step size. In other words, with the compatible order-of-accuracy accuracy for temporal and spatial sche schemes, mes, CFL=1.0 is a critical value. value When the CFL is less than one, the spatial error is the dominant truncation term in the solution; when the CFL is larger than one, the temporal accuracy dominates. One can also see that the log log-to-log slope of L2 error with CFL number is 2, which verifies the 2nd order accuracy of the BDF2 temporal scheme. For the 3rd order spatial schemes, a much smaller CFL number is required to ensure time-step step size independence. For the 3rd order upwind scheme, the critical CFL rd number is about 0.3, and for the 3 central scheme, which has 4th order accuracy, the CFL number must be less than 0.04 before the spatial accuracy becomes important. In other words, using the existing 2nd order BDF2 temporal scheme, one has to reduce the time-step step size to a value much smaller than 1 in order to take full advantage of a 3rd order spatial scheme. For the pressure--based based solver, as almost all equations are solved implicitly, reducing the time step to accommodate high-order order spatial accuracy implies an increase in computational time and a decrease in computational efficiency. The situation becomes even worse for the 4th order scheme. For example, for the 4th order upwind scheme, the critical CFL number is about CFL=0.1, and for the 4th order central scheme, which is 5th order accurate, the critical CFL number is around CFL=0.01. Even at this low CFL number, one can still see the sensitivity of the time-step size. Before we discuss the results from the BDF3 sc scheme, heme, it should be first noted that the multistep method using the Backward Differentiation Formulae is not self self-starting for high-orders. orders. At each time step, the solutions from several previous time steps are required, but these solutions do not exist at the first time step. Choices for starting these methods are to use lower order accuracy for the first few time steps, or to start the schemes with a multistage method, which is self-starting. starting. We have derived a variable time time-step sized and self-starting BDF DF formula. With this algorithm, for the first few time steps the low order scheme is used but with user specified smaller time steps. Our test showed that it can significantly reduce the effect of the initial low low-order order scheme and hence achieve the targeted higher order-of-accuracy. 12 American Institute of Aeronautics and Astronautics 092407

1.E-01

CFL Number Effect of 2nd Order Temporal Scheme

1.E-02 L2 Error

1.E-03 1.E-04

2nd_up 2nd_cn 3rd_up 3rd_cn 4th_up

1.E-05 1.E-06


1.E-07 1.E-03

1.E 1.E-02

1.E-01 CFL Number

1.E+00

1.E+01

Figure 8.. Time step dependence study for the 2nd order temporal scheme. Figure 9 shows the difference in error reduction for the 4th order spatial schemes using the BDF2, BDF3 and BDF4 temporal schemes. One can see the different error convergence rates for the BDF2, BDF3 and BDF4 schemes, verifying ifying the correct implementation of the 3rd and 4th order temporal schemes. Also, using the same CFL number (time step size) both the BDF3 and BDF4 scheme exhibit lower error than the BDF2 scheme. It also shows that the 4th order temporal scheme will al allow effective utilization of the 4th order central spatial accuracy with a reasonable CFL number of 1.0.

Figure 9.. Convergence rate of BDF2 and BDF3 for the 4th order spatial schemes

C. 2D Supersonic Flow Over a Bump in a Channel Consider ider a 2D domain as shown in Fi Figure 10, which corresponds to a rectangular channel of height L and length 3L. L. The channel has a circular arc bump (in the middle) of length L and maximum imum thickness equal to 0.1 L along the bottom wall. A grid comprised entirely of quadrilateral cells is utilized. The inlet Mach number is M=1.64, and upper and lower walls are specified as symmetry planes.



Figure 110. Supersonic flow over a bump in a channel. Given in Figure 11 are pressure distributions using different spatial schemes for a steady state solution. It is generally perceived that the pressure-based based method is only applicable for low speed flows, and it lacks good shock capturing capability. This example prove proves the capabilities of the current pressure-based based method in the application of a wide range of Mach numbers. Due to the existence of shock waves in the flow domain, the central schemes give raise to high oscillations in the solution. On the other hand, upwind schemes are stable. e. As one can see from Figure 16, the first-order order upwind method is very diffusive, and the shock waves have been smeared out. The 2nd order upwind scheme shows some oscillation ne near the shock wave, while the 3rd order er upwind scheme has significantly improved the solution accuracy. Using the 4th order scheme the result is very sharp and the shocks are well preserved.

Figure 11.. Pressure distribution over bump in channel using different spatial schemes. Figure 12 shows the application of adaptive mesh to the solution of this supersonic flow problem using 3rd order upwind.. First the initial condition is set as uniform flow, the same as the inlet condition. A large time step size of 500 seconds is applied. The grid iss adapted based on pressure gradient from an initial coarse grid solution. As time progresses, the flow field is converging toward its final solution and the shock becomes progressively sharper.



Figure 12.. Adaptive grid refinement for the supersonic pr problem using 3rd order upwind scheme D. Incompressible Flow Over a Generic Micro Micro-Munitions Wing The Generic Micro-Munitions (GENMM) GENMM) wing con configuration is shown in Figure 13.. Even though the wing is not rectangular in the spanwise direction, it has fairly cconstant onstant chord for most of the span. This allows the modeling of fluid dynamics from a 2D model. The chord length of the airfoil is chosen as 5”, and the incoming flow velocity is specified as U ∞ = 7 . 87 m / s (17.7 mi/hr). The fluid density is taken aas ρ = 1.0 kg/m3 and the dynamic viscosity of the fluid is 1.666x10-5 (kg/m-s). s). The Reynolds number is 60,000.

Figure 13.. GENMM wing configuration used for the present DNS study A quadrilateral grid is generated as shown Figure 114 [28], and the baseline grid has a total number of 100K cells. The solution adaptive refinement strategy will be applied here to resolve the fine scale vortical flow features. Second-order accurate spatial and temporal mporal schemes are employed employed. 15 American Institute of Aeronautics and Astronautics 092407


Figure 14. Computational mputational grid for the baseline 2D GENMM for CFD simulations

Figure 15.. The generation and transportation of leading edge vortex captured using adaptive mesh. (b) and(c) are the zoom zoom-in view of the flow field. Figure 15 shows ws the 2D results of generated vorticity around the GENMM at Re = 60,000 with zero angle-ofangle attack. At initialization, the entire flowfield is set to the freestream velocity. Unlike other types of airfoils, this GENMM has a curved mean camber line, which lleads eads to flow separation at the lower surface even at zero angle-ofangle attack (AoA). To efficiently and accurately capture and preserve the leading edge vortices, grid adaptation based on local vorticity is applied as shown in Figure 15. The vorticity field around ound the leading edge depicts a nice roll up of the leading vortex due to flow separation. This roll up in turn induces a secondary vortex with an opposite sense of 16 American Institute of Aeronautics and Astronautics 092407


rotation (marked in blue color) (Figure 15(b)). The generated vortices will then convect downstream wnstream with the mean flow. After that a new leading edge vortex is created again (Figure 15(c)). Figure 16 shows an instantaneous vortical flow field around the airfoil. A von Karman stre streett is evident at the trailing edge of the airfoil, and at the same time one can also see a Tollmien–Schlichting Schlichting (T (T-S) S) wave developing on the suction surface due to a streamwise instability in the viscous boundary layer. In the wake region, there are multiple shed vortices interacting with each other, leading to the unsteady dy nature of the wake flow. From the numerical simulation point of view, the vorticity strength is well preserved at the trailing edge indicating a low numerical dissipation of the present solution using adaptive mesh refinement. It is also observed that tthe he solution adaptation approach can well capture the formation and separation of a vortex. Once the transient flow sets in, the flowfield becomes highly unsteady.

Figure 16. Vortex shedding from us using ing adaptive mesh for GENMM at AOA=0 deg., and Re=60,000. E. 3D Transition Flow over SD7003 03 Airfoil The flow transition to turbulence is typically three dimensional. Hain et al. [29]] carried out water tunnel experiments for the SD7003 airfoil at a critical Reynolds number of 66,000, which is very close to the present consideration of 60,000. They hey employed time time-resolved resolved particle image velocimetry (PIV) for visualizing the developing flow structure. At 4 degree Angle of Attack (AoA) (AoA), a laminar separation bubble was observed. It was found that Kelvin-Helmholtz Helmholtz instabilities result in the generation and amplification mplification of spanwise vortices in the shear layer above the bubble, which then lead to three dimensional (3 (3-D) breakdown to turbulence. The simulation tion model is shown in Figure 17. The spanwise extent is 0.2c, and the present base-line base model has 50 cells in the spanwise direction. There are total of 600 grid points wrapping around the airfoil airfoil,, resulting in a total of 3.2 million cells.

Figure 17. 3-D D grid for SD7003 airfoil DNS. The spanwise extension =0.2C, with 50 cells. Total mesh size: 3.2 million.



3-D D Vortical Structure of the Simulated Flow Figure 18 shows the vortical flow field computed by the present DNS at a Reynolds number of 60,000 and angleof-attack attack (AoA) of 4 degree. At the leading edge, there is a shear layer developed due to the adverse pressure gradient. This shear layer shedss vortices as illustrated in Figure 18. Unlike the 2-D D case where the flow is constrained in the spanwise direction, however, the shear layer in 33-D D rolls up and experiences vortex burst. This vortex breakdown in 3-D D transitions the flow into a true three dimensional field.. As one can see, the flow eventually becomes turbulent after laminar to turbulence transition. This is substantially different from the 2D simulation in which the results remain essentially tially laminar.

Figure 18. 3-D D Vortical flow field computed by the present DNS at Re=60,000, α=4o.

study Figure 19. Computational grids around a SD7003 airfoil used for spatial sensitivity study. To appreciate the capabilities of grid adaption adaption,, grids with different resolutions were generated first to model the same physical problem. Figure 19 displays the four different grids without adaptation used in this study. The 18 American Institute of Aeronautics and Astronautics 092407


coarsest grid has about 1 million cells,, and the finest has about 9.7 million cells. It is expected that with the increase of the grid resolution, the capability in resolving subgrid subgrid-scale eddies will improve.. This is indeed the case as illustrated in Figure 20,, where the 3D vvorticity field represented in terms of Q factor is plotted.

D instantaneous vortical flow field computed by the present DNS at Re=60,000, α=4o using Figure 20. 3-D different grid resolutions and 2nd order spatial upwind scheme. Vorticity-Based Adaptation Q surface above the SD7003 wing with 1 million cells using adaptation based on Shown in Figure 21 is the iso-Q vorticity magnitude. Here the vorticity vvector is computed from the velocity field as:

Ω ij =

1  ∂u i ∂u j − 2  ∂x j ∂xi

   

(40)

In comparison to the results from the original grid in Figure 20,, much finer structures of vortical flow and its transient nature are captured with the adaptive mesh. It is also apparent that due to the unsteadiness, adaptation is required for each time step. Unlike the 22D case where the vorticity is coherent and well-structured, structured, in 3D the initial vortices from the shear layer break down into multiple sub sub-scale scale vortices and flow starts transition to turbulence. It is because of this turbulent nature that th the vorticity distribution in the transition region is very much random in nature.. The transition region, where Reynolds stress is non non-zero, zero, requires continuous adaptation. It is also noticed from the cross-sectional grid in Figure 221 that using vorticity as the adaption criterion leads to adaptation in i the near-wall wall region, even on the pressure side.

Figure 21.. Adaptation based on vorticity field for SD7003 wing, Re=60,000. AoA=4 deg using 3rd order upwind.



Q-Criterion-Based Adaptation As seen from Figure 21,, using vorticity as the adaption criteri criterion results in significant amount of grid cells being placed near the laminar boundaryy layer region on both the suction and nd pressure sides. From a flow physics point of view, it is the transition tion zone on the upper surface that contributes to the drag increase. se. The refinement is crucial in the non-zero mean Reynolds-stress stress zone instead of in the boundary layer. Actually, y, the mean Reynolds-stress Reynolds is zero within the laminar boundary layers next to wall. Shown in Figure 22 is the Q=0 iso-surface. surface. This surface is detached from the airfoil surface, and adaptation around this surface can be potentially more efficient. An adaptation scheme using the iso-Q Q criterion surface is developed and the results are shown in Figure 22. It is evident that using the Q=0 surface as the adaption criterion produces a grid that is only refined near the vortical core region. However, as the flow is unsteady, adaptation must be carried out at each time step. Mean Reynolds-Stress Stress Based Adaptation When the flow is fully ly developed, the m mean Reynolds-stress is invariant with time. Using mean Reynolds-stress Reynolds based adaptation can potentially save significant time as only a few adaptations are requir required. ed. To illustrate this, Figure 38 shows a 2D SD7003 airfoil rfoil at 8 degree AoA. Figure 223a displays the mean Reynolds-stress Reynolds distribution which is essentially independent dependent of time, and Figu Figure 23b displays the instantaneous vorticity distribution. The figures clearly show that the adaptation re region based on the mean Reynolds-stress stress covers the vorticity generation and transport regions.. The larger size of the adaptation zone is easily offset by the savings in CPU time. Also, one can see that the important physics occurs in this non-zero mean Reynolds-stress zone.

Figure 22. Grid adaptation based on Q=0 surface.



Figure 23. Distribution of mean Reynolds Reynolds-stress and grid id adaption using mean Reynolds-stress Reynolds for flow over a 2D SD7003 airfoil (a), and the vortical flow field (b) . Figure 24 shows the corresponding case for 3D SD7003 wing. Again, we se see that at the adaptation region encompasses the vortex breakdown physics. The nice feature of this adaptation is that the adaptation is static even for the unsteady flow field, and only a few adaptations are required. It is clear that the mean turbulence ReynoldsReynolds stress based grid adaption is the best approach for transition flow as it takes the turbulent flow physics into consideration, enabling the grid refinement to the most critical regions of the flow.

Figure 24.. 3D grid adaptation based on the mean Reyno Reynolds-Stress Stress distribution for a SD7003 wing. Re=60,000. AOA=4 deg. using 3rd order upwind

V. Conclusions High-order schemes of up to 4th-order order spatial and temporal accuracy for unstructured grids are combined with adaptive grid refinement to study unsteady eady incompressible vortical flow and compressible flow with shock waves. A novel high-order order scheme for unstructured grids is presented that, unlike nlike the popular discontinuous Galerkin method, does not introduce any additional degree degrees-of-freedom in each grid cell. Instead, the high-order order accuracy is achieved by using Taylor series expansions, and by computing high-order order derivatives as a function of lower order ones. This approach adds only very small overhead and can be easily implemented into existing production produ CFD codes. Several verification problems demonstrate that the scheme can achieve the designed formal order of accuracy. Unlike many of the existing adaptive mesh refinement methods based on pressure or velocity gradient, this approach utilizes first-principles physics-based quantities as the refinement criterion. These quantities include mean Reynoldsstress, vorticity magnitude and Q-criterion criterion. Transitional flow over a SD7003 airfoil is used to evaluate the most efficient adaptation method. It is shown own that the mean Reynolds-stress-based grid adaptation ion is the best approach for transitional flow as it takes the turbulent flow physics into consideration and focuses the grid refinement on the most critical regions of the flow. 21 American Institute of Aeronautics and Astronautics 092407

VI. Acknowledgements This study was supported under an Air Force SBIR Phase II contract number Contract: FA8651-11-C-0139. The contributions of Drs. X.G. Tan, A. Zhou, Robert Harris, Ravi Kanann and Sami Habchi of CFDRC are greatly appreciated.

VII. References


1

Marviplis, D. J. “Grid Resolution Study of a Drag Prediction Workshop Configuration using NSU3D Unstructured Mesh Solver,” 23rd AIAA Applied Aerodynamics Conference, Toronto, AIAA Paper 2005-4729, June 2005. 2Fraysse, F., Valero, E., and Ponsin, J., “Comparison of Mesh Adaptation Using the Adjoint Methodology and Truncation Error Estimates”, AIAA JOURNAL, Vol. 50, No. 9, September, 2012 3. Gao, H., and Wang, Z. J., “A Residual-Based Procedure for hp Adaptation on 2D Hybrid Meshes,” 49th AIAA Aerospace Sciences Meeting, Orlando, FL, AIAA Paper 2011-492, Jan. 2011. 4. Kamkar, S., Jameson, A., Wissink, A., and Sankaran, V., “Automated Off-Body Cartesian Mesh Adaption for Rotorcraft Simulations”, 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 4 - 7 January 2011, Orlando, Florida. AIAA-2011-1269. 5. Bin, J., Uzun, A., and Hussani, M.Y., “Adaptive Mesh Redistribution Method for Domains with complex boundaries”, J. of Computational Physics, Vol. 230, No. 8, pp.3178-3204, 2011 6 http://en.wikipedia.org/wiki/Overflow_(software) 7 Frink,N.T., Pirzadeh, S.Z, Parikh, P.C., Pandya, M.J., and Bhat, M.K., "The NASA Tetrahedral Unstructured Software System" , The Aeronautical Journal, Vol. 104, No. 1040, October 2000, pp. 491-499. 8 Luke, E., "On Robust and Accurate Arbitrary Polytope CFD Solvers (Invited)," AIAA 2007-3956, AIAA Computational Fluid Dynamics Conference, Miami, 2007. 9 http://fun3d.larc.nasa.gov/ 10 http://www.esi-group.com/products/Fluid-Dynamics/cfd-fastran 11 http://www.metacomptech.com/ 12 http://www.esi-group.com/products/multiphysics/ace-multiphysics-suite 13 http://www.ansys.com/Products/Simulation+Technology/Fluid+Dynamics/Fluid+Dynamics+Products/ANSYS +Fluent 14 http://www.cd-adapco.com/ 15 Luke, E., Thakur, S., Thomson, D., Wright, J., and Shyy, W., “Recent Progress Towards a Rule-Based Computational Tool for Liquid Rocket Combustion”, 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit 9 - 12 July 2006, Sacramento, California. AIAA 2006-5043. 16 Yang, H. Q., Habchi, S. D., and Przekwas, A. J., “General Strong Conservation Formulation of Navier-Stokes Equations in Nonothorthogonal Curvilinear Coordinates”, AIAA Journal, Vol. 32, No. 5, pp. 936-940, 1994. 17 Jiang, Y., Przekwas, A. J., “Implicit, Pressure-Based Incompressible Navier-Stokes Equations Solver for Unstructured Mesh”, AIAA-94-0303, 1994. 18 Chen, Z. J., Przekwas, A. J., “A Coupled Pressure-based Computational Method for Incompressible/Compressible Flows”, J. of Computational Physics, Vol. 229, pp. 9150-9165, 2010. 19 Yang, H. Q., Sheta, Essam, Habchi, S. D., Przekwas, A. J., Huang, A., “Novel Volume of Solid Technology for Nonlinear Aeroelastic Stability Analysis”, AIAA-2008-667. 20 Barth, T.J. and Frederickson, P.O., "Higher-Order Solution of the Euler Equations on Unstructured Grids Using Quadratic Reconstruction," AIAA-90-0013,1990. 21 Venkatakrishnan, V. and Mavriplis, D.J., "Implicit Solvers for Unstructured Meshes," AIAA-91-1537-CP, 1991. 22 Yang, H. Q., Chen, Z. J., and Dudley, J. G., “Development of High-Order Scheme in Unstructured Mesh for Direct Numerical Simulations”, Proceedings of ASME 2013 Fluid Engineering Division Summer Meeting, FEDSM2013-16385, July 07-11, 2013, Incline Village, Nevada, USA. 23 Yang. H. Q., Chen, Z. J., A. J. Przekwas, and Dudley, J. G., “High-Order Pressure-based Solver for Aeroacoustic Simulations”, 19th AIAA/CEAS Aeroacoustics Conference, Berlin, Germany, May 27-30, 2013, AIAA-2013-2021 24 Berger, M., Oliger, J. “Adaptive Mesh Refinement for Hyperbolic Partial Differential Equations”, J. of Computational Physics, V 53, 484–512, 1984.


25


Berger MJ., Colella P. “ Local Adaptive Mesh Refinement for Shock Hydrodynamics”, J of Computational Physics, V82, 6484 , 1989 26 Przekwas, A. J., Wilkerson, P., Chen, Z. J., Reeves, D., and Zhou, A., “Blast Wave and Acoustic Hearing Loss Injury Prevention and Protection”, Annual Report, 29 September 2010. CFDRC Project Number: 8943/Annual Year 1, 2010. 27 Fidkowski, K. J., “Review of Output-Based Error Estimation and Mesh Adaptation in Computational Fluid Dynamics,” AIAA Journal, Vol. 49, No. 4, 2011, pp. 673–694. 28 Yang. H. Q., Dudley, J. G. and David, J.,, “Aerodynamic Investigation of Generic SUAS”, 42nd AIAA Fluid Dynamic Conference and Exhibit, 25-28 June, 2012, New Orleans, LS, 2012, AIAA-2012-2969. 29 Hain, R, Kahler, C. J. and Radesoiel, R., “Dynamics of Laminar Separation Bubbles at Low-Reynolds Number Aerofoils”, Journal of Fluid Mechanics, Vol. 630, pp. 129-153, 2009.
