10.2514/6.2014-1277 AIAA SciTech Forum 13-17 January 2014, National Harbor, Maryland 52nd Aerospace Sciences Meeting A Fourth-Order Boundary Treatment for Viscous Fluxes on Cartesian Grid Finite-Volume Methods 1 1 Xinfeng Gao , Stephen Guzik 1 Computational Fluid Dynamics and Propulsion Laboratory, Colorado State University 2 Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 2 and Phillip Colella , Appied Numerical Algorithms Group, Lawrence Berkeley National Laboratory This study focuses on a fourth-order boundary treatment for nite-volume schemes to solve the compressible Navier-Stokes equations on a Cartesian grid. A fourth-order nite-volume stencil is derived for the viscous stress tensor operator and the modied fourth-order stencil near the physical boundary is developed. Fourier error analysis and stability analysis are performed for the fourth-order elliptic operator. For time integration, we use the fourth-order Runge-Kutta method. The fourth-order scheme was applied to the transient Couette ow and the solution accuracy was veried. Nomenclature Greek Symbols Notation ODE ordinary dierential equation PDE partial dierential equation RHS right-hand side kD k-dimensional space h ∆x, ∆y t time step ρ µ τ λ γ φ ξ spatial step time I. uid density uid viscosity uid stress tensor eigenvalues of the diusion model term coecients of four-stage Rung-Kutta a scalar quantity a length variable Introduction High order methods are required for applications such as compressible ows and combustion problems requiring long-time integration. In addition, high order methods can reduce impact of loss of accuracy at boundaries where mesh is not smooth (e.g., boundaries of mesh renement, boundaries between dierent mappings, etc.). In this study, we focus on a fourth-order boundary treatment for nite-volume schemes to solve the compressible Navier-Stokes equations on a Cartesian grid. For a compressible gas, the system of governing equations include the continuity, momentum, and energy equations given by. ∂ρ ~ + ∇·(ρ~u) = 0 ∂t ∂ ~ ρ~u~u +pI~~ = ∇· ~ ~~τ +ρf~ (ρ~u)+ ∇· ∂t h i ∂ ~ ρ~u e+ p ~ ~~τ ·~u − ∇·~ ~ q +ρf~·~u (ρe)+ ∇· = ∇· ∂t ρ where ρ is the density, ~u is the velocity, is the total specic energy with molecular stress tensor, and ~q h p is the pressure, being the enthalpy. ~ I~ (1) (2) (3) e = |~u|2 /2 + h − p/ρ molecular viscosity, ~ ~τ is the is the identity tensor, and In addition, µ is the is the molecular heat ux vector. The pressure is given by the ideal gas law. The molecular uid stress tensor, ~ ~τ , is dened as ~~ 1 ~~ ~ ~~τ = 2µ(S − I ∇ · ~u) 3 1 of 13 American Institute of Aeronautics and Astronautics Copyright © 2014 by Xinfeng Gao. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission. where ~ ~ S is the strain rate tensor dened by 1 ∂ui 2 ∂xj + ∂uj ∂xi . f~ is the body force. The molecular heat ux is modeled using Fourier's law. II. Fourth-Order Finite-Volume Approximations in Space The system of governing equations (Eqns.(13)) can be written in a conservative form as ∂U ~ ~ ~ = S. +∇· F−G ∂t (4) From the general divergence form, the integral form of the governing equations can be derived as Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 ∂ ∂t Z U dV Z + Vi ~ ~ ~ ∇ · (F − G) − S dV = 0 . Vi We consider Cartesian-grid nite-volume methods where the cell centers are marked by the points (i0 , ..., iD−1 ) = i ∈ ZD . D where x0 ∈ R is some Vi take the form Vi = [x0 + (i − 21 u)∆x, x0 + (i + 12 u)∆x] coordinates, ∆x is the mesh spacing, and u is the vector whose com- Thus, the control volumes xed origin of ponents are all equal to one. Representing the integrals as averages and applying the divergence theorem of Gauss, we arrive at a semi-discrete form D d 1 X ~ ~ ~ ~ hUii = − hFii+ 12 ed − hFii− 21 ed − hGii+ 12 ed − hGii− 21 ed +hSii , dt ∆xd (5) d=1 with h·i representing the cell-averaged or face-averaged quantities. The vectors of the conserved and the primitive solution variables for the three-dimensional case are given by W = [ρ, ux , uy , uz , T ] T T and , respectively. The primitive solution variables will be used in the next section. For completeness, the inviscid and viscous ux dyads ~ = F ~ G= U = [ρ, ρux , ρuy , ρuz , ρe] ρux ρu2x + p ρux uy ρux uz (ρe + p)ux 0 τxx τxy τxz ux τxx + uy τxy + uz τxz −κ ∂T ∂x , , ~ and G ~ are given in detail as F ρuy ρuy ux ρu2y + p ρuy uz (ρe + p)uy , ρuz ρuz ux ρuz uy ρu2z + p (ρe + p)uz 0 τyx τyy τyz ux τyx + uy τyy + uz τyz −κ ∂T ∂y , , 0 τzx τzy τzz ux τzx + uy τzy + uz τzz −κ ∂T ∂z respectively. The viscous stress tensor in 3D Cartesian coordinates is given by τxx τ = τxy τxz τyx τyy τyz ∂uy µ ∂uz x − 2( + ) 3 4 ∂u ∂x ∂y ∂z τzx ∂u x µ( ∂xy + ∂u τzy = ∂y ) τzz ∂ux z µ( ∂u ∂x + ∂z ) The source vector, S, x µ( ∂u ∂y + µ 3 4 ∂uy ∂y ∂uy ∂x ) x − 2( ∂u ∂x z µ( ∂u ∂y + z + ∂u ) ∂z ∂uy ∂z ) x µ( ∂u ∂z + ∂uz ∂x ) ∂uy ∂z ∂uz ∂y ) µ( µ 3 + ∂ux z 4 ∂u ∂z − 2( ∂x + ∂uy ∂y ) contains terms related to the body force, and the nite rate chemistry if a combustion problem is solved, and have the form of h i S = 0, ρfx , ρfy , ρfz , ρf~ · ~u . 2 of 13 American Institute of Aeronautics and Astronautics . II.A. The 4th -Order Hyperbolic Face-Averaged Flux on a Single-Level Grid Our approach to compute the the hyperbolic face-averaged ux, hFidi+ 1 ed , follows the procedures described 2 1 by McCorquodale and Colella. Some care was taken in transforming from conservative to primitive variables in order to preserve fourth-order accuracy. For completeness, we briey repeat the procedure here. 2 (2) hUii . Note hUii to hWii through the cell-centered value Ui by using Ui = hUii − ∆x 24 ∆ P 2 (2) (2) that ∆ is a second order accurate Laplacian (for a scalar quantity, q, ∆ qi = d 1/∆x (qi−ed −2qi + qi+ed ). The transformation between the conservative and primitive variables is given by Wi = W(Ui ) and Wi = W(hUii ). The latter is used to help compute cell-averaged primitive value, hWii = 2 (2) Wi + ∆x Wi . Note that the Laplacian is applied to Wi instead of Wi in order to minimize the 24 ∆ 1. Convert from size of stencil required, this substitution does not compromise the fourth-order accuracy. Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 2. Interpolate from hWii to hWii+ 21 ed hWii+ 12 ed = 3. Calculate face-centered Laplacian is using 1 7 (hWii + hWii+ed ) − (hWii−ed + hWii+2ed ) . 12 12 Wi+ 1 ed by Wi+ 1 ed 2 2 ∆(d,2) qi+ 12 ed = X d0 6=d 4. Compute face-averaged uxes hFidi+ 1 ed d , where the transverse 1 0 0 1 − 2q . q + q 1 d 1 d i+ 2 ed i+ 2 ed +ed ∆x2 i+ 2 e −e with 2 The superscript 2 (d,2) = hWii+ 12 ed − ∆x hWii+ 21 ed 24 ∆ (6) hFidi+ 1 ed = Fd (Wdi+ 1 ed ) + 2 2 ∆x2 (d,2) d ) . F (hWid 24 ∆ i+ 21 ed is used to denote the direction of the ux vector. Additional information for the modied stencils for the hyperbolic spatial discretization near physical boundaries, including issues on limiters and the dissipation mechanisms for strong shocks can be found in McCorquodale and Colella. II.B. The 1 4th -Order Elliptic Face-Averaged Flux on a Single-Level Grid hGidi+ 1 ed , we need to evaluate the face-averaged stress 2 ~ u + ∇~ ~ uT − 2 I~~∇ ~ · ~u with a varying viscosity. τ = µ(~x) ∇~ 3 In order to compute the elliptic face-averaged ux, tensor operator. The viscous tensor is dened by The components of the stress tensor operator consist of the normal and tangential velocity gradients. The 4th -order velocity gradient eld is constructed using the stencils as follows 1 + h~ u i − 15h~ u i − h~ u i 15h~ u i d d d i i−e i+2e i+e 12∆xd0 1 7 8h~uii+ed0 − 8h~uii−ed0 − h~uii+2ed0 + h~uii−2ed0 12∆xd0 12 ∂~u i 1 h +8h~uii+ed0 +ed − 8h~uii−ed0 +ed − h~uii+2ed0 +ed + h~uii−2ed0 +ed ∂xd0 i+ 2 ed 1 − 8h~uii+ed0−ed −8h~uii−ed0−ed −h~uii+2ed0−ed +h~uii−2ed0−ed 12 ! +8h~uii+ed0 +2ed −8h~uii−ed0+2ed −h~uii+2ed0+2ed +h~uii−2ed0+2ed 0 d =d 0 d 6= d (7a) (7b) The computation of the velocity gradient normal to the face, Eqn.(7a), is straightforward from a fourthorder centered stencil. The tangential velocity gradient calculation, Eqn.(7a), is formulated based on the following three steps: 3 of 13 American Institute of Aeronautics and Astronautics 1. In the stencil of four cells adjacent to face i + 12 ed (i.e., cells j = {i − ed , i, i + ed , i + 2ed }), evaluate the face-averaged solution variables 1 7 h~uij + h~uij+ed0 − huij−ed0 + huij+2ed0 + O(∆x4 ) 12 12 1 7 h~uij + h~uij−ed0 − huij−2ed0 + huij+ed0 + O(∆x4 ) . = 12 12 h~uij+ 1 ed0 = 2 h~uij− 1 ed0 2 (8) 2. Compute the cell-averaged tangential derivative of solution variables (specically, where the face direction is orthogonal to the derivative direction.) h ∂~u 1 0 0 ij = h~ u i − h~ u i . 0 j+ 12 ed j− 21 ed ∂xd0 ∆xd (9) Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 Therefore, the cell-averaged tangential derivative can be expressed by h 1 ∂~u 0 0 + h~ 0 − h~ 0 − 8h~ + O(∆x4 ) . u i u i u i 8h~ u i ij = 0 d d d d j−2e j+2e j−e j+e ∂xd0 12∆xd (10) 3. Calculate the face-averaged tangential derivative of solution variables using the centered interpolation scheme given by Eqn.(6) h ∂~u 7 i 1 = ∂xd0 i+ 2 ed 12 h ∂~u ∂~u i +h i ∂xd0 i ∂xd0 i+ed − 1 12 h ∂~u ∂~u i +h i ∂xd0 i−ed ∂xd0 i+2ed + O(∆x4 ) . (11) Now, substituting the cell-averaged tangential derivatives for the required cells into the equation, we arrive the complete stencil for calculating the face-averaged tangential velocity gradient. The viscous ux vector G can be expressed as a column vector of τ, i.e., G ≡ τ (:, c). The face-averaged elliptic ux is then evaluated by ∂~u hGidi+ 1 ed = Gd hµii+ 21 ed , h ii+ 1 ed . 2 2 ∂xd0 II.B.1. Fourier Error Analysis of the 4th -Order Spatial Operator We now conrm the 4th -order accuracy of the spatial discretization scheme using Fourier error analysis by applying the spatial dierencing operator to approximate the derivatives of the analytical Fourier function. The Fourier error analysis is performed for the cell-averaged tangential derivative using the stencil given by Eqn.(10) since its accuracy determines that of the stress tensor operator. Although the analysis is illustrated in two dimensions for simplicity, the extension to three dimensions is straightforward. Assume an arbitrary periodic function and express this function by the complex exponential part, i.e. W ~κ (~x) = e2πi~κ·~x , as f (~x) = ∞ X fˆ~κ W ~κ (~x) = −∞ where ~κ ∞ X fˆ~κ e2πi~κ·~x , −∞ is a vector of wavenumber and runs from −∞ to ∞ in each of the vector direction. First, let us f δ , i.e., discretization of f δ on a periodic domain by f = fˆ~κδ ω (~κ) , where δ is the spatial spacing, ∆x or ∆y . Simply, it is the −M 2 +1 case in this study that ∆x = ∆y on the same grid level. We use integer vector ~ j to index the grid points, since i is used for complex value notation; for example, in two dimensions, ~ j has two components (jx , jy ). review some denitions useful to the Fourier analysis. We dene a grid function, P M2 Then the discrete complex component is written as (~ κ) ω~j ≡ W ~κ (x~j ) ~ = e2πi~κ·jδ (jx = 0, 1, · · · , M −1; jy = 0, 1, · · · , M −1) . 4 of 13 American Institute of Aeronautics and Astronautics Also, an interpolation/spectral function is introduced as M 2 X fs (~x) = fˆ~κδ W ~κ (~x) , −M 2 +1 with its Fourier coecient taken from the grid function and its complex component taken from the analytical fs (~xj ) = f~j , function. The relation, can be derived from above. In addition, the analysis will use the shift operator and the cell-averaged value. The shift operator is dened as M M δ f~j+~ s 2 X = ~ e2πi~κ·(j+~s)δ fˆ~κδ = −M 2 +1 where ~s = m~x + n~y , ~e = ~x + ~y , and 2 X ~ eiβ fˆ~κδ e2πi~κ·jδ , −M 2 +1 β ≡ 2π(~κ · ~s)δ . For an arbitrary cell (jx , jy + n), where n is an integer, Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 the cell-averaged value takes the form of M 2 X hf ijx ,jy +n = ei2πκy n∆y sin β2x sin −M 2 +1 βx 2 βy 2 βy 2 ~ fˆ~κ∆y e2πi~κ·j∆y . Next, we illustrate the Fourier analysis on the cell-averaged tangential derivative calculation using stencil in Eqn.(10) for cell(jx , jy ). Substitute the cell-averaged values of the Fourier function for all the cells used in the stencil, we obtain h ∂f 1 ijx ,jy = 8hf ijx ,jy +1 − 8hf ijx ,jy −1 − hf ijx ,jy +2 + hf ijx ,jy −2 ∂y 12∆y M βy βx 2 X 1 i2πκy ∆y −i2πκy ∆y i2πκy 2∆y −i2πκy 2∆y sin 2 sin 2 ˆ∆y 2πi~ κ·~j∆y = 8e −8e −e +e f~κ e βx βy 12∆y M 2 − 2 2 +1 M 2 = X sin β2x sin β2y ∆y 2πi~κ·~j∆y 1 ˆ e f . 8eiβy − 8e−iβy − ei2βy + e−i2βy ~ κ βx βy 12∆y M 2 − 2 2 +1 Further, we use the Taylor series for eix = P∞ 0 (ix)k k! 2 3 4 5 = 1 + ix − x2! + i x3! + x4! + i x5! · · · to simplify the above equation as ∂f h ijx ,jy ∂y M 2 X = i2πκy ∆y + O(∆y 4 ) sin β2x sin βx 2 −M 2 +1 βy 2 βy 2 ~ fˆ~κ∆y e2πi~κ·j∆y M β βy βx 2 X sin β2x sin 2y βy ˆ∆y 2πi~κ·~j∆y 4 sin 2 sin 2 ˆ∆y 2πi~ κ·~j∆y i f e + O(∆y ) βx f~κ e . βx βy βy ∆y ~κ M 2 2 = − 2 2 +1 (12) 2 In fact, the rst term on the right-hand side of Eqn.(12) is the analytical derivative of the Fourier function. We can verify it by calculating an analytical average value of the derivative on an interval [~ xo h ∂f i~x ∂y = = = 1 h2 Z 1 h2 Z (jyo + 12 )h Z (jyo − 21 )h (jxo + 12 )h (jxo − 12 )h (jyo + 12 )h Z (jyo − 21 )h (jxo + 12 )h (jxo − 12 )h − h2 ~e, ~xo + h2 ~e] ∂f (~x) dxdy ∂y ∞ X i2πκy fˆ~κ e2πi(~κ·~x) dxdy −∞ ∞ β X sin β2x sin 2y βy ˆ 2πi~κ·~x i f~κ e . βx βy h −∞ 2 2 Further, we can verify that the second term on the right-hand side of Eqn.(12) is sin( β2x ) βx 2 = 1− 2 βx 24 4 βx + 1920 + · · · . Indeed, we have a th 4 ∂f -order spatial operator h ∂y ijx ,jy O(∆y 4 ) = h ∂f x ∂y i~ with the aid of + O(∆y 4 ). The ∂~ u analysis was also performed for h ∂xd0 ii+ 12 ed given by Eqn.(11) and a fourth-order accuracy was achieved. 5 of 13 American Institute of Aeronautics and Astronautics II.B.2. Modied Stencils near Physical Boundaries for the Face-Averaged Gradients Our strategy for computing uxes is to rst extrapolate the solution to two layers of ghost cells at the boundary and then use the centered interior scheme everywhere. Next, at physical boundaries where more information is known, the uxes are corrected to accommodate the extra information. For example, in a two- x-direction, the additional information is ~u = (0, 0) on the top wall and ~u = (U0 , 0) on the bottom wall. We generally ∂~ u ∂~ u know that h ∂x i = 0 on the boundary faces but need to correct h ∂y i and modify the y -direction uxes near ∂~ u the boundary accordingly. We can compute h ∂y i on the wall boundaries based on two approaches. The rst dimensional Couette ow with the innite walls aligned with the the velocities known on the wall, approach is a one-sided scheme by employing one layer of ghost cells near the physical boundary. The second approach is to employ two layers of ghost cells at the boundary and the centered interior scheme will be applied everywhere. An objective in both approaches is to avoid the necessity of nding information across corners of complex geometries. This implies that the boundary update scheme is necessarily one-dimensional in a direction orthogonal to the boundary. Details are presented for both schemes. Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 For the rst approach, we employ one layer of ghost cells. The values for the ghost cells are computed rst based on the known information on the physical boundary using a fourth-order interpolator/extrapolator. We illustrate the procedure in a one-dimensional case as shown in Figs. 1(a) 1(c). j j j+ 1 2 j+ 7 2 j− 1 2 j+ 5 2 j− 3 2 j+ 3 2 j− 5 2 j+ 1 2 j− 7 2 j− 1 2 j−1 j−2 j−3 (a) Ghost cell on j (b) High side Consider a moving (c) Low side high side Figure 1. Modied stencils near physical boundaries. wall boundary on the top side where the wall velocity is known as j be the ghost cell as shown in Fig. 1(a) and computing huij ∆y huiw , a face-averaged quantity. Let cell is the spatial spacing. The mathematical procedures for are described as follows. We start with the assumption that the face average can be interpolated from the indenite integrals of u in each cell U= R ξ u(ξ)dξ such that U (ξi+ 21 ) = Ui+ 12 = X huik ∆ξk , (13) k>i 1 2 is used to represent any face in Fig. 1(a) (i = {j, . . . , j − 4}). Note here that ξ and y are dened to increase moving away from the wall towards the interior. The interpolating function U (ξ) = a + bξ + cξ 2 + dξ 3 + eξ 4 is dened from the ve faces, starting at the location of j + 21 in the stencil shown where index i+ in Figs. 1(a) and 1(b): Uj+ 12 Uj− 12 Uj− 32 Uj− 52 Uj− 72 =0 = huij ∆y = huij + huij−1 ∆y . (14) = huij + huij−1 + huij−2 ∆y = huij + huij−1 + huij−2 + huij−3 ∆y Replacing the left-hand side of above equation with appropriate values of interpolating function and rearranging the equations, we can solve a linear system, Ax = b, for the coecients in the interpolating 6 of 13 American Institute of Aeronautics and Astronautics function 1 0 1 ∆y 1 2∆y 1 3∆y 1 4∆y 0 (∆y)2 (2∆y)2 (3∆y)2 (4∆y)2 0 (∆y)3 (2∆y)3 (3∆y)3 (4∆y)3 0 (∆y)4 (2∆y)4 (3∆y)4 (4∆y)4 a b c d e = 0 huij ∆y huij + huij−1 ∆y huij + huij−1 + huij−2 ∆y huij + huij−1 + huij−2 + huij−3 ∆y . (15) huif = b + Uj− 12 = huiw and using Dierentiating the resulting interpolant yields the formula for computing face-averaged values, 2cξ + 3dξ 2 + 4eξ 3 . This is the approach used to derive Eqn.(6). However, knowing the formula, we can instead compute the cell-averaged value for the ghost cell Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 at j− j from the face-averaged value 1 2 and the cell-averaged values involved in the stencil by 1 5 13 huij = 4 huiw − huij−3 + huij−2 − huij−1 . 12 12 12 (16) j , one will need to update the uxes for the face j − 23 in Fig. 1(b), or 3 j + 2 in Fig. 1(c)), since the cell j is involved in the fourth-order centered stencil for that face. It is likely also necessary to update the faces orthogonal to the y -direction near the boundary, probably starting with With the value of the ghost cell Eqn.(9), using the known velocity on the boundary face and recomputing the velocity on the next interior face with one of Eqns.(8) using the recently determined value in ghost cell have only analyzed Couette ow which is uniform in the x-direction j. Note that in our results we and thus does not permit investigating this particular detail. Knowing the value in the ghost cell, the next step is to use this value in a one-sided fourth-order stencil for computation of the velocity gradients on physical boundaries. We now describe the procedures for ∂U (ξ) = hu(ξ)i resulting from ∂ξ − 12 ]∆y as shown in Fig.(1(b)), obtaining the modied stencils near physical boundaries. Recall the relation the denition of the indenite integrals. For huij over the interval [j + 1 2, j we have δUj (y) = Uj− 21 − Uj+ 12 = −∆yhuij , where now y is dened to increase in the same direction as j . We perform Uj− 52 , Uj− 27 at the location of (j − 12 )∆y in Fig. 1(b) and arrive at Uj+ 12 = Uj− 12 U 3 =U 1 j− 2 j− 2 Uj− 5 = Uj− 1 2 2 Uj− 27 = Uj− 12 Taylor expansions of Uj+ 21 , Uj− 32 , (∆y)2 ∂ 2 U (∆y)3 ∂ 3 U (∆y)4 ∂ 4 U 5 2 ∂y 2 + 3! ∂y 3 + 4! ∂y 4 + O(∆y ) 2 3 4 2 3 4 (−∆y) ∂ U (−∆y) ∂ U (−∆y) ∂ U 5 (−∆y) ∂U ∂y + 2 ∂y 2 + 3! ∂y 3 + 4! ∂y 4 + O(∆y ) 2 3 4 2 3 4 (−2∆y) ∂ U (−2∆y) ∂ U (−2∆y) ∂ U 5 (−2∆y) ∂U ∂y + 2 ∂y 2 + 3! ∂y 3 + 4! ∂y 4 + O(∆y ) 2 3 4 2 3 4 (−3∆y) ∂ U (−3∆y) ∂ U (−3∆y) ∂ U 5 (−3∆y) ∂U ∂y + 2 ∂y 2 + 3! ∂y 3 + 4! ∂y 4 + O(∆y ) + ∆y ∂U ∂y + + + + . (17) We introduce a simple mathematical manipulation to correlate the integral value with the cell-averaged value as follows Uj+ 12 U 3 j− 2 Uj− 25 U 7 j− 2 − Uj− 12 = ∆yhuij − Uj− 21 = −∆yhuij−1 (18) − Uj− 21 = −∆yhuij−2 − ∆yhuij−1 − Uj− 21 = −∆yhuij−3 − ∆yhuij−2 − ∆yhuij−1 Now, rewrite Eqn.(17) in a matrix form and use the correlations expressed in Eqn.(18), we have huij −huij−1 −huij−2 − huij−1 −huij−3 − huij−2 − huij−1 1 −1 = −2 −3 1 2 1 2 4 2 9 2 1 6 − 61 − 68 − 27 6 1 24 1 24 16 24 81 24 ∂U ∂y 2 ∆y ∂∂yU2 3 (∆y)2 ∂∂yU3 4 (∆y)3 ∂∂yU4 7 of 13 American Institute of Aeronautics and Astronautics j− 12 . (19) Again, the indenite integrals are dierentiable and result in h ∂U ∂y 2 3 iT 4 , ∆y ∂∂yU2 , (∆y)2 ∂∂yU3 , (∆y)3 ∂∂yU4 j− 12 = h 2 3 2 ∂ u 3 ∂ u hu(y)i , ∆yh ∂u ∂y i , (∆y) h ∂y 2 i , (∆y) h ∂y 3 i iT j− 12 . Consequently, Eqn.(19) can be solved for the following face-averaged derivatives hu(y)i ∆yh ∂u ∂y i (∆y)2 h ∂ 2 u i ∂y 2 3 (∆y)3 h ∂∂yu3 i = 1 4 11 12 3 2 − 32 1 j− 12 1 2 1 2 1 3 6 6 −3 −4 1 − 12 1 − 12 1 2 1 huij −huij−1 −huij−2 − huij−1 −huij−3 − huij−2 − huij−1 . (20) Therefore, the modied stencil for the rst derivative at a physical boundary near the top side in Fig. 1(b)) Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 takes the form of h ∂u 1 i 1 = ∂y j− 2 ∆y 11 9 3 1 huij − huij−1 − huij−2 + huij−3 12 12 12 12 (21) For the physical boundaries in the low side of the domain, we follow the same procedure as described above to obtain the corresponding modied stencil. Uj+ 32 , Uj+ 52 , Uj+ 72 at the location of (j hu(y)i ∆yh ∂u ∂y i (∆y)2 h ∂ 2 u i ∂y 2 3 (∆y)3 h ∂∂yu3 i + But, the Taylor expansions are performed for faces 1 2 )∆y . And the nal linear system is − 14 11 = 12 − 32 1 j+ 21 3 2 1 2 − 12 −6 6 3 −4 1 12 1 − 12 1 −2 1 3 1 −huij huij+1 huij+2 + huij+1 huij+3 + huij+2 + huij+1 Uj− 12 , . (22) Accordingly, the modied stencil for a physical boundary near the low side is h ∂u 1 i 1 = ∂y j+ 2 ∆y − 11 9 3 1 huij + huij+1 + huij+2 − huij+3 . 12 12 12 12 (23) The second approach is to use two layers of ghost cells. We compute the values for the two-layer ghost cells based on the boundary conditions and then apply the center scheme for calculating the viscous ten- j+1 j tion II.B.2 to compute cell-averaged values for the two layers of ghost cells. For the rst layer of cells shown in 1 5 13 huij = 4 huiw − huij−3 + huij−2 − huij−1 . 12 12 12 j−2 For the second layer of cells, we can use the same ap- j−3 proach and extrapolate a bit further away. However, recall that the polynomial we derive is for the indenite integral. When we extrapolate further, what we are cal- (a) Ghost cells on culating is the average solution in a larger rectangle of Figure 2. aries. Uj+ 23 = −huij+1 ∆y , j− 1 2 j− 3 2 j− 5 2 j− 7 2 j (∆y)2 −(∆y)3 (b) Faces on high side Two-layer ghost cells near physical bound- and the rst row of Eqn.(15) with −∆y 1 2 high side cells. One can replace the rst row of Eqn.(14) with 1 j+ j−1 Fig.2, we can again use Eqn.(16), h 3 2 j+1 sor operators on the interior cell faces and the physical boundaries. We follow the procedures described in Sec- j+ (∆y)4 ih i h i a = −huij+1 ∆y , 8 of 13 American Institute of Aeronautics and Astronautics to obtain 2 37 83 huij+1 + huij = 4 5huiw − huij−3 + (huij−2 + huij−1 ) . 3 12 12 Since we just calculated huij , we can use this to now nd the average solution in cell (24) j + 1. Finally, the centered scheme is applied to compute the velocity gradients and viscous uxes on interior cell faces and physical boundaries. III. III.A. Stability Analysis Eigensystem Analysis of the Matrix Dierence Operators ∂φ x, t) is ∂t = −a∇φ + ν∆φ, where φ(~ is the propagation speed, a real constant which may be positive or negative, and ν is Let us consider an advection-diusion equation on a periodic domain a scalar quantity, a Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 the diusivity, a constant. Applying the fourth-order nite-volume formulation procedure as described in Section II to the advection-diusion equation, we obtain a set of coupled, homogeneous, linear, rst-order ODE's with constant coecients. Represent them by the equation D D X d ~ a X ~ d+ ν ~ d hφi = − Bp (1, −8, 0, 8, −1)hφie Bp (−1, 16, −30, 16, −1)hφie 2 dt 12∆x 12∆x d=1 d=1 (25) ∆x is the uniform spatial spacing. The periodic ma1 B 12 p (1, −8, 0, 8, −1) is the advection matrix operator and 1 B (−1, 16, −30, 16, −1) is the diusion matrix operator. The p 12 where trix Im eigenvalues of the advection matrix operator, λd λa 's, are all pure imaginary and those of the diusion matrix operator, 1.372 real and negative. We nd that λa |λa | ∈ (0, 1.372) and λd 's, are |λd | ∈ (0, 5.333) as illustrated in Fig. 3. The eigenvalues play important Re role in the stability contour that will be discussed next. -5.333 In practical application, the governing equations are typically -1.372 nonlinear. This is often used as a guide for estimating the worthiness of a method for more general problems and it serves as a Figure 3. λd . III.B. Distribution of eigenvalues λa and The fairly reliable necessary stability condition, but it is by no means a sucient one. 4th -order Runge-Kutta Scheme Having information on the eigensystems of the semi-discrete linear forms, we can decouple a complete system of ODE's into a set of linear, rst-order equations. We can derive the stability function of a time integration scheme. In this study, we apply the standard 4th -order Runge-Kutta scheme to a linear convection model ODE problem. For convenience, we describe the procedure for stability function derivation using the standard RK4, written in the form k i = h F φn + s X αij kj , tn + γi h j=1 φn+1 = φn + s X β j kj (26) βj are parameters given in a Butcher tableau j=1 where φ is a solution quantity, (Table 1), and s h is the time step, αij , γi and is four, the total number of stages. th Application of the standard 4 -order Runge-Kutta to a model ODE for a linear convection problem, Ps dφ ~ e + λhA~k . Rearrange the j=1 αij kj . Express it in vector form by k = λhφn~ dt = λφ, yields ki = λh φn + equation as ~k = I − λhA −1 λhφn~e , 9 of 13 American Institute of Aeronautics and Astronautics (27) Table 1. the Butcher tableau for the standard Runge-Kutta scheme. γ1 γ2 γ3 γ4 0 α21 α31 α41 β1 where = α32 α42 β2 α43 β3 1 k4 0 1 0 0 1 1/6 1/3 1/3 0 α A = 21 α31 α41 1/2 0 0 α32 α42 1/6 0 0 . 0 0 0 0 α43 0 Substitute Eqn.(27) into Eqn.(26), one arrives h −1 i φn+1 = φn + ~b T ~k = 1 + ~b T λh(I − λhA ~e φn , ~b T = [β1 , β2 , β3 , β4 ] k R(λh) kp ≤ 1. That is where and let R(λh) = 1 + ~b T λh(I − λhA k R(λh) kp −1 ~e. (28) A sucient condition for stability is det(I − λhA + λh~e~b T ) ≤ 1. det(I − λhA) = dφ dt = λφ + iχφ, where φ is a scalar are eigenvalues associated with the diusion and convection terms, we obtain the Apply the same procedure to a linear convection-diusion ODE problem, solution quantity, λ and χ stability function of det I − (λh + iχh)A + (λh + iχh)~e~b T k R(λh + iχh) kp = ≤ 1. det I − (λh + iχh)A Stability Region for Standard RK4 3 2 2 1 1 0 −1 −2 −3 −3 (29) Stability Region for Standard RK4 3 Im(λ h + i χ h) Im(λ h + i χ h) Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 1/2 1/2 β4 1 1 ~e = , 1 k1 k2 ~k = , k3 1/2 0 −1 −2 −2.5 −2 −1.5 −1 Re(λ h + i χ h) −0.5 0 −3 −0.25 −0.2 0.5 −0.15 −0.1 −0.05 0 0.05 Re(λ h + i χ h) 0.1 0.15 0.2 0.25 Figure 4. Stability contour (white region) for the standard Runge-Kutta methods generated from Eqn. (29). The plot on the right is a close-up view near the origin of the contour plot on the left. Figure 4 indicates the stability region of the standard RK4 method. This was used for time step calculations in this study. For problems dominated by strong diusion, we can derive the time constraint from the limiting value of |λh| ≤ 2 ν 2.5∆x λh = λd ∆x . For problems 2 hD and the time spacing h ≤ ν|λ | d max D limiting |χh| ≤ 2.937, where χh = λa D CFL with CFL = ah/∆x. 2.785, where dominated by convection, we see the Accordingly, the CFL number is limited by 2.937 D|λa |max . 10 of 13 American Institute of Aeronautics and Astronautics Although it is out of the scope of this paper, it is worth mentioning that we are concerned with the numerical integration of a sti system of ODE's (Eqn.(5)) due to diusion and combustion, since this is of interest to our future study. The use of explicit time schemes is often inecient. To cope with the stiness, we will adopt a linearly implicit Runge-Kutta method. Linearly implicit Runge-Kutta methods are a particular instance of additive Runge-Kutta (ARK) methods. We are particularly interested in the ARK4(3)6L[2]SA method, where the number is the number of stages, the 4 is the order of the main method, 3 is the order of the embedded method, 6 L represents a L-stable method, and SA denotes stiy accurate. The stability function is k R(λh + iχh) kp = det(I − λAI − iχAE + (λ + iχ)~e~b T ) det(I − λAI − iχAE ) λ and χ are eigenvalues from the matrix dierence operators for diusion and convection, respectively, AI , AE , and ~b are coecients for the ARK4 scheme. Zhang, Johansen and Colella2 investigated the where and stability region for a convection-diusion system in detail. More details on ARK schemes can be found in 3, 4 Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 Calvo and Kennedy. IV. Results We apply our algorithm to solve a transient Couette ow and compare the numerical solution to the exact analytical solution for numerical error convergence study. The analytical transient solution of a Couette ow between two innite plates may be found in the texts by Batchelor, the form ( u(y, t) = U0 where n 5 6 Sherman, or White, 7 and is often in 8 H is the channel height, ν ∞ X y sin(nπy/H) − n2 π22 νt (1 − ) − 2 e H , H nπ n=1 ) is the kinematic viscosity, U0 is the moving velocity of the bottom wall, is an integer, y is location normal to the wall moving direction, and t is time. In this paper, the uid is 80 ◦ C, air at 2 3 with a kinematic viscosity of 2.09e-05 m /s and a density of 1.0 kg/m . The channel height is 0.1m and the moving wall velocity is 0.0418 m/s. The Reynolds number is 200 based on the channel height and the moving wall velocity. The problem has an exact analytical so- 1 lution allowing for error to be Exact described 0.9 by a particular norm Computed given by R Lp = |u − uexact |p dΩ Ω dΩ where 0.8 p1 , 0.7 is the area 0.6 associated with a discrete By error measure. evaluating a se- 5 0 0 y H 10 0 0.5 s quence of grids, the er- 0.4 ror can be t to the relation size β ND in 32 0.3 Lp = αND β , where 8 D, 0.1 denotes the spatial −4) and the absolute α 0 β ≈ describes 0 0.1 0.2 0.3 magni- tude of the error. s s 1 order of accuracy (for the fourth-order, s 16 0.2 is the grid dimension s 0.4 s 0.5 u U0 0.6 0.7 0.8 0.9 1 Figure 5. Comparison of the computed and exact solution for the transient Couette ow at time = 1.0s, 8.0s, 16.0s, 32.0s, 100.0s, and 500s. The steady-state exact solution is u(y) U0 =1− y h. 11 of 13 American Institute of Aeronautics and Astronautics Our algorithm solves for cell- ∞ averaged values of the solution variables. For the error calculation, the −4 2 10 exact solution should be integrated for 1 each cell as ∞ ( 10−5 2 y )+ hu(y, t)iyi = U0 (y − 2H )yi + ∆y 2 ∞ X cos(nπy/H) − n2 π22 νt H e 2H . ∆y n2 π 2 yi − Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 determined once the addition of extra terms through a cycle of 2π in the sine operator no longer aected the precision. Note that we found it impracti- cal to solve the above equation analytically for t ≤ 0.1, 1 1 10−6 10−7 2 The number of terms in the series was 2 2 Solution Error n=1 ∞ ∞ ∞ 2 2 dy = −5 dx 1 1 ∞ ∞ 2 10−8 2 1 1 ∞ 10−9 ∞ 2 since a large number 1 2 of terms in the sum are required. The computed transient L∞ t0 = 0 s L2 t0 = 0 s L1 t0 = 0 s L∞ t0 = 0.48 s L2 t0 = 0.48 s L1 t0 = 0.48 s ∞ 1 10−10 Couette dy = −4 dx ow solution on the 64 grid and the an- 2 1 alytical solution are compared for time at 1.0, 8.0, 16.0, 32.0, 100.0, and 500.0 ∞ ∞ 10−11 2 seconds and plotted in Fig. 5. The com- 1 2 puted and the exact solutions overlap. The solution reaches steady state after 1 10−12 about 500 seconds. Numerical values of the transient Couette ow solution errors are measured with the L∞ -, L1 -, and 10−13 L2 -norms 16 32 64 for solution at 32 seconds. The solution 128 256 ND errors and the convergence rates beFigure 6. Solution accuracy for transient Couette ow. tween consecutive grid resolutions are tabulated in Tables 2 and 3. Table 2. Numerical values of the transient Couette ow solution errors measured with the L∞ -, L1 -, and L2 -norms at 32 seconds and the convergence rates between consecutive grid resolutions. The ow eld was initialized with zero. L norm L∞ L1 L2 32 rate 64 rate 128 rate 256 rate 4.711e-08 3.859 3.119e-09 3.917 2.007e-10 3.958 1.273e-11 3.979 6.621e-09 3.951 4.210e-10 3.975 2.657e-11 3.986 1.669e-12 3.987 1.550e-08 3.933 9.995e-10 3.955 6.355e-11 3.975 1.669e-12 3.992 Table 3. Numerical values of the transient Couette ow solution errors measured with the L∞ -, L1 -, and L2 -norms at 32 seconds and the convergence rates between consecutive grid resolutions. The ow eld was initialized with an exact solution at 0.48 seconds. L norm L∞ L1 L2 32 rate 64 rate 128 rate 256 rate 7.571e-07 4.110 2.881e-08 4.716 9.207e-10 4.968 2.880e-11 4.998 1.158e-07 4.109 4.457e-09 4.699 1.444e-10 4.947 1.669e-12 5.060 2.580e-07 4.110 9.893e-09 4.705 3.189e-10 4.955 9.676e-12 5.042 The dierence between the two tables is the ow initialization. Unlike the cases in Table 2 where the ow eld was initialized with zero velocity, the cases in Table 3 were initialized with an exact solution at 0.48 seconds. Figure 6 shows the norms of the solution error norms and the slopes for all the norms fall between -5 and -4. Therefore, the viscous operators are at least fourth-order accurate. 12 of 13 American Institute of Aeronautics and Astronautics The algorithm exhibits a super convergence rate of 5 when the ow led was initialized with an exact solution at 0.48 seconds. It is believed that a smooth initial ow eld contributed to the super convergence. V. Concluding Remarks and Future Work The fourth-order boundary treatment for a single-level algorithm to solve for the compressible NavierStokes equations is described in detail including the mathematical modeling, the nite-volume formulation, and the modied stencils near the physical boundary, and the numerical stability analysis. The single-level algorithm is veried with a fourth-order accuracy in both time and space. In future work we will explore multidimensional ows that permit full investigation of our proposed corrections for tangential derivatives on faces perpendicular to boundaries. Additionally, we will extend the AMR algorithm for viscous and combustion problems. Downloaded by UNIVERSITY OF ADELAIDE on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2014-1277 References 1 McCorquodale, P. and Colella, P., A high-order nite-volume method for conservation laws on locally rened grids, Comm. App. Math. Comput. Sci., Vol. 6, No. 1, 2011, pp. 125. 2 Zhang, Q., Johansen, H., and Colella, P., A fourth-order accurate nite-volume method with structured adaptive mesh renement for solving the advection-diusion equation, SIAM J. Sci. Comput., Vol. 34, 2012, pp. B179B201. 3 Calvo, M. and Novo, J., Linearly implicit Runge-Kutta methods for advection-reaction-diusion equations, Applied Numerical Mathematics , Vol. 37, 2001, pp. 535549. 4 Kennedy, C. and Carpenter, M., Additive Runge-Kutta schemes for convection-diusion-reaction equations, Applied Numerical Mathematics , Vol. 44, 2003, pp. 139181. 5 Batchelor, G. K., Introduction to Fluid Dynamics , Cambridge University Press, 1970, pp. 193e195. 6 Sherman, F., Viscous Flows , McGraw-Hill, 1991, pp. 335336. 7 White, F. M., Viscous Fluid Flow , McGraw-Hill, 1991, pp. 132143. 8 Muzychka, Y. and Yovanovich, M., Unsteady viscous ows and Stoke's rst problem, International Journal of Thermal Sciences , Vol. 49, 2010, pp. 820828. 9 Schlichting, H., Boundary-Layer Theory , McGraw-Hill, Toronto, 7th ed., 1979. 10 Gao, X., Northrup, S., and Groth, C., Parallel Solution-Adaptive Method for Two-Dimensional Non-Premixed Combusting Flows, Progress in Computational Fluid Dynamics, An International Journal , Vol. 11, No. 2, 2011, pp. 7695. 11 Guzik, S. M. J., McCorquodale, P., and Colella, P., A Freestream-Preserving High-Order Finite-Volume Method for Mapped Grids with Adaptive-Mesh Renement, AIAA 2012-0574, 50th AIAA Aerospace Sciences Meeting, 2012. 13 of 13 American Institute of Aeronautics and Astronautics

1/--страниц