Amano R.S., Sunden B. (Eds.) Computational Fluid Dynamics and Heat Transfer: Emerging Topics
Подождите немного. Документ загружается.
Sunden CH005.tex 10/9/2010 15: 0 Page 176
176 Computational Fluid Dynamics and Heat Transfer
r
s
4
3
1
2
s 1
s 1
r 1
r 1
Figure 5.1. Two-dimensional quadrilateral element.
laminar flow: U, V, P, T and six degrees of freedom for turbulent flows: U, V, P,
T, K, ω. Quadrilateral element should have nine nodes if midside nodes are used
(Figure 5.1).
Take this element as an example for using shape function and Jacobian matrix,
the velocity component u
i
and temperature T are approximated by using bilinear
interpolation function,
ϕ = ϑ =
⎡
⎢
⎢
⎢
⎢
⎣
1
4
(1 −r)(1 − s)
1
4
(1 +r)(1 − s)
1
4
(1 +r)(1 + s)
1
4
(1 −r)(1 − s)
⎤
⎥
⎥
⎥
⎥
⎦
(1)
Two pressurediscretizations arepossiblein thiskindof element:a bilinearcon-
tinuousapproximation, ψ ∈Q
1
, withthepressure degreesof freedomlocatedat the
four corner nodes, or a piecewise constant discontinuous pressure approximation,
ψ =1∈Q
0
,withthepressuredegreesoffreedomassociatedwithelementcentroid.
The construction of the finite-element matrices requires the computation of
various derivatives and integrals of the bilinear inter polation function. Since the
basis functions are given in terms of the normalized coordinates r and s and the
derivatives andintegrals are interms ofthe physicalx, y coordinates,the following
relations need to be defined:
⎡
⎢
⎢
⎣
∂ϕ
∂r
∂ϕ
∂s
⎤
⎥
⎥
⎦
=
⎡
⎢
⎢
⎣
∂x
∂r
∂y
∂r
∂x
∂s
∂y
∂s
⎤
⎥
⎥
⎦
⎡
⎢
⎢
⎣
∂ϕ
∂x
∂ϕ
∂y
⎤
⎥
⎥
⎦
=
⎡
⎢
⎢
⎣
∂N
T
∂r
x
∂N
T
∂s
y
∂N
T
∂s
x
∂N
T
∂s
y
⎤
⎥
⎥
⎦
⎡
⎢
⎢
⎣
∂ϕ
∂x
∂ϕ
∂y
⎤
⎥
⎥
⎦
= J
⎡
⎢
⎢
⎣
∂ϕ
∂x
∂ϕ
∂y
⎤
⎥
⎥
⎦
(2)
whereJ is the Jacobian matrix; by inverting it, equation (2) provides the necessary
relation for the derivatives of the basis function as below:
⎡
⎢
⎢
⎣
∂ϕ
∂x
∂ϕ
∂y
⎤
⎥
⎥
⎦
= J
−1
⎡
⎢
⎢
⎣
∂ϕ
∂r
∂ϕ
∂s
⎤
⎥
⎥
⎦
(3)
Sunden CH005.tex 10/9/2010 15: 0 Page 177
Equal-order segregated finite-element method 177
In thecaseof integral evaluation, to completethe transformation fromphysical
coordinate to normalized coordinates, an elemental area for quadrilateral element
could be evaluated as follows:
dx dy =|J|dr ds (4)
where |J|=Determine of J.
Equations (3) and (4) allow any integral of functions of x and y to be expressed
asintegralsofrational functionsin ther and s coordinate system.Thisis extremely
important in evaluating element matrices, which is the basic step in finite-element
method.
The details for other elements are ignored in this chapter.
Triangular element
The three-node triangular element can be used to model either two-dimensional
Cartesian or axisymmetric geometries. Each node has four degrees of freedom
for laminar flow: U, V, P, T and six degrees of freedom for turbulent flows: U, V,
P, T, K, ω.
Triangular element should have six nodes if midside nodes are used.
5.2.2 Three-dimensional elements
Brick (hexahedral) element
Theeight-node brickelement isused tomodel three-dimensionalgeometries asthe
general base. Each node has five degrees of freedom for laminar flow, i.e.,V, W, P,
T and six degrees of freedom for turbulent flows, i.e.,V,W, P,T, K, ω.
Brick element has total 21 nodes if midside nodes are used; this is also used as
the base element for other degenerated elements with midside nodes.
Tetrahedral element
The four-node tetrahedral elementcan be used tomodel three-dimensional geome-
tries as a degenerated form from three-dimensional brick element; the mapping
relation is listed in Table 5.1. Each node has five degrees of freedom for laminar
flow, i.e., U,V,W, P,T and seven degrees of freedom for turbulent flows, i.e., U,V,
W, P, T, K, ω.
Tetrahedral element has 10 nodes if midside nodes are used.
Pyramid element
Thefive-nodepyramid elementcanbeusedtomodelthree-dimensionalgeometries
as adegenerated form from three-dimensional brickelement; the mapping relation
islistedinTable5.1.Eachnodehasfivedegreesoffreedomforlaminarflow, i.e.,U,
V, W, P,T and seven degrees of freedom for turbulent flows, i.e., U,V,W, P,T, K, ω.
Pyramid element has 13 nodes if midside nodes are used.
Sunden CH005.tex 10/9/2010 15: 0 Page 178
178 Computational Fluid Dynamics and Heat Transfer
Table 5.1. Mapping from degenerated three-dimensional elements
to computational nodal configuration
3
1
2
4
5
1
2
3
3
4
5
6
6
1
2
3
4
1
2
3
4
5
6
1
2
3
3
4
5
6
6
1
2
3
3
4
5
3
3
Degenerated elements Computational nodal configuration
3D tetrahedral element
12345678
12334433
12345678
12334533
12345678
12334533
3D pyramid element
3D wedge element
Wedge element
The six-node wedge (prism) element can be used to model three-dimensional
geometries as a degenerated form from three-dimensional brick element; the map-
pingrelationislistedinTable5.1.Eachnodehasfivedegreesoffreedomforlaminar
Sunden CH005.tex 10/9/2010 15: 0 Page 179
Equal-order segregated finite-element method 179
flow, i.e., U,V,W, P,T and seven degrees of freedom for turbulent flows, i.e., U,V,
W, P, T, K, .
Wedge element has 15 nodes if midside nodes are used.
5.2.3 Degenerated elements
In general finite-element code development, the three-dimensional six-node
wedges, five-node pyramids, and four-node tetrahedral are considered as degen-
erate form of eight-node brick element by collapsing nodes. Table 5.1 attempts to
demonstrate how degenerate elements are formed by collapsing nodes if midside
nodes are not used.The first row is node ID in brick element and the second row is
fordegeneratednodes; thetworowsmakeamappingbetweendegenerated element
nodes and brick element nodes. The high-order elements with midside nodes are
notcommonlyused forfinite element–basedfluid flowand heattransfer codes,and
this is discussed in later sections of this chapter.
In two-dimensional case, the three-node triangle element is degenerated from
four-node quadrilateral element in a similar manner.
5.2.4 Special elements (rod and shell)
In heat transfer of some finite-element codes, arbitrary thickness rod and shell
elements may be defined as those which allow heat conduction along the rod or in
the shell of the element, but offer no resistance to heat flow across the element.
Rod elements have only length and may be used in either two-dimensional or
three-dimensionalsimulations. Shellelementshaveboth lengthand widthandmay
only be used in three-dimensional simulations.
5.3 Governing Equations for Fluid Flow and Heat Transfer
Problems
The numerical modeling process begins with a physical model that is based on
a number of simplifying physical assumptions. These assumptions can be made
in light of the understanding of the physical processes that are involved. A math-
ematical model is built from the physical model and then a series of PDEs are
developed.
For most of the practical physical processes, these PDEs cannot be solved
analytically but must be solved numerically. Analytic solutions have been found
only for a small subset of all possible fluid flow problems. Nonlinear nature of
the governing equations can produce extremely complex flow fields and multiple
time-dependent solutions, even for simplest geometries. Normally, the nature of
themathematical modelwould dictatethetypeofnumerical modelemployedanda
reviewofthecur rentliteraturewouldsuggestthemostappropriatenumericalmodel
to use. A combination of the mathematical model and numerical model leads to
a set of governing equations. Before these governing equations can be solved,
Sunden CH005.tex 10/9/2010 15: 0 Page 180
180 Computational Fluid Dynamics and Heat Transfer
asolution algorithmshould bedevelopedthatcansolvethe equationscorrectlyand
efficiently.
In this section, the gover ning equations and auxiliary equations are focused
uponforthefluid flowandheattransferphenomenaand dynamicmovingboundary
problems.The details are discussed as following aspects:
(a) Formation of the set of nonlinear coupled PDEs for mass and momentum
balances and various types of boundary conditions.
(b) Handling method for four types of terms in the governing equations, including
transient, advection, diffusion, and source.
(c) Uniform stabilized GLS formulation, including SUPG method and PSPG
method, which demonstrates viable equal-order velocity and pressure inter-
polations.
5.3.1 General form of governing equations
In this section, thestrong form ofthefluid flow problemis presented using a linear
constitutive relationship, leading to a set of nonlinear coupled PDEs for mass and
momentum balances. Various types of boundary conditions are also introduced.
Let and (0,T) denote the spatial and temporal domains, with x and t repre-
sentingthe coordinatesassociatedwith and(0,T).Theboundary ofthedomain
may involve several internal boundaries.
Conservation of mass
In solving fluid flow problems, it is important that the mass is conserved in the
system.The continuity equation may be written as
∂ρ
∂t
+∇·(ρ
V) = 0 (5)
where the term ∂ρ/∂t stands for the change in density over time; it also describes
the shrinkage-induced flow in heat transfer solidification.
Conservation of momentum
The differential equation governing the conservation of momentum in a given
direction for a Newtonian fluid can be written as
∂(ρ
V)
∂t
+∇·(ρ
V
V) =∇·(µ∇
V) −∇P +
B +S
m
(6)
where µ is the dynamic viscosity, P is the pressure,
B is the body force per unit
volume, and S
m
is momentum source term.
The component form of equation (6) in 3D is
ρ
∂u
1
∂t
+ u
1
∂u
1
∂x
1
+ u
2
∂u
1
∂x
2
+ u
3
∂u
1
∂x
3
−µ
∂
2
u
1
∂x
2
1
+
∂
2
u
1
∂x
2
2
+
∂
2
u
1
∂x
2
3
+
∂p
∂x
1
−f
1
= 0
(7)
Sunden CH005.tex 10/9/2010 15: 0 Page 181
Equal-order segregated finite-element method 181
ρ
∂u
2
∂t
+ u
1
∂u
2
∂x
1
+ u
2
∂u
2
∂x
2
+ u
3
∂u
2
∂x
3
−µ
∂
2
u
2
∂x
2
1
+
∂
2
u
2
∂x
2
2
+
∂
2
u
2
∂x
2
3
+
∂p
∂x
2
−f
2
= 0
(8)
ρ
∂u
3
∂t
+ u
1
∂u
3
∂x
1
+ u
2
∂u
3
∂x
2
+ u
3
∂u
3
∂x
3
−µ
∂
2
u
3
∂x
2
1
+
∂
2
u
3
∂x
2
2
+
∂
2
u
3
∂x
2
3
+
∂p
∂x
3
−f
3
= 0
(9)
If the flow is turbulent, the time-averaged N–S equations have the same for m
as equation (6), except that the fluid viscosity is replaced by an effective viscosity,
µ
eff
, which is defined as
µ
eff
= µ +µ
t
(10)
Here µ is still the fluid dynamic viscosity while µ
t
is the turbulence viscosity.
Conservation of energy
The equation governing the conservation of energy in ter ms of enthalpy can be
written in general as
∂(ρH)
∂t
+∇·(ρ
VH) =∇·(k∇T) +S
h
(11)
whereH isthe specificenthalpy, k isthethermalconductivity, T isthe temperature,
and S
h
is the volumetric rate of heat generation. The first term on the right-hand
sideof equation(11) representstheinfluence ofconduction heattransfer withinthe
fluid, according to the Fourier’s law of conduction.
The energy equation could also be written in terms of static temperature; how-
ever, the enthalpy form of energy equation generally considers the multiphase
flow such as steam/water, moist gas flow, and phase change like solidification and
melting.
Transport equations
The general transport equation for advection–diffusion is described as
∂
∂t
(ρφ) +∇(ρ
V · φ) =∇(ρ ·∇φ) + S (12)
Here,φ isthephasequantityingeneralexpression, isthe diffusioncoefficient
in general expression, the first left-hand term is the temporal term, and the second
left-handterm isthe convectiveterm.Thefirst right-handterm isthe diffusive term
and the second right-hand term is the source term.
Boundary conditions
Different types of boundary conditions can be encountered in fluid flow problems.
The simplest one could be the Dirichlet boundary condition, which is given by
u = u
g
(13)
Sunden CH005.tex 10/9/2010 15: 0 Page 182
182 Computational Fluid Dynamics and Heat Transfer
Another boundary conditionis traction-freeboundary condition,which isoften
applied for the far downstream boundary of the flow domain and is given by
σ · n = 0 (14)
For some special cases, traction force can also be applied as boundary
conditions.
For solid boundaries or symmetric planes, the normal component of velocity
must be specified as “no penetration” condition, i.e.,
u ·n = u
n
= 0 (15)
Theshearstressatthesymmetricplanesmustbespecifiedasanaturalboundary
condition as follows:
µ
∂u
τ
∂n
= 0 (16)
To be more specific, for a Newtonian fluid and a straight edge, equation (16)
can be rewritten in terms of its normal and tangential components:
σ
n
= (n · σ) · n =−p +2µ
∂u
n
∂n
(17)
σ
τ
= (n · σ) · τ = µ
∂u
τ
∂n
+
∂u
n
∂τ
(18)
If the symmetry line is aligned with a Cartesian axis, equation (15) becomes
a Dirichlet boundary condition while equation (16) is satisfied by setting σ
τ
=0.
Thus, a symmetry line involves both Dirichlet and Neumann type of boundaries.
5.3.2 Discretized equations and solution algorithm
Thissectionis usedtoreducethegoverningPDEsto asetofalgebraic equations.In
thismethod,thedependentvariablesarerepresentedbypolynomialshapefunctions
over a small area or volume (element in finite-element method). These represen-
tations are substituted into the governing PDEs and then the weighted integral of
these equations over the element is taken where the weight function is chosen to
be the same as the shape function.The result is a set of algebraic equations for the
dependent variable at discrete points or nodes on every element.
General form of discretization equations
The momentum, energy, temperature, and turbulent equations all have the form of
a scalar transport equation.There arefourtypes of terms in the equation, transient,
advection, diffusion, and source. For the purpose of describing the discretization
Sunden CH005.tex 10/9/2010 15: 0 Page 183
Equal-order segregated finite-element method 183
Table 5.2. General meaning of transportation representation
φ Meaning C
φ
φ
S
φ
u X velocity 1 µ
e
ρg
x
−
∂p
∂x
+ Rx
v X velocity 1 µ
e
ρg
x
−
∂p
∂x
+ Rx
w X velocity 1 µ
e
ρg
x
−
∂p
∂x
+ Rx
T Temperature C
p
K Q
v
+ E
k
+ W
v
+ µ +∂p/∂t
k Turbulence kinetic 1 µ
t
/σ
k
µ
t
/µ − ρε +C
4
βg
i
∂T
∂x
i
/σ
t
energy
ε Turbulence dissipation 1 µ
t
/σ
e
C
1
µ
t
ε/k − C
2
ρε
2
/k+
rate
C
1
C
µ
C
3
βkg
i
∂T
∂x
i
/σ
t
methods, the variable is referred to as φ. Then the form of the scalar transport
equation is
∂
∂t
(ρC
φ
φ) +
∂
∂t
(ρuC
φ
φ) +
∂
∂t
(ρvC
φ
φ)
∂
∂t
(ρwC
φ
φ)
(19)
=
∂
∂x
φ
∂φ
∂x
+
∂
∂y
φ
∂φ
∂y
+
∂
∂z
φ
∂φ
∂z
+ S
φ
where
C
φ
is the transient and advection coefficient,
φ
is the diffusion coefficient,
and S
φ
are source terms.
The pressure equation can be derived from continuity equation by segregated
solution scheme, which is discussed in the next chapter. Equation (19) is a general
form of velocity, temperature, turbulent variables k,ε, and species. Table 5.2 lists
the corresponding equation representation for C
φ
,
φ
, and S
φ
.
The discretization process, therefore, consists of deriving the element matrices
to put together the matrix equation in simplified matrix for m as follows:
([A
Transient
e
] +[A
Advection
e
] +[A
Diffusion
e
]){φ
e
}={S
φ
e
} (20)
Galerkin method of weighted residuals is used to form the element integrals,
denoted by N
e
the weighting function for the element, which is also the shape
function.
Sunden CH005.tex 10/9/2010 15: 0 Page 184
184 Computational Fluid Dynamics and Heat Transfer
Temporal term
For transientanalyses, thetransient terms are discretizedusing animplicitor back-
ward difference method. Using the matrix algebra notation, a typical steady-state
scalar transport equation with momentum, energy, and turbulence variables can be
written:
A
ij
u
j
= F
i
where A
ij
contains the discretized advection and diffusion terms from the gov-
erning equations, u
j
is the solution vector or values of the dependent variable
(u,v,w,T,K,. ...) and F
j
contains the source terms.
The transient terms in the governing equations took the form ρ∂φ/∂t, where
φ represents the dependent variable (u,v,w,....). This term is discretized using a
backward difference method:
∂φ
∂t
=
φ
new
− φ
old
t
(21)
Add this term to the matrix equation above
(A
ij
+ B
ij
)u
new
j
= F
i
+ B
ii
u
old
j
(22)
where B
ij
is a diagonal matrix composed of terms like:
B
ii
=
1
t
&
N
i
ρd
e
(23)
This discretized transient equation must be solved iteratively at each time step to
determine all of the new variables (variable values at the latest time).
Advection term
There are two approaches to discretize this most challengeable term in fluid flow:
the monotone streamline-upwind (MSU) approach is the first-order accurate and
produces smooth and monotone solutions, and the SUPG approach is the second-
order accurate.
Monotonestreamlineupwind(MSU)approach. In theMSUmethod,theadvec-
tiontermishandledthroughamonotonestreamlineapproachbasedontheideathat
advectiontransportisalongcharacteristiclines. Itisusefultothinkoftheadvection
transport formulation in terms of a quantity being transported in a known velocity
field.Thevelocityfielditselfcouldbeenvisionedasasetofstreamlineseverywhere
tangent to the velocity vectors. The advection terms can be expressed in terms of
the streamline velocities.
In advection-dominated transport phenomena, one assumes that no trans-
port occurs across characteristic lines, i.e., all transfer occurs along streamlines.
Sunden CH005.tex 10/9/2010 15: 0 Page 185
Equal-order segregated finite-element method 185
Therefore, if one assumes that the advection term is expressed in streamline
direction:
∂(ρc
φ
v
x
φ)
∂x
+
∂(ρc
φ
v
y
φ)
∂y
+
∂(ρc
φ
v
z
φ)
∂z
=
∂(ρc
φ
v
s
φ)
∂s
(24)
then the advection term is constant throughout an element:
[A
advection
e
] =
∂(ρc
φ
v
s
φ)
∂s
&
Nd
e
(25)
Figure 5.2 illustrates the streamline across a brick element; this formulation
could be used for every element, each of which will have one node which gets
contributions from inside the element; the derivative is calculated using a simple
difference as follows:
∂(ρc
φ
v
s
φ)
∂s
=
(ρc
φ
v
s
φ)
up
− (ρc
φ
v
s
φ)
down
s
(26)
where
down=Subscript for value at downstream node
up=Subscript for value taken at the location at which the streamline through the
streamline through the upwind node enters the element. The value is interpolated
in terms of the nodes in between.
s=Distance from the upstream point to the downstream node.
The whole process goes through all the elements and identifies the downwind
nodes.The calculation is madebased on the velocities tofindwhere the streamline
through thedownwindnode came from.Weight factors arecalculated based on the
proximity of the upwind location to the neighboring nodes.
Streamline-upwind/Petrov–Galerkin (SUPG) approach. This method is used
to discretize the advection term and an additional diffusion-like perturbation term,
which acts only in the advection direction.
Forthree-dimensionalproblem, theperturbed interpolation functionisgiven by
W
i
= N
i
+
β
¯
h
2|
V|
u
∂N
i
∂x
+ v
∂N
i
∂y
+ w
∂N
i
∂z
(27)
Here,
V =(u,v,w) is the velocity field in three-dimensional domain,
h is the
averaged size of an element, which, for eight-node brick element shown in Figure
5.2, can be calculated as
¯
h =
1
|
V|
(
|h
1
|+|h
2
|+|h
3
|
)
(28)