The 3D Fourier equation with the Robin’s boundary condition usingthe finite difference metod

(1)

THE 3D FOURIER EQUATION WITH THE ROBIN’S BOUNDARY CONDITION USING THE FINITE DIFFERENCE METHOD

Grzegorz Biernat, Sylwia Lara-Dziembek, Edyta Pawlak

Institute of Mathematics, Czestochowa University of Technology

Czestochowa, Poland

grzegorz.biernat@im.pcz.pl, sylwia.lara@im.pcz.pl, edyta.pawlak@im.pcz.pl

Abstract. This article is a continuation of the discussion on the two-dimensional Fourier equation with Robin’s boundary condition. In the paper the Finite Difference Method is applied.

Keywords: block matrices, determinant, Fourier equation, boundary condition

Introduction

In this paper we consider the heat flow in the three-dimensional area described by the Fourier equation, supplemented by the third boundary condition.

The paper presents the exact solution of the problem.

Solution of the problem

We consider the 3D Fourier equation [1-3]

( ) ( ) ( ) ( )

2 2 2

, , , , , , , , , , , ,

T x y z t T x y z t T x y z t T x y z t

c t

x y z

λ ρ

 ∂ ∂ ∂  ∂

+ + =

 

 ∂ ∂ ∂  ∂

 

(1)

where λ is a thermal conductivity [W/mK], c is a specific heat [J/kgK], ρ is a mass density [kg/m

³

], T is temperature [K], x, y, z denote the geometrical co-ordinates and t is time [s].

The second differential partial derivatives using FDM are as follows:

( )

2

1, , , , , , 1, , ,

2 2

2

, 1, , , , , , 1, ,

2 2

2

, , 1, , , , , , 1,

2 2

2 , 1 1

i j k l i j k l i j k l

T T T

T i m

x x

T T T

T j n

y y

T T T

T k p

z z

− +

∆ = ≤ ≤ −

∆ ∆

− +

∆ = ≤ ≤ −

∆ ∆

− +

∆ = ≤ ≤ −

∆ ∆

(2)

and the first differential derivative of the time has a form

, , , , , , 1

i j k l i j k l

, 1

T T

T l q

t t

−

∆ = ≤ ≤

∆ ∆ (3)

The initial condition has the following form

, , ,0

, 0

i j k ini

T = T ≤ ≤ l q (4)

where T

_ini

is the initial temperature.

The system of equations containing the Fourier equation with the Robin conditions [2] is presented

( )

1, , ,

0, , ,

,1, ,

,0, ,

, ,1,

, ,0,

2 2 2

2

j k l env

i k l env

i j l env

T T

x

T T

y

T T

z

T T T T

c t

x y z

λ α

λ ρ

 −

= −

 ∆

  −

 = −

 ∆

 −

 = −

 ∆

  

 ∆ ∆ ∆ ∆

+ + =

 

  ∆ ∆ ∆  ∆



(5)

where α is the heat transfer coefficient and T

_env

is the ambient temperature of envi- ronment.

Now we consider the systems of equations corresponding to the eight vertices and six lateral planes of our cube.

In the first vertex of the cube (for i = 0, j = 0, k = 0) we have

( )

( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

1,0,0,

0,0,0,

0,1,0,

0,0,0,

0,0,1,

0,0,0,

0,0, 0, 1, 0,0, 0, 0,0,

2 2 2 2 2

0,1, 0, 0,0,0, 0, 0

2 2 2 2

2

2 2 2

2

l env

env l l env l

l env l

T T

x

T T

y

T T

z

T T T T T

x x x y y

T T T T

y z z z

λ α

λ λ λ λ λ

λ λ λ λ

− = −

∆

− = −

∆

− = −

∆

− + + − +

∆ ∆ ∆ ∆ ∆

+ + − +

∆ ∆ ∆ ∆

0,0,0, 0,0,0, 1 ,1,

l l

l

T T

c t

ρ

⁻

 

 

 

 

  −

 =

 ∆



(6)

(3)

Therefore we obtain the system of equations

( ) ( ) ( )

( )

0,0,0, 1, 0,0,

0,0,0, 0,1, 0,

0,0,0, 0, 0,1,

0, 0,0,

2 2 2

1,0,0, 0,1,0, 0,0,1,

2 2 2

2

2 2

2 2 2

l l env

l

l l l

T T T

x x

T T T

y y

T T T

z z

c T

x y z t

T T T

x y z

x

λ λ

α α

λ λ

α α

λ λ

α α

λ λ λ ρ

λ λ λ

λ

 

− =  − 

∆  ∆ 

 

− ∆ =   − ∆  

 

− ∆ =   − ∆  

 

 + + + 

 ∆ ∆ ∆ ∆ 

 

− − − =

∆ ∆ ∆

= ∆ ( )

²

( )

² ^env ^0,0,0,^l ¹

T c T

y z t

λ λ ρ

−

 

 

 

 

 

   

  + +  +

  ∆ ∆  ∆

  



(7)

The main matrix system above has the form

( )

²

( )

²

( )

²

( )

²

( )

²

( )

²

0 0

2 0 0

2 2 2 2

x

y A

z c

x y z t x y z

α λ α λ α λ

λ λ λ ρ λ λ λ

 − 

 ∆ 

 

 − 

 ∆ 

 

=  − 

 ∆ 

 

+ + + − − −

 

∆ ∆ ∆ ∆ ∆ ∆ ∆

 

 

(8)

Using the determinant of the matrix A we obtain the limitations on the steps of the differential grid [4]

0 2

det 0

x y z A α λ α λ α λ

 − >

 ∆

 

− >

 ∆

 

 − >

 ∆

 >



(9)

(4)

So

2 2 2

x y z

λ λ

α α

λ λ

α α

λ λ

α α

 < ∆ <

 

 < ∆ <

 

 < ∆ <

 

(10)

The determinant of the matrix A is given by the following formula

( ) ( ) ( )

3

2 2 2

1 1 1

det

4 2

A c

x y z x y z t x y z

λ λ λ λ ρ

α

   

 

= ∆ ∆ ∆   ∆ + ∆ + ∆ + ∆ −   ∆ + ∆ + ∆     (11)

Next, we calculate the determinant corresponding to the variable T

0,0,0,l

( ) ( ) ( )

0, 0, 0,

3

0,0,0, 1

2 2 2

det

1 1 1

4 2

T l

env l

A

T c T

x y z x y z x y z t

λ λ λ λ ρ

α

−

=

     

  + + − + +  + 

   

∆ ∆ ∆     ∆ ∆ ∆  ∆ ∆ ∆    ∆  

(12) The temperature of node (0, 0, 0) at the time l is as follows

( ) ( ) ( )

0,0,0, 1

2 2 2

0,0,0,

2 2 2

1 1 1

2 1 1 1

2

env l

l

T c T

x y z t

x y z

T c

x y z t

x y z

λ λ λ ρ

α

λ λ λ ρ

α

−

   

 + + −  ∆ + ∆ + ∆   + ∆

 ∆ ∆ ∆   

 

=  

+ + −  ∆ + ∆ + ∆  + ∆

∆ ∆ ∆  

(13)

According to following notation

( )

²

( )

²

( )

²

1 1 1

, , , ,

2 a b c d e c

x y z t

x y z

λ λ λ ρ

α  

= ∆ = ∆ = ∆ =   ∆ + ∆ + ∆   = ∆ (14)

we have

( )

0, 0,0, 1

0,0,0,

env l

l

a b c d T eT

T

a b c d e

+ + − +

−

= + + − +

(15) or

0,0,0,l env 0,0,0,l 1

a b d

T T T

a b d a b d

⁻

= − +

− + − +

(16)

(5)

We proceed for the other vertices, similarly.

Then we repeat the procedure for each edge.

And so, on the edge (i, 0, 0), where 1 ≤ ≤ i m − 1, j = 0, k = we receive 0

( )

( ) ( ) ( )

( ) ( )

1,0, 0, 1,0,0,

,0,0,

,1,0,

,0, 0,

,0,1,

,0, 0,

1, 0,0, , 0,0, 1,0, 0,

2 2 2

,0, 0, ,1,0,

2 2 2

2 2

2

i l i l

i l env

i l i l i l

env i l i l

env

T T

x

T T

y

T T

z

T T T

x x x

T T T

y y y

T

z z

λ α

λ λ λ

λ λ

+ −

− +

− = −

∆

− = −

∆

− = −

∆

− + +

∆ ∆ ∆

+ − + +

∆ ∆ ∆

+ −

∆ ∆ ( )

,0, 0, ,0, 0, 1

,0,0, 2 , 0,1,

i l i l

T T

T T c

z t

λ ρ

⁻

 

 

 

 

 −

+ =

 ∆ ∆



(17)

Therefore

( ) ( ) ( )

3

2 2 2

1 1 1

det

4 2

A c

x y z x y z t x y z

λ λ λ λ ρ

α

   

 

= ∆ ∆ ∆   ∆ + ∆ + ∆ + ∆ −   ∆ + ∆ + ∆     (18)

( ) ( ) ( )

, 0, 0,

3

1, 0,0, , 0,0, 1

2 2 2

det

4 1 1 1

2

i l

T

env i l i l

A x y z

T T c T

x y z t

y z x

λ

λ λ λ ρ

α

− −

= ⋅

∆ ∆ ∆

     

   

⋅      ∆ + ∆ −   ∆ + ∆ + ∆     + ∆ + ∆    (19)

Ultimately, the temperature on the edge (i, 0, 0) 1 ≤ ≤ i m − 1, j = 0, k = at the 0 moment l is given by the formula

( ) ( ) ( )

,0,0,

1,0,0, ,0,0, 1

2 2 2

1 1 1

2 1 1 1

2

i l

env i l i l

T

T T c T

x y z t

y z x

c

x y z t

x y z

λ λ λ ρ

α

λ λ λ ρ

α

− −

=

   

 + −  ∆ + ∆ + ∆   + + ∆

 ∆ ∆    ∆

 

=  

+ + −  ∆ + ∆ + ∆  + ∆

∆ ∆ ∆  

(20)

(6)

Assuming markings (14) we obtain

,0, 0, 1,0, ,0, 1

i l env i l i l

b c d a e

T T T T

a b c d e a b c d e

⁻

a b c d e

⁻

= + − + +

+ + − + + + − + + + − +

(21) Analogously we calculate the temperature for the other edge points.

Next, on the plane (i, j, 0), where 1 ≤ ≤ i m − 1, 1 ≤ ≤ − j n 1, k = we receive 0

( ) ( ) ( )

1, ,0, , ,0, 1, ,0,

2 2 2

, 1,0, , ,0, , 1,0,

2 2 2

, ,0, , ,1,

2 2 2

, ,0, , ,0, 1

2

i j l i j l i j l

env i j l i j l

i j l i j l

T T T

x x x

T T T

y y y

T T T

z z z

c c

T T

t t

λ λ λ

ρ ρ

− +

−

− + +

∆ ∆ ∆

+ − + +

∆ ∆ ∆

+ − + =

∆ ∆ ∆

= −

∆ ∆

(22)

In the case of other planes of the cube we proceed similarly.

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, , , , , , 1,

2 2 2

, , , , , , 1

2

i j k l i j k l i j k l

i j k l i j k l

T T T

x x x

T T T

y y y

T T T

z z z

c c

T T

t t

λ λ λ

ρ ρ

− +

−

− + +

∆ ∆ ∆

+ − + +

∆ ∆ ∆

+ − + =

∆ ∆ ∆

= −

∆ ∆

(23)

Conclusions

Both this article and the previous one [3] are the introduction to the discussion

on the boundary-initial problem for n-dimensional heat flow described by the

Fourier equation.

(7)

References

[1] Biernat G., Lara-Dziembek S., Pawlak E., The determinants of the block matrices in the 3D Fourier equation, Scientific Research of the Institute of Mathematics and Computer Science 2012, 4(11), 5-10.

[2] Mochnacki B., Suchy J.S., Modelowanie i symulacja krzepnięcia odlewów, WN PWN, Warszawa 1993.

[3] Biernat G., Lara-Dziembek S., Pawlak E., The finite difference method in the 2D Fourier equation with Robin’s boundary condition, Journal of Applied Mathematics and Computational Mechanics 2013, 4(12), 5-11.

[4] Biernat G., Mazur J., Finite Difference Method in the Fourier equation with Newton’s boundary conditions direct formulas, Scientific Research of the Institute of Mathematics and Computer Science 2008, 2(7), 123-128.

(8)