INTERVAL BOUNDARY ELEMENT METHOD FOR 1D TRANSIENT DIFFUSION PROBLEM
Alicja Piasecka Belkhayat
Department for Strength of Materials and Computational Mechanics Silesian University of Technology, Gliwice, Poland
Abstract. In this paper the description of an unsteady heat transfer for one-dimensional problem is presented. It is assumed that all thermophysical parameters (specific heat, mass density and heat conduction) are given as intervals. The problem discussed has been solved using the 1st scheme of the boundary element method. The interval Gauss elimination method with the decomposition procedure has been applied to solve the obtained interval system of equations. In the final part of the paper the results of numerical computations are shown.
1. Interval arithmetic
Let us consider an interval ,~x which can be defined as a set of the form [4]
{ }
, :
x%≡ x x = ∈x R x≤ ≤x x (1)
where x and x denote the lower and the upper bounds, respectively. An interval is called thin if x=x and thick if x<x.
The sum of two intervals a%= a a, and b%= b b, can be written as ,
c% %= + =a b% a+b a+b (2) The difference is of the form
,
c% %= − =a b% a−b a−b (3) The product of the intervals is described by the formula
( ) ( )
min , , , max , , ,
c% %= ⋅ =a b% a b a b a b a b⋅ ⋅ ⋅ ⋅ a b a b a b a b⋅ ⋅ ⋅ ⋅ (4) The inversion of the interval c% can be expressed as
1 1 , 1 , 0 ,
c%= b%= b b ∉ b b (5)
The quotient of two intervals can be written as
1 , 0 ,
c% %=a b%= ⋅a% b% ∉ b b (6)
2. Formulation of the problem
Let us consider a plate domain of thickness L. The transient temperature field is described by the following linear energy equation [1, 2]
2
2
( , ) ( , )
0 : , T x t , T x t ( , )
x L c c Q x t
t x
∂ ∂
< < = λ λ +
∂ ∂ (7)
where λ λ, , c c denote the interval values of the thermal conductivity and the , volumetric specific heat, respectively, T is the temperature, Q is the heat source, x is the spatial co-ordinate and t is the time.
The equation (7) can be expressed as follows
2
2
( , ) ( , ) 1
0 : , ( , )
,
T x t T x t
x L a a Q x t
t x c c
∂ ∂
< < = +
∂ ∂ (8)
where a a, = λ λ, c c, is the interval diffusion coefficient and its lower and upper bounds can be defined according to the rules of the interval arithmetic as
( )
( )
min / , / , / , / max / , / , / , /
a c c c c
a c c c c
= λ λ λ λ
= λ λ λ λ (9)
The above equation is supplemented by the following boundary-initial conditions
1
2
0
0: ( , ) ( , )
: ( , ) ( , )
0: ( , ) ( )
b
b
x T x t T x t
x L T x t T x t
t T x t T x
= =
= =
= =
(10)
where Tb1 and Tb2 are the known boundary temperatures and T0 is the initial temperature.
3. Interval boundary integral equation Let us introduce the time grid
0 1 2 1
0= < < < <t t t K tf− <tf < <K tF < ∞ (11)
with a constant time step ∆t=tf −tf−1.
Taking into account the criterion of the weighted residual method (WRM) the following interval boundary integral equation is obtained [1-3]
1
2
2 0
( , ) ( , ) 1
, ( , ) ( , , , ) d d 0
,
f
f
t L
f
t
T x t T x t
a a Q x t T x t t x t
x t c c
−
∂ −∂ + ∗ ξ =
∂ ∂
∫ ∫
% (12)where T%∗( , ,ξ x tf, )t is the interval fundamental solution, ξ is the point where the concentrated heat source is applied. The interval function TT%( , ,∗( , ,ξ x tx tff, ), )t is t expressed as
1 ( )2
( , , , ) , exp
4 , ( )
2 , ( )
f
f f
T x t t T T x
a a t t a a t t
∗ ξ = ∗ ∗ = − − ξ
−
π −
% (13)
where
2 2
2 2
2
1 ( ) 1 ( )
min exp , exp ,
4 ( ) 4 ( )
2 ( ) 2 ( )
1 ( ) 1 ( )
exp , exp
4 ( ) 4 ( )
2 ( ) 2 ( )
1 ( ) 1
max exp ,
4 ( )
2 ( ) 2
f f
f f
f f
f f
f f
x x
T a t t a t t a t t a t t
x x
a t t a t t
a t t a t t
T x
a t t
a t t a
∗
∗
− ξ − ξ
= − −
− −
π − π −
− ξ − ξ
− −
− −
π − π −
− ξ
= −
π − − π
2
2 2
( )
exp ,
4 ( )
( )
1 ( ) 1 ( )
exp , exp
4 ( ) 4 ( )
2 ( ) 2 ( )
f f
f f
f f
x a t t t t
x x
a t t a t t
a t t a t t
− ξ
−
− −
− ξ − ξ
− −
− −
π − π −
(14)
The fundamental solution TT%( , ,∗( , ,ξ x tx tff, ), )t fulfils the properties t
2
2
( , , , ) ( , , , )
, ( ) ( )
lim ( , , , ) ( )
f
f f
f
f t t
T x t t T x t t
a a x t t
x t
T x t t x
∗ ∗
∗
→
∂ ξ +∂ ξ = −δ − ξ δ −
∂ ∂
ξ = δ − ξ
% %
% %
% %
(15)
where ( )δ ⋅% is the interval Dirac’s delta function.
The interval heat flux corresponding to the fundamental solution results from the following formula
2
2
3 / 2
( , , , )
( , , , ) , ,
1 2( ) ( )
, exp
4 , ( ) 4 , ( )
2 , ( )
, ( ) ( )
exp 4 , ( )
4 , ( )
f f
f f
f
f f
T x t t
q x t t q q
x
x x
a a t t a a t t
a a t t
x x
a a t t a a t t
∗ ξ = ∗ ∗ = − λ λ ∂ ∗ ξ =
∂
− ξ − ξ
λ λ − =
− −
π −
λ λ − ξ − − ξ
−
π −
% %
(16)
and the lower and the upper bounds of ~q*(ξ,x,tf,t) we can define as
2
3 / 2
2
3 / 2
, ( ) ( )
min exp
4 , ( )
4 , ( )
, ( ) ( )
max exp
4 , ( )
4 , ( )
f f
f f
x x
q a a t t a a t t
x x
q a a t t a a t t
∗
∗
λ λ − ξ − ξ
= −
−
π −
λ λ − ξ − ξ
= −
−
π −
(17)
In this way we obtain the following criterion
1 1
1
2
2
0 0
0
( , ) ( , )
, ( , , , ) d d ( , , , ) d d
1 ( , ) ( , , , ) d d 0
,
f f
f f
f
f
t L t L
f f
t t
t L
f
t
T x t T x t
a a T x t t x t T x t t x t
x t
Q x t T x t t x t c c
− −
−
∗ ∗
∗
∂ ξ − ∂ ξ +
∂ ∂
ξ =
∫ ∫ ∫ ∫
∫ ∫
% %
%
(18)
Integrating twice by parts the first component of the equation (18) and integrating by parts with respect to time the second component of this equation and then using the properties of the fundamental solution finally we obtain the interval boundary integral equation
1
1
1
0
1 1
0 0
0
( , ) 1 ( , , , ) ( , ) d
,
1 ( , , , ) ( , ) d ( , , , ) ( , ) d
,
1 ( , ) ( , , , ) d d
,
f
f
f
f
f
f
t x L
f f
t x
t x L L
f f f f
t x
t L
f
t
T t T x t t q x t t
c c
q x t t T x t t T x t t T x t x
c c
Q x t T x t t x t c c
−
−
−
=
∗
=
=
∗ ∗ − −
=
∗
ξ + ξ =
ξ + ξ +
ξ
∫
∫ ∫
∫ ∫
% % %
% % %
%
%
(19)
where q%= − λ λ ∂ ∂, T% x and ( ,T%ξ tf) is the interval function of the temperature.
4. Numerical approximation
Let us consider the constant elements with respect to time, which can be defined as follows
1 ( , ) ( , )
, :
( , ) ( , )
f
f f
f
T x t T x t
t t t
q x t q x t
− =
∈
=
% %
% % (20)
The equation (20) can be written as
1
1
1
0
1 1
0 0
1
0
( , ) 1 ( , ) ( , , , ) d
,
1 ( , ) ( , , , ) d ( , , , ) ( , ) d
,
1 ( , ) ( , , , ) d d
,
f
f
f
f
f
f
t x L
f f f
t x
t x L L
f f f f f
t x
L t
f f
t
T t q x t T x t t t
c c
T x t q x t t t T x t t T x t x
c c
Q x t T x t t t x
c c
−
−
−
=
∗
=
=
∗ ∗ − −
=
− ∗
ξ + ξ =
ξ + ξ +
ξ
∫
∫ ∫
∫ ∫
% % %
% % % %
%
(21)
At first, we calculate the integrals
1
1 sgn( )
( , ) ( , , , ) d erfc
, 2 2 ,
f
f
t
f
t
x x
h x q x t t t
c c − a a t
∗ − ξ − ξ
ξ = ξ =
∆
∫
% % (22)
where
sgn( ) sgn( )
min erfc , erfc
2 2 2 2
sgn( ) sgn( )
max erfc , erfc
2 2 2 2
x x
x x
h
a t a t
x x
x x
h a t a t
− ξ − ξ − ξ − ξ
= ∆ ∆
− ξ − ξ − ξ − ξ
= ∆ ∆
(23)
and
1
2
( , ) 1 ( , , , ) d
,
( )
exp erfc
2 , 2 ,
4 , ( )
, ,
f
f
t
f
t
f
g x T x t t t
c c
x x
t x
a a t a a t t
c c
−
ξ = ∗ ξ =
− ξ − ξ
∆ − − ξ −
− λ λ ∆
λ λ
∫
%%
(24)
where
2
2
( )
min exp
4 , ( )
, ,
2 , erfc 2 ,
( )
max exp
4 , ( )
, ,
2 , erfc 2 ,
f
f
t x
g
a a t t c c
x x
a a t
t x
g
a a t t c c
x x
a a t
∆ − ξ
= − −
−
λ λ
− ξ − ξ
λ λ ∆
∆ − ξ
= − −
−
λ λ
− ξ − ξ
λ λ ∆
(25)
After introducing the obtained results (22) and (24) to the equation (21) one obtains
0 0
1 1
( , ) ( , ) ( , ) ( , ) ( , )
( , ) ( , )
x L x L
f f f
x x
f f
T t g x q x t h x T x t
P t Z t
= =
= =
− −
ξ + ξ = ξ +
ξ + ξ
% % % % %
% % (26)
or
1 1
( , ) ( , ) ( , ) ( , 0) (0, )
( , ) ( , ) ( , 0) (0, ) ( , ) ( , )
f f f
f f f f
T t g L q L t g q t
h L T L t h T t P t − Z t −
ξ + ξ − ξ =
ξ − ξ + ξ + ξ
% % % % %
% % % % % % (27)
where
1 1 1
0
2
1
0
( , ) ( , , , ) ( , ) d
1 ( )
exp ( , ) d
4 ,
2 ,
L
f f f f
L
f
P t T x t t T x t x
x T x t x
a a t
a a t
− ∗ − −
−
+ −
+ −
ξ = ξ =
− ξ
−
∆
π ∆
∫
∫
% % %
%
(28)
with
2
1
0
2
1
0
1 ( )
min exp ( , ) d
4 ,
2 ,
1 ( )
max exp ( , ) d
4 ,
2 ,
L
f
L
f
P x T x t x
a a t
a a t
P x T x t x
a a t
a a t
−
+ −
+ −
−
+ −
+ −
− ξ
= −
∆
π ∆
− ξ
= −
∆
π ∆
∫
∫
(29)
and
1 1
0
( , ) ( , ) ( , ) d
L
f f
Z%ξ t − =
∫
Q x t − g%ξ x x (30) with1 1
0 0
1 1
0 0
min ( , ) ( , ) d , ( , ) ( , ) d
max ( , ) ( , ) d , ( , ) ( , ) d
L L
f f
L L
f f
Z Q x t g x x Q x t g x x
Z Q x t g x x Q x t g x x
− −
− −
= ξ ξ
= ξ ξ
∫ ∫
∫ ∫
(31)
For ξ→0+ and ξ→L− one obtains a system of two equations. Taking into account the defined boundary conditions (10) we get the following system of equations
( ) ( )
( ) ( )
1
2
11 11 12 12
21 21 2 2 2 2
11 11 12 12
21 21 2 2 2 2
1 1 1 1
1 1
0, , 0,
, ,
, , , , ,
(0, )
, ,
( , )
, ,
(0, ), (0, ) (0, ), (0, )
( , ), ( , )
f f
f f
f b
f b
f f f f
f f
q t q t
G G G G
G G G G q L t q L t
T t
H H H H
T L t
H H H H
P t P t Z t Z t
P L t P L t Z
− − − −
− −
=
⋅ +
+
( ,L tf−1),Z L t( , f−1)
(32)
where the elements of interval matrices G% and H% are defined as
11 12
21 2 2
(0, 0) (0, )
( , 0) ( , )
G g G g L
G g L G g L L
= − =
= − = −
% % % %
% % % % (33)
and
11 12
21 2 2
0.5 (0, )
( , 0) 0.5
H H h L
H h L H
= − =
= − = −
%
% %
%
% % (34)
The interval Gauss elimination method with the decomposition procedure [4, 5] has been used to solve the system of equations (32). After determining the
’missing’ boundary values the functions ( ,T%ξ tf) at internal nodes of the domain considered are calculated using the formula
2 1
1 1
( , ) ( , ) ( , ) ( , 0) (0, )
( , ) ( , ) ( , 0) (0, ) ( , ) ( , )
f f f
b b
f f f f
T t h L T L t h T t
g L q L t g q t P t − Z t −
ξ = ξ − ξ −
ξ + ξ + ξ + ξ
% %
%
% %
% % % %
(35)
5. Results of computations
Let us consider a plate domain of dimension L = 0.1 m. The following input data have been introduced: λ = 35 W/(m ⋅ K), λ = λ – 0.02 ⋅ λ W/(m ⋅ K),
λ = λ + 0.02 ⋅ λ W/(m ⋅ K), c = 7500 ⋅ 650 J/(m3 ⋅ K), c = c – 0.02 ⋅ c J/(m3 ⋅ K), c = c + 0.02 ⋅ c J/(m3 ⋅ K), Q = 10 000 W/m3, initial temperature T0 = 500°C, boundary values Tb1 = 100°C, Tb2 = 200°C.
The domain considered has been divided into 20 linear elements, the time step
∆t = 0.5 s.
100 200 300 400 500
0,00 0,02 0,04 0,06 0,08 0,10 x[m]
T [oC]
Tem_L Tem_R
25s 50s 75s
Fig. 1. Temperature distribution
300 350 400 450 500 550
0 5 10 15 20 25 t[s]
T [oC]
Tem_L Tem_R
2
1 3
Fig. 2. Cooling curves
In Figure 1 the distribution of temperature for the time 25, 50 and 75 s is shown (Tem_L, Tem_R denote the first and the second endpoints of the interval).
Figure 2 illustrates the cooling curves obtained at the points x1 = 0.02 m, x2 = 0.05 m and x3 = 0.08 m in the domain considered.
References
[1] Majchrzak E., Boundary element method in heat transfer, Publ. of the Technological University of Czestochowa, Czestochowa 2001 (in Polish).
[2] Majchrzak E., Mochnacki B., The BEM application for numerical solution of non-steady and non-linear thermal diffusion problems, Computer Assisted Mechanics and Engineering Sciences 1996, 3, 4, 327-346.
[3] Mochnacki B., Suchy J.S., Numerical methods in computations of foundry processes, PFTA, Cracow 1995.
[4] Neumaier A., Interval methods for system of equations, Cambridge University Press, Cambridge, New York, Port Chester, Melbourne, Sydney 1990.
[5] Piasecka-Belkhayat A., Application of the interval methods for solving linear thermal diffusion problems, Scientific Research of the Institute of Mathematics and Computer Science 2005, 1(4), 204-210.
