Multiplicity Analysis of a Thermistor Problem—A Possible Explanation of Delamination Fracture

Abstract

In the present study, a numerical bifurcation analysis of a PTC thermistor problem is carried out, considering a realistic heat dissipation mechanism due to conduction, nonlinear temperature-dependent natural convection, and radiation. The electric conductivity is modeled as a strongly nonlinear and smooth function of the temperature between two limiting values, based on measurements. The temperature field has been resolved for both cases were either the current or the voltage (nonlocal problem) is the controlling parameter. With the aid of an efficient continuation algorithm, multiple steady-state solutions that do not depend on the external circuit have been identified as a result of the inherent nonlinearities. The analysis reveals that the conduction–convection parameter and the type of the imposed boundary conditions have a profound effect on the solution structure and the temperature profiles. For the case of current control, depending on the boundary conditions, a complex and interesting multiplicity pattern appears either as a series of nested cusp points or as enclosed branches emanating from pitchfork bifurcation points. The stability analysis reveals that when the device edges are insulated, only the uniform solutions are stable, namely, one “cold” and one “hot”. A key feature of the “hot” state is that the corresponding temperature is proportional to the input power and its magnitude could be one or even two orders of magnitude higher than the “cold” one. Therefore, the change over from the “cold” to the “hot” state induces a thermal shock and could perhaps be the reason for the mechanical failure (delamination fracture) of PTC thermistors.

1. Introduction

Thermistors are thermally sensitive resistors and have either a negative (NTC) or positive (PTC) resistance/temperature coefficient. We will only discuss the latter, albeit certain similarities and differences, especially for the flash sintering phenomenon, will be addressed in the appropriate sections. Its characteristic feature is the strongly dependent electric conductivity. PTCs are manufactured from silicon, barium, lead, and strontium titanates with the addition of yttrium, manganese, tantalium, and silica [1,2,3,4,5]. PTC thermistors are widely used as current-limiting devices, that is, as non-destructible (resettable) fuses for electric circuit protection, sensing excessive currents. They can also be encountered as a micro self-heating thermostat for microelectronic, biomedical, and chemical applications. Common geometrical configurations are the washer, the disk, and the rod type [1,2,3]. Although PTCs and in general electroceramics are, in principle, loaded electrically, a significant number of mechanical failures are recorded annually. This may be explained on the basis of the Joule self-heating effect, which causes temperature differences, thermal strains, and excessive thermo-mechanical stresses that may cause the failure of the device. As the current technological trends point towards greater device miniaturization while operating at higher power densities, there is an increasing demand for a thorough analysis and understanding of the underlying coupled electrothermal phenomena in electroceramic devices (Dewitte et al. [6], Supancic [7], Danzer et al. [8]).

The thermistor, as a coupled thermo-electric problem, has attracted significant attention. An early result (1900) described by Diesselhorst [9] for the steady-state problem shows that the temperature may be expressed as a function of the electric potential provided that certain types of boundary conditions are imposed. Cimatti [10] extended the result of Diesselhorst to obtain existence and uniqueness conditions for the steady-state problem. Further results on existence and uniqueness were obtained by Cimatti [11], Cimatti and Prodi [12], Xie and Allegretto [13], and Antontsev and Chipot [14]. Bahadir [15,16] and Çatal [17] employed finite element numerical techniques to solve the thermistor problem assuming a step electric conductivity function. Kutluay et al. [18] obtained a heat balance integral solution of the same problem considering a modified electric conductivity function. Ammi and Torres [19] numerically solved a nonlocal parabolic equation in time and space domains resulting from the thermistor problem. Golosnoy and Sykulski [20] compared various computational techniques for coupled nonlinear thermo-electric problems. The thermistor problem was used as a test case with an electric conductivity being a nonlinear function of the temperature and the electric field intensity. Karpov [21] demonstrated the bistability conditions, the switching autowave properties, and the emergence of dissipative structures of essentially a thermistor problem under current control. Apart from the temperature-dependent electric conductivity, convective heat dissipation with a constant heat transfer coefficient over the lateral surface was assumed. It is worth noting that the thermistor problem is closely associated with the flash sintering of ceramics [22] such as, for instance, yttria-stabilized zirconia, magnesia-doped alumina, and strontium titanate, among others [23]. The essence of the process is that when an operating parameter such as the furnace temperature exceeds the corresponding limit point, established by the applied voltage that separates the stable from the unstable steady states—the Joule heating greatly exceeds the heat dissipation mechanism due to radiation and the temperature increases significantly. The process controller is then switched from voltage to current control to maintain the temperature within the specified limits [24,25].

From the literature review above, it appears that the thermistor problem has been studied with various assumptions and/or restrictions related primarily to the form of electrical conductivity, the heat dissipation mechanism, and, in certain cases, with the influence of the external electric circuit. The latter is also associated with the existence of multiple steady-state solutions, up to three, as determined from the number of intersection points between the current–voltage characteristic curves of the external (linear) circuit and the thermistor (Fowler et al. [26], Howison et al. [27], Zhou and Westbrook [28], Cimatti [29]). The aim of the present study is to provide a reasonable explanation for the most common reason PTC thermistors fail—namely, delamination fracture due to excessive thermal loading (shock). To this end, a one-dimensional PTC device model based on nonlinear electric resistivity combined with a distributed heat dissipation mechanism due to conduction, nonlinear temperature-dependent natural convection, and radiation has been developed. The problem formulated in this way admits multiple steady-state solutions that do not depend on the external circuit. The numerical bifurcation analysis reveals that the PTC thermistor is a bistable system that can operate between a “cold” and a “hot” stable state. The changeover from the “cold” to the “hot” state induces a kind of thermal shock that could be a reason for potential failure, since the “hot” state is associated with a very high temperature. This hypothesis is further supported by a similar behavior that is encountered in other bistable systems such as superconductors and boiling wires. Depending on the boundary conditions, a complex and interesting multiplicity pattern appears either as a series of nested cusp points or as enclosed branches emanating from pitchfork bifurcation points.

2. Analysis

2.1. Energy Balance

Consider a horizontal cylindrical segment of a conductor of uniform material density γ with constant thermal conductivity k. The segment has a diameter D and length L, as 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 partial 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)

Here, 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. 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 [30], Lupi [31], Lupi et al. [32]). Thus, when the current is controlled, Equation (1) reads:

$$\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]+\widehat{\rho }{J}^2.$$

(2)

On the other hand, in many practical applications, the applied voltage across the device is the controlling parameter. This may be taken into account by introducing the electric potential $\Phi$, as in [33]:

$$E=-\frac{d\Phi }{dX}=\widehat{\rho }J.$$

(3)

Integrating the above relationship and considering that the current remains constant, albeit still unknown, yields:

$$J{\displaystyle {\int }_{ 0}^{ L}\widehat{\rho }(T)dX}=-\left[\Phi (L)-\Phi (0)\right]=V,$$

(4)

where V is the voltage drop across the device, as shown in Figure 1. Substituting the current density J from Equation (4) into Equation (2), the energy balance for voltage control takes the form:

$$\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]+\frac{\widehat{\rho }{V}^2}{{\left[{\displaystyle {\int }_{ 0}^{ L}\widehat{\rho }(T)dX}\right]}^2}.$$

(5)

Compared with the current control problem, this is a nonlocal problem, since the solution depends on the resistivity integral over the device.

2.2. Electric Resistivity

A characteristic feature of a ceramic PTC device is the strongly nonlinear 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 [21], 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)]},$$

(6)

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 (6) becomes significant. The form of Equation (6) is supported by a significant volume of measurements, as, for instance, those reported by Brzozowski and Castro [34], Wang et al. [35], Wang et al. [36], Luo et al. [37], Takeda et al. [38], and Rowlands and Vaidhyanathan [39].

2.3. Heat Transfer Model

The heat generated in the device due to the current flow (Metaxas [30], Lupi [31], Lupi et al. [32]) 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 [40] 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},$$

(7)

where $f(\mathrm{Pr})$ is a weak function of the Prandtl number Pr. Equation (7) covers a very wide range of Rayleigh numbers, namely, in the range from 10⁻⁶ to 10⁹, while it maintains a simple and compact mathematical form. Equation (7) is applied locally in the evaluation of the convective heat transfer coefficient along the device axis (see Figure 1), in a similar manner as it was utilized by Faghri and Sparrow [41]. 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 ,$$

(8)

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.

2.4. Boundary Conditions: Problems P1 and P2

As 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,$$

(9)

whereas for problem P2, the imposed symmetrical boundary conditions are:

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

(10)

2.5. Electrothermal Model in Dimensionless Form

Introducing dimensionless variables

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

(11)

the temperature distribution along the device takes the form of

$$\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],$$

(12)

for current control and

$$\frac{\partial \Theta }{\partial \tau }=\frac{{\partial }^2\Theta }{\partial x^2}-u^2\left[\mathrm{Nu}(\Theta -1)+C_h({\Theta }^4-1)-{\left(\frac{v}{u{\displaystyle {\int }_{ 0}^{ 1}\rho dx}}\right)}^2\rho \right],$$

(13)

for voltage control. In the above equations, u is the conduction–convection parameter (CCP), which is extensively used in conjugate heat transfer and electro-thermal problems associated with superconductors and metallic conductors as well [42,43,44], and is defined as:

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

(14)

where the reference heat transfer coefficient $h_{\mathrm{ref}}$ is conveniently 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}},$$

(15)

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}}}.$$

(16)

Under steady-state conditions, the partial differential Equation (12) reduces to a two-point boundary value problem for current control:

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

(17)

and, similarly, Equation (13) for the voltage control:

$${\Theta }^{″}-u^2\left[\mathrm{Nu}(\Theta -1)+C_h({\Theta }^4-1)-{\left(\frac{v}{u\overline{\rho }}\right)}^2\rho \right]=0 0<x<1,$$

(18)

where

$$\overline{\rho }={\displaystyle {\int }_{ 0}^{ 1}\rho dx}.$$

The voltage–current relationship becomes

$$uj{\displaystyle {\int }_{ 0}^{ 1}\rho (\Theta )dx}-v=0,$$

(19)

and the boundary conditions are

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

(20)

for problem P1 and

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

(21)

for problem P2. Convenient reference values have been adopted for the electric resistivity ${\widehat{\rho }}_{\mathrm{ref}}={\widehat{\rho }}_h$ and the applied voltage $V_{\mathrm{ref}}^2=k{T}_{\infty }{\widehat{\rho }}_{\mathrm{ref}}$

2.6. 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 },$$

(22)

is determined by the eigenvalues λ of the corresponding Sturm–Liouville problem, after substituting Equation (22) into the original partial differential equation, Equation (12):

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

(23)

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 (23), 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.

3. Results and Discussion

The second-order, two-point boundary value problem described by Equation (17) for current control and Equation (18) for voltage control are 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. [45]) and the second one uses a spline collocation method described by Ascher et al. [46]. Continuation along the various branches has been carried out along the lines suggested by Seydel [47]. For the computation of the singular points, an extended problem is formed from the partial derivatives of Equations (17) and (18) with respect to the parameters, according to Witmer et al. [48].

Before we analyze the complete numerical solution, it is instructive to first discuss the uniform solutions of Equation (17), 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.$$

(24)

A geometrical (graphical) solution is depicted in Figure 2, where the heat generation and the heat dissipation curves are 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 (24) may be obtained from the number of the intersection points between the heat generation and the heat dissipation curves, as shown in Figure 3. Three uniform solutions exist: two are stable, namely, the “cold” and the “hot”, whereas the third one at an intermediate temperature is unstable. From this simplified analysis, a few important observations may be pointed out. The “cold” temperature remains practically constant and very close to the ambient temperature for a wide range of current densities. On the other hand, the “hot” temperature appears proportional to the input power (i.e., $\sim j^2$), since at this temperature range, the electric resistivity approaches its maximum value, i.e., $\rho \sim 1$. Consequently, the “hot” temperature can increase one and even two orders of magnitude compared to the “cold” one. Thus, when a current instability or disturbance gradually propagates through the circuit and the limit point to right is exceeded, as schematically shown in Figure 3, then the “hot” state prevails, and depending on the operating conditions, a significantly higher temperature will be established even if the radiation contribution is significant. This is also known as a hysteresis loop, and under these circumstances, it closely resembles thermal shock. However, this is not a thermal runaway, as is the case during flash sintering, since the temperature could be excessive, but it is still bounded. Therefore, it is reasonable to assume that the bistability (“hot”–“cold”) and the associated hysteresis loop could be a potential reason for the mechanical failure (i.e., delamination fracture; Dewitte et al. [6], Supancic [7], Danzer et al. [8]) of PTC thermistors. In practice, however, the conditions prevailing could worsen either due to the dense packing of the electronic devices or because of poor cooling. This inevitably will shift the “hot” temperature even higher.

3.1. Current Control 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 shown in Figure 4, where the edge temperature ${\Theta }_e=\Theta (0)$ is plotted against $j^2$. For lower values of the conduction–convection parameter $(u=0.5)$, the three uniform solutions of Equation (24) are recovered (Figure 4a). As u increases, the number of solutions increases as well. Five solutions exist for $u=2$ in Figure 4b, seven solutions for $u=3$ in Figure 4c, and eleven for $u=5$ in Figure 4d. The corresponding temperature profiles are shown in Figure 5a–d. In every case, the solution structure consists 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 the boundary conditions in Equation (20) have two salient features. The first one is that as long as the nonuniform solutions are unstable, only the uniform ones are physically encountered, with the “cold” one being far more preferable from the operating point of view, since it results in the reduced thermal stresses and thermal loading (reduced fracture probability) of the device, as already explained in the previous paragraph. The second feature is a direct consequence of the first one, since the stable temperature profiles may be obtained from the solution of the much simpler algebraic problem in Equation (24) instead of solving the complete boundary value problem of Equation (17). Hence, the results concerning the bistability and the magnitude of the “hot” temperature previously obtained from the lumped (zero-dimensional) model are verified from the one-dimensional model as well. The hypothesis put forward of the induced thermal shock during the changeover from the “cold” to the “hot” state is further supported from the similar behavior encountered in other bistable systems, as in, for example, superconductors, when, for instance, a current lead feeding a superconductor magnet made from a high temperature superconductor switches from the superconducting (“cold”) state to the normal one (“hot”), say because of a loss of coolant flow accident, and then the temperature level could be in excess of 3000 K (see Krikkis [44] for a numerical calculation and Dresner [49] in paragraphs 10.3 to 10.8 for an analytical one under certain simplifying assumptions). The reason for such a similarity is the form of the electric resistivity, where, within a short temperature range, the resistivity changes by several orders of magnitude. The same bistability and hysteresis loops occur in boiling systems through the nonlinear and nonmonotonic boiling curve describing the three boiling modes, the stable nucleate and film, and the unstable transition one. When the wall temperature exceeds the critical heat flux, the system experiences a jump-like behavior, as shown in Figure 3, from the nucleate (“cold”) regime to the film (“hot”) regime and vice versa. The latter entails a significantly higher temperature (Zhukov and Barelko [50], Krikkis [51]). This nonlinear behavior is well understood and documented both theoretically and experimentally.

3.2. Problem P2

For the case where the edge temperature is fixed, the bifurcation analysis reveals a rich and interesting pattern, as shown in Figure 6. 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\dots )$ in Figure 6a, where the regions with unique, three, five, and seven (and so on) solutions are designated accordingly. 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 7. Again, the CCP has a profound effect on the solution structure, since the number of steady states calculated depends on its magnitude. Indeed, starting from $u=0.5$ and a unique solution in Figure 6b, three solutions emerge for $u=1$ in Figure 6c, five for $u=1.5$ in Figure 6d, and seven in Figure 6e 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 8a, the center is maintained at a lower temperature, while the temperature increases towards the edge $({\Theta }_c<{\Theta }_e)$. When the CCP increases and three solutions are present, as in Figure 8b, 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 8c for $u=1.5$ and seven in Figure 8d when $u=2$. The stable ones remain as 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 additional solutions emerge. Furthermore, as u increases, a temperature plateau around the center is formed and the temperature gradient along the device, although it cannot be entirely eliminated as was the case in problem P1, 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 [52] was the first to publish imbedded cusps in the study of first-order reactions in a 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. [53]). Another example of similar multiplicity behavior is encountered when the three boiling modes (nucleate, transition, and film) are excited by a uniform heat-generating source along a non-isothermal extended surface immersed in a boiling liquid (Krikkis [51]).

3.3. Voltage Control Problems P1 and P2

For the range of parameters considered when the voltage is the controlling parameter, the nonlocal problem in Equation (18) admits up to three solutions, as shown from the multiplicity pattern on the $(v,{\Theta }_e)$ plane in Figure 9 and from current–voltage characteristic in Figure 10 for problem P1. Two singular points exist and the temperature remains bounded for both problems P1 and P2. Clearly, this a key feature of PTC resistivity. In comparison, the multiplicity pattern for an NTC resistivity characteristic is simpler as it is restricted to a single limit point, but it is associated with an unbounded temperature (thermal runaway) once the voltage exceeds this limit point, leading to the so-called flash sintering phenomenon. Moreover, from a practical point of view, Figure 10 demonstrates the device potential for use as a current limiter, since the current density decreases with the applied voltage both in the regions with a unique solution and within the multiplicity region. When the heat transfer from the device edges (problem P2) is taken into account, again, up to three solutions exist, but the solution structure is qualitatively different. A bifurcation diagram on the $(j,{\Theta }_e)$ plane with the applied voltage v as a parameter is depicted in Figure 11. Along a line of constant voltage, the edge temperature shows only a weak dependence on the applied current, and in effect, the temperature is mostly affected by the voltage. Projecting the limit points of Figure 11 onto the $(v,j)$, plane we obtain the multiplicity region of the voltage–current characteristic shown in Figure 12. For this case, the multiplicity domain is confined between two cusp points, C₁ and C₂, respectively. For a better understanding of the solution structure, a geometrical representation of the associated cusp point is shown in Figure 13.

4. Conclusions and Future Work

An electrothermal model for a barium titanate-based PTC thermistor has been set up, featuring a nonlinear heat rejection mechanism consisting of natural convection combined with radiation through the device’s lateral surface. A smooth and strongly temperature-dependent electric conductivity function based on experimental data has been adopted. The temperature field has been resolved for both cases, where either the current or the voltage (nonlocal problem) is the controlling parameter. With the aid of an efficient continuation algorithm, multiple steady-state solutions that do not depend on the external circuit have been identified as a result of the inherent nonlinearities. Important findings are the profound effect of the conduction–convection parameter and of the boundary conditions on the multiplicity structure.

For the current control case, regardless of the type of boundary conditions imposed, the number of 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”. A key feature of this bistability is that the “cold” solution remains constant and close to the ambient temperature for a wide range of applied currents, but on the other hand, the “hot” steady state is proportional to the input power and the temperature is substantially higher. Therefore, a change in the operating point from the “cold” to the “hot” state closely resembles thermal shock, and could perhaps be a reason for the mechanical failure (delamination fracture) of PTC thermistors. Therefore, the control and effective reduction of the heat transmitted through the edges results in the reduced thermal loading of the device both in the steady-state operation and during transients, since the temperature gradient across the device vanishes. On the other hand, when the voltage is controlled, the multiple solutions are reduced to three, while the maximum temperature appears lower compared with the current controlled case.

As for future work, it would be interesting to extend the analysis to 2D/3D geometries and include alternating voltage and current as the controlling parameters. Another research direction could be the potential multiplicity features of ceramic devices with more complicated electrical resistivity functions resembling a combination of positive–negative temperature coefficients, as, for instance, those reported by Takeda et al. [38] and Fang et al. [54]. This particular form of the resistivity admits self-sustained oscillations (Hopf bifurcation), introducing a challenging dynamic behavior, as shown by Elmer [55].