- Email: [email protected]
Numerical solutions of the thermistor problem by spline finite elements
Applied Mathematics and Computation 162 (2005) 475–489 www.elsevier.com/locate/amc
Numerical solutions of the thermistor problem by spline finite elements 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 (PTC) thermistor problem, having a step function electrical conductivity that is a highly non-linear function of the temperature, using subdomain collocation and Petrov–Galerkin methods based on spline finite elements. The resulting system of ordinary differential equations is solved by the usual Crank–Nicolson finite difference method using a variant of Thomas algorithm. It is shown that the numerical solutions obtained by the present methods exhibit the correct physical characteristics of the problem and, they are in very good agreement with the exact solution. A Fourier stability analysis of each method is also investigated. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Thermistor; Step electrical conductivity; Subdomain collocation; Petrov–Galerkin
1. Introduction Thermistor is a thermo-electric device made from a ceramic material whose electrical conductivity varies strongly with temperature which can be used as a switching device in many electronic circuits. The study of thermistor problems in heat and current flow has a long history of applications in several areas of electronics and its related industries [1,2]. There are generally two kinds of thermistors; one is a positive temperature coefficient (PTC) thermistor which its electrical conductivity decreases with
*
Corresponding author. E-mail address: [email protected] (S. Kutluay).
0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.125
476
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
increasing temperature, and the other is a negative temperature coefficient thermistor which its electrical conductivity increases with increasing temperature. These problems are highly nonlinear since each problem consists of a heat equation with temperature dependent electrical conductivity. The steady-state thermistor problem has been studied by several authors [3– 9] where existence and uniqueness of solutions were given. Westbrook [10] obtained steady-state solutions of a one-dimensional thermistor problem with a step function electrical conductivity by using an iterative method based on finite elements. Preis et al. [11] investigated a two-dimensional thermistor problem by the finite element method. Wood and Kutluay [12] presented a simplified model of one particular thermistor and gave an approximate functional solution for the evolving thermistor problem by an application of the heat balance integral. Kutluay et al. [13] obtained approximate steady-state solutions of the problem under various assumptions on the step function electrical conductivity using explicit finite difference, linear Galerkin finite element and heat balance integral methods. More recently, Kutluay and Wood [14] gave approximate steady-state solutions of a PTC thermistor problem, having a ramp electrical conductivity that is a higly nonlinear function of the temperature, using a standard explicit finite difference method. In this study, we have applied subdomain collocation and Petrov–Galerkin methods based on spline finite elements to the PTC thermistor problem with a temperature dependent step function electrical conductivity to obtain its approximate steady-state solutions.
2. Statement of the problem This paper deals with a one-dimensional PTC thermistor problem which can be modelled as follows: The temperature U ðx; tÞ is governed by the heat-flow problem as 2 oU o2 U o/ ¼ 2 þ ar ; 0 < x < 1; t > 0; ð1Þ ot ox ox subject to the boundary conditions oU ¼ 0; ox
x ¼ 0; t > 0;
oU þ bU ¼ 0; ox
x ¼ 1; t > 0
ð2Þ ð3Þ
and the initial condition U ðx; 0Þ ¼ 0;
0 6 x 6 1;
ð4Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
477
in which a is the ratio of electric heating to heat diffusion and b is a positive heat transfer coefficient. The electrical potential /ðx; tÞ is governed by the current-conservation problem as o o/ r ¼ 0; 0 < x < 1; t > 0; ð5Þ ox ox subject to the boundary conditions /ð0; tÞ ¼ 0;
t > 0;
/ð1; tÞ ¼ 1;
t>0
and the initial condition /ðx; 0Þ ¼ x;
0 6 x 6 1;
in which r is the electrical conductivity. In many cases of practical interest, r is a function of the temperature which has been taken by several authors [10,12,15] as 1; 0 6 U 6 1; rðU Þ ¼ ð6Þ d; U > 1; which is mathematically equivalent to a step function as illustrated in Fig. 1 and where, typically, d ¼ 105 . Since U ðx; 0Þ ¼ 0 and b > 0, we assume monotonicity of the temperature profile such that the point x ¼ 0 will always be the hottest and will also be the first point to reach the critical temperature Uc ¼ 1 above which the conductivity drops. As a result of the decreasing electrical conductivity the rate of heat loss at x ¼ 1 will be equal to the internal heat generation and an equilibrium will be attained. Westbrook [10] showed that three distinct phases are observable, as the cold, warm and hot phases. Hence, by the conductivity (6), the σ 1
δ 0
1
U
Fig. 1. Behaviour of step electrical conductivity r against U.
478
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489 U
U
U
Uc =1
Uc =1
Uc =1
0
x
1
0
(a) Cold
s
1
x
(b) Warm
0
1
x
(c) Hot
Fig. 2. Steady-state configurations.
steady-state temperature may be any one of the three scenarios shown in Fig. 2 [13]. As seen in Fig. 2(b) the warm phase has two distinct regions: a hot region 0 6 x 6 s (above Uc ¼ 1) of low electrical conductivity separated by a free boundary s from a cold region s 6 x 6 1 (below or equal to Uc ¼ 1) of high electrical conductivity. By the electrical conductivity defined by Eq. (6) the evolving problem given by Eqs. (1)–(4) will move sequentially through the cold, warm and hot phases (if at all). Kutluay et al. [16] showed that, in the warm phase, the sharp decreasing of the electrical conductivity from 1 to d with increasing temperature can cause oscillation in the predicted temperature when the finite difference methods are applied to the thermistor problem. It is not able to be recovered during the process since the temperature domain in the cold and warm regions of the warm phase are not symmetrical. As a result of this, they studied with the modified electrical conductivity defined by rðU Þ ¼ 1 þ ðd 1Þs:
ð7Þ
It is clearly seen that in the cold phase, s ¼ 0 and so r ¼ 1 which is the same with the conductivity (6), in the hot phase, s ¼ 1 and so r ¼ d which is again the same with the conductivity (6), and but in the warm phase, 0 < s < 1 and so r ¼ 1 þ ðd 1Þs: An oscillation also occurs when the finite element methods apply to the problem. So, in this study, we will also consider the electrical conductivity given by Eq. (7). By virtue of the conductivity (7) the exact solution of the current-conservation problem given by Eq. (5) can be easily obtained as /ðx; tÞ ¼ x;
0 6 x 6 1; t P 0
and so the thermistor problem is reduced to the heat conduction problem oU o2 U ¼ 2 þ ar; ot ox
0 < x < 1; t > 0;
ð8Þ
subject to the conditions (2)–(4). In this paper, we solve this heat conduction problem for each particular steady-state phase by using subdomain collocation and Petrov–Galerkin methods.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
479
3. Exact steady-state solutions of the problem An exact steady-state solution of the problem given by Eq. (8) with the conditions (2)–(4) can be easily obtained by the relative values of a, b and d for each particular phase. In this section, the warm phase solution will be particularly examined under the electrical conductivity given by Eq. (7) since the solutions of the cold and hot phases are the same as in [12]. The asterisk * denotes exact steady-state solutions in this study. 3.1. Cold phase (s ¼ 0, r ¼ 1) The steady-state solution to this phase is a 2 1 þ x2 ; 0 6 x 6 1: U ðxÞ ¼ 2 b
ð9Þ
In this phase, 0 6 U ð1Þ 6 U ð0Þ 6 1 provides the inequality a 6 2b=ð2 þ bÞ
ð10Þ
giving a cold steady-state solution [10,12]. 3.2. Warm phase (0 < s < 1, r ¼ 1 þ (d 1)s ) The warm exact steady-state solution can be easily obtained as a½1 þ ðd 1Þs 2 1 þ x2 ; 0 6 x 6 1: U ðxÞ ¼ 2 b
ð11Þ
Since U ðs Þ ¼ 1, Eq. (11) leads to the following equation for the steadystate values s of sðtÞ: 2 2 a½1 þ ðd 1Þs 1 þ s 2 ¼ 0; ð12Þ b which is solved numerically by using the Newton method and used in Eq. (11) to determine the exact steady-state temperature. Since s 2 ð0; 1, for all time, from Eq. (12) we obtain the inequality að2 þ bÞ > 2b P 2ad
ð13Þ
implies a warm steady-state which is identical to that obtained by Westbrook [10] and Wood and Kutluay [12].
480
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
3.3. Hot phase (s ¼ 1, r ¼ d) The steady-state solution to this phase is ad 2 2 1þ x : U ðxÞ ¼ 2 b
ð14Þ
Since U ð1Þ P 1 in this phase, the values of a, b and d satisfying the inequality ð15Þ
b 6 ad; giving a hot steady-state [10,12].
4. Methods of solution The interval 0 ¼ x0 < x1 < < xN ¼ 1 is divided into N equal finite elements by h ¼ 1=N . The quadratic B-splines Qm ðxÞ, at the knots xm which form a basis for functions over ½0; 1 are defined as [17] 8 2 2 2 ðxmþ2 xÞ 3ðxmþ1 xÞ þ 3ðxm xÞ ; ½xm1 ; xm ; > > < 2 2 1 ½xm ; xmþ1 ; Qm ðxÞ ¼ 2 ðxmþ2 xÞ2 3ðxmþ1 xÞ ; h > ðx xÞ ; ½x > mþ1 ; xmþ2 ; : mþ2 0; otherwise: ð16Þ So the global approximation UN ðx; tÞ to the function U ðx; tÞ can be written in terms of the splines as UN ðx; tÞ ¼
N X
Qj ðxÞdj ðtÞ;
ð17Þ
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 ;
ð18Þ
0 6 n 6 h:
Eq. (16) can be expressed as 8 2 Qm1 1 < ðh2 nÞ ; Qm ¼ 2 h þ 2hn 2n2 ; h : 2 Qmþ1 n;
0 6 n 6 h:
ð19Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
481
Since all other splines are zero over the element ½xm ; xmþ1 , the approximation (17) over this element can be written in terms of the basis functions given by Eq. (19) as UN ðx; tÞ ¼
mþ1 X
Qj ðxÞdj ðtÞ:
ð20Þ
j¼m1
Using splines (19) and the approximation (20), the nodal values Um and Um0 at the node xm are given in terms of element parameters dm by Um ¼ U ðxm Þ ¼ dm1 þ dm ;
ð21Þ
2 Um0 ¼ U 0 ðxm Þ ¼ ðdm dm1 Þ; h
ð22Þ
where the prime denotes differentiation with respect to x. Applying a weighted residual method to Eq. (8) with weighting functions Wm produces the weak form Z xN Wm ðUt Uxx arÞ dx ¼ 0; ð23Þ x0
where m ¼ 0; 1; . . . ; N 1. 4.1. Subdomain collocation method In the subdomain collocation (SC) method [18], the weighting functions Wm are 1; xm 6 x 6 xmþ1 ; Wm ¼ 0; otherwise and so, the weak form (23) over each finite element ½xm ; xmþ1 yields the integral equation Z xmþ1 ðUt Uxx arÞ dx ¼ 0: ð24Þ xm
Integrating (24) and using Eqs. (20) and (22), we obtain h _ 2 ðdm1 þ 4d_ m þ d_ mþ1 Þ ðdm1 2dm þ dmþ1 Þ har ¼ 0; 3 h
ð25Þ
where the dot denotes differentiation with respect to t. Substituting the following forward difference for d_ and the weight average for d d_ ¼ ðdnþ1 dn Þ=Dt;
d ¼ ð1 hÞdn þ hdnþ1 ;
in Eq. (25) results in a system of N 1 linear algebraic equations as
ð26Þ
482
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489 nþ1 ð1 6hrÞdm1 þ ð4 þ 12hrÞdmnþ1 þ ð1 6hrÞdnþ1 mþ1
¼ ð1 þ 6ð1 hÞrÞdnm1 þ ð4 12ð1 hÞrÞdnm þ ð1 þ 6ð1 hÞrÞdnmþ1 þ 3kar;
ð27Þ
where (m ¼ 1; 2; . . . ; N 2), Dtð kÞ is the time step, r ¼ k=h2 and 0 6 h 6 1 is the weight factor. The cases h ¼ 0, h ¼ 1=2 and h ¼ 1 correspond, respectively, to the explicit, Crank–Nicolson and fully implicit finite difference methods. Using Eqs. (21) and (22), the boundary conditions (2) and (3) become dn0 ¼ dn1 ;
ð28Þ
nþ1 d0nþ1 ¼ d1
and dnN ¼ cdnN 1 ;
ð29Þ
dNnþ1 ¼ cdNnþ1 1 ;
where c ¼ ð2 hbÞ=ð2 þ hbÞ. Substituting the expressions (28) and (29) in Eq. (27) and using a Crank– Nicolson approach we obtain the resulting system of equations ð5 þ 3rÞdnþ1 þ ð1 3rÞd1nþ1 ¼ ð5 3rÞdn0 þ ð1 þ 3rÞdn1 þ 3kar; 0 3rÞdnþ1 m1
6rÞdnþ1 m
ð1 þ ð4 þ þ ð1 ¼ ð1 þ 3rÞdnm1 þ ð4 6rÞdnm þ ð1 þ 3rÞdnmþ1 þ 3kar; 3rÞdnþ1 N2
m ¼ 0;
3rÞdnþ1 mþ1
m ¼ 1; 2; . . . ; N 2;
cÞÞdnþ1 N1
ð1 þ ð4 þ c þ 3rð2 ¼ ð1 þ 3rÞdnN 2 þ ð4 þ c 3rð2 cÞÞdnN1 þ 3kar;
m ¼ N 1:
ð30Þ
This system can be easily and efficiently solved by using a variant of Thomas algorithm. 4.1.1. Stability analysis of the method For stability analysis it is convenient to use the Fourier series method pffiffiffiffiffiffiffi (see, e.g., [19], Section 2). Substituting the Fourier mode dnm ¼ gn eiumh (i ¼ 1) into Eq. (30) leads to the growth factor g of the form g¼
2ð1 þ 3rÞ cos2 ðuh=2Þ þ 1 6r : 2ð1 3rÞ cos2 ðuh=2Þ þ 1 þ 6r
For the stability, g must be jgj 6 1 which is satisfied for all positive values of stability parameter r only. So, the system (30) is unconditionally stable.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
483
4.2. Petrov–Galerkin method Applying a Petrov–Galerkin (PG) method with the weighting functions Wm defined by linear B-spline basis functions [17] 8 ðx xÞ 2ðxm xÞ; ½xm1 ; xm ; 1 < mþ1 Wm ðxÞ ¼ xmþ1 x; ½xm ; xmþ1 ; ð31Þ h: 0; otherwise; the weak form (23), over each finite element ½xm ; xmþ1 , with some manipulation, leads to the integral equation Z xmþ1 ðWUt þ Wx Ux W arÞ dx ¼ bW ðxN ÞU ðxN ; tÞ: ð32Þ xm
Using the local coordinate transformation (18), from Eq. (31) we obtain the linear B-spline functions for the finite element ½xm ; xmþ1 Wm ¼ 1 n=h;
Wmþ1 ¼ n=h:
Substituting Eq. (20) in Eq. (32) leads to mþ1 Z xmþ1 mþ1 Z xmþ1 X X Wi Qj dx d_ ej þ Wi 0 Q0j dx dej j¼m1
xm
xm
j¼m1
Z mþ1 X e þb Wi ðxN ÞQj ðxN Þdj ¼ ar j¼m1
xmþ1
Wi dx ;
xm
which also be written in a matrix form as Ae d_ e þ ðBe þ C e Þde ¼ d e ;
ð33Þ
where 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 Wi Qj dx; Aeij ¼ xm Z xmþ1 Wi 0 Q0j dx; Beij ¼ ð34Þ xm e Cij ¼ bWi ðxN ÞQj ðxN Þ; Z xmþ1 die ¼ ar Wi dx; xm
where i and j takes only the values, respectively m, m þ 1 and m 1, m, m þ 1 for the element ½xm ; xmþ1 . The element matrices (34) can be obtained, with some manipulation, as
484
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
Aeij
h 3 ¼ 12 1
8 8
1 ; 3
Beij
1 1 ¼ h 1
0 0
1 ; 1
die
ar 1 ¼ 2 1
and Cije is a zero matrix but 0 0 0 N 1 : Cij ¼ 0 b b Assembling all contributions from all elements, Eq. (33) leads to the matrix differential equation Ad_ þ ðB þ CÞd ¼ d;
ð35Þ T
where d ¼ ðd1 ; d0 ; . . . ; dN Þ and A, B and C are ðN þ 1Þ ðN þ 2Þ matrices and d is an ðN þ 1Þ 1 matrix which are defined from the element matrices Ae , Be , C e and d e , respectively. Using a Crank–Nicolson approach for d and forward finite difference for d_ given by Eq. (26), Eq. (35) takes the form Dt Dt A þ ðB þ CÞ dnþ1 ¼ A ðB þ CÞ dn þ dDt: ð36Þ 2 2 Using one of the boundary conditions (2) and (3) in the system (36) leads to an ðN þ 1Þ ðN þ 1Þ matrix system which can be easily solved by the Thomas algorithm. A typical member of Eq. (36) in terms of the nodal parameters dnm is nþ1 nþ1 nþ1 þ a2 dm1 þ a2 dnþ1 a1 dm2 m þ a1 dmþ1
¼ a3 dnm2 þ a4 dnm1 þ a4 dnm þ a3 dnmþ1 ;
ð37Þ
where a1 ¼
h k ; 12 2h
a2 ¼
11h k þ ; 12 2h
a3 ¼
h k þ ; 12 2h
a4 ¼
11h k : 12 2h
4.2.1. Stability analysis of the method For stability analysis, again substituting the Fourier mode dnm ¼ gn eiumh , pffiffiffiffiffiffi ffi (i ¼ 1) in Eq. (37) leads to the growth factor g of the form g¼
a ib ; c id
where a ¼ a3 cos 2uh þ h cos uh þ a4 ;
b ¼ a3 sin 2uh ða3 a4 Þ sin uh;
c ¼ a1 cos 2uh þ h cos uh þ a2 ;
d ¼ a1 sin 2uh ða1 a2 Þ sin uh:
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
485
For the stability, the magnitude of the growth factor g must be jgj 6 1 which is mathematically equivalent to c2 þ d 2 a2 b2 P 0. By virtue of any algebraic package programme it is seen that c2 þ d 2 a2 b2 P 0 for h > 0, k > 0 and u > 0. Thus, the scheme (37) is unconditionally stable.
5. Numerical results and conclusion In order to evaluate the temperature profiles we use the values of a, b and d satisfying inequalities (10), (13) and (15) for the cold, warm and hot phases, respectively. We also use the steady-state test
1 max jU0new U0old j; jUNnew UNold j 6 0:5 108 ; k
Nh ¼ 1
ð38Þ
for any one of the three phases. All calculations were performed in Fortran code under a Pentium 4 processor using double precision arithmetic. To test the accuracy of both numerical methods we have calculated the kek1 norm defined as N 1 X U ðxm Þ ; kek1 ¼ 1 N þ 1 m¼0 U ðxm Þ
e ¼ ½e0 ; e1 ; . . . ; eN T :
We insistly take the same values of a, b and d for each particular phase in case the reader easily compares the results obtained using the present methods with earlier works [12,13,16]. 5.1. Cold steady-state The values a ¼ 0:1 and b ¼ 0:2 satisfying inequality (10) are sufficient to obtain a cold steady-state solution. Such a choice of the values of a and b is only typical of a problem evolving to a cold steady-state. Alternatively, different values of a and b, still satisfying the inequality (10), simple serve to alter the time taken to attain steady-state. The results obtained by the present methods (for k ¼ 0:01 and h ¼ 0:1) were reported when the steady-state test (38) is satisfied, and presented in Table 1 with the exact solution (9). It is observed that the numerical solutions are seen to be in exact agreement with the exact one. The only remarkable difference is the time to attain the steady-state. The subdomain collocation solution is obtained after something like 89.9 units of CPU time whereas the Petrov– Galerkin solution is obtained after some 65.32 units of CPU time.
486
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
Table 1 Solutions in cold steady-state x
U SC, PG 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.5500 0.5495 0.5480 0.5455 0.5420 0.5375 0.5320 0.5255 0.5180 0.5095 0.5000
5.2. Warm steady-state If the predicted temperature at a node is such that at time level n, U ðxm Þ > 1 and U ðxmþ1 Þ < 1 in the warm phase, the moving boundary s lies between adjacent values U ðxm Þ and U ðxmþ1 Þ. Since U ðx; tÞ ¼ 1 at x ¼ sðtÞ, the location of the interface s at time level n can be easily estimated by using a linear interpolation as [16] sn ¼ xm þ
hð1 U ðxm ÞÞ : U ðxmþ1 Þ U ðxm Þ
ð39Þ
Thus, the motion of sðtÞ can be obtained from s_ ¼
snþ1 sn : Dt
For a warm steady-state solution, the values a ¼ 0:5 and b ¼ 0:2 satisfying inequality (13) are sufficient. The numerical results obtained by the both methods (for k ¼ 0:01 and h ¼ 0:1, 0.05, 0.025, 0.0125) are compared with the exact solution (11) in Table 2. It is seen that numerical solutions are in good agreement with the exact one, and exhibit the expected convergence as the mesh size h is refined. The predicted location of the interface s is obtained by Eq. (39) which is in good agreement with the exact one given by Eq. (12). Fig. 3 shows the movement (location and speed) of the interface s. In both numerical methods, the boundary s smoothly reaches its steady-state value of about 0.62 in some 5.2 units of CPU time.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
487
Table 2 Solutions in warm steady-state x
U
Exact
SC and PG h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:0125
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.036749 1.035806 1.032979 1.028266 1.021669 1.013186 1.002819 0.990566 0.976429 0.960406 0.942499
1.036641 1.035699 1.032871 1.028160 1.021563 1.013081 1.002715 0.990464 0.976328 0.960307 0.942401
1.036587 1.035645 1.032818 1.028106 1.021510 1.013028 1.002663 0.990412 0.976277 0.960257 0.942352
1.036585 1.035643 1.032816 1.028104 1.021507 1.013026 1.002660 0.990410 0.976275 0.960255 0.942350
kek1
0.000160
0.000056
0.000004
0.000001
1.036583 1.035641 1.032814 1.028102 1.021506 1.013024 1.002659 0.990408 0.976273 0.960253 0.942348
Fig. 3. The movement of the interface sðtÞ.
5.3. Hot steady-state Both numerical solutions to this phase are easily obtained for the values a ¼ 105 and b ¼ 0:2 satisfying inequality (15). The agreement between our numerical results (for k ¼ 0:01 and h ¼ 0:1) and the exact solution (14) match each other, shown in Table 3. Again, only difference is the time to attain steady-state. The subdomain collocation solution is obtained after something
488
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
Table 3 Solutions in hot steady-state x
U SC, PG and exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
5.500 5.495 5.480 5.455 5.420 5.375 5.320 5.255 5.180 5.095 5.000
like 102.19 units of CPU time whereas the Petrov–Galerkin solution is obtained after some 77.61 units of CPU time. In order to show how good the numerical predictions exhibit the correct physical behaviour of the problem, we only give graph in Fig. 4. The numerical solutions (for k ¼ 0:01 and h ¼ 0:025) obtained by the present methods in the cold, warm and hot phases during the steady-state are drawn on the same diagram, but curves cannot be distinguishable since they are very close to each other.
6.0 5.5 5.0
Temperature, U
4.5
Hot
4.0 3.5 3.0 2.5 2.0 1.5 1.0
Warm Cold
0.5 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Distance, x Fig. 4. Temperature profiles in the cold, warm and hot phases at steady-state.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
489
In this paper, we have presented approximate steady-state solutions of a one-dimensional PTC thermistor problem with a temperature dependent step function electrical conductivity using subdomain collocation and Petrov– Galerkin methods. It is shown that the numerical solutions obtained by the present methods satisfy all the physical characteristics of the problem. As a conclusion, both methods are capable of producing numerical solutions for the PTC thermistor problem. Further, the present methods can be applied to a two-dimensional analogue of the problem.
References [1] F.J. Hyde, Thermistors, Butterwoth Group (ILIFFE books), London, 1971. [2] E.D. Macklen, Thermistors, Electrochemical Publications, Scotland, 1979. [3] G. Cimatti, Remark on existence and uniqueness for the thermistor problem under mixed boundary conditions, Q. Appl. Math. 47 (1989) 117–121. [4] 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. € [5] H. Diesselhorst, Uber das probleme eines elektrisch erwarmten leiters, Annal. Phys. 1 (1900) 312–325. [6] A. Fowler, I. Frigaard, S.D. Howison, Temperature surges in thermistors, SIAM J. Appl. Math. 52 (1992) 998–1011. [7] S.D. Howison, J.F. Rodrigues, M. Shillor, Stationary solutions to the thermistor problem, J. Math. Anal. Appl. 174 (1993) 573–588. [8] J. Wiedmann, The thermistor problem, Nonlinear Differ. Equ. Appl. 4 (1997) 133–148. [9] S. Zhou, D.R. Westbrook, Numerical solutions of the thermistor equations, J. Comput. Appl. Math. 79 (1997) 101–118. [10] D.R. Westbrook, The thermistor: a problem in heat and current flow, Numer. Meth. PDEs 5 (1989) 259–273. [11] 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. [12] A.S. Wood, S. Kutluay, A heat balance integral model of the thermistor, Int. J. Heat Mass Transfer 38 (10) (1995) 1831–1840. [13] 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. [14] S. Kutluay, A.S. Wood, Numerical solutions of the thermistor problem with a ramp electrical conductivity, Appl. Math. Comput. 148 (2004) 145–162. [15] A.S. Wood, Modelling the thermistor, in: M. Heili€ o (Ed.), Proceedings of the 5th European Conference on Mathematics in Industry, Kluwer Academic, Stuttgart, 1991, pp. 397–400. [16] 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. [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.
Numerical solutions of the thermistor problem by spline finite elements 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 (PTC) thermistor problem, having a step function electrical conductivity that is a highly non-linear function of the temperature, using subdomain collocation and Petrov–Galerkin methods based on spline finite elements. The resulting system of ordinary differential equations is solved by the usual Crank–Nicolson finite difference method using a variant of Thomas algorithm. It is shown that the numerical solutions obtained by the present methods exhibit the correct physical characteristics of the problem and, they are in very good agreement with the exact solution. A Fourier stability analysis of each method is also investigated. Ó 2004 Elsevier Inc. All rights reserved. Keywords: Thermistor; Step electrical conductivity; Subdomain collocation; Petrov–Galerkin
1. Introduction Thermistor is a thermo-electric device made from a ceramic material whose electrical conductivity varies strongly with temperature which can be used as a switching device in many electronic circuits. The study of thermistor problems in heat and current flow has a long history of applications in several areas of electronics and its related industries [1,2]. There are generally two kinds of thermistors; one is a positive temperature coefficient (PTC) thermistor which its electrical conductivity decreases with
*
Corresponding author. E-mail address: [email protected] (S. Kutluay).
0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.125
476
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
increasing temperature, and the other is a negative temperature coefficient thermistor which its electrical conductivity increases with increasing temperature. These problems are highly nonlinear since each problem consists of a heat equation with temperature dependent electrical conductivity. The steady-state thermistor problem has been studied by several authors [3– 9] where existence and uniqueness of solutions were given. Westbrook [10] obtained steady-state solutions of a one-dimensional thermistor problem with a step function electrical conductivity by using an iterative method based on finite elements. Preis et al. [11] investigated a two-dimensional thermistor problem by the finite element method. Wood and Kutluay [12] presented a simplified model of one particular thermistor and gave an approximate functional solution for the evolving thermistor problem by an application of the heat balance integral. Kutluay et al. [13] obtained approximate steady-state solutions of the problem under various assumptions on the step function electrical conductivity using explicit finite difference, linear Galerkin finite element and heat balance integral methods. More recently, Kutluay and Wood [14] gave approximate steady-state solutions of a PTC thermistor problem, having a ramp electrical conductivity that is a higly nonlinear function of the temperature, using a standard explicit finite difference method. In this study, we have applied subdomain collocation and Petrov–Galerkin methods based on spline finite elements to the PTC thermistor problem with a temperature dependent step function electrical conductivity to obtain its approximate steady-state solutions.
2. Statement of the problem This paper deals with a one-dimensional PTC thermistor problem which can be modelled as follows: The temperature U ðx; tÞ is governed by the heat-flow problem as 2 oU o2 U o/ ¼ 2 þ ar ; 0 < x < 1; t > 0; ð1Þ ot ox ox subject to the boundary conditions oU ¼ 0; ox
x ¼ 0; t > 0;
oU þ bU ¼ 0; ox
x ¼ 1; t > 0
ð2Þ ð3Þ
and the initial condition U ðx; 0Þ ¼ 0;
0 6 x 6 1;
ð4Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
477
in which a is the ratio of electric heating to heat diffusion and b is a positive heat transfer coefficient. The electrical potential /ðx; tÞ is governed by the current-conservation problem as o o/ r ¼ 0; 0 < x < 1; t > 0; ð5Þ ox ox subject to the boundary conditions /ð0; tÞ ¼ 0;
t > 0;
/ð1; tÞ ¼ 1;
t>0
and the initial condition /ðx; 0Þ ¼ x;
0 6 x 6 1;
in which r is the electrical conductivity. In many cases of practical interest, r is a function of the temperature which has been taken by several authors [10,12,15] as 1; 0 6 U 6 1; rðU Þ ¼ ð6Þ d; U > 1; which is mathematically equivalent to a step function as illustrated in Fig. 1 and where, typically, d ¼ 105 . Since U ðx; 0Þ ¼ 0 and b > 0, we assume monotonicity of the temperature profile such that the point x ¼ 0 will always be the hottest and will also be the first point to reach the critical temperature Uc ¼ 1 above which the conductivity drops. As a result of the decreasing electrical conductivity the rate of heat loss at x ¼ 1 will be equal to the internal heat generation and an equilibrium will be attained. Westbrook [10] showed that three distinct phases are observable, as the cold, warm and hot phases. Hence, by the conductivity (6), the σ 1
δ 0
1
U
Fig. 1. Behaviour of step electrical conductivity r against U.
478
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489 U
U
U
Uc =1
Uc =1
Uc =1
0
x
1
0
(a) Cold
s
1
x
(b) Warm
0
1
x
(c) Hot
Fig. 2. Steady-state configurations.
steady-state temperature may be any one of the three scenarios shown in Fig. 2 [13]. As seen in Fig. 2(b) the warm phase has two distinct regions: a hot region 0 6 x 6 s (above Uc ¼ 1) of low electrical conductivity separated by a free boundary s from a cold region s 6 x 6 1 (below or equal to Uc ¼ 1) of high electrical conductivity. By the electrical conductivity defined by Eq. (6) the evolving problem given by Eqs. (1)–(4) will move sequentially through the cold, warm and hot phases (if at all). Kutluay et al. [16] showed that, in the warm phase, the sharp decreasing of the electrical conductivity from 1 to d with increasing temperature can cause oscillation in the predicted temperature when the finite difference methods are applied to the thermistor problem. It is not able to be recovered during the process since the temperature domain in the cold and warm regions of the warm phase are not symmetrical. As a result of this, they studied with the modified electrical conductivity defined by rðU Þ ¼ 1 þ ðd 1Þs:
ð7Þ
It is clearly seen that in the cold phase, s ¼ 0 and so r ¼ 1 which is the same with the conductivity (6), in the hot phase, s ¼ 1 and so r ¼ d which is again the same with the conductivity (6), and but in the warm phase, 0 < s < 1 and so r ¼ 1 þ ðd 1Þs: An oscillation also occurs when the finite element methods apply to the problem. So, in this study, we will also consider the electrical conductivity given by Eq. (7). By virtue of the conductivity (7) the exact solution of the current-conservation problem given by Eq. (5) can be easily obtained as /ðx; tÞ ¼ x;
0 6 x 6 1; t P 0
and so the thermistor problem is reduced to the heat conduction problem oU o2 U ¼ 2 þ ar; ot ox
0 < x < 1; t > 0;
ð8Þ
subject to the conditions (2)–(4). In this paper, we solve this heat conduction problem for each particular steady-state phase by using subdomain collocation and Petrov–Galerkin methods.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
479
3. Exact steady-state solutions of the problem An exact steady-state solution of the problem given by Eq. (8) with the conditions (2)–(4) can be easily obtained by the relative values of a, b and d for each particular phase. In this section, the warm phase solution will be particularly examined under the electrical conductivity given by Eq. (7) since the solutions of the cold and hot phases are the same as in [12]. The asterisk * denotes exact steady-state solutions in this study. 3.1. Cold phase (s ¼ 0, r ¼ 1) The steady-state solution to this phase is a 2 1 þ x2 ; 0 6 x 6 1: U ðxÞ ¼ 2 b
ð9Þ
In this phase, 0 6 U ð1Þ 6 U ð0Þ 6 1 provides the inequality a 6 2b=ð2 þ bÞ
ð10Þ
giving a cold steady-state solution [10,12]. 3.2. Warm phase (0 < s < 1, r ¼ 1 þ (d 1)s ) The warm exact steady-state solution can be easily obtained as a½1 þ ðd 1Þs 2 1 þ x2 ; 0 6 x 6 1: U ðxÞ ¼ 2 b
ð11Þ
Since U ðs Þ ¼ 1, Eq. (11) leads to the following equation for the steadystate values s of sðtÞ: 2 2 a½1 þ ðd 1Þs 1 þ s 2 ¼ 0; ð12Þ b which is solved numerically by using the Newton method and used in Eq. (11) to determine the exact steady-state temperature. Since s 2 ð0; 1, for all time, from Eq. (12) we obtain the inequality að2 þ bÞ > 2b P 2ad
ð13Þ
implies a warm steady-state which is identical to that obtained by Westbrook [10] and Wood and Kutluay [12].
480
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
3.3. Hot phase (s ¼ 1, r ¼ d) The steady-state solution to this phase is ad 2 2 1þ x : U ðxÞ ¼ 2 b
ð14Þ
Since U ð1Þ P 1 in this phase, the values of a, b and d satisfying the inequality ð15Þ
b 6 ad; giving a hot steady-state [10,12].
4. Methods of solution The interval 0 ¼ x0 < x1 < < xN ¼ 1 is divided into N equal finite elements by h ¼ 1=N . The quadratic B-splines Qm ðxÞ, at the knots xm which form a basis for functions over ½0; 1 are defined as [17] 8 2 2 2 ðxmþ2 xÞ 3ðxmþ1 xÞ þ 3ðxm xÞ ; ½xm1 ; xm ; > > < 2 2 1 ½xm ; xmþ1 ; Qm ðxÞ ¼ 2 ðxmþ2 xÞ2 3ðxmþ1 xÞ ; h > ðx xÞ ; ½x > mþ1 ; xmþ2 ; : mþ2 0; otherwise: ð16Þ So the global approximation UN ðx; tÞ to the function U ðx; tÞ can be written in terms of the splines as UN ðx; tÞ ¼
N X
Qj ðxÞdj ðtÞ;
ð17Þ
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 ;
ð18Þ
0 6 n 6 h:
Eq. (16) can be expressed as 8 2 Qm1 1 < ðh2 nÞ ; Qm ¼ 2 h þ 2hn 2n2 ; h : 2 Qmþ1 n;
0 6 n 6 h:
ð19Þ
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
481
Since all other splines are zero over the element ½xm ; xmþ1 , the approximation (17) over this element can be written in terms of the basis functions given by Eq. (19) as UN ðx; tÞ ¼
mþ1 X
Qj ðxÞdj ðtÞ:
ð20Þ
j¼m1
Using splines (19) and the approximation (20), the nodal values Um and Um0 at the node xm are given in terms of element parameters dm by Um ¼ U ðxm Þ ¼ dm1 þ dm ;
ð21Þ
2 Um0 ¼ U 0 ðxm Þ ¼ ðdm dm1 Þ; h
ð22Þ
where the prime denotes differentiation with respect to x. Applying a weighted residual method to Eq. (8) with weighting functions Wm produces the weak form Z xN Wm ðUt Uxx arÞ dx ¼ 0; ð23Þ x0
where m ¼ 0; 1; . . . ; N 1. 4.1. Subdomain collocation method In the subdomain collocation (SC) method [18], the weighting functions Wm are 1; xm 6 x 6 xmþ1 ; Wm ¼ 0; otherwise and so, the weak form (23) over each finite element ½xm ; xmþ1 yields the integral equation Z xmþ1 ðUt Uxx arÞ dx ¼ 0: ð24Þ xm
Integrating (24) and using Eqs. (20) and (22), we obtain h _ 2 ðdm1 þ 4d_ m þ d_ mþ1 Þ ðdm1 2dm þ dmþ1 Þ har ¼ 0; 3 h
ð25Þ
where the dot denotes differentiation with respect to t. Substituting the following forward difference for d_ and the weight average for d d_ ¼ ðdnþ1 dn Þ=Dt;
d ¼ ð1 hÞdn þ hdnþ1 ;
in Eq. (25) results in a system of N 1 linear algebraic equations as
ð26Þ
482
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489 nþ1 ð1 6hrÞdm1 þ ð4 þ 12hrÞdmnþ1 þ ð1 6hrÞdnþ1 mþ1
¼ ð1 þ 6ð1 hÞrÞdnm1 þ ð4 12ð1 hÞrÞdnm þ ð1 þ 6ð1 hÞrÞdnmþ1 þ 3kar;
ð27Þ
where (m ¼ 1; 2; . . . ; N 2), Dtð kÞ is the time step, r ¼ k=h2 and 0 6 h 6 1 is the weight factor. The cases h ¼ 0, h ¼ 1=2 and h ¼ 1 correspond, respectively, to the explicit, Crank–Nicolson and fully implicit finite difference methods. Using Eqs. (21) and (22), the boundary conditions (2) and (3) become dn0 ¼ dn1 ;
ð28Þ
nþ1 d0nþ1 ¼ d1
and dnN ¼ cdnN 1 ;
ð29Þ
dNnþ1 ¼ cdNnþ1 1 ;
where c ¼ ð2 hbÞ=ð2 þ hbÞ. Substituting the expressions (28) and (29) in Eq. (27) and using a Crank– Nicolson approach we obtain the resulting system of equations ð5 þ 3rÞdnþ1 þ ð1 3rÞd1nþ1 ¼ ð5 3rÞdn0 þ ð1 þ 3rÞdn1 þ 3kar; 0 3rÞdnþ1 m1
6rÞdnþ1 m
ð1 þ ð4 þ þ ð1 ¼ ð1 þ 3rÞdnm1 þ ð4 6rÞdnm þ ð1 þ 3rÞdnmþ1 þ 3kar; 3rÞdnþ1 N2
m ¼ 0;
3rÞdnþ1 mþ1
m ¼ 1; 2; . . . ; N 2;
cÞÞdnþ1 N1
ð1 þ ð4 þ c þ 3rð2 ¼ ð1 þ 3rÞdnN 2 þ ð4 þ c 3rð2 cÞÞdnN1 þ 3kar;
m ¼ N 1:
ð30Þ
This system can be easily and efficiently solved by using a variant of Thomas algorithm. 4.1.1. Stability analysis of the method For stability analysis it is convenient to use the Fourier series method pffiffiffiffiffiffiffi (see, e.g., [19], Section 2). Substituting the Fourier mode dnm ¼ gn eiumh (i ¼ 1) into Eq. (30) leads to the growth factor g of the form g¼
2ð1 þ 3rÞ cos2 ðuh=2Þ þ 1 6r : 2ð1 3rÞ cos2 ðuh=2Þ þ 1 þ 6r
For the stability, g must be jgj 6 1 which is satisfied for all positive values of stability parameter r only. So, the system (30) is unconditionally stable.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
483
4.2. Petrov–Galerkin method Applying a Petrov–Galerkin (PG) method with the weighting functions Wm defined by linear B-spline basis functions [17] 8 ðx xÞ 2ðxm xÞ; ½xm1 ; xm ; 1 < mþ1 Wm ðxÞ ¼ xmþ1 x; ½xm ; xmþ1 ; ð31Þ h: 0; otherwise; the weak form (23), over each finite element ½xm ; xmþ1 , with some manipulation, leads to the integral equation Z xmþ1 ðWUt þ Wx Ux W arÞ dx ¼ bW ðxN ÞU ðxN ; tÞ: ð32Þ xm
Using the local coordinate transformation (18), from Eq. (31) we obtain the linear B-spline functions for the finite element ½xm ; xmþ1 Wm ¼ 1 n=h;
Wmþ1 ¼ n=h:
Substituting Eq. (20) in Eq. (32) leads to mþ1 Z xmþ1 mþ1 Z xmþ1 X X Wi Qj dx d_ ej þ Wi 0 Q0j dx dej j¼m1
xm
xm
j¼m1
Z mþ1 X e þb Wi ðxN ÞQj ðxN Þdj ¼ ar j¼m1
xmþ1
Wi dx ;
xm
which also be written in a matrix form as Ae d_ e þ ðBe þ C e Þde ¼ d e ;
ð33Þ
where 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 Wi Qj dx; Aeij ¼ xm Z xmþ1 Wi 0 Q0j dx; Beij ¼ ð34Þ xm e Cij ¼ bWi ðxN ÞQj ðxN Þ; Z xmþ1 die ¼ ar Wi dx; xm
where i and j takes only the values, respectively m, m þ 1 and m 1, m, m þ 1 for the element ½xm ; xmþ1 . The element matrices (34) can be obtained, with some manipulation, as
484
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
Aeij
h 3 ¼ 12 1
8 8
1 ; 3
Beij
1 1 ¼ h 1
0 0
1 ; 1
die
ar 1 ¼ 2 1
and Cije is a zero matrix but 0 0 0 N 1 : Cij ¼ 0 b b Assembling all contributions from all elements, Eq. (33) leads to the matrix differential equation Ad_ þ ðB þ CÞd ¼ d;
ð35Þ T
where d ¼ ðd1 ; d0 ; . . . ; dN Þ and A, B and C are ðN þ 1Þ ðN þ 2Þ matrices and d is an ðN þ 1Þ 1 matrix which are defined from the element matrices Ae , Be , C e and d e , respectively. Using a Crank–Nicolson approach for d and forward finite difference for d_ given by Eq. (26), Eq. (35) takes the form Dt Dt A þ ðB þ CÞ dnþ1 ¼ A ðB þ CÞ dn þ dDt: ð36Þ 2 2 Using one of the boundary conditions (2) and (3) in the system (36) leads to an ðN þ 1Þ ðN þ 1Þ matrix system which can be easily solved by the Thomas algorithm. A typical member of Eq. (36) in terms of the nodal parameters dnm is nþ1 nþ1 nþ1 þ a2 dm1 þ a2 dnþ1 a1 dm2 m þ a1 dmþ1
¼ a3 dnm2 þ a4 dnm1 þ a4 dnm þ a3 dnmþ1 ;
ð37Þ
where a1 ¼
h k ; 12 2h
a2 ¼
11h k þ ; 12 2h
a3 ¼
h k þ ; 12 2h
a4 ¼
11h k : 12 2h
4.2.1. Stability analysis of the method For stability analysis, again substituting the Fourier mode dnm ¼ gn eiumh , pffiffiffiffiffiffi ffi (i ¼ 1) in Eq. (37) leads to the growth factor g of the form g¼
a ib ; c id
where a ¼ a3 cos 2uh þ h cos uh þ a4 ;
b ¼ a3 sin 2uh ða3 a4 Þ sin uh;
c ¼ a1 cos 2uh þ h cos uh þ a2 ;
d ¼ a1 sin 2uh ða1 a2 Þ sin uh:
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
485
For the stability, the magnitude of the growth factor g must be jgj 6 1 which is mathematically equivalent to c2 þ d 2 a2 b2 P 0. By virtue of any algebraic package programme it is seen that c2 þ d 2 a2 b2 P 0 for h > 0, k > 0 and u > 0. Thus, the scheme (37) is unconditionally stable.
5. Numerical results and conclusion In order to evaluate the temperature profiles we use the values of a, b and d satisfying inequalities (10), (13) and (15) for the cold, warm and hot phases, respectively. We also use the steady-state test
1 max jU0new U0old j; jUNnew UNold j 6 0:5 108 ; k
Nh ¼ 1
ð38Þ
for any one of the three phases. All calculations were performed in Fortran code under a Pentium 4 processor using double precision arithmetic. To test the accuracy of both numerical methods we have calculated the kek1 norm defined as N 1 X U ðxm Þ ; kek1 ¼ 1 N þ 1 m¼0 U ðxm Þ
e ¼ ½e0 ; e1 ; . . . ; eN T :
We insistly take the same values of a, b and d for each particular phase in case the reader easily compares the results obtained using the present methods with earlier works [12,13,16]. 5.1. Cold steady-state The values a ¼ 0:1 and b ¼ 0:2 satisfying inequality (10) are sufficient to obtain a cold steady-state solution. Such a choice of the values of a and b is only typical of a problem evolving to a cold steady-state. Alternatively, different values of a and b, still satisfying the inequality (10), simple serve to alter the time taken to attain steady-state. The results obtained by the present methods (for k ¼ 0:01 and h ¼ 0:1) were reported when the steady-state test (38) is satisfied, and presented in Table 1 with the exact solution (9). It is observed that the numerical solutions are seen to be in exact agreement with the exact one. The only remarkable difference is the time to attain the steady-state. The subdomain collocation solution is obtained after something like 89.9 units of CPU time whereas the Petrov– Galerkin solution is obtained after some 65.32 units of CPU time.
486
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
Table 1 Solutions in cold steady-state x
U SC, PG 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.5500 0.5495 0.5480 0.5455 0.5420 0.5375 0.5320 0.5255 0.5180 0.5095 0.5000
5.2. Warm steady-state If the predicted temperature at a node is such that at time level n, U ðxm Þ > 1 and U ðxmþ1 Þ < 1 in the warm phase, the moving boundary s lies between adjacent values U ðxm Þ and U ðxmþ1 Þ. Since U ðx; tÞ ¼ 1 at x ¼ sðtÞ, the location of the interface s at time level n can be easily estimated by using a linear interpolation as [16] sn ¼ xm þ
hð1 U ðxm ÞÞ : U ðxmþ1 Þ U ðxm Þ
ð39Þ
Thus, the motion of sðtÞ can be obtained from s_ ¼
snþ1 sn : Dt
For a warm steady-state solution, the values a ¼ 0:5 and b ¼ 0:2 satisfying inequality (13) are sufficient. The numerical results obtained by the both methods (for k ¼ 0:01 and h ¼ 0:1, 0.05, 0.025, 0.0125) are compared with the exact solution (11) in Table 2. It is seen that numerical solutions are in good agreement with the exact one, and exhibit the expected convergence as the mesh size h is refined. The predicted location of the interface s is obtained by Eq. (39) which is in good agreement with the exact one given by Eq. (12). Fig. 3 shows the movement (location and speed) of the interface s. In both numerical methods, the boundary s smoothly reaches its steady-state value of about 0.62 in some 5.2 units of CPU time.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
487
Table 2 Solutions in warm steady-state x
U
Exact
SC and PG h ¼ 0:1
h ¼ 0:05
h ¼ 0:025
h ¼ 0:0125
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.036749 1.035806 1.032979 1.028266 1.021669 1.013186 1.002819 0.990566 0.976429 0.960406 0.942499
1.036641 1.035699 1.032871 1.028160 1.021563 1.013081 1.002715 0.990464 0.976328 0.960307 0.942401
1.036587 1.035645 1.032818 1.028106 1.021510 1.013028 1.002663 0.990412 0.976277 0.960257 0.942352
1.036585 1.035643 1.032816 1.028104 1.021507 1.013026 1.002660 0.990410 0.976275 0.960255 0.942350
kek1
0.000160
0.000056
0.000004
0.000001
1.036583 1.035641 1.032814 1.028102 1.021506 1.013024 1.002659 0.990408 0.976273 0.960253 0.942348
Fig. 3. The movement of the interface sðtÞ.
5.3. Hot steady-state Both numerical solutions to this phase are easily obtained for the values a ¼ 105 and b ¼ 0:2 satisfying inequality (15). The agreement between our numerical results (for k ¼ 0:01 and h ¼ 0:1) and the exact solution (14) match each other, shown in Table 3. Again, only difference is the time to attain steady-state. The subdomain collocation solution is obtained after something
488
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
Table 3 Solutions in hot steady-state x
U SC, PG and exact
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
5.500 5.495 5.480 5.455 5.420 5.375 5.320 5.255 5.180 5.095 5.000
like 102.19 units of CPU time whereas the Petrov–Galerkin solution is obtained after some 77.61 units of CPU time. In order to show how good the numerical predictions exhibit the correct physical behaviour of the problem, we only give graph in Fig. 4. The numerical solutions (for k ¼ 0:01 and h ¼ 0:025) obtained by the present methods in the cold, warm and hot phases during the steady-state are drawn on the same diagram, but curves cannot be distinguishable since they are very close to each other.
6.0 5.5 5.0
Temperature, U
4.5
Hot
4.0 3.5 3.0 2.5 2.0 1.5 1.0
Warm Cold
0.5 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Distance, x Fig. 4. Temperature profiles in the cold, warm and hot phases at steady-state.
S. Kutluay, A. Esen / Appl. Math. Comput. 162 (2005) 475–489
489
In this paper, we have presented approximate steady-state solutions of a one-dimensional PTC thermistor problem with a temperature dependent step function electrical conductivity using subdomain collocation and Petrov– Galerkin methods. It is shown that the numerical solutions obtained by the present methods satisfy all the physical characteristics of the problem. As a conclusion, both methods are capable of producing numerical solutions for the PTC thermistor problem. Further, the present methods can be applied to a two-dimensional analogue of the problem.
References [1] F.J. Hyde, Thermistors, Butterwoth Group (ILIFFE books), London, 1971. [2] E.D. Macklen, Thermistors, Electrochemical Publications, Scotland, 1979. [3] G. Cimatti, Remark on existence and uniqueness for the thermistor problem under mixed boundary conditions, Q. Appl. Math. 47 (1989) 117–121. [4] 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. € [5] H. Diesselhorst, Uber das probleme eines elektrisch erwarmten leiters, Annal. Phys. 1 (1900) 312–325. [6] A. Fowler, I. Frigaard, S.D. Howison, Temperature surges in thermistors, SIAM J. Appl. Math. 52 (1992) 998–1011. [7] S.D. Howison, J.F. Rodrigues, M. Shillor, Stationary solutions to the thermistor problem, J. Math. Anal. Appl. 174 (1993) 573–588. [8] J. Wiedmann, The thermistor problem, Nonlinear Differ. Equ. Appl. 4 (1997) 133–148. [9] S. Zhou, D.R. Westbrook, Numerical solutions of the thermistor equations, J. Comput. Appl. Math. 79 (1997) 101–118. [10] D.R. Westbrook, The thermistor: a problem in heat and current flow, Numer. Meth. PDEs 5 (1989) 259–273. [11] 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. [12] A.S. Wood, S. Kutluay, A heat balance integral model of the thermistor, Int. J. Heat Mass Transfer 38 (10) (1995) 1831–1840. [13] 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. [14] S. Kutluay, A.S. Wood, Numerical solutions of the thermistor problem with a ramp electrical conductivity, Appl. Math. Comput. 148 (2004) 145–162. [15] A.S. Wood, Modelling the thermistor, in: M. Heili€ o (Ed.), Proceedings of the 5th European Conference on Mathematics in Industry, Kluwer Academic, Stuttgart, 1991, pp. 397–400. [16] 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. [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.