Prerpints.org logo
Preprint
Article

This version is not peer-reviewed.

Electrothermal Instabilities in Barium Titanate Based Ceramics

A peer-reviewed article of this preprint also exists.

Submitted:

07 March 2024

Posted:

08 March 2024

You are already at the latest version

Abstract
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 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 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 experimental flash sintering results obtained for 3% and 8% yttria stabilized zirconia.
Keywords: 
BaTiO3 and titanates
;  
thermistor
;  
Joule heating
;  
nonlocal problem
;  
bifurcation analysis
;  
delamination fracture
Subject: 
Physical Sciences  -   Applied Physics

I. 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 nonlinear 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 widely used as current limiting devices, that is as nondestructible (resettable) fuses for electric circuit protection, sensing excessive currents or constant temperature heating elements. Technical and theoretical aspects may be found in the review papers [2,3]. Although PTCs and in general electroceramics are in principle loaded electrically a significant number of mechanical failures is being 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 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 thorough analysis and understanding of the underlying coupled electrothermal phenomena in electroceramic devices [4,5,6,7].
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,15,16,17,18,19]) attention. From the literature above it appears that the thermistor problem has been studied with various assumptions and/or restrictions related primarily with the form of the electric conductivity (in essence an idealization of the PTC effect), the heat dissipation mechanism and in certain cases with the influence of the external electric circuit. The later is also associated with the existence of multiple steady state solutions, up to three, as determined from the number of the intersection points between the current-voltage characteristic curves of the external (linear) circuit and the thermistor [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 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 thermistors 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].

II. Analysis

II.1. Energy Balance

Consider a horizontal cylindrical segment of a conductor of uniform material density with constant thermal conductivity k. The specimen has diameter D and length L as it is schematically shown in Figure 1a. 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 ∞ . An energy balance along the longitudinal direction Z yields the following differential equation for the device temperature T:
C ∂ T ∂ t = ∂ ∂ Z ( k ∂ T ∂ Z ) − P A [ h c ( T − T ∞ ) + ϵ σ S B ( T 4 − T ∞ 4 ) ] + E J
In the equation above, C is the volumetric 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, σ S B 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 = ρ ^ ( T ) J , ρ ^ being the electric resistivity [22]. 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 Φ as [23]
E = − d Φ d Z = ρ ^ J
Integrating the above relationship and considering that the current remains constant, albeit still unknown, yields:
J ∫   0   L ρ ^ ( T ) d Z = − [ Φ ( L ) − Φ ( 0 ) ] = V ,
where V is the voltage drop across the device, as shown schematically in Figure 1c. Substituting the current density J from Eq. (3) into Eq.(1) the energy balance for voltage control takes the form:
C ∂ T ∂ t = ∂ ∂ Z ( k ∂ T ∂ Z ) − P A [ h c ( T − T ∞ ) + ϵ σ S B ( T 4 − T ∞ 4 ) ] + ρ ^ V 2 [ ∫   0   L ρ ^ ( T ) d Z ] 2 .
In contrast to the corresponding current control problem, this is a nonlocal problem since the solution depends on the resistivity integral over the device. Neumann type boundary conditions are imposed on the device ends:
∂ T ∂ Z | Z = 0 = ∂ T ∂ Z | Z = L = 0 ,
implying that heat is dissipated only through the lateral surface of the conductor.

II.2. Heat Transfer Model

The heat generated in the device due to the current flow (Joule heating) 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 [24] is employed:
Nu = 0.36 + 0.518 Ra 1 / 4 f ( Pr ) ,   f ( Pr ) = [ 1 + ( 0.559 Pr ) 9 / 16 ] 4 / 9 ,
where f ( Pr ) is a weak function of the Prandtl number Pr. Equation (6) covers a very wide range of Rayleigh numbers namely in the range from 10-6 to 109, while it maintains a simple and compact mathematical form. Equation (6) 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 [25]. Consequently the local Rayleigh number Ra is evaluated from the local temperature difference as:
Ra = g β D 3 [ T ( Z ) − T ∞ ] / α ν ,
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.

II.3. Electric Resistivity. The Heywang–Jonker Model

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 (PTC effect). After Heywang’s [26] fundamental research, the PTC effect has been explained in terms of a double-Schottky-barrier model by numerous authors [26,27,28,29]. The majority charge carriers (electrons in the case of donor–doped BaTiO3) are trapped by acceptor states in the grain boundary core region, 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 Curie temperature (ferroelectric–paraelectric phase transition point), the potential, φ , in the space charge layer can be calculated from Poisson’s equation
d 2 φ d x 2 = − e ( n D − n A ) ε 0 ε ,
where ε 0 , is the permittivity of vacuum, ε is the relative permittivity, e is the elementary charge, whereas n D , and n A are the and densities of donor and acceptor states in the bulk region, respectively. The solution of Eq.(8) with boundary conditions φ ( 0 ) = φ 0 and φ ( w 0 ) = 0 yields the depletion zone width w 0
w 0 = N s 2 ( n D − n A ) ,
and the potential barrier height φ 0 of the space charge layers
φ 0 = e N s 2 8 ε 0 ε ( n D − n A ) ,
with N s being the density of occupied acceptor states in the boundary layer. The temperature dependence of the relative permittivity is given by the Curie–Weiss law
ε = c T − T C ,
where c = 1.5 × 10 5 K is the Curie–Weiss constant and T C corresponds to the Curie–Weiss temperature. The density of occupied acceptor states can be expressed as a function of temperature in terms of Fermi–Dirac statistics:
N s = N s 0 1 + exp [ ( ε φ 0 + ℰ F − ℰ A ) / k B T ] ,
with N s 0 being the acceptor state density in the boundary layer and ℰ A is the energy level of the acceptor states relative to the conduction band edge. The energy level ℰ F corresponds to the Fermi level defined as
ℰ F = k B T ln [ N c / ( n D − n A ) ] ,
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].
ρ ^ = ρ ^ g [ 1 + ζ w 0 ψ exp ( ψ ) ] ,   ψ = e φ 0 k B T ,
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 encountered, 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 Eq.(14) well within the NTC regime and up to 900°C has been experimentally verified [31].

II.4. 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:
z = Z / L ,   τ = α t / L 2 ,   Θ = T / T ref ,   ρ = ρ ^ / ρ ^ ref ,   v = V / V ref ,   V ref 2 = k ( T ρ ^ ) ref .
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 Eq.(4) the temperature distribution along the device takes the form:
∂ Θ ∂ τ = ∂ 2 Θ ∂ z 2 − u 2 [ Nu ( Θ − 1 ) + C h ( Θ 4 − 1 ) − ( v u ∫   0   1 ρ d z ) 2 ρ ] .
In the above equation, the conduction–convection parameter (CCP) is defined as:
u 2 = h ref L 2 k ( A / P ) ,
where the reference heat transfer coefficient h ref is defined through the Nusselt number:
Nu = h c k ∞ / D = h c h ref .
In terms of the dimensionless variables defined above the local Rayleigh number becomes:
Ra ( z ) = Ra ∞ [ Θ ( z ) − 1 ] ,   Ra ∞ = g β D 3 T ∞ α ν .
The current density parameter is related to the current density as below
j 2 = J 2 ( A / P ) T ∞ ( ρ ^ h ) ref ,
and C h is the ratio of the radiative heat transfer coefficient to the reference heat transfer coefficient
C h = ϵ σ S B T ∞ 3 h ref = h r h ref .
Under steady state conditions the partial differential equation Eq.(16) reduces to a two point boundary value problem for the device temperature Θ ( z )
Θ ″ − u 2 [ Nu ( Θ − 1 ) + C h ( Θ 4 − 1 ) − ( v u 〈 ρ 〉 ) 2 ρ ] = 0 ,
where the primes denote differentiation with respect to z, i.e Θ ′ = d Θ / d z and
〈 ρ 〉 = ∫   0   1 ρ ( Θ ) d z .
The dimensionless voltage-current relationship (the global constraint) becomes
v − u j 〈 ρ 〉 = 0 ,
with boundary conditions
Θ ′ ( z = 0 ) = Θ ′ ( z = 1 ) = 0 .

II.5. Stability

The linear stability of a certain steady state j s , Θ s ( z ) to small perturbations is determined by substituting the below expansions into Eqs.(16, 21)
Θ ( z , τ ) = Θ s ( z ) + δ   Θ ( z , τ ) ,
j ( τ ) = j s + δ ​ j ( τ ) .
Retaining only linear terms the current density perturbation reads:
δ ​ j = − u j s 2 v ∫   0   1 ρ Θ δ   Θ d z .
Assuming the below form of the temperature perturbation in terms of the eigenvalue λ and the eigenfunction ϑ ( z )
Δ Θ ( z , τ ) = e λ τ ϑ ( z ) ,
we get an integro-differential problem for the eigenvalue λ
λ ϑ = ϑ ″ − [ u 2 Δ Q Θ ϑ − 2 u ( j 3 ρ ) s v ∫   0   1 ρ Θ ϑ d z ] = 0 ,
where
Δ Q Θ = ∂ ∂ Θ [ Nu ( Θ − 1 ) + C h ( Θ 4 − 1 ) − j s 2 ρ ( Θ ) ] Θ = Θ s ,
and
ρ Θ = ∂ ρ ∂ Θ | Θ = Θ s ,   ρ s = ρ ( Θ s ( z ) ) .
The corresponding boundary conditions are ϑ ′ ( 0 ) = ϑ ′ ( 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 microwave heating of ceramic materials [32].

III. Numerical Methods

An efficient way to solve the second order integro–differential two–point boundary value problem in Eq.(20) is first to transform it into a standard form that can be handled by ordinary differential equations (ODE) solvers. Employing the new variables Θ 1 = Θ , Θ 2 = Θ ′ and Θ ′ 3 = ρ ( Θ 1 ) , Eq.(20) is transformed into a system of first order equations
Θ 1 ′ = Θ 2 Θ 2 ′ = u 2 [ Nu ( Θ 1 − 1 ) + C h ( Θ 1 4 − 1 ) − j 2 ρ ( Θ 1 ) ] Θ 3 ′ = ρ ( Θ 1 ) ,
which can be integrated by standard solvers say from z = 0 to z = 1 since the voltage–current constraint is now an implicit boundary condition in j at z = 1 i.e:
v − u j Θ 3 ( j ; z = 1 ) = 0 .

IV. Results and Discussion

Before we analyze the complete numerical solution, it is instructive to discuss first the uniform solutions of Eq.(20), which reduces to an algebraic one for a constant temperature profile:
Nu ( Θ − 1 ) + C h ( Θ 4 − 1 ) − ( v / u ) 2 ρ ( Θ ) = 0 .
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 Eq.(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 be also understood from the voltage–current relationship shown in Figure 4 which resembles an unusual lobe–like curve. Such non–linear and non–monotonic 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.
Interestingly, when the conduction term Θ ″ is taken into consideration, a far more complicated solution structure and multiplicity pattern emerges as it is 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. Comparing 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].

IV.1. 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 delamination fracture. As a possible explanation, it has been postulated that excessive stress amplitudes may be 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 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 later 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 is also applicable to metal–oxide (i.e. ZnO) varistors. The ZnO based varistor is highly non–linear two terminal polycrystalline device commonly known as multi–component metal oxide varistors (MOVs) [38,39,40]. These MOVs exhibit excellent surge-withstanding capability via highly symmetric nonlinear current–voltage behavior with respect to the polarity of the applied voltage. This behavior is attributed to the electrical potential barriers (Schottky barriers) formed involving thin insulating layers around the successive ZnO grains [41,42,43]. These features allowed varistor applications in the power systems’ protection, suppression of internally generated spikes in electronic circuits, valve surge absorption and surge protection of electric and electronic circuitry among others. As a surge arrester 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].

IV.2. Comparison with Experiments

Above Θ max in the NTC regime where the potential barrier assumes an Arrhenius type temperature dependence, as in Eq.(33), the thermistor problem and especially the thermal runaway phenomenon described above is closely associated with the flash sintering of ceramics as for instance yttria stabilized zirconia, magnesia doped alumina and strontium titanate among others. The essence of the process is that when an operating parameter such as the furnance 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 blows up. The process controller is then switching from voltage to current control to maintain the temperature within the specified limits [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:
ρ ^ = ρ ^ 0 exp ( ℰ a / R T ) ,
where ℰ a is the activation energy, R is the universal gas constant, and ρ ^ 0 the pre–exponential factor. Now T ∞ ( Θ ∞ ) stands for the furnance temperature which is no longer 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:
ρ ^ ref = ρ ^ 0 exp ( ℰ a / R T ref ) = ρ ^ 0 exp ( γ ) .
with the above reference values the dimensionless electric resistivity yields:
ρ = exp [ 1 γ ( 1 Θ − 1 ) ] .
Utilizing the Frank–Kamenetskii approximation for the exponential term of the electric resistivity Hewitt et al. [53] obtained closed form solutions of Eq.(20) for plane and radial geometries under certain assumptions. In the present study the bifurcation and stability analysis is carried out numerically, therefore no simplified assumptions for the form of the electric resistivity have been introduced so the same model described by Eqs.(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 Eq.(20) together with the monotonically decreasing function of the electric resistivity in Eq.(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 Eq.(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 furnance 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 a NTC device. Similar instability criteria have been proposed for the uniform temperature model [55]. The singular points calculated from Eq.(32) are in good agreement with the experimental data for 3YSZ [49], shown in Figure 8 and for 8YSZ [48,54] shown in Figure 9. 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 Eq.(32) is quite successful and consequently has been widely used for correlating experimental flash sintering data.

V. 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 steady state 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 experimental flash sintering results obtained for 3% and 8% yttria stabilized zirconia. Important findings are:
  • 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 to the distributed model, only uniform solutions are admitted one stable and one unstable.

References

  1. Haayman, P.; Dam, R.; Klasens, H. Verfahren zur Herstellung halbleitenden Materials. German Patent 929350.
  2. Huybrechts, B.; Ishizaki, K.; Takata, M. The positive temperature coefficient of resistivity in barium titanate. J. Mater. Sci. 1995, 30, 2463–2474. [Google Scholar] [CrossRef]
  3. Preis, W.; Sitte, W. Electrical properties of grain boundaries in interfacially controlled functional ceramics. J. Electroceram. 2015, 34, 185–206. [Google Scholar] [CrossRef]
  4. Danzer, R.; Kaufmann, B.; Supancic, P. Failure of high power varistor ceramic components. J. Eur. Ceram Soc. 2020, 40, 3766–3770. [Google Scholar] [CrossRef]
  5. Dewitte, C.; Elst, R.; Delannay, F. On the mechanism of delamination fracture of BaTiO3–based PTC thermistors. J. Eur. Ceram. Soc. 1994, 14, 481–492. [Google Scholar] [CrossRef]
  6. Supancic, P. Mechanical stability of BaTiO3–based PTC thermistors components: experimental investigation and theoretical modeling. J. Eur. Ceram. Soc. 2000, 20, 2009–2024. [Google Scholar] [CrossRef]
  7. Chlup, Z.; Drdlík, D.; Hadraba, H. Temperature effect on elastic and fracture behaviour of lead-free piezoceramic BaTiO3. J. Eur. Ceram. Soc. 2023, 43, 1509–1522. [Google Scholar] [CrossRef]
  8. Smith, D.; Ghayoub, N.; Charissou, I.; Bellon, O.; Abélard, P.; Edwards, A. Transient Thermal Gradients in Barium Titanate Positive Temperature Coefficient (PTC) Thermistors. J. Am. Ceram Soc. 1998, 81, 1789–1796. [Google Scholar] [CrossRef]
  9. Mader, G.; Meixner, H.; Kleinschmidt, P. Study of microscopic heat sources in semiconducting barium titanate ceramics. J. Appl. Phys. 1984, 56, 2832–2836. [Google Scholar] [CrossRef]
  10. Heinen, B.; Waser, R. Influence of the thickness and area of NiCr/Ag electrodes on the characteristics of BaTiO3–ceramic based positive–temperature–coefficient thermistor. J. Mat. Sci. 1998, 33, 4603–4608. [Google Scholar]
  11. Karpov, E. Bistability, autowaves and dissipative structures in semiconductor fibers with anomalous resistivity properties. Philosophical Magazine. 2012, 92, 1300–1316. [Google Scholar] [CrossRef]
  12. Platzer, A.; Supancic, P.; Danzer, R. Thermography and Simulation of Advanced Nonlinear Electroceramic Components. Adv. Engin. Mater. 2006, 8, 1134–1139. [Google Scholar] [CrossRef]
  13. Krikkis, R. Multiplicity Analysis of a Thermistor problem – A Possible Explanation of Delamination Fracture. J. 2023, 6, 517–535. [Google Scholar] [CrossRef]
  14. Cimatti, G. Remark on existence and uniqueness for the thermistor problem under mixed boundary conditions. Quart. Appl. Math. 1989, 47, 117–121. [Google Scholar] [CrossRef]
  15. Fowler, A.; Frigaard, I.; Howison, S. Temperature surges in current-limiting circuit devices. SIAM J. Appl. Math. 1992, 52, 998–1011. [Google Scholar] [CrossRef]
  16. Howison, S.; Rodrigues, J.; Shilor, M. Stationary solutions to the thermistor problem. J. Math. Anal. Appl. 1993, 174, 573–588. [Google Scholar] [CrossRef]
  17. Zhou, S.; Westbrook, D. Numerical solutions of the thermistor equations. J. Comput. Appl. Math.; 1997, 79, 101–118. [Google Scholar] [CrossRef]
  18. Golosnoy, I.; Sykulski, J. Numerical modeling of non-linear coupled thermo-electric problems. A comparative study. COMPEL. 2009, 28, 639–655. [Google Scholar]
  19. Liero, M.; Koprucki, T.; Fisher, A.; Scholz, R.; Glitzky, A. p–Laplace thermistor modeling of electrothermal feedback in organic semiconductor devices. Z. Angew. Math. Phys. 2009, 66, 639–655. [Google Scholar] [CrossRef]
  20. Romanovskii, V.; Watanabe, K. Operating modes of high- composite superconductors and thermal runaway conditions under current charging. Supercond. Sci. Technol. 2006, 19, 541–550. [Google Scholar] [CrossRef]
  21. Krikkis, R. On the Thermal Dynamics of Metallic and Superconducting Wires. Bifurcations, Quench, the Destruction of Bistability and Temperature Blowup. J 2021, 4, 803–823. [Google Scholar] [CrossRef]
  22. Metaxas, A. Foundations of Electroheat. A Unified Approach; Wiley, 1996. [Google Scholar]
  23. Lupi, S. Fundamentals of Electroheat. Electrical Technologies for Process Heating; Springer, 2017. [Google Scholar]
  24. Churchill, S.; Chu, H.H. Correlating equations for laminar and turbulent free convection from a horizontal cylinder. Int. J. Heat Mass Transfer. 1975, 18, 1049–1053. [Google Scholar] [CrossRef]
  25. Faghri, M.; Sparrow, E. Forced convection in a horizontal pipe subjected to nonlinear external natural convection and to external radiation. Int. J. Heat Mass Transfer. 1986, 23, 861–872. [Google Scholar] [CrossRef]
  26. Heywang, W. Resistivity Anomaly in Doped Barium Titanate. J. Am. Ceram. Soc. 1964, 47, 484–490. [Google Scholar] [CrossRef]
  27. Jonker, G. Some aspects of semiconducting barium titanate. Solid State Electron. 1964, 7, 895–903. [Google Scholar] [CrossRef]
  28. Kulwicki, B.; Purdes, A. Diffusion potentials in BaTiO3 and the theory of PTC materials. Ferroelectrics. 1970, 1, 253–263. [Google Scholar] [CrossRef]
  29. Ihrig, H.; Puschert, W. A systematic experimental and theoretical investigation of the grain boundary resistivities of n–doped BaTiO3 ceramics. J. Appl. Phys. 1977, 48, 3081–3088. [Google Scholar] [CrossRef]
  30. Mallick, G.; Emtage, P. Current–Voltage Characteristics of Semiconducting Barium Titanate Ceramic. J. Appl. Phys. 1968, 39, 3088–3094. [Google Scholar] [CrossRef]
  31. Preis, W.; Sitte, W. Modeling of transport properties of interfacially controlled electroceramics: Application to n–conducting barium titanate. J. Electroceram. 2011, 27, 83–88. [Google Scholar] [CrossRef]
  32. Kriegsmann, G. Pattern formation in microwave heated ceramics: cylinders and slabs. IMA Journal of Applied Mathematics. 2001, 66, 1–32. [Google Scholar] [CrossRef]
  33. Wacker, A.; Schöll, E. Criteria for stability in bistable electrical devices with S– or Z–shaped current voltage characteristic. J. Appl. Phys. 1995, 78, 7352–7357. [Google Scholar] [CrossRef]
  34. Romanowskii, V.; Watanabe, K. Multi-Stable Static Resistive States of High Tc Superconductors With Temperature-Decreasing Power Exponent of Their Voltage-Current Characteristic. IEEE Transactions on Applied Superconductivity. 2010, 20, 2134–2137. [Google Scholar] [CrossRef]
  35. Fisher, A.; Pfalz, M.; Vandewal, K.; Lenk, S.; Liero, M.; Glitzky, A.; Reineke, S. Full Electrothermal OLED Model Including Nonlinera Self–heating Effects. Phys. Rev. Applied. 2018, 10, 014023. [Google Scholar] [CrossRef]
  36. Elmer, F. Limit cycles of the ballast resistor caused by intrinsic instabilities. Z. Phys. B–Cond. Mat. 1992, 87, 377–386. [Google Scholar] [CrossRef]
  37. Krikkis, R. Analysis of High-Temperature Superconducting Current Leads: Multiple Solutions, Thermal Runaway, and Protection. J. 2023, 6, 302–317. [Google Scholar] [CrossRef]
  38. Gupta, T. Application of Zinc Oxide Varistors. J. Amer. Ceram Soc. 1990, 73, 1817–1840. [Google Scholar] [CrossRef]
  39. Clarke, D. Varistor Ceramics. J. Am. Ceram Soc. 1999, 82, 485–502. [Google Scholar] [CrossRef]
  40. Li, J.; Li, S.; Cheng, P.; Alim, M. Advances in ZnO–Bi2O3 based varistors. J Mater Sci: Mater Electron. 2015, 26, 4782–4809. [Google Scholar] [CrossRef]
  41. Mahan, G.; Philipp, H. Theory of conduction in ZnO varistors. J. Appl. Phys. 1979, 50, 2799–2812. [Google Scholar] [CrossRef]
  42. Blatter, G.; Greuter, F. Carrier transport through grain boundaries in semiconductors. Phys. Rev B. 1986, 33, 3952–3966. [Google Scholar] [CrossRef] [PubMed]
  43. Greuter, F.; Blatter, G. Electrical properties of grain boundaries in polycrystalline compound semiconductors. Semicond. Sci. Tech. 1990, 5, 111–137. [Google Scholar] [CrossRef]
  44. Eda, K. Destruction mechanism of ZnO varistors due to high currents. J. Appl. Phys. 1984, 56, 2948–2955. [Google Scholar] [CrossRef]
  45. Späck-Leigsnering, Y.; Gjonaj, E.; De Gersem, H.; Weiland, T.; Gießel, M.; Hinrichsen, V. Electroquasistatic-Thermal Modeling and Simulation of Station Class Surge Arresters. IEEE Trans. Magn. 2016, 42. [Google Scholar] [CrossRef]
  46. Seyyedbarzegar, S.; Mirzaie, M. Thermal balance diagram modelling of surge arrester for thermal stability analysis considering ZnO varistor degradation effect. IET Gener. Transm. Distrib. 2016, 10, 1570–1581. [Google Scholar] [CrossRef]
  47. Likitha, S.; Kanyakumari, M.; Pauly, P.J.; Aradhya, R.S.; Vasudev, N. Estimation of Critical Resistive Leakage Current of polymer housed ZnO surge arrester by electro-thermal modeling. J. Elec. Sys. Inf. Tech. 2018, 5, 861–873. [Google Scholar] [CrossRef]
  48. Cologna, M.; Raskova, B.; Raj, R. Flash Sintering of Nanograin Zirconia in < 5 s at 850°C. J. Am. Ceram Soc. 2010, 93, 3556–3559. [Google Scholar]
  49. Todd, R.; Zapata-Solvas, E.; Bonilla, R.; Sneddon, T.; Wilshaw, P. Electrical characteristics of flash sintering: thermal runaway of Joule heating. J. Eur. Ceram. Soc. 2015, 35, 1865–1877. [Google Scholar] [CrossRef]
  50. Yu, M.; Grasso, S.; Mckinnon, R.; Saunders, T.; Reece, M. Review of flash sintering: materials, mechanisms and modeling. Advances in Applied Ceramics. 2016. [Google Scholar] [CrossRef]
  51. Biesuz, M.; Sglavo, V. Flash sintering of ceramics. J. Eur. Ceram. Soc. 2019, 39, 115–145. [Google Scholar] [CrossRef]
  52. Zhou, H.; Li, X.; Zhu, Y.; et al. Review of flash sintering with strong electric field. High Voltage. 2022, 7, 1–11. [Google Scholar] [CrossRef]
  53. Hewitt, I.; Lacey, A.; Todd, R. A mathematical model for flash sintering. Math. Model. Nat. Phenom. 2015, 10, 77–89. [Google Scholar] [CrossRef]
  54. Downs, J.; Sglavo, V. Electric Field Assisted Sintering of Cubic Zirconia at 390°C. J. Am. Ceram. Soc. 2013, 96, 1342–1344. [Google Scholar] [CrossRef]
  55. Dong, Y.; Chen, I.W. Predicting the Onset of Flash Sintering. J. Am. Ceram. Soc. 2015, 98, 2333–2335. [Google Scholar] [CrossRef]
Preprints 100820 g001
Figure 1. (a) conductor geometry (b) energy balance on an elementary volume element (c) simplified electric circuit.
Figure 1. (a) conductor geometry (b) energy balance on an elementary volume element (c) simplified electric circuit.
Preprints 100820 g001
Preprints 100820 g002
Figure 2. Energy balance and uniform solutions.
Figure 2. Energy balance and uniform solutions.
Preprints 100820 g002
Preprints 100820 g003
Figure 3. Solution structure projected on ( v , Θ ) plane and instability threshold.
Figure 3. Solution structure projected on ( v , Θ ) plane and instability threshold.
Preprints 100820 g003
Preprints 100820 g004
Figure 4. Solution structure projected on ( v , j ) plane.
Figure 4. Solution structure projected on ( v , j ) plane.
Preprints 100820 g004
Preprints 100820 g005
Figure 5. Solution structure for the distributed model. Θ b = Θ ( 0 ) .
Figure 5. Solution structure for the distributed model. Θ b = Θ ( 0 ) .
Preprints 100820 g005
Preprints 100820 g006
Figure 6. Temperature profiles.
Figure 6. Temperature profiles.
Preprints 100820 g006
Preprints 100820 g007
Figure 7. Instability threshold for flash sintering.
Figure 7. Instability threshold for flash sintering.
Preprints 100820 g007
Preprints 100820 g008
Figure 8. Instability threshold for 3YSZ. Comparison with experimental data [49].
Figure 8. Instability threshold for 3YSZ. Comparison with experimental data [49].
Preprints 100820 g008
Preprints 100820 g009
Figure 9. Instability threshold for 8YSZ. Comparison with experimental data [48,54].
Figure 9. Instability threshold for 8YSZ. Comparison with experimental data [48,54].
Preprints 100820 g009
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
bsky
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2025 MDPI (Basel, Switzerland) unless otherwise stated