CHE-614
Computational Fluid Dynamics
Descretisation, Interface Conductivity,
Nonlinearity and Source Term Linearization
Dr. Muhammad Nadeem
[email protected] Room 15, . : 3774
Types of Discretisation
Cell-centered type of grid Vertex-based type of grid
1 2 3 1 2 3 4
• most common type of grid for • most common type of grid for
the finite volume method the finite difference method
• the numerical flux is the focus • the flow variable is the focus
entity in this case entity in this case
2
Location of Control Volume Faces
Practice A: Faces located midways between grids points
Practice B: Grid points placed at the center of control volumes
N N
n
n
W w P e E
s W w P e E
S
s
S
y y
x x
Practice A Practice B
For uniform grid, both practices are identical
The “midway” faces in Practice A do provide greater accuracy in calculating
the heat flux across the face.
Since grid point P is not located at the center of control volume faces in
Practice A, assumption of uniform value is erroneous. 3
Location of Control Volume Faces
Material discontinuity and discontinuous boundary conditions can be
handled using Practice B
4
Location of Control Volume Faces
Half control volume at boundary Control volume face at coincide
with boundary
Practice A Practice B
Practice A produces half-control volume at boundaries while in Practice
B, control volume face lie on the boundary.
Treatment of Half Control Volume
If the boundary temperature is given, no particular difficulty arises, and
no additional equations are required.
If the value of 𝑞𝐵 itself is given, (δx)i
𝑞𝐵
𝑑 𝑑𝑇
𝑘 +𝑆 =0
𝑑𝑥 𝑑𝑥 B i I
𝑑𝑇 𝑑𝑇 Δx
ҧ =0
𝑘 ቤ − 𝑘 ቤ + 𝑆Δ𝑥
𝑑𝑥 𝑖 𝑑𝑥 𝐵
𝑑𝑇
𝑞𝐵 − 𝑞𝑖 + 𝑆𝑢 + 𝑆𝑃 𝑇𝑃 Δ𝑥 = 0 𝑞 = −𝑘
𝑑𝑥
𝑘𝑖 𝑇𝐼 − 𝑇𝐵
𝑞𝐵 + − + 𝑆𝑢 + 𝑆𝑃 𝑇𝐵 Δ𝑥 = 0
𝛿𝑥 𝑖
The required equation for 𝑇𝐵 becomes 𝑎𝐵 𝑇𝐵 = 𝑎𝐼 𝑇𝐼 + 𝑏
where 𝑘
𝑎𝐼 = 𝑎𝐵 = 𝑎𝐼 − 𝑆𝑝 ∆𝑥 𝑏 = 𝑆𝑢 ∆𝑥 + 𝑞𝐵
𝛿𝑥 𝑖
6
Treatment of Half Control Volume
If the heat flux 𝑞𝐵 is specified in terms of a heat transfer coefficient h
and a surrounding-fluid temperature 𝑇𝑓 such that
𝑞𝐵 = −ℎ 𝑇𝑓 − 𝑇𝐵 (δx)i
𝑞𝐵
𝑑 𝑑𝑇
𝑘 +𝑆 =0
𝑑𝑥 𝑑𝑥 B i I
𝑑𝑇 𝑑𝑇 Δx
ҧ =0
𝑘 ቤ − 𝑘 ቤ + 𝑆Δ𝑥
𝑑𝑥 𝑖 𝑑𝑥 𝐵
𝑞𝑖 − 𝑞𝐵 + 𝑆𝑢 + 𝑆𝑃 𝑇𝑃 Δ𝑥 = 0
𝑘𝑖 𝑇𝐼 − 𝑇𝐵
+ ℎ 𝑇𝑓 − 𝑇𝐵 + 𝑆𝑢 + 𝑆𝑃 𝑇𝐵 Δ𝑥 = 0
𝛿𝑥 𝑖
The required equation for 𝑇𝐵 becomes 𝑎𝐵 𝑇𝐵 = 𝑎𝐼 𝑇𝐼 + 𝑏
where 𝑘
𝑎𝐼 = 𝑎𝐵 = 𝑎𝐼 − 𝑆𝑝 ∆𝑥 + ℎ 𝑏 = 𝑆𝑢 ∆𝑥 + ℎ𝑇𝑓
𝛿𝑥 𝑖
7
The Interface Conductivity
So far, we have considered the situation with uniform conductivity k.
In this case, the conductivity at δxe
control-volume face, say e, is P e E
represented by ke.
Significant variations of k are frequently ΔxP ΔxE
encountered x
We shall often know the value of k only at the grid points W, P, E, and so on.
The most straightforward procedure for obtaining the interface diffusion
coefficient ke is to assume a linear variation of k between points P and E
𝑘𝑒 = 𝑓𝑒 𝑘𝑃 + 1 − 𝑓𝑒 𝑘𝐸
where the interpolation factor fe is a ratio defined in terms of the distances
shown in Figure 0.5∆𝑥𝐸 0.5∆𝑥𝑃
𝑓𝑒 = and 1 − 𝑓𝑒 =
𝛿𝑥𝑒 𝛿𝑥𝑒
If the interface e were midway between the grid points, fe would be 0.5, and ke
would be the arithmetic mean of kP and kE
9
The Interface Conductivity …
As long as ke is smoothly varying, this is a perfectly adequate interpolation.
This simple approach cannot accurately handle the abrupt changes of
conductivity that may occur in composite materials.
It is not the local value of conductivity at the interface e that concerns us
primarily, but the main objective is to obtain a good representation for the heat
flux 𝑞𝑒 at the interface via
𝑘𝑒 𝑇𝑃 − 𝑇𝐸
𝑞𝑒 =
𝛿𝑥𝑒
The desired expression for 𝑘𝑒 is the one that leads to a “correct” 𝑞𝑒 .
Let us consider that the control volume 𝑘𝑃 δxe 𝑘𝐸
surrounding the grid point P is filled P e E
with a material of uniform conductivity
𝑘𝑃 , and the one around E with a ΔxP ΔxE
material of conductivity 𝑘𝐸 .
10
The Interface Conductivity …
Surface flux at interface using 𝑘𝑒 𝑇𝑃 − 𝑇𝐸
one-dimensional steady-state 𝑞𝑒 =
𝛿𝑥𝑒
analysis as
𝑘𝑃 δxe 𝑘𝐸
P e E
𝑘𝑃 𝑇𝑃 − 𝑇𝑒 𝑘𝐸 𝑇𝑒 − 𝑇𝐸
𝑞𝑃 = ΔxP ΔxE 𝑞𝐸 =
0.5∆𝑥𝑃 0.5∆𝑥𝐸
0.5∆𝑥𝑃 𝑞𝑃 = 𝑞𝑒 = 𝑞𝐸 0.5∆𝑥𝐸
𝑇𝑒 = 𝑇𝑃 − 𝑞𝑃 𝑇𝑒 = 𝑇𝐸 + 𝑞𝐸
𝑘𝑃 𝑘𝐸
0.5∆𝑥𝑃 0.5∆𝑥𝐸
𝑇𝑃 − 𝑞 = 𝑇𝐸 + 𝑞
𝑘𝑃 𝑘𝐸
0.5∆𝑥𝑃 0.5∆𝑥𝐸
𝑇𝑃 − 𝑇𝐸 = 𝑞 +
𝑘𝑃 𝑘𝐸
𝑇𝑃 − 𝑇𝐸
𝑞=
0.5∆𝑥𝑃 Τ𝑘𝑃 + 0.5∆𝑥𝐸 Τ𝑘𝐸
11
The Interface Conductivity …
For the composite slab between points P and E, a steady one-dimensional
analysis (without sources) leads to
𝑇𝑃 − 𝑇𝐸
𝑞𝑒 =
0.5∆𝑥𝑃 Τ𝑘𝑃 + 0.5∆𝑥𝐸 Τ𝑘𝐸
Thus, an equivalent interface conductivity may be defined as
𝛿𝑥𝑒 0.5∆𝑥𝑃 0.5∆𝑥𝐸
= +
𝑘𝑒 𝑘𝑃 𝑘𝐸
The term 𝛿𝑥𝑒 /𝑘𝑒 may be seen as a resistance to the diffusion transfer between
P and E. −1
1 − 𝑓𝑒 𝑓𝑒
𝑘𝑒 = +
𝑘𝑃 𝑘𝐸
When the interface e is placed midway between P and E, we have 𝑓𝑒 = 0.5;
then 2𝑘𝑃 𝑘𝐸
𝑘𝑒−1 = 0.5 𝑘𝑃−1 + 𝑘𝐸−1 𝑘𝑒 =
𝑘𝑃 + 𝑘𝐸
𝑘𝑒 is the harmonic mean of 𝑘𝑃 and 𝑘𝐸 , rather than the arithmetic mean
12
The Interface Conductivity …
This definition of 𝑘𝑒 leads to following expression for coefficient 𝑎𝐸 .
−1
0.5∆𝑥𝑃 0.5∆𝑥𝐸
𝑎𝐸 = +
𝑘𝑃 𝑘𝐸
The effectiveness of this formulation can be quickly seen in the
following two limiting cases:
1. Let 𝑘𝐸 → 0. Then 𝑘𝑒 = 0
This implies that the heat flux at the face of an insulator becomes zero,
as it should.
𝑘𝐸
2. Let 𝑘𝑃 ≫ 𝑘𝐸 . Then 𝑘𝑒 =
𝑓𝑒
This is to be expected because the high-conductivity material around
point P would offer negligible resistance in comparison with the
material around E
The other implication is that 𝑘𝑒 is not equal to 𝑘𝐸 , but rather 1Τ𝑓𝑒 times it
14
The Interface Conductivity …
Our purpose is to get a correct value of 𝑞𝑒 via
𝑘𝑒 𝑇𝑃 − 𝑇𝐸
𝑞𝑒 =
𝛿𝑥𝑒
The use of 𝑘𝑒 = 𝑘𝐸 Τ𝑓𝑒 yields
𝑘𝐸 𝑇𝑃 − 𝑇𝐸
𝑞𝑒 =
0.5∆𝑥𝐸
When 𝑘𝑃 ≫ 𝑘𝐸 , the temperature 𝑇𝑃 will prevail right up to the interface e,
and the temperature drop 𝑇𝑃 − 𝑇𝐸 will actually take place over the distance
0.5∆𝑥𝐸 . Thus, the correct heat flux will be as given by this equation.
This interface-conductivity formula is based on the steady, no-source, one-
dimensional situation in which the conductivity varies in a stepwise
fashion from one control volume to the next.
Even in situations with nonzero sources or with continuous variation of
conductivity, it performs much better than the arithmetic-mean formula.
15
Nonlinearity
The discretization equation aPTP = aETE + aWTW + Su is a linear
algebraic equation
The set of such equations is solved by the methods for linear algebraic
equations
Nonlinear situations are frequently encountered even in heat
conduction
The conductivity k may depend on T, or the source S may be a
nonlinear function of T.
Then, the coefficients in the discretization equation will
themselves depend on T
Such situations are handled by iteration.
16
Nonlinearity
This process involves the following steps:
1. Start with a guess or estimate for the values of T at all grid points.
2. From these guessed T’s, calculate tentative values of the coefficients in
the discretization equation.
3. Solve the nominally linear set of algebraic equations to get new values of
T.
4. With these T’s as better guesses, return to step 2 and repeat the process
until further repetitions (called iterations) cease to produce any
significant changes in the values of T.
This final unchanging state is called the convergence of the iterations
It is possible that successive iterations would not ever converge to a
solution. The values of T may steadily drift or oscillate with increasing
amplitude.
This process, is called divergence.
17
Source-Term Linearization
When the source S depends on T, we express the dependence in a linear
form given by S = Su + SpTP
This is done because
1) our nominally linear framework would allow only a formally
linear dependence, and
2) the incorporation of linear dependence is better than treating S as a
constant.
When S is a nonlinear function of T, we must linearize it, i.e., specify
the values of Su and Sp, which may themselves depend on T.
During each iteration cycle, Su and Sp would then be recalculated from
the new values of T.
The linearization of S should be a good representation of the S ~ T
relationship.
Further, the basic rule about nonpositive Sp must be obeyed
18
Source-Term Linearization
There are many ways of splitting a given expression for S into Su and
SpTP.
Some of these are illustrated by the following examples
Example 1 Given: S = 5 – 4T.
Su = 5, SP = – 4. This is the most obvious form and is recommended.
𝑆𝑢 = 5 − 4𝑇𝑃∗ , SP = 0.
This is the approach of the lazy person who throws the entire S into Su and
sets SP equal to zero.
𝑆𝑢 = 5 + 7𝑇𝑃∗ , SP = – 11.
This proposes a steeper S ~ T relationship than the one actually given.
The result will be that the convergence of the iterations will slow down.
19
Source-Term Linearization
Example 2 Given: S = 3 + 7T.
Su = 3, SP = 7.
In general this is not acceptable, as it makes SP positive. If the problem could
be solved without iteration, this linearization would give the correct solution,
but if iteration is employed for some reason the presence of a positive SP may
cause divergence
𝑆𝑢 = 3 + 7𝑇𝑃∗ , SP = 0.
This is the practice one should follow when a negative SP is not naturally
forthcoming.
𝑆𝑢 = 3 + 9 𝑇𝑃∗ , SP = – 2.
This is an artificial creation of a negative SP. It will, in general, slow down
the convergence.
20
Source-Term Linearization
Example 3 Given: S = 4 – 5T 3.
𝑆𝑢 = 4 − 5𝑇𝑃∗3 , SP = 0.
This is the lazy-person approach, which fails to take advantage of the known
dependence of S on T.
𝑆𝑢 = 4, SP = −5𝑇𝑃∗2 .
This looks like the correct linearization, but the given S ~ T curve is steeper
than this implies. ∗
𝑑𝑆
𝑆 = 𝑆∗ + 𝑇𝑃 − 𝑇𝑃∗ = 4 − 5𝑇𝑃∗3 − 15𝑇𝑃∗2 𝑇𝑃 − 𝑇𝑃∗
𝑑𝑇
Thus, 𝑆𝑢 = 4 + 10𝑇𝑃∗3 , SP = −15𝑇𝑃∗2 .
This linearization represents the tangent to the S ~ T curve at Tp
𝑆𝑢 = 4 + 20𝑇𝑃∗3 , SP = – 25𝑇𝑃∗2 .
This linearization, which is steeper than the given S ~ T curve, would slow
down convergence. 21
Source-Term Linearization
1 𝑆𝑢 = 4 − 5𝑇𝑃∗3 , SP = 0.
S Given curve
𝑆𝑢 = 4 − 5𝑇 3
2 𝑆𝑢 = 4, SP = −5𝑇𝑃∗2 .
∗3
3 𝑆𝑢 = 4 + 10𝑇𝑃 , SP = −15𝑇𝑃∗2 .
4 𝑆𝑢 = 4 + 20𝑇𝑃∗3 , SP = −25𝑇𝑃∗2 .
𝑇𝑃∗
T
22
Assignment
An exterior wall of a building is a composite wall made of four layers
as shown below:
Consider the layers as homogeneous with following data:
a1 = 2 cm, a2 = 1 cm,
ho hin
a3 = 5 cm, a4 = 1 cm, A B C D
To Tin
a1 a2 a3 a4
kA = 0.72 W/m K, kB = 0.12 W/m K,
kC = 0.026 W/m K, and kD = 0.22 W/m K.
Selecting a suitable discretization scheme, calculate the temperature
distribution within the composite wall and evaluate the heat flow per
unit width based assuming a winter design indoor and outdoor
conditions as: Tin = 25°C, hin = 10 W/m2 °C and
Tout = −25°C, ho = 50 W/m2 C.
23
Assignment
Also, compare these numerical solutions to the analytical solution
Temperature Distribution:
𝑥 − 𝑥𝑛−1 𝑎𝑚
𝑇𝑖𝑛 − 𝑇𝑜 ℎ𝑜 + σ𝑛−1 +1
𝑘𝑛−1 𝑚=1
𝑘𝑚
𝑇𝑛 = 𝑇𝑜 +
𝑎𝑚 ℎ𝑖𝑛
ℎ𝑖𝑛 σ𝑁 + +1
𝑚=1
𝑘𝑚 ℎ𝑜
Heat Flux:
𝑇𝑖𝑛 − 𝑇𝑜
𝑞" =
1 𝑁 𝑎𝑚 1
+ σ𝑚=1 +
ℎ𝑖𝑛 𝑘𝑚 ℎ𝑜
24
