Multiplicity Analysis of a Thermistor Problem

1. Analysis

1.2. Energy Balance

Consider a horizontal cylindrical segment of a conductor of uniform material density γ with constant thermal conductivity k. The segment has diameter D and length L as it is schematically shown in Figure 1. Heat is being dissipated by conduction through the core of the device and by natural convection and radiation through the lateral surface area, in an ambient environment of constant temperature $T_{\infty }$. An energy balance along the longitudinal direction X yields the following differential equation for the device temperature T:

$$\gamma C\frac{\partial T}{\partial t}=\frac{\partial }{\partial X}\left(k\frac{\partial T}{\partial X}\right)-\frac{P}{A}\left[h_c(T-T_{\infty })+\epsilon \sigma (T^4-T_{\infty }^4)\right]+EJ$$

(1)

In the equation above, C is the specific heat capacity, A is the cross sectional area, P is the perimeter, $h_c$ is the convective heat transfer coefficient, ε is the surface emissivity, σ is the Stefan-Boltzmann constant, E is the electric field intensity and J is the current density through the device.

1.3. Electric Resistivity

A characteristic feature of a ceramic PTC device is the strongly non-linear dependence of its resistivity with respect to temperature. Driven by a transition of the ferroelectric PTC material the resistance increases several orders of magnitude in a relatively small temperature interval, for instance between 100°C and 200°C. The same smooth curve with continuous derivatives, with respect to temperature for the subsequent numerical bifurcation analysis has been adopted from Karpov [18], which represents a barium titanate (BaTiO₃) based device:

$$\widehat{\rho }(T)={\widehat{\rho }}_l+\frac{1}{{({\widehat{\rho }}_h-{\widehat{\rho }}_l)}^{-1}+\exp [-0.12(T-95)]}$$

(2)

where ${\widehat{\rho }}_l=2$Ωm and ${\widehat{\rho }}_h={10}^4$Ωm are the asymptotic resistivity values corresponding to the low and high operating temperatures respectively. The temperature of 95°C signifies the onset of the transition from the low to the high resistivity value as the contribution of the exponential term in Equation (2) becomes significant. The form of Equation (2) is supported by a significant volume of measurements as for instance reported by Brzozowski and Castro [27], Wang et al. [28], Wang et al. [29], Luo et al. [30], Takeda et al. [31], Rowlands and Vaidhyanathan [32].

1.4. Heat Transfer Model

The heat generated in the device due to the current flow (Joule heating, Metaxas [33], Lupi [34], Lupi et al. [35]) is dissipated to the surrounding environment through natural convection and radiation exchange. For the circumferential average Nusselt number Nu, the correlation of Churchill and Chu [36] is employed:

$$\mathrm{Nu}=0.36+0.518\frac{{\mathrm{Ra}}^{1/4}}{f(\mathrm{Pr})},f(\mathrm{Pr})={\left[1+{\left(\frac{0.559}{\mathrm{Pr}}\right)}^{9/16}\right]}^{4/9}$$

(3)

where $f(\mathrm{Pr})$ is a weak function of the Prandtl number Pr. Equation (3) covers a very wide range of Rayleigh numbers namely in the range from 10^-6 to 10⁹, while it maintains a simple and compact mathematical form. Equation (3) is applied locally in the evaluation of the convective heat transfer coefficient along the device axis, in a similar manner as it was utilized by Faghri and Sparrow [37]. Consequently the local Rayleigh number Ra is evaluated from the local temperature difference as:

$$\mathrm{Ra}=g\beta D^3[T(X)-T_{\infty }]/\alpha \nu$$

(4)

where g is the acceleration due to gravity, β is the thermal expansivity, D is the device diameter, α is the thermal diffusivity and ν is the kinematic viscosity.

1.5. Boundary Conditions. Problems P1 and P2.

As it will be described in the next paragraphs the boundary conditions have a profound effect on the bifurcation structure in general and on the temperature distribution in particular. Thus for problem P1 the edges of the device are considered adiabatic:

$${\frac{\partial T}{\partial X}|}_{X=0}={\frac{\partial T}{\partial X}|}_{X=L}=0$$

(5)

whereas for problem P2 the imposed boundary conditions are:

$${T|}_{X=L}=T_e,{\frac{\partial T}{\partial X}|}_{X=L/2}=0$$

(6)

1.6. The Electrothermal Model in Dimensionless Form

Considering a constant (dc) current flowing through the device, the electric field intensity is related to the current density through Ohm’s law, $E=\widehat{\rho }(T)J$, Metaxas [33], Lupi [34], Lupi et al. [35]. Introducing dimensionless variables

$$x=X/L,\tau ={\alpha t/L}^2,\Theta =T/T_{\infty },\rho =\widehat{\rho }/{\widehat{\rho }}_{\mathrm{ref}}$$

(7)

the temperature distribution along the device takes the form:

$$\frac{\partial \Theta }{\partial \tau }=\frac{{\partial }^2\Theta }{\partial x^2}-u^2\left[\mathrm{Nu}(\Theta -1)+C_h({\Theta }^4-1)-j^2\rho \right]$$

(8)

In the above equation, the conduction- convection parameter (CCP) is defined as:

$$u^2=\frac{h_{\mathrm{ref}}{L}^2}{k(A/P)}$$

(9)

where the reference heat transfer coefficient $h_{\mathrm{ref}}$ is defined through the Nusselt number:

$$\mathrm{Nu}=\frac{h_c}{k_{\infty }/D}=\frac{h_c}{h_{\mathrm{ref}}}$$

In terms of the dimensionless variables defined above the local Rayleigh number becomes:

$$\mathrm{Ra}={\mathrm{Ra}}_{\infty }(\Theta -1),{\mathrm{Ra}}_{\infty }=\frac{g\beta D^3{T}_{\infty }}{\alpha \nu }$$

The current density parameter is related to the current density as below

$$j^2=\frac{J^2(A/P)}{T_{\infty }}{\left(\frac{\widehat{\rho }}{h}\right)}_{\mathrm{ref}}$$

(10)

and $C_h$ is the ratio of the radiative heat transfer coefficient to the reference heat transfer coefficient

$$C_h=\frac{\epsilon \sigma T_{\infty }^3}{h_{\mathrm{ref}}}=\frac{h_r}{h_{\mathrm{ref}}}$$

(11)

Under steady state conditions the partial differential equation Equation (8) reduces to a two point boundary value problem for the device temperature $\Theta (x)$

$${\Theta }^{″}-u^2\left[\mathrm{Nu}\right(\Theta -1)+C_h({\Theta }^4-1)-j^2\rho ]=0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0<x<1$$

(12)

with boundary conditions

$${\Theta }^{\prime }(x=0)={\Theta }^{\prime }(x=1)=0$$

(13)

for problem P1 and

$$\Theta (x=1)={\Theta }_e,{\Theta }^{\prime }(x=0.5)=0$$

(14)

for problem P2.

1.7. Stability

The stability of a certain steady state ${\Theta }_s(x)$ to small perturbations $\vartheta (x)$ i.e.

$$\Theta (x,\tau )={\Theta }_s(x)+\vartheta (x)e^{\lambda \tau }$$

(15)

is determined by the eigenvalues λ of the corresponding Sturm-Liouville problem, after substituting Equation (15) into the original partial differential equation, Equation (8),

$${\vartheta }^{″}-\left[u^2\Delta Q_{\Theta }+\lambda \right]\vartheta =0,0<x<1$$

(16)

where

$$\Delta Q_{\Theta }=\frac{\partial }{\partial \Theta }{\left[\mathrm{Nu}(\Theta -1)+C_h({\Theta }^4-1)-j^2\rho (\Theta )\right]}_{\Theta ={\Theta }_s}$$

The corresponding boundary conditions for problem P1 are ${\vartheta }^{\prime }(0)={\vartheta }^{\prime }(1)=0$ and for problem P2 are ${\vartheta }^{\prime }(0.5)=\vartheta (1)=0$. During branch tracing, for every steady state that has been calculated from Equation (12), the associated eigenvalue problem, Equation (16), is subsequently numerically solved and a sufficient number of eigenvalues is determined. Stable solutions are characterized by negative eigenvalues whereas positive ones correspond to unstable temperature distributions.

2. Results and Discussion

The second order, two-point boundary value problem described by Equation (12) is transformed into a system of first order equations through the transformation ${\Theta }_1=\Theta$, ${\Theta }_2={\Theta }^{\prime }$ and solved numerically. In order to ensure an accurate and reliable numerical solution, two different methods have been employed resulting in identical results under the strict tolerances imposed. The first one utilizes a multi-shooting Runge-Kutta formula pair of order 8(7), Hairer et al. [38] and the second one a spline collocation method described by Ascher et al. [39]. Continuation along the various branches has been carried out along the lines suggested by Seydel [40]. For the computation of the singular points an extended problem is being formed from the partial derivatives of Equation (12) with respect to the parameters, according to Witmer et al. [41].

Before we analyze the complete numerical solution, it is instructive to discuss first the uniform solutions of Equation (12), which reduces to an algebraic one for a constant temperature profile:

$$\mathrm{Nu}(\Theta -1)+C_h({\Theta }^4-1)-j^2\rho (\Theta )=0$$

(17)

A geometrical (graphical) solution is depicted in Figure 2 where the heat generation and the heat dissipation curves are being plotted for a variety of current parameters $j^2$ and reference Rayleigh numbers ${\mathrm{Ra}}_{\infty }$ (Figure 2a). In Figure 2b the effect of radiation through $C_h$ on the heat rejection rate is demonstrated. Depending on the combination of $j^2$, ${\mathrm{Ra}}_{\infty }$ and $C_h$ up to three solutions of Equation (17) may be obtained from the number of the intersection points.

2.1. Problem P1

Interestingly, when the conduction term and the associated conduction-convection parameter u is taken into consideration, a far more complicated solution structure and multiplicity pattern emerges as it is shown in Figure 3, where the edge temperature ${\Theta }_e=\Theta (0)$ is plotted against $j^2$. For lower values of the conduction-convection parameter $(u=1)$ the three uniform solutions of Equation (17) are being recovered (Figure 3a). As u increases the number of solutions increases as well. Five for $u=2$ in Figure 3b, seven for $u=3$ in Figure 3c and eleven for $u=5$ in Figure 3d. The corresponding temperature profiles are shown in Figure 4 (a to d). In every case the solution structure is consisting of the three uniform solutions, one stable “cold”, one stable “hot” and one unstable at an intermediate temperature. Additional unstable solutions emerge in the form of standing waves as u gradually increases. It is worth pointing out that the solution obtained by imposing boundary conditions Equation (13) have two salient features. The first one is that as long as the nonuniform solutions are unstable only the uniform ones are physically encountered, the “cold” one being far more preferable from the operating point of view since it results in reduced thermal stresses and thermal loading (reduced fracture probability) of the device. The second one stems from the first and the associated stability analysis, namely the temperature profiles may be obtained from the solution of the simpler algebraic Equation (17) instead of solving the complete boundary value problem, Equation (12). In other words effective control of the boundary conditions (and the edge temperature, say through the application of an insulating layer) appears to have a profound effect on the operating and performance characteristics of the device.

2.2. Problem P2

For the case where the edge temperature is fixed, the bifurcation analysis reveals a rich and interesting pattern, shown in Figure 5. Selecting the center temperature ${\Theta }_c=\Theta (0.5)$ as the bifurcation parameter, the projection of the limit points on the $({\Theta }_c,j^2)$ plane forms a pattern of nested cusp points $(C_1{C}_2{C}_3\mathrm{...})$ in Figure 5a, where the regions with unique, three, five, seven (and so on) solutions are designated. For a better understanding of the complexity of the solution structure and the clarity of the exposition a geometrical perspective of the nested cusps is depicted in Figure 6. Again the CCP has a profound effect on the solution structure since the number of the steady states calculated depend on its magnitude. Indeed starting from $u=0.5$ and a unique solution in Figure 5b, three solutions emerge for $u=1$ in Figure 5c, five for $u=1.5$ in Figure 5d and seven in Figure 5e for $u=2$. The temperature profiles for problem P2 are symmetrical with respect to the center of the device. For the unique solution calculated in Figure 7a the center is maintained at a lower temperature, while it increases towards the edge $({\Theta }_c<{\Theta }_e)$. When the CCP increases and three solutions are present as in Figure 7b, for the same edge temperature ${\Theta }_e$ two stable solutions exist one “cold” $({\Theta }_c<{\Theta }_e)$ and one “hot” $({\Theta }_c>{\Theta }_e)$. In contrast to the “cold” solution described earlier, the peak temperature for the “hot” one now appears around the center and gradually decreases to ${\Theta }_e$. As CCP further increases more solutions emerge, five in Figure 7c for $u=1.5$ and seven in Figure 7d when $u=2$. Yet the stable ones remain the “cold” and “hot” set described earlier whereas the ones appearing as standing waves are unstable. In other words the bistability is retained as the CCP increases and several solutions emerge. Furthermore, as u increases a temperature plateau around the center is being formed and the temperature gradient along the device although it cannot be entirely eliminated as it was the case in problem P1, it is definitely reduced. It is worth pointing out that a multiplicity structure consisting of nested cusps is not a unique feature of this particular electrothermal problem. Aris [42] was the first to publish imbedded cusps in the study of first order reaction in spherical pellet. The same solution structure emerged in the analysis of another chemical reacting system, namely an isothermal Langmuir-Hinshelwood reaction in a cylindrical or spherical pellet, Witmer et al. [43]. Another example of similar multiplicity behavior is encountered when the three boiling modes (nucleate, transition and film) are being excited by a uniform heat generating source along a non-isothermal extended surface immersed in a boiling liquid, Krikkis [44].

3. Conclusions

An electrothermal model for a barium titanate based thermistor has been set up, featuring a nonlinear temperature dependent natural convection combined with radiation heat rejection mechanism through the device lateral surface. A smooth and strongly temperature dependent electric conductivity function based on experimental data has been adopted. As a result the problem formulated in this way admits multiple steady state solutions which do not depend on the external circuit.

Important findings are the profound effect of the conduction-convection parameter and of the boundary conditions on the multiplicity structure. Regardless of the type of the boundary conditions imposed, the number of the multiple steady states depends on the magnitude of the conduction-convection parameter, as for instance up to eleven solutions have been calculated when the edges of the device are insulated (problem P1) and up to seven when the edge temperature is being fixed (problem P2). The stability analysis reveals that for problem P1 only the uniform solutions are stable, namely one “cold” and one “hot”. Therefore, control and effective reduction of the heat transmitted through the edges results in a reduced thermal loading of the device both in the steady state operation and during transients since the temperature gradient across the device is vanishing.