Theme 3 Finite Volume Method (FVM)

MKM 411 – Computational Fluid Dynamics (CFD) Theme 3
Finite Volume Method (FVM)

THEME 3
Overview
→ Basic Methodology
→ Discretization of the Domain
→ Integrate the DE over the CV’s
→ Approximate the Integral
→ Assembling the Discrete System
→ Solve Iteratively or Simultaneously
→ Diffusion and Convective Problems
→ Central Difference, Upwind & QUICK Methods
Use integration
→ Staggered Arrangement
→ Pressure-Velocity Coupled Problems
In FVM We Take Integration from the Equations

→ Therefore, we have to change the format of our five equations (Mass, momentum, and Energy)
to new formats (if needed) so that the terms are easy to integrate.
→ In this way, all different terms of our equations must be like:
∂
(A) + ⋯
∂x, or y or z or t
→ With this idea, let’s get back to our equations and find out which parts must be changed.
∂
We are looking for → (A) +⋯
∂x, or y or z or t
The problem starts with the format of the momentum and energy equations.
For example, if we consider the x-direction of the momentum equation:

2 2 2
∂v ∂vx ∂vx ∂vx ∂p ∂ v ∂ vx ∂ vx Right-hand side terms are
x-direction: ρ ( ∂tx + vx ∂x
+ vy ∂y
+ vz ∂z
) =- ∂x
+ ρfx + μ ( ∂x2x + ∂y2
+ ∂z2
) all the correct form
We can’t take vx out (it is not constant)
∂ ∂vx
∂t
(vx ) ✔ vx ∂x
For Incompressible and Constant Properties

∂vx ∂vy ∂vz
Continuity Equation: ∇∙v=0 → ∂x
+ ∂y
+ ∂z
=0
2
∂ vx ∂
Dv = (v )
Momentum Equation: ρ Dt = - ∇p + ρf + μ∇2 v ∂x 2 ∂x x
2 2 2
∂v ∂vx ∂vx ∂vx ∂p ∂ v ∂ vx ∂ vx
x-direction: ρ ( ∂tx + vx ∂x
+ vy ∂y
+ vz ∂z
) =- ∂x
+ ρfx + μ ( ∂x2x + ∂y2
+ ∂z2
)
2 2 2
∂v ∂vy ∂vy ∂vy ∂p ∂ v ∂ vy ∂ vy
y-direction: ρ ( ∂ty + vx ∂x
+ vy ∂y
+ vz ∂z
) =- ∂y
+ ρfy + μ ( ∂x2y + ∂y2
+ ∂z2
)
2 2 2
∂v ∂vz ∂vz ∂vz ∂p ∂ v ∂ vz ∂ vz
z-direction: ρ ( ∂tz + vx ∂x
+ vy ∂y
+ vz ∂z
) =- ∂z
+ ρfz + μ ( ∂x2z + ∂y2
+ ∂z2
)
1
DT
Energy Equation: ρc = Φ + k∇2 T + ġ
Dt
2 2 2
∂T ∂T ∂T ∂T ∂ T ∂ T ∂ T
ρc ( +u +v + w ) = Φ + k( + + 2 ) + ġ
∂t ∂x ∂y ∂z ∂x2 ∂y2 ∂z
∂u 2 ∂v 2 ∂w 2 ∂v ∂u 2 ∂w ∂v 2 ∂u ∂w 2
Φ = 2μ [( ∂x ) + (∂y) + ( ∂z ) ] + μ [∂x + ∂y
] + μ [ ∂y
+ ∂z
] + μ [ ∂z + ∂x
]
We need to change the format of the momentum and energy equations

2 2 2
∂u ∂u ∂u ∂u ∂p ∂ u ∂ u ∂ u
ρ ( ∂t + u ∂x + v ∂y + w ∂z ) = - ∂x
+ ρfx + μ ( ∂x2 + ∂y2
+ ∂z2
) u = ui⃗ + vj⃗ + wk
⃗⃗
∂u ∂u ∂u
Let’s see if we can write: div(uu) = u ∂x + v ∂y + w ∂z
⃗⃗) = ( ∂ ⃗i +
div(uu) = ∇ ∙ (u2⃗i + uvj⃗ + uwk
∂ ⃗ ∂
j + ∂z ⃗⃗k) ∙ (u2⃗i + uvj⃗ + uwk ⃗⃗) = ∂ (u2 ) + ∂ (uv) + ∂ (uw)
∂x ∂y ∂x ∂y ∂z
∂u ∂u ∂v ∂w ∂u ∂u ∂v ∂w ∂u ∂u ∂u
= 2u ∂x + (v ∂y + u ∂y) + (u ∂z + w ∂z ) = (u ∂x + u ∂y + u ∂z ) + (u ∂x + v ∂y + w ∂z )
∂u ∂v ∂w ∂u ∂u ∂u ∂u ∂u ∂u
= u ( ∂x + + ∂z ) + (u ∂x + v ∂y + w ∂z ) = (u ∂x + v ∂y + w ∂z ) ✔
⏟ ∂y
Continuity
for incompressible
=0
Du ∂u ∂u ∂u ∂u ∂u The equation will be given

Therefore, we can write: ρ = ρ( + div(uu)) = ρ ( +u +v +w ) on the formula sheet
Dt ∂t ∂t ∂x ∂y ∂z
For Constant Properties We usually keep the div/∇ in the equation for this section of work
Du ∂p
x- momentum in the CFD textbook: ρ Dt = - ∂x
+ μ div(grad u) + SMx
Du ∂p
x-momentum as already determined: ρ =- + μ∇2 u + ρfx Laplacian derivative: ∇2 u
Dt ∂x
2 2 2
∂ ∂ ⃗ ∂ ⃗⃗ ∂u ∂u ⃗ ∂u ⃗⃗ ∂ u ∂ u ∂ u
div(grad u) = ∇ ∙ (∇u) = ∇2 u = (∂x ⃗i + ∂y
j + ∂z k ) ∙ ( ∂x ⃗i + ∂y
j + ∂z k ) = ∂x2
+ ∂y2
+ ∂z2
= ∇2 u ✔
SMx is the source term (e.g., body force):

Source term: SMx = ρfx S – source term; M – equation (momentum in this instance); x – direction (x-direction in this instance)
Governing Equations in the CFD Textbook
In the textbook, the velocity vector is indicated as u: u = ui⃗ + vj⃗ + wk

⃗⃗
∂ρ ∂(ρu) ∂(ρv) ∂(ρw) ∂ρ
Mass conservation: ∂t
+ ∂x
+ ∂y
+ ∂z
=0 → ∂t
+ div(ρu) = 0
For FVM, we only work with incompressible, constant property, and Newtonian fluids.
Continuity: div(u) = 0 ∂u ∂v ∂w
∂x
+ ∂y + ∂z = 0 ✔
⃗ ⃗ ⃗⃗
u = ui + vj + wk
For FVM, we only work with incompressible, constant property, and Newtonian fluids.
Therefore, the Navier-Stokes equation is valid for momentum Source term
Du ∂p
x- momentum in the CFD textbook: ρ Dt = - ∂x
We already had:
2 2 2
∂u ∂u ∂u ∂u ∂p ∂ u ∂ u ∂ u Du ∂p
ρ ( ∂t + u ∂x + v ∂y + w ∂z ) = - ∂x
+ ρfx + μ ( ∂x2 + ∂y2
+ ∂z2
) ← ρ Dt = - ∂x
+ μ∇2 u + ρfx
2
For Constant Properties We usually keep the div/∇ in the equation for this section of work
Du ∂p
x- momentum in the CFD textbook: ρ =- + μ div(grad u) + SMx
Dt ∂x
Du ∂p
x-momentum as already determined: ρ Dt = - ∂x
+ μ∇2 u + ρfx Laplacian derivative: ∇2 u
2 2 2
∂ ∂ ⃗ ∂ ∂u ∂u ⃗ ∂u ∂ u ∂ u ∂ u
div(grad u) = ∇ ∙ (∇u) = ∇2 u = (∂x ⃗i + ∂y
j + ∂z ⃗⃗k) ∙ ( ∂x ⃗i + ∂y
j + ∂z ⃗⃗k) = ∂x2
+ ∂y2
+ ∂z2
= ∇2 u ✔
SMx is the source term (e.g., body force):

Source term: SMx = ρfx S – source term; M – equation (momentum in this instance); x – direction (x-direction in this instance)
Du ∂p
ρ Dt
=- ∂x
Du ∂t
We can write: = + div(uu)
Dt ∂x
⃗⃗) = ( ∂ ⃗i +
div(uu) = ∇ ∙ (u2⃗i + uvj⃗ + uwk
∂ ⃗ ∂
j + ⃗⃗k) ∙ (u2⃗i + uvj⃗ + uwk⃗⃗) = ∂ (u2 ) + ∂ (uv) + ∂ (uw)
∂x ∂y ∂z ∂x ∂y ∂z
∂u ∂u ∂v ∂w ∂u ∂u ∂v ∂w ∂u ∂u ∂u
= 2u ∂x + (v ∂y + u ∂y) + (u ∂z + w ∂z ) = (u ∂x + u ∂y + u ∂z ) + (u ∂x + v ∂y + w ∂z )
∂u ∂v ∂w ∂u ∂u ∂u ∂u ∂u ∂u
= u ( ∂x + + ) + (u ∂x + v ∂y + w ∂z ) = (u ∂x + v ∂y + w ∂z ) ✔
⏟ ∂y ∂z
Continuity
for incompressible
=0
Du ∂u ∂u ∂u ∂u ∂u The equation will be given

Therefore, we can write: ρ Dt = ρ ( ∂t + div(uu)) = ρ ( ∂t + u ∂x + v ∂y + w ∂z ) on the formula sheet
For Incompressible, Constant Property, and Newtonian Fluids
Du ∂p ∂u ∂p
x-momentum: ρ Dt = - ∂x
+ μ div(grad u) + SMx or ρ ( ∂t + div(uu)) = - ∂x
x-momentum as already determined:
2 2 2
Du ∂p ∂u ∂u ∂u ∂u ∂p ∂ u ∂ u ∂ u
ρ =- + ρfx + μ∇2 u → ρ( +u +v +w ) =- + ρfx + μ ( + + )
Dt ∂x ∂t ∂x ∂y ∂z ∂x ∂x2 ∂y2 ∂z2
Dv ∂p ∂v ∂p
y-momentum: ρ Dt = - ∂y
+ μ div(grad v) + SMy or ρ ( ∂t + div(vu)) = - ∂y
+ μ div(grad v) + SMy
y-momentum as already determined:
2 2 2
Dv ∂p ∂v ∂v ∂v ∂v ∂p ∂ v ∂ v ∂ v
ρ Dt = - ∂y
+ ρfy + μ∇2 v → ρ ( ∂t + u ∂x + v ∂y + w ∂z ) = - ∂y
+ ρfy + μ (∂x2 + ∂y2
+ ∂z2
)
Dw ∂p ∂w ∂p
z-momentum: ρ =- + μ div(grad w) + SMz or ρ ( + div(wu)) = - + μ div(grad w) + SMz
Dt ∂z ∂t ∂z
z-momentum as already determined:

2 2 2
Dw ∂p ∂w ∂w ∂w ∂w ∂p ∂ w ∂ w ∂ w
ρ Dt
=- ∂z
+ ρfz + μ∇2 w → ρ ( ∂t + u ∂x + v ∂y + w ∂z ) = - ∂z
+ ρfz + μ ( ∂x2 + ∂y2
+ ∂z2
)
We only consider steady cases in FVM in MKM 411, therefore, for steady, incompressible, constant
property, and Newtonian fluids:
x-momentum:
∂p
New form for FVM: ρ div(uu) = - ∂x
2 2 2
∂u ∂u ∂u ∂p ∂ u ∂ u ∂ u
Previous form: ρ (u ∂x + v ∂y + w ∂z ) = - ∂x
+ ρfx + μ ( ∂x2 + ∂y2
+ ∂z2
)
3
y-momentum:
∂p
New form for FVM: ρ div(vu) = - ∂y
2 2 2
∂v ∂v ∂v ∂p ∂ v ∂ v ∂ v
Previous form: ρ (u ∂x + v ∂y + w ∂z ) = - ∂y
+ ρfy + μ (∂x2 + ∂y2
+ ∂z2
)
z-momentum:
∂p
New form for FVM: ρ div(wu) = - + μ div(grad w) + SMz
∂z
2 2 2
∂w ∂w ∂w ∂p ∂ w ∂ w ∂ w
Previous form: ρ (u ∂x + v ∂y + w ∂z
) =- ∂z
+ ρfz + μ ( ∂x2 + ∂y2
+ ∂z2
)
Conservation of Energy (in the CFD Textbook)

Di
ρ = - p div u(k grad T) + Φ + Si Si is the source term for heat generation
Dt
For incompressible, constant property, and Newtonian fluids i = cT

DT
ρc Dt = k div(grad T) + Φ + Si Laplacian of Temperature: div(grad T) = ∇ ∙ (∆T) = ∇2 T
DT
We already had: ρc Dt = Φ + k∇2 T + ġ
DT ∂T
ρc Dt = k div(grad T) + Φ + Si or ρc (∂t + div(Tu)) = k div(grad T) + Φ + Si
⃗⃗) = ( ∂ ⃗i +
div(Tu) = ∇ ∙ (Tui⃗ + Tvj⃗ + Twk
∂ ⃗ ∂ ⃗⃗
j + ∂z k ) ∙ (Tui⃗ + Tvj⃗ + Twk⃗⃗) = ∂ (Tu) + ∂ (Tv) +
∂
(Tw)
∂x ∂y ∂x ∂y ∂z
∂T ∂u ∂T ∂v ∂T ∂w ∂u ∂v ∂w ∂T ∂T ∂T
= (u +T ) + (v + T ) + (w + T ) = (T + T + T ) + (u + v + w )
∂x ∂x ∂y ∂y ∂z ∂z ∂x ∂y ∂z ∂x ∂y ∂z
∂u ∂v ∂w ∂T ∂T ∂T ∂T ∂T ∂T
= T( + + ) + (u + v + w ) = (u + v + w ) ✔
⏟∂x ∂y ∂z ∂x ∂y ∂z ∂x ∂y ∂z
Continuity
for incompressible
=0
∂T
ρc (∂t + div(Tu)) = k div(grad T) + Φ + Si
We only consider steady cases in FVM in MKM 411, therefore, for steady, incompressible, constant
ρc(div(Tu)) = k div(grad T) + Φ + Si
For solids: k div(grad T) + Si = 0
Therefore, for steady, incompressible, constant property, and Newtonian fluids:
Continuity: div(u) = 0 ρfx

∂p
x-momentum: ρ div(uu) = - ∂x
+ μ div(grad u) + SMx ρfy
∂p
y-momentum: ρ div(vu) = - + μ div(grad v) + SMy ρfz These equations will be on the formula sheet
∂y
z-momentum: ρ div(wu) = -
∂p
+ μ div(grad w) + SMz ė g
∂z
Energy Equation: ρc(div(Tu)) = k div(grad T) + Φ + Si
4
Finite Volume Method

What is the Finite Volume Method in simple words? What are we looking for in the Finite Volume
Method?
General Methodology: Finite Volume = Basic Methodology
In general, the FVM involves the following steps:
1. Divide the Domain into Control Volumes (CVs)

The union of al CVs should cover the whole
problem domain.
→ Each node in a control
volume will be at the centre
This is different to the Finite Difference Method where the nodes
of the control volume were on the corners of the divisions
Dimensions Description
One
Dimensional
→ We use capital letters for nodes North

→ We usually call the interface between two small
Control Volumes a face (surface)
→ We usually indicate the faces with lower case West East
letters, e.g., the east face of node P is e, the west
face of point P is w
→ Therefore, we will indicate the information of the South
nodes by capital letters and the information of
Two the faces (interface between two CVs) by lower ∆x and ∆y are not going to change
Dimensional case letters between nodes for MKM 411.
→ Notation:
If the notation is for a face, the
→ TP : temperature at node P value will need to be used to find
→ uP : velocity in the x-direction at node P the parameter at the node.
→ vP : velocity in the y-direction at node P
→ ue : velocity in the x-direction at the east face
→ uw : velocity in the x-direction at the west face
→ vn : velocity in the y-direction at the north face
→ vs : velocity in the y-direction at the south face
2. Finding the Governing Differential Equations and Simplifying them for the Problem (CVs)
For example: One-dimensional steady heat conduction without heat generation
∂T
General equation: ρc (∂t + div(Tu)) = div(k grad T) + Φ + Si
Simplifying the equation for steady one-dimensional heat conduction without heat generation:
d dT
ρc(0 + 0) = ∇ ∙ (k ∇ T) + 0 + 0 → dx
(k dx) =0 ∇ →
d
dx
and one-dimensional ∴ k∇T → k
dT
dx
5
3. Integrate the Differential Equation Over the CVs

Over the CV → ΔV
d dT
∫ (k ) dV = 0
∆V dx dx
Each CV gives one equation and finally, the final number of

equations much be equal to the number of the CVs
4. Applying the Divergence Theorem for Changing the Integral Over the Volume of the Surface
∫V ∇ ∙ Φ dV = ∮A n
⃗⃗ Φ dA (Divergence theorem) V → A
d dT
∫∆V dx (k dx) dV = 0 We usually omit n
⃗⃗ because we use rectangular shapes in MKM 411
dT
Φ = (k dx)
d dT dT
∫∆V dx (k dx) dV = ∮A (k dx) ⃗n⃗ dA = 0
That’s why we use these forms; for steady, incompressible, constant

ρfx
Continuity: div(u) = 0
∂p ρfy
These equations will be on the formula sheet.
∂p
y-momentum: ρ div(vu) = - + μ div(grad v) + SMy ρfz These equations used the divergence theorem
∂y for the surface integral.
∂p ė g
z-momentum: ρ div(wu) = - + μ div(grad w) + SMz
∂z
5. Approximation of the Surface Integrals

dT
∮A (k dx) n
⃗⃗ dA = 0
dT dT dT
∮A (k dx) n
⃗⃗ dA = (kA )
dx
- (kA )
dx w
=0 (face after – face before)
e
d dT dT dT
∫ΔV dx (k dx) dV = (kA dx) - (kA dx) = 0
e w
6. Approximation of Function Values and Derivatives by Interpolation with Nodal Values

dT TE - TP T -T δx : δx – for different Δx’s, PE – between P and E → usually leave as Δx
(kA ) = ke Ae = k Ae E P k isPEconstant ∴ k → k
dx e δxPE Δx e
dT TP - TW TP - TW
(kA dx) = kw Aw δx = k Aw Δx
w WP
7. Assembling and Solution of the Discrete Algebraic System

dT dT TE - TP TP - TW
(kA ) - (kA ) =0 → (kAe )- (kAw )=0
dx e dx w Δx Δx
kA kA kAw kA These equations apply for all
→ ( Δxe - Δxw ) TP = T + Δxe TE
Δx W
internal nodes (boundary nodes
slightly different)
6
→ The result is a set of linear algebraic equations

→ There is one equation for each CV
For each control volume: ap TP = aWTW + aE TE
→ Solve iteratively or simultaneously
aP TP = aW TW + aE TE
[A][T] = [C] Matrix form
Gaussian elimination or [T] = [A]-1 [C]
Finite Volume Method (Fast Review Up to Now)

1. Divide the domain into control volumes
2. Finding (or simplifying) the differential equation(s) for the CVs

DT 1D, conduction, without heat generation d dT
ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) =0
DT 1D, conduction, with heat generation d dT
ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) + ė g = 0
3. Integrate the differential equations(s) over the CVs
d dT
1D, conduction, without heat generation: ∫∆V dx (k dx) dV = 0
d dT
1D, conduction, with heat generation: ∫∆V dx (k dx) dV + ∫V ė g dV = 0
4. Using mathematical relations and applying approximation (all internal nodes)
Don’t need to write the divergence theorem in the exam
∫V ∇ ∙ Φ dV = ∮A n
⃗⃗ Φ dA (Divergence theorem) d
→ → A (final form)
dx
dT dT dT T -T dT TP - TW
1D, conduction, without heat generation: (kA dx) - (kA dx) = 0, | = E∆x P, dx|
dx e
= ∆x
e w w
dT dT ̅̇ g V = dT T -T dT TP - TW
1D, conduction, with heat generation: (kA ) - (kA ) + e 0, | = E P, | =
dx e dx w dx e ∆x dx w ∆x
The above equations are the general equations
5. Finding the discretised equation(s) for each CV
1D, conduction, without heat generation: aP TP = aW TW + aE TE
1D, conduction, with heat generation: aP TP = aW TW + aE TE + Su
6. Assembling and solution of a discrete system of equations
[A][T] = [C]
Example 4.1
Source-free → no heat generation
Consider the problem of source-free heat
conduction in an insulated rod whose ends are
maintained at constant temperatures of 100°C and
500°C respectively. The one-dimensional problem is
sketched in the Figure. Calculate the steady state
temperature distribution in the rod. The thermal conductivity, k, equals 1000 W/m/K, the cross-
-3
sectional area, A, is 10 × 10 m2 .
Solution
d dT
Energy Equation: ρc(div(Tu)) = k div(grad T) + Φ + Si → (k ) =0 No convection ∴ no div(Tu)
dx dx
Let us divide the length of the rod into five equal control volumes
There will be 5 nodes, 5 CVs, 5 equations, and 5 unknowns
Nodes 1 and 5 are boundary nodes; boundary nodes have at least one less node than internal nodes
7
d dT
Simplified equation: dx
(k dx) =0
d dT dT dT
General equation that works for all nodes: ∫ΔV dx (k dx) dV = (kA dx) - (kA dx) =0
e w
For internal nodes:

The grid consists of five nodes. For each one of nodes 2, 3, and 4, the temperature values to the
east and west are available as nodal values.
d dT dT dT
∫ΔV dx (k dx) dV = (kA dx) - (kA dx) =0
e w
dT dT dT TE - TP dT TP - TW
(kA ) - (kA ) =0 → (kA ) = ke Ae and (kA ) = kw Aw
dx e dx w dx e δxPE dx w δxWP
k kw k k
After rearranging: (δxe Ae + A ) TP
δxWP w
= (δxw Aw ) TW + (δxe Ae ) TE
PE WP PE
The thermal conductivity (ke = kw = k), node spacing (δx), and cross-sectional area (Ae = Aw = A)
are constants. Therefore, the discretised equation for nodal points 2, 3, and 4 is aP TP = aWTW + aE TE
aW aE aP
k k
A A aW + aE
δx δx
Nodes 1 and 5 are boundary nodes
δx/2 at boundary
dT dT
(kA dx) - (kA dx) =0
e w
TE - TP TP - TA k 2k k 2k
→ kA ( ) - kA ( ) =0 → ( A+ A) TP = 0 ∙ TW + ( A) TE + ( A) TA TA is known
δx e δx/2 δx δx δx δx
w
Source term
Equation for boundary node 1: aP TP = aW TW + aE TE + Su
aW aE aP Su
kA 2kA 2kA
0 aW + aE + T
δx δx δx A
For boundary node 5:
δx/2 at boundary
8
dT dT
(kA dx) - (kA dx) =0
e w
TB - TP TP - TW k 2k k 2k
→ kA ( ) - kA ( ) =0 → ( A+ A) TP = ( A) TW + 0 ∙ TE + ( A) TB TB is known
δx/2 δx w δx δx δx δx
e
Equation for boundary node 5: aP TP = aW TW + aE TE + Su Source term
aW aE aP Su
kA 2kA 2kA
0 aW + aE + T
δx δx δx B
Node aW aE Su aP
1 0 100 200TA 300
2 100 100 0 200
3 100 100 0 200
4 100 100 0 200
5 100 0 200TB 300
The resulting set of algebraic equations for this example is:

300T1 = 100T2 + 200TA 300 -100 0 0 0 T1 200TA TA = 100
200T2 = 100T1 + 100T3 -100 200 -100 0 0 T2 0
Matrix
200T3 = 100T2 + 100T4 → 0 -100 200 -100 0 T3 = 0
200T4 = 100T3 + 100T5 0 0 -100 200 -100 T4 0
300T5 = 100T4 + 200TB [ 0 0 0 -100 300] [T5 ] [ 200TB ] TB = 500
T1 140 If more temperature values are required, the problem needs to be divided into more
T2 220 control volumes.
Numerical solution: T3 = 300 A smaller Δx is usually better (more accurate), but the computation time will be more.
The numerical solution has now been found; the analytical solution needs to be found
T4 380 next.
[T5 ] [460]
k can be taken out because it is a constant
d dT d dT dT
Analytical solution: dx
(k dx) =0 → ( )
dx dx
=0 → (dx) = C1 → T(x) = C1 x + C 2
T(0) = 100: T(0) = C1 (0) + C2 = 100 ∴ C2 = 100
T(0.5) = 500: T(0.5) = C1 (0.5) + 100 = 500 ∴ C1 = 800
∴ Exact solution: T(x) = 800x + 100
9


ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) =0
ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) + ė g = 0
d dT
d dT
∫V ∇ ∙ Φ dV = ∮A ⃗n⃗ Φ dA (Divergence theorem) → → A (final form)
d
dx
dx e
= ∆x
e w w
dT dT ̅ dT TE - TP dT TP - TW
1D, conduction, with heat generation: (kA ) - (kA ) + ė g V = 0, | = ∆x , dx| =
dx e dx w dx e w ∆x
[A][T] = [C]
Central Difference Scheme

No temperature value is given at face e → there is no node at the face ∴ we
should use the average between node E and node P
TE + TP
Ae Te = ? → Ae Te = ∆y ( )
2
The Central Difference Scheme applies in FVM because the assumption is that the
CV is small enough that the centre of the cell is equal to the average of the cell.
Do not mix Ae Te with:

Node after - before in the direction
dT dT T⏞
E - TP
Ae (dx) =? → Ae (dx) = ∆y ( ∆x
)
e e
d dT d dT dT dT
(k ) =0 → ∫ΔV dx (k dx) dV = (kA dx) - (kA dx) =0
dx dx e w
dT TE - TP dT TP - TW
If A is constant: (kA dx) = ke A Δx
and (kA dx) = kw A Δx
e w
If k is constant: kw = ke = k
If k is changing with temperature: use the Central Difference Scheme (average from nodal values)
kP + kW kP + kE
→ kw = 2
and ke = 2
10
Let’s Review the Regulation on Faces (up to now)

Aw Tw = ? → Aw = ∆y × ⏟t = ∆y and Tw = ?
1
The thickness (t) is usually 1
We do not have any nodes at the faces; therefore, we must

decide on a scheme for these faces.
The Central Difference Scheme gives us the average value
of the nodes that the face is in between
TW + TP The Central Difference Scheme for an unknown at a face is the average value of
Tw =
2 the nodes that the face is in between
Central Difference Scheme

The Central Difference Scheme can be applied to any unknown Γ (gamma → k)
and any unknown ϕ (phi → T); these equations can also be applied to heat and
mass transfer
dϕ dϕ dϕ ϕE - ϕP
(ΓA dx ) = Γe Ae ( dx ) → (ΓA dx ) = Γe Ae ( δxPE
)
e e e
dϕ dϕ dϕ ϕP - ϕW
(ΓA dx ) = Γw Aw ( dx ) → (ΓA dx ) = Γw Aw (δxWP
)
w w w
ΓP + ΓE ΓW + ΓP
Γe = 2
and Γw = 2
← This part will be different for an Upwind Scheme
An Tn = ? → An = ∆x × ⏟t = ∆x and Tn = ?
1
The thickness (t) is usually 1 4 internal nodes = a length of 4Δx

The Central Difference Scheme gives us the average value
of the nodes that the face is in between
TN + TP The Central Difference Scheme for an unknown at a face is the average value of
Tn =
Integrate the Differential Equation (s) over the CVs
∫ΔV D.E. dV Incorrect integration impacts the final answer greatly; there are four different cases of differential equations
Differential equation(s) can include:

2
∂ ϕ
1. ← second order DE
∂F2
∂ϕ
2. ∂F
← first order DE
We usually don’t use third order DEs
3. ϕ ← unknown variable
4. C ← constant
ϕ can be T, u, v, w. P
F can be x, y, z
11
Examples of different D.E.s:

2
∂ ϕ d dϕ d dT
1. ⇒ ( ) ⇒ like ( )
∂η2 dF dF dx dx
d dT dT dT
∫ΔV dx (k dx) dV = (kA dx) - (kA dx)
⏟ e w
First order approximation
dT TE - TP dT T -T
(dx) = Δx
and (dx) = PΔx W
e w
∂ϕ d d
2. ∂F
⇒ dx
(ϕ) ⇒ like dx (u) Take ϕ out
d
∫ΔV dx (u) dV = (Au)e - (Au)w
d
If ϕ = uT and F = x → ∫ΔV dx (uT) dV = (AuT)e - (AuT)w = Ae ue Te - Aw uw Tw
∂ϕ uE + uP T +T uW + uP T +T
Central Difference Scheme → = Δy ( ) ( E P) - Δy ( ) ( W P)
∂F 2 2 2 2
3. ϕ
∫ΔV T dV = TP ΔV The Central Difference Scheme applies in FVM because the assumption is that the CV is small enough
that the centre of the cell is equal to the average of the cell.
∫ΔV u dV = uP ΔV
∫ΔV v dV = vP ΔV
For any unknowns at the face, we can choose one of the following:
1. Central Differencing Scheme (CDS)
2. Upwind Differencing Scheme (UDS) Will be told in a test/exam which method to use
3. Flux-Blending Technique (FBT)
4. QUICK Scheme
An Example with a Source Term

d dT
∇ ∙ (k∇T) + ė g = 0 or for 1D dx
(k dx) + ė g = 0
ė g = 200T + 1000 W/m3

nd Unknown
2 order
variable Constant
⏞
d dT
dx
(k dx) + ⏞
200T + ⏞
1000 = 0
d dϕ dϕ dϕ
∫∆V dx (Γ dx ) dV + ∫∆V S dV = (ΓA dx ) - (ΓA dx ) + S̅ ΔV = 0
e w
dT dT
(kA dx) - (kA dx) + 200TP ΔV + 1000ΔV = 0
e w
The general transport equation: div(Γ grad ϕ) + Sϕ = 0
d dT
Like: ∇ ∙ (k∇T) + ė g = 0 or for 1D dx
(k dx) + ė g = 0
The control volume integration:

The first term is a second order term;
∫CV div(Γ grad ϕ) dV + ∫CV Sϕ dV = ∫A n ∙ (Γ grad ϕ) dA + ∫CV Sϕ dV = 0 the second term is a first order term
12
The following steps apply for all the nodes
Step 1: Grid Generation
In general: δxWP ≠ δxPE ≠ δxwe
In this course: δxWP = δxPE = δxwe = Δx Δx is fixed
Step 2: Discretisation
(Integration of the governing equation to yield a
discretised equation)
div(Γ grad ϕ) + Sϕ = 0
∫ div(Γ grad ϕ) dV + ∫ Sϕ dV = ∫ n ∙ (Γ grad ϕ) dA + ∫ Sϕ dV = 0

CV CV A CV
↓
1D
d dϕ dϕ dϕ
∫ (Γ ) dV + ∫ S dV = (ΓA ) - (ΓA ) + S̅ ΔV = 0
∆V dx dx ∆V dx e dx w
(ΓA ) = Γe Ae ( ) → (ΓA ) = Γe Ae ( )
dx e dx e dx e δxPE
(ΓA dx ) = Γw Aw ( dx ) → (ΓA dx ) = Γw Aw ( δxWP
)
w w w
Central Difference Method

⏞ ΓP + ΓE ΓW + ΓP
Γe = and Γw = ← This part will be different for an Upwind Scheme
2 2
dϕ dϕ ϕE - ϕP ϕP - ϕW
(ΓA dx ) - (ΓA dx ) + S̅ ΔV = 0 → Γe Ae ( δxPE
) - Γw Aw (δxWP
) + S̅ ΔV = 0
e w
In general: the source term S may be a function of the dependent variable → S̅ ΔV = Su + Sp ϕP

ϕE - ϕP ϕP - ϕW
→ Γe Ae ( δxPE
) - Γw Aw ( δxWP
) + (Su + Sp ϕP ) = 0
Γ Γw Γ Γ
This can be rearranged as: (δxe Ae + A
δxWP w
- Sp ) ϕP = (δxw Aw ) ϕW + (δxe Ae ) ϕE + Su
PE WP PE
aP ϕP = aW ϕW + aE ϕE + Su
aW aE aP
Γw Γe
A A aW + aE - Sp
δxWP w δxPE e


ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) =0
ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) + ė g = 0
13

d dT
d dT
∫V ∇ ∙ Φ dV = ∮A n
⃗⃗ Φ dA (Divergence theorem) → → A (final form)
d
dx
dx e
= ∆x
e w w
[A][T] = [C]
General Form: Diffusion Problems

Consider the steady state diffusion of a property ϕ in a one-dimensional
domain.
d dϕ
dx
(Γ dx ) +S=0
Where Γ is the diffusion coefficient and S is the source term.

d dϕ dϕ dϕ
∫∆V dx (Γ dx ) dV + ∫∆V S dV = (ΓA dx ) - (ΓA dx ) + S̅ ΔV = 0
e w
In a uniform grid, linearly interpolated values for Γe and Γw are given (using the Central Difference
Scheme) as
ΓW + ΓP ΓP + ΓE
Γw = 2
and Γe = 2
dϕ ϕP - ϕW dϕ ϕE - ϕP
(ΓA ) = Γw Aw ( ) and (ΓA ) = Γe Ae ( ) Remember @ the boundary one side has a length of δx/2
dx w δxWP dx e δxPE
The source term S may be a function of the dependent variable: S̅ ΔV = Su + Sp ϕP

ϕE - ϕP ϕP - ϕW
Γe Ae ( δxPE
) - Γw Aw ( δxWP
) + (Su + Sp ϕP ) = 0
The table is used for if the solution needs to be

aP ϕP = aW ϕW + aE ϕE + Su
coded
aW aE aP
Γw Γe
A A aW + aE - Sp
δxWP w δxPE e
dϕ dϕ
(ΓA ) - (ΓA ) + S̅ ΔV = 0
dx e dx w
(ΓA dx ) = Γe Ae ( dx ) → (ΓA dx ) = Γe Ae ( δxPE
)
e e e
(ΓA dx ) = Γw Aw ( dx ) → (ΓA dx ) = Γw Aw ( δxWP
)
w w w
Central Difference Method

⏞ ΓP + ΓE ΓW + ΓP
Γe = 2
and Γw = 2
← This part will be different for an Upwind Scheme
14
Example 4.2
The figure alongside shows a large plate of thickness L = 2 cm with
constant thermal conductivity k = 0.5 W/mK and uniform heat
generation q = 1000 kW/m3 . (q → ė g ) ė g is written as q in the textbook
This is a large plate with a small thickness value ∴ the 1D assumption is valid
Assuming that the dimensions in the y- and z-directions are so large that
temperature gradients are significant in the x-direction only, calculate
the steady state temperature distribution. Compare the numerical result
with the analytical solution.
Therefore: Solid, Steady, and 1D
Solution
Steady: (Solid, Steady, and 1D)
d dT d dT
Energy Equation: ρc(div(Tu)) = k div(grad T) + Φ + Si → (k ) + ė g = 0 → (k ) +q=0
dx dx dx dx
The domain is divided into five control

volumes
There will be 5 nodes, 5 CVs, 5 equations, and 5
unknowns; A and B have a constant temperature
boundary condition
d dT
∫∆V dx (k dx) dV + ∫∆V q dV = 0
dT dT
The general equation: [(kA dx) - (kA dx) ] + q∆V = 0
e w
The general equation is obtained exactly after integration
Therefore, the discretised equation for nodal points 2, 3, and 4 is (2, 3, and 4 have nodes before and after)
TE - TP T -T ke A kw A ke A kw A
[ke A ( ) - kw A ( P W )] + qAδx = 0 → ( + ) TP = ( ) TW + ( ) TE + qAδx
δx δx δx δx δx δx
A will get cancelled out
Equation for nodes 2, 3, and 4: aP TP = aW TW + aE TE + Su The equation will apply even if there are 1000 nodes
aW aE aP Sp Su
kA kA
aW + aE - Sp 0 qAδx
δx δx
Nodes 1 and 5 are boundary nodes

dT dT
[(kA dx) - (kA dx) ] + q∆V = 0
e w
TE - TP TP - TA
[ke A ( δx
) - kA A ( )] + qAδx = 0
δx/2
aW aE aP Sp Su
kA 2kA 2kA
0 aW + aE - Sp - qAδx + T
δx δx δx A
There is more than one source term for node 1 because the value for TA is known
15

dT dT
[(kA dx) - (kA dx) ] + q∆V = 0
e w
TB - TP TP - TW
[kB A ( ) - kw A ( δx
)] + qAδx = 0
δx/2
aW aE aP Sp Su
kA 2kA 2kA
0 aW + aE - Sp - qAδx + T
δx δx δx B
There is more than one source term for node 5 because the value for TB is known
There are now 5 equations and 5 unknowns
Node aW aE Su Sp aP = aW + aE - Sp
1 0 125 4000 + 250TA - 250 375
2 125 125 4000 0 250
3 125 125 4000 0 250
4 125 125 4000 0 250
5 125 0 4000 + 250TB - 250 375
[A][T] = [C]
375 -125 0 0 0 T1 29000 T1 150

-125 250 -125 0 0 T2 4000 T2 218
0 -125 250 -125 0 T3 = 4000 → T3 = 254
0 0 -125 250 -125 T4 4000 T4 258
[ 0 0 0 -125 375] [T5 ] [54000] [T5 ] [230]
The numerical solution has now been
determined
The analytical solution: T=

T -T q
[BL A + 2k
(L - x)] x + TA
Usually have a function plot for the

exact/analytical solution and points for the
numerical solution
Example 4.3
The cylindrical fin (alongside) has a uniform cross-sectional
area A. Calculate the temperature distribution along the fin.
heat conduction
convection
⏞
d dT
Heat transfer equation: dx
(kA dx) - ⏞
hP(T - T∞ ) = 0
P is the fin perimeter
As seen in the equation above, there is no heat generation
T - T∞ cosh[n(L - x)]
Analytical solution: TB - T∞
= cosh(nL)
A fin is added to a surface to improve heat transfer; an insulated endpoint is not practical in real life application but is used in this example
for simplicity purposes.
16
Solution
Now, FVM:
d dT
dx
(kA dx) - hP(T - T∞ ) = 0
kA is constant (divide by kA on both sides)

d dT hP
( )
dx dx
- n2 (T - T∞ ) = 0 n2 = kA Known temperature:
T∞ = 20°C (Ambient temperature)
Divide the length into five control volumes
If there are 5 CVs, there will be 5 nodes; nodes 1 and 5 are boundary
nodes; nodes 2, 3, and 4 are internal nodes
d dT
∫∆V (dx) dV - ∫∆V n2 (T - T∞ ) dV = 0 Assume CDS (average)
dx
dT dT
[( A ) -
dx e
(A ) ]
dx w
- [n2 (TP - T∞ ) A δx] = 0
General equation for all CVs: A cancels out

dT dT 1 General Equation for all
[(dx) - (dx) ] - [n2 (TP - T∞ )δx] = 0 nodes
e w
TE - TP TP - TW
For internal nodes: [( ) - ( )] - [n2 (TP - T∞ )δx] = 0 1 General Equation for internal nodes
δx δx
1 1 1 1 2 2
The internal node equation can be re-arranged as: ( + )T = ( ) TW + ( ) TE + n
⏟ δxT∞ - n
⏟ δxTP
δx δx P δx δx
Su Sp
Equation for internal nodal points 2, 3, and 4: aP TP = aWTW + aE TE + Su
aW aE aP Sp Su
1 1
aW + aE - Sp - n2 δx n2 δxT∞
δx δx

Using the general equation:
dT dT
[(dx) - (dx) ] - [n2 (TP - T∞ )δx] = 0
e w
The equation can be switched to node numbers
T -T TP - TB T -T T - 100
→ [( Eδx P ) - ( )] - [n2 (TP - T∞ )δx] = 0 → [( 2δx 1 ) - ( 1δx/2 )] - [n2 (T1 - 20)δx] = 0
δx/2
This can be re-arranged as: aP TP = aW TW + aE TE + Su
aW aE aP Sp Su
1 2 2
0 aW + aE - Sp - n2 δx - n2 δxT∞ + T
δx δx δx B

dT dT
[(dx) - (dx) ] - [n2 (TP - T∞ )δx] = 0
e w
T -T T -T The equation for boundary node
→ [0 - ( P δx W)] - [n2 (TP - T∞ )δx] = 0 → [0 - ( 5δx 4)] - [n2 (T5 - T∞ )δx] = 0 5 is easier than that of boundary
node 1 due to the insulated end
having a temperature gradient
This can be re-arranged as: aP TP = aW TW + aE TE + Su of 0 on the east face
17
aW aE aP Sp Su
1
0 aW + aE - Sp - n2 δx n2 δxT∞
δx
There are now 5 equations and 5 unknowns
hP
L = 1 m, = 25 m-2
kA
Node aW aE Su Sp aP = aW + aE - Sp
1 0 5 100 + 10TB - 15 20
2 5 5 100 -5 15
3 5 5 100 -5 15
4 5 5 100 -5 15
5 5 0 100 -5 10
20 -5 0 0 0 T1 1100 T1 64.22
-5 15 -5 0 0 T2 100 T2 36.91
0 -5 15 -5 0 T3 = 100 → T3 = 26.50
0 0 -5 15 -5 T4 100 T4 22.60
[ 0 0 0 -5 10 [T5 ]
] [ 100 ] [T5 ] [21.30]
For the numerical solution:
For 5 CVs (5 nodes): δx = 0.2 m
For 10 CVs (10 nodes): δx = 0.1 m (this matches the exact solution)
∴ The solution is more accurate with more CVs
In general Example 4.3 is not a practical example!
Insulation = axis of symmetry: Insulation is usually linked to symmetry
The two situations

alongside are
mathematically equal
The problem on the left has
been solved, therefore, the
problem on the right has been
solved
Example 4.3 with a convection boundary condition at the end

There is no insulation or constant temperature boundary condition at the right end
heat conduction
convection
⏞
d dT
Heat transfer equation: dx
(kA dx) - ⏞
hP(T - T∞ ) = 0
P is the fin perimeter
18
For an insulated boundary condition:

dT dT
[(dx) - (dx) ] - [n2 (TP - T∞ )δx] = 0
e w
T -T T -T
→ [0 - ( P δx W)] - [n2 (TP - T∞ )δx] = 0 → [0 - ( 5δx 4)] - [n2 (T5 - T∞ )δx] = 0
The equation for boundary node 5 is easier than that of boundary node 1 due to the insulated end having a temperature gradient of 0
on the east face
For a convection boundary condition: Now we have convection; TA > 20°C at the boundary

dT dT
[(dx) - (dx) ] - [n2 (TP - T∞ )δx] = 0
e w
T -T T5 - T4
→ A
( Δx/25
) - ( Δx
) - [n2 (TP - T∞ )Δx] = 0 Now we have 5 equations, but 6 unknowns
e w There is no temperature value @ A
Therefore, we need one extra equation: A boundary condition equation for convection can be used
2kT5
dT T -T ( - T∞ )
q̇ e = h(TA - T∞ ) = - k( ) → h(TA - T∞ ) = - k ( A 5) → TA = ∆x
2k
dx e Δx/2 (h + )
∆x
T -T T5 - T4
A
( Δx/25
) - ( Δx
) - [n2 (TP - T∞ )Δx] = 0 The 6 unknowns can now be solved for
e w
Governing Equations
For steady, incompressible, constant property, and Newtonian fluids:
Continuity: div(u) = 0
ρfx
∂p
ρfy
∂p These equations will be on the formula sheet.
y-momentum: ρ div(vu) = - ∂y
These equations used the divergence theorem
ρfz for the surface integral.
∂p
z-momentum: ρ div(wu) = - + μ div(grad w) + SMz
∂z ė g
∂ ∂ ∂
Continuity: div(u) = 0 → ∇∙u=0 → u = ui⃗ + vj⃗ + wk
⃗⃗ → ∇∙u= ∂x
(u) + ∂y
(v) + ∂z
(w)
∂p ∂
+ μ div(grad u) + SMx → 𝜌∇ ∙ (uu) = - ∂x
(p) + μ∇ ∙ (∇u) + SMx
∂ ∂ ∂ ∂ ∂ ∂u ∂ ∂u ∂ ∂u
ρ (∂x (uu) + ∂y
(uv) + ∂z
(uw)) =- ∂x
(p) + μ (∂x ( ∂x ) + ( )
∂y ∂y
+ ( ))
∂z ∂z
+ SMx
Energy Equation: ρc(div(Tu)) = k div(grad T) + Φ + Si → ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T) + Φ + Si

∂ ∂ ∂ ∂ ∂T ∂ ∂T ∂ ∂T
ρ (∂x (uT) + ∂y
(vT) + ∂z
(wT)) = k (∂x (∂x) + ( )
∂y ∂y
+ ( ))
∂z ∂z
+ Φ + Si
19


ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) =0
ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) + ė g = 0
d dT
d dT
∫V ∇ ∙ Φ dV = ∮A ⃗n⃗ Φ dA (Divergence theorem) d
→ → A (final form)
dx
dx e
= ∆x
e w w
[A][T] = [C]
Aw Tw = ? → Aw = ∆y × ⏟t = ∆y and Tw = ?
1

The Central Difference Scheme gives us the average value of the
nodes that the face is in between
Tw =
An Tn = ? → An = ∆x × ⏟t = ∆x and Tn = ?
1
The thickness (t) is usually 1 4 internal nodes = a length of 4Δx

TN + TP The Central Difference Scheme for an unknown at a face is the average value of
Tn =
20

dT dT
Aw (dx) =? → Aw = ∆y × ⏟t = ∆y and (dx) =?
w w
1
dT TThe node after the face in the x-direction - TThe node before the face in the x-direction
( ) =
dx face ∆x
dT TP - TW
(dx) = ∆x
w
dT dT
Ae (dx) = ? → Ae = ∆y × ⏟t = ∆y and (dx) = ?
e e
1
(dx) = ∆x
face
dT TE - TP
(dx) = ∆x
e
∂T ∂T
As (∂y) = ? → As = ∆x × ⏟t = ∆x and (∂y) = ?
s 1 s
∂T TThe node after the face in the y-direction - TThe node before the face in the y-direction
(∂y) = ∆y
face
∂T TP - TS
As ( ) = ∆x ( )
∂y s ∆y
∂T ∂T
An ( ) = ? → An = ∆x × ⏟t = ∆x and ( ) =?
∂y n ∂x n
1
∂T TThe node after the face in the y-direction - TThe node before the face in the y-direction
(∂y) = ∆y
face
∂T TN - TP
An (∂y) = ∆x ( ∆y
)
n
The Finite Volume Method for Two-Dimensional Diffusion Problems

∂ ∂ϕ ∂ ∂ϕ ∂ ∂T ∂ ∂T
∂x
(Γ ∂x ) + ∂y
(Γ ∂y ) +S=0 → ∂x
(k ∂x) + ∂y
(k ∂y) + ė g = 0
∂ ∂ϕ ∂ ∂ϕ
∫∆V ∂x
(Γ ∂x ) dx ∙ dy + ∫∆V ∂y
(Γ ∂y ) dx ∙ dy + ∫∆V Sϕ dV = 0
∂ϕ ∂ϕ ∂ϕ ∂ϕ
[Γe Ae ( ∂x ) - Γw Aw ( ∂x ) ] + [Γn An ( ∂y ) - Γs As ( ∂y ) ] + S̅ ΔV = 0
e w n s
∂ϕ ϕP - ϕW
Flux across the west face = Γw Aw ∂x | = Γw Aw ( δxWP
)
w
∂ϕ ϕE - ϕP
Flux across the east face = Γe Ae ∂x | = Γe Ae ( δx )
e PE
∂ϕ ϕP - ϕS These equations apply for all internal nodes
Flux across the south face = Γs As ∂y | = Γs As ( δy )
s SP
∂ϕ ϕN - ϕP
Flux across the north face = Γn An ∂y | = Γn An ( δyPN
)
n
21
∂ϕ ∂ϕ ∂ϕ ∂ϕ
e w n s
Ae = Aw = Δy
An = As = Δx
Using the Central Difference Method:

ΓP + ΓE ΓW + ΓP ΓP + ΓN ΓP + ΓS
Γe = 2
, Γw = 2
, Γn = 2
, and Γs = 2
∂ϕ ∂ϕ ∂ϕ ∂ϕ
e w n s
ϕE - ϕP ϕP - ϕW ϕN - ϕP ϕP - ϕS
→ Γe Ae ( δxPE
) - Γw Aw ( δxWP
) + Γn An ( δyPN
) - Γs As ( δySP
) + S̅ ΔV = 0
ϕE - ϕP ϕP - ϕW ϕN - ϕP ϕP - ϕS
→ Γe Ae ( ) - Γw Aw ( ) + Γn An ( ) - Γs As ( ) + S̅ ΔV = 0
δx δx δy δy
The source term S may be a function of the dependent variable: S̅ ΔV = Su + Sp ϕP
Therefore, the equation can be rearranged as

Γ Aw
w Γe A e Γs A s Γn A n w Γ Aw e e Γ A ΓA Γ A
( δx + δxPE
+ δySP
+ δyPN
- Sp ) ϕP = ( δx ) ϕW + ( δx ) ϕE + (δys s ) ϕS + (δyn n) ϕN + Su
WP WP PE SP PN
This can be re-arranged as: aP ϕP = aW ϕW + aE ϕE + aS ϕS + aN ϕN + Su
aW aE aS aN aP
Γw Aw Γe Ae Γs As Γn An
aW + aE + aS + aN - Sp
δxWP δxPE δySP δyPN
The results of the process above show that the 1D and 2D cases yield similar results
Finite Volume Method for Three-Dimensional Diffusion Problems

∂ ∂ϕ ∂ ∂ϕ ∂ ∂ϕ
∂x
(Γ ∂x ) + ∂y
(Γ ∂y ) + ∂z
(Γ ∂z ) +S=0 The equations are steady state (there is no time dependence)
∂ ∂T ∂ ∂T ∂ ∂T
For Example: ∂x
(k ∂x) + ∂y
(k ∂y) + ∂z
(k ∂z) + ė g = 0
∂ϕ ∂ϕ ∂ϕ ∂ϕ ∂ϕ ∂ϕ
[Γe Ae ( ∂x ) - Γw Aw ( ∂x ) ] + [Γn An ( ∂y ) - Γs As ( ∂y ) ] + [Γt At ( ∂z ) - Γb Ab ( ∂z ) ] + S̅ ΔV = 0
e w n s t b
ϕE - ϕP ϕP - ϕW ϕ -ϕ ϕ -ϕ
[Γe Ae ( δx
) - Γw Aw (
δx
)] + [Γn An ( Nδy P ) - Γs As ( P S )]
δy
ϕT - ϕP ϕP - ϕB
+ [Γt At ( δz
) - Γb Ab ( δz
)] + (Su + Sp ϕP ) = 0
This can be re-arranged as: aP ϕP = aW ϕW + aE ϕE + aS ϕS + aN ϕN + aB ϕB + aT ϕT + Su
aW aE aS aN aB aT aP
Γw Aw Γe Ae Γs As Γn An Γb Ab Γt At
aW + aE + aS + aN + aB + aT - Sp
δxWP δxPE δySP δyPN δxBP δxPT
Summary of Discretised Equations for Diffusion Problems
aP ϕP = ∑ anb ϕnb + Su
aP = ∑ anb - Sp
22
Dimension aW aE aS aN aB aT
Γw Aw Γe Ae
1D - - - -
δxWP δxPE
Γw Aw Γe Ae Γs As Γn An
2D - -
δxWP δxPE δySP δyPN
Γw Aw Γe Ae Γs As Γn An Γb Ab Γt At
3D
δxWP δxPE δySP δyPN δxBP δxPT
The results of the process above show that the 1D, 2D, and 3D cases yield the same results
∂ ∂ϕ ∂ ∂ϕ
∫∆V ∂x
(Γ ∂x ) dx ∙ dy + ∫∆V ∂y
(Γ ∂y ) dx ∙ dy + ∫∆V Sϕ dV = 0
∂ϕ ∂ϕ
[Γe Ae ( ∂x ) - Γw Aw ( ∂x ) ]
e w
The 3rd step is to integrate the differential equation(s) over the CVs
∫ΔV D.E. dV Incorrect integration impacts the final answer greatly; there are four different cases of differential equations
Differential equation(s) can include:

2
∂ ϕ ϕ can be 5 unknowns and the equations can be made up of 2nd order derivatives,1st order
1. 2 ← second order DE
∂h derivatives and unknown terms.
∂ϕ The ϕ term or the velocity term can be inside the derivative in the energy equation.
2. ∂h
← first order DE
We usually work up to 2D in MKM 411.
3. ϕ ← unknown variable
ϕ can be T, u, v, w. P
h can be x, y, z
Examples of different D.E.s:

2
∂ ϕ d dϕ d dT ∂ ∂T
1. 2 ⇒ ( )
dh dh
⇒ like dx (dx) or ( )
∂y ∂y
∂h
d dT dT dT
∫ΔV dx (k dx) dV = (kA dx) - (kA dx)
⏟ e w
∂ ∂T ∂T ∂T
or ∫ΔV ∂y (k ∂y) dV = (kA ∂y) - (kA ∂y)
⏟ n s
∂ϕ d d
2. ∂h
⇒ dx
(ϕ) ⇒ like dx (u) Take ϕ out
d
∫ΔV dx (u) dV = (Au)e - (Au)w
d
If ϕ = v and h = y → ∫ΔV dy (v) dV = (Av)n - (Av)s
d
If ϕ = uT and h = x → ∫ΔV dx (uT) dV = (AuT)e - (AuT)w
3. ϕ ∫ΔV ϕ dV = ? ϕ can be T, u, v, w. P
The Central Difference Scheme applies in FVM because the assumption is that the CV is small enough
∫ΔV T dV = TP ΔV that the centre of the cell is equal to the average of the cell.
Central Differencing is not a good scheme when there is convection
∫ΔV u dV = uP ΔV
∫ΔV v dV = vP ΔV
23


ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) =0
ρc Dt = Φ + ∇ ∙ ∇(kT) + ġ → dx
(k dx) + ė g = 0
d dT
d dT
∫V ∇ ∙ Φ dV = ∮A ⃗n⃗ Φ dA (Divergence theorem) → → A (final form)
d
dx
dx e
= ∆x
e w w
[A][T] = [C]
Aw Tw = ? → Aw = ∆y × ⏟t = ∆y and Tw = ?
1

Tw =
dT dT
Aw (dx) =? → Aw = ∆y × ⏟t = ∆y and (dx) =?
w w
1
(dx) = ∆x
face
dT TP - TW
( ) =
dx w ∆x
Check yourself-1 Time: 3 minutes
Write the equivalent of the expressions in terms of nodal values and Δx and Δy. Use the Central
dT
Difference Method when needed. An Tn = ?, An (dy) = ?, As vs = ?
n
24
Solution
TP + TN
An Tn = ∆x ( 2
)
dT TN - TP
An (dy) = ∆x ( ∆y
)
n
vP + vs
As vs = ∆x ( 2
)
Next, we will deal with convection
The Finite Volume Method for Convection-Diffusion Problems
div(ρuϕ) = div(Γ grad ϕ) + Sϕ Energy equation
IF ϕ = cp T, Γ = k/cp , ρ = constant: Change to steady state

2
ρcp ∇ ∙ (uT) = ∇ ∙ (k∇T) + ė gen → ρcp (T∇u + u∇T) = k∇ T + ė gen → ρcp (u∇T) = k∇2 T + ė gen ✔
The equation is steady with no heat dissipation (Φ)
Steady convection-diffusion problems: div(ρuϕ) = div(Γ grad ϕ) + Sϕ
Steady continuity equation: div(ρu) = 0
Steady One-Dimensional Convection and Diffusion
Steady convection-diffusion equation:

d d dϕ
(ρuϕ) = (Γ ) 1D without a source term, Δx is fixed
dx dx dx
Steady continuity (mass) equation:

d(ρu)
dx
=0
Integration over the control volume:

d d dϕ
dx
(ρuϕ) = dx
(Γ dx )
∂ϕ ∂ϕ
→ (ρuAϕ)e - (ρuAϕ)w = (ΓA ) - (ΓA )
∂x e ∂x w
Integration over the continuity equation: Velocity does not necessarily move in the positive direction
d(ρu)
dx
=0 → (ρuA)e - (ρuA)w = 0
It is convenient to define two variables, F and D. In general:

Convective mass flux: F = ρuA If velocity is already given, software can now use the velocity in F = ρuA
Γ
Diffusion conductance: D = δx A
Γ
For 1D: F = ρu and D = δx A is omitted
Steady convection-diffusion equation:

d d dϕ ∂ϕ ∂ϕ
dx
(ρuϕ) = dx
(Γ dx ) → (ρuAϕ)e - (ρuAϕ)w = (ΓA ) - (ΓA )
∂x ∂x e w
Steady continuity (mass) equation:

d(ρu)
dx
=0 → (ρuA)e - (ρuA)w = 0
25
F = ρuA Fe = (ρuA)e and Fw = (ρuA)w
Γ Γe Γw
D = A De = A and Dw = A
δx δx e δx w
Continuity equation: (ρuA)e - (ρuA)w = 0 → Fe = (ρuA)e and Fw = (ρuA)w → Fe - Fw = 0
Therefore, for 1D and steady:
∂ϕ ∂ϕ
Convection-diffusion equation: (ρuAϕ)e - (ρuAϕ)w = (ΓA ) - (ΓA )
∂x e ∂x w
ϕE - ϕP ϕP - ϕW
(ρuAϕ)e - (ρuAϕ)w = (ΓA ) - (ΓA )
δx e δx w
Γ Γw
Fe = (ρuA)e Fw = (ρuA)w De = δxe Ae Dw = A
δx w
Fe ϕe - Fw ϕw = De (ϕE - ϕP ) - Dw (ϕP - ϕW)
The Central Differencing Scheme: Used to make sure that all equations have nodal terms and not face terms
ϕe = (ϕP + ϕE )/2 and ϕw = (ϕW + ϕP )/2
Fe Fw
2
(ϕP + ϕE ) - 2
(ϕW + ϕP ) = De (ϕE - ϕP ) - Dw (ϕP - ϕW)
Fw Fe Fw Fe
This can be re-arranged to give: [(Dw - ) + (De + )] ϕP = (Dw + ) ϕW + (De - ) ϕE
2 2 2 2
aP ϕP = aW ϕW + aE ϕE
We would like to find: aP = f(aW , aE )

Fw Fe Fw Fe
[(Dw + 2
) + (De -
2
) + (Fe - Fw )] ϕP = (Dw +
2
) ϕW + (De -
2
) ϕE v
aP ϕP = aW ϕW + aE ϕE
aW aE aP
Fw Fe
Dw + De - aW + aE + (Fe - Fw )
2 2
For 1D: Fe - Fw = 0 Fe = 0 for 1D, Fe ≠ 0 for 2D
Example
Consider the following energy equation in 2D
ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T)
If the velocity field is: The velocity is not unknown in this instance – it is given
The velocity vector is 2D; the
u(u, v) = 2i⃗ + 3j⃗ velocity field is not a function
of x and y
By applying the Finite Volume Method (FVM), find the general
algebraic equation for internal nodes. Use the Central Difference
Scheme where needed. Assume the fluid is incompressible, and k, c are constants. Δx ≠ Δy
26
Solution
What should we do first? We have the velocity field, therefore, we should check the continuity
equation:
∂u ∂v ∂(2) ∂(3)
∂x
+ ∂y
=0 → ∂x
+ ∂y
=0 → 0-0=0✔
∂ ∂ ∂ ∂T ∂ ∂T
ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T) → ρc (∂x (uT) + ∂y
(vT)) = k (∂x (∂x) + ( ))
∂y ∂y
∂ ∂ ∂ ∂T ∂ ∂T
→ ∫ ρc (∂x (uT) + ∂y
(vT)) dV = ∫ k (∂x (∂x) + ( )) dV
∂y ∂y
The general equation for all nodes:

∂T ∂T ∂T ∂T
ρc((AuT)e - (AuT)w + (AvT)n - (AvT)s ) = k ((A ∂x) - (A ∂x) + (A ∂y) - (A ∂y) )
e w n s
∂T ∂T ∂T ∂T
→ ρc(Ae ue Te - Aw uw Tw + An vn Tn - As vs Ts ) = k (Ae (∂x) - Aw (∂x) + An (∂y) - As (∂y) )
e w n s
T +T T +T T +T T +T
→ ρc (∆yue E P - ∆yuw W P + ∆xvn N P - ∆xvs S P )
2 2 2 2
TE - TP TP - TW TN - TP TP - TS
= k (∆y ( ∆x
) - ∆y ( ∆x
) + ∆x ( ∆y
) - ∆x ( ∆y
))
u(u, v) = 2i⃗ + 3j⃗

TE + TP TW + TP TN + TP TS + TP
ρc (2∆y 2
- 2∆y 2
+ 3∆x 2
- 3∆x 2
)
→ TE - TP TP - TW TN - TP TP - TS
= k (∆y ( ) - ∆y ( ) + ∆x ( ) - ∆x ( ))
∆x ∆x ∆y ∆y
Example
ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T)
If the velocity field is: The velocity is not unknown in this instance – it is given
The velocity field is now

u(u, v) = (2x + y)i⃗ + (x - 2y)j⃗ a function of x and y

Scheme where needed. Assume the fluid is incompressible, and k, c are constants. Δx = Δy
Solution
What should we do first? We have the velocity field, therefore, we should check the continuity
equation:
∂u ∂v ∂(2x + y) ∂(x - 2y)
∂x
+ ∂y
=0 → ∂x
+ ∂y
=0 → 2-2=0✔
∂ ∂ ∂ ∂T ∂ ∂T
ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T) → ρc ( (uT) + (vT)) = k( ( ) + ( ))
∂x ∂y ∂x ∂x ∂y ∂y
∂ ∂ ∂ ∂T ∂ ∂T
→ ∫ ρc (∂x (uT) + ∂y
(vT)) dV = ∫ k (∂x (∂x) + ( )) dV
∂y ∂y
27

∂T ∂T ∂T ∂T
e w n s
∂T ∂T ∂T ∂T
e w n s
ue = (2x + y)e
Ae = Aw = ∆y
uw = (2x + y)w
u(u, v) = (2x + y)i⃗ + (x - 2y)j⃗ →
vn = (x - 2y)n
An = As = ∆x
vs = (x - 2y)s
ρc(∆yue Te - ∆yuw Tw + ∆xvn Tn - ∆xvs Ts )
∂T ∂T ∂T ∂T
= k (∆y ( ) - ∆y ( ) + ∆x ( ) - ∆x ( ) )
∂x e ∂x w ∂y n ∂y s
TE + TP T +T T +T TS + TP
→ ρc (∆yue - ∆yuw W P + ∆xvn N P - ∆xvs )
2 2 2 2
= k (∆y ( ) - ∆y ( ) + ∆x ( ) - ∆x ( ))
∆x ∆x ∆y ∆y
Δx = Δy, ue = (2x + y)e , uw = (2x + y)w , vn = (x - 2y)n , vs = (x - 2y)s

→ ρc∆x (ue - uw + vn - vs ) = k((TE - 2TP + TW ) + (TN - 2TP + TS ))
2 2 2 2
→ ρc∆x ((2x + y)e - (2x + y)w + (x - 2y)n - (x - 2y)s ) = k(TE + TW + TN + TS - 4TP )
2 2 2 2
ρc (∆yue 2
- ∆yuw 2
+ ∆xvn 2
- ∆xvs 2
)
ue = ?, uw = ?, vn = ?, vs = ?
= k (∆y ( ∆x
) - ∆y ( ∆x
) + ∆x ( ∆y
) - ∆x ( ∆y
))
Assuming: ∆x = ∆y = 10 cm with u(u, v) = (2x + y)i⃗ + (x - 2y)j⃗

ue = (2x + y)e = 2xe + ye = 2(0.15) + 0.1 = 0.4 m/s
uw = (2x + y)w = 2xw + yw = 2(0.05) + 0.1 = 0.2 m/s
vn = (x - 2y)n = xn - 2yn = 0.1 - 2(0.15) = - 0.2 m/s
vs = (x - 2y)s = xs - 2ys = 0.1 - 2(0.05) = 0 m/s
ρc∆x ((2x + y)e 2
- (2x + y)w 2
+ (x - 2y)n 2
- (x - 2y)s 2
) = k(TE + TW + TN + TS - 4TP )
ue = 0.4 m/s, uw = 0.2 m/s, vn = - 0.2 m/s, vs = 0 m/s

TE + TP TW + TP TN + TP
→ ρc(0.1) (0.4 2
- 0.2 2
- 0.2 2
) = k(TE + TW + TN + TS - 4TP )
aP TP = aW TW + aE TE + aN TN + aS TS
Example
ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T)

Scheme where needed. Assume the fluid is incompressible, and
k, c are constants. Δx = Δy
If only one node is asked, the algebraic equation only needs to be determined for that
node
28
Solution
∂ ∂ ∂ ∂T ∂ ∂T
ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T) → ρc (∂x (uT) + ∂y
(vT)) = k (∂x (∂x) + ( ))
∂y ∂y
∂ ∂ ∂ ∂T ∂ ∂T
→ ∫ ρc (∂x (uT) + ∂y
(vT)) dV = ∫ k (∂x (∂x) + ( )) dV
∂y ∂y

∂T ∂T ∂T ∂T
e w n s
∂T ∂T ∂T ∂T
e w n s
T +T T +T T +T T +T
ρc (∆yue E P - ∆yuw W P + ∆xvn N P - ∆xvs S P)
2 2 2 2
→ TE - TP TP - TW TN - TP TP - TS
= k (∆y ( ∆x
) - ∆y ( ∆x
) + ∆x ( ∆y
) - ∆x ( ∆y
))
ue = ?, uw = ?, vn = ?, vs = ?
uE + uP uW + uP vN + vP vS + vP
ue = 2
, uw = 2
, vn = 2
, vs = 2
Example
By applying the Finite Volume Method (FVM) and the Central
Difference Scheme, where needed, find the algebraic equation
for the energy equation for nodes 1 and 3.
Assume: Steady, 2D, air is incompressible and has constant

properties, no heat generation, no slip condition, Δx ≠ Δy, and
the viscous heat dissipation function is given as:
∂v
Φ = μ ∂y + B(Tv) + D
Incompressible → the standard energy equation can be used

No slip → practically means that there is viscosity, μ ≠ 0, velocity at the wall can be zero
Heat dissipation → B and D are constants in this instance (not always constants in other
instances); μ is constant
Solution
∂T
Energy equation for incompressible with constant properties: ρc ( + ∇ ∙ (Tu)) = k∇ ∙ (∇T) + Φ + ė gen
∂t
In the case of steady and no heat generation: ρc(∇ ∙ (Tu)) = k∇ ∙ (∇T) + Φ
∂ ∂ ∂ ∂T ∂ ∂T ∂v
2D with the given Φ: ρc (∂x (uT) + ∂y
(vT)) = k (∂x (∂x) + ( ))
∂y ∂y
+ μ ∂y + B(Tv) + D
∂ ∂ ∂ ∂T ∂ ∂T ∂v
∫ ρc (∂x (uT) + ∂y
(vT)) dV = ∫ k (∂x (∂x) + ( )) dV
∂y ∂y
+ ∫ (μ ∂y + B(Tv) + D) dV

ρc((AuT)e - (AuT)w + (AvT)n - (AvT)s )
∂T ∂T ∂T ∂T
= k ((A ∂x) - (A ∂x) + (A ∂y) - (A ∂y) ) + μ((Av)n - (Av)s ) + B(TP vP )∆V + D∆V
e w n s
∆V is usually equal to ∆x∆y∆t = ∆x∆y(1) = ∆x∆y
29
Node 1
uair = 4 cos θ m/s It is very important the write the information of the specific node
∂T
∂T ( ) = 0 because
Tw = Tair = 30°C, uw = 4 cos θ m/s, ( ) = 0, vs = 0 ∂y s
∂y s of the insulation
at the south face
ρc((AuT)e - (AuT)w + (AvT)n - (AvT)s )
∂T ∂T ∂T ∂T
= k ((A ∂x) - (A ∂x) + (A ∂y) - (A ∂y) ) Assume
∂T
≠
∂T
n
e w n s ∂x ∂y
e
+ μ((Av)n - (Av)s ) + B(TP vP )∆V + D∆V
w
Insulation does not mean that the temperature at the face is zero s
uE + uP T +T vN + vP T +T
ρc (∆y ( ) ( E P) - ∆y(4 cos θ)(30) + ∆x ( ) ( N P) - ∆x(0)Ts )
2 2 2 2
TE - TP TP - TW TN - TP vN + vP
= k (∆y ( ∆x
) - ∆y ( ∆x/2
) + ∆x ( ∆y
) - ∆x(0)) + μ∆x (( 2
) - (0)) + B(TP vP )∆x∆y + D∆x∆y
The face w is at Δx/2 (halfway between nodes that are Δx apart)
u2 + u1 T +T v4 + v1 T +T
ρc (∆y ( 2
) ( 2 2 1) - ∆y(4 cos θ)(30) + ∆x ( 2
) ( 2 2 1 ))
T2 - T1 1 T - 30 T4 - T1 v4 + v1
= k (∆y ( ∆x
) - ∆y ( ∆x/2 ) + ∆x ( ∆y
)) + μ∆x ( 2
) + B(T1 v1 )∆x∆y + D∆x∆y
Node 3 Different to node 1; equation will be different
vair = - 4 sin θ m/s
Te = 0°C, ue = 0 m/s, Tn = 30°C, vn = - 4 sin θ m/s
ρc((AuT)e - (AuT)w + (AvT)n - (AvT)s ) n

∂T ∂T ∂T ∂T e
= k ((A ) - (A ) + (A ) - (A ) )
∂x e ∂x w ∂y n ∂y s w
+ μ((Av)n - (Av)s ) + B(TP vP )∆V + D∆V s
Face e is important: T = 0
Must switch face terms to nodal terms
uW + uP T +T vS + vP T +T
ρc (∆y(0)Te - ∆y ( ) ( W P) + ∆x(- 4 sin θ)(30) - ∆x ( ) ( S P ))
2 2 2 2
ET -T
P TP - TW N P T -T TP - TS vP + vS
= k (∆y ( ∆x/2 ) - ∆y ( ∆x
) + ∆x ( ∆y/2 ) - ∆x ( ∆y
)) + μ∆x (- 4 sin θ - 2
) + B(TP vP )∆x∆y + D∆x∆y
u4 + u3 T +T v2 + v3 T +T
ρc (- ∆y ( 2
) ( 4 2 3) + ∆x(- 4 sin θ)(30) - ∆x ( 2
) ( 2 2 3 ))
0 - T3 T3 - T4 30 - T3 T3 - T2 v3 + v2
= k (∆y ( ) - ∆y ( ) + ∆x ( ) - ∆x ( )) + μ∆x (- 4 sin θ - ) + B(T3 v3 )∆x∆y + D∆x∆y
∆x/2 ∆x ∆y/2 ∆y 2
30

Theme 3 Finite Volume Method (FVM)

Uploaded by

Document Informationclick to expand document information

Copyright:

Available Formats

Theme 3 Finite Volume Method (FVM)

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Theme 3 Finite Volume Method (FVM)

Uploaded by

Copyright:

Available Formats

MKM 411 – Computational Fluid Dynamics (CFD) Theme 3

Finite Volume Method (FVM)

In FVM We Take Integration from the Equations

For example, if we consider the x-direction of the momentum equation:

We can’t take vx out (it is not constant)

For Incompressible and Constant Properties

We need to change the format of the momentum and energy equations

Du ∂u ∂u ∂u ∂u ∂u The equation will be given

SMx is the source term (e.g., body force):

Governing Equations in the CFD Textbook

In the textbook, the velocity vector is indicated as u: u = ui⃗ + vj⃗ + wk

SMx is the source term (e.g., body force):

Du ∂u ∂u ∂u ∂u ∂u The equation will be given

For Incompressible, Constant Property, and Newtonian Fluids

z-momentum as already determined:

Conservation of Energy (in the CFD Textbook)

For incompressible, constant property, and Newtonian fluids i = cT

For solids: k div(grad T) + Si = 0

Therefore, for steady, incompressible, constant property, and Newtonian fluids:

Continuity: div(u) = 0 ρfx

Energy Equation: ρc(div(Tu)) = k div(grad T) + Φ + Si

Finite Volume Method

General Methodology: Finite Volume = Basic Methodology

In general, the FVM involves the following steps:

1. Divide the Domain into Control Volumes (CVs)

→ We use capital letters for nodes North

For example: One-dimensional steady heat conduction without heat generation

3. Integrate the Differential Equation Over the CVs

Each CV gives one equation and finally, the final number of

That’s why we use these forms; for steady, incompressible, constant

Energy Equation: ρc(div(Tu)) = k div(grad T) + Φ + Si

5. Approximation of the Surface Integrals

6. Approximation of Function Values and Derivatives by Interpolation with Nodal Values

7. Assembling and Solution of the Discrete Algebraic System

→ The result is a set of linear algebraic equations

Finite Volume Method (Fast Review Up to Now)

2. Finding (or simplifying) the differential equation(s) for the CVs

For internal nodes:

Nodes 1 and 5 are boundary nodes

For boundary node 5:

Equation for boundary node 5: aP TP = aW TW + aE TE + Su Source term

The resulting set of algebraic equations for this example is:

T(0) = 100: T(0) = C1 (0) + C2 = 100 ∴ C2 = 100

T(0.5) = 500: T(0.5) = C1 (0.5) + 100 = 500 ∴ C1 = 800

∴ Exact solution: T(x) = 800x + 100

Finite Volume Method (Fast Review Up to Now)

2. Finding (or simplifying) the differential equation(s) for the CVs

Central Difference Scheme

Do not mix Ae Te with:

Let’s Review the Regulation on Faces (up to now)

We do not have any nodes at the faces; therefore, we must

Central Difference Scheme

Let’s Review the Regulation on Faces (up to now)

We do not have any nodes at the faces; therefore, we must

Integrate the Differential Equation (s) over the CVs

Differential equation(s) can include:

Examples of different D.E.s:

An Example with a Source Term

ė g = 200T + 1000 W/m3