The n-dimensional Fourier equation with the Robin’s boundary conditionusing the finite difference method

(1)

THE n-DIMENSIONAL 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

Częstochowa, Poland

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

Abstract. The article is a continuation of the considerations on the three-dimensional Fourier equation supplemented by the third boundary condition (the Robin’s condition).

The paper is based on the Finite Difference Method.

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

Introduction

The subject of our discussion is the heat flow in the n-dimensional area described by the Fourier equation with Robin’s boundary condition. In this paper we present the exact solution to this problem.

Solution of the problem

We consider the n-dimensional Fourier equation [1, 2]

( ) ( ) ( )

( )

1 2 2 2

1 2 1 2 1 2

2 2 2

2 1 2

, , ..., , , , ..., , , , ..., ,

...

, , ..., , λ

ρ

 ∂ ∂ ∂ 

 + + +  =

 ∂ ∂ ∂ 

 

= ∂

∂

n

n n n

n

T x x x t T x x x t T x x x t

x x x

T x x x t

c

t

(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], ₁ , ₂ , ...,

n

x x x denote the geometrical co-ordinates and t is time [s].

Assuming the following difference quotients we obtain the differential approxima-

tion of the second derivatives appearing in equation (1)

(2)

( )

1 2 1 2 1 2

1 1 2 1 2 1 2

2 1 2 1 2 1 2

2 1, ,..., , , ,..., , 1, ,..., ,

1 1

2 2

1 2 , 1,..., , , ,..., , , 1,..., ,

2 2

2 2 , ,..., 1, , ,..., , , ,.

2

2 , 1 1

2 − +

− +

−

− +

∆ = ≤ ≤ −

∆ ∆

− +

∆ = ≤ ≤ −

∆ ∆

− +

∆ =

∆ M

n n n

n n

i i i l i i i l i i i l

i i i l i i i l i i

n

T T T

T

i m

x x

T T T

T i m

x x

T T T

T

x ( )

.., 1,

2 ⁺ , 1 ≤ ≤ − 1

∆

i n l

n n

n

i m

x

(2)

and the approximation of the first derivative of the time

1 2 , ,..., , 1 2 , ,..., , 1

− , 1

∆ −

= ≤ ≤

∆ ∆

n n

i i i l i i i l

T T

T l q

t t (3)

The initial condition has the form

1 2 , ,..., ,

, 0

= ≤ ≤

i i i l n 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-4] is presented

( )

2 0, 2 ,..., ,

,1,..., , 1

, 0,..., , 1

, ,...,1,

1 2

, ,..., 0,

1 2

1 1, ,..., ,

1

2 2 2 2

2 2 2

2

2 ...

λ α

λ ρ

 −

 = −

 ∆

 −

  = −

 ∆

 

 −

 = −

 ∆

   ∂ ∂ ∂  ∆

  + + +  =

 ∂ ∂ ∂  ∆

 

 M

n

i i n l

i i l

n

i i l env

T env

n

T T

x

T T

x

T T

x

T T T T

c

x x x t



(5)

where α is the heat transfer coefficient and T _env is the ambient temperature of environment.

Now we consider the systems of equations corresponding to the exemplary vertices,

edges and lateral planes of our n-cube.

(3)

In the first vertex of the n-cube (for ₁ , ₂ , ...,

i i i n ) we have

( )

( ) ( ) ( )

( ) ( )

1,0,..., 0,

0, 0,...,0, 1

0,1,0, ...,

0,0,...,0, 2

0, 0,...,1,

0,0,...,0,

0,0,0,..., 1,0,0,...0,

2 2 2

1 1 1

0,0,0,...,

2 2

2

2 λ α

λ α

λ λ λ

λ λ

− = −

∆

− = −

∆

− = −

∆

− + +

∆ ∆ ∆

+ −

∆ ∆

M

l env

n

env l l

env

T T

x

T T

x

T T

x

T T T

x x x

T T

x x ( )

( ) ⁽ ⁾ ⁽ ⁾

0, 1, 0,...0, 2

2 0,0, 0,...0, 0,0, 0,...0, 1

0,0,0,..., 0, 0,0,...1,

2 2 2

...

2 λ

λ λ λ

ρ ⁻

+ + +

∆

+ − + = −

∆ ∆ ∆

∆

l l

env l l

n

T x

T T

T T T c

z z t

x

(6)

Thus, the system of equations is as follows:

( ) ( )

0,0,...,0, 1,0,...,0,

1 1

0,0,...,0, 0,1,...0,

2 2

0,0,...,0, 0,0,...,1,

0,0,..., 0, 1,

2 2

1 1

2 2

2 λ λ

α α

λ λ

α α

λ λ

α α

λ ρ λ

=

 

− ∆ =    − ∆   

 

− ∆ =    − ∆   

 

− ∆ =    − ∆   

 

 +  −

 ∆ ∆  ∆

 

 

∑

M

l l env

n n

n

l

p p

T T T

x x

T T T

x x

T T T

x x

c T T

x t x ( ) ( )

( )

0,..., 0, 2 0,1,...,0, 2 0,0,...,1,

2 0,0,...,0, 1 2

1 λ ... λ

λ ρ

−

=

 

 

 

 

 − − − =

 ∆ ∆

 

= +

 ∆ ∆



∑

l l l

n n

env l

p p

T T

x x

T c T

x t

(7)

(4)

The main matrix system above has the form

( ) ( ) ( ) ( )

1

2 2 2 2 2

1 1 2

0 0

2 0 0

2 2

α λ

λ ρ λ λ λ

=

 − 

 ∆ 

 

 − 

 

 ∆ 

 

=  

−

 

 ∆ 

 

 + − − − 

 ∆ ∆ ∆ ∆ ∆ 

 

∑

L

M M M O M

L

n n

p p n

x

x A

x c

x t x x x

(8)

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

1

2 0 2

0 2

det 0

α λ

 − >

 ∆

 

− >

 ∆

 

 

 

− >

 ∆

  >

 M

n

x

x A

(9)

So

1

2 2 2

2 λ λ

α α

λ λ

α α

λ λ

α α

 < ∆ <

 

 < ∆ <

 

 

 < ∆ <



M

n

x x

x

(10)

(5)

The determinant of the matrix A can be written in the following form:

( )

3

2 1 1

1 det 1 4 2

λ λ ρ

α

= =

=

 

 

= ⋅ + −

 ∆ ∆ ∆ 

∆    

∑ ∑

∏

n n

n

p p p

c A

t x

x x

(11)

Now, we calculate the determinant corresponding to the variable T 0, 0,…, 0, l

( )

0,0,..., 0,

3 0,0,..., 0, 1 2

1 1

1 det 1

4 2

λ λ ρ

α ₋

= =

=

   

   

= ∆ ⋅       ∆ − ∆    + ∆   

∑ ∑

l ∏

n n

T n env l

p p p

c

A T T

x t

x x

(12)

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

( )

0,0,...,0, 1 2

1 1

0,0,...,0,

2 1 1

1 2 1 2

λ ρ

α

λ ρ

α

−

= =

 

 −  +

 ∆  ∆

 ∆ 

 

=

+ −

∆ ∆

∆

∑ ∑

n n

env l

p p p

p

l n n

p p p

p

T c T

x t

x T

c

t x

x

(13)

Assuming the following markings

( ) ² ¹

, 1,..., , 1 ,

2 λ ρ

α

=

= = = =

∆ ∆

∆

∑

n p

p p

p

a p n b d c

x t

x

(14)

we get

0,0,...,0, 1 1

0,0,..., 0,

1 −

=

 

− +

 

 

=

− +

∑

n

p env l

p

l n

p p

a b T d T

T

a b d

(15)

or in another notation

1 0,0,...,0, 0,0,...,0, 1

1 1

=

−

= =

−

= +

− + − +

∑

∑ ∑

n p p

l n env n l

p p

a b

T T d T

a b d a b d

(16)

Similarly, we proceed for the remaining vertices.

(6)

Subsequently we repeat the procedure for each edge .

And so, on the edge (i 1 ,0,…,0), where 1 ≤ i ₁ ≤ m ₁ − 1, i ₂ = 0, ..., i _n = 0 we receive

( )

( ) ( ) ( )

1 1

1

1 1 1 1

1,0,...,0, 1, 0,..., 0,

,0,,...,0, 1

,1,0,..., 0,

, 0,..., 0, 2

, 0,..., 1,

, 0,..., 0,

1, 0,..., 0, ,0,...,0, 1, 0,

2 2 2

1 1 1

2

2 λ α

λ α

λ λ λ

+ −

− +

− = −

∆

− = −

∆

− = −

∆

− +

∆ ∆ ∆

M

i l i l

i l env

n

i l i l i

T T

x

T T

x

T T

x

T T T

x x x

( ) ( ) ( )

1 1

...,0,

,0,..., 0, ,1,...,0,

2 2 2

,0,...,0, ,0,...,0, 1

,0,...,0, ,0,..., 1,

2 2 2

2 ...

2 λ λ λ

λ λ λ

ρ ⁻

+

+ − + + +

∆ ∆ ∆

+ − + = −

∆ ∆ ∆ ∆

l

env i l i l

i l i l

env i l i l

n n n

T T T

x x x

T T

T T T c

x x x t

(17)

Thus, the respective determinants take the following forms:

¹ ( ) ² ¹

1 det 1 4 2

λ λ ρ

α

= =

=

 

 

= ∆ ⋅    ∆ + ∆ − ∆   

∑ ∑

∏

n n n

n

p p p

c A

t x

x x

(18)

( ) ( )

, 0,..., 0,

1,0,...,0, ,0,...,0, 1

2 2

2 1

1 1

det

1 4 2

λ λ λ ρ

α − −

= =

=

   

   

= ∆ ⋅       ∆ − ∆    + ∆ + ∆   

∑ ∑

∏

i l

T

n n n

env i l i l

n

p p p

A

c

T T T

x t

x x

x

(19)

Ultimately, the temperature on the edge (i, 0, 0) 1 ≤ i ₁ ≤ m ₁ − 1, i ₂ = 0, ..., i _n = 0

at the moment l is given by the formula

(7)

( ) ( )

( )

1,0,..., 0, ,0,..., 0, 1

2 2

2 1

1 ,0,...,0,

2 1 1

1 2 1

2 λ λ ρ

α

λ ρ

α

− −

= =

 

 −  + +

 ∆  ∆

 

=

− +

∆ ∆

∆

∑ ∑

n n

env i l i l

p p p

p

i l n n

p p p

p

c

T T T

x t

x x

T

c

x t

x

(20)

Using designations (14), we obtain

1 2

,0,...,0, 1,0,...,, ,0,...,0, 1

1 1 1

=

− −

= = =

−

= + +

− + − + − +

∑

∑ ∑ ∑

n p p

i l n env n i l n i l

p p p

a b

a d

T T T T

a b d a b d a b d

(21)

Analogously we calculate the temperature for the other edge points.

On the plane (i 1 , i 2 , …, 0), 1 ≤ i ₁ ≤ m ₁ − 1, 1 ≤ i ₂ ≤ m ₂ − 1,

3 = 0,..., _n = 0

i i we receive

( ) ( ) ( )

1 2 1 2 1 2

1 2 1 2

1, , 0,...,0, , ,0,..., 0, 1, ,0,...,0,

2 2 2

1 1 1

, 1,0,..., 0, , ,0,...,0, , 1, 0,...,0,

2 2 2

, , 0,...,0, , ,0,...,1

2 2 2

2 2 ...

2 λ λ λ

λ λ λ

− +

− + +

∆ ∆ ∆

+ − + + +

∆ ∆ ∆

+ − +

∆ ∆ ∆

i i l i i l i i l

env i i l i i

n n n

T T T

x x x

T T T

x x x

T T T

x x x , ( 1 2 , , 0,...,0, 1 2 , ,0,...,0, 1 )

ρ

= − −

l ∆ i i l i i l

c 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 2 1 2 1 2

1, ,..., , , ,..., , 1, ,..., ,

2 2 2

1 1 1

, 1,..., , , ,..., , , 1,..., ,

2 2 2

, ,..., 1, , ,..., , , ,.

2 2 2

2 2 ...

2 λ λ λ

λ λ λ

− +

−

− + +

∆ ∆ ∆

+ − + + +

∆ ∆ ∆

+ − +

∆ ∆ ∆

n n n

n n

i i i l i i i l i i i l

i i i l i i i l i i

n n n

T T T

x x x

T T T

x x x

T T T

x x x

( 1 2 1 2 )

.., 1,

, ,..., , , ,..., , 1

ρ

+

−

=

= −

∆

n

n n

i l

i i i l i i i l

c T T

t

(23)

(8)

Conclusion

This article is the continuation of the discussion on the boundary-initial problem for the heat flow described by the Fourier equation [3, 4].

References

[1] Biernat G., Lara-Dziembek S., Pawlak E., The determinants of the block band matrices based on the n-dimensional Fourier equation. Part 1, Journal of Applied Mathematics and Computational Mechanics 2013, 3(12), 5-13.

[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-12.

[4] Biernat G., Lara-Dziembek S., Pawlak E., The 3D Fourier equation with the Robin’s boundary condition using the Finite Difference Method, Journal of Applied Mathematics and Computa- tional Mechanics 2014, 1(13), 5-12.

[5] 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.

The n-dimensional Fourier equation with the Robin’s boundary conditionusing the finite difference method