Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

49, 1, pp. 265-276, Warsaw 2011

INTERVAL BOUNDARY ELEMENT METHOD FOR

TRANSIENT DIFFUSION PROBLEM IN TWO-LAYERED

DOMAIN

Alicja Piasecka-Belkhayat

Department of Strength of Materials and Computational Mechanics, Silesian University of

Technology, Gliwice, Poland; e-mail: alicja.piasecka@polsl.pl

In the paper, the description of an unsteady heat transfer for one-
dimensional problemproceeding in a two-layereddomain is presented. It
is assumed that all thermophysical parameters appearing in the mathe-
matical descriptionof the problemanalysedare givenasdirected interval
values. The problem discussed has been solved using the 1st scheme of
the interval boundary element method. The interval Gauss elimination
method has been applied to solve the obtained interval system of equ-
ations. In the final part of the paper, results of numerical computations
are shown.

Key words: directed interval arithmetic, interval boundary element me-
thod, heat transfer

1. Introduction

Heat transfer problems are usually solved using equations with deterministic
parameters.However, inmost cases of the engineering practice, values of these
parameters cannotbedefinedwithahighprecision, and in suchcases it ismuch
more convenient to define these parameters as intervals.

In the available literature we can find examples of papers using the inte-
rval arithmetic (Neumaier, 1990) and the theory of fuzzy sets (Zadeh, 1965)
allowing one to solve problems taking into account ”uncertainties” in thema-
thematical model. We can also find papers dedicated to interval boundary
element method (Burczyński and Skrzypczyk, 1997) and the interval finite
element method (Muhanna et al., 2005). However, most of these papers are
related to boundary problems, andwe can hardly find any examples of papers
dealing with boundary-initial problems.



266 A. Piasecka-Belkhayat

In this paper, the Interval Boundary ElementMethod (IBEM) for solving
non-steady heat transfer problems with directed interval thermophysical pa-
rameters of both sub-domains has been presented with the approach of the
directed interval arithmetic (Markov, 1995, Popova, 2001). This assumption
is closer to real physical conditions of the process considered because it is
difficult to estimate the thermophysical parameters appearing in the mathe-
matical model. The main advantage of the directed interval arithmetic upon
the classical interval arithmetic is that the obtained temperature intervals are
much narrower.
In theory as well as in practice, it is valuable to develop the Interval Bo-

undary ElementMethod (IBEM).

2. Directed interval arithmetic

Let us consider a directed interval ã which can be defined as a set D of all
directed pairs of real numbers defined as follows (Kużelewski, 2008; Markov,
1995; Popova, 2001)

ã= 〈a−,a+〉 := {ã∈D| a−,a+ ∈R} (2.1)
where a− and a+ denote the beginning and the end of the interval, respecti-
vely.
The left or the right endpoint of the interval ã can be denoted as as,

s∈ {+,−}, where s is a binary variable. This variable can be expressed as a
product of two binary variables and is defined as

++=−−=+
+−=−+=− (2.2)

An interval is called proper if a− ¬ a+, improper if a− ­ a+ and dege-
nerate if a− = a+. The set of all directed interval numbers can be written
as D = P ∪I, where P denotes the set of all directed proper intervals and
I denotes the set of all improper intervals.
Additionally, a subset Z=ZP ∪ZI ∈D should be defined, where
ZP = {ã∈P | a− ¬ 0¬ a+} ZI = {ã∈ I| a+ ¬ 0¬ a−} (2.3)
For directed interval numbers two binary variables are defined. The first

of them is the direction variable

τ(ã)=

{
+ if a− ¬ a+

− if a− >a+
(2.4)



Interval BEM for transient diffusion problem... 267

and the other is the sign variable

σ(ã)=

{
+ if a− > 0, a+ > 0

− if a− < 0, a+ < 0
ã∈D\Z (2.5)

For the zero argument σ(〈0,0〉) = σ(0) = +, for all intervals including the
zero element ã∈Z, σ(ã) is not defined.
The sum of two directed intervals ã = 〈a−,a+〉 and b̃ = 〈b−,b+〉 can be

written as
ã+ b̃= 〈a−+ b−,a++ b+〉 ã, b̃∈D (2.6)

The difference is of the form

ã− b̃= 〈a−− b+,a+− b−〉 ã, b̃∈D (2.7)

The product of the directed intervals is described by the formula

ã · b̃=






〈
a−σ(̃b) · b−σ(̃a),aσ(̃b) · bσ(̃a)

〉
ã, b̃∈D\Z

〈
aσ(̃a)τ (̃b) · b−σ(̃a),aσ(̃a)τ (̃b) · bσ(̃a)

〉
ã∈D\Z, b̃∈Z

〈
a−σ(̃b) · bσ(̃b)τ (̃a),aσ(̃b) · bσ(̃b)τ (̃a)

〉
ã∈Z, b̃∈D\Z

〈min(a− · b+,a+ · b−),max(a− · b−,a+ · b+)〉 ã, b̃∈ZP
〈max(a− · b−,a+ · b+),min(a− · b+,a+ · b−)〉 ã, b̃∈ZI
0 (ã∈ZP , b̃∈ZI)∪ (ã∈ZI, b̃∈ZP)

(2.8)
The quotient of two directed intervals can be written as

ã/b̃=






〈
a−σ(̃b)/bσ(̃a),aσ(̃b)/b−σ(̃a)

〉
ã, b̃∈D\Z

〈
a−σ(̃b)/b−σ(̃b)τ (̃a),aσ(̃b)/b−σ(̃b)τ (̃a)

〉
ã∈Z, b̃∈D\Z

(2.9)

In the directed interval arithmetic, two extra operators are defined – inversion
of summation

−
D
ã= 〈−a−,−a+〉 ã∈D (2.10)

and inversion of multiplication

1/
D
ã= 〈1/a−,1/a+〉 ã∈D\Z (2.11)

So, two additional mathematical operations can be defined as follows

ã−
D
b̃= 〈a−− b−,a+− b+〉 ã, b̃∈D (2.12)



268 A. Piasecka-Belkhayat

and

ã/
D
b̃=





〈
a−σ(̃b)/b−σ(̃a),aσ(̃b)/bσ(̃a)

〉
ã, b̃∈D\Z

〈
a−σ(̃b)/bσ(̃b),aσ(̃b)/bσ(̃b)

〉
ã∈Z, b̃∈D\Z

(2.13)

Now, it is possible to obtain the number zero by subtraction of two identical
intervals ã−

D
ã=0and the number one as the result of the division ã/

D
ã=1,

which was impossible when applying classical interval arithmetic (Neumayer,
1990).

3. Heat transfer model in two-layered domain

Let us consider a two-layered domain of dimension L = L1 +L2. The heat
conduction process in the first sub-domain is describedby the following energy
equation (Mochnacki and Suchy, 1995; Majchrzak, 2001)

x∈ (0,L1) : 〈c−1 ,c
+
1 〉
∂T1(x,t)

∂t
= 〈λ−1 ,λ

+
1 〉
∂2T1(x,t)

∂x2
(3.1)

where 〈c−1 ,c
+
1 〉 is the directed interval volumetric specific heat for the first

sub-domain, 〈λ−1 ,λ
+
1 〉 is the directed interval thermal conductivity, T1, x, t

denote temperature, spatial co-ordinate and time, respectively.
Equation (3.1) can be expressed as follows

x∈ (0,L1) :
∂T1(x,t)

∂t
= 〈a−1 ,a

+
1 〉
∂2T1(x,t)

∂x2
(3.2)

where 〈a−1 ,a
+
1 〉 = 〈λ

−

1 ,λ
+
1 〉/〈c

−

1 ,c
+
1 〉 is the directed interval diffusion coeffi-

cient, and its beginning and end can be defined according to the rules of the
directed interval arithmetic (Markov, 1995).
Taking into account the assumption that λ̃1, c̃1 ∈ D\Z. one obtains the

following formula

λ̃1/D c̃1 =
〈
λ
−σ(̃c1)
1 /c

−σ(λ̃1)
1 ,λ

σ(̃c1)
1 /c

σ(λ̃1)
1

〉
(3.3)

For example, for the interval coefficients λ̃1 = 〈34,35〉 and c̃1 =
= 〈4900000,5400000〉 the sign variables are σ(λ̃1) = +, σ(c̃1) = +, so the
quotient of λ̃1 and c̃1 can be calculated according to the formula

λ̃1/D c̃1 = 〈λ
−+
1 /c

−+
1 ,λ

+
1 /c
+
1 〉= 〈λ

−

1 /c
−

1 ,λ
+
1 /c
+
1 〉 (3.4)



Interval BEM for transient diffusion problem... 269

and the directed interval diffusion coefficient ã1 is computed as follows

ã1 =
λ̃1
c̃1
=

〈34,36〉
〈4900000,5400000〉

= 〈34,36〉/
D
〈4900000,5400000〉 =

(3.5)

=
〈 34
4900000

,
36

5400000

〉
≈〈0.0000069,0.0000066〉

As a result, the interval obtained is improper.

The temperature field in the other sub-domain is determined by the energy
equation

x∈ (L1,L2) : 〈c−2 ,c
+
2 〉
∂T2(x,t)

∂t
= 〈λ−2 ,λ

+
2 〉
∂2T2(x,t)

∂x2
(3.6)

where 〈c−2 ,c
+
2 〉, 〈λ

−

2 ,λ
+
2 〉 are the directed interval values of volumetric specific

heat and thermal conductivity, respectively, and T2 denotes temperature for
the second sub-domain.

The above equation, (3.6), can be transformed as follows

x∈ (L1,L2) :
∂T2(x,t)

∂t
= 〈a−2 ,a

+
2 〉
∂2T2(x,t)

∂x2
(3.7)

where 〈a−2 ,a
+
2 〉= 〈λ

−

2 ,λ
+
2 〉/〈c

−

2 ,c
+
2 〉 is thedirected interval diffusion coefficient

for the second layer.

Equations (3.2) and (3.7) must be supplemented by the boundary-initial
conditions of the following form

x=0 : q̃(0, t) =−〈λ−1 ,λ
+
1 〉
∂T1(x,t)

∂x
= q̃L

x=L2 : q̃(L2, t)=−〈λ−2 ,λ
+
2 〉
∂T2(x,t)

∂x
= q̃R

t=0 : T1(x,0)=T10(x) T2(x,0)=T20(x)

(3.8)

and the continuity condition on the contact surface between two layers

x=L1 :





−〈λ−1 ,λ

+
1 〉
∂T1(x,t)

∂x
=−〈λ−2 ,λ

+
2 〉
∂T2(x,t)

∂x

T1(x,t)=T2(x,t)

(3.9)

where q̃L, q̃R are the given interval boundary heat fluxes, T10 and T20 are the
initial temperatures for the first and second layer, respectively.



270 A. Piasecka-Belkhayat

4. Interval boundary element method

In this paper, the 1st scheme of the interval boundary elementmethod is used
(Brebbia et al., 1984; Majchrzak, 1998, 2001). At first, the time grid must be
introduced

0= t0 <t1 <t2 < ... < tf−1 <tf < ... < tF <∞ (4.1)

with a constant time step ∆t= tf − tf−1.
Let us consider the constant elements with respect to time (Majchrzak,

2001, Majchrzak andMochnacki, 1996). The boundary integral equation cor-
responding to the transition tf−1 → tf for the first layer is following

T̃1(ξ,t
f)+

[
1

c̃1
q̃1(x,t

f)

tf∫

tf−1

T̃∗1(ξ,x,t
f, t) dt

]x=L1

x=0

=

(4.2)

=

[
1

c̃1
T̃1(x,t

f)

tf∫

tf−1

q̃∗1(ξ,x,t
f , t) dt

]x=L1

x=0

+

L1∫

0

T̃∗1(ξ,x,t
f , tf−1)T̃1(x,t

f−1) dx

where ξ is the observation point, q̃1(x,t
f) is the directed interval heat flux.

The directed interval fundamental solution T∗1 (ξ,x,t
f, t) is of the following

form (Brebbia et al., 1984; Majchrzak, 2001)

T̃∗1 (ξ,x,t
f, t)=

1

2
√
πã1(t

f − t)
exp
[
− (x− ξ)

2

4ã1(t
f − t)

]
(4.3)

and should be calculated according to the formula

〈T∗−,T∗+〉=exp
{[
− (x− ξ)

2

4(tf − t)
]
/
D
〈a−1 ,a

+
1 〉
}
/
D

[
2
√
π〈a−1 ,a

+
1 〉(tf − t)

]
(4.4)

Because σ{−(x− ξ)2/[4(tf − t)]}=−, σ(ã1) =+, the interval fundamental
solution can be calculated as follows (see Eq. (2,13)) (Markov, 1995; Moore
and Bierbaum, 1979)

〈T∗−,T∗+〉=
(4.5)

=
〈
exp
[
−
(x− ξ)2
4a+1 (t

f − t)
]
,exp
[
−
(x− ξ)2
4a−1 (t

f − t)
]〉
/
D

[
2
√
π〈a−1 ,a

+
1 〉(tf − t)

]



Interval BEM for transient diffusion problem... 271

Thedirected interval heatfluxresulting fromthe interval fundamental solution
should be found in an analytical way, and then

q̃∗1(ξ,x,t
f , t)=−〈λ−1 ,λ

+
1 〉
∂T̃∗1 (ξ,x,t

f, t)

∂n
= 〈λ−1 ,λ

+
1 〉(x−ξ) ·

(4.6)

·
〈
exp
[
−
(x−ξ)2
4a+1 (t

f − t)
]
,exp
[
−
(x− ξ)2
4a−1 (t

f − t)
]〉
/
D
{4
√
π[〈a−1 ,a

+
1 〉(t

f − t)]3/2}

Theboundary integral equation corresponding to the transition tf−1 → tf
for the other layer can be expressed as follows

T̃2(ξ,t
f)+

[
1

c̃2
q̃2(x,t

f)

tf∫

tf−1

T̃∗2(ξ,x,t
f, t) dt

]x=L2

x=L1

=

(4.7)

=

[
1

c̃2
T̃2(x,t

f)

tf∫

tf−1

q̃∗2(ξ,x,t
f , t) dt

]x=L2

x=L1

+

L2∫

L1

T̃∗2(ξ,x,t
f , tf−1)T̃2(x,t

f−1) dx

where T̃∗2 (ξ,x,t
f, t) is thedirected interval fundamental solution for the second

sub-domain.
The numerical approximation of interval equations (4.2) and (4.7) leads to

the system of interval equations




−H̃111 −H̃112 G̃112 0
−H̃121 −H̃122 G̃122 0
0 −H̃211 G̃211 −H̃212
0 −H̃221 G̃221 −H̃222







T̃1(0, t
f)

T̃1(L1, t
f)

q̃(L1, t
f)

T̃2(L2, t
f)



=




−G̃111 · q̃L
−G̃121 · q̃L
−G̃212 · q̃R
−G̃222 · q̃R



+




P̃1(0, t
f−1)

P̃1(L1, t
f−1)

P̃2(L1, t
f−1)

P̃2(L2, t
f−1)




(4.8)
where

G̃e11 =−G̃e22 =−
√
∆t

√
λ̃ec̃eπ

L0 =0

(4.9)

G̃e12 =−G̃e21 =
√
∆t

√
λ̃ec̃eπ

exp
[
−
(Le−Le−1)2
4ãe∆t

]
−
Le−Le−1
2λ̃e

erfc
(Le−Le−1
2
√
ãe∆t

)

and

H̃e11 = H̃
e
22 =−

1

2
H̃e12 = H̃

e
21 =
1

2
erfc
(Le−Le−1
2
√
ãe∆t

)
(4.10)



272 A. Piasecka-Belkhayat

while

P̃1(0, t
f−1)=

1

2
√
πã1∆t

L1∫

0

exp
[
− x

2

4ã1∆t

]
T̃1(x,t

f−1) dx

P̃1(L1, t
f−1)=

1

2
√
πã1∆t

L1∫

0

exp
[
−(x−L1)

2

4ã1∆t

]
T̃1(x,t

f−1) dx

(4.11)

P̃2(L1, t
f−1)=

1

2
√
πã2∆t

L2∫

L1

exp
[
−(x−L1)

2

4ã2∆t

]
T̃2(x,t

f−1) dx

P̃2(L2, t
f−1)=

1

2
√
πã2∆t

L2∫

L1

exp
[
−
(x−L2)2
4ã2∆t

]
T̃2(x,t

f−1) dx

where e=1 denotes the first layer, e=2 denotes the other one.

For example, for the interval coefficients λ̃1 = 〈34,36〉 and c̃1 =
= 〈4900000,5400000〉, the product of λ̃1 and c̃1 is calculated according to
the formula

λ̃1 · c̃1 = 〈λ−+1 · c
−+
1 ,λ

+
1 · c

+
1 〉= 〈λ

−

1 · c
−

1 ,λ
+
1 · c

+
1 〉=

(4.12)
= 〈34 ·4900000,36 ·5400000〉= 〈1.666 ·108,1.944 ·108〉

The interval Gauss elimination method (Neumayer, 1990; Piasecka-
Belkhayat, 2008) has been used to solve the interval system of equations (4.8).
After determining the ’missing’ boundary values for both layers, the interval
temperatures T̃fe at the internal points ξ

i can be calculated using the formu-
las:

— for the first layer (ξ∈ (0,L1))

T̃1(ξ,t
f)=

1

2
exp
(
−
L1− ξ√
ã1∆t

)
T̃1(L1, t

f)+
1

2
exp
(
−
ξ√
ã1∆t

)
T̃1(0, t

f)+

+

√
∆t

2
√
λ̃1c̃1

exp
(
−
L1−ξ√
ã1∆t

)
q̃(L1, t

f)+ (4.13)

+

√
∆t

2
√
λ̃1c̃1

exp
(
− ξ√
ã1∆t

)
qL(0, t

f)+ P̃1(ξ,t
f−1)



Interval BEM for transient diffusion problem... 273

— for the other layer (ξ∈ (L1,L2))

T̃2(ξ,t
f)=

1

2
exp
(
−
L2− ξ√
ã2∆t

)
T̃2(L2, t

f)+
1

2
exp
(
−
ξ−L1√
ã2∆t

)
T̃2(L1, t

f)+

+

√
∆t

2
√
λ̃2c̃2

exp
(
−L2− ξ√
ã2∆t

)
qR(L2, t

f)+ (4.14)

+

√
∆t

2
√
λ̃2c̃2

exp
(
−
ξ−L1√
ã2∆t

)
q̃(L1, t

f)+ P̃2(ξ,t
f−1)

5. Numerical example

In the paper, an example of one-dimensional heat transient transfer in a
two-layered domain of dimensions L1 = 0.02m and L2 = 0.02m is pre-
sented. On both sides the boundary condition of the 2nd type of the form
q̃L = 〈9800,10200〉W/m2 and q̃R =0W/m2 has been assumed. The first and
other layer have been divided into 20 constant internal cells, respectively.

The following input data have been introduced: λ1 = 90W/(mK),
c1 = 3.916MJ/(m

3K), λ2 = 35W/(mK), c2 = 5.175MJ/(m
3K). The ini-

tial temperature of the first and other sub-domain is T01 =T02 =20
◦C, time

step ∆t=1s.

All thermophysical parameters of the two-layered domain are assumed to
be directed interval values:

λ̃1 = 〈λ1−0.05λ1,λ1+0.05λ1〉
c̃1 = 〈c1−0.05c1,c1+0.05c1〉
λ̃2 = 〈λ2−0.05λ2,λ2+0.05λ2〉
c̃2 = 〈c2−0.05c2,c2+0.05 · c2〉

Figure 1 illustrates the temperature distribution in the domain analysed
obtained for chosen times. The dashed and solid lines denote the lower and
the upper bounds of the temperature intervals, respectively.

Figure 2 shows a comparison between the temperature distribution obta-
ined for the time 80s and the results obtainedwith classical BEM for thermo-
physical parameters defined without intervals (dotted line).



274 A. Piasecka-Belkhayat

Fig. 1. Temperature distribution for chosen times

Fig. 2. Comparison of temperature distribution for the time 80s

6. Conclusions

In this paper, the description of an unsteady heat transfer for 1D problem for
the two-layered domain has been presented. All the thermophysical parame-
ters appearing in the mathematical model of the domain analysed have been
considered as directed interval values. The problem discussed has been solved
using the 1st scheme of the interval boundary element method according to
the rules of the directed interval arithmetic.

Themain advantage of the directed interval arithmetic upon the classical
interval arithmetic is that the obtained temperature intervals are much nar-
rower (Piasecka-Belkhayat, 2007). The problem analysed can be extended to
multi-layered domains.



Interval BEM for transient diffusion problem... 275

References

1. Brebbia C.A., Telles J.C.F., Wrobel L.C., 1984, Boundary Element
Techniques, Springer-Verlag

2. BurczyńskiT., SkrzypczykJ., 1997,Fuzzyaspects of theboundaryelement
method,Engineering Analysis with Boundary Elements, 19, 209-216

3. KużelewskiA., 2008, Implementationof amodeling and simulationalgorithm
for inaccuratelydefinedboundaryproblems,Ph.D.Thesis,TechnicalUniversity,
Białystok [in Polish]

4. Majchrzak E., 1998, Numerical modelling of bio-heat transfer using the bo-
undary element method, Journal of Theoretical and Applied Mechanics, 36, 2,
437-455

5. Majchrzak E., 2001, Boundary Element Method in Heat Transfer, Publ. of
the Technological University of Czestochowa, Częstochowa [in Polish]

6. Majchrzak E., Mochnacki B., 1996, The BEM application for numerical
solution of non-steady and non-linear thermal diffusion problems, Computer
Assisted Mechanics and Engineering Sciences, 3, 4, 327-346

7. Markov S.M., 1995, On directed interval arithmetic and its applications, Jo-
urnal of Universal Computer Science, 1, 514-526

8. Mochnacki B., Suchy J.S., 1995, Numerical Methods in Computations of
Foundry Processes, PFTA, Cracow

9. MooreR.E., BierbaumF., 1979,Methods andApplications of Interval Ana-
lysis, SIAM

10. Muhann R.L. Mullen R.L., Zhang H., 2005, Penalty-based solution for
the interval finite elementmethods, Journal of Engineering Mechanics, ASCE,
131, 1102-1111

11. Neumaier A., 1990, Interval Methods for System of Equations, Cambridge
University Press, Cambridge, NewYork, Port Chester, Melbourne, Sydney

12. Piasecka-BelkhayatA., 2007, Solution of energy equationusing the interval
boundary element method, Scientific Research of the Institute of Mathematics
and Computer Science, 6, 1, Częstochowa, 219-226

13. Piasecka-Belkhayat A., 2008, Interval boundary element method for 2D
transient diffusion problem,EngineeringAnalysis with Boundary Elements,32,
5, Elsevier, 424-430

14. Popova E.D., 2001, Multiplication distributivity of proper and improper in-
tervals,Reliable Computing, 7, 129-140

15. Zadeh L.A., 1965, Fuzzy sets, Information and Control, 8, 338-353



276 A. Piasecka-Belkhayat

Przedziałowa metoda elementów brzegowych dla zadań dyfuzji

w obszarach dwuwarstwowych

Streszczenie

W pracy przedstawiono opis nieustalonego przepływu ciepła dla zadań jednowy-
miarowychwobszarachdwuwarstwowych.Założono, że wszystkie parametry termofi-
zyczne pojawiające się w opisiematematycznymanalizowanego zadania sąwartościa-
mi przedziałowymi skierowanymi. Omawiane zagadnienie zostało rozwiązane za po-
mocą pierwszego schematu przedziałowejmetody elementów brzegowych.Do rozwią-
zania otrzymanego interwałowego układu równań zastosowano przedziałową metodę
eliminacji Gaussa.Wkońcowej części pracy pokazanowyniki obliczeń numerycznych.

Manuscript received April 19, 2010; accepted for print September 27, 2010