Electrothermal Instabilities in Barium-Titanate-Based Ceramics

: An electrothermal analysis for barium-titanate-based ceramics is presented, combining the Heywang–Jonker model for the electric resistivity with a heat dissipation mechanism based on natural convection and radiation in a one-dimensional model on the device level with voltage as the control parameter. Both positive-temperature-coefficient (PTC) and negative temperature coefficient (NTC) effects are accounted for through the double Schottky barriers at the grain boundaries of the material. The problem formulated in this way admits uniform and non-uniform multiple-steady-state solutions that do not depend on the external circuit. The numerical bifurcation analysis reveals that the PTC effect gives rise to several multiplicites above the Curie point, whereas the NTC effect is responsible for the thermal runaway (temperature blowup). The thermal runaway phenomenon as a potential thermal shock could be among the possible reasons for the observed thermomechanical failures (delamination fracture). The theoretical results for the NTC regime and the thermal runaway are in agreement with the experimental flash sintering results obtained for barium titanate, and 3% and 8% yttria-stabilized zirconia.


Introduction
Since the first invention, back in 1955 [1], of the basic materials and processing technologies of positive-temperature-coefficient (PTC) ceramics, the field of interfacially controlled electroceramics is now of paramount technological importance, as well as high scientific interest.Two different categories of materials may be identified: (i) PTC resistors (thermistors) based on barium titanate, and (ii) ZnO varistors.A salient feature for both types is the highly non-linear voltage-current relationship caused by the formation of double Schottky barriers at the grain boundaries which are affected by the local composition of grain boundary regions' (doping and segregation) oxygen partial pressure and the applied voltage.We will discuss only the former, albeit certain similarities regarding the existence of multiple solutions and the thermal runaway phenomenon will be addressed in the appropriate sections.PTC thermistors are manufactured from silicon, barium, lead, and strontium titanates with the addition of yttrium, manganese, tantalium, and silica.They are extensively used for electric circuit protection, sensing excessive currents, and constant-temperature heating elements, and as current limiting devices, that is, as nondestructible (resettable) fuses.The technical and theoretical aspects may be found in the review papers [2,3].Although PTCs and, in general, electroceramics are, in principle, loaded electrically as components of electrical and/or electronic circuits, a significant number of mechanical failures is being recorded annually.Although the reasons for the failure are not fully understood, several studies are based on the Joule self-heating effect which causes temperature differences, thermal strains, and excessive thermo-mechanical stresses that may cause the failure of the device [4][5][6][7].Today, there is an increasing demand for a thorough analysis and understanding of the underlying non-linearly coupled electrothermal phenomena in electroceramic devices, since the current technological trends point towards designs with a reduced size and increasingly higher power densities.
The thermistor as a strongly non-linear coupled electrothermal problem (i.e., the exponential increase of the electric resistivity above the Curie temperature) has attracted significant engineering (theoretical and experimental [5][6][7][8][9][10][11][12][13]), and mathematical (applied and numerical [14-19]) attention.The above literature review suggests that the thermistor problem has been studied with various assumptions and/or restrictions.These usually include a constant heat transfer coefficient, an idealization of the PTC effect through a simplified electric conductivity function and, in certain cases, the influence of an external linear electric circuit.The latter admits up to three solutions as determined from the number of the intersection points between the current-voltage characteristic curves of the thermistor and the external circuit [15][16][17].
Yet, a complete analysis of the inherent electrothermal instabilities on the device level, induced by the combined positive and negative temperature coefficient (NTC) behavior, called thereafter the P-NTC effect, is still missing.The aim of the present study is to provide insight and, hopefully, a reasonable explanation to the most common reason of thermistor failure, namely, the delamination fracture due to excessive thermal loading (shock).To this end, the Heywang-Jonker model for the electric resistivity of the device is adopted, which describes both the positive and the negative temperature coefficient regimes.This non-linear and non-monotonic electric resistivity function of the temperature is combined with a non-linear temperature-dependent natural convection and radiation heat dissipation mechanism to form a one-dimensional distributed device model based on voltage control.The problem formulated in this way admits multiple-steady-state solutions that can be either uniform or non-uniform and do not depend on the external circuit.The numerical bifurcation analysis reveals that the PTC effect gives rise to several multiplicites above the Curie point, whereas the NTC effect is responsible for the temperature blowup (thermal runaway) which, unless detected and prevented, will lead to the destruction of the device.This result is further supported from a similar behavior that is encountered in other bistable systems such as superconductors [20] and boiling wires [21].

Energy Balance
Following the energy analysis in [13] for the device shown in Figure 1a, the temperature variation T along the longitudinal direction Z is describe by the following partial differential equation: This is the balance between the heat generated by the Joule effect and the heat dissipated by conduction, radiation, and natural convection.In the equation above, A is the cross-sectional area, T ∞ is the ambient temperature, C is the volumetric specific heat capacity, h c is the convective heat transfer coefficient, P is the perimeter, ε is the surface emissivity, σ SB is the Stefan-Boltzmann constant, J is the current density through the device, and E is the electric field intensity.For a constant DC current, the electric field intensity will be proportionally related to the current density through Ohm's law, E = ρ(T)J although the electric resistivity ρ will be a function of the local temperature [22].Since, in many practical applications, the controlling parameter is the applied voltage, it is instructive to introduce the electric potential Φ as [23]: Integrating the above relationship and considering that the current remains constant, albeit still unknown, yields: where V is the voltage drop across the device, as shown schematically in Figure 1c.Substituting the current density J from Equation (3) into Equation ( 1), the energy balance for voltage control takes the form: In contrast to the corresponding current control problem, this is a non-local problem since the solution depends on the resistivity integral over the device.Neumann-type boundary conditions are imposed on the device ends: implying that heat is dissipated only through the lateral surface of the conductor.

Heat Transfer Model
For the heat dissipation by natural convection, the correlation of Churchill and Chu [24] for the circumferentially average Nusselt number Nu is adopted: where f (Pr) is a weak function of the Prandtl number Pr. Equation ( 6) is valid for a very wide range of Rayleigh numbers, from 10 −6 to 10 9 , while it maintains a simple and compact mathematical form, which facilitates the calculation of the partial derivatives during the continuation on the stable and unstable branches of the solution.Equation ( 6) is applied locally in the evaluation of the convective heat flux along the device axis, in a similar manner as it was applied by Faghri and Sparrow for the calculation of the non-linear external heat flux in a horizontal pipe [25].Hence, the local Rayleigh number Ra is evaluated from the local temperature difference as: 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.

Electric Resistivity: The Heywang-Jonker Model
A salient 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 narrow temperature interval (PTC effect).After Heywang's [26] fundamental research, several authors [26][27][28][29] have explained the PTC effect in terms of a double-Schottky-barrier model.The acceptor states in the grain boundary core region acts as a trap for the majority charge carriers (electrons in the case of donor-doped BaTiO 3 ) which gives rise to a negative net charge of the grain boundary core and adjacent space charge layers in the grains where the majority charge carriers are depleted.Above the ferroelectric-paraelectric phase transition point (Curie temperature), the potential, ϕ, in the space charge layer can be obtained from Poisson's equation: where 0 is the permittivity of vacuum, is the relative permittivity, and e is the elementary charge, whereas n D and n A are the densities of the donor and acceptor states in the bulk region, respectively.The solution of Equation ( 8) with boundary conditions ϕ(0) = ϕ 0 and ϕ(w 0 ) = 0 yields the depletion zone width w 0 and the potential barrier height ϕ 0 of the space charge layers with N s being the density of the occupied acceptor states in the boundary layer.The temperature dependence of the relative permittivity is given by the Curie-Weiss law where c = 1.5 × 10 5 K is the Curie-Weiss constant and T C corresponds to the Curie-Weiss temperature.The density of the occupied acceptor states can be expressed as a function of temperature in terms of Fermi-Dirac statistics: with N s0 being the acceptor state density in the boundary layer and ε A the energy level of the acceptor states relative to the conduction band edge.The energy level ε F corresponds to the Fermi level defined as with N c being the density of states of the conduction band, which corresponds to the density of titanium ions in the lattice (N c = 1.56 × 10 22 cm −3 ) [27,29], owing to a small band width of the order of 0.1 eV.In the linear response regime, Jonker showed that the total resistivity is given by [27,30]: where ρg and ζ denote the grain resistivity and the grain boundary density, respectively.Equation ( 14) describes the characteristic increase in the resistivity by several orders of magnitude above the ferroelectric-paraelectric transition temperature T C = 110-120 • C.However, after exceeding T max = 250-270 • C where the maximum resistivity is encoun-tered, the resistivity is entering the negative temperature coefficient (NTC) regime where its magnitude is rapidly decreasing with increasing temperature.It is exactly this combined non-linear and non-monotonic P-NTC effect that gives rise to multiple steady states and eventually to thermal runaway as it will be shown in the following paragraphs.The validity of Equation ( 14) well within the NTC regime and up to 900 • C has been experimentally verified [31].

The Electrothermal Model in Dimensionless Form
The electrothermal model developed in the previous paragraph may be recast in dimensionless form through the below variables' transformation: To simplify the convective and radiative terms a convenient choice for the reference temperature is T ref = T ∞ , whereas the magnitude of the maximum resistivity is taken as the reference value, ρref = ρmax .Introducing the new variables into Equation ( 4), the temperature distribution along the device takes the form: In the above equation, the conduction-convection parameter (CCP) is defined as: where k is the thermal conductivity and L is the device length.The reference heat transfer coefficient h ref is defined through the Nusselt number: In terms of the dimensionless variables defined above, the local Rayleigh number becomes: The current density parameter is related to the current density as below and C h is the ratio of the radiative heat transfer coefficient to the reference heat transfer coefficient Under steady-state conditions, the partial differential equation, Equation ( 16), reduces to a two-point boundary value problem for the device temperature Θ(z) where Θ z = dΘ/dz, Θ zz = d 2 Θ/dz 2 , and The dimensionless voltage-current relationship (the global constraint) becomes with boundary conditions

Stability
The linear stability of a certain steady state j s , Θ s (z) to small perturbations is deter- mined by substituting the below expansions into Equations ( 16) and ( 21) Retaining only linear terms, the current density perturbation reads: Assuming the below form of the temperature perturbation in terms of the eigenvalue λ and the eigenfunction ϑ(z) δΘ(z, τ) = e λτ ϑ(z), (26) we obtain an integro-differential problem for the eigenvalue λ where and The corresponding boundary conditions are ϑ z (0) = ϑ z (1) = 0. Stable solutions are characterized by negative eigenvalues, whereas positive ones correspond to unstable temperature distributions.It is worth pointing out that similar non-local eigenvalue problems appear, for instance, in the theory of the microwave heating of ceramic materials [32].

Numerical Methods
An efficient way to solve the second-order integro-differential two-point boundary value problem in Equation ( 20) is, first, to transform it into a standard form that can be handled by ordinary differential equation (ODE) solvers.Employing the new variables 20) is transformed into a system of first-order equations which can be integrated by standard solvers, say, from z = 0 to z = 1, since the voltagecurrent constraint is now an implicit boundary condition in j at z = 1, i.e., v − ujΘ 3 (j; z = 1) = 0. (31)

Results and Discussion
Before we analyze the complete numerical solution, it is instructive to discuss, first, the uniform solutions of Equation (20), which reduces to an algebraic one for a constant temperature profile: A geometrical (graphical) solution is depicted in Figure 2 where the heat generation Q g and the heat dissipation Q d curves are being plotted.Depending on the magnitude of the applied voltage v, up to two solutions of Equation ( 32) may be obtained from the number of the intersection points shown in Figure 2. The solutions are projected in the (v, Θ) plane in Figure 3 with u as a parameter.From the two solutions, the one corresponding to a lower temperature, i.e., Θ < Θ max is stable (solid line), whereas the other one is unstable (dashed line).The two branches, the stable and the unstable one, are connected through a singular point with a characteristic voltage magnitude.Any applied voltage that exceeds this threshold will lead to an instability, namely, a thermal runaway (temperature blowup).Hence, the curve connecting the singular points forms the instability threshold which depends only on the conduction-convection parameter u since the ambient temperature and the donor and acceptor densities remain fixed.The associated instabilities may also be understood from the voltage-current relationship shown in Figure 4 which resembles an unusual lobe-like curve.Such non-linear and nonmonotonic v-j characteristic curves are encountered, for instance, in semiconductors (S-, N-, or Z-shaped) [33], superconductors [34], and organic LEDs [35], leading to interesting and complex bifurcations.

Results and Discussion
Before we analyze the complete numerical solution, it is instructive to discuss, first, the uniform solutions of Equation (20), which reduces to an algebraic one for a constant temperature profile: A geometrical (graphical) solution is depicted in Figure 2 where the heat generation  and the heat dissipation  curves are being plotted.Depending on the magnitude of the applied voltage , up to two solutions of Equation ( 32) may be obtained from the number of the intersection points shown in Figure 2. The solutions are projected in the (, Θ) plane in Figure 3 with u as a parameter.From the two solutions, the one corresponding to a lower temperature, i.e., Θ < Θ is stable (solid line), whereas the other one is unstable (dashed line).The two branches, the stable and the unstable one, are connected through a singular point with a characteristic voltage magnitude.Any applied voltage that exceeds this threshold will lead to an instability, namely, a thermal runaway (temperature blowup).Hence, the curve connecting the singular points forms the instability threshold which depends only on the conduction-convection parameter u since the ambient temperature and the donor and acceptor densities remain fixed.The associated instabilities may also be understood from the voltage-current relationship shown in Figure 4 which resembles an unusual lobe-like curve.Such non-linear and non-monotonic vj characteristic curves are encountered, for instance, in semiconductors (S-, N-, or Zshaped) [33], superconductors [34], and organic LEDs [35], leading to interesting and complex bifurcations.
Energy balance and uniform solutions of Equation (32).Heat generation (  ) and dissipation curves ( ) are plotted against temperature.For a certain range of the applied voltage , two intersection points exist, one below and one above Θ , corresponding to one stable and one unstable solution, respectively.Above a certain threshold (singular point at  = 2.672), i.e., for  = 3, no intersection point exists and a thermal runaway is triggered.(32).Heat generation (Q g ) and dissipation curves (Q d ) are plotted against temperature.For a certain range of the applied voltage v, two intersection points exist, one below and one above Θ max , corresponding to one stable and one unstable solution, respectively.Above a certain threshold (singular point at v = 2.672), i.e., for v = 3, no intersection point exists and a thermal runaway is triggered.32) projected on the (v, j) plane.Increasing the voltage beyond the singular points, the current blows up, resulting in temperature blowup.
Interestingly, when the conduction term Θ zz is taken into consideration, a far more complicated solution structure and multiplicity pattern emerges, as shown in Figure 5. There, the base temperature Θ b = Θ(0) has been selected as the bifurcation parameter and it is plotted against the dimensionless applied voltage v. Compared with Figure 3, both uniform solutions are recovered, and, in addition, several non-uniform solutions appear mostly within the zone Θ C < Θ < Θ max .The temperature profiles that correspond to the multiplicity pattern of Figure 5 are shown in Figure 6 for v = 0.1.Our linear stability analysis presented in paragraph II.5 reveals that the spatially periodic profiles (specifically 2 and 4 in Figure 6) are unstable, while profiles 1 and 5, which are antisymmetric with respect to the center of the device, are stable.This is in agreement with the stability analysis results of Elmer [36] for an idealized resistivity-temperature curve consisting of a step function below the transition point and a linearly decreasing segment above the transition point.This means that, in contrast to the uniform case where a single path for the instability existed, now, the path to instability may come from the stable non-uniform solution as indicated by the red arrows in Figure 6.An important conclusion that can drawn from Figure 6 is that the abrupt change in the resistivity (PTC effect) is responsible for the multiplicity below approximately the temperature of maximum resistivity while the NTC effect leads to a thermal runaway since the heat generation rate exceeds the heat that can be dissipated by natural convection and radiation.A similar thermal runaway phenomenon is observed in metallic conductors and high-temperature superconductors as a result of the non-linear relationship between the temperature and the electric resistivity [21,37].
Interestingly, when the conduction term Θ is taken into consideration, a far more complicated solution structure and multiplicity pattern emerges, as shown in Figure 5. There, the base temperature Θ = Θ(0) has been selected as the bifurcation parameter and it is plotted against the dimensionless applied voltage v. Compared with Figure 3, both uniform solutions are recovered, and, in addition, several non-uniform solutions appear mostly within the zone Θ < Θ < Θ .The temperature profiles that correspond to the multiplicity pattern of Figure 5 are shown in Figure 6 for  = 0.1.Our linear stability analysis presented in paragraph II.5 reveals that the spatially periodic profiles (specifically 2 and 4 in Figure 6) are unstable, while profiles 1 and 5, which are antisymmetric with respect to the center of the device, are stable.This is in agreement with the stability analysis results of Elmer [36] for an idealized resistivity-temperature curve consisting of a step function below the transition point and a linearly decreasing segment above the transition point.This means that, in contrast to the uniform case where a single path for the instability existed, now, the path to instability may come from the stable non-uniform solution as indicated by the red arrows in Figure 6.An important conclusion that can drawn from Figure 6 is that the abrupt change in the resistivity (PTC effect) is responsible for the multiplicity below approximately the temperature of maximum resistivity while the NTC effect leads to a thermal runaway since the heat generation rate exceeds the heat that can be dissipated by natural convection and radiation.A similar thermal runaway phenomenon is observed in metallic conductors and high-temperature superconductors as a result of the non-linear relationship between the temperature and the electric resistivity [21,37].20) with bifurcation parameter Θ  = Θ(0).Apart from the uniform solutions that are recovered, as in Figure 3, stable (continuous) and unstable (dashed) non-uniform solutions exist as well.The instability mechanism is the same as before, but, now, the path to thermal runway could also be through the stable non-uniform branches as indicated by the red arrows.20) with bifurcation parameter Θ b = Θ(0).Apart from the uniform solutions that are recovered, as in Figure 3, stable (continuous) and unstable (dashed) non-uniform solutions exist as well.The instability mechanism is the same as before, but, now, the path to thermal runway could also be through the stable non-uniform branches as indicated by the red arrows.20).Solutions designated as 1, 2, 4, and 5 correspond to non-uniform profiles, whereas 3 and 6 correspond to uniform profiles.

The Relationship between Temperature and Mechanical Failure
In general, the mechanical failures of electroceramic components are attributed to thermal stresses induced by the temperature gradients developed.Dewitte et al. [5] carried out a thermo-elastic analysis which verified that a temperature difference between the core and surface exists during transient switching processes and the amplitudes of the resulting thermal stresses are sensitive to the applied boundary conditions on the edges of the device.However, the calculated mechanical stresses in a homogeneous PTC component are too low to explain the delamination fracture.As a possible explanation, it has been postulated that excessive stress amplitudes may be a consequence of inhomogeneities in the resistance field within the ceramic material.Supancic [6] extended the resistivity model to include the varistor effect, that is, the resistivity dependence on the electric field intensity, also known as the varistor effect (voltage-dependent resistor).The analysis and the subsequent experimental verification showed that the varistor effect significantly changes the thermo-electrical response of the device since higher temperature gradients are prevailing, which give rise to substantial thermo-mechanical stresses.In any case, with or without the varistor effect, the calculated temperature remains bounded [5,6].This phenomenon may be explained from a different perspective, as soon as it is recognized that, because of the PTC effect, the problem is a bistable one; that is, it admits multiple solutions (states), one "cold" and one "hot" [13,15].The latter state is inevitably associated with temperatures of sufficient magnitude either under voltage or under current control.However, the present study shows that, when the whole P-NTC effect is taken into consideration, on one hand, the stable non-uniform steady state may exceed the temperature of maximum resistivity at high voltages as shown in Figure 5.This temperature profile is not symmetric with respect to the center and a significant temperature gradient may develop across the device.On the other hand, an electrothermal instability, which is a  20).Solutions designated as 1, 2, 4, and 5 correspond to non-uniform profiles, whereas 3 and 6 correspond to uniform profiles.

The Relationship between Temperature and Mechanical Failure
In general, the mechanical failures of electroceramic components are attributed to thermal stresses induced by the temperature gradients developed.Dewitte et al. [5] carried out a thermo-elastic analysis which verified that a temperature difference between the core and surface exists during transient switching processes and the amplitudes of the resulting thermal stresses are sensitive to the applied boundary conditions on the edges of the device.However, the calculated mechanical stresses in a homogeneous PTC component are too low to explain the delamination fracture.As a possible explanation, it has been postulated that excessive stress amplitudes may be a consequence of inhomogeneities in the resistance field within the ceramic material.Supancic [6] extended the resistivity model to include the varistor effect, that is, the resistivity dependence on the electric field intensity, also known as the varistor effect (voltage-dependent resistor).The analysis and the subsequent experimental verification showed that the varistor effect significantly changes the thermoelectrical response of the device since higher temperature gradients are prevailing, which give rise to substantial thermo-mechanical stresses.In any case, with or without the varistor effect, the calculated temperature remains bounded [5,6].This phenomenon may be explained from a different perspective, as soon as it is recognized that, because of the PTC effect, the problem is a bistable one; that is, it admits multiple solutions (states), one "cold" and one "hot" [13,15].The latter state is inevitably associated with temperatures of sufficient magnitude either under voltage or under current control.However, the present study shows that, when the whole P-NTC effect is taken into consideration, on one hand, the stable non-uniform steady state may exceed the temperature of maximum resistivity at high voltages as shown in Figure 5.This temperature profile is not symmetric with respect to the center and a significant temperature gradient may develop across the device.On the other hand, an electrothermal instability, which is a consequence of the NTC effect at higher temperatures, may trigger a temperature blowup when the runaway voltage is exceeded, which could be another reason for the thermo-mechanical failures observed in the electroceramic devices.
The preceding multiplicity analysis and the instability leading to thermal runaway is also applicable to metal oxide (i.e., ZnO) varistors.The ZnO-based varistor is a highly non-linear two-terminal polycrystalline device, commonly known as multi-component metal oxide varistors (MOVs) [38][39][40].Because of their highly symmetric non-linear currentvoltage characteristic curve with respect to the polarity of the applied voltage, these MOVs demonstrate an excellent surge-withstanding capability.Similarly to the barium-based ceramics, these electric characteristics can be explained with the formation of potential barriers (Schottky barriers) involving thin insulating layers around the successive ZnO grains [41][42][43].Because of their special characteristics, varistiors are widely used in the power systems' protection, for suppressing internally generated spikes in electronic circuits, as surge absorbers, and as surge protectors of electric and electronic circuitry, among others.For the surge arrester application in particular, the energy absorption capability while maintaining thermal stability is of paramount importance.Three failure modes have been identified thus far: (i) electrical puncture, (ii) physical cracking, and (iii) thermal runaway [39,44].Similar energy balances between heat generation and heat dissipation leading to thermal runaway have been recently reported [38,[45][46][47].

Comparison with Experiments
Above Θ max in the NTC regime where the potential barrier assumes an Arrheniustype temperature dependence, as in Equation ( 33), the thermistor problem, and, especially, the thermal runaway phenomenon described above, is closely associated with the flash sintering of ceramics as in, for instance, barium and strontium titanates, and 3% and 8% yttria-stabilized zirconia, just to name a few.The key feature of the process is to exercise control when an operating parameter, usually the furnace temperature, exceeds the corresponding limit point, established by the voltage applied to the specimen.Under these conditions, the heat dissipation mechanism can no longer balance the Joule heating and the temperature blows up.To maintain the temperature below a certain threshold, the process controller will switch from voltage to current control [48,49].Detailed reviews may be found in the papers [50][51][52].In this case, the electric resistivity is given by a simpler and monotonic Arrhenius relationship: where ε a is the activation energy, R is the universal gas constant, and ρ0 the pre-exponential factor.Now, T ∞ (Θ ∞ ) stands for the furnace temperature, which is no longer a constant but rather a parameter.Thus, a different reference temperature has been selected, namely, T ref = 1000 K.A convenient choice for the reference electric resistivity is: With the above reference values, the dimensionless electric resistivity yields: Utilizing the Frank-Kamenetskii approximation for the exponential term of the electric resistivity, Hewitt et al. [53] obtained closed-form solutions of Equation (20) for plane and radial geometries under certain assumptions.In the present study, the bifurcation and stability analysis are carried out numerically; therefore, no simplified assumptions for the form of the electric resistivity have been introduced, so the same model described by Equation (30) and Neumann boundary conditions (insulated edges) has been used for the numerical solution of the flash sintering phenomenon.The solution structure and the singular points are presented in Figure 7, which looks very similar to Figure 3, which describes the corresponding uniform solutions in the whole P-NTC regime.Hence, it is worth pointing out that the boundary value problem in Equation ( 20), together with the monotonically decreasing function of the electric resistivity in Equation (35) and Neumann boundary conditions, Θ z (0) = Θ z (1) = 0, admits only uniform solutions, i.e., Θ z = Θ zz = 0.In fact, under these circumstances the solution can be obtained by solving the much simpler zero-dimensional model described by the algebraic relationship in Equation (32).This kind of boundary conditions, i.e., insulated edges, may be justified by the presence of the connecting electrodes (see Figure 1c) which diminish cooling through the edges acting like insulators.As long as only uniform solutions exist, there is a single singular point that connects the stable and the unstable branches.As soon as the voltage or the furnace temperature limits specified by the singular point, a limit point in this case, are exceeded, a thermal runaway is under way and the process controller must switch from voltage to current control in order to maintain the temperature within limits.This is because, for the current control mode, a unique solution exists for an NTC device.Similar instability criteria have been proposed for the uniform temperature model [54].The singular points calculated from Equation (32) are in good agreement with the experimental data for 3YSZ [49] shown in Figure 8, for 8YSZ [48,55] shown in Figure 9, and for barium titanate [56] shown in Figure 10.The absence of non-uniform solutions when the edges of the specimen are insulated or, in general, when the heat transfer is poor could be a reason why a single energy balance between Joule heating and cooling through radiation and natural convection in Equation ( 32) is quite successful and, consequently, has been widely used for correlating experimental flash sintering data.
worth pointing out that the boundary value problem in Equation ( 20), together with the monotonically decreasing function of the electric resistivity in Equation (35) and Neumann boundary conditions, Θ (0) = Θ (1) = 0, admits only uniform solutions, i.e., Θ = Θ = 0 .In fact, under these circumstances the solution can be obtained by solving the much simpler zero-dimensional model described by the algebraic relationship in Equation (32).This kind of boundary conditions, i.e., insulated edges, may be justified by the presence of the connecting electrodes (see Figure 1c) which diminish cooling through the edges acting like insulators.As long as only uniform solutions exist, there is a single singular point that connects the stable and the unstable branches.As soon as the voltage or the furnace temperature limits specified by the singular point, a limit point in this case, are exceeded, a thermal runaway is under way and the process controller must switch from voltage to current control in order to maintain the temperature within limits.This is because, for the current control mode, a unique solution exists for an NTC device.Similar instability criteria have been proposed for the uniform temperature model [54].The singular points calculated from Equation (32) are in good agreement with the experimental data for 3YSZ [49] shown in Figure 8, for 8YSZ [48,55] shown in Figure 9, and for barium titanate [56] shown in Figure 10.The absence of non-uniform solutions when the edges of the specimen are insulated or, in general, when the heat transfer is poor could be a reason why a single energy balance between Joule heating and cooling through radiation and natural convection in Equation ( 32) is quite successful and, consequently, has been widely used for correlating experimental flash sintering data.

Conclusions
An electrothermal model for barium-titanate-based ceramics has been developed, combining the Heywang-Jonker model for the electric resistivity with a heat dissipation mechanism based on natural convection and radiation in a one-dimensional model on the device level with voltage as the control parameter.Both PTC and NTC effects are accounted for through the double Schottky barriers at the grain boundaries of the material.The problem formulated in this way admits uniform and non-uniform multiple-steadystate solutions that do not depend on the external circuit.The instability thresholds calculated for the NTC regime and the associated thermal runaway are in agreement with the experimental flash sintering results obtained for barium titanate, and 3% and 8% yttria-stabilized zirconia.The important findings are as follows:

•
The PTC effect gives rise to multiple solutions mainly in the temperature range between the Curie and the maximum resistivity points.

•
Thermal runaway is due to the NTC effect.The runaway voltage depends on the conduction-convection parameter u.

•
Thermal runaway as a thermal shock is a potential reason for the thermo-mechanical failures observed (delamination fracture).

•
For the NTC regime (flash sintering) when Neumann boundary conditions are imposed on the distributed model, only uniform solutions are admitted, one stable and one unstable.

Figure 1 .
Figure 1.(a) Conductor geometry; (b) energy balance on an elementary volume element; and (c) simplified electric circuit.

Figure 2 .
Figure 2. Energy balance and uniform solutions of Equation(32).Heat generation (Q g ) and dissipation curves (Q d ) are plotted against temperature.For a certain range of the applied voltage v, two intersection points exist, one below and one above Θ max , corresponding to one stable and one unstable solution, respectively.Above a certain threshold (singular point at v = 2.672), i.e., for v = 3, no intersection point exists and a thermal runaway is triggered.

Figure 3 .
Figure 3. Solution structure of Equation (32) projected on the (, Θ) plane and instability threshold.Stable solutions are plotted wih a continuous line, whereas the unstable ones are plotted with a dashed line.The thermal runaway triggered by an increase in the applied voltage is shown schematically with red arrows.The instability threshold is formed by the singular points (•) and depends on the conduction-convection parameter .

Figure 4 .
Figure 4. Solution structure of Equation (32) projected on the (, ) plane.Increasing the voltage beyond the singular points, the current blows up, resulting in temperature blowup.

Figure 3 .
Figure 3. Solution structure of Equation (32) projected on the (v, Θ) plane and instability threshold.Stable solutions are plotted wih a continuous line, whereas the unstable ones are plotted with a dashed line.The thermal runaway triggered by an increase in the applied voltage is shown schematically with red arrows.The instability threshold is formed by the singular points (•) and depends on the conduction-convection parameter u.

Figure 3 .
Figure 3. Solution structure of Equation (32) projected on the (, Θ) plane and instability thres Stable solutions are plotted wih a continuous line, whereas the unstable ones are plotted w dashed line.The thermal runaway triggered by an increase in the applied voltage is sh schematically with red arrows.The instability threshold is formed by the singular points (• depends on the conduction-convection parameter .

Figure 4 .
Figure 4. Solution structure of Equation (32) projected on the (, ) plane.Increasing the vo beyond the singular points, the current blows up, resulting in temperature blowup.

Figure 4 .
Figure 4. Solution structure of Equation(32) projected on the (v, j) plane.Increasing the voltage beyond the singular points, the current blows up, resulting in temperature blowup.

Figure 5 .
Figure 5. Solution structure for the distributed model in Equation (20) with bifurcation parameter Θ  = Θ(0).Apart from the uniform solutions that are recovered, as in Figure3, stable (continuous) and unstable (dashed) non-uniform solutions exist as well.The instability mechanism is the same as before, but, now, the path to thermal runway could also be through the stable non-uniform branches as indicated by the red arrows.

Figure 5 .
Figure 5. Solution structure for the distributed model in Equation (20) with bifurcation parameter Θ b = Θ(0).Apart from the uniform solutions that are recovered, as in Figure3, stable (continuous) and unstable (dashed) non-uniform solutions exist as well.The instability mechanism is the same as before, but, now, the path to thermal runway could also be through the stable non-uniform branches as indicated by the red arrows.

Figure 6 .
Figure 6.Numerically computed temperature profiles for the distributed model of Equation(20).Solutions designated as 1, 2, 4, and 5 correspond to non-uniform profiles, whereas 3 and 6 correspond to uniform profiles.

Figure 6 .
Figure 6.Numerically computed temperature profiles for the distributed model of Equation(20).Solutions designated as 1, 2, 4, and 5 correspond to non-uniform profiles, whereas 3 and 6 correspond to uniform profiles.

Figure 7 .
Figure 7.Instability threshold formed by the singular points (•) for flash sintering.Θ now stands for the furnace temperature which is no longer constant.The electric resistivity is given by the Arrhenius form in Equation(35).When the NTC effect is combined with Neumann boundary conditions, Equation(20) and Equation(32) have the same uniform solutions.

Figure 7 .
Figure 7.Instability threshold formed by the singular points (•) for flash sintering.Θ ∞ now stands for the furnace temperature which is no longer constant.The electric resistivity is given by the Arrhenius form in Equation(35).When the NTC effect is combined with Neumann boundary conditions, Equations(20) and(32) have the same uniform solutions.