- Email: [email protected]
Finite element approaches to the PTC thermistor problem
Applied Mathematics and Computation 163 (2005) 147–162 www.elsevier.com/locate/amc
Finite element approaches to the PTC thermistor problem S. Kutluay, A. Esen Department of Mathematics, Faculty of Arts and Science, In€on€u University, 44069 Malatya, Turkey
Abstract In this paper, subdomain collocation (SC) and Petrov–Galerkin (PG) methods based on spline finite elements have been applied to the one-dimensional positive temperature coefficient (PTC) thermistor problem with a temperature dependent ramp function electrical conductivity to obtain the predicted temperature distributions and the locations of the free boundaries in the steady-state case. The numerical results obtained by the present methods have been compared with the exact solution and also those obtained by earlier authors. A Fourier stability analysis of each method is investigated. 2004 Elsevier Inc. All rights reserved. Keywords: PTC thermistor; Ramp electrical conductivity; Subdomain collocation; Petrov– Galerkin; B-spline functions
1. Introduction In this paper we consider an initial boundary value problem of coupled nonlinear partial differential equations which consists of a heat-flow equation with the joule heating as a source and a current flow equation with a temperature dependent electrical conductivity. This problem is a thermo-electric problem called as a PTC thermistor problem which is mathematically modeled as follows: Heat-flow problem: 2
Ut ¼ Uxx þ arð/x Þ ;
0 < x < 1; t > 0;
E-mail addresses: [email protected] (S. Kutluay), [email protected] (A. Esen). 0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.01.024
ð1Þ
148
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
Ux ð0; tÞ ¼ 0;
t > 0;
Ux ð1; tÞ þ bU ð1; tÞ ¼ 0; U ðx; 0Þ ¼ 0;
t > 0;
0 6 x 6 1:
ð2Þ ð3Þ
Current-flow problem: ðr/x Þx ¼ 0;
0 < x < 1; t > 0;
ð4Þ
/ð0; tÞ ¼ 0; /ð1; tÞ ¼ 1;
t > 0; t > 0;
ð5Þ
/ðx; 0Þ ¼ x;
0 6 x 6 1;
ð6Þ
where U ðx; tÞ is the temperature, /ðx; tÞ is the electrical potential, b is a surface heat transfer coefficient, a is the ratio of electrical conductivity and r is a temperature dependent electrical conductivity which decreases as the temperature increases. The existence, uniqueness and behaviour of solutions to this problem were given in [1–6]. Numerical solutions of the PTC thermistor problem have been undertaken by employing various forms of the finite element [7–9], finite difference [10–13] and heat balance integral methods [14,15]. In this paper, SC and PG methods have been applied to the PTC thermistor problem with an electrical conductivity proposed in the next section to obtain its approximate steady-state solutions. The numerical results obtained by both methods are compared with exact ones and also earlier work.
2. Electrical conductivities for the problem Many authors [7,11,15] considered the above given thermistor problem with a step function electrical conductivity of the form 1; 0 6 U 6 1; rðU Þ ¼ d 105 ; U > 1: In a postscript, Howison [16] makes the point that a step function is a crude representation for rðU Þ. Recently, Kutluay and Wood [12] introduced a slightly more realistic model for rðU Þ of the form 8 0 6 U 6 1; < 1; rðU Þ ¼ 1 þ ð1 U Þð1 dÞ; 1 < U < 2; ð7Þ : d 105 ; U P 2; which is mathematically equivalent to a ramp function. They have applied a standard explicit finite difference (EFD) method to obtain approximate steadystate solutions of the thermistor problem with the electrical conductivity (7).
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
149
They also showed that under the ramp electrical conductivity (7) the evolving thermistor problem defined by Eqs. (1)–(6) moves sequentially through cold, warm/cold, warm, hot/warm and hot phases in turn. When SC and PG methods are applied to the model problem, some complications or difficulties may arise in setting up weak formulations of warm/ cold and hot/warm phases since each of both phases consists of two regions with different conductivities separated from a free boundary as seen in Fig. 1. To over come difficulties in applications of SC and PG methods to the model problem, in the present paper we will use the following conductivities rðU Þ ¼ 1 þ ð1 U Þð1 dÞs1
ð8Þ
rðU Þ ¼ 1 þ ð1 U Þð1 dÞ þ ½d ð1 þ ð1 U Þð1 dÞÞs2
ð9Þ
and
for the cold/warm and hot/warm phases, respectively. Each of the proposed conductivities can be easily obtained using a linear interpolation by virtue of the conductivity (7). In the above equations, s1 is a free boundary separating the cold region from the warm region of the warm/cold phase and s2 is also a free boundary but it separates the hot region from the warm region of the hot/ warm phase. The locations of free boundaries s1 and s2 can be estimated by using a linear interpolation as sn1 ¼ xm þ
hð1 U ðxm ÞÞ U ðxmþ1 Þ U ðxm Þ
ð10Þ
sn2 ¼ xm þ
hð2 U ðxm ÞÞ ; U ðxmþ1 Þ U ðxm Þ
ð11Þ
and
respectively, where U ðxm Þ is an approximation to U ðx; tÞ at a node xm and a time level n. With the choice of the conductivities given by Eqs. (8) and (9), the U
(a)
U
2
2
1
1
Hot Warm
0
(b)
Warm
Cold s1
1
x
0
s2
1 x
Fig. 1. Steady-state configurations of (a) warm/cold and (b) hot/warm phases.
150
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
cold (s1 ¼ 0), warm (s1 ¼ 1 or s2 ¼ 0) and hot (s2 ¼ 1) phases are exactly the same with the conductivity (7) unlike the warm/cold and hot/warm phases known as transition phases should be reorganized. It means that each of nontransition phases (cold, warm and hot phases) satisfies exactly the same physical requirements obtained in [12]. We assume that /ðx; tÞ ¼ x in 0 6 x 6 1 for t P 0 and then the thermistor problem is reduced to the heat flow problem Ut ¼ Uxx þ ar;
0 < x < 1; t > 0;
ð12Þ
subject to the conditions (2) and (3). In this paper, we will only investigate transition phases (warm/cold and hot/ warm phases) in detail.
3. Exact solutions of the transition phases 3.1. Warm/cold solution e ðxÞ to this phase by using a We obtain the exact steady-state solution U standard analytical method as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 e ðxÞ ¼ 2c1 cosh U að1 dÞs1 x þ þ 1; ð13Þ ð1 dÞs1 where i pffiffiffiffiffiffiffiffiffiffiffiffiffi h 1 be að1dÞs1 ðd1Þs 1 1 ffi: c1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pffiffiffiffiffiffiffiffiffiffiffiffi b að1 dÞs1 þ b þ að1 dÞs1 e að1dÞs1 e ðs1 Þ ¼ 1 into Eq. (13) leads to a non-linear equation for the free Inserting U interface s1 : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c1 ð1 dÞs1 cosh s1 að1 dÞs1 þ 1 ¼ 0; ð14Þ which can be solved numerically. If s1 becomes equal to 1, Eq. (14) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ b¼ ð1 dÞ and so the problem enters a warm phase. Since 0 < s1 6 1 in this phase, we obtain an inequality, in terms of a; b and d, of the form pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ 6 b; ð15Þ ð1 dÞ
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
151
which is used to get warm/cold steady-state solutions. Inequality (15) giving a sub-bound for b is exactly the same with the one in [12] which is evident that the conductivity (8) can be used as a good representation for r to model the warm/cold phase. 3.2. Hot/warm solution e ðxÞ to this phase can be obtained by using The exact steady-state solution U a standard analytical method as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2s2 2d e ðxÞ ¼ 2c2 cosh U ; ð16Þ að1 dÞð1 s2 Þx þ 1 s2 ð1 dÞð1 s2 Þ where i pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h 2s2 2d be að1dÞð1s2 Þ 1s ð1dÞð1s2 Þ 2 c2 ¼ : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b að1 dÞð1 s2 Þ þ b þ að1 dÞð1 s2 Þ e að1dÞð1s2 Þ e ðs2 Þ ¼ 2, from Eq. (16) we obtain an equation for the interface s2 : Since U pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c2 ð1 dÞð1 s2 Þ cosh s2 að1 dÞð1 s2 Þ þ d ¼ 0; ð17Þ which can be solved numerically. For s2 ¼ 0, from Eq. (17) we obtain pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ b¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
: 2 d 1 þ cosh að1 dÞ Since 0 < s2 6 1 for all t, we obtain an inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ b< pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
; 2 d 1 þ cosh að1 dÞ
ð18Þ
which is used to find hot/warm steady-state solutions. This inequality giving an upper-bound for b is exactly the same with the one in [12] which is evident that the conductivity (9) can be used as a good representation for r to model the hot/warm phase.
4. The finite element solutions The region ½0; 1 is partitioned into N finite elements of equal length h by the knots xi such that 0 ¼ x0 < x1 < < xN 1 < xN ¼ 1. Over the element ½xm ; xmþ1 , an approximation to the exact solution U ðx; tÞ in terms of quadratic B-splines Qm with the required properties can be given as [17]
152
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
UN ¼ Qm1 dm1 þ Qm dm þ Qmþ1 dmþ1 ; where the dm are time dependent element parameters to be determined from the boundary and weighted residual conditions. In terms of a local coordinate system n given by n ¼ x xm ð0 6 n 6 hÞ the base functions Qm are 2 n 2n 2n2 n2 Qm1 ¼ 1 ; Qm ¼ 1 þ 2 ; Qmþ1 ¼ 2 : h h h h The nodal values Um and Um0 at the node in terms of element parameters dm can be obtained as Um ¼ U ðxm Þ ¼ dm1 þ dm
ð19Þ
hUm0 ¼ hU 0 ðxm Þ ¼ 2ðdm dm1 Þ;
ð20Þ
and
respectively. In the above equation and throughout this paper, the prime denotes derivative with respect to x. Applying a weighted residual technique we obtain a weighted residual form of Eq. (12) as Z 1 Wm ðUt Uxx arÞ dx ¼ 0; ð21Þ 0
where the Wm ðm ¼ 0ð1ÞN 1Þ are the weight functions. 4.1. Subdomain collocation (SC) method: In the SC method [18], the weight functions Wm are 1; xm 6 x 6 xmþ1 ; Wm ¼ 0; otherwise and so, the weighted residual form (21) over each finite element ½xm ; xmþ1 leads to the integral equation Z xmþ1 ðUt Uxx arÞ dx ¼ 0: ð22Þ xm
Substituting the electrical conductivities (8) and (9) into Eq. (22), and utilizing Eqs. (19) and (20), after a simple arrangement, we respectively obtain the following integral equations: 2 h _ dm1 þ 4d_ m þ d_ mþ1 ðdm1 2dm þ dmþ1 Þ 3 h ahð1 dÞs1 ðdm1 þ 4dm þ dmþ1 Þ ¼ hað1 þ ð1 dÞs1 Þ þ ð23Þ 3
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
and
153
2 h _ dm1 þ 4d_ m þ d_ mþ1 ðdm1 2dm þ dmþ1 Þ 3 h ahð1 dÞð1 s2 Þ þ ðdm1 þ 4dm þ dmþ1 Þ ¼ hað2 dÞ 2hað1 dÞs2 : 3 ð24Þ
In the above equation and throughout this paper, the dot denotes differentiation with respect to t. Inserting Eqs. (19) and (20) into Eqs. (23) and (24), and using the following forward difference for d_ and the Crank–Nicolson approach for d [19] d_ ¼ ðdnþ1 dn Þ=Dt;
d ¼ ðdn þ dnþ1 Þ=2;
ð25Þ
we obtain the following systems of N linear algebraic equations with N þ 2 unknown as nþ1 nþ1 ð1 3r þ c1 Þdnþ1 m1 þ ð4 þ 6r þ 4c1 Þdm þ ð1 3r þ c1 Þdmþ1
¼ ð1 þ 3r c1 Þdnm1 þ ð4 6r 4c1 Þdnm þ ð1 þ 3r c1 Þdnmþ1 þ 3kað1 þ ð1 dÞs1 Þ
ð26Þ
and nþ1 nþ1 ð1 3r þ c2 Þdnþ1 m1 þ ð4 þ 6r þ 4c2 Þdm þ ð1 3r þ c2 Þdmþ1
¼ ð1 þ 3r c2 Þdnm1 þ ð4 6r 4c2 Þdnm þ ð1 þ 3r c2 Þdnmþ1 þ 3kað2 dÞ 6kað1 dÞs2 ;
ð27Þ
respectively, where Dtð kÞ is the time step, r ¼ k=h2 , c1 ¼ kað1 dÞs1 =2, and c2 ¼ kað1 dÞð1 s2 Þ=2. The element parameters d1 and dN appearing in Eqs. (26) and (27) are eliminated by using the boundary conditions (2). Thus, each one of the recurrence relationships (26) and (27) can be written as nþ1 ð5 þ 3r þ 5c1 Þdmnþ1 þ ð1 3r þ c1 Þdmþ1 ¼ ð5 3r 5c1 Þdnm
þ ð1 þ 3r c1 Þdnmþ1 þ 3kað1 þ ð1 dÞs1 Þ;
m ¼ 0;
nþ1 nþ1 ð1 3r þ c1 Þdm1 þ ð4 þ 6r þ 4c1 Þdmnþ1 þ ð1 3r þ c1 Þdmþ1
¼ ð1 þ 3r c1 Þdnm1 þ ð4 6r 4c1 Þdnm þ ð1 þ 3r c1 Þdnmþ1 þ 3kað1 þ ð1 dÞs1 Þ;
m ¼ 1ð1ÞN 2;
nþ1 ð1 3r þ c1 Þdm1 þ ðð4 þ kÞð1 þ c1 Þ þ 3rð2 kÞÞdmnþ1
¼ ð1 þ 3r c1 Þdnm1 þ ðð4 þ kÞð1 c1 Þ 3rð2 kÞÞdnm þ 3kað1 þ ð1 dÞs1 Þ;
m¼N 1
ð28Þ
154
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
and nþ1 ð5 þ 3r þ 5c2 Þdmnþ1 þ ð1 3r þ c2 Þdmþ1 ¼ ð5 3r 5c2 Þdnm
þ ð1 þ 3r c2 Þdnmþ1 þ 3kað2 dÞ 6kað1 dÞs2 ; ð1 3r þ
nþ1 c2 Þdm1
þ ð4 þ 6r þ
4c2 Þdmnþ1
m ¼ 0;
nþ1 þ ð1 3r þ c2 Þdmþ1
¼ ð1 þ 3r c2 Þdnm1 þ ð4 6r 4c2 Þdnm þ ð1 þ 3r c2 Þdnmþ1 þ 3kað2 dÞ 6kað1 dÞs2 ; ð1 3r þ
nþ1 c2 Þdm1
þ ðð4 þ kÞð1 þ c2 Þ þ 3rð2
m ¼ 1ð1ÞN 2;
kÞÞdmnþ1
¼ ð1 þ 3r c2 Þdnm1 þ ðð4 þ kÞð1 c2 Þ 3rð2 kÞÞdnm þ 3kað2 dÞ 6kað1 dÞs2 ;
m ¼ N 1; ð29Þ
respectively, where k ¼ ð2 hbÞ=ð2 þ hbÞ. Each of the above systems can be easily and efficiently solved by using the Thomas algorithm. A Fourier stability analysis shows that the schemes (28) and (29) are unconditionally stable since the Crank–Nicolson approach is used for the element parameters d. 4.2. Petrov–Galerkin (PG) method The weighted residual form (21) over each finite element ½xm ; xmþ1 , with integration by parts, leads to the weak form Z xmþ1 ðWUt þ Wx Ux W arÞ dx ¼ bW ðxN ÞU ðxN ; tÞ: ð30Þ xm
The weight functions Wm are taken as linear B-spline basis functions for the finite element ½xm ; xmþ1 of the form [17] Wm ¼ 1 n=h;
Wmþ1 ¼ n=h:
Substituting the electrical conductivities (8) and (9) into Eq. (30), after a simple arrangement we respectively obtain the following equations: Z xmþ1 ðWUt þ Wx Ux þ að1 dÞs1 WU að1 þ ð1 dÞs1 ÞW Þ dx xm
¼ bW ðxN ÞU ðxN ; tÞ
ð31Þ
and Z
xmþ1
ðWUt þ Wx Ux þ að1 dÞð1 s2 ÞWU aðð2 dÞ 2ð1 dÞs2 ÞW Þ dx
xm
¼ bW ðxN ÞU ðxN ; tÞ:
ð32Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
155
Inserting Eqs. (19) and (20) into Eqs. (31) and (32) lead to mþ1 Z xmþ1 mþ1 Z xmþ1 X X e 0 0 _ Wi Qj dx dj þ Wi Qj dx dej xm
j¼m1
xm
j¼m1
þ að1 dÞs1
mþ1 Z X j¼m1
¼ að1 þ ð1 dÞs1 Þ
xmþ1
Wi Qj dx dej þ b
xm
Z
Wi ðxN ÞQj ðxN Þdej
j¼m1
xmþ1
mþ1 X
Wi dx xm
and mþ1 Z X
mþ1 Z X e _ Wi Qj dx dj þ
xm
j¼m1
xmþ1
mþ1 X j¼m1
xmþ1
0
Wi Q0j dx
dej þ að1 dÞð1 s2 Þ
xm
j¼m1
Z
xmþ1
mþ1 X Wi Qj dx dej þ b Wi ðxN ÞQj ðxN Þdej
xm
¼ aðð2 dÞ 2ð1 dÞs2 Þ
Z
j¼m1 xmþ1
Wi dx ; xm
which can also be written in matrix forms as Ae d_ e þ ðc3 Ae þ Be þ C e Þde ¼ ða þ c3 Þd e
ð33Þ
Ae d_ e þ ðc4 Ae þ Be þ C e Þde ¼ ðað1 ð1 dÞs2 Þ þ c4 Þd e ;
ð34Þ
and
respectively, where c3 ¼ að1 dÞs1 , c4 ¼ að1 dÞð1 s2 Þ, de ¼ ðdm1 ; dm ; dmþ1 ÞT are the element parameters, and Ae ; Be ; C e and d e are the element matrices given by the following integrals: Z xmþ1 Aeij ¼ Wi Qj dx; xm Z xmþ1 Wi 0 Q0j dx; Beij ¼ ð35Þ xm e Cij ¼ bWi ðxN ÞQj ðxN Þ; Z xmþ1 die ¼ Wi dx; xm
where i ¼ m; m þ 1 and j ¼ m 1; m; m þ 1. The element matrices (35) are obtained as h 3 8 1 1 1 0 1 h 1 Ae ¼ ; Be ¼ ; de ¼ 12 1 8 3 h 1 0 1 2 1
156
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
and C e is a zero matrix but 0 0 0 : C N 1 ¼ 0 b b Assembling contributions from all elements, Eqs. (33) and (34) respectively yield the following matrix equations: Ad_ þ ðc3 A þ B þ CÞd ¼ ða þ c3 Þd
ð36Þ
Ad_ þ ðc4 A þ B þ CÞd ¼ ðað1 ð1 dÞs2 Þ þ c4 Þd;
ð37Þ
and
T
where d ¼ ðd1 ; d0 ; . . . ; dN Þ and A; B, C are ðN þ 1Þ ðN þ 2Þ matrices, and d is an ðN þ 1Þ 1 matrix which are derived from the corresponding element matrices Ae ; Be ; C e and d e . Substituting the approximations (25) into Eqs. (36) and (37) lead to equations ½ð2 þ c3 kÞA þ k ð B þ CÞdnþ1 ¼ ½ð2 c3 kÞA k ð B þ CÞdn þ 2k ða þ c3 Þd ð38Þ and ½ð2 þ c4 kÞA þ kðB þ CÞdnþ1 ¼ ½ð2 c4 kÞA kðB þ CÞdn þ 2kðað1 ð1 dÞs2 Þ þ c4 Þd;
ð39Þ
respectively. Utilizing one of the boundary conditions (2) in Eqs. (38) and (39) yield an ðN þ 1Þ ðN þ 1Þ pentadiagonal matrix system and so is easily and efficiently solved with a variant of the Thomas algorithm. Again a Fourier stability analysis shows that schemes (38) and (39) are unconditionally stable since the Crank–Nicolson approach is used for the element parameters de .
5. Numerical results and conclusion In order to obtain the steady-state temperature distribution U ðx; tÞ and the locations of the free boundaries s1 and s2 , respectively, in the warm/cold and hot/warm phases constructed in the previous sections we used criteria 1 max U0new U0old ; UNnew UNold 6 0:5 108 k with k ¼ 0:01 for both the SC and PG solutions. In order to show how good the numerical solutions of the problem in comparison with the exact ones we used the error norm L2 defined by
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
L2 ¼
N 2 1 X e ðxm Þ U ðxm Þ U N þ 1 m¼0
157
!1=2 :
To compare our results obtained by using the SC and PG methods with the previous numerical results of Kutluay and Wood [12] using the EFD method we used values for a ¼ 0:5; 102 with b ¼ 0:5 satisfying inequalities (15) and (18) to get numerical solutions in the warm/cold and hot/warm phases, respectively. All the numerical computations were executed on a Pentium 4 PC in the Fortran code using double-precision arithmetic to reduce the round-off error to a minimum. The SC solutions to the warm/cold phase (for a ¼ 0:5 and b ¼ 0:5) are obtained by using Eq. (28) and are given in Table 1 with the norm L2 for a variety values of the mesh size h. The norm L2 decreases from 0.319602 · 104 (for h ¼ 0:1) to 0.020591 · 104 (for h ¼ 0:025). It is clearly seen that the predicted temperature values converge to the exact ones as the mesh size h is refined. Table 2 displays the PG solutions of the warm/cold phase obtained by using Eq. (38) with the norm L2 . It is seen from the table that the agreement of the numerical solutions with the exact one is remarkably good. The results obtained by SC and PG methods are presented in Table 3 with the corresponding phase solutions given by Kutluay and Wood [12]. The solutions obtained by the present two methods (for h ¼ 0:05) coincide to five significant digits in the most cases and agree with those of Kutluay and Wood [12]. Table 4 shows a comparison of the predicted motion (location and speed) of the free boundary s1 obtained by virtue of Eq. (10) for k ¼ 0:01 and various mesh sizes. It is seen that the estimated values of s1 and s_ 1 obtained by the present methods Table 1 Values of the temperature distribution as predicted by the numerical (SC method) and exact solutions in the warm/cold phase x
Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.156290 1.154116 1.147582 1.136662 1.121311 1.101463 1.077037 1.047929 1.014018 0.975162 0.931199
1.156191 1.154017 1.147486 1.136570 1.121224 1.101383 1.076965 1.047867 1.013967 0.975123 0.931173
1.156170 1.153996 1.147466 1.136551 1.121206 1.101366 1.076950 1.047854 1.013956 0.975116 0.931169
L2 104
0.319602
0.063870
0.020591
1.156154 1.153980 1.147450 1.136535 1.121191 1.101353 1.076937 1.047843 1.013946 0.975107 0.931162
158
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
Table 2 Values of the temperature distribution as predicted by the numerical (PG method) and exact solutions in the warm/cold phase Numerical
x
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.156242 1.154067 1.147536 1.136619 1.121273 1.101431 1.077012 1.047913 1.014013 0.975170 0.931221
1.156179 1.154005 1.147474 1.136559 1.121214 1.101375 1.076958 1.047863 1.013965 0.975125 0.931178
1.156168 1.153994 1.147463 1.136549 1.121204 1.101365 1.076949 1.047854 1.013957 0.975117 0.931171
L2 104
0.231638
0.047248
0.018844
1.156154 1.153980 1.147450 1.136535 1.121191 1.101353 1.076937 1.047843 1.013946 0.975107 0.931162
Table 3 Comparison of numerical (SC and PG) solutions of the warm/cold phase with results from [12] for h ¼ 0:05 and k ¼ 0:01 x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Present
[12]
SC
PG
EFD
Exact
1.156191 1.154017 1.147486 1.136570 1.121224 1.101383 1.076965 1.047867 1.013967 0.975123 0.931173
1.156179 1.154005 1.147474 1.136559 1.121214 1.101375 1.076958 1.047863 1.013965 0.975125 0.931178
1.144337 1.142197 1.135767 1.125014 1.109885 1.090304 1.066174 1.037373 1.003757 0.965177 0.921597
1.145311 1.143173 1.136749 1.126008 1.110894 1.091334 1.067228 1.038456 1.004874 0.966255 0.922624
converge to the exact ones as the mesh size h is refined. For h ¼ 0:025, the norm L2 is 0.020591 · 104 in the SC method whereas it is 0.018844 · 104 in the PG method. Observe that the PG method needs less steady-state time ðt1 Þ and generates slightly more accurate solutions than the SC method. The SC solutions to the hot/warm phase (for a ¼ 102 and b ¼ 0:5) are obtained by using Eq. (29) and are displayed in Table 5 with the norm L2 for various values of h. The norm L2 decreases from 4.725202 ·104 (for h ¼ 0:1) to 0.125614 · 104 (for h ¼ 0:025). It is seen that the SC solutions are in good agreement with the exact solution and better accuracy can be achieved for smaller values of h. Table 6 presents the PG solutions of the hot/warm phase
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
159
Table 4 Comparison of the predicted location and speed of the free boundary s1 for k ¼ 0:01 and various mesh sizes at steady-state time ðt1 Þ Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
SC
s1 s_ 1 107 t1
0.836077 0.112404 22.17
0.837157 0.115755 22.13
0.837342 0.114231 22.11
0.837581 0.0 –
PG
s1 s_ 1 107 t1
0.836076 2.290229 19.29
0.837156 2.178131 19.35
0.837344 1.802939 19.63
0.837581 0.0 –
Table 5 Values of the temperature distribution as predicted by the numerical (SC method) and exact solutions in the hot/warm phase (e1 ¼ 4 107 , e2 ¼ 3 107 ) x
Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
2.000003 1.999998 1.999980 1.999926 1.999779 1.999372 1.998250 1.995160 1.986642 1.963170 1.898485
2.000001 1.999996 1.999975 1.999916 1.999757 1.999328 1.998168 1.995038 1.986590 1.963790 1.902255
2 þ e1 1.999995 1.999974 1.999914 1.999753 1.999319 1.998153 1.995021 1.986605 1.963993 1.903238
L2 104
4.725202
0.727451
0.125614
2 þ e2 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903601
obtained by using Eq. (39) with the norm L2 . We can see from the table that the agreement of the PG solutions with the exact solution is good enough. The results obtained by SC and PG methods are shown in Table 7 with the corresponding phase solutions given by Kutluay and Wood [12]. They are all in very good agreement with each other. Table 8 presents a comparison of the steady-state values of the numerical motion (location and speed) of the free boundary s2 obtained by virtue of Eq. (11) for k ¼ 0:01 and various mesh sizes. It is seen from Table 8 that the predicted values of s2 and s_ 2 obtained by SC and PG methods are reasonably in good agreement and also exhibit the expected convergence as h is refined. For h ¼ 0:025, the norm L2 is 0.125614 · 104 in the SC method whereas it is 0.000815 · 104 in the PG method. It is clearly seen that the PG method produces slightly more accurate solutions than the SC
160
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
Table 6 Values of the temperature distribution as predicted by the numerical (PG method) and exact solutions in the hot/warm phase (e1 ¼ 4 107 ; e2 ¼ 3 107 ) Numerical
x
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
2.000001 1.999996 1.999976 1.999918 1.999763 1.999342 1.998208 1.995143 1.986870 1.964535 1.904234
2.000001 1.999995 1.999974 1.999916 1.999756 1.999327 1.998173 1.995067 1.986715 1.964248 1.903815
2+e1 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964089 1.903602
L2 104
0.748935
0.178411
0.000815
2+e2 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903601
Table 7 Comparison of numerical (SC and PG) solutions of the hot/warm phase with results from [12] for h ¼ 0:05 and k ¼ 0:01 x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Present
[12]
SC
PG
EFD
Exact
2.000001 1.999996 1.999975 1.999916 1.999757 1.999328 1.998168 1.995038 1.986590 1.963790 1.902255
2.000001 1.999995 1.999974 1.999916 1.999756 1.999327 1.998173 1.995067 1.986715 1.964248 1.903815
2.000008 2.000003 1.999985 1.999929 1.999776 1.999364 1.998255 1.995271 1.987241 1.965634 1.907492
2.000002 1.999997 1.999978 1.999923 1.999774 1.999368 1.998266 1.995268 1.987121 1.964973 1.904771
method. However, the SC method needs less steady-state time ðt1 Þ than the PG method. In this paper, subdomain collocation and Petrov–Galerkin B-spline finite element methods were efficiently applied to the PTC thermistor problem with the electrical conductivities defined by expressions (8) and (9) in Section 1 in order to obtain steady-state temperature distributions and the motion (location and speed) of the free interfaces. It is shown that all the numerical results obtained by using SC and PG methods display the correct physical phenomena of the problem, and are reasonably in good agreement with each other and
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
161
Table 8 Comparison of the predicted location and speed of the free boundary s2 for k ¼ 0:01 and various mesh sizes at steady-state time ðt1 Þ Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
SC
s2 s_ 2 107 t1
0.057120 1.550098 0.96
0.034289 1.497480 1.40
0.028094 0.026169 4.63
0.025334 0.0 –
PG
s2 s_ 2 107 t1
0.017410 2.348426 0.68
0.021160 0.018321 2.01
0.025326 0.011836 6.91
0.025334 0.0 –
exact ones. In general, more accurate numerical solutions will result from smaller mesh size h. However, more computer storage and computations will be required for smaller mesh size. As a conclusion, the electrical conductivities proposed in this paper can be used as a good representation for r to model the PTC thermistor problem.
References [1] G. Cimatti, Remark on existence and uniqueness for the thermistor problem under mixed boundary conditions, Q. Appl. Math. 47 (1989) 117–121. [2] G. Cimatti, G. Prodi, Existence results for a non-linear elliptic system modelling a temperature dependent electrical resistor, Ann. Mat. Pura Appl. 162 (1992) 33–42. [3] A. Fowler, I. Frigaard, S.D. Howison, Temperature surges in thermistors, SIAM J. Appl. Math. 52 (1992) 998–1011. [4] S.D. Howison, J.F. Rodrigues, M. Shillor, Stationary solutions to the thermistor problem, J. Math. Anal. Appl. 174 (1993) 573–588. [5] J. Wiedmann, The thermistor problem, Nonlinear Differ. Equations Appl. 4 (1997) 133–148. [6] S. Zhou, D.R. Westbrook, Numerical solutions of the thermistor equations, J. Comput. Appl. Math. 79 (1997) 101–118. [7] D.R. Westbrook, The thermistor: a problem in heat and current flow, Numer. Meth. PDEs 5 (1989) 259–273. [8] K. Preis, O. Biro, R. Dyczij-Edlinger, K.R. Richter, Zs. Badics, H. Riedler, H. St€ onger, Application of FEM to coupled electric, thermal and mechanical problems, IEEE Trans. Magn. 30 (5) (1994) 3316–3319. [9] S. Kutluay, A. Esen, A B-spline finite element method for the thermistor with the modified electrical conductivity, Appl. Math. Comput., AMC 8384, in press. [10] S. Kutluay, A.R. Bahadir, A. Ozdes, A variety of finite difference methods to the thermistor with a new modified electrical conductivity, Appl. Math. Comput. 106 (1999) 205–213. [11] A.S. Wood, Modelling the thermistor, in: M. Heili€ o (Ed.), Proceedings of 5th European Conference on Maths. Industry, Kluwer Academic, Stuttgart, 1991, pp. 397–400. [12] S. Kutluay, A.S. Wood, Numerical solutions of the thermistor problem with a ramp electrical conductivity, Appl. Math. Comput. 148 (2004) 145–162. [13] S. Kutluay, A.R. Bahadir, A. Ozdes, Various methods to the thermistor problem with a bulk electrical conductivity, Int. J. Numer. Meth. Eng. 45 (1999) 1–12.
162
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
[14] S. Kutluay, A.S. Wood, A. Esen, A heat balance integral approach to the thermistor problem with a modified electrical conductivity, submitted for publication. [15] A.S. Wood, S. Kutluay, A heat balance integral model of the thermistor, Int. J. Heat Mass Transfer 38 (10) (1995) 1831–1840. [16] S. Howison, Complex variables in industrial mathematics, in: H. Neunzert (Ed.), Proceedings of 2nd European Symposium Maths. Industry, Oberwolfach, Kluwer Academic, Stuttgart, 1988, pp. 153–166. [17] P.M. Prenter, Splines and Variational Methods, John-Wiley, New York, 1975. [18] O.C. Zienkiewics, R.L. Taylor, The Finite Element Method, Butterworth-Heinemann, Oxford, 2002. [19] G.D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Clarendon, Oxford, 1987.
Finite element approaches to the PTC thermistor problem S. Kutluay, A. Esen Department of Mathematics, Faculty of Arts and Science, In€on€u University, 44069 Malatya, Turkey
Abstract In this paper, subdomain collocation (SC) and Petrov–Galerkin (PG) methods based on spline finite elements have been applied to the one-dimensional positive temperature coefficient (PTC) thermistor problem with a temperature dependent ramp function electrical conductivity to obtain the predicted temperature distributions and the locations of the free boundaries in the steady-state case. The numerical results obtained by the present methods have been compared with the exact solution and also those obtained by earlier authors. A Fourier stability analysis of each method is investigated. 2004 Elsevier Inc. All rights reserved. Keywords: PTC thermistor; Ramp electrical conductivity; Subdomain collocation; Petrov– Galerkin; B-spline functions
1. Introduction In this paper we consider an initial boundary value problem of coupled nonlinear partial differential equations which consists of a heat-flow equation with the joule heating as a source and a current flow equation with a temperature dependent electrical conductivity. This problem is a thermo-electric problem called as a PTC thermistor problem which is mathematically modeled as follows: Heat-flow problem: 2
Ut ¼ Uxx þ arð/x Þ ;
0 < x < 1; t > 0;
E-mail addresses: [email protected] (S. Kutluay), [email protected] (A. Esen). 0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2004.01.024
ð1Þ
148
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
Ux ð0; tÞ ¼ 0;
t > 0;
Ux ð1; tÞ þ bU ð1; tÞ ¼ 0; U ðx; 0Þ ¼ 0;
t > 0;
0 6 x 6 1:
ð2Þ ð3Þ
Current-flow problem: ðr/x Þx ¼ 0;
0 < x < 1; t > 0;
ð4Þ
/ð0; tÞ ¼ 0; /ð1; tÞ ¼ 1;
t > 0; t > 0;
ð5Þ
/ðx; 0Þ ¼ x;
0 6 x 6 1;
ð6Þ
where U ðx; tÞ is the temperature, /ðx; tÞ is the electrical potential, b is a surface heat transfer coefficient, a is the ratio of electrical conductivity and r is a temperature dependent electrical conductivity which decreases as the temperature increases. The existence, uniqueness and behaviour of solutions to this problem were given in [1–6]. Numerical solutions of the PTC thermistor problem have been undertaken by employing various forms of the finite element [7–9], finite difference [10–13] and heat balance integral methods [14,15]. In this paper, SC and PG methods have been applied to the PTC thermistor problem with an electrical conductivity proposed in the next section to obtain its approximate steady-state solutions. The numerical results obtained by both methods are compared with exact ones and also earlier work.
2. Electrical conductivities for the problem Many authors [7,11,15] considered the above given thermistor problem with a step function electrical conductivity of the form 1; 0 6 U 6 1; rðU Þ ¼ d 105 ; U > 1: In a postscript, Howison [16] makes the point that a step function is a crude representation for rðU Þ. Recently, Kutluay and Wood [12] introduced a slightly more realistic model for rðU Þ of the form 8 0 6 U 6 1; < 1; rðU Þ ¼ 1 þ ð1 U Þð1 dÞ; 1 < U < 2; ð7Þ : d 105 ; U P 2; which is mathematically equivalent to a ramp function. They have applied a standard explicit finite difference (EFD) method to obtain approximate steadystate solutions of the thermistor problem with the electrical conductivity (7).
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
149
They also showed that under the ramp electrical conductivity (7) the evolving thermistor problem defined by Eqs. (1)–(6) moves sequentially through cold, warm/cold, warm, hot/warm and hot phases in turn. When SC and PG methods are applied to the model problem, some complications or difficulties may arise in setting up weak formulations of warm/ cold and hot/warm phases since each of both phases consists of two regions with different conductivities separated from a free boundary as seen in Fig. 1. To over come difficulties in applications of SC and PG methods to the model problem, in the present paper we will use the following conductivities rðU Þ ¼ 1 þ ð1 U Þð1 dÞs1
ð8Þ
rðU Þ ¼ 1 þ ð1 U Þð1 dÞ þ ½d ð1 þ ð1 U Þð1 dÞÞs2
ð9Þ
and
for the cold/warm and hot/warm phases, respectively. Each of the proposed conductivities can be easily obtained using a linear interpolation by virtue of the conductivity (7). In the above equations, s1 is a free boundary separating the cold region from the warm region of the warm/cold phase and s2 is also a free boundary but it separates the hot region from the warm region of the hot/ warm phase. The locations of free boundaries s1 and s2 can be estimated by using a linear interpolation as sn1 ¼ xm þ
hð1 U ðxm ÞÞ U ðxmþ1 Þ U ðxm Þ
ð10Þ
sn2 ¼ xm þ
hð2 U ðxm ÞÞ ; U ðxmþ1 Þ U ðxm Þ
ð11Þ
and
respectively, where U ðxm Þ is an approximation to U ðx; tÞ at a node xm and a time level n. With the choice of the conductivities given by Eqs. (8) and (9), the U
(a)
U
2
2
1
1
Hot Warm
0
(b)
Warm
Cold s1
1
x
0
s2
1 x
Fig. 1. Steady-state configurations of (a) warm/cold and (b) hot/warm phases.
150
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
cold (s1 ¼ 0), warm (s1 ¼ 1 or s2 ¼ 0) and hot (s2 ¼ 1) phases are exactly the same with the conductivity (7) unlike the warm/cold and hot/warm phases known as transition phases should be reorganized. It means that each of nontransition phases (cold, warm and hot phases) satisfies exactly the same physical requirements obtained in [12]. We assume that /ðx; tÞ ¼ x in 0 6 x 6 1 for t P 0 and then the thermistor problem is reduced to the heat flow problem Ut ¼ Uxx þ ar;
0 < x < 1; t > 0;
ð12Þ
subject to the conditions (2) and (3). In this paper, we will only investigate transition phases (warm/cold and hot/ warm phases) in detail.
3. Exact solutions of the transition phases 3.1. Warm/cold solution e ðxÞ to this phase by using a We obtain the exact steady-state solution U standard analytical method as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 e ðxÞ ¼ 2c1 cosh U að1 dÞs1 x þ þ 1; ð13Þ ð1 dÞs1 where i pffiffiffiffiffiffiffiffiffiffiffiffiffi h 1 be að1dÞs1 ðd1Þs 1 1 ffi: c1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pffiffiffiffiffiffiffiffiffiffiffiffi b að1 dÞs1 þ b þ að1 dÞs1 e að1dÞs1 e ðs1 Þ ¼ 1 into Eq. (13) leads to a non-linear equation for the free Inserting U interface s1 : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c1 ð1 dÞs1 cosh s1 að1 dÞs1 þ 1 ¼ 0; ð14Þ which can be solved numerically. If s1 becomes equal to 1, Eq. (14) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ b¼ ð1 dÞ and so the problem enters a warm phase. Since 0 < s1 6 1 in this phase, we obtain an inequality, in terms of a; b and d, of the form pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ 6 b; ð15Þ ð1 dÞ
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
151
which is used to get warm/cold steady-state solutions. Inequality (15) giving a sub-bound for b is exactly the same with the one in [12] which is evident that the conductivity (8) can be used as a good representation for r to model the warm/cold phase. 3.2. Hot/warm solution e ðxÞ to this phase can be obtained by using The exact steady-state solution U a standard analytical method as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2s2 2d e ðxÞ ¼ 2c2 cosh U ; ð16Þ að1 dÞð1 s2 Þx þ 1 s2 ð1 dÞð1 s2 Þ where i pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h 2s2 2d be að1dÞð1s2 Þ 1s ð1dÞð1s2 Þ 2 c2 ¼ : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b að1 dÞð1 s2 Þ þ b þ að1 dÞð1 s2 Þ e að1dÞð1s2 Þ e ðs2 Þ ¼ 2, from Eq. (16) we obtain an equation for the interface s2 : Since U pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c2 ð1 dÞð1 s2 Þ cosh s2 að1 dÞð1 s2 Þ þ d ¼ 0; ð17Þ which can be solved numerically. For s2 ¼ 0, from Eq. (17) we obtain pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ b¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
: 2 d 1 þ cosh að1 dÞ Since 0 < s2 6 1 for all t, we obtain an inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ b< pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
; 2 d 1 þ cosh að1 dÞ
ð18Þ
which is used to find hot/warm steady-state solutions. This inequality giving an upper-bound for b is exactly the same with the one in [12] which is evident that the conductivity (9) can be used as a good representation for r to model the hot/warm phase.
4. The finite element solutions The region ½0; 1 is partitioned into N finite elements of equal length h by the knots xi such that 0 ¼ x0 < x1 < < xN 1 < xN ¼ 1. Over the element ½xm ; xmþ1 , an approximation to the exact solution U ðx; tÞ in terms of quadratic B-splines Qm with the required properties can be given as [17]
152
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
UN ¼ Qm1 dm1 þ Qm dm þ Qmþ1 dmþ1 ; where the dm are time dependent element parameters to be determined from the boundary and weighted residual conditions. In terms of a local coordinate system n given by n ¼ x xm ð0 6 n 6 hÞ the base functions Qm are 2 n 2n 2n2 n2 Qm1 ¼ 1 ; Qm ¼ 1 þ 2 ; Qmþ1 ¼ 2 : h h h h The nodal values Um and Um0 at the node in terms of element parameters dm can be obtained as Um ¼ U ðxm Þ ¼ dm1 þ dm
ð19Þ
hUm0 ¼ hU 0 ðxm Þ ¼ 2ðdm dm1 Þ;
ð20Þ
and
respectively. In the above equation and throughout this paper, the prime denotes derivative with respect to x. Applying a weighted residual technique we obtain a weighted residual form of Eq. (12) as Z 1 Wm ðUt Uxx arÞ dx ¼ 0; ð21Þ 0
where the Wm ðm ¼ 0ð1ÞN 1Þ are the weight functions. 4.1. Subdomain collocation (SC) method: In the SC method [18], the weight functions Wm are 1; xm 6 x 6 xmþ1 ; Wm ¼ 0; otherwise and so, the weighted residual form (21) over each finite element ½xm ; xmþ1 leads to the integral equation Z xmþ1 ðUt Uxx arÞ dx ¼ 0: ð22Þ xm
Substituting the electrical conductivities (8) and (9) into Eq. (22), and utilizing Eqs. (19) and (20), after a simple arrangement, we respectively obtain the following integral equations: 2 h _ dm1 þ 4d_ m þ d_ mþ1 ðdm1 2dm þ dmþ1 Þ 3 h ahð1 dÞs1 ðdm1 þ 4dm þ dmþ1 Þ ¼ hað1 þ ð1 dÞs1 Þ þ ð23Þ 3
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
and
153
2 h _ dm1 þ 4d_ m þ d_ mþ1 ðdm1 2dm þ dmþ1 Þ 3 h ahð1 dÞð1 s2 Þ þ ðdm1 þ 4dm þ dmþ1 Þ ¼ hað2 dÞ 2hað1 dÞs2 : 3 ð24Þ
In the above equation and throughout this paper, the dot denotes differentiation with respect to t. Inserting Eqs. (19) and (20) into Eqs. (23) and (24), and using the following forward difference for d_ and the Crank–Nicolson approach for d [19] d_ ¼ ðdnþ1 dn Þ=Dt;
d ¼ ðdn þ dnþ1 Þ=2;
ð25Þ
we obtain the following systems of N linear algebraic equations with N þ 2 unknown as nþ1 nþ1 ð1 3r þ c1 Þdnþ1 m1 þ ð4 þ 6r þ 4c1 Þdm þ ð1 3r þ c1 Þdmþ1
¼ ð1 þ 3r c1 Þdnm1 þ ð4 6r 4c1 Þdnm þ ð1 þ 3r c1 Þdnmþ1 þ 3kað1 þ ð1 dÞs1 Þ
ð26Þ
and nþ1 nþ1 ð1 3r þ c2 Þdnþ1 m1 þ ð4 þ 6r þ 4c2 Þdm þ ð1 3r þ c2 Þdmþ1
¼ ð1 þ 3r c2 Þdnm1 þ ð4 6r 4c2 Þdnm þ ð1 þ 3r c2 Þdnmþ1 þ 3kað2 dÞ 6kað1 dÞs2 ;
ð27Þ
respectively, where Dtð kÞ is the time step, r ¼ k=h2 , c1 ¼ kað1 dÞs1 =2, and c2 ¼ kað1 dÞð1 s2 Þ=2. The element parameters d1 and dN appearing in Eqs. (26) and (27) are eliminated by using the boundary conditions (2). Thus, each one of the recurrence relationships (26) and (27) can be written as nþ1 ð5 þ 3r þ 5c1 Þdmnþ1 þ ð1 3r þ c1 Þdmþ1 ¼ ð5 3r 5c1 Þdnm
þ ð1 þ 3r c1 Þdnmþ1 þ 3kað1 þ ð1 dÞs1 Þ;
m ¼ 0;
nþ1 nþ1 ð1 3r þ c1 Þdm1 þ ð4 þ 6r þ 4c1 Þdmnþ1 þ ð1 3r þ c1 Þdmþ1
¼ ð1 þ 3r c1 Þdnm1 þ ð4 6r 4c1 Þdnm þ ð1 þ 3r c1 Þdnmþ1 þ 3kað1 þ ð1 dÞs1 Þ;
m ¼ 1ð1ÞN 2;
nþ1 ð1 3r þ c1 Þdm1 þ ðð4 þ kÞð1 þ c1 Þ þ 3rð2 kÞÞdmnþ1
¼ ð1 þ 3r c1 Þdnm1 þ ðð4 þ kÞð1 c1 Þ 3rð2 kÞÞdnm þ 3kað1 þ ð1 dÞs1 Þ;
m¼N 1
ð28Þ
154
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
and nþ1 ð5 þ 3r þ 5c2 Þdmnþ1 þ ð1 3r þ c2 Þdmþ1 ¼ ð5 3r 5c2 Þdnm
þ ð1 þ 3r c2 Þdnmþ1 þ 3kað2 dÞ 6kað1 dÞs2 ; ð1 3r þ
nþ1 c2 Þdm1
þ ð4 þ 6r þ
4c2 Þdmnþ1
m ¼ 0;
nþ1 þ ð1 3r þ c2 Þdmþ1
¼ ð1 þ 3r c2 Þdnm1 þ ð4 6r 4c2 Þdnm þ ð1 þ 3r c2 Þdnmþ1 þ 3kað2 dÞ 6kað1 dÞs2 ; ð1 3r þ
nþ1 c2 Þdm1
þ ðð4 þ kÞð1 þ c2 Þ þ 3rð2
m ¼ 1ð1ÞN 2;
kÞÞdmnþ1
¼ ð1 þ 3r c2 Þdnm1 þ ðð4 þ kÞð1 c2 Þ 3rð2 kÞÞdnm þ 3kað2 dÞ 6kað1 dÞs2 ;
m ¼ N 1; ð29Þ
respectively, where k ¼ ð2 hbÞ=ð2 þ hbÞ. Each of the above systems can be easily and efficiently solved by using the Thomas algorithm. A Fourier stability analysis shows that the schemes (28) and (29) are unconditionally stable since the Crank–Nicolson approach is used for the element parameters d. 4.2. Petrov–Galerkin (PG) method The weighted residual form (21) over each finite element ½xm ; xmþ1 , with integration by parts, leads to the weak form Z xmþ1 ðWUt þ Wx Ux W arÞ dx ¼ bW ðxN ÞU ðxN ; tÞ: ð30Þ xm
The weight functions Wm are taken as linear B-spline basis functions for the finite element ½xm ; xmþ1 of the form [17] Wm ¼ 1 n=h;
Wmþ1 ¼ n=h:
Substituting the electrical conductivities (8) and (9) into Eq. (30), after a simple arrangement we respectively obtain the following equations: Z xmþ1 ðWUt þ Wx Ux þ að1 dÞs1 WU að1 þ ð1 dÞs1 ÞW Þ dx xm
¼ bW ðxN ÞU ðxN ; tÞ
ð31Þ
and Z
xmþ1
ðWUt þ Wx Ux þ að1 dÞð1 s2 ÞWU aðð2 dÞ 2ð1 dÞs2 ÞW Þ dx
xm
¼ bW ðxN ÞU ðxN ; tÞ:
ð32Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
155
Inserting Eqs. (19) and (20) into Eqs. (31) and (32) lead to mþ1 Z xmþ1 mþ1 Z xmþ1 X X e 0 0 _ Wi Qj dx dj þ Wi Qj dx dej xm
j¼m1
xm
j¼m1
þ að1 dÞs1
mþ1 Z X j¼m1
¼ að1 þ ð1 dÞs1 Þ
xmþ1
Wi Qj dx dej þ b
xm
Z
Wi ðxN ÞQj ðxN Þdej
j¼m1
xmþ1
mþ1 X
Wi dx xm
and mþ1 Z X
mþ1 Z X e _ Wi Qj dx dj þ
xm
j¼m1
xmþ1
mþ1 X j¼m1
xmþ1
0
Wi Q0j dx
dej þ að1 dÞð1 s2 Þ
xm
j¼m1
Z
xmþ1
mþ1 X Wi Qj dx dej þ b Wi ðxN ÞQj ðxN Þdej
xm
¼ aðð2 dÞ 2ð1 dÞs2 Þ
Z
j¼m1 xmþ1
Wi dx ; xm
which can also be written in matrix forms as Ae d_ e þ ðc3 Ae þ Be þ C e Þde ¼ ða þ c3 Þd e
ð33Þ
Ae d_ e þ ðc4 Ae þ Be þ C e Þde ¼ ðað1 ð1 dÞs2 Þ þ c4 Þd e ;
ð34Þ
and
respectively, where c3 ¼ að1 dÞs1 , c4 ¼ að1 dÞð1 s2 Þ, de ¼ ðdm1 ; dm ; dmþ1 ÞT are the element parameters, and Ae ; Be ; C e and d e are the element matrices given by the following integrals: Z xmþ1 Aeij ¼ Wi Qj dx; xm Z xmþ1 Wi 0 Q0j dx; Beij ¼ ð35Þ xm e Cij ¼ bWi ðxN ÞQj ðxN Þ; Z xmþ1 die ¼ Wi dx; xm
where i ¼ m; m þ 1 and j ¼ m 1; m; m þ 1. The element matrices (35) are obtained as h 3 8 1 1 1 0 1 h 1 Ae ¼ ; Be ¼ ; de ¼ 12 1 8 3 h 1 0 1 2 1
156
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
and C e is a zero matrix but 0 0 0 : C N 1 ¼ 0 b b Assembling contributions from all elements, Eqs. (33) and (34) respectively yield the following matrix equations: Ad_ þ ðc3 A þ B þ CÞd ¼ ða þ c3 Þd
ð36Þ
Ad_ þ ðc4 A þ B þ CÞd ¼ ðað1 ð1 dÞs2 Þ þ c4 Þd;
ð37Þ
and
T
where d ¼ ðd1 ; d0 ; . . . ; dN Þ and A; B, C are ðN þ 1Þ ðN þ 2Þ matrices, and d is an ðN þ 1Þ 1 matrix which are derived from the corresponding element matrices Ae ; Be ; C e and d e . Substituting the approximations (25) into Eqs. (36) and (37) lead to equations ½ð2 þ c3 kÞA þ k ð B þ CÞdnþ1 ¼ ½ð2 c3 kÞA k ð B þ CÞdn þ 2k ða þ c3 Þd ð38Þ and ½ð2 þ c4 kÞA þ kðB þ CÞdnþ1 ¼ ½ð2 c4 kÞA kðB þ CÞdn þ 2kðað1 ð1 dÞs2 Þ þ c4 Þd;
ð39Þ
respectively. Utilizing one of the boundary conditions (2) in Eqs. (38) and (39) yield an ðN þ 1Þ ðN þ 1Þ pentadiagonal matrix system and so is easily and efficiently solved with a variant of the Thomas algorithm. Again a Fourier stability analysis shows that schemes (38) and (39) are unconditionally stable since the Crank–Nicolson approach is used for the element parameters de .
5. Numerical results and conclusion In order to obtain the steady-state temperature distribution U ðx; tÞ and the locations of the free boundaries s1 and s2 , respectively, in the warm/cold and hot/warm phases constructed in the previous sections we used criteria 1 max U0new U0old ; UNnew UNold 6 0:5 108 k with k ¼ 0:01 for both the SC and PG solutions. In order to show how good the numerical solutions of the problem in comparison with the exact ones we used the error norm L2 defined by
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
L2 ¼
N 2 1 X e ðxm Þ U ðxm Þ U N þ 1 m¼0
157
!1=2 :
To compare our results obtained by using the SC and PG methods with the previous numerical results of Kutluay and Wood [12] using the EFD method we used values for a ¼ 0:5; 102 with b ¼ 0:5 satisfying inequalities (15) and (18) to get numerical solutions in the warm/cold and hot/warm phases, respectively. All the numerical computations were executed on a Pentium 4 PC in the Fortran code using double-precision arithmetic to reduce the round-off error to a minimum. The SC solutions to the warm/cold phase (for a ¼ 0:5 and b ¼ 0:5) are obtained by using Eq. (28) and are given in Table 1 with the norm L2 for a variety values of the mesh size h. The norm L2 decreases from 0.319602 · 104 (for h ¼ 0:1) to 0.020591 · 104 (for h ¼ 0:025). It is clearly seen that the predicted temperature values converge to the exact ones as the mesh size h is refined. Table 2 displays the PG solutions of the warm/cold phase obtained by using Eq. (38) with the norm L2 . It is seen from the table that the agreement of the numerical solutions with the exact one is remarkably good. The results obtained by SC and PG methods are presented in Table 3 with the corresponding phase solutions given by Kutluay and Wood [12]. The solutions obtained by the present two methods (for h ¼ 0:05) coincide to five significant digits in the most cases and agree with those of Kutluay and Wood [12]. Table 4 shows a comparison of the predicted motion (location and speed) of the free boundary s1 obtained by virtue of Eq. (10) for k ¼ 0:01 and various mesh sizes. It is seen that the estimated values of s1 and s_ 1 obtained by the present methods Table 1 Values of the temperature distribution as predicted by the numerical (SC method) and exact solutions in the warm/cold phase x
Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.156290 1.154116 1.147582 1.136662 1.121311 1.101463 1.077037 1.047929 1.014018 0.975162 0.931199
1.156191 1.154017 1.147486 1.136570 1.121224 1.101383 1.076965 1.047867 1.013967 0.975123 0.931173
1.156170 1.153996 1.147466 1.136551 1.121206 1.101366 1.076950 1.047854 1.013956 0.975116 0.931169
L2 104
0.319602
0.063870
0.020591
1.156154 1.153980 1.147450 1.136535 1.121191 1.101353 1.076937 1.047843 1.013946 0.975107 0.931162
158
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
Table 2 Values of the temperature distribution as predicted by the numerical (PG method) and exact solutions in the warm/cold phase Numerical
x
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.156242 1.154067 1.147536 1.136619 1.121273 1.101431 1.077012 1.047913 1.014013 0.975170 0.931221
1.156179 1.154005 1.147474 1.136559 1.121214 1.101375 1.076958 1.047863 1.013965 0.975125 0.931178
1.156168 1.153994 1.147463 1.136549 1.121204 1.101365 1.076949 1.047854 1.013957 0.975117 0.931171
L2 104
0.231638
0.047248
0.018844
1.156154 1.153980 1.147450 1.136535 1.121191 1.101353 1.076937 1.047843 1.013946 0.975107 0.931162
Table 3 Comparison of numerical (SC and PG) solutions of the warm/cold phase with results from [12] for h ¼ 0:05 and k ¼ 0:01 x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Present
[12]
SC
PG
EFD
Exact
1.156191 1.154017 1.147486 1.136570 1.121224 1.101383 1.076965 1.047867 1.013967 0.975123 0.931173
1.156179 1.154005 1.147474 1.136559 1.121214 1.101375 1.076958 1.047863 1.013965 0.975125 0.931178
1.144337 1.142197 1.135767 1.125014 1.109885 1.090304 1.066174 1.037373 1.003757 0.965177 0.921597
1.145311 1.143173 1.136749 1.126008 1.110894 1.091334 1.067228 1.038456 1.004874 0.966255 0.922624
converge to the exact ones as the mesh size h is refined. For h ¼ 0:025, the norm L2 is 0.020591 · 104 in the SC method whereas it is 0.018844 · 104 in the PG method. Observe that the PG method needs less steady-state time ðt1 Þ and generates slightly more accurate solutions than the SC method. The SC solutions to the hot/warm phase (for a ¼ 102 and b ¼ 0:5) are obtained by using Eq. (29) and are displayed in Table 5 with the norm L2 for various values of h. The norm L2 decreases from 4.725202 ·104 (for h ¼ 0:1) to 0.125614 · 104 (for h ¼ 0:025). It is seen that the SC solutions are in good agreement with the exact solution and better accuracy can be achieved for smaller values of h. Table 6 presents the PG solutions of the hot/warm phase
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
159
Table 4 Comparison of the predicted location and speed of the free boundary s1 for k ¼ 0:01 and various mesh sizes at steady-state time ðt1 Þ Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
SC
s1 s_ 1 107 t1
0.836077 0.112404 22.17
0.837157 0.115755 22.13
0.837342 0.114231 22.11
0.837581 0.0 –
PG
s1 s_ 1 107 t1
0.836076 2.290229 19.29
0.837156 2.178131 19.35
0.837344 1.802939 19.63
0.837581 0.0 –
Table 5 Values of the temperature distribution as predicted by the numerical (SC method) and exact solutions in the hot/warm phase (e1 ¼ 4 107 , e2 ¼ 3 107 ) x
Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
2.000003 1.999998 1.999980 1.999926 1.999779 1.999372 1.998250 1.995160 1.986642 1.963170 1.898485
2.000001 1.999996 1.999975 1.999916 1.999757 1.999328 1.998168 1.995038 1.986590 1.963790 1.902255
2 þ e1 1.999995 1.999974 1.999914 1.999753 1.999319 1.998153 1.995021 1.986605 1.963993 1.903238
L2 104
4.725202
0.727451
0.125614
2 þ e2 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903601
obtained by using Eq. (39) with the norm L2 . We can see from the table that the agreement of the PG solutions with the exact solution is good enough. The results obtained by SC and PG methods are shown in Table 7 with the corresponding phase solutions given by Kutluay and Wood [12]. They are all in very good agreement with each other. Table 8 presents a comparison of the steady-state values of the numerical motion (location and speed) of the free boundary s2 obtained by virtue of Eq. (11) for k ¼ 0:01 and various mesh sizes. It is seen from Table 8 that the predicted values of s2 and s_ 2 obtained by SC and PG methods are reasonably in good agreement and also exhibit the expected convergence as h is refined. For h ¼ 0:025, the norm L2 is 0.125614 · 104 in the SC method whereas it is 0.000815 · 104 in the PG method. It is clearly seen that the PG method produces slightly more accurate solutions than the SC
160
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
Table 6 Values of the temperature distribution as predicted by the numerical (PG method) and exact solutions in the hot/warm phase (e1 ¼ 4 107 ; e2 ¼ 3 107 ) Numerical
x
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
2.000001 1.999996 1.999976 1.999918 1.999763 1.999342 1.998208 1.995143 1.986870 1.964535 1.904234
2.000001 1.999995 1.999974 1.999916 1.999756 1.999327 1.998173 1.995067 1.986715 1.964248 1.903815
2+e1 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964089 1.903602
L2 104
0.748935
0.178411
0.000815
2+e2 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903601
Table 7 Comparison of numerical (SC and PG) solutions of the hot/warm phase with results from [12] for h ¼ 0:05 and k ¼ 0:01 x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Present
[12]
SC
PG
EFD
Exact
2.000001 1.999996 1.999975 1.999916 1.999757 1.999328 1.998168 1.995038 1.986590 1.963790 1.902255
2.000001 1.999995 1.999974 1.999916 1.999756 1.999327 1.998173 1.995067 1.986715 1.964248 1.903815
2.000008 2.000003 1.999985 1.999929 1.999776 1.999364 1.998255 1.995271 1.987241 1.965634 1.907492
2.000002 1.999997 1.999978 1.999923 1.999774 1.999368 1.998266 1.995268 1.987121 1.964973 1.904771
method. However, the SC method needs less steady-state time ðt1 Þ than the PG method. In this paper, subdomain collocation and Petrov–Galerkin B-spline finite element methods were efficiently applied to the PTC thermistor problem with the electrical conductivities defined by expressions (8) and (9) in Section 1 in order to obtain steady-state temperature distributions and the motion (location and speed) of the free interfaces. It is shown that all the numerical results obtained by using SC and PG methods display the correct physical phenomena of the problem, and are reasonably in good agreement with each other and
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
161
Table 8 Comparison of the predicted location and speed of the free boundary s2 for k ¼ 0:01 and various mesh sizes at steady-state time ðt1 Þ Numerical
Exact
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
SC
s2 s_ 2 107 t1
0.057120 1.550098 0.96
0.034289 1.497480 1.40
0.028094 0.026169 4.63
0.025334 0.0 –
PG
s2 s_ 2 107 t1
0.017410 2.348426 0.68
0.021160 0.018321 2.01
0.025326 0.011836 6.91
0.025334 0.0 –
exact ones. In general, more accurate numerical solutions will result from smaller mesh size h. However, more computer storage and computations will be required for smaller mesh size. As a conclusion, the electrical conductivities proposed in this paper can be used as a good representation for r to model the PTC thermistor problem.
References [1] G. Cimatti, Remark on existence and uniqueness for the thermistor problem under mixed boundary conditions, Q. Appl. Math. 47 (1989) 117–121. [2] G. Cimatti, G. Prodi, Existence results for a non-linear elliptic system modelling a temperature dependent electrical resistor, Ann. Mat. Pura Appl. 162 (1992) 33–42. [3] A. Fowler, I. Frigaard, S.D. Howison, Temperature surges in thermistors, SIAM J. Appl. Math. 52 (1992) 998–1011. [4] S.D. Howison, J.F. Rodrigues, M. Shillor, Stationary solutions to the thermistor problem, J. Math. Anal. Appl. 174 (1993) 573–588. [5] J. Wiedmann, The thermistor problem, Nonlinear Differ. Equations Appl. 4 (1997) 133–148. [6] S. Zhou, D.R. Westbrook, Numerical solutions of the thermistor equations, J. Comput. Appl. Math. 79 (1997) 101–118. [7] D.R. Westbrook, The thermistor: a problem in heat and current flow, Numer. Meth. PDEs 5 (1989) 259–273. [8] K. Preis, O. Biro, R. Dyczij-Edlinger, K.R. Richter, Zs. Badics, H. Riedler, H. St€ onger, Application of FEM to coupled electric, thermal and mechanical problems, IEEE Trans. Magn. 30 (5) (1994) 3316–3319. [9] S. Kutluay, A. Esen, A B-spline finite element method for the thermistor with the modified electrical conductivity, Appl. Math. Comput., AMC 8384, in press. [10] S. Kutluay, A.R. Bahadir, A. Ozdes, A variety of finite difference methods to the thermistor with a new modified electrical conductivity, Appl. Math. Comput. 106 (1999) 205–213. [11] A.S. Wood, Modelling the thermistor, in: M. Heili€ o (Ed.), Proceedings of 5th European Conference on Maths. Industry, Kluwer Academic, Stuttgart, 1991, pp. 397–400. [12] S. Kutluay, A.S. Wood, Numerical solutions of the thermistor problem with a ramp electrical conductivity, Appl. Math. Comput. 148 (2004) 145–162. [13] S. Kutluay, A.R. Bahadir, A. Ozdes, Various methods to the thermistor problem with a bulk electrical conductivity, Int. J. Numer. Meth. Eng. 45 (1999) 1–12.
162
S. Kutluay, A. Esen / Appl. Math. Comput. 163 (2005) 147–162
[14] S. Kutluay, A.S. Wood, A. Esen, A heat balance integral approach to the thermistor problem with a modified electrical conductivity, submitted for publication. [15] A.S. Wood, S. Kutluay, A heat balance integral model of the thermistor, Int. J. Heat Mass Transfer 38 (10) (1995) 1831–1840. [16] S. Howison, Complex variables in industrial mathematics, in: H. Neunzert (Ed.), Proceedings of 2nd European Symposium Maths. Industry, Oberwolfach, Kluwer Academic, Stuttgart, 1988, pp. 153–166. [17] P.M. Prenter, Splines and Variational Methods, John-Wiley, New York, 1975. [18] O.C. Zienkiewics, R.L. Taylor, The Finite Element Method, Butterworth-Heinemann, Oxford, 2002. [19] G.D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Clarendon, Oxford, 1987.