- Email: [email protected]
Finite element solution of the thermistor problem with a ramp electrical conductivity
Applied Mathematics and Computation 161 (2005) 897–913 www.elsevier.com/locate/amc
Finite element solution of the thermistor problem with a ramp electrical conductivity S. Kutluay *, A. Esen Department of Mathematics, Faculty of Arts and Science, In€on€u University, 44069 Malatya, Turkey
Abstract This paper presents approximate steady-state solutions of a one-dimensional positive temperature coefficient thermistor problem with a ramp function electrical conductivity using the Galerkin cubic B-spline finite element method. The numerical results obtained by the present method have been compared with the exact one and also those obtained by earlier authors, and are found to be in very good agreement with each other. Further, it is shown that the numerical solutions satisfy the correct physical phenomena of the problem. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Thermistor problem; Ramp electrical conductivity; Galerkin finite element; Cubic Bsplines
1. Introduction This paper deals with an initial-boundary problem consisting of a heat and current flow in a circuit device called a thermistor which is used a switching device in electronic circuits. A typical thermistor is a cylinder made from a ceramic material whose electrical conductivity is a function of the temperature.
*
Corresponding author. 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.2003.12.059
898
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
Let U ðx; tÞ and /ðx; tÞ be the temperature and electrical potential, respectively, and rðU Þ the electrical conductivity. Then, the thermistor problem is modelled in a dimensionless form as follows. Heat-flow problem: 2 oU o2 U o/ ¼ 2 þ ar ; 0 < x < 1; t > 0; ð1Þ ot ox ox oU ¼ 0; ox
x ¼ 0; t > 0;
oU þ bU ¼ 0; ox U ðx; 0Þ ¼ 0;
ð2Þ
x ¼ 1; t > 0;
ð3Þ ð4Þ
0 6 x 6 1:
Current-flow problem: o o/ r ¼ 0; 0 < x < 1; t > 0; ox ox
ð5Þ
/ð0; tÞ ¼ 0; /ð1; tÞ ¼ 1;
t > 0; t > 0;
ð6Þ
/ðx; 0Þ ¼ x;
0 6 x 6 1:
ð7Þ
In the above equations, a is the ratio of the rate of heat production, b is a positive surface heat transfer coefficient in the thermistor and r, the electrical conductivity, is of the form [1] 8 0 6 U 6 1; < 1; rðU Þ ¼ 1 þ ð1 U Þð1 dÞ; 1 < U < 2; ð8Þ : d 105 ; U P2 which is mathematically equivalent to a ramp function, as illustrated in Fig. 1. σ 1
δ 0
1
2
Fig. 1. Variation of r with U .
U
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
899
The above constructed model problem is known as a positive temperature coefficient (PTC) thermistor since its conductivity r given by Eq. (8) decreases with an increase of temperature U . The variation of r with U makes the thermistor useful in several areas such as thermometry, the measurement of microwave-frequency power, thermal relays and control devices activated by changes in temperature such as a switching device or as a constant temperature heater [2–4]. The PTC thermistor problem under various assumptions on r has been considered by many authors in the steady-state case. We refer the reader to [5–9] and references therein where various aspects of modeling, existence, uniqueness and behaviour of solutions are presented. A variety of techniques including finite difference [10–12], finite element [12–15] and approximate analytical methods [16] have been applied to the PTC thermistor problem particularly with a step function electrical conductivity. More recently, Kutluay and Wood [1] presented approximate steady-state solutions to the problem with a ramp electrical conductivity using the standard explicit finite difference method. They showed that a ramp electrical conductivity seems to be satisfactory in explaining the physical properties of the problem. In this paper we consider the one-dimensional PTC thermistor problem with a ramp electrical conductivity and use the Galerkin finite element (GFE) method based on the cubic B-spline functions in order to determine the evolving temperature distributions and the position of the free boundary during the steady-state case. Further, the numerical results obtained by the GFE have been compared with those obtained by the explicit finite difference (EFD) method in [1].
2. Determination of the particular phases Since U ðx; 0Þ ¼ 0 and b > 0 in the thermistor problem, we assume monotonicity of the temperature profile then the point x ¼ 0 will always be the hottest and will the first point to reach a value of 1 and 2 above each of which the behaviour of r changes. Owing to the decrease in r, the rate of heat loss at x ¼ 1 will be equal to the internal heat generation and so a steady-state (or equilibrium) temperature will be reached. The steady-state may be any one of the following five phases as seen in Fig. 2. Cold phase ð0 6 t 6 t0 Þ: In this phase, 0 6 U ðx; tÞ 6 1 and so rðU Þ ¼ 1. The physical requirements U ð0; tÞ P U ð1; tÞ for all t 2 ½0; t0 and U ð0; tÞ ¼ 1 at t ¼ t0 should be satisfied. Warm/cold phase ðt0 < t 6 t1 Þ: In this phase, U ðx; tÞ > 1 and so rðU Þ ¼ 1 þ ð1 U Þð1 dÞ. U ðx; tÞ 6 1 and so rðU Þ ¼ 1. The physical requirements U ð0; tÞ P U ð1; tÞ for all t 2 ðt0 ; t1 and U ð1; tÞ ¼ 1 at t ¼ t1 should be satisfied.
900
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
U
U
U
2
2
2
1
1
1
1
0
x
0
s1
x
1
U
U
2
2
1
1 s2
0
(d)
1
x
(c)
(b)
(a)
0
1
x
0
1
x
(e)
Fig. 2. Steady-state configurations: (a) cold, (b) warm/cold, (c) warm, (d) hot/warm and (e) hot.
Further, this phase has a free boundary s1 , on which U ðx; tÞ ¼ 1, separating a warm region ðU > 1Þ from a cold region ðU 6 1Þ. Warm phase ðt1 < t 6 t2 Þ: In this phase, 1 < U ðx; tÞ < 2 and so rðU Þ ¼ 1 þ ð1 U Þð1 dÞ. The physical requirement U ð0; tÞ P U ð1; tÞ for all t 2 ðt1 ; t2 should be satisfied. Hot/warm phase ðt2 < t 6 t3 Þ: In this phase, U ðx; tÞ P 2 and so r ¼ d. U ðx; tÞ < 2 and so r ¼ 1 þ ð1 U Þð1 dÞ. The physical requirements U ð0; tÞ P U ð1; tÞ for all t 2 ðt2 ; t3 and U ð1; tÞ ¼ 2 at t ¼ t3 should be satisfied. Further, this phase has a free boundary s2 , on which U ðx; tÞ ¼ 2, separating a hot region ðU > 2Þ from a warm region ðU 6 2Þ. Hot phase ðt > t3 Þ: In this phase, U ðx; tÞ P 2 and so r ¼ d. The physical requirement U ð0; tÞ P U ð1; tÞ > 2 for all t > t3 should be satisfied. It is seen from Eq. (8) that the conductivity r in the warm/cold decreases from 1 at the temperature 0 6 U 6 1 to 1 þ ð1 U Þð1 dÞ at the temperature 1 < U < 2, in the hot/warm phase decreases from 1 þ ð1 U Þð1 dÞ at the temperature 1 < U < 2 to d at the temperature U P 2. This decrease results in an oscillation in the temperatures when the GFE method is applied to the above mentioned two phases. It is not able to be recovered during the process since neither the temperature domain in the cold and warm regions of the warm/cold phase nor the temperature domain in the hot and warm regions of the hot/warm phase are symmetrical. In order to avoid the non-diminishing oscillations in the numerical computations, in this study we define the conductivities
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
rðU Þ ¼ 1 þ ð1 U Þð1 dÞs1
901
ð9Þ
and rðU Þ ¼ 1 þ ð1 U Þð1 dÞ þ ½d ð1 þ ð1 U Þð1 dÞÞ s2
ð10Þ
for warm/cold and hot/warm phases, respectively. Each of (9) and (10) is obtained by using a linear interpolation technique. Eqs. (9) and (10) can be used to determine the conductivity of each particular phase regarding the values of s1 and s2 . For example, in the cold phase, s1 ¼ 0 and so Eq. (9) gives r ¼ 1 which is exactly the same with the conductivity (8). In the warm/cold phase, s1 2 ð0; 1 (see Fig. 2(b)) and so the conductivity is as given in Eq. (9). In the warm phase, s1 ¼ 1 and so Eq. (9) gives r ¼ 1 þ ð1 U Þð1 dÞ which is also exactly the same with the conductivity (8). In the hot/warm phase, s2 2 ð0; 1 (see Fig. 2(d)) and so the conductivity is as given in Eq. (10). In hot phase, s2 ¼ 1 and so Eq. (10) gives r ¼ d which is again the same with the conductivity (8). In this paper, we assume that /ðx; tÞ ¼ x everywhere in 0 6 x 6 1 for all time t. The thermistor problem given by Eqs. (1)–(7) is then reduced to the heat conduction problem oU o2 U ¼ 2 þ ar; ot ox
0 < x < 1; t > 0;
ð11Þ
supplemented by the conditions (2)–(4).
3. Exact steady-state solutions The time partial derivative in Eq. (11) vanishes since the temperature does not vary with time t at steady-state. The partial differential equation is then reduced to an ordinary differential equation that can be easily solved one of standard analytical and numerical solution methods. Under the conductivity considered in the above section, the cold, warm and hot phases remain the same as in [1] but unlike the warm/cold and hot/warm phases will need to reformulate. Thus we will only examine exact steady-state solutions of the warm/cold phase using the conductivity (9) and hot/warm phase using the conductivity (10). In order to distinguish exact steady-state values from time dependent ones we will use the asterisk * as a superscript for exact steady-state values. Cold steady-state: The cold solution is [1] a 2 U ðxÞ ¼ 1 þ x2 ; 0 6 x 6 1 2 b
902
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
subject to the inequality a 6 2b=ð2 þ bÞ:
ð12Þ
Warm/cold steady-state: We easily obtain the exact steady-state solution to this phase as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 U ðxÞ ¼ 2c1 cosh að1 dÞs 1 x þ þ 1; ð13Þ ð1 dÞs 1 where pffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 be að1dÞs1 ðd1Þs
1 1 c1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi ffi : b að1 dÞs 1 þ ðb þ að1 dÞs 1 Þe2 að1dÞs1 Since U ðs 1 Þ ¼ 1, from Eq. (13) we obtain an equation for the interface s 1 : qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c1 ð1 dÞs 1 cosh s 1 að1 dÞs 1 þ 1 ¼ 0 ð14Þ which can be solved numerically. For s 1 ¼ 1, Eq. (14) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ : b¼ 1d Since s 1 2 ð0; 1 for all time, we obtain the values of a, b and d satisfying the inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ 6b ð15Þ 1d implies a warm/cold steady-state. The resulting inequality is identical to the sub-bound of b obtained by Kutluay and Wood [1] which is evident that Eq. (9) can be used as a good representation for the electrical conductivity of the warm/cold phase. Warm steady-state: The warm solution is [1] 1 U ðxÞ ¼ 1 þ 1d ! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b cosh að1 dÞx 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ sinh að1 dÞ þ b cosh að1 dÞ subject to satisfying the inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ að1 dÞ tanh að1 dÞ : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 b 6 1 d 2 d 1 þ cosh að1 dÞ
ð16Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
903
Hot/warm steady-state: We obtain the exact steady-state solution to this phase as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2s 2 U ðxÞ ¼ 2c2 cosh að1 dÞð1 s 2 Þx 1 s 2 2d ; ð17Þ þ ð1 dÞð1 s 2 Þ where
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2s 2d be að1dÞð1s2 Þ 1s2 ð1dÞð1s
Þ 2 2 c2 ¼ : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b að1 dÞð1 s2 Þ þ b þ að1 dÞð1 s 2 Þ e að1dÞð1s2 Þ Since U ðs 2 Þ ¼ 2, from Eq. (17) we obtain an equation for the interface s 2 : qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c2 ð1 dÞð1 s 2 Þ cosh s 2 að1 dÞð1 s 2 Þ þ d ¼ 0 ð18Þ which can be solved numerically. For s 2 ¼ 0, Eq. (18) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ b¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 d 1 þ cosh að1 dÞ Since s 2 2 ð0; 1 for all time, we obtain the values of a, b and d satisfying the inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ ð19Þ b< pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 d 1 þ cosh að1 dÞ for a hot/warm steady-state solution. The resulting inequality is identical to the upper bound of b obtained by Kutluay and Wood [1] which is evident that Eq. (10) can be used as a good representation for the electrical conductivity of the hot/warm phase. Hot steady-state: The hot solution is [1] ad 2 U ðxÞ ¼ 1 þ x2 2 b subject to the inequality 2b 6 ad:
ð20Þ
904
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
4. The method of solution In this section, we examine approximate steady-state solutions of each particular phase, having the ramp electrical conductivity, by using the Galerkin cubic B-spline finite element method. The interval 0 ¼ x0 < x1 < < xN 1 < xN ¼ 1 is divided into N equal finite elements by h ¼ 1=N . The cubic B-splines wm ðxÞ, at the knots xm which form a basis for functions over ½0; 1 , are defined as [17] 8 3 ðx xm2 Þ ; > > > 3 2 3 2 > 1 < h3 þ 3h2 ðx xm1 Þ þ 3hðx xm1 Þ2 3ðx xm1 Þ3 ; wm ðxÞ ¼ 3 h þ 3h ðxmþ1 xÞ þ 3hðxmþ1 xÞ 3ðxmþ1 xÞ ; h > > 3 > > : ðxmþ2 xÞ ; 0;
½xm2 ; xm1 ; ½xm1 ; xm ; ½xm ; xmþ1 ; ½xmþ1 ; xmþ2 ; otherwise:
ð21Þ Let UN ðx; tÞ be an approximation to the exact solution U ðx; tÞ. Then the global approximation UN ðx; tÞ can be written in terms of the splines as UN ðx; tÞ ¼
N þ1 X
wj ðxÞdj ðtÞ;
ð22Þ
j¼1
where dj are time dependent parameters to be determined from the boundary and weighted residual conditions. In terms of a local coordinate transformation n ¼ x xm ð0 6 n 6 hÞ, Eq. (21) can be expressed as 8 3 ðh nÞ > > wm1 < 3 2 1 3hðh nÞ2 3ðh nÞ3 ; 0 6 n 6 h: wm ¼ 3 h3 þ 3h2 ðh nÞ þ 2 3 h > > wmþ1 : h3 þ 3h n þ 3hn 3n n ð23Þ Since all other splines are 0 over the element ½xm ; xmþ1 , the approximation (22) over this element can be written in terms of the basis functions (23) as UN ¼
mþ2 X
w j dj :
ð24Þ
j¼m1
Using splines (23) and the approximation (24), the nodal values Um and Um0 at the node xm in terms of element parameters dm can be obtained as Um ¼ U ðxm Þ ¼ dm1 þ 4dm þ dmþ1 ; hUm0 ¼ hU 0 ðxm Þ ¼ 3ðdmþ1 dm1 Þ:
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
905
In the above equation and throughout this paper, the prime denotes derivative with respect to x. In the Galerkin method, basis functions (23) are used as weighting functions Wm . So, the weak form of the model problem given by Eq. (11) is Z
1
Wm ðUt Uxx arÞ dx ¼ 0:
ð25Þ
0
As mentioned in Section 2, it should be sufficient to work on the electrical conductivities given by Eqs. (9) and (10) only. Using the electrical conductivities (9) and (10) in Eq. (25), we respectively obtain the following integral equations over each finite element ½xm ; xmþ1 : Z xmþ1 ðWUt þ Wx Ux þ að1 dÞs1 WU að1 þ ð1 dÞs1 ÞW Þ dx xm
¼ bW ðxN ÞU ðxN ; tÞ
ð26Þ
and Z
xmþ1
ðWUt þ Wx Ux þ að1 dÞð1 s2 ÞWU aðð2 dÞ
xm
2ð1 dÞWs2 ÞÞ dx ¼ bW ðxN ÞU ðxN ; tÞ:
ð27Þ
Inserting Eqs. (23) and (24) into Eqs. (26) and (27) lead to mþ2 Z X j¼m1
xmþ1
mþ2 Z X e _ wi wj dx dj þ
xm mþ2 Z X
¼ að1 þ ð1 dÞs1 Þ
dej
mþ2 X wi wj dx dej þ b wi ðxN Þwj ðxN Þdej
xm
j¼m1
Z
xmþ1
w0i w0j dx
xm
j¼m1
þ að1 dÞs1
xmþ1
j¼m1
xmþ1
wi dx xm
and mþ2 Z X j¼m1
xmþ1
mþ2 Z X e _ wi wj dx dj þ
xm
þ að1 dÞð1 s2 Þ
j¼m1
¼ aðð2 dÞ 2ð1 dÞs2 Þ
Z
xmþ1
xm xmþ1
xm
w0i w0j dx
dej
xm
j¼m1 mþ2 Z X
xmþ1
mþ2 X wi wj dx dej þ b wi ðxN Þwj ðxN Þdej
wi dx ;
j¼m1
906
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
which are also be written in matrix forms as Ae d_ e þ ðc1 Ae þ Be þ C e Þde ¼ ða þ c1 Þd e
ð28Þ
and Ae d_ e þ ðc2 Ae þ Be þ C e Þde ¼ ðað1 ð1 dÞs2 Þ þ c2 Þd e ;
ð29Þ
respectively, where c1 ¼ að1 dÞs1 and c2 ¼ að1 dÞð1 s2 Þ. In the above equations and throughout this paper, the dot denotes derivative with respect to t, de ¼ ðdm1 ; dm ; dmþ1 ; dmþ2 Þ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 wj dx; xm Z xmþ1 w0i w0j dx; Beij ¼ ð30Þ xm Cije ¼ bwi ðxN Þwj ðxN Þ; Z xmþ1 die ¼ wi dx; xm
where i; j ¼ m 1; m; m þ 1; m þ 2. After a simple arrangement, the element matrices (30) can be obtained as 2 3 20 129 60 1 60 7 h 6 6 129 1188 933 7 Aeij ¼ 6 7; 140 4 60 933 1188 129 5 1 18 1 6 6 21 Beij ¼ 6 10h 4 36
60 21 102
129 36 87
87
102
3
36
21
2
20 3 3 36 7 7 7; 21 5
2
3 1 7 h6 6 11 7 die ¼ 6 7 4 4 11 5
18
1
and Cije is a zero matrix but 2 3 0 0 0 0 6 0 b 4b b 7 7 CijN 1 ¼ 6 4 0 4b 16b 4b 5: 0 b 4b b Assembling all contributions from all elements, Eqs. (28) and (29) lead to matrix differential equations Ad_ þ ðc1 A þ B þ CÞd ¼ ða þ c1 Þd ð31Þ and Ad_ þ ðc2 A þ B þ CÞd ¼ ðað1 ð1 dÞs2 Þ þ c2 Þd;
ð32Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
907
T
respectively, where d ¼ ðd1 ; d0 ; . . . ; dN þ1 Þ and A, B and C are ðN þ 3Þ ðN þ 3Þ matrices and d is an ðN þ 3Þ 1 matrix which are defined from the element matrices Ae , Be , C e and d e , respectively. Substituting the Crank–Nicolson approach d ¼ ðdn þ dnþ1 Þ=2 and the forward finite difference d_ ¼ ðdnþ1 dn Þ=k into Eqs. (31) and (32) yields the system of equations ½ð2 þ c1 kÞA þ kðB þ CÞ dnþ1 ¼ ½ð2 c1 kÞA kðB þ CÞ dn þ 2kða þ c1 Þd ð33Þ and ½ð2 þ c2 kÞA þ kðB þ CÞ dnþ1 ¼ ½ð2 c2 kÞA kðB þ CÞ dn þ 2kðað1 ð1 dÞs2 Þ þ c2 Þd;
ð34Þ
respectively, where kð DtÞ is the time step. Each resulting system is an ðN þ 3Þ ðN þ 3Þ septadiagonal equation system and so is easily and efficiently solved with a variant of the Thomas algorithm. A Fourier stability analysis shows the schemes (33) and (34) to be unconditionally stable since the Crank–Nicolson approach is used for the element parameters de .
5. Numerical results and conclusion To make a comparison with earlier work [1] we consistently use a choice of values for a (0.1, 0.5, 1, 102 and 105 ) with b ¼ 0:5 satisfying inequalities (12), (15), (16), (19) and (20) to obtain steady-state solutions for the cold, warm/ cold, warm, hot/warm and hot phases, respectively. We use the steady-state criteria 1 maxfjU0new U0old j; jUNnew UNold jg 6 0:5 108 k to get the Galerkin finite element solution of each particular phase. All calculations were performed on a Pentium 4 processor using Microsoft Fortran Compiler. In order to show how good the numerical solutions of the problem in comparison with the exact ones we shall use the error norm L2 defined by " #1=2 N 1 X 2 L2 ¼ jU ðxm Þ U ðxm Þj : N þ 1 m¼0 Cold steady-state solutions: The values a ¼ 0:1 and b ¼ 0:5 satisfy inequality (12) and so give a cold steady-state. Replacing s1 by 0, Eq. (33) is solved using the Thomas algorithm. The numerical results obtained by the GFE method
908
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
Table 1 Steady-state solutions in the cold phase x
GFE, [1] and exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.2500 0.2495 0.2480 0.2455 0.2420 0.2375 0.2320 0.2255 0.2180 0.2095 0.2000 t1 (GFE) ¼ 39.56, t1 ([1]) ¼ 39.54
(for k ¼ 0:01, h ¼ 0:05) are compared with exact one and also those obtained by the EFD method in [1] as shown in Table 1. It is seen that they are all in very good agreement with each other. Warm/cold steady-state solutions: The values a ¼ 0:5 and b ¼ 0:5 satisfy inequality (15) and then give a warm/cold steady-state. This phase has a free boundary s1 ðtÞ 2 ð0; 1 (see Fig. 2(b)) which is an unknown function of time and must be determined as a part of solution. If the predicted temperature is U ðxm Þ > 1 and U ðxmþ1 Þ < 1 at a time level n, the boundary s1 is between adjacent values U ðxm Þ and U ðxmþ1 Þ. Since U ðs1 ; tÞ ¼ 1, the location of the free boundary s1 can be estimated by using a linear interpolation as [1] sn1 ¼ xm þ
hð1 U ðxm ÞÞ : U ðxmþ1 Þ U ðxm Þ
ð35Þ
Thus, the motion of sðtÞ can be obtained from s_ 1 ¼
s1nþ1 sn1 : Dt
Eq. (33) is again solved using the Thomas algorithm. The results obtained by the GFE method (for k ¼ 0:01, h ¼ 0:1, 0.05, 0.025, 0.0125) are presented in Table 2. It is observed that the numerical solutions are in good agreement with the exact one, and also exhibit the expected convergence as the mesh size h is refined. The results obtained by the GFE method (for k ¼ 0:01, h ¼ 0:05) are also compared with those obtained by the EFD method [1]. It is seen that they agree with each other, as shown in Table 3. The steady-state value of the interface s1 (Eq. (35)) obtained by the GFE method using Eq. (33) is about 0.84 in some 22.25 units of steady-state time
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
909
Table 2 Steady-state solutions in the warm/cold phase
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:0125
Exact solution
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.076959 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
1.156163 1.153989 1.147458 1.136544 1.121199 1.101360 1.076945 1.047849 1.013953 0.975113 0.931167
1.156154 1.153980 1.147450 1.136535 1.121191 1.101353 1.076937 1.047843 1.013946 0.975107 0.931162
L2 105
0.231657
0.047376
0.018918
0.008275
x
Numerical solution
Table 3 Comparison of steady-state solutions in the warm/cold phase with results from [1] at h ¼ 0:05 and k ¼ 0:01 x
GFE
[1] EFD
Exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.156179 1.154005 1.147474 1.136559 1.121214 1.101375 1.076959 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
t1 ¼ 22:25
t1 ¼ 23:82
ðt1 Þ whereas it is obtained by the EFD method in [1] as about 0.81 in something 23.82 units of steady-state time. Table 4 exhibits the approximate and exact values of s1 and s_ 1 with the steady-state times at various mesh sizes. It is clearly seen that as the mesh size h is reduced the predicted steady-state values of s1 and s_ 1 smoothly and rapidly reach their exact steady-state values. Warm steady-state solutions: A warm steady-state solution can be easily obtained by using the values a ¼ 1 and b ¼ 0:5. Replacing s1 by 1, Eq. (33) is again solved using the Thomas algorithm. The results obtained by the GFE method (for k ¼ 0:01, h ¼ 0:1, 0.05) are displayed in Table 5. There is a good
910
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
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 times ðt1 Þ
s1 s_ 1 107 t1
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:0125
Exact
0.836076 0.596046 22.23
0.837156 0.596024 22.25
0.837344 0.596006 22.27
0.837435 0.596002 32.23
0.837581 0.0
agreement between numerical and exact solutions. It is seen that GFE solutions exhibit the expected convergence as the mesh size h is reduced. As shown in Table 5, the numerical solutions (for h ¼ 0:05) hit the exact steady-state ones, and also agree with those given in [1]. Thus, it is not necessary to take smaller values of h. Hot/warm steady-state solutions: It is sufficient to use the values a ¼ 102 and b ¼ 0:5 to get a hot/warm steady-state solution. As shown in Fig. 2(d), this phase has a free boundary s2 2 ð0; 1 which is an unknown function of time and so must be determined as a part of solution with a similar manner to s1 . If the predicted temperature is U ðxm Þ > 2 and U ðxmþ1 Þ < 2 at a time level n, the boundary s2 is between adjacent values U ðxm Þ and U ðxmþ1 Þ. Since U ðs2 ; tÞ ¼ 2, the location of the free boundary s2 can be estimated by a linear interpolation as [1]
sn2 ¼ xm þ
hð2 U ðxm ÞÞ : U ðxmþ1 Þ U ðxm Þ
ð36Þ
Table 5 Comparison of steady-state solutions in the warm phase with results from [1] x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
GFE
EFD [1]
h ¼ 0:1
h ¼ 0:05
h ¼ 0:05
1.486325 1.483754 1.476017 1.463035 1.444679 1.420766 1.391055 1.355249 1.312991 1.263857 1.207354
1.486324 1.483754 1.476016 1.463035 1.444679 1.420765 1.391054 1.355249 1.312991 1.263856 1.207354
1.486363 1.483793 1.476058 1.463080 1.444729 1.420822 1.391120 1.355325 1.313080 1.263961 1.207476
t1 ¼ 13:94
t1 ¼ 14:35
Exact solution 1.486324 1.483754 1.476016 1.463035 1.444679 1.420765 1.391054 1.355249 1.312991 1.263856 1.207354
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
911
Thus, the motion of s2 ðtÞ can be easily obtained by using the usual forward difference formulae as s_ 2 ¼
s2nþ1 sn2 : Dt
Eq. (34) is solved using a variant of Thomas algorithm. The obtained GFE solutions (for k ¼ 0:01, h ¼ 0:1, 0.05, 0.025, 0.0125) are presented in Table 6. They are in good agreement with each other and also exact one. It is seen that numerical approximations show a high degree of accuracy when the mesh size h is reduced. Again, good agreement with the exact values is evident, as in convergence. Table 7 displays a comparison of approximate steady-state solutions obtained by the GFE method (for k ¼ 0:01, h ¼ 0:05) with those obtained by the EFD method in [1] which are reasonably in good agreement with each other. The steady-state values of the free boundary s2 (Eq. (36)) obtained by the GFE method using Eq. (34) is about 0.03 in something 6.94 units of steadystate time whereas it is obtained by using EFD method in [1] as about 0.12 in some 0.69 units of steady-state time. Table 8 displays the approximate and exact values of s2 and s_ 2 with the steady-state times at various mesh sizes. It is observed that as the mesh size h is reduced the predicted values of s2 and s_ 2 smoothly and rapidly approach to their exact values. Hot steady-state solutions: The values a ¼ 105 and b ¼ 0:5 is sufficient to obtain a hot steady-state solution. Replacing s2 by 1, Eq. (34) is solved using the Thomas algorithm. The GFE solutions (for k ¼ 0:01, h ¼ 0:05) are compared with EFD solutions (for h ¼ 0:05) given in [1] and also the exact solution, in Table 9. It is seen that they are all in full agreement with each other. Table 6 Steady-state solutions in the hot/warm phase ðe ¼ 3:0 107 Þ x
Numerical solution
Exact solution
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.999975 1.999918 1.999761 1.999338 1.998197 1.995120 1.986815 1.964410 1.904081
2.000001 1.999995 1.999974 1.999916 1.999756 1.999327 1.998172 1.995066 1.986712 1.964241 1.903806
2þe 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903602
L2 105
5.615767
1.710306
0.004061
2þe 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903601
912
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
Table 7 Comparison of steady-state solutions in the hot/warm phase with results from [1] at h ¼ 0:05 and k ¼ 0:01 x
GFE
[1] EFD
Exact
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.999995 1.999974 1.999916 1.999756 1.999327 1.998172 1.995066 1.986712 1.964241 1.903806
2.000008 2.000003 1.999985 1.999929 1.999976 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
t1 ¼ 6:94
t1 ¼ 0:69
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 times ðt1 Þ
s2 s_ 2 107 t1
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
Exact
0.015077 0.785326 3.36
0.020950 0.091881 6.94
0.025315 0.072008 23.26
0.025334 0.0
Table 9 Steady-state solutions in the hot phase x
GFE, [1] and exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
2.500 2.495 2.480 2.455 2.420 2.375 2.320 2.255 2.180 2.095 2.000 t1 (GFE) ¼ 44.96, t1 ([1]) ¼ 54.18
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
913
In this paper, the approximate steady-state solutions of the PTC thermistor problem with temperature dependent electrical conductivities described in Section 2 using the cubic GFE method have been presented. It is shown that the numerical results obtained by the GFE method exhibit the correct physical phenomena of the PTC thermistor problem. As a conclusion, we have demonstrated that the proposed electrical conductivities can be used effectively to obtain GFE solutions in transition phases (warm/cold and hot/warm) of the thermistor. References [1] S. Kutluay, A.S. Wood, Numerical solutions of the thermistor problem with a ramp electrical conductivity, Appl. Math. Comput. 148 (2004) 145–162. [2] R.P. Turner, ABCÕs of Thermistors, Foulsham-Sams, Slough, 1970. [3] F.J. Hyde, Thermistors, Butterworth (ILIFFE Books), London, 1971. [4] E.D. Macklen, Thermistors, Electrochemical Publications, Scotland, 1979. [5] G. Cimatti, Remark on existence and uniqueness for the thermistor problem under mixed boundary conditions, Q. Appl. Math. 47 (1989) 117–121. [6] 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. [7] A. Fowler, I. Frigaard, S.D. Howison, Temperature surges in thermistors, SIAM J. Appl. Math. 52 (1992) 998–1011. [8] S.D. Howison, J.F. Rodrigues, M. Shillor, Stationary solutions to the thermistor problem, J. Math. Anal. Appl. 174 (1993) 573–588. [9] J. Wiedmann, The thermistor problem, Nonlin. Differ. Equ. Appl. 4 (1997) 133–148. [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.), Proc. 5th Eur. Conf. Math. Industry, Kluwer Academic, Stuttgart, 1991, pp. 397–400. [12] S. Kutluay, A.R. Bahadir, A. Ozdes, Various methods to the thermistor problem with a bulk electrical conductivity, Int. J. Numer. Meth. Engng. 45 (1999) 1–12. [13] S. Kutluay, A. Esen, A B-spline finite element method for the thermistor with the modified electrical conductivity, Appl. Math. Comput. (in press). [14] D.R. Westbrook, The thermistor: a problem in heat and current flow, Numer. Meth. PDEs 5 (1989) 259–273. [15] 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. [16] A.S. Wood, S. Kutluay, A heat balance integral model of the thermistor, Int. J. Heat Mass Transfer 38 (10) (1995) 1831–1840. [17] P.M. Prenter, Splines and Variational Methods, John-Wiley, New York, 1975.
Finite element solution of the thermistor problem with a ramp electrical conductivity S. Kutluay *, A. Esen Department of Mathematics, Faculty of Arts and Science, In€on€u University, 44069 Malatya, Turkey
Abstract This paper presents approximate steady-state solutions of a one-dimensional positive temperature coefficient thermistor problem with a ramp function electrical conductivity using the Galerkin cubic B-spline finite element method. The numerical results obtained by the present method have been compared with the exact one and also those obtained by earlier authors, and are found to be in very good agreement with each other. Further, it is shown that the numerical solutions satisfy the correct physical phenomena of the problem. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Thermistor problem; Ramp electrical conductivity; Galerkin finite element; Cubic Bsplines
1. Introduction This paper deals with an initial-boundary problem consisting of a heat and current flow in a circuit device called a thermistor which is used a switching device in electronic circuits. A typical thermistor is a cylinder made from a ceramic material whose electrical conductivity is a function of the temperature.
*
Corresponding author. 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.2003.12.059
898
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
Let U ðx; tÞ and /ðx; tÞ be the temperature and electrical potential, respectively, and rðU Þ the electrical conductivity. Then, the thermistor problem is modelled in a dimensionless form as follows. Heat-flow problem: 2 oU o2 U o/ ¼ 2 þ ar ; 0 < x < 1; t > 0; ð1Þ ot ox ox oU ¼ 0; ox
x ¼ 0; t > 0;
oU þ bU ¼ 0; ox U ðx; 0Þ ¼ 0;
ð2Þ
x ¼ 1; t > 0;
ð3Þ ð4Þ
0 6 x 6 1:
Current-flow problem: o o/ r ¼ 0; 0 < x < 1; t > 0; ox ox
ð5Þ
/ð0; tÞ ¼ 0; /ð1; tÞ ¼ 1;
t > 0; t > 0;
ð6Þ
/ðx; 0Þ ¼ x;
0 6 x 6 1:
ð7Þ
In the above equations, a is the ratio of the rate of heat production, b is a positive surface heat transfer coefficient in the thermistor and r, the electrical conductivity, is of the form [1] 8 0 6 U 6 1; < 1; rðU Þ ¼ 1 þ ð1 U Þð1 dÞ; 1 < U < 2; ð8Þ : d 105 ; U P2 which is mathematically equivalent to a ramp function, as illustrated in Fig. 1. σ 1
δ 0
1
2
Fig. 1. Variation of r with U .
U
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
899
The above constructed model problem is known as a positive temperature coefficient (PTC) thermistor since its conductivity r given by Eq. (8) decreases with an increase of temperature U . The variation of r with U makes the thermistor useful in several areas such as thermometry, the measurement of microwave-frequency power, thermal relays and control devices activated by changes in temperature such as a switching device or as a constant temperature heater [2–4]. The PTC thermistor problem under various assumptions on r has been considered by many authors in the steady-state case. We refer the reader to [5–9] and references therein where various aspects of modeling, existence, uniqueness and behaviour of solutions are presented. A variety of techniques including finite difference [10–12], finite element [12–15] and approximate analytical methods [16] have been applied to the PTC thermistor problem particularly with a step function electrical conductivity. More recently, Kutluay and Wood [1] presented approximate steady-state solutions to the problem with a ramp electrical conductivity using the standard explicit finite difference method. They showed that a ramp electrical conductivity seems to be satisfactory in explaining the physical properties of the problem. In this paper we consider the one-dimensional PTC thermistor problem with a ramp electrical conductivity and use the Galerkin finite element (GFE) method based on the cubic B-spline functions in order to determine the evolving temperature distributions and the position of the free boundary during the steady-state case. Further, the numerical results obtained by the GFE have been compared with those obtained by the explicit finite difference (EFD) method in [1].
2. Determination of the particular phases Since U ðx; 0Þ ¼ 0 and b > 0 in the thermistor problem, we assume monotonicity of the temperature profile then the point x ¼ 0 will always be the hottest and will the first point to reach a value of 1 and 2 above each of which the behaviour of r changes. Owing to the decrease in r, the rate of heat loss at x ¼ 1 will be equal to the internal heat generation and so a steady-state (or equilibrium) temperature will be reached. The steady-state may be any one of the following five phases as seen in Fig. 2. Cold phase ð0 6 t 6 t0 Þ: In this phase, 0 6 U ðx; tÞ 6 1 and so rðU Þ ¼ 1. The physical requirements U ð0; tÞ P U ð1; tÞ for all t 2 ½0; t0 and U ð0; tÞ ¼ 1 at t ¼ t0 should be satisfied. Warm/cold phase ðt0 < t 6 t1 Þ: In this phase, U ðx; tÞ > 1 and so rðU Þ ¼ 1 þ ð1 U Þð1 dÞ. U ðx; tÞ 6 1 and so rðU Þ ¼ 1. The physical requirements U ð0; tÞ P U ð1; tÞ for all t 2 ðt0 ; t1 and U ð1; tÞ ¼ 1 at t ¼ t1 should be satisfied.
900
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
U
U
U
2
2
2
1
1
1
1
0
x
0
s1
x
1
U
U
2
2
1
1 s2
0
(d)
1
x
(c)
(b)
(a)
0
1
x
0
1
x
(e)
Fig. 2. Steady-state configurations: (a) cold, (b) warm/cold, (c) warm, (d) hot/warm and (e) hot.
Further, this phase has a free boundary s1 , on which U ðx; tÞ ¼ 1, separating a warm region ðU > 1Þ from a cold region ðU 6 1Þ. Warm phase ðt1 < t 6 t2 Þ: In this phase, 1 < U ðx; tÞ < 2 and so rðU Þ ¼ 1 þ ð1 U Þð1 dÞ. The physical requirement U ð0; tÞ P U ð1; tÞ for all t 2 ðt1 ; t2 should be satisfied. Hot/warm phase ðt2 < t 6 t3 Þ: In this phase, U ðx; tÞ P 2 and so r ¼ d. U ðx; tÞ < 2 and so r ¼ 1 þ ð1 U Þð1 dÞ. The physical requirements U ð0; tÞ P U ð1; tÞ for all t 2 ðt2 ; t3 and U ð1; tÞ ¼ 2 at t ¼ t3 should be satisfied. Further, this phase has a free boundary s2 , on which U ðx; tÞ ¼ 2, separating a hot region ðU > 2Þ from a warm region ðU 6 2Þ. Hot phase ðt > t3 Þ: In this phase, U ðx; tÞ P 2 and so r ¼ d. The physical requirement U ð0; tÞ P U ð1; tÞ > 2 for all t > t3 should be satisfied. It is seen from Eq. (8) that the conductivity r in the warm/cold decreases from 1 at the temperature 0 6 U 6 1 to 1 þ ð1 U Þð1 dÞ at the temperature 1 < U < 2, in the hot/warm phase decreases from 1 þ ð1 U Þð1 dÞ at the temperature 1 < U < 2 to d at the temperature U P 2. This decrease results in an oscillation in the temperatures when the GFE method is applied to the above mentioned two phases. It is not able to be recovered during the process since neither the temperature domain in the cold and warm regions of the warm/cold phase nor the temperature domain in the hot and warm regions of the hot/warm phase are symmetrical. In order to avoid the non-diminishing oscillations in the numerical computations, in this study we define the conductivities
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
rðU Þ ¼ 1 þ ð1 U Þð1 dÞs1
901
ð9Þ
and rðU Þ ¼ 1 þ ð1 U Þð1 dÞ þ ½d ð1 þ ð1 U Þð1 dÞÞ s2
ð10Þ
for warm/cold and hot/warm phases, respectively. Each of (9) and (10) is obtained by using a linear interpolation technique. Eqs. (9) and (10) can be used to determine the conductivity of each particular phase regarding the values of s1 and s2 . For example, in the cold phase, s1 ¼ 0 and so Eq. (9) gives r ¼ 1 which is exactly the same with the conductivity (8). In the warm/cold phase, s1 2 ð0; 1 (see Fig. 2(b)) and so the conductivity is as given in Eq. (9). In the warm phase, s1 ¼ 1 and so Eq. (9) gives r ¼ 1 þ ð1 U Þð1 dÞ which is also exactly the same with the conductivity (8). In the hot/warm phase, s2 2 ð0; 1 (see Fig. 2(d)) and so the conductivity is as given in Eq. (10). In hot phase, s2 ¼ 1 and so Eq. (10) gives r ¼ d which is again the same with the conductivity (8). In this paper, we assume that /ðx; tÞ ¼ x everywhere in 0 6 x 6 1 for all time t. The thermistor problem given by Eqs. (1)–(7) is then reduced to the heat conduction problem oU o2 U ¼ 2 þ ar; ot ox
0 < x < 1; t > 0;
ð11Þ
supplemented by the conditions (2)–(4).
3. Exact steady-state solutions The time partial derivative in Eq. (11) vanishes since the temperature does not vary with time t at steady-state. The partial differential equation is then reduced to an ordinary differential equation that can be easily solved one of standard analytical and numerical solution methods. Under the conductivity considered in the above section, the cold, warm and hot phases remain the same as in [1] but unlike the warm/cold and hot/warm phases will need to reformulate. Thus we will only examine exact steady-state solutions of the warm/cold phase using the conductivity (9) and hot/warm phase using the conductivity (10). In order to distinguish exact steady-state values from time dependent ones we will use the asterisk * as a superscript for exact steady-state values. Cold steady-state: The cold solution is [1] a 2 U ðxÞ ¼ 1 þ x2 ; 0 6 x 6 1 2 b
902
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
subject to the inequality a 6 2b=ð2 þ bÞ:
ð12Þ
Warm/cold steady-state: We easily obtain the exact steady-state solution to this phase as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 U ðxÞ ¼ 2c1 cosh að1 dÞs 1 x þ þ 1; ð13Þ ð1 dÞs 1 where pffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 be að1dÞs1 ðd1Þs
1 1 c1 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi ffi : b að1 dÞs 1 þ ðb þ að1 dÞs 1 Þe2 að1dÞs1 Since U ðs 1 Þ ¼ 1, from Eq. (13) we obtain an equation for the interface s 1 : qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c1 ð1 dÞs 1 cosh s 1 að1 dÞs 1 þ 1 ¼ 0 ð14Þ which can be solved numerically. For s 1 ¼ 1, Eq. (14) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ : b¼ 1d Since s 1 2 ð0; 1 for all time, we obtain the values of a, b and d satisfying the inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ tanh að1 dÞ 6b ð15Þ 1d implies a warm/cold steady-state. The resulting inequality is identical to the sub-bound of b obtained by Kutluay and Wood [1] which is evident that Eq. (9) can be used as a good representation for the electrical conductivity of the warm/cold phase. Warm steady-state: The warm solution is [1] 1 U ðxÞ ¼ 1 þ 1d ! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b cosh að1 dÞx 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1 dÞ sinh að1 dÞ þ b cosh að1 dÞ subject to satisfying the inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ að1 dÞ tanh að1 dÞ : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 b 6 1 d 2 d 1 þ cosh að1 dÞ
ð16Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
903
Hot/warm steady-state: We obtain the exact steady-state solution to this phase as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2s 2 U ðxÞ ¼ 2c2 cosh að1 dÞð1 s 2 Þx 1 s 2 2d ; ð17Þ þ ð1 dÞð1 s 2 Þ where
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2s 2d be að1dÞð1s2 Þ 1s2 ð1dÞð1s
Þ 2 2 c2 ¼ : pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b að1 dÞð1 s2 Þ þ b þ að1 dÞð1 s 2 Þ e að1dÞð1s2 Þ Since U ðs 2 Þ ¼ 2, from Eq. (17) we obtain an equation for the interface s 2 : qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2c2 ð1 dÞð1 s 2 Þ cosh s 2 að1 dÞð1 s 2 Þ þ d ¼ 0 ð18Þ which can be solved numerically. For s 2 ¼ 0, Eq. (18) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ b¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 d 1 þ cosh að1 dÞ Since s 2 2 ð0; 1 for all time, we obtain the values of a, b and d satisfying the inequality pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d að1 dÞ sinh að1 dÞ ð19Þ b< pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 d 1 þ cosh að1 dÞ for a hot/warm steady-state solution. The resulting inequality is identical to the upper bound of b obtained by Kutluay and Wood [1] which is evident that Eq. (10) can be used as a good representation for the electrical conductivity of the hot/warm phase. Hot steady-state: The hot solution is [1] ad 2 U ðxÞ ¼ 1 þ x2 2 b subject to the inequality 2b 6 ad:
ð20Þ
904
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
4. The method of solution In this section, we examine approximate steady-state solutions of each particular phase, having the ramp electrical conductivity, by using the Galerkin cubic B-spline finite element method. The interval 0 ¼ x0 < x1 < < xN 1 < xN ¼ 1 is divided into N equal finite elements by h ¼ 1=N . The cubic B-splines wm ðxÞ, at the knots xm which form a basis for functions over ½0; 1 , are defined as [17] 8 3 ðx xm2 Þ ; > > > 3 2 3 2 > 1 < h3 þ 3h2 ðx xm1 Þ þ 3hðx xm1 Þ2 3ðx xm1 Þ3 ; wm ðxÞ ¼ 3 h þ 3h ðxmþ1 xÞ þ 3hðxmþ1 xÞ 3ðxmþ1 xÞ ; h > > 3 > > : ðxmþ2 xÞ ; 0;
½xm2 ; xm1 ; ½xm1 ; xm ; ½xm ; xmþ1 ; ½xmþ1 ; xmþ2 ; otherwise:
ð21Þ Let UN ðx; tÞ be an approximation to the exact solution U ðx; tÞ. Then the global approximation UN ðx; tÞ can be written in terms of the splines as UN ðx; tÞ ¼
N þ1 X
wj ðxÞdj ðtÞ;
ð22Þ
j¼1
where dj are time dependent parameters to be determined from the boundary and weighted residual conditions. In terms of a local coordinate transformation n ¼ x xm ð0 6 n 6 hÞ, Eq. (21) can be expressed as 8 3 ðh nÞ > > wm1 < 3 2 1 3hðh nÞ2 3ðh nÞ3 ; 0 6 n 6 h: wm ¼ 3 h3 þ 3h2 ðh nÞ þ 2 3 h > > wmþ1 : h3 þ 3h n þ 3hn 3n n ð23Þ Since all other splines are 0 over the element ½xm ; xmþ1 , the approximation (22) over this element can be written in terms of the basis functions (23) as UN ¼
mþ2 X
w j dj :
ð24Þ
j¼m1
Using splines (23) and the approximation (24), the nodal values Um and Um0 at the node xm in terms of element parameters dm can be obtained as Um ¼ U ðxm Þ ¼ dm1 þ 4dm þ dmþ1 ; hUm0 ¼ hU 0 ðxm Þ ¼ 3ðdmþ1 dm1 Þ:
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
905
In the above equation and throughout this paper, the prime denotes derivative with respect to x. In the Galerkin method, basis functions (23) are used as weighting functions Wm . So, the weak form of the model problem given by Eq. (11) is Z
1
Wm ðUt Uxx arÞ dx ¼ 0:
ð25Þ
0
As mentioned in Section 2, it should be sufficient to work on the electrical conductivities given by Eqs. (9) and (10) only. Using the electrical conductivities (9) and (10) in Eq. (25), we respectively obtain the following integral equations over each finite element ½xm ; xmþ1 : Z xmþ1 ðWUt þ Wx Ux þ að1 dÞs1 WU að1 þ ð1 dÞs1 ÞW Þ dx xm
¼ bW ðxN ÞU ðxN ; tÞ
ð26Þ
and Z
xmþ1
ðWUt þ Wx Ux þ að1 dÞð1 s2 ÞWU aðð2 dÞ
xm
2ð1 dÞWs2 ÞÞ dx ¼ bW ðxN ÞU ðxN ; tÞ:
ð27Þ
Inserting Eqs. (23) and (24) into Eqs. (26) and (27) lead to mþ2 Z X j¼m1
xmþ1
mþ2 Z X e _ wi wj dx dj þ
xm mþ2 Z X
¼ að1 þ ð1 dÞs1 Þ
dej
mþ2 X wi wj dx dej þ b wi ðxN Þwj ðxN Þdej
xm
j¼m1
Z
xmþ1
w0i w0j dx
xm
j¼m1
þ að1 dÞs1
xmþ1
j¼m1
xmþ1
wi dx xm
and mþ2 Z X j¼m1
xmþ1
mþ2 Z X e _ wi wj dx dj þ
xm
þ að1 dÞð1 s2 Þ
j¼m1
¼ aðð2 dÞ 2ð1 dÞs2 Þ
Z
xmþ1
xm xmþ1
xm
w0i w0j dx
dej
xm
j¼m1 mþ2 Z X
xmþ1
mþ2 X wi wj dx dej þ b wi ðxN Þwj ðxN Þdej
wi dx ;
j¼m1
906
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
which are also be written in matrix forms as Ae d_ e þ ðc1 Ae þ Be þ C e Þde ¼ ða þ c1 Þd e
ð28Þ
and Ae d_ e þ ðc2 Ae þ Be þ C e Þde ¼ ðað1 ð1 dÞs2 Þ þ c2 Þd e ;
ð29Þ
respectively, where c1 ¼ að1 dÞs1 and c2 ¼ að1 dÞð1 s2 Þ. In the above equations and throughout this paper, the dot denotes derivative with respect to t, de ¼ ðdm1 ; dm ; dmþ1 ; dmþ2 Þ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 wj dx; xm Z xmþ1 w0i w0j dx; Beij ¼ ð30Þ xm Cije ¼ bwi ðxN Þwj ðxN Þ; Z xmþ1 die ¼ wi dx; xm
where i; j ¼ m 1; m; m þ 1; m þ 2. After a simple arrangement, the element matrices (30) can be obtained as 2 3 20 129 60 1 60 7 h 6 6 129 1188 933 7 Aeij ¼ 6 7; 140 4 60 933 1188 129 5 1 18 1 6 6 21 Beij ¼ 6 10h 4 36
60 21 102
129 36 87
87
102
3
36
21
2
20 3 3 36 7 7 7; 21 5
2
3 1 7 h6 6 11 7 die ¼ 6 7 4 4 11 5
18
1
and Cije is a zero matrix but 2 3 0 0 0 0 6 0 b 4b b 7 7 CijN 1 ¼ 6 4 0 4b 16b 4b 5: 0 b 4b b Assembling all contributions from all elements, Eqs. (28) and (29) lead to matrix differential equations Ad_ þ ðc1 A þ B þ CÞd ¼ ða þ c1 Þd ð31Þ and Ad_ þ ðc2 A þ B þ CÞd ¼ ðað1 ð1 dÞs2 Þ þ c2 Þd;
ð32Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
907
T
respectively, where d ¼ ðd1 ; d0 ; . . . ; dN þ1 Þ and A, B and C are ðN þ 3Þ ðN þ 3Þ matrices and d is an ðN þ 3Þ 1 matrix which are defined from the element matrices Ae , Be , C e and d e , respectively. Substituting the Crank–Nicolson approach d ¼ ðdn þ dnþ1 Þ=2 and the forward finite difference d_ ¼ ðdnþ1 dn Þ=k into Eqs. (31) and (32) yields the system of equations ½ð2 þ c1 kÞA þ kðB þ CÞ dnþ1 ¼ ½ð2 c1 kÞA kðB þ CÞ dn þ 2kða þ c1 Þd ð33Þ and ½ð2 þ c2 kÞA þ kðB þ CÞ dnþ1 ¼ ½ð2 c2 kÞA kðB þ CÞ dn þ 2kðað1 ð1 dÞs2 Þ þ c2 Þd;
ð34Þ
respectively, where kð DtÞ is the time step. Each resulting system is an ðN þ 3Þ ðN þ 3Þ septadiagonal equation system and so is easily and efficiently solved with a variant of the Thomas algorithm. A Fourier stability analysis shows the schemes (33) and (34) to be unconditionally stable since the Crank–Nicolson approach is used for the element parameters de .
5. Numerical results and conclusion To make a comparison with earlier work [1] we consistently use a choice of values for a (0.1, 0.5, 1, 102 and 105 ) with b ¼ 0:5 satisfying inequalities (12), (15), (16), (19) and (20) to obtain steady-state solutions for the cold, warm/ cold, warm, hot/warm and hot phases, respectively. We use the steady-state criteria 1 maxfjU0new U0old j; jUNnew UNold jg 6 0:5 108 k to get the Galerkin finite element solution of each particular phase. All calculations were performed on a Pentium 4 processor using Microsoft Fortran Compiler. In order to show how good the numerical solutions of the problem in comparison with the exact ones we shall use the error norm L2 defined by " #1=2 N 1 X 2 L2 ¼ jU ðxm Þ U ðxm Þj : N þ 1 m¼0 Cold steady-state solutions: The values a ¼ 0:1 and b ¼ 0:5 satisfy inequality (12) and so give a cold steady-state. Replacing s1 by 0, Eq. (33) is solved using the Thomas algorithm. The numerical results obtained by the GFE method
908
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
Table 1 Steady-state solutions in the cold phase x
GFE, [1] and exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.2500 0.2495 0.2480 0.2455 0.2420 0.2375 0.2320 0.2255 0.2180 0.2095 0.2000 t1 (GFE) ¼ 39.56, t1 ([1]) ¼ 39.54
(for k ¼ 0:01, h ¼ 0:05) are compared with exact one and also those obtained by the EFD method in [1] as shown in Table 1. It is seen that they are all in very good agreement with each other. Warm/cold steady-state solutions: The values a ¼ 0:5 and b ¼ 0:5 satisfy inequality (15) and then give a warm/cold steady-state. This phase has a free boundary s1 ðtÞ 2 ð0; 1 (see Fig. 2(b)) which is an unknown function of time and must be determined as a part of solution. If the predicted temperature is U ðxm Þ > 1 and U ðxmþ1 Þ < 1 at a time level n, the boundary s1 is between adjacent values U ðxm Þ and U ðxmþ1 Þ. Since U ðs1 ; tÞ ¼ 1, the location of the free boundary s1 can be estimated by using a linear interpolation as [1] sn1 ¼ xm þ
hð1 U ðxm ÞÞ : U ðxmþ1 Þ U ðxm Þ
ð35Þ
Thus, the motion of sðtÞ can be obtained from s_ 1 ¼
s1nþ1 sn1 : Dt
Eq. (33) is again solved using the Thomas algorithm. The results obtained by the GFE method (for k ¼ 0:01, h ¼ 0:1, 0.05, 0.025, 0.0125) are presented in Table 2. It is observed that the numerical solutions are in good agreement with the exact one, and also exhibit the expected convergence as the mesh size h is refined. The results obtained by the GFE method (for k ¼ 0:01, h ¼ 0:05) are also compared with those obtained by the EFD method [1]. It is seen that they agree with each other, as shown in Table 3. The steady-state value of the interface s1 (Eq. (35)) obtained by the GFE method using Eq. (33) is about 0.84 in some 22.25 units of steady-state time
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
909
Table 2 Steady-state solutions in the warm/cold phase
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:0125
Exact solution
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.076959 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
1.156163 1.153989 1.147458 1.136544 1.121199 1.101360 1.076945 1.047849 1.013953 0.975113 0.931167
1.156154 1.153980 1.147450 1.136535 1.121191 1.101353 1.076937 1.047843 1.013946 0.975107 0.931162
L2 105
0.231657
0.047376
0.018918
0.008275
x
Numerical solution
Table 3 Comparison of steady-state solutions in the warm/cold phase with results from [1] at h ¼ 0:05 and k ¼ 0:01 x
GFE
[1] EFD
Exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.156179 1.154005 1.147474 1.136559 1.121214 1.101375 1.076959 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
t1 ¼ 22:25
t1 ¼ 23:82
ðt1 Þ whereas it is obtained by the EFD method in [1] as about 0.81 in something 23.82 units of steady-state time. Table 4 exhibits the approximate and exact values of s1 and s_ 1 with the steady-state times at various mesh sizes. It is clearly seen that as the mesh size h is reduced the predicted steady-state values of s1 and s_ 1 smoothly and rapidly reach their exact steady-state values. Warm steady-state solutions: A warm steady-state solution can be easily obtained by using the values a ¼ 1 and b ¼ 0:5. Replacing s1 by 1, Eq. (33) is again solved using the Thomas algorithm. The results obtained by the GFE method (for k ¼ 0:01, h ¼ 0:1, 0.05) are displayed in Table 5. There is a good
910
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
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 times ðt1 Þ
s1 s_ 1 107 t1
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:0125
Exact
0.836076 0.596046 22.23
0.837156 0.596024 22.25
0.837344 0.596006 22.27
0.837435 0.596002 32.23
0.837581 0.0
agreement between numerical and exact solutions. It is seen that GFE solutions exhibit the expected convergence as the mesh size h is reduced. As shown in Table 5, the numerical solutions (for h ¼ 0:05) hit the exact steady-state ones, and also agree with those given in [1]. Thus, it is not necessary to take smaller values of h. Hot/warm steady-state solutions: It is sufficient to use the values a ¼ 102 and b ¼ 0:5 to get a hot/warm steady-state solution. As shown in Fig. 2(d), this phase has a free boundary s2 2 ð0; 1 which is an unknown function of time and so must be determined as a part of solution with a similar manner to s1 . If the predicted temperature is U ðxm Þ > 2 and U ðxmþ1 Þ < 2 at a time level n, the boundary s2 is between adjacent values U ðxm Þ and U ðxmþ1 Þ. Since U ðs2 ; tÞ ¼ 2, the location of the free boundary s2 can be estimated by a linear interpolation as [1]
sn2 ¼ xm þ
hð2 U ðxm ÞÞ : U ðxmþ1 Þ U ðxm Þ
ð36Þ
Table 5 Comparison of steady-state solutions in the warm phase with results from [1] x 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
GFE
EFD [1]
h ¼ 0:1
h ¼ 0:05
h ¼ 0:05
1.486325 1.483754 1.476017 1.463035 1.444679 1.420766 1.391055 1.355249 1.312991 1.263857 1.207354
1.486324 1.483754 1.476016 1.463035 1.444679 1.420765 1.391054 1.355249 1.312991 1.263856 1.207354
1.486363 1.483793 1.476058 1.463080 1.444729 1.420822 1.391120 1.355325 1.313080 1.263961 1.207476
t1 ¼ 13:94
t1 ¼ 14:35
Exact solution 1.486324 1.483754 1.476016 1.463035 1.444679 1.420765 1.391054 1.355249 1.312991 1.263856 1.207354
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
911
Thus, the motion of s2 ðtÞ can be easily obtained by using the usual forward difference formulae as s_ 2 ¼
s2nþ1 sn2 : Dt
Eq. (34) is solved using a variant of Thomas algorithm. The obtained GFE solutions (for k ¼ 0:01, h ¼ 0:1, 0.05, 0.025, 0.0125) are presented in Table 6. They are in good agreement with each other and also exact one. It is seen that numerical approximations show a high degree of accuracy when the mesh size h is reduced. Again, good agreement with the exact values is evident, as in convergence. Table 7 displays a comparison of approximate steady-state solutions obtained by the GFE method (for k ¼ 0:01, h ¼ 0:05) with those obtained by the EFD method in [1] which are reasonably in good agreement with each other. The steady-state values of the free boundary s2 (Eq. (36)) obtained by the GFE method using Eq. (34) is about 0.03 in something 6.94 units of steadystate time whereas it is obtained by using EFD method in [1] as about 0.12 in some 0.69 units of steady-state time. Table 8 displays the approximate and exact values of s2 and s_ 2 with the steady-state times at various mesh sizes. It is observed that as the mesh size h is reduced the predicted values of s2 and s_ 2 smoothly and rapidly approach to their exact values. Hot steady-state solutions: The values a ¼ 105 and b ¼ 0:5 is sufficient to obtain a hot steady-state solution. Replacing s2 by 1, Eq. (34) is solved using the Thomas algorithm. The GFE solutions (for k ¼ 0:01, h ¼ 0:05) are compared with EFD solutions (for h ¼ 0:05) given in [1] and also the exact solution, in Table 9. It is seen that they are all in full agreement with each other. Table 6 Steady-state solutions in the hot/warm phase ðe ¼ 3:0 107 Þ x
Numerical solution
Exact solution
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.999975 1.999918 1.999761 1.999338 1.998197 1.995120 1.986815 1.964410 1.904081
2.000001 1.999995 1.999974 1.999916 1.999756 1.999327 1.998172 1.995066 1.986712 1.964241 1.903806
2þe 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903602
L2 105
5.615767
1.710306
0.004061
2þe 1.999995 1.999974 1.999914 1.999752 1.999318 1.998152 1.995023 1.986626 1.964088 1.903601
912
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
Table 7 Comparison of steady-state solutions in the hot/warm phase with results from [1] at h ¼ 0:05 and k ¼ 0:01 x
GFE
[1] EFD
Exact
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.999995 1.999974 1.999916 1.999756 1.999327 1.998172 1.995066 1.986712 1.964241 1.903806
2.000008 2.000003 1.999985 1.999929 1.999976 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
t1 ¼ 6:94
t1 ¼ 0:69
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 times ðt1 Þ
s2 s_ 2 107 t1
h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
Exact
0.015077 0.785326 3.36
0.020950 0.091881 6.94
0.025315 0.072008 23.26
0.025334 0.0
Table 9 Steady-state solutions in the hot phase x
GFE, [1] and exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
2.500 2.495 2.480 2.455 2.420 2.375 2.320 2.255 2.180 2.095 2.000 t1 (GFE) ¼ 44.96, t1 ([1]) ¼ 54.18
S. Kutluay, A. Esen / Appl. Math. Comput. 161 (2005) 897–913
913
In this paper, the approximate steady-state solutions of the PTC thermistor problem with temperature dependent electrical conductivities described in Section 2 using the cubic GFE method have been presented. It is shown that the numerical results obtained by the GFE method exhibit the correct physical phenomena of the PTC thermistor problem. As a conclusion, we have demonstrated that the proposed electrical conductivities can be used effectively to obtain GFE solutions in transition phases (warm/cold and hot/warm) of the thermistor. References [1] S. Kutluay, A.S. Wood, Numerical solutions of the thermistor problem with a ramp electrical conductivity, Appl. Math. Comput. 148 (2004) 145–162. [2] R.P. Turner, ABCÕs of Thermistors, Foulsham-Sams, Slough, 1970. [3] F.J. Hyde, Thermistors, Butterworth (ILIFFE Books), London, 1971. [4] E.D. Macklen, Thermistors, Electrochemical Publications, Scotland, 1979. [5] G. Cimatti, Remark on existence and uniqueness for the thermistor problem under mixed boundary conditions, Q. Appl. Math. 47 (1989) 117–121. [6] 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. [7] A. Fowler, I. Frigaard, S.D. Howison, Temperature surges in thermistors, SIAM J. Appl. Math. 52 (1992) 998–1011. [8] S.D. Howison, J.F. Rodrigues, M. Shillor, Stationary solutions to the thermistor problem, J. Math. Anal. Appl. 174 (1993) 573–588. [9] J. Wiedmann, The thermistor problem, Nonlin. Differ. Equ. Appl. 4 (1997) 133–148. [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.), Proc. 5th Eur. Conf. Math. Industry, Kluwer Academic, Stuttgart, 1991, pp. 397–400. [12] S. Kutluay, A.R. Bahadir, A. Ozdes, Various methods to the thermistor problem with a bulk electrical conductivity, Int. J. Numer. Meth. Engng. 45 (1999) 1–12. [13] S. Kutluay, A. Esen, A B-spline finite element method for the thermistor with the modified electrical conductivity, Appl. Math. Comput. (in press). [14] D.R. Westbrook, The thermistor: a problem in heat and current flow, Numer. Meth. PDEs 5 (1989) 259–273. [15] 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. [16] A.S. Wood, S. Kutluay, A heat balance integral model of the thermistor, Int. J. Heat Mass Transfer 38 (10) (1995) 1831–1840. [17] P.M. Prenter, Splines and Variational Methods, John-Wiley, New York, 1975.