Amano R.S., Sunden B. (Eds.) Computational Fluid Dynamics and Heat Transfer: Emerging Topics
Подождите немного. Документ загружается.
Sunden CH005.tex 10/9/2010 15: 0 Page 196
196 Computational Fluid Dynamics and Heat Transfer
Construct a function,
&
⎡
⎢
⎢
⎣
∂
∂u
∂x
N
∂x
+
∂
∂u
∂y
N
∂y
+
∂
∂u
∂z
N
∂z
⎤
⎥
⎥
⎦
d
(46)
=
&
∂u
∂x
∂N
∂x
+
∂u
∂y
∂N
∂y
+
∂u
∂z
∂N
∂z
d +
&
N
∂
2
u
∂x
2
+
∂
2
u
∂y
2
+
∂
2
u
∂z
2
d
Using Green–Gauss theorem,
&
⎡
⎢
⎢
⎣
∂
∂u
∂x
N
∂x
+
∂
∂u
∂y
N
∂y
+
∂
∂u
∂z
N
∂z
⎤
⎥
⎥
⎦
d =
&
N
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d
(47)
Equations (45), (46), and (47) are combined to get the new form of the
discretized x-momentum equation:
&
ρN
u
∂u
∂x
+ v
∂u
∂y
+ w
∂u
∂z
d =−
&
N
∂p
∂x
d +
&
NF
x
d
+
&
N
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d −
&
∂u
∂x
∂N
∂x
+
∂u
∂y
∂N
∂y
+
∂u
∂z
∂N
∂z
d
Move the last item in the right-hand side to the left-hand side, then
LHS = [A
u
]{u}=
&
ρN
u
∂u
∂x
+ v
∂u
∂y
+ w
∂u
∂z
+
µ
∂u
∂x
∂N
∂x
+ µ
∂u
∂y
∂N
∂y
+ µ
∂u
∂z
∂N
∂z
d
= ρ
&
N
i
u
∂N
j
∂x
+ v
∂N
j
∂y
+ w
∂N
j
∂z
d{u}+
&
µ
∂N
i
∂x
∂N
j
∂x
d
{u}
+
&
µ
∂N
i
∂y
∂N
j
∂y
d
{u}+
&
µ
∂N
i
∂z
∂N
j
∂z
d
{u}
= [C(u,v,w) + K
11
+ K
22
+ K
33
]{u} (48a)
where, C(U) is the advection item and K
11
,K
22
,K
33
are the diffusion items.
To obtain C(U), velocities u,v,w are considered the average values in an element.
RHS =−
&
N
i
∂p
∂x
d +
&
N
i
F
x
d +
&
N
i
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d =−M
14
{p}+{F
xA
}+{F
x
} (48b)
[A
v
] and [A
w
] could be achieved by similar approaches.
Sunden CH005.tex 10/9/2010 15: 0 Page 197
Equal-order segregated finite-element method 197
SUPG itemsThe SUPG weighting function proposed by Brooks is defined by:
W
k
=N
k
+¯p
k
Here the streamline-upwind contribution ¯p
k
is defined as:
¯p
k
=
˜
k
e
(u
e
1
)
2
+ (u
e
2
)
2
u
e
j
∂N
i
∂x
j
(49)
Only convective items are weighted with the SUPG functions for the inconsistent
approach,whichisbelievedtobemorestablethantheconsistentSUPGformulation
[10]. The extra items from ¯p
k
should be added into K
11
,K
22
,K
33
in equation (49)
in the implementation.
Equation of pseudo velocity (prediction phase)
In matrix form of the seg regated momentum equations:
[M
11
]{u}=−[M
14
]{p}+{F
xA
}+{F
x
} (50a)
[M
22
]{v}=−[M
24
]{p}+{F
yA
}+{F
y
} (50b)
[M
33
]{w}=−[M
34
]{p}+{F
zA
}+{F
z
} (50c)
where
[M
11
] = [C(u,v,w) + K
11
+ K
22
+ K
33
] +[SUPG(u,v,w)]
[M
22
] = [C(u,v,w) + K
11
+ K
22
+ K
33
] +[SUPG(u,v,w)]
[M
33
] = [C(u,v,w) + K
11
+ K
22
+ K
33
] +[SUPG(u,v,w)]
[M
14
] =
&
N
i
∂p
∂x
d
[M
24
] =
&
N
i
∂p
∂y
d
[M
34
] =
&
N
i
∂p
∂z
d
{F
xA
}=
&
N
i
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d
{F
yA
}=
&
N
i
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d
{F
zA
}=
&
N
i
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d
{F
x
}=
&
N
i
F
x
d
{F
y
}=
&
N
i
F
y
d
{F
z
}=
&
N
i
F
z
d
where is the fluid domain (volume); and is the flow boundaries (surface).
There are no SUPG items in {F
x
}{F
y
}{F
z
} for the inconsistent approach.
Sunden CH005.tex 10/9/2010 15: 0 Page 198
198 Computational Fluid Dynamics and Heat Transfer
The above three equations can be solved separately to obtain the velocity
estimate (prediction phase).
HatVelocity (V
hat
) Rewrite equations (50a), (50b), (50c):
a
uii
u
i
=−
i=j
a
uij
u
j
−
&
N
i
∂p
∂x
d + µ
&
N
i
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d +{F
x
}
i
(51a)
a
vii
v
i
=−
i=j
a
vij
v
j
−
&
N
i
∂p
∂y
d + µ
&
N
i
∂v
∂x
n
x
+
∂v
∂y
n
y
+
∂v
∂z
n
z
d +{F
y
}
i
(51b)
a
wii
w
i
=−
i=j
a
wij
u
j
−
&
N
i
∂p
∂z
d + µ
&
N
i
∂w
∂x
n
x
+
∂w
∂y
n
y
+
∂w
∂z
n
z
d +{F
z
}
i
(51c)
For constant pressure gradient in an element, rewrite equation (51a):
u
i
=
1
a
uii
⎡
⎣
−
i=j
a
uij
u
j
+ µ
&
N
i
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d +{F
x
}
i
⎤
⎦
−
1
a
uii
&
N
i
d
∂p
∂x
i.e.,
u
i
=ˆu
i
− k
p
i
∂p
∂x
(52a)
where
k
p
ui
=
1
a
uii
&
N
i
d (53a)
ˆu
i
=
1
a
uii
⎡
⎣
−
i=j
a
uij
u
j
+ µ
&
N
i
∂u
∂x
n
x
+
∂u
∂y
n
y
+
∂u
∂z
n
z
d +{F
x
}
i
⎤
⎦
(54a)
Rewrite equation (51b):
v
i
=
1
a
vii
⎡
⎣
−
i=j
a
vij
v
j
+ µ
&
N
i
∂v
∂x
n
x
+
∂v
∂y
n
y
+
∂v
∂z
n
z
d +{F
y
}
i
⎤
⎦
−
1
a
vii
&
N
i
d
∂p
∂y
(50)
i.e.,
v
i
=ˆv
i
− k
p
i
∂p
∂y
(52b)
Sunden CH005.tex 10/9/2010 15: 0 Page 199
Equal-order segregated finite-element method 199
where
k
p
vi
=
1
a
vii
&
N
i
d (53b)
ˆv
i
=
1
a
vii
⎡
⎣
−
i=j
a
vij
v
j
+ µ
&
N
i
∂v
∂x
n
x
+
∂v
∂y
n
y
+
∂v
∂z
n
z
d +{F
y
}
i
⎤
⎦
(54b)
Rewrite equation (51c):
w
i
=
1
a
wii
⎡
⎣
−
i=j
a
wij
w
j
+ µ
&
N
i
∂w
∂x
n
x
+
∂w
∂y
n
y
+
∂w
∂z
n
z
d +{F
z
}
i
⎤
⎦
−
1
a
wii
&
N
i
d
∂p
∂z
i.e.,
w =ˆw
i
− k
p
i
∂p
∂z
(52c)
where
k
p
wi
=
1
a
wii
&
N
i
d (53c)
ˆw
i
=
1
a
wii
⎡
⎣
−
i=j
a
wij
w
j
+ µ
&
N
i
∂w
∂x
n
x
+
∂w
∂y
n
y
+
∂w
∂z
n
z
d +{F
z
}
i
⎤
⎦
(54c)
To ensure mass continuity, the velocity should be updated (corrected) at each
iterationbywhich velocity correction isreformulatedfrom the new pressurevalue,
refer to equation 62a, 62b and 62c for velocity correction details.
V
Corrected
= V
hat
+ V
Correction_Term_By_Pressure
(55)
Pressure Poisson equation
On applying Green–Gauss theorem and integration to the continuity equation,
the integral of the divergence in a domain equals the net flux across the domain
boundary:
&
∂(uN)
∂x
+
∂(vN)
∂y
+
∂(wN)
∂z
d =
&
Nv ·nd (56)
The continuity equation becomes:
&
N
∂u
∂x
+
∂v
∂y
+
∂w
∂z
d =−
&
∂N
∂x
u +
∂N
∂y
v +
∂N
∂z
w
d
(57)
+
&
N(un
x
+ vn
y
+ wn
z
)d = 0
Sunden CH005.tex 10/9/2010 15: 0 Page 200
200 Computational Fluid Dynamics and Heat Transfer
The modified continuity equation (57) using green function transformation is
targeted to solve pressure.
Step 1: Substitute interpolated velocity from equations (44a), (44b), and (44c)
to continuity equation (57)
−
&
∂N
i
∂x
(N
j
u
j
) +
∂N
i
∂y
(N
j
v
j
) +
∂N
i
∂z
(N
j
w
j
)
d +
&
N
i
(un
x
+vn
y
+wn
z
)d = 0
(58)
Step 2: Substitute velocity prediction from equations (53a), (53b), and (53c)
into equation (58)
&
∂N
i
∂x
N
j
K
p
j
∂p
∂x
+
∂N
i
∂y
N
j
K
p
j
∂p
∂y
+
∂N
i
∂z
N
j
K
p
j
∂p
∂z
d
=
&
∂N
i
∂x
(N
j
ˆu
j
) +
∂N
i
∂y
(N
j
ˆv
j
) +
∂N
i
∂z
(N
j
ˆw
j
)
d −
&
N
i
(un
x
+ vn
y
+ wn
z
)d
(59)
Step 3: Rewrite to final form
Notice that p=N
i
p
i
, so substituting
∂p
∂x
=
∂(N
i
p
i
)
∂x
=p
i
∂N
i
∂x
(assume that p is con-
stant in an element here) in equation (59) (using similar way for and ) we obtain
∂p
∂y
and
∂p
∂z
.
$
&
∂N
i
∂x
(N
k
K
p
k
)
∂N
j
∂x
d +
&
∂N
i
∂y
(N
k
K
p
k
)
∂N
j
∂y
d
+
&
∂N
i
∂z
(N
k
K
p
k
)
∂N
j
∂z
d
%
·{p}
(60)
=
&
∂N
i
∂x
(N
j
ˆu
j
)d +
&
∂N
i
∂y
(N
j
ˆv
j
)d +
&
∂N
i
∂z
(N
j
ˆw
j
)d
−
&
N
i
(un
x
+ vn
y
+ wn
z
)d
Moving V
hat
out of the integration in equation (60), we obtain:
$
&
∂N
i
∂x
(N
k
K
p
k
)
∂N
j
∂x
d +
&
∂N
i
∂y
(N
k
K
p
k
)
∂N
j
∂y
d
+
&
∂N
i
∂z
(N
k
K
p
k
)
∂N
j
∂z
)
d
%
·{p}
(61)
=
&
∂N
i
∂x
N
j
d{ˆu}+
&
∂N
i
∂y
N
j
d{ˆv}+
&
∂N
i
∂z
N
j
d{ˆw}
−
&
N
i
(un
x
+ vn
y
+ wn
z
)d
Now equation (61) is the final form of pressure Poisson equation.
Sunden CH005.tex 10/9/2010 15: 0 Page 201
Equal-order segregated finite-element method 201
It is interesting to note that the element pressure matrices are identical to
those obtained in classical diffusion-type problems, with the term K replacing
the diffusion coefficient.
Here are more discussions on segregated solution scheme.
(1) The coefficient matrix for the pressure equation is similar to that obtained for
thediffusiontermintheconventionalfinite-elementformulationiftheviscosity
is replaced by N
i
K
p
j
the formulation developed by Rice [5] andVellando [12];
this fact is indicative of the stability and robust nature of the resulting pressure
field.
(2) The LHS and RHS of equation (69) could be constructed and evaluated in
element level and then assembled into global matrix LHS and RHS; this is
typical finite-element behavior.
(3) In equation (60), k
p
i
is a node-based proper ty, it should be assembled from
the saved element-level matrix,
/
N
i
d, a
uii
is also a global property from
assembly.
(4) Equation (60) is in the form of Poisson equation; the LHS matrix has the
properties of symmetric and positive definite as indicated by Shaw [7].
(5) Specialboundaryintegrationsonallinflowandoutflowboundariesarerequired
to solve pressure, which corresponds to the fourth item in equation (61).
Velocity correction
Write equations(51a), (51b), and (51c) usingweightingfunctionp(x,y)=N
i
p
i
and
chain law, and discarding velocity gradient items on boundaries. ˆu
i
, ˆv
i
, ˆw
i
are used
for the velocity prediction as “hat” velocity.
u
i
=
1
a
uii
⎡
⎣
−
i=j
a
uij
u
j
+{F
x
}
i
⎤
⎦
−
1
a
uii
&
N
i
∂N
j
∂x
d
{p}
Or, u
i
=ˆu
i
−
1
a
uii
&
N
i
∂N
j
∂x
d
{p}=ˆu
i
−
1
a
uii
{M
14
}{p} (62a)
Get corrected velocity v, w in the same way as u,
v
i
=ˆv
i
−
1
a
vii
&
N
i
∂N
j
∂y
d
{p}=ˆv
i
−
1
a
vii
{M
24
}{p} (62b)
w
i
=ˆw
i
−
1
a
wii
&
N
i
∂N
j
∂z
d
{p}=ˆw
i
−
1
a
wii
{M
34
}{p} (62c)
At this step, u,v,w are directly corrected from the predicted and with trivial
effort, the second item is the pressure coupling effect on the velocity.
Special treatments on boundary conditions
There are some special treatments in the segregated formulation.
Wall boundary: For wall boundary, the surface integral or natural boundary
conditionforthepressureequation iszero.Toincorporatethe knownvelocitiesinto
Sunden CH005.tex 10/9/2010 15: 0 Page 202
202 Computational Fluid Dynamics and Heat Transfer
pressure equation, it is required that “hat” velocity be zero. In implementation, the
pressure coefficients, k
p
ui
, k
p
vi
, k
p
wi
at the wall nodes should be set as zero.
Specified velocity boundary: The surface integration of velocity boundary
conditionsshouldbeaddedintothepressureequation(61)’srighthandsideasinthe
fourth item in the right hand side; for specified velocityboundary, the contribution
is fixed.
Inflow/OutflowboundaryZerogradientvelocityboundaryconditionisapplied
at the inflow/outflow boundaries (including inlet/outlet boundary and pressure
boundary), i.e., no special treatment in the momentum equation is required at
this kind of boundary. However, the contribution from the surface integration of
natural boundary conditions should be added to the right-hand side of the pressure
equation as done in the fourth item on the right-hand side of equation (60).
Segregatedschemenotonlysimplydividesthebigcoupling matrixintosmaller
matrices, but also sets up a pressure equation (Poisson equation) with positive-
definite symmetric properties, which can greatly expedite the pressure solution
with better accuracy. Most of the fluid flow problems are driven by pressure, and
thereforethevelocitycomponentswillfollowagoodsolutionofpressureinpractice.
Solution stability control
Inertial relaxation. Iner tial relaxation is used in the governing equations to slow
downtheconvergencerateinthesamemannerastransittermsareusedinnonsteady
problems, so the inertial relaxation factor is also called virtual time step.
The inertial relaxation is added to the governing equations in the following
manner:
A
ii
+
ρ
i
/
Nd
t
inertia
φ
i
+
j=i
A
ij
φ
j
= F
i
+
ρ
i
/
Nd
t
inertia
φ
iold
(63)
The second term in parentheses and the last term on the right-hand side are the
inertial relaxation terms. t
inertia
can be automatically adjusted according to the
global matrix condition number.
•
Small t
inertia
value will enhance the solution stability and slow down the
convergence rate.
•
t
inertia
is applied to momentum equations only, and it will indirectly affect the
pressure equation because of the slow marching of the velocity change.
•
For convergence solution, the inertial term has no affect on the final solution
because it was eliminated in both LHS and RHS.
Explicit relaxation for momentum and pressure equations. Explicit relaxation
formomentum andpressure equationsisthe regulartechnique, whichisalsocalled
under-relaxation. In this form, the new solution is weighted by the old solution
using the formula:
φ = (1 − α)φ
old
+ αφ
new
(64)
where φ
new
is the updated current solution and φ
old
is the previous value. α is the
relaxationfactor.Ifα =1,allpreviousvaluesareignored,andallupdatedvaluesare
Sunden CH005.tex 10/9/2010 15: 0 Page 203
Equal-order segregated finite-element method 203
used. For convergence difficulties, low value of α will help. Automatic parameter
control was adaptedto eliminate theuser’s input andadjustment effort.Thecontrol
method based on the convergence history identifies the convergence pattern via an
intelligence control theory as proposed by Xie [16].
Solution strategy in fluid flow and heat transfer problems
Because the governing equations are nonlinear, they must be solved iteratively.
Picard successive substitution or Newton–Raphson iterative methods are used. In
Picard successive substitution method, estimates of the solution variables (U, V,
W, P, K, E,T, etc.) are substituted into the governing equations. The equations are
solvedfornewvaluesthatarethenusedastheestimatesforthenextpass.Theglobal
solver controller will either perform a fixed number of these global iterations or
checkfortheconvergencecriterion.Thecontrollerwillstopwheneitherisreached.
Theconvergencecriterionisthelevelatwhichthespecifiedvariable’sresidualnorm
must reach.
MethodA: For weak (one way) coupled flow and heat transfer problem
Loop (nonlinear) until convergence
Solve momentum equations by segregated scheme
Solve pressure equations and velocity correction by segregated scheme
End Loop
Solve energy equations
Method B: For fully coupled flow and heat transfer problem
Loop (nonlinear) until convergence
Solve momentum equations by segregated scheme
Solve pressure equations and velocity correction by segregated scheme
Solve energy equations
End Loop
Method C: For weak (one way) coupled flow and heat transfer problem
Loop (nonlinear) until convergence
Solve coupled continuity, momentum equations
End Loop
Solve energy equations
Method D: For fully coupled flow and heat transfer problem
Loop (nonlinear) until convergence
Solve coupled continuity, momentum equations
Solve energy equations
End Loop
Sunden CH005.tex 10/9/2010 15: 0 Page 204
204 Computational Fluid Dynamics and Heat Transfer
Method E: For fully coupled flow and heat transfer problem
Loop (nonlinear) until convergence
Solve wholly coupled continuity, momentum and energy equations
End Loop
Owing to the importance of the segregated pressure–velocity solution scheme
versuscoupledschemeinmostoftheconvectiveheattransferproblems,themajority
ofcomputationaleffortsinheattransferproblemslieinsolvingtheflowandpressure
fields.
Scenario I demonstrates the decoupled flow and heat transfer, i.e., there are
no temperature-dependent properties in the flow. We can first solve the isothermal
flow (energy equation turned off) to yield a converged flow-field solution and then
solve the energy transport equation alone; this is a naturally segregated approach.
Scenario IIdepicts the coupledflowand heat transfer(typically natural-convection
problem); the common solution practice is to realize the flow–thermal coupling
using thenonlinear iteration untilconvergence.Thisis exactly theadvantage ofthe
velocity–pressure “segregated scheme.”
Figure 5.4 shows the five solution schemes. Methods A and B are segregated
schemes described in this chapter, methods C and D are coupled momentum
and pressure schemes, and method E is completely coupled velocity–pressure–
temperaturesolutionscheme.MethodEisrarelyusedbecauseofthebulkproperties
in global matrix by nature especially for large industry application. In any case, all
methods include solving energy equation.
For multiphysics system of equations to be solved by segregated manner, the
solution scheme is illustrated in Figure 5.4:
1. Velocity–pressure sequential solver (momentum-continuity)
2. Turbulence subiteration
3. Nonlinear iterations
4. Time stepping iteration
At each time step, the nonlinear iterations are performed for the momentum-
continuity, turbulence, front-tracking equation, and temperature. Subiterations of
turbulence transport equations are also used to accelerate the overall convergence
of the iterative process.
Inmoregeneral andwider rangemultiphysicsapplications withstructure, elec-
tromagnetic,acousticsdisciplinescouplingwithfluidflowandheattransfer,typical
application such as fluid–structure interaction (FSI), two-tier approaches may be
applied for the architecture with one unified simulation environment (workbench).
Intierone,CFDpackagetightlyintegratefluidflow,heattransfer,masstransfer,and
acoustics components, and structure package tightly integrate linear dynamics and
explicit dynamics analysis components. In top level tier two, customized interpro-
cesscommunicationtechnologyorgeneralinterprocesscoupling toolslikeMPCCI
Sunden CH005.tex 10/9/2010 15: 0 Page 205
Equal-order segregated finite-element method 205
(I)
(IV)
(II)
Solve u momentum
Solve v momentum
Solve w momentum
Solve pressure and
correct velocity
Solve flow
Initial condition(IC)
Solve turbulence kinetic
energy
Solve eddy dissipation
Update m
t
Solve T (heat transfer)
Solve front tracking
(free surface problem)
While t t
max
(II)
(III)
Figure 5.4. General solution procedure of the segregated solver.
could be used for the communication among different independent packages such
as CFD, structure, electromagnetic analysis.
5.4.3 Data storage and block I/O process
Data storage
The FEM results in a system of linear equations containing a large number of
equations and unknown degree of freedom. A big amount of information requires
a clever data-keeping strategy to avoid whole matrix storage. Knowing that the
sparse properties of global matrix, we are using data storing which is the so-called
skyline or column profile storage. Instead of storing every single matrix element,
we could think of storing only the nonzero hits and corresponding position index
in the matrix assembly phase (Figure 5.5).
However,evenwiththeuseofsuchareducedstoragemode, thematrixisusually
so large that it cannot be stored in core, so that it must be segmented in blocks and
stored on low-speed disk storage.