26/02/2020
MAT6669 Finite difference 1
MAT6669
Numerical modelling of heat transfer
Part 1
An explicit finite difference scheme
Magnus Anderson and Karl Travis
Units and nomenclature
𝑘
𝜌
𝑐𝑝 𝛼 𝑚 𝑄 𝜀
𝜎 𝐴
thermal conductivity density
specific heat Diffusivity
Mass
Heat flux
Emissivity Stefan-Boltzmann constant Area
distance
Temperature
Environment temperature Heat transfer coefficient
Latent heat
Volume fraction of transformed phase
(𝐽𝑠−1𝑚−1𝐾−1) (𝑘𝑔 𝑚−3)
( 𝐽 𝑘𝑔−1 𝐾−1) ( 𝑚2𝑠−1)
(𝑘𝑔)
(𝐽𝑠−1) (𝑑𝑖𝑚𝑒𝑛𝑠𝑖𝑜𝑛𝑙𝑒𝑠𝑠) (𝐽𝑠−1𝑚−2𝐾−4) (𝑚2)
(𝑚)
(𝐾)
(𝐾) (𝐽𝑠−1𝑚−2𝐾−1)
(𝐽 𝑘𝑔−1) (𝑑𝑖𝑚𝑒𝑛𝑠𝑖𝑜𝑛𝑙𝑒𝑠𝑠)
𝛥𝑥 𝑇 𝑇
∞ h
𝛼 𝐿 𝛽
𝑓𝛼
26/02/2020
MAT6669 Heat and Materials with Application 2
Heat transfer
We are interested in the temperature evolution within a volume that has known initial temperature and boundary conditions. The exterior of the volume loses heat with the environment through convection and radiation.
𝑅𝑎𝑑𝑖𝑎𝑡𝑖𝑜𝑛
𝑄=𝜀𝜎𝐴 𝑇4 −𝑇4 ∞
𝑄 = −𝑘 𝐴 𝐶𝑜𝑛𝑑𝑢𝑐𝑡𝑖𝑜𝑛
𝛥𝑇 𝛥𝑥
C𝑜𝑛𝑣𝑒𝑐𝑡𝑖𝑜𝑛 𝑄=h𝐴𝑇 −𝑇
∞
𝐿𝑎𝑡𝑒𝑛𝑡 h𝑒𝑎𝑡 𝑓𝑟𝑜𝑚 𝑝h𝑎𝑠𝑒 𝑡𝑟𝑎𝑛𝑠𝑓𝑜𝑟𝑚𝑎𝑡𝑖𝑜𝑛𝑠
𝛼 𝑑𝑓𝛼 𝑄=𝑚𝐿𝛽 𝑑𝑡
26/02/2020
MAT6669 Heat and Materials with Application
3
Temperature evolution
Interior
𝜌𝑐 𝜕𝑇− 𝜕 𝑘𝜕𝑇 − 𝜕 𝑘𝜕𝑇 =0 𝑝𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦
Conduction, heat generation through latent heat or chemical reactions.
𝑠𝑢𝑚 𝑜𝑓h𝑒𝑎𝑡 𝑓𝑙𝑢𝑥𝑒𝑠
Surfaces
𝑛
𝑑𝑇
𝑖 𝑝𝑑𝑡
Heat fluxes from conduction, convection, radiation, latent heat, chemical reactions.
=𝑄 =𝑚𝐶
26/02/2020
MAT6669 Heat and Materials with Application 4
𝑎𝑐𝑡𝑖𝑛𝑔 𝑤𝑖𝑡h𝑖𝑛 𝑐𝑜𝑛𝑡𝑟𝑜𝑙 𝑣𝑜𝑙𝑢𝑚𝑒
𝑖=1
The finite difference method
The finite difference method is derived through the manipulation of Taylor series to obtain approximations of derivatives. A Taylor series is a method of approximating a function using a polynomial. The following YouTube video provides a good description of the Taylor series;
The following is a Taylor series considering a forward difference of h
𝑓𝑥+h =𝑓𝑥 +𝜕𝑓h+1𝜕2𝑓h2+1𝜕3𝑓h3+1𝜕4𝑓h4+⋯
𝜕𝑥 2! 𝜕𝑥2 3! 𝜕𝑥3 4! 𝜕𝑥4
We can rearrange this equation for 𝜕𝑓 𝜕𝑥
𝜕𝑓=𝑓𝑥+h −𝑓𝑥 −1𝜕2𝑓h−1𝜕3𝑓h2−1𝜕4𝑓h3-… 𝜕𝑥 h 2!𝜕𝑥2 3!𝜕𝑥3 4!𝜕𝑥4
We can group the remaining terms into a term descriptive of the error
𝜕𝑓=𝑓 𝑥+h −𝑓 𝑥 +𝒪(h) 𝜕𝑥 h
The following YouTube videos describe the derivation of a central finite difference scheme
https://www.youtube.com/watch?v=f67jpwJu-d0 https://www.youtube.com/watch?v=r419a4VYi8o
26/02/2020 MAT6669 Heat and Materials with Application 5
Finite difference equations First derivative
Two point, forward
𝜕𝑢𝑖 = 𝑢𝑖+1 − 𝑢𝑖 + 𝒪 h 𝜕𝑥 h
Two point, backward
Second derivative
Three point, forward
𝜕2𝑢𝑖 = 𝑢𝑖+2 − 2 𝑢𝑖+1 + 𝑢𝑖 + 𝒪 h 𝜕𝑥2 h2
Three point, backward
𝜕𝑢𝑖 = 𝑢𝑖 − 𝑢𝑖−1 + 𝒪 h
𝜕𝑥h 𝜕𝑥2h2
𝜕2𝑢𝑖 = 𝑢𝑖 − 2 𝑢𝑖−1 + 𝑢𝑖−2 + 𝒪 h
Three point, forward
𝜕𝑢𝑖 =−3𝑢𝑖 +4𝑢𝑖+1 −𝑢𝑖+2 +𝒪 h2 𝜕𝑥 2h
Three point, backward
Three point, central
𝜕2𝑢𝑖 =𝑢𝑖+1 −2𝑢𝑖 +𝑢𝑖−1 +𝒪 h2 𝜕𝑥2 h2
Four point forward
Two point, central Four point, backward
𝜕𝑢𝑖 =3𝑢𝑖 −4𝑢𝑖−1 −𝑢𝑖−2 +𝒪 h2
𝜕2𝑢𝑖 =−𝑢𝑖+3 +4𝑢𝑖+2 −5𝑢𝑖+1 +2𝑢𝑖 +𝒪 h2 𝜕𝑥2h 𝜕𝑥2 h2
𝜕𝑢𝑖 = 𝑢𝑖+1 −𝑢𝑖−1 +𝒪 h2 𝜕2𝑢𝑖 −𝑢𝑖−3 +4𝑢𝑖−2 −5𝑢𝑖−1 +2𝑢𝑖 2 𝜕𝑥2h 𝜕𝑥2=h2+𝒪h
Four point, central
𝜕𝑢𝑖 =𝑢𝑖−2 −8𝑢𝑖−1 +8𝑢𝑖+1 −𝑢𝑖+2 +𝒪 h4 𝜕𝑥 12h
26/02/2020
MAT6669 Heat and Materials with Application 6
Temperature representation
It is possible to discretize the domain into small sections of widths Δ𝑥 and Δ𝑦 and then describe the temperature using an array or a vector. For plotting the temperature in matlab, the temperature must be expressed as a 2D array. To manipulate the problem to the generic form 𝐴 𝑥 = 𝑏, the temperature must be expressed in a vector.
Array – matrix form
𝑇𝑇𝑇 1,1 1,2 1,3
𝑇=𝑇𝑇𝑇 2,1 2,2 2,3
Array – vector form
𝑇 1,1
𝑇𝑥1 𝑦1 1,2 Position 𝑥2 𝑦1
𝑇𝑇𝑇
𝑇 array 1,3
𝑇 2,1
𝑥3 𝑦1
𝑥1 𝑦2
3,1
3,2 3,3
𝑇=
𝑇𝑥𝑦 2,2 𝑙𝑜𝑐𝑎𝑡𝑖𝑜𝑛 = 2 2
𝑇 𝑥3 𝑦2
Position vectors
𝑥𝑥𝑥 2,3 𝑥𝑦 𝑥=123𝑇13
𝑦1
𝑇𝑇𝑇 1,1 1,2 1,3
3,1
𝑇 3,2
𝑇 3,3
𝑥2 𝑦3
𝑥𝑦 3 3
𝑦𝑇𝑇𝑇 𝑦= 2 2,1 2,2 2,3
𝑦𝑇𝑇𝑇
3 3,1 3,2 3,3
26/02/2020
MAT6669 Heat and Materials with Application
7
Temperature discretization
Consider an array matrix form representation of the temperature field with 𝑛 many rows and 𝑚 many columns. Special conditions are needed to define the kinetics the boundaries to account for heat transfer through radiation and convection. For a 3 point finite difference scheme, the following special cases exist
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 1,1 1,2 1,3 1,… 1,𝑗 1,… 1,𝑚−2 1,𝑚−1 1,𝑚
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 2,1 2,2 2,3 2,… 2,𝑗 2,… 2,𝑚−2 2,𝑚−1 2,𝑚
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 3,1 3,2 3,3 3,… 3,𝑗 3,… 3,𝑚−2 3,𝑚−1 3,𝑚
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 …,1 …,2 …,3 …,… …,𝑗 …,… …,𝑚−2 …,𝑚−1 …,𝑚
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 𝑖,1 𝑖,2 𝑖,3 𝑖,… 𝑖,𝑗 𝑖,… 𝑖,𝑚−2 𝑖,𝑚−1 𝑖,𝑚
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 …,1 …,2 …,3 …,… …,𝑗 …,… …,𝑚−2 …,𝑚−1 …,𝑚
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 𝑛−2,1 𝑛−2,2 𝑛−2,3 𝑛−2,… 𝑛−2,𝑗 𝑛−2,… 𝑛−2,𝑚−2 𝑛−2,𝑚−1 𝑛−2,𝑚
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 𝑛−1,1 𝑛−1,2 𝑛−1,3 𝑛−1,… 𝑛−1,𝑗 𝑛−1,… 𝑛−1,𝑚−2 𝑛−1,𝑚−1 𝑛−1,𝑚
𝑻=
𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇
𝑛,1 𝑛,2
𝑛,3 𝑛,…
𝑛,𝑗 𝑛,…
𝑛,𝑚−2 𝑛,𝑚−1 𝑛,𝑚
Convection Radiation Conduction Heat generation
Corners
Special cases
Edges
Interior nodes
Conduction Heat generation
26/02/2020
MAT6669 Heat and Materials with Application
8
Finite difference schemes: Interior points
Let us distinguish the temperature fields at different points in time in addition to space. Let 𝑇𝑝 refer to 𝑖,𝑗
the temperature at the current time and 𝑇p+1 the temperature after a time of Δ𝑡. 𝑖,𝑗
Next time step
Current time step
𝜕𝑇 𝑇p+1−𝑇𝑝
𝜕𝑡
= 𝑖,𝑗
Δ𝑡
𝑖,𝑗+𝒪(Δ𝑡) Row below
Location of interest
𝑇p 𝑖,𝑗−1
𝑇p 𝑖−1,𝑗
𝑝 𝑇𝑝
𝑇 𝑖,𝑗+1 𝑖,𝑗
𝑇𝑝 𝑖+1,𝑗
𝜕2𝑇 𝜕𝑥2
=
𝑇𝑝 𝑖+1,𝑗
Row above
−2𝑇𝑝 +𝑇𝑝
𝑖,𝑗 𝑖−1,𝑗 + 𝒪(Δ𝑥2)
Δ𝑥2
Location of interest
Column to the right
Column to the left
𝜕2𝑇 𝜕𝑦2
=
𝑇𝑝 𝑖,𝑗+1
−2𝑇𝑝 +𝑇𝑝
𝑖,𝑗 𝑖,𝑗−1 + 𝒪(Δ𝑦2)
Δ𝑦2
26/02/2020
MAT6669 Heat and Materials with Application
9
Interior nodes: Finite difference scheme
The 2D heat equation
𝜌𝐶 𝜕𝑇− 𝜕 𝑘𝜕𝑇 − 𝜕 𝑘𝜕𝑇 =0 𝑝𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦
Apply the chain rule:
𝜌𝐶 𝜕𝑇−𝑘𝜕2𝑇−𝑘𝜕2𝑇−𝜕𝑘𝜕𝑇−𝜕𝑘𝜕𝑇=0 𝑝 𝜕𝑡 𝜕𝑥2 𝜕𝑦2 𝜕𝑥𝜕𝑥 𝜕𝑦𝜕𝑦
Assume negligble spatial gradients in thermal conductivity
𝜌𝐶 𝜕𝑇−𝑘𝜕2𝑇−𝑘𝜕2𝑇=0 𝑝 𝜕𝑡 𝜕𝑥2 𝜕𝑦2
Let𝛼= 𝑘 𝜌𝑐𝑝
𝜕𝑇 − 𝛼 𝜕2𝑇 − 𝛼 𝜕2𝑇 = 0 𝜕𝑡 𝜕𝑥2 𝜕𝑦2
26/02/2020
MAT6669 Heat and Materials with Application 10
Interior nodes: Finite difference scheme
𝜕𝑇 − 𝛼 𝜕2𝑇 − 𝛼 𝜕2𝑇 = 0 𝜕𝑡 𝜕𝑥2 𝜕𝑦2
Use 1st order forward FD scheme for temporal derivative and 2nd order central FD scheme for spatial derivative
𝑇𝑝+1−𝑇𝑝 𝑇𝑝
𝑖,𝑗
−2𝑇𝑝 +𝑇𝑝 𝑇𝑝
𝑖,𝑗 𝑖−1,𝑗 −𝛼 𝑖,𝑗+1
−2𝑇𝑝 +𝑇𝑝
𝑖,𝑗 𝑖,𝑗−1 =0
Δ𝑦2
𝑖,𝑗 −𝛼 𝑖+1,𝑗 𝑇p+1−𝑇𝑝 𝑇𝑝
Δ𝑥2
Assume Δ𝑥 = Δ𝑦
−2𝑇𝑝 +𝑇𝑝 𝑇𝑝
𝑖,𝑗 𝑖−1,𝑗 −𝛼 𝑖,𝑗+1
Rearrange for 𝑇p+1 𝑖,𝑗
𝑖,𝑗
𝑖,𝑗 −𝛼 𝑖+1,𝑗
−2𝑇𝑝 +𝑇𝑝
Δ𝑡
Δ𝑡
𝑖,𝑗 Δ𝑥2 Δ𝑥2
𝑖,𝑗−1 =0
𝑇p+1=𝑇𝑝 +F 𝑇𝑝
𝑖,𝑗 𝑖,𝑗 o 𝑖+1,𝑗
+𝑇𝑝 +𝑇𝑝 +𝑇𝑝 −4𝑇𝑝
Fourier number
𝐹 =𝛼 Δ𝑡
𝑜 Δ𝑥 2
𝑖−1,𝑗 𝑖,𝑗+1 𝑖,𝑗−1
𝑖,𝑗
26/02/2020 MAT6669 Heat and Materials with Application
11
Boundary elements
Consider the
temperature at the
boundary node 𝑇𝑝 𝑖,𝑚
𝑇𝑝 𝑖,𝑚
26/02/2020
MAT6669 Heat and Materials with Application
12
Boundary elements
𝑇𝑝 𝑖,𝑚−1
𝑇𝑝 𝑖−1,𝑚
𝑇𝑝 𝑖,𝑚
𝑇𝑝 𝑖+1,𝑚
We can use finite
difference to
determine the change
in temperature of 𝑇𝑝 𝑖,𝑛
considering its neighbouring nodes.
26/02/2020
MAT6669 Heat and Materials with Application 13
Boundary elements
𝑇𝑝 𝑖,𝑚−1
𝑇𝑝 𝑖−1,𝑚
𝑇p 𝑖,𝑚
𝑇p 𝑖+1,𝑚
What about heat loss or gain through convection and radiation?
26/02/2020
MAT6669 Heat and Materials with Application 14
Boundary elements
1Δ𝑥 2
𝑇𝑝 𝑖,𝑚−1
𝑇𝑝 𝑖−1,𝑚
𝑇𝑝 𝑖,𝑚
2 𝑖+1,𝑚
The width of the geometry is given by 𝜔
Consider the heat flowing through the following volume.
Note the dimensions (1 Δ𝑥 by Δ𝑦 )
𝑇𝑝
26/02/2020
MAT6669 Heat and Materials with Application 15
Δ𝑦
Boundary elements
The diffusive fluxes entering and leaving the node 𝑇 are given by
𝑄3
𝑄 ,𝑄 ,𝑄 ,and𝑄 . 123 4
Conservation of energy dictates that 𝑑𝑇
𝑄 +𝑄 +𝑄 +𝑄 =𝑚𝐶 1234 𝑝𝑑𝑡
Where 𝑚 is the mass, and 𝐶 is the 𝑝
specific heat capacity.
An additional flux may be added for heat generated through chemical reactions / latent heat.
𝑚,𝑛
𝑄2
𝑄4
𝑄1
26/02/2020
MAT6669 Heat and Materials with Application 16
Boundary elements
𝑑𝑇 1234 𝑝𝑑𝑡
𝑄 +𝑄 +𝑄 +𝑄 =𝑚𝐶
𝑇𝑝 𝑖−1,𝑚
The diffusive fluxes are given by Fourier’s law
𝑄1 = k
𝑝𝑝
𝑄3
Δ𝑥 𝑇 −𝑇 𝜔 𝑖−1,𝑚 𝑖,𝑚
2 Δ𝑦
𝑇𝑝 −𝑇𝑝 𝑄2 = k Δ𝑦𝜔 𝑖,𝑚−1 𝑖,𝑚
𝑇𝑝 𝑄
𝑖,𝑚−12𝑝4 Δ𝑥
𝑄
𝑇𝑝𝑝
𝑖,𝑚 Δ𝑥 𝑇 −𝑇 𝑄3 = k 𝜔 𝑖−1,𝑚 𝑖,𝑚
𝑄1
𝑖+1,𝑚
𝑇𝑝
Area
Where the bracketed term on the left is the area, with 𝜔 being the width of the geometry.
2 Δ𝑦
26/02/2020
MAT6669 Heat and Materials with Application 17
Boundary elements
𝑇𝑝
𝑄 𝑖−1,𝑚
3
𝑇𝑝𝑄2 𝑄4
𝑇𝑝 𝑖,𝑚
The heat lost or gained through convection is given by
𝑄4 = 𝑄𝑟𝑎𝑑 + 𝑄𝑐𝑜𝑛
𝑄 =hΔ𝑦𝜔𝑇−𝑇𝑝 𝑐𝑜𝑛 ∞ 𝑚,𝑛
4 𝑟𝑎𝑑 ∞ 𝑚,𝑛
𝑖,𝑚−1
𝑄 =𝜀𝜎Δ𝑦𝜔𝑇4−𝑇𝑝
𝑄1
𝑖+1,𝑚
Where 𝑇 is the environment ∞
temperature.
h is the heat transfer coefficient
(𝑊 𝑚−2𝐾−1) .
The terms 𝜀 and 𝜎 refer to the
emissivity and the Stefan-Boltzmann
constant
𝑇𝑝
26/02/2020
MAT6669 Heat and Materials with Application 18
Explicit scheme
𝑄 +𝑄 +𝑄 +𝑄 =𝑚𝐶
Let us ignore radiation, and only consider convection heat losses/gains, 𝑄4 = 𝑄𝑐𝑜𝑛
𝑇𝑝 −𝑇𝑝 Δ𝑥 𝑇𝑝 −𝑇𝑝 𝑖,𝑚−1 𝑖,𝑚 + k 𝜔 𝑖−1,𝑚 𝑖,𝑚
+hΔ𝑦𝜔 𝑇 −𝑇𝑝
∞ 𝑚,𝑛 𝑝
Δ𝑡
𝑑𝑇 1 2 3 4 𝑝𝑑𝑡
Δ𝑥 𝑇𝑝 −𝑇𝑝
k 𝜔 𝑖−1,𝑚 2 Δ𝑦
𝑖,𝑚
+ k Δ𝑦𝜔
𝑇p+1 − 𝑇𝑝
𝑖,𝑚
The mass is given by 𝑚 = 1 Δ𝑥Δ𝑦𝜔𝜌 2
=𝑚𝐶 𝑖,𝑚
Δ𝑥 2 Δ𝑦
Δ𝑥 𝑇p −𝑇𝑝
𝑇p −𝑇𝑝
𝑖,𝑚−1 𝑖,𝑚 +k
Δ𝑥 𝑇𝑝 −𝑇𝑝
k
𝜔 𝑖−1,𝑚 2 Δ𝑦
𝑖,𝑚
+k Δ𝑦𝜔 1
𝜔 𝑖−1,𝑚 Δ𝑥 2 Δ𝑦
𝑖,𝑚
+hΔ𝑦𝜔 𝑇 −𝑇𝑝
∞ 𝑚,𝑛 2 𝑝 Δ𝑡
𝑇p+1 − 𝑇𝑝 = Δ𝑥Δ𝑦𝜔𝜌𝐶 𝑖,𝑚 𝑖,𝑚
26/02/2020 MAT6669 Heat and Materials with Application
19
Explicit scheme
Δ𝑥 𝑇p −𝑇𝑝
𝑇𝑝 −𝑇𝑝 Δ𝑥 𝑇𝑝 −𝑇𝑝
Divide by 𝜔 and assume Δx = Δ𝑦
+k Δ𝑥 1
k 𝑖−1,𝑚 2 Δ𝑥
𝑖,𝑚
𝑖,𝑚−1 𝑖,𝑚 +k 𝑖−1,𝑚
𝑖,𝑚
= Δ𝑥2𝜌𝐶 𝑖,𝑚 𝑖,𝑚 ∞ 𝑚,𝑛 2 𝑝 Δ𝑡
+hΔ𝑥 𝑇 −𝑇𝑝
𝑇p+1=𝑇𝑝 +𝐹 2𝐵𝑇 +𝑇𝑝 +2𝑇𝑝
Δ𝑥 2 Δ𝑥 𝑇p+1 − 𝑇𝑝
Simplify and rearrange for 𝑇𝑛+1 𝑖,𝑚
𝑖,𝑚
𝑖,𝑚 0
𝑖 ∞
𝑖−1,𝑚 𝑖,𝑚−1
+𝑇𝑝 𝑖−1,𝑚
− 4+2𝐵 𝑇𝑝 𝑖 i,m
𝛼=
𝑘 𝜌𝑐𝑝
Fourier number
𝐹=𝛼 Δ𝑡 𝑜 Δ𝑥2
Biot number
𝐵=Δ𝑥h 𝑖 𝑘
26/02/2020
MAT6669 Heat and Materials with Application
20
Explicit scheme: Left edge
𝑇1 , 1 𝑇2,1 𝑇3,1 𝑇…,1
𝑻= 𝑇𝑖,1 𝑇…,1
𝑇𝑛 − 2 , 1 𝑇𝑛 − 1 , 1 𝑇𝑛 , 1
𝑇1,2 𝑇2,2 𝑇3,2 𝑇…,2 𝑇𝑖,2 𝑇…,2
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , 𝑗 𝑇2,𝑗 𝑇3,𝑗 𝑇…,𝑗 𝑇𝑖,𝑗 𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄4 𝑖,1
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 2
𝑄1
Δ𝑥 𝑇𝑝 −𝑇𝑝
1Δ𝑥
2 𝑇𝑝
𝑖,2
Δ𝑦
p
𝑇 2Δ𝑦
𝑖−1,1
𝑄1 = k
𝜔 𝑖+1,1 𝑖,1
𝑄 =hΔ𝑦𝜔 𝑇 −𝑇𝑝
𝑇𝑝 𝑖,1
𝑇𝑝 −𝑇𝑝 𝑄3=k Δ𝑥𝜔 𝑖−1,1 𝑖,1
2 Δ𝑦
2∞
𝑖,1
𝑇𝑝 𝑝𝑝 𝑖+1,1 𝑇 −𝑇
𝑄4 =k Δ𝑦𝜔
𝑖,2
𝑖,1
Δ𝑥
26/02/2020 MAT6669 Heat and Materials with Application
21
Explicit scheme: Right edge
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , 𝑗 𝑇2,𝑗 𝑇3,𝑗 𝑇…,𝑗 𝑇𝑖,𝑗 𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑖,𝑚 𝑄1
4
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 2
Δ𝑥 𝑇𝑝 −𝑇𝑝 𝜔 𝑖+1,𝑚 𝑖,𝑚
𝑇𝑝 −𝑇𝑝 𝑖,𝑚−1 𝑖,𝑚
1Δ𝑥
𝑇𝑝 𝑖,𝑚−1
2 Δ𝑦
𝑇𝑝 𝑖−1𝑚
𝑇𝑝 𝑖,𝑚
𝑇𝑝 𝑖+1,𝑚
𝑄1 = k
𝑄2 = k Δ𝑦𝜔
2 Δ𝑦
Δ𝑥
Δ𝑥 𝑇𝑝 −𝑇𝑝
𝜔 𝑖−1,𝑚
2 Δ𝑦
𝑖,𝑚
𝑄3 = k
𝑄 =hΔ𝑦𝜔 𝑇 −𝑇𝑝
4∞
𝑖,𝑚
26/02/2020
MAT6669 Heat and Materials with Application
22
Explicit scheme: Top edge
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄 1,𝑗 4
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 2
𝑇1 , 𝑗
𝑇2,𝑗
𝑇3,𝑗
𝑇…,𝑗
𝑇𝑖,𝑗 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1
𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑄1
𝑝𝑝𝑝 𝑇 𝑇 𝑇
𝑇𝑝 −𝑇𝑝 2,𝑗 1,𝑗
Δ𝑦
Δ𝑦 𝑇𝑝 −𝑇𝑝 𝑄2 = k 𝜔 1,𝑗−1 1,𝑗
2 Δ𝑥 𝑄 =hΔ𝑥𝜔 𝑇 −𝑇𝑝
Δ𝑦 𝑇𝑝 −𝑇𝑝 𝑄4 = k 𝜔 1,𝑗+1 1,𝑗
2 Δ𝑥
1,𝑗−1 1,𝑗
1,𝑗+1
𝑄1 = k Δ𝑥𝜔
1Δ𝑦 2
𝑇𝑝 2,𝑗
Δ𝑥
3∞
1,𝑗
26/02/2020
MAT6669 Heat and Materials with Application
23
Explicit scheme: Bottom edge
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , 𝑗 𝑇2,𝑗 𝑇3,𝑗 𝑇…,𝑗 𝑇𝑖,𝑗 𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑛,𝑗 𝑄1
4
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 2
𝑇𝑝 𝑛−1,𝑗
1Δ𝑦 2
Δ𝑥
𝑄 =hΔ𝑥𝜔 𝑇 −𝑇𝑝
1∞
𝑛,𝑗
𝑄2 = k
𝑄3 =k Δ𝑥𝜔
Δ𝑦 𝑇𝑝 −𝑇𝑝 𝜔 𝑛,𝑗−1 𝑛,𝑗
𝑇𝑝 −𝑇𝑝 𝑛−1,𝑗 𝑛,𝑗
2 Δ𝑥
Δ𝑦
Δ𝑦 𝑇𝑝 −𝑇𝑝
26/02/2020
MAT6669 Heat and Materials with Application
24
𝑇𝑝 𝑛,𝑗+1
𝑇𝑝 𝑛,𝑗
𝑇𝑝 𝑛,𝑗−1
𝑄4 = k
2
𝜔
𝑛,𝑗+1 𝑛,𝑗 Δ𝑥
Explicit scheme: Edge summary
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑇1 , 𝑗
𝑇2,𝑗
𝑇3,𝑗
𝑇…,𝑗
𝑇𝑖,𝑗 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1
𝑘
𝐹=𝛼 Δ𝑡
𝐵=Δ𝑥h
𝛼=
𝑝 Δ𝑥 𝑘
Left edge
𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝜌𝑐
𝑜𝑖 2
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑇p+1 =𝑇𝑝 +𝐹 2𝐵𝑇 +𝑇𝑝 +𝑇𝑝 +2𝑇𝑝 − 4+2𝐵 𝑇𝑝
𝑖,1 𝑖,1 𝑜 𝑖∞ 𝑖+1,1 𝑖−1,1 𝑖,2
Right edge
𝑖 𝑖,1
𝑇p+1 =𝑇𝑝 +𝐹 2𝐵𝑇 +𝑇p +2𝑇𝑝 +𝑇𝑝
𝑖,𝑚 𝑖,𝑚 𝑜 𝑖 ∞ 𝑖+1,𝑚 𝑖,𝑚−1 𝑖−1,𝑚
Top edge
𝑇p+1 =𝑇𝑝 +𝐹 2𝐵𝑇 +2𝑇p +𝑇𝑝
1,𝑗 1,𝑗 𝑜 𝑖∞ 2,𝑗 1,𝑗−1 1,𝑗+1
Bottom edge
𝑇p+1 =𝑇𝑝 +𝐹 2𝐵𝑇 +𝑇𝑝 +2𝑇p +𝑇𝑝
𝑛,𝑗 𝑛,𝑗 𝑜 𝑖 ∞ 𝑛,𝑗−1 𝑛−1,𝑗 𝑛,𝑗+1
− 4+2𝐵 𝑇𝑝
𝑖 𝑖,𝑚
+𝑇p
− 4+2𝐵 𝑇𝑝 𝑖 1,𝑗
− 4+2𝐵 𝑇𝑝 𝑖 𝑛,𝑗
26/02/2020
MAT6669 Heat and Materials with Application
25
Top left corner
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇𝑝 1,1
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , 𝑗 𝑇2,𝑗 𝑇3,𝑗 𝑇…,𝑗 𝑇𝑖,𝑗 𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄 1,1 4
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 4
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑄1
𝑇𝑝 𝑝𝑝
1,2
Δ𝑥 𝑇−𝑇 𝑄1 = k 𝜔 2,1 1,1
1Δ𝑦 2
1Δ𝑥 2
2 Δ𝑦 𝑄=hΔ𝑦𝜔 𝑇−𝑇𝑝
𝑇𝑝 2,1
22∞
𝑄 =h Δ𝑥𝜔 𝑇 −𝑇𝑝
32∞
1,1
Δ𝑦 𝑇𝑝 −𝑇𝑝 𝑄4 = k 𝜔 1,2 1,1
2 Δ𝑥
1,1
26/02/2020
MAT6669 Heat and Materials with Application
26
Top right corner
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , 𝑗 𝑇2,𝑗 𝑇3,𝑗 𝑇…,𝑗 𝑇𝑖,𝑗 𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄 1,𝑚 4
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 4
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑇𝑝 1,𝑚−1
𝑇𝑝 1,𝑚
Δ𝑥 𝑇𝑝−𝑇𝑝 𝜔 2,𝑚 1,𝑚
2 Δ𝑦
𝑄1
1Δ𝑦 2
1Δ𝑥 2
𝑄1 = k
Δ𝑦 𝑇𝑝 −𝑇𝑝 𝜔 1,𝑚−1 1,𝑚
𝑇𝑝 2,𝑚
2 Δ𝑥 32∞
𝑄2 = k
𝑄 =h Δ𝑥𝜔 𝑇 −𝑇𝑝
1,𝑚
𝑄 =h Δ𝑦𝜔 𝑇 −𝑇𝑝
26/02/2020
MAT6669 Heat and Materials with Application
27
42∞
1,𝑚
Bottom right corner
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , 𝑗 𝑇2,𝑗 𝑇3,𝑗 𝑇…,𝑗 𝑇𝑖,𝑗 𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄 𝑛,𝑚 4
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 4
𝑄=hΔ𝑥𝜔 𝑇−𝑇𝑝
1 2 ∞ 𝑛,𝑚
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑄1
1Δ𝑥 2
1Δ𝑦 2
𝑇𝑝 𝑛−1,𝑚
𝑄2 = k
𝑄3 = k
Δ𝑦 𝑇𝑝 −𝑇𝑝 𝜔 𝑛,𝑚−1 𝑛,𝑚
2 Δ𝑥
Δ𝑥 𝑇𝑝 −𝑇𝑝 𝜔 𝑛−1,𝑚 𝑛,𝑚
2 Δ𝑦
26/02/2020
MAT6669 Heat and Materials with Application 28
𝑇𝑝 𝑛,𝑚−1
𝑇p 𝑛,𝑚
𝑄 =h Δ𝑦𝜔 𝑇 −𝑇𝑝
4 2 ∞ 𝑛,𝑚
Bottom left corner
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , 𝑗 𝑇2,𝑗 𝑇3,𝑗 𝑇…,𝑗 𝑇𝑖,𝑗 𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑄3
𝑄2 𝑇𝑝 𝑄 𝑛,1 4
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑚 = 1 Δ𝑥Δ𝑦𝜔 𝜌 4
𝑄1
𝑄 =h Δ𝑥𝜔 𝑇 −𝑇𝑝
1Δ𝑥 2
1Δ𝑦 2
𝑇𝑝 𝑛−1,1
Δ𝑥 𝑇𝑝 −𝑇𝑝 𝜔 𝑛−1,1 𝑛,1
12∞
𝑄 =h Δ𝑦𝜔 𝑇 −𝑇𝑝
22∞
𝑛,1
𝑛,1
𝑄3 = k 𝑄4 = k
2 Δ𝑦 Δ𝑦 𝑇𝑝 −𝑇𝑝
𝜔 𝑛,2
2 Δ𝑥
𝑛,1
𝑇𝑝 𝑛,1
p 𝑛,2
26/02/2020
MAT6669 Heat and Materials with Application
29
𝑇
Explicit scheme: Corner summary
𝑇1 , 1 𝑇1,2 𝑇2,1 𝑇2,2 𝑇3,1 𝑇3,2 𝑇…,1 𝑇…,2
𝑻= 𝑇𝑖,1 𝑇𝑖,2 𝑇…,1 𝑇…,2 𝑇𝑛 − 2 , 1 𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 1 𝑇𝑛 − 1 , 𝑇𝑛 , 1 𝑇𝑛 , 2
𝑇1,3 𝑇2,3 𝑇3,3 𝑇…,3 𝑇𝑖 , 3 𝑇…,3
𝑇𝑛 − 2 , 𝑇𝑛 − 1 , 𝑇𝑛 , 3
𝑇1,… 𝑇2,… 𝑇3,… 𝑇…,… 𝑇𝑖,… 𝑇…,… 𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇1 , … 𝑇1 , 𝑚 − 2 𝑇1 , 𝑚 − 1 𝑇2,… 𝑇2,𝑚−2 𝑇2,𝑚−1 𝑇3,… 𝑇3,𝑚−2 𝑇3,𝑚−1 𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝑇1,𝑚 𝑇2,𝑚 𝑇3,𝑚 𝑇…,𝑚 𝑇𝑖,𝑚 𝑇…,𝑚 𝑇𝑛 − 2 , 𝑚 𝑇𝑛 − 1 , 𝑚 𝑇𝑛 , 𝑚
𝑇1 , 𝑗
𝑇2,𝑗
𝑇3,𝑗
𝑇…,𝑗
𝑇𝑖,𝑗 𝑇𝑖,… 𝑇𝑖,𝑚−2 𝑇𝑖,𝑚−1
𝑘
𝐹=𝛼 Δ𝑡
𝐵=Δ𝑥h
𝛼=
𝑝 Δ𝑥 𝑘
Top left corner
𝑇…,𝑗 𝑇𝑛 − 2 , 𝑗 𝑇𝑛 − 1 , 𝑗 𝑇𝑛 , 𝑗
𝑇…,… 𝑇…,𝑚−2 𝑇…,𝑚−1
𝜌𝑐
𝑜𝑖 2
2 2
3 3
𝑇𝑛 − 2 , … 𝑇𝑛 − 1 , … 𝑇𝑛 , …
𝑇𝑛 − 2 , 𝑚 − 2 𝑇𝑛 − 1 , 𝑚 − 2 𝑇𝑛 , 𝑚 − 2
𝑇𝑛 − 2 , 𝑚 − 1 𝑇𝑛 − 1 , 𝑚 − 1 𝑇𝑛 , 𝑚 − 1
𝑇𝑝+1 =𝑇𝑝 +𝐹 4𝐵𝑇 +2𝑇𝑝 +2𝑇𝑝 −4 1+𝐵 𝑇𝑝
1,1 1,1 𝑜 𝑖 ∞ 2,1 1,2
Top right corner
𝑇𝑝+1=𝑇𝑝 +F 4𝐵𝑇 +2𝑇𝑝 +2𝑇𝑝
1,𝑚 1,𝑚 o 𝑖∞ 2,𝑚 1,𝑚−1
Bottom right corner
𝑖 1,1
−41+𝐵 𝑇𝑝
𝑖 1,𝑚
𝑇𝑝+1=𝑇𝑝 +F4𝐵𝑇+2𝑇𝑝 +2𝑇𝑝
𝑛,𝑚 𝑛,𝑚 o 𝑖∞ 𝑛−1,𝑚 𝑛,𝑚−1
Bottom left corner
−41+𝐵𝑇𝑝
𝑖 𝑛,𝑚
𝑇𝑝+1=𝑇𝑝 +F 4𝐵𝑇 +2𝑇𝑝 +2𝑇𝑝 −41+𝐵 𝑇𝑝 𝑛,1 𝑛,1 o 𝑖∞ 𝑛−1,1 𝑛,2 𝑖 𝑛,1
26/02/2020
MAT6669 Heat and Materials with Application
30
Time step
We are using a first order accurate finite difference approximation for the temporal derivative.
Significant error will occur when the following criteria are not met
Δ𝑥 2
𝛼Δ𝑡 ≥4
This may be rearranged to determine an appropriate time step as a function of
the material properties:
Δ𝑥 2 4𝛼
Δ𝑡 =
This does not account for the heat transfer coefficient. The time step may need to be reduced further depending upon the Biot number, 𝐵𝑖.
26/02/2020 MAT3820 Introduction to FEM 31