THE FINITE DIFFERENCE METHOD IN THE 2D FOURIER EQUATION WITH ROBIN’S BOUNDARY CONDITION
Grzegorz Biernat, Sylwia Lara-Dziembek, Edyta Pawlak Institute of Mathematics, Czestochowa University of Technology
Częstochowa, Poland
grzegorz.biernat@im.pcz.pl, sylwia.lara@im.pcz.pl, edyta.pawlak@im.pcz.pl
Abstract. This paper contains the application of the Finite Difference Method in the two-dimensional Fourier equation using Robin’s boundary condition (the third boundary condition).
Keywords: block matrices, n-band matrices, determinant, Fourier equation, boundary condition
Introduction
In this paper we consider the Fourier equation describing the heat transfer in a two-dimensional domain. This equation is supplemented by the third boundary condition, which is often called Robin’s condition.
This condition is the most “natural” condition, which can be assumed on the parts of area adjacent to the environment. It constitutes a mathematical notation of the heat flux's continuity and it is based on the law of Newton.
In the paper we give the direct FDM formula for the solution of the Fourier equation with Robin’s boundary condition.
1. Solution of the problem
We consider the 2D Fourier equation [1-3]
( ) ( ) ( )
2 2
2 2
, , , , , ,
T x y t T x y t T x y t
x y c t
λ ∂ ∂ ρ ∂
+ =
∂ ∂ ∂
(1)
where λ is a thermal conductivity [W/mK], c is a specific heat [J/kgK], ρ is a mass density [kg/m
3], T is temperature [K], , x y denote the geometrical co-ordinates and t is time [s].
Then approximations of the second order partial derivatives using FDM are as follows
( )
( )
2
1, , , , 1, ,
2 2
2
, 1, , , , 1,
2 2
2 , 1 1
2 , 1 1
i j l i j l i j l
i j l i j l i j l
T T T
T i m
x x
T T T
T j n
y y
− +
− +
− +
∂ = ≤ ≤ −
∂ ∆
− +
∂ = ≤ ≤ −
∂ ∆
(2)
The initial condition has the form
, , 0
, 0
i j ini
T = T ≤ ≤ l q (3)
where T
iniis the initial temperature.
Below we present a system of equations containing the Fourier equation with the third boundary conditions [2]
( )
( )
1, ,
0, ,
,1,
, 0,
2 2
2 2
2
2
j l env
j l env
i l env
i l env
T T
T T
x
T T
T T
y
T T T
x y c t
λ α
λ α
λ ρ
−
= −
∆
−
= −
∆
∆ ∆ ∆
∆ + ∆ = ∆
(4)
where α is the heat transfer coefficient and T
envis the ambient temperature (the temperature of environment).
Successively we consider the systems of equations corresponding to the four verti- ces and lateral edges of our domain, which is a rectangle.
For i = 0, j = 0 we have
( )
( )
( ) ( ) ( ) ( ) ( ) ( )
1, 0,
0, 0,
0,1,
0, 0,
0, 0, 1, 0, 0, 0, 0,1,
2 2 2 2 2 2
0, 0, 0, 0, 1
2
2
2 2
l env
l env
l env
l env
env l l env l l
l l
T T
T T
x
T T
T T
y
T T T T T T
x x x y y y
T T
c t
λ α
λ α
λ λ λ λ λ λ
ρ
− −
= −
∆
−
= −
∆
− + + − + =
∆ ∆ ∆ ∆ ∆ ∆
−
=
∆
(5) After transformations we get the following form
( ) ( ) ( ) ( )
( ) ( )
0, 0, 1, 0,
0, 0, 0,1,
0, 0, 1, 0, 0,1,
2 2 2 2
0, 0, 1
2 2
2 2
2 2
2 2
l l env
l l env
l l l
env l
T T T
x x
T T T
y y
c T T T
x y t x y
T c T
x y t
λ λ
α α
λ λ
α α
λ λ ρ λ λ
λ λ ρ
−
− = −
∆ ∆
− = −
∆ ∆
+ + − − =
∆ ∆ ∆ ∆ ∆
= + +
∆ ∆ ∆
(6)
The main matrix above the system of equations is as follows
( )
2( )
2( )
2( )
20 2
0
2
2 2
x
A y
c
x y t x y
α λ α λ
λ λ ρ λ λ
−
∆
−
= ∆
+ + − −
∆ ∆ ∆ ∆ ∆
(7)
On the basis of the determinant of the matrix A we obtain the limitations on the
steps of the differential grid [4]
0 2
0 2
det 0
x y A α λ α λ
− >
∆
− >
∆
>
(8)
Hence
2 2
x y
λ λ
α α
λ λ
α α
< ∆ <
< ∆ <
(9)
The determinant of the main matrix is expressed by the following formula
( ) ( )
2
2 2
1 1
det
2 2
A c
x y x y t x y
λ λ λ ρ
α
= ∆ ∆ ∆ + ∆ + ∆ − ∆ + ∆
(10)
Then we have the determinant corresponding to the variable T
0,0,l( ) ( )
0 ,0 ,
2
0, 0, 1
2 2
1 1
det
2 2
T l env l
A T c T
x y x y x y t
λ λ λ ρ
α
−
= ∆ ∆ ∆ + ∆ − ∆ + ∆ + ∆
(11) Thus, the temperature of node (0,0) at the time l is as follows
( ) ( )
( ) ( )
0, 0, 1
2 2
0, 0,
2 2
1 1
2
1 1
2
env l
l
T c T
x y t
x y
T c
x y t
x y
λ λ ρ
α
λ λ ρ
α
−
+ − + +
∆ ∆ ∆ ∆ ∆
=
+ − ∆ + ∆ + ∆
∆ ∆
(12)
Taking the following notation
( )
2( )
21 1
, , ,
2
a b c d c
x y t
x y
λ λ ρ
α
= ∆ = ∆ = ∆ + ∆ = ∆ (13)
we have
0, 0,l env 0, 0,l1
a b c d
T T T
a b c d a b c d
−= + − +
+ − + + − +
(14)
Similarly, we proceed for the other vertices: (m, 0), (m, n), (0, n).
Then we repeat the procedure for each edge.
And so, on the edge (i, 0), where 1 ≤ i ≤ m − 1 , j = 0 we receive
( )
( )
( ) ( ) ( ) ( ) ( ) ( )
1, 0, 1, 0,
, 0,
,1,
, 0,
1, 0, , 0, 1, 0, , 0, ,1,
2 2 2 2 2 2
, 0, , 0, 1
2
2
2 2
i l i l
i l env
i l env
i l env
i l i l i l env i l i l
i l i l
T T
T T
x
T T
T T
y
T T T T T T
x x x y y y
T T
c t
λ α
λ α
λ λ λ λ λ λ
ρ
+ −
− +
−
−
= −
∆
−
= −
∆
− + + − + =
∆ ∆ ∆ ∆ ∆ ∆
−
=
∆
(15) Therefore
( ) ( )
2
2 2
1 1
det
2 2
A c
x y x y t x y
λ λ λ ρ
α
= ∆ ∆ ∆ + ∆ + ∆ − ∆ + ∆
(16)
( ) ( )
,0 ,
2
1, 0, , 0, 1
2 2
1 1
det
2 2
i l
T env i l i l
A T T c T
x y y x y x t
λ λ λ ρ
α
− −
= ∆ ∆ ∆ − ∆ + ∆ + ∆ + ∆
(17)
Thus, the temperature on the edge (i, 0) 1 ≤ i ≤ m − 1 , j = 0 at the moment l is given
by the formula
( ) ( )
( ) ( )
1, 0, , 0, 1
2 2
, 0,
2 2
1 1
2
1 1
2
env i l i l
i l
T T c T
x y t
y x
T c
x y t
x y
λ λ ρ
α
λ λ ρ
α
− −
− + + +
∆ ∆ ∆ ∆ ∆
=
+ − ∆ + ∆ + ∆
∆ ∆
(18)
Assuming markings (13) we obtain
, 0, 1, 0, , 0, 1
i l env i l i l
b c a d
T T T T
a b c d a b c d
−a b c d
−= − + +
+ − + + − + + − +
(19)
Analogously we calculate the temperature for the other edge points:
– 1 ≤ ≤ i m − 1, j = n – i = 0, 1 ≤ ≤ − j n 1 – i = m , 1 ≤ ≤ − . j n 1
The system of equations in the internal points of the area takes the following form
( ) ( ) ( )
( ) ( ) ( )
1, , , , 1, ,
2 2 2
, 1, , , , 1,
2 2 2
, , , , 1
2
2
i j l i j l i j l
i j l i j l i j l
i j l i j l
T T T
x x x
T T T
y y y
c c
T T
t t
λ λ λ
λ λ λ
ρ ρ
− +
− +
−
