close

Вход

Забыли?

вход по аккаунту

?

6.2014-1277

код для вставкиСкачать
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
Документ
Категория
Без категории
Просмотров
2
Размер файла
458 Кб
Теги
2014, 1277
1/--страниц
Пожаловаться на содержимое документа