*Article* **Compact Thermal Model of the Pulse Transformer Taking into Account Nonlinearity of Heat Transfer** †

#### **Krzysztof Górecki 1,\*, Kalina Detka <sup>1</sup> and Krzysztof Górski <sup>2</sup>**


Received: 14 May 2020; Accepted: 29 May 2020; Published: 1 June 2020

**Abstract:** This paper presents a compact nonlinear thermal model of pulse transformers. The proposed model takes into account differentiation in values of the temperatures of a ferromagnetic core and each winding. The model is formulated in the form of an electric network realising electrothermal analogy. It consists of current sources representing power dissipated in the core and in each of the windings, capacitors representing thermal capacitances and controlled current sources modelling the influence of dissipated power on the thermal resistances in the proposed model. Both self-heating phenomena in each component of the transformer and mutual thermal couplings between each pair of these components are taken into account. A description of the elaborated model is presented, and the process to estimate the model parameters is proposed. The proposed model was verified experimentally for different transformers. Good agreement between the calculated and measured waveforms of each component temperature of the tested pulse transformers was obtained. Differences between the results of measurements and calculations did not exceed 9% for transformers with a toroidal core and 13% for planar transformers.

**Keywords:** nonlinear thermal model; SPICE; pulse transformer; thermal phenomena; self-heating; modelling; measurements

#### **1. Introduction**

Pulse transformers are an important component of switched-mode power converters [1–3]. These transformers have simple structure, and they consist of two kinds of components, i.e., a ferromagnetic core and at least two windings. During the operation of the considered device, an increase in temperature of each transformer component is observed [4–6]. This increase is a result of thermal phenomena occurring in the pulse transformer, such as self-heating in each component of the transformer and mutual thermal interaction between each pair of these components [7–9].

Knowing the core temperature and the windings temperatures is important from the point of view of electrical and magnetic properties of a pulse transformer. As is shown in [10,11], the temperature significantly changes the characteristics of ferromagnetic materials used to make a transformer core and causes a change in the resistances of the windings. In particular, an excessive increase in temperature can lead to damage in the insulation of the windings or can reduce the magnetic permeability of the core [8,11,12]. Additionally, an increase in the temperatures of electronic components causes a decrease in their lifetime [13,14].

In order to calculate the temperature waveforms of electronic components with thermal phenomena taken into account, a thermal model of these components is indispensable [15,16]. In many papers [10,11,17–25], thermal models of transformers are described, but they have disadvantages. One group of thermal models is microscopic models, which make it possible to calculate temperature distribution in the considered component. For example, in [11,18], the finite element method is used to determine the temperature distribution in a transformer, but at the same time the distribution of the wasted power per unit of volume in the transformer is assumed. Another group of thermal models is compact thermal models, which take into account only one temperature characterising the whole device [17].

Reference [26] presents three-dimensional (3-D) numerical compact thermal models of planar transformers. The DC thermal network of the cited model is inspired by the Delphi method, which allows obtaining shorter time of calculations than the finite volume method. The presented results do not illustrate the influence of dissipated power on the temperatures of the core and the windings. Additionally, this model is dedicated to the ANSYS software, and it is difficult to implement it in other software, e.g., in SPICE.

The similar approach to modelling thermal properties of electronic components is presented in [27,28]. The model presented in [27] uses homogenisation techniques to reduce calculation requirements. The advantage of this model is the reduction of the number of thermal factors to 6–8, and good agreement between the results of calculations and measurements can still be obtained. Unfortunately, the mentioned model could be used in the software dedicated only to a 3-D thermal analysis. In turn, the approach to the thermal modelling of components of electric machines presented in [28] requires time-consuming measurements of the tested prototype. In the paper [28], some temperature waveforms of the tested machine are presented.

Reference [29] describes a simplified form of an electrothermal model dedicated to planar magnetic components (inductors and transformers), which operate in the space industry. An important part of this model is the thermal model dedicated to the ANSYS program. Due to the fact that a lot of factors are taken into account, the network representation of the proposed model is complex and contains three subcircuits representing the thermal network of the transformer windings, the thermal network of the transformer core and the thermal network of the connection between the windings and the Printed Circuit Board (PCB).

References [30,31] are dedicated to parameters estimation of thermal models of electronics devices in the Delphi-inspired form. To this end, genetic algorithms are used. Unfortunately, for multiple heat sources in such models, calculations are time-consuming, and the total simulation time can reach 800 h.

Reference [32] presents a thermal model of a planar transformer. The form of this model is obtained on the basis of computation fluid dynamics. The structure of this model is adequate for planar transformers only. Unfortunately, in the paper [32], no results of experimental verification of this model are presented.

Compact thermal models of transformers are described in [10,11,17,18,26]; however, differences between the core and windings temperatures are typically not taken into account in the mentioned models. In addition, in the models described in [24,25], the dependence of the dissipation efficiency of heat generated in the device on the power dissipated in the transformer is omitted [33].

It is widely known [7,15,16,34–36] that some factors, such as ambient temperature, cooling systems or power dissipated in electronic devices, influence the efficiency of heat removal from the devices. Reference [37] presents a nonlinear thermal model of a planar transformer. In this model, the influence of power dissipated in particular components of the transformer on the efficiency of heat removal is taken into account.

This paper, which is an extended version of Reference [38], proposes a compact nonlinear thermal model (CNTM) of a pulse transformer based on a thermal model of a planar transformer described in [19]. In comparison to [19], which presents a linear thermal model, the nonlinearity of the heat transfer process is taken into account. In comparison to [38], a detailed description of the nonlinear thermal model of the transformer and the new results of measurements and calculations illustrating the effectiveness of this model for transformers including different ferromagnetic cores are presented in this paper. The model proposed by the authors' takes into account self-heating phenomena in all components of the transformer and mutual thermal couplings between each pair of these components.

The selected thermal models of transformers given in the literature are discussed in Section 2. The form of a nonlinear thermal model is presented in Section 3. The obtained results of measurements and calculations proving the effectiveness of the elaborated model are shown in Section 4. The advantage of this nonlinear thermal model over a linear thermal model is experimentally confirmed for the selected transformers containing cores with different shapes and made of different ferromagnetic materials. For all models described in the following sections, thermal resistance and thermal capacitance are given in K/W and J/K, respectively, whereas all temperatures are expressed in Celsius degrees.

#### **2. Selected Thermal Models of Transformers in the Literature**

A network representation of classical linear thermal models of transformers [10,39] is shown in Figure 1. It can be seen that only one internal temperature of the whole transformer is used. In this model, differences in the temperatures of the core and the windings are not taken into account.

**Figure 1.** Network representation of the thermal model described in [10,39].

In the presented model, the controlled current source G1 represents the sum of the power dissipated in the core and in the windings of the transformer. Thermal capacitance is represented by capacitor Cth, while Rth1 is the thermal resistance characterising heat convection and Rth2 is thermal resistance characterising the heat radiation from the surface of the examined device. Voltage source Ta is the ambient temperature, and the voltage in node Tj is the temperature of the transformer.

In [17,33], a thermal model of magnetic components (transformers and inductors) is proposed. This model enables calculating the temperature difference between the core and ambient temperature (ΔTC) and the temperature difference between the windings and ambient temperature (ΔTW) by taking into account self-heating and mutual thermal couplings between the core and the windings. The network representation of this model is shown in Figure 2.

However, in this model, only single thermal time constants for the windings (RthW and CthW) and for the core (RthC and CthC) are taken into account. Current sources GPC and GPW describe power losses in the core and in the windings, respectively; whereas current sources GPC1 and GPW1 model the influence of mutual thermal couplings between the core and the winding on the temperature differences ΔTW and ΔTC.

**Figure 2.** Network representation of the thermal model of magnetic components from [18].

The thermal model of a transformer proposed in [17], of which the diagram is presented in Figure 3, makes it possible to calculate the temperature of the core TC and the temperatures of both windings TW by taking into account self-heating and mutual thermal couplings between the core and the windings. This model has the form of RC Foster networks excited by current sources representing the powers dissipated in the core (PthC and PthWC1) and in the windings (PthW and PthWC1).

**Figure 3.** Network representation of the thermal model of a transformer proposed in [17].

In order to take into account thermal couplings between the core and the windings, controlled current sources PthWC1 and PthCW1 are applied. The voltage source Ta represents ambient temperature. In the described model, common RC networks are used to model self-heating and mutual thermal couplings between the components of the transformer.

A disadvantage of thermal models of transformers described in this section is they do not take into account differentiations of windings temperatures, nonlinearities of thermal properties and thermal couplings occurring between the components of the transformer. Therefore, in the following Section, the authors' nonlinear thermal model of the transformer is proposed. In this model, nonlinearities of thermal phenomena and thermal couplings between the components of the transformer are taken into account.

#### **3. Proposed Nonlinear Thermal Model of Pulse Transformers**

As is shown in [17,19], temperature distributions on each of the windings and on the core in pulse transformers are practically uniform. Therefore, their thermal properties can be described with the use of a compact thermal model [7,19]. On the other hand, the temperatures of the core and the windings can be significantly different from each other [7]. Therefore, in such models, the differences of the core and windings temperatures should be taken into account. In each transformer component, a self-healing phenomenon occurs, and additionally, thermal interaction between each pair of transformer components is observed.

The nonlinear thermal model of the transformer worked out by the authors has the form of a subcircuit dedicated for SPICE. This form is inspired by the linear thermal model of a planar transformer proposed in [19] and by the nonlinear thermal model of semiconductor devices described in [37,40]. In the new model, the differences in the temperatures of components in the transformer (the core and all the windings) are taken into account. These temperatures result from self-heating phenomena in every component and mutual thermal couplings between each pair of the components of the transformer. The presented model makes it possible to calculate the temperature waveforms of the core and each winding by taking into account both the mentioned phenomena. When formulating the considered model, based on the results of measurements shown in [7], changes in the powers dissipated in each component influence only the values of transformer thermal resistances in the thermal model, whereas thermal capacitances do not depend on power. The network representation of the proposed model is presented in Figure 4.

**Figure 4.** Network representation of a compact nonlinear thermal model (CNTM) of a pulse transformer.

In this network, nine circuits can be distinguished. Three of them located on the left-hand side of Figure 4 are used to calculate the temperature waveforms of particular components of the transformer—the primary winding TW1, the secondary winding TW2 and the core TC. These temperatures (given in ◦C) correspond to the voltages (given in V) in the nodes and are denoted as TW1, TW2 and TC, respectively. Current sources IC, IW1 and IW2 model power losses in the core and the two windings, respectively.

Circuits including capacitors and controlled current sources model self-transient thermal impedances of the primary winding (CW11, ... , CW1n and GW11, ... , GW1n), of the secondary winding (CW21, ... , CW2n and GW21, ... , GW2n) and of the core (CC1, ... , CCn and GC1, ... , GCn), respectively. The voltages on these circuits (given in V) correspond to differences (given in ◦C) between the temperatures of the transformer components and ambient temperature, resulting from self-heating phenomena.

The controlled voltage sources E1, E2 and E3 represent the influence of mutual thermal couplings between the components of the transformer on the temperatures of these components. Capacitors represent the thermal capacitances of all the elements in the heat flow path, whereas the controlled current sources represent nonlinear thermal resistances between the elements in the heat flow path.

The output voltage of the controlled voltage source E1 is equal to the sum of voltages in nodes TW11 and TWC1, the output voltage of the controlled voltage source E2 is equal to the sum of voltages in nodes TW21 and TWC2, and the output voltage of the controlled voltage source E3 is equal to the sum of voltages in nodes TCW1 and TCW2. Voltage sources V1, V2 and V3 represent ambient temperature.

Six other subcircuits are used to model mutual thermal couplings between each pair of the transformer components. Current sources in these subcircuits represent power (given in W) dissipated in particular components of the transformer (IW12 and IW1C in the primary winding, IW21 and IW2C in the secondary winding and IC1 and IC2 in the core). Networks containing capacitors and controlled current sources connected to the above-mentioned current sources model mutual transient thermal impedances between each pair of components of the transformer.

All self-transient and mutual transient thermal impedances can be described by Equation (1) [7,15]:

$$Z\_{th}(t) = R\_{th} \cdot \left[1 - \sum\_{i=1}^{N} a\_i \cdot \exp\left(-\frac{t}{\tau\_{thi}}\right)\right] \tag{1}$$

where Rth is the thermal resistance, ai is the coefficient (without unit) corresponding to thermal time constant τthi (given in s), and N is the number of thermal time constants.

The dependence of Rth on the dissipated power is described as:

$$R\_{th} = R\_{th0} \cdot \left[1 + \alpha \cdot \exp\left(-\frac{p - p\_0}{b}\right)\right] \tag{2}$$

where Rth0 denotes the minimum value of the thermal resistance, p denotes the power dissipated in a heating component of the transformer, α is the parameter without unit, and p0 and b are model parameters given in W.

Changes in values of thermal resistance are modelled by the controlled current source Gi. The output current of the controlled source is described as:

$$G\_{\bar{l}} = V\_{\bar{G}\bar{l}} / (a\_{\bar{l}} \cdot \mathbb{R}\_{th}),\tag{3}$$

where VGi denotes the voltage on the current source Gi.

Thermal capacitances are given as:

$$C\_i = \tau\_{thi} / \left( a\_i \cdot R\_{th} \right). \tag{4}$$

The values of the model parameters are estimated using the results of measurements of self- and mutual transient thermal impedances existing in the transformer thermal model. Measurements of such parameters at different powers dissipated in the core and in the windings of the tested transformer are realised by the method described in [7]. The values of parameters in Equation (1) are estimated for each transient thermal impedance using the method described in [15,41]. Next, parameters α, p0 and b in Equation (2) are estimated for self-thermal and mutual thermal resistance in the transformer model by using local estimation [15].

#### **4. Results**

In order to verify the presented model and its practical use, measurements and calculations of the temperature waveform of each component of the tested transformers, which contained ferromagnetic cores with different shapes and sizes and were made of different materials, were performed. A planar transformer with a ferrite core and transformers with ring cores made of different materials were measured and modelled as examples.

The planar transformer, of which the cuboidal core dimension was 22 mm × 16 mm × 9 mm, was made of ferrite material 3F3 and contained windings in the form of printed paths on laminate FR-4, which was 1 mm thick. The primary winding contained three turns with a width of 2.5 mm, and the secondary winding contained four turns with a width of 1 mm. Transformers with a toroidal core contained identical primary and secondary windings. On each of them, 20 turns of copper wire in the enamel with a diameter of 0.8 mm were wound. The toroidal core had an external diameter equal

to 26 mm, the internal diameter equal to 16 mm and a width equal to 11 mm. Cores made of toroidal powdered iron (RTP), ferrites (RTF) and nanocrystals (RTN) were used for toroidal transformers.

The values of the parameters in the thermal model were estimated for transformers containing ferromagnetic cores with different shapes and dimensions and made of different ferromagnetic materials.

The values of the parameters in Equation (2), which represents the proposed compact nonlinear thermal model of the transformer (CNMT), are summarized in Table 1 (for a toroidal transformer) and Table 2 (for a planar transformer). The values of the parameters presented in Tables 1 and 2 describe the model of the network representation shown in Figure 4.

**Table 1.** Values of the parameters in Equation (2) for a toroidal transformer with a toroidal powdered iron core (RTP).



**Table 2.** Values of the parameters in Equation (2) for a planar transformer.

Comparing the values of the parameters describing thermal resistances in the thermal model of the toroidal transformer, it is obvious that the values of parameters α, p0 and b were nearly the same for all considered thermal resistances. Differences were observed in the values of the parameter Rth0. This meant that courses RthW1(pW1), RthW1W2(pW1) and RthW1C(pW1) were nearly parallel. In contrast, big differences were observed between the values of the parameters describing the dependences of thermal resistances in the thermal model of the planar transformer on powers PW1 and PC.

The values of thermal capacitances were estimated with the use of the method described in [15]. As an example, in Table 3, the values of these parameters obtained for a planar transformer are summarised.


**Table 3.** Values of the parameters ai and τthi of the selected self-transient and mutual transient thermal impedances in a thermal model of a planar transformer.

It can be clearly seen that, in the considered self-transient and mutual transient thermal impedances, different numbers of thermal time constants occurred. In self-transient thermal impedances, three or four thermal time constants were used, whereas mutual transient thermal impedances were described

with two or three thermal time constants. The values of the considered thermal time constants were in a range from 40 μs to over 1000 s.

For example, by means of Equation (2), the measured (indicated by points) and modelled (indicated by lines) dependences of thermal resistances in the thermal model of the transformer on the power dissipated in one of the components of the transformer are shown in Figures 5 and 6.

Figure 5 illustrates the influence of power PW1 dissipated in the primary winding of the transformer containing a toroidal core made of powdered iron on the thermal resistance of the winding RthW1 (blue colour) and on mutual thermal resistances between this winding and the secondary winding RthW1W2 (black colour), as well as those between the core and the primary winding RthW1C (red colour).

Figure 6a illustrates the influence of power dissipated in the primary winding of the planar transformer on the thermal resistances RthW1 (blue colour), RthW1W2 (black colour) and RthW1C (red colour). Figure 6b illustrates the influence of power dissipated in the core of the transformer on the thermal resistance of the core RthC (blue colour) and the mutual thermal resistances between the core and both windings RthCW1 (red colour) and RthCW2 (green colour).

As it is obvious in Figures 5 and 6, it is possible to accurately model the measured dependences of the considered thermal resistance on the dissipated power using Equation (2). Visible differences between self-thermal and mutual thermal resistances were observed for both considered transformers. These thermal resistances differed from one another. It is worth noticing that the considered differences were bigger for the planar transformer than for the toroidal transformer. Changes in dissipated power can cause changes in thermal resistance even by 25%.

**Figure 5.** Measured (indicated by points) and modelled (indicated by lines) dependences of the selected thermal resistances in the thermal model of the transformer with a toroidal core made of powdered iron on the power dissipated in the primary winding.

**Figure 6.** Measured (indicated by points) and modelled (indicated by lines) dependences of the selected thermal resistances in the thermal model of the planar transformer on the power dissipated in the primary winding (**a**) and in the core (**b**).

Calculations and measurements were performed at power dissipation in only one of the components of the tested transformers and at simultaneous dissipation of power in different components of the tested devices. Power dissipated in each component of the transformer always had a shape of a single rectangular impulse with a long duration time. The results of calculations obtained by means of the nonlinear thermal model were compared to the results of calculations performed by means of the linear thermal model described in [19] and the results of measurements performed with the use of a pyrometer. The windings and the core were excited with different powers. The successive figures (Figures 7–12) present the calculated and measured waveforms of the temperatures of the primary winding TW1, the secondary winding TW2 and the core TC of the tested transformers. In these figures, the results of measurements is represented by points, the results of calculations performed using the CNTM is represented by solid lines, and the results of calculations using the compact linear thermal model (CLTM) is represented by dashed lines [19].

Calculations were performed for the values of parameters describing the nonlinear thermal model of transformers according to the principles shown in Section 3. On the other hand, for the model from [19], the values of parameters estimated at the lowest measured values of power dissipated in each component of the transformers were used.

Figure 7 presents the measured and calculated waveforms of temperatures TW1 and TC in the toroidal transformer with the RTP core at a dissipated power PW1 of 2.83 W in the primary winding.

It can be observed that the nonlinear thermal model makes it possible to obtain very good agreement between the results of calculations and measurements. In contrast, the values of the considered temperatures obtained with the use of calculations performed with the linear thermal model were higher than the results of measurements by even 40 ◦C.

**Figure 7.** Measured and calculated temperature waveforms of the primary winding TW1 and the core TC of a toroidal transformer with an RTP at a dissipated power PW1 of 2.83 W in the primary winding.

Figure 8 shows the measured and calculated temperature waveforms of the primary winding TW1, the secondary winding TW2 and the core TC for a transformer containing an RTN. These waveforms were obtained, while the primary winding of the tested transformer was excited by a single rectangular pulse with duration times equal to 4000 s (Figure 8a) and 5000 s (Figure 8b). PW1 was 1.05 W in the case presented in Figure 8a, and PW1 was equal to 2.82 W in the case presented in Figure 8b.

As it is seen, very good agreement between the results of calculations performed with the use of the CNTM and measurements was achieved for both the values of power dissipated in this transformer. In contrast, for the linear thermal model, an excess of temperature of each component of the transformer over ambient temperature was overestimated by 10% at lower values of the considered powers and even 50% overestimated for higher values of the considered dissipated power.

**Figure 8.** Measured and calculated temperature waveforms of the components of a transformer with a toroidal nanocrystalline core (RTN) at dissipated powers PW1 of 1.05 W (**a**) and 2.82 W (**b**) in the primary winding.

Figure 9 illustrates the measured and calculated temperature waveforms of the primary winding TW1, the secondary winding TW2 and the core TC for the transformer with an RTF. The primary winding was stimulated by a single rectangular impulse with a duration time equal to 7000 s at PW1 of 0.77 W.

By analysing the temperature waveforms of the core TC, the primary winding TW1 and the secondary winding TW2 presented in Figure 9, it can be observed that the results of calculations performed by means of the nonlinear thermal model assured better agreement with the results of measurements than the results of calculations obtained by means of the model described in [19]. The obtained difference between the results of calculations performed with the use of the linear thermal model and the measurements results was 10%.

Figure 10 presents the measured and calculated waveforms of the temperature of the primary winding TW1, the temperature of the secondary winding TW2 and the core temperature TC of the planar transformer by stimulating the primary winding with a single rectangular impulse having a duration time equal to 5500 s and PW1 of 1 W (Figure 10a) and with a single rectangular impulse having a duration time of 5000 s and PC of 2 W (Figure 10b).

As one can notice, the results of calculations by means of the nonlinear thermal model assured better agreement with the results of measurements than the results of calculations performed by means of the linear thermal model. In Figure 10a, it is visible that the difference between the waveforms of temperatures TW1 and TC obtained with the use of the considered models was more than 2 ◦C. The biggest difference between the results of calculations was observed for temperature TW2. In Figure 10b, the considered differences were bigger than in the case presented in Figure 10a. These differences reached even 20 ◦C. When power was dissipated in the core only, the temperatures of both windings were nearly the same.

**Figure 9.** Measured and calculated waveforms of the temperature of the primary winding TW1, the temperature of the secondary winding TW2 and the temperature of the core TC at a dissipated power PW1 of 0.77 W in the primary winding.

**Figure 10.** Measured and calculated waveforms of the temperature of the primary winding TW1, the temperature of the secondary winding TW2 and the temperature of the core TC of the planar transformer at PW1 of 1 W in the primary winding (**a**) and PC of 2 W in the core (**b**).

The results of measurements and calculations presented above corresponded to untypical situations, when power was dissipated in one of the transformer components only. Figures 11 and 12 present the results of measurements and calculations of temperature waveforms of components of the selected transformers obtained in the case when power was dissipated in both the components of the selected transformers.

Figure 11 shows the measured and calculated temperature waveforms of the primary winding TW1, the secondary winding TW2 and the core TC for the transformer containing the toroidal ferrite core (RTF) when it was excited by power dissipated simultaneously in the primary and secondary windings with a single rectangular pulse having duration times equal to 4500 s (Figure 11a) and 5500 s (Figure 11b). PW1 was equal to 0.88 W and PW2 was equal to 1 W in the case presented in Figure 11a, whereas PW1 was equal to 2.4 W and PW2 was equal to 2.7 W in the case presented in Figure 11b.

By analysing temperature waveforms of the core and the windings of the transformer containing an RTF (Figure 11), it is easy to observe that the difference between the results of calculations performed using the linear thermal model and the results of measurements was up to 20 ◦C. In contrast, the results of calculations made using the nonlinear thermal model showed very good agreement with the results of measurements. Differences between the results of measurements and calculations performed with the linear thermal model were bigger at higher powers dissipated in the components of the transformer.

Figure 12 presents the calculated and measured waveforms of the temperatures of each component of the planar transformer by simultaneously stimulating the core and the primary winding with a single rectangular impulse having a duration time equal to 4500 s, PW1 of 2.3 W and PC of 2 W.

**Figure 11.** Measured and calculated temperature waveforms of the primary winding TW1, the secondary winding TW2 and the core TC with simultaneous power dissipation in the windings: (**a**) PW1 = 0.88 W and PW2 = 1.05 W; (**b**) PW1 = 2.4W and PW2 = 2.7 W.

**Figure 12.** Measured and calculated waveforms of the temperature of the primary winding TW1, the temperature of the secondary winding TW2 and the temperature of the core TC of the planar transformer with an impulse power PC of 2 W in the core and an impulse power PW1 of 2.3 W in the primary winding.

As one can notice, the results of calculations by means of the nonlinear thermal model assured better agreement with the results of measurements than the results of calculations performed by means of the linear thermal model. The observed differences between the results of measurements and the results of calculations performed by means of the linear thermal model exceeded even 40 ◦C, whereas the difference between the results of calculations performed by means of the nonlinear thermal model and the results of measurements was no more than about 5 ◦C. It is also worth noticing that in the considered operation conditions the temperatures of the transformer components were low and did not exceed 10 ◦C.

Table 4 contains the values of the maximum errors of the temperatures of transformers components δTW1, δTW2 and δTC obtained using the CLTM and the CNTM of the transformer analysed in this section.


**Table 4.** Values of the maximum relative error of the temperatures of transformers components calculated with the CNTM and the compact linear thermal model (CLTM).

It is clearly seen in Table 4 that the considered relative error of calculations of temperatures obtained using the CLTM achieved 48.8% for the primary winding temperature of the transformer with an RTP core whereas such an error of calculations of the primary winding temperature obtained using the CNTM was equal to 8.16% for the same transformer. It is also worth noticing that the value of the relative error of the temperatures of transformers components obtained using the CNTM did not exceed 9% for all the considered transformers containing a toroidal core made of different materials. It was also observed that an increase in power dissipated in the transformer components caused an increase of the value of the relative error of calculations for both the considered models.

#### **5. Conclusions**

This paper describes a new CNTM of a pulse transformer elaborated by the authors. It allows determining the waveforms of the temperatures of the core and each winding. Calculations were performed by taking into account a self-heating phenomenon and mutual thermal interaction between the transformer components. In our model, the dependences of self-thermal and mutual thermal resistances depend on power dissipated in the transformer components. The proposed model was verified experimentally for transformers including ferromagnetic cores with different shapes and made of different materials. High accuracy of this model and its advantage over the CLTM were also demonstrated. The results of calculations and measurements presented in this paper also confirmed that power dissipated in the transformer components influenced mostly thermal resistance and its influence on thermal capacity was neglected.

It is also worth noticing that the obtained differences between the values of the temperatures of the transformer components can be equal to 50 ◦C. Such big differences justify the use different temperatures of the windings and the core in the thermal model. The linear thermal model makes it possible to obtain good agreement between the results of calculations and measurements only for very low powers dissipated in each component of the examined transformers. The nonlinear thermal model makes it possible to obtain good agreement between the results of calculations and measurements over a wide range of power dissipated in the windings or in the core of the examined transformers by considering the dependence of Rth on power.

The presented results of calculations and measurements proved that the nonlinear thermal model of a pulse transformer proposed in this paper is able to accurately describe dynamic thermal properties of the transformer including cores with different dimensions and shapes and made of different ferromagnetic materials. It was also shown that using this model one can obtain accurate results of calculations at different conditions of power dissipation in the tested transformers.

It is also worth noticing that the proposed model is universal and it can be useful for different type of transformers such as planar transformers or transformers containing a toroidal core. Additionally, this model takes into account properties of materials used to make the core. The biggest advantage of the proposed model is its simple structure and that easy implementation, e.g., in the SPICE program, which is widely used by designers of electronic circuits. In addition, analyses made with the proposed model are not time-consuming, which is very important from the economic point of view. The proposed model can be used to formulate electrothermal models of transformers, which take into account the nonlinearity of a heat removal process.

The CNTM of the pulse transformer proposed in this paper can be useful in designing switch-mode power supplies, and it can be used as a component of an electrothermal model of transformers.

**Author Contributions:** Conceptualization (K.G. (Krzysztof Górecki)); methodology (K.G. (Krzysztof Górecki), K.G. (Krzysztof Górski)); validation (K.G. (Krzysztof Górecki), K.D., K.G. (Krzysztof Górski)); investigation (K.G. (Krzysztof Górecki), K.G. (Krzysztof Górski)); resources (K.D.); writing—original draft preparation (K.G. (Krzysztof Górecki), K.D.); writing—review and editing (K.G. (Krzysztof Górecki), K.D.); visualization (K.G. (Krzysztof Górecki), K.G. (Krzysztof Górski)); supervision (K.G. (Krzysztof Górecki)). All authors have read and agreed to the published version of the manuscript.

**Funding:** The project is financed in the framework of the program by the Ministry of Science and Higher Education called "Regionalna Inicjatywa Doskonało´sci" in the years 2019–2022 (project number: 006/RID/2018/19; the sum of financing: 11,870,000 PLN).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations and Notations**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Electro-Thermal Simulation of Vertical VO2 Thermal-Electronic Circuit Elements**

#### **Mahmoud Darwish, Péter Neumann, and János Mizsei and László Pohl \***

Department of Electron Devices, Budapest University of Technology and Economics, 1117 Budapest, Hungary; mahmoud@eet.bme.hu (M.D.); neumann@eet.bme.hu (P.N.); mizsei@eet.bme.hu (J.M.)

**\*** Correspondence: pohl@eet.bme.hu

Received: 15 May 2020; Accepted: 1 July 2020; Published: 3 July 2020

**Abstract:** Advancement of classical silicon-based circuit technology is approaching maturity and saturation. The worldwide research is now focusing wide range of potential technologies for the "More than Moore" era. One of these technologies is thermal-electronic logic circuits based on the semiconductor-to-metal phase transition of vanadium dioxide, a possible future logic circuits to replace the conventional circuits. In thermal-electronic circuits, information flows in a combination of thermal and electronic signals. Design of these circuits will be possible once appropriate device models become available. Characteristics of vanadium dioxide are under research by preparing structures in laboratory and their validation by simulation models. Modeling and simulation of these devices is challenging due to several nonlinearities, discussed in this article. Introduction of custom finite volumes method simulator has however improved handling of special properties of vanadium dioxide. This paper presents modeling and electro-thermal simulation of vertically structured devices of different dimensions, 10 nm to 300 nm layer thicknesses and 200 nm to 30 μm radii. Results of this research will facilitate determination of sample sizes in the next phase of device modeling.

**Keywords:** beyond CMOS; VO2; thermal-electronic circuits; electro-thermal simulation; vertical structure

#### **1. Introduction**

Miniaturization of electronic devices' feature size is the primary driver of computing development. Scaling down has been described by Moore's law for many years. Recently, complimentary metal–oxide–semiconductor (CMOS) scaling down has been slower than ever before owing to scaling limits that appear on very small dimensions. The technological advancement described by Moore's Law for CMOS technology since the past 50 years is expected to flatten out completely by 2025 [1,2]. Finding "Beyond CMOS" solutions has become more essential to keep the computing advancement flourishing. It takes about 10 years for a new technology to move from the laboratory to mass production, for example, FinFET [1]. There are several potential "More than Moore" devices in laboratories [3], but currently, none is in a phase to state that this will be the successor.

A report commissioned by the Defense Technical Information Center (DTIC) [4] proposes four basic computational models: classical digital computing, analog computing, neuro-inspired computing and quantum computing. The latter three areas are suitable for solving special problems much more effectively than before; however, they are generally not able to replace traditional tasks where classical digital computing performs well, such as 3D graphics, mobile devices and similar others.

Architecture specialization can help to increase computing power for a while [5], for example, the use of tensor processing units (TPUs) developed for artificial intelligence (AI) [6] or the wider use of field-programmable gate arrays (FPGAs) in data centers [7]. The degree of specialization is well illustrated by the fact that there were 4 separate accelerators in the Apple A4 system on a chip (SoC), 28 in the A8, and more than 40 in the A12 [8]. Another way to increase the performance of CMOS

technology-based systems is to use design and technological developments such as chip stacking in 3D using through-silicon vias (TSVs) [9], advanced energy management, near-threshold voltage (NTV) operation and application of an increased number of metal layers [3].

Following are some more promising future devices compatible with CMOS technology [10]:


In this article, we discuss a potential device that is compatible with CMOS integration and in which information transmits by two physical processes: electrical and thermal. This device is the phonsistor (a portmanteau of "phonon transistor") [18]. A digital circuit can be built from this device, which is called thermal-electronic logic circuits (TELC) [19,20]. TELC can be implemented using materials that have semiconductor-to-metal transition (SMT). Vanadium dioxide (VO2) has shown properties as a good candidate for implementation of TELC structures with its SMT at around 67 ◦C [21].

Vanadium dioxide has some special properties whose investigation is not yet complete. Its optical properties [22], the relationship between SMT and structural phase transition [23,24], the nature of SMT [20,25] or the monoclinic metallic phase [26] are studied. Its applications cover many areas such as thermal rectification [27,28], high performance electromechanical switches [29], smart window coatings [30], neuromorphic devices [31,32] and non-volatile memory arrays [32]. Phonsistor and TELC are not the only transistor-type applications, we also find examples of purely thermal transistors in the literature [33,34].

VO2 encounters electrical resistivity change when it exposes to heat. The material shows much lesser ability to conduct electrical current at room temperature than at higher temperatures. Moreover, the resistivity drops steeply around 67 ◦C (340 K), changing the material electrical properties to a conductor. The transition from a low-conducting phase to a high-conducting phase is called "semiconductor-to-metal transition" (SMT). Phase transition occurs due to the sudden change in the structure of VO2 from the tetragonal structure (at low temperatures) to the monoclinic structure (at high temperatures) [19]. Figure 1 illustrates the electric resistivity change by 3 to 4 orders of magnitude versus temperature. The electro-thermal behavior of VO2 is highly technology-dependent. It also depends on the manufacturing process and the substrate [35–37].

**Figure 1.** Temperature dependence of the resistivity of vanadium dioxide for (100) VO2; based on [20].

A simple phonsistor consists of two components: a heating resistor and an SMT resistor. The heat generated by Joule heating in the resistor transports by phonons to the SMT resistor by diffusion. This process is similar to the diffusion of minority charge carriers at the base of a bipolar transistor. The switching voltage and holding current of the SMT resistor decrease due to external heating; this is the role of the heating resistor [38]. The abrupt variation in the VO2 electric resistivity makes it capable of switching on and off electrically depending on its thermal state. Hence the term "thermal-electronic logic circuits" (TELC) is used.

Figure 2 shows a sample TELC structure where R1, R2, R3, R4, and R5 are conventional resistors and V1, V2, and V3 are VO2 SMT resistors. For instance, V1 switches ON if R1 AND R2 are both ON. V2 switches ON if R3 OR R4 is ON. Combinations that are more sophisticated can be implemented by altering the ambient temperature, driving current, elements sizes, distances between the resistors and their locations. When the driving current flows in the conventional resistors, the Joule-heating will heat them and the generated heat will diffuse to the surrounding by convection. With enough heating delivered, VO2 SMT resistors will reach the desired temperature to drop the resistivity value and transition to the conducting phase. This action allows current to pass through VO2 SMT resistors. Furthermore, heat dissipation layers inserted between each section of the structure is important to prevent heat overlapping.

**Figure 2.** Implementation of VO2 thermal switching in a TELC circuit.

Many problems still need solutions to implement the circuit shown. For example, how the dimensions of an SMT resistor affect the electro-thermal characteristics of the device needs clarification. This paper investigates that how the characteristics of vertical VO2 resistors depend on their sizes by simulation. Previously, we studied lateral structures by performing high-resolution electro-thermal simulations using a range of different dimensions [39]. Measurement results for lateral TELC structures are also available in [40].

The goal of our long-term study is designing and creating real TELC structures that require a comprehensive study of VO2 thin films' behavior along with how these structures act at different dimensions, arrangements and ambient temperatures. In this article, we examine VO2 based vertical structures at different dimensions.

#### **2. Experimentation**

Standard CMOS technology steps were used to prepare the samples for vertical thermo-electrical device test in Figure 3. The substrate was an n-type Si with 3.5 × <sup>10</sup><sup>17</sup> cm−<sup>3</sup> dopant concentration for better spreading resistance of the vertical structure. The 100 nm SiO2 insulator layer was grown by thermal oxidation. A window was etched into the insulator layer that will be the active conductive channel shape of the vertical structure. The VO2 layer (approximately 50 nm thin) was deposited by RF sputtering at 650 ◦C high temperature. The radii of the window changed between 10 μm up to 100 μm.

**Figure 3.** Cross-sectional view of the vertical VO2 structure prepared in the lab.

For the tailoring of SMT layer, wet chemistry was applied to remove the surround cover area of the window. The top and collector Pt contacts (approx. 10 nm thin) of devices were also prepared by RF sputtering. Standard optical lithography was applied to transfer the patterns during technology steps. During the vertical structure preparation, lateral structures were prepared to check the SMT layer quality and functional operation. The dimensions of VO2 channels in lateral structures were 100 μm × 700 μm.

Probe station was used for electrical measurement of samples. The electrical measurement was carried out by Keithley SourceMeter and samples were elevated at different environmental temperatures by a Cole-Parmer Thermostat System.

#### **3. Electro-Thermal Simulation**

To study thermal-electrical circuits by simulation, a distributed parametric simulator is required, typically with some finite algorithm (FDM, FEM, FVM). The difficulty of simulation is illustrated by the fact that the resistivity of VO2 changes rapidly as a function of temperature (see Figure 1) and the characteristic has a hysteresis.

Steep characteristic is the bigger problem. One consequence of this is developing of more stable working points. Figure 4 shows a VO2 resistor. If a small current is applied to the resistor, whose temperature has increased (because of Joule heating), it does not reach the level required for the phase change, the current density will be of same magnitude throughout the resistor volume. At higher

driving currents, the temperature begins to rise. In the part of resistor that first reaches the required temperature for phase change, the resistivity suddenly decreases and the current density increases. In the remaining resistor, however, the current density decreases and thus the temperature also decreases. The positive feedback results in vast majority of current flowing in a hot channel between the anode and cathode. If the material of resistor, structure and environment are homogeneous, the channel forms in the centerline of resistor. However, if this is not the case, the channel may create elsewhere, or even more channels may create. In consequence, a steady-state simulation that takes into account the process of current and temperature distribution must be performed. This example is achieved by time-domain simulation.

(**c**) **Figure 4.** A 200 μm × 20 μm × 0.5 μm VO2 resistor: (**a**) Structure; (**b**) Temperature distribution at low current (2 mA); (**c**) Temperature distribution at high current (8 mA).

The structure of the device shown in Figure 4 and the boundary conditions of the simulation are the same as those of the device shown in [41]. The 200 μm × 20 μm × 0.5 μm VO2 resistor is directly connected to the platinum anode and cathode. There was 25 ◦C air above the device, and the bottom was connected to a cold plate at 50 ◦C. In the simulation, the driving current was increased in 2 mA increments up to 20 mA and then decreased in the same increments to 0. The temperature distribution is shown in Figure 4b is for the increasing 2 mA, while the distribution is shown in Figure 4c is for the decreasing 8 mA. In [41], the process of channel formation and termination can be viewed in a series of images.

Another problem arising from the steep characteristic is that the simulation is difficult to converge. In the literature, we find that either electro-thermal simulation results for very specific VO2 structures [29,42], or they achieve the convergence by applying special boundary conditions [43]. In the case of thermal-electrical circuits, these solutions cannot be applied. Our attempts with ANSYS ended in failure, so we made the custom, finite volumes method based simulator developed at the department suitable for handling the special properties of vanadium dioxide [41].

#### *3.1. Computational Method (SUNRED)*

In the custom simulator, the examined structure divides by a grid into elementary rectangular cells and each cell is filled with a homogeneous material see Figure 5a. The real structure is three-dimensional as the figure shows a layer of it for clarity. The elementary cells are modeled by a circuit. Adjacent elementary cells are connected by an electrical and a thermal line (the ground is not marked), as shown in Figure 5b. The model obtained from the finite volume method gives the internal structure of elementary cells [44].

A typical structure is shown in Figure 5c. Dissipation caused by current flowing in the electrical part is realized by the dependent power source of the thermal model. Electrical and thermal circuit elements are temperature dependent. The computational task is to determine the voltages and temperatures with knowledge of the excitations and boundary conditions. The simulator performs the solution by successive network reduction method (SUNRED). Figure 5d, which merges cell pairs in each step, knocks out the internal nodes from the equations until finally, a single cell remains. The boundary conditions are applied to this cell, and then the voltages and temperatures of the internal nodes are calculated in a series of successive backward substitution steps.

**Figure 5.** Operation of the custom solver (for clarity in 2D): (**a**) Dividing the model into elementary cells; (**b**) Model network built from elementary cells; (**c**) Electro-thermal elementary cell, red: thermal, blue: electrical; (**d**) Successive network reduction.

In the time domain, the solution is done by the backward Euler method, and the nonlinear calculation is by successive approximation. A heuristic method is used to deal with the extreme nonlinearity of VO2, details of which are described in [41]. Having more than one value of resistivity at one temperature value, in addition to the memory effect in VO2, make the hysteresis loop essential in the materials resistivity function. The hysteresis of the resistivity of VO2 is modeled by a parallelogram, see Figure 6.

#### *3.2. Hysteresis Model*

The resistivity function consists of two domains in SUNRED, normal domain and hysteresis domain. In the hysteresis domain, resistivity is multi-valued and any resistance-temperature combination can occur, not just the edges. For instance, when the structure temperature is raised starting from the room temperature 25 ◦C, resistivity drops linearly in the temperature range (T0 = 25 ◦C, T1 = 58 ◦C). In the following temperature range (T1 = 58 ◦C, T2 = 64 ◦C), resistivity remains constant in the hysteresis domain. The material phase change begins above T2 and hence triggers the semiconductor-to-metal transition. SMT leads to a sharp drop in resistivity in the temperature range (T2 = 64 ◦C, T4 = 73 ◦C). At 73 ◦C, the material encounters a full phase transition and resistivity above this temperature is always constant. If the temperature between T2 and T4 starts to decrease, the resistivity remains constant until it reaches the T3-T1 section. It then goes to section T3-T1 to T1 and then T0. During the cooling phase, resistivity function does not follow the same pattern as the heating phase because of the hysteresis. It remains constant until reaching T3 = 67 ◦C. A reverse phase change occurs in the temperature range (T3 = 67 ◦C, T1 = 58 ◦C), bringing VO2 back to the nonconductive state. Resistivity moves up linearly on further reduction of temperature [41].

**Figure 6.** Hysteresis model in the custom simulator.

#### *3.3. Geometry of the Simulated Structure*

An equivalent geometry of our vertical structure in Figure 7 is built in the SUNRED simulator (*r* = 10 μm, *h* = 50 nm). The VO2 layer is sandwiched between platinum and silicon layers. The silicon substrate layer is 6.2 mm × 6.2 mm and 0.3 mm thickness. For faster simulations and simplicity, the quarter of the entire geometry is simulated. In this case, electric current must be multiplied by four while describing the full structure and the cavity length (*r*) represents the radius of the whole cavity. The bottom of the substrate (the bottom boundary condition) is a heat transfer coefficient (HTC) specified in Section 4. All other sides have a boundary condition of 10W/m2 K *HTC* (the value for air).

**Figure 7.** SUNRED model for VO2 vertical resistor.

#### *3.4. Excitation Required for the Simulation*

SUNRED capabilities allow different methods for creating excitation in the simulated structure. Electrical current, power (heat flow), voltage and temperature are the available options. For our study, we use electric current on the anode as excitation for the structure. The principal is to assign a different current value in each simulation step. Setting the cathode as the electric ground of the geometry (electric

potential = zero) initiates the direction of the electric current from anode to cathode. Optimal path for the current is to pass through VO2 and Si, reaching the cathode and creating the current flow because air and SiO2 have high resistivity (refer to Figure 3).

The simulation starts with a low excitation current and gradual increase in each step. Driving current warms up VO2 layer above the ambient temperature (50 ◦C) due to the generation of Joule-heat. Higher excitation current leads to a higher temperature of VO2 resistor and further elevation of temperature is required to simulate the phase-change of VO2. A similar excitation current profile is used for the second half of the simulation but with a falling direction instead of rising, Figure 8. The peak value and difference of excitation current in each step can be modified to improve the resolution. Smaller current steps give higher resolution.

**Figure 8.** Excitation current per simulation step.

#### **4. Results and Discussion**

Each flow of current accompanies an electric voltage developed between anode and cathode. (Voltage, Current) pair in each step is gathered to draw the V-I characteristic curve of our structure. For instance, the V-I curve for one of the performed simulations on vertical VO2 resistors is shown in Figure 9a. The simulated geometry is (*r* = 10 μm, *h* = 60 nm) and the ambient temperature is 50 ◦C. After analyzing several V-I characteristics, a general visualization of the V-I curve of vertical VO2 geometries can be formalized as illustrated in Figure 9b.

**Figure 9.** V-I curve of VO2 vertical structures: (**a**) simulation example, *r* = 10 μm, *h* = 60 nm, and ambient temperature is 50 ◦C and (**b**) the general shape for all other simulations.

The starting point of the study is to simulate the geometry in Figure 7 using the real size of produced sample in laboratory (*r* = 10 μm, *h* = 50 nm). The low-current section of the measured curve is strongly nonlinear, which did not occur in case of lateral SMT resistors [40,41]. The main difference between the current vertical structure and the previously studied lateral structures is that in lateral structures the platinum anode and cathode were directly attached to VO2, while in the vertical structure the silicon substrate was inserted between the cathode and VO2.

More than 90% of the voltage applied to the entire structure drops in the region between the cathode and the VO2 resistor even when the VO2 resistor is OFF. The three main parts of the region are the platinum-silicon transition at the cathode, the substrate, and the silicon-VO2 transition. Since the doping of the substrate is high, its resistivity is around 0.05 Ω cm, the voltage on it is only a fraction of the voltage applied to the region. In addition, the temperature-dependent nonlinearity of the substrate is very low [45], it cannot be responsible for the measured nonlinearity. The Pt-Si interface can show a Shottky like barrier which will depend on the Si dopant concentration [46]. Annealing could create a platinum-silicide transition that would give an ohmic contact, however, annealing cannot be applied due to VO2. Based on the band structure of VO2 [47], it is likely that a barrier will also be created at the junction of the two semiconductor materials. Most of the voltage in the region, therefore, drops at the two transitions.

Based on the measured results, the current dependence on resistivity of the transitions is mainly responsible for the nonlinearity, the temperature dependence is not significant, because at higher currents the curve is almost linear. Since the purpose is the investigation of size dependence of VO2 resistor at higher currents, the simulation is greatly simplified if nonlinearity is not taken into account and the transitions are modeled with constant resistivity. Accordingly, in the simulated model, a 100 nm thick resistor layer has been placed between both the Pt-Si and Si-VO2 interfacing surfaces.

Three cases were examined: (a) if a barrier is assumed only at the Pt-Si transition, (b) if a barrier is assumed only at the Si-VO2 transition and (c) if barriers are assumed at both transitions with same voltage drops. In case (a) the resistivity of the layer between Pt-Si is 1.33 × <sup>10</sup><sup>5</sup> <sup>Ω</sup> cm and the *HTC* at the bottom of the substrate is 46 W/m2 K; in case (b), the resistivity of the layer between Si-VO2 is 377 Ω cm and the *HTC* at the bottom of the substrate is 100W/m2 K; and in case (c), the resistivity of the layer between Pt-Si is 6.7 × <sup>10</sup><sup>4</sup> <sup>Ω</sup> cm, between Si-VO2 it is <sup>189</sup> <sup>Ω</sup> cm and the HTC is 60 W/m2 K. Results of the three simulations are shown in Figure 10b.

Both curves show an acceptable degree of validity. However, it is difficult to reach a near-absolute accuracy while neglecting the nonlinearity of the transitions. The location of the resistor defines the dissipation point. From circuitry point of view, the three cases are equivalent. We assume that there is a barrier at both transitions, so in the rest of the article, we perform the simulations with the model corresponding to case (c).

The experimental data in Figure 10a show a "noisy" behavior in the ascending current phase. When the current increases through the SMT layer, some conductive channels are switched ON, which as a consequense gives the noisy path. When the current decreases, the SMT conductive channel collapse faster. This phenomenon can also be observed with decreasing current, but it is much more moderated than with the increasing current case. The origin of the noisy signal is that a lot of parallel conductive channels are in the SMT layer (depends on the crystal structure). The real VO2 layer is not a homogenous crystal, it has a spiky structure [48]. While input power increases, some of the conductive channels are "switched ON", which are connected in parallel with the original SMT phase. If the dissipated power decreases (the resistance decreases) then it can switch back, which shows a jump back to the original path. This phenomenon also appears in other structures, for example in [40], the lateral structure with two Pt contacts are directly attached to VO2, no Si substrate is interposed.

**Figure 10.** V-I curve of (*r* = 10 μm, *h* = 50 nm) vertical VO2 resistor based on: (**a**) Laboratory measurement. (**b**) Simulation results.

#### *4.1. Dependence on Thickness*

In the first set of simulations, we investigated how the VO2 structure behaves on different thicknesses and what is the effect of changing *h* on the characteristics of our structure. The range of *h* used in the simulations is (10 nm to 300 nm) which is investigated at different radii. Figure 11 shows how the resistance (*R*) of VO2 layer changes with thickness. Resistance during the nonconductive and conductive phases is the point of interest here (refer to Figure 9b). Resistance is calculated as the 1/*slope* = (Δ*V*/Δ*I*) of the V-I curve. It can be concluded from Figure 11a that the nonconductive phase resistance increases as thickness increases. Furthermore, *R* dependence on *h* becomes lower at larger radii. On the other hand, conductive phase resistance shows no direct relation with thickness at these *h* and *r* ranges, Figure 11b.

**Figure 11.** VO2 resistance during (**a**) the nonconductive phase and (**b**) during the conductive phase at different thicknesses.

The current at which VO2 resistor is turned on (when resistivity drops significantly) is called the opening current. Similarly, the current at which the VO2 resistor is turned off (during the cooling phase) is called the closing current. Voltages associated with these currents are opening voltage and closing voltage, respectively. Figure 12 illustrates the dependence of opening and closing currents on the variation of thickness. Since the opening current is associated with the nonconductive phase resistance and the closing current is associated with the conductive phase resistance, the opening current falls as

thickness increases due to the elevation of *R* versus thickness (in the nonconductive phase), Figure 12a. On the other hand, the current remains unchanged as thickness increases because resistance is not affected by thickness (in the conductive phase), Figure 12b. A high degree of compatibility can be noticed between Figures 11 and 12.

**Figure 12.** The dependence of the (**a**) opening and (**b**) closing currents on the thickness of the VO2 resistor.

#### *4.2. Dependence on Radius*

The second set of simulations is to study how the structure responds to change in radius. Simulations performed with different radii *r* in the range of (200 nm to 30 μm) and are investigated at different thicknesses. Figure 13 concludes that radius has a strong effect on the VO2 resistance, and the resistance drops nonlinearly for both the nonconductive phase (Figure 13a) and the conductive phase (Figure 13b) as radius increases. It is explained with the relation *R* = *ρ <sup>l</sup> <sup>A</sup>* , *ρ* is the resistivity, *l* is length (which is *h* here) and *A* is the cross-sectional area (which is here *πr*2).

**Figure 13.** VO2 resistance during (**a**) the nonconductive phase and (**b**) during the conductive phase at different radii.

As a consequence of the resistance dependence on radius, opening and closing currents increase as radius increases. Figure 14 illustrates the effect of changing the radius on the opening current (Figure 14a) and the closing current (Figure 14b). For low current applications, it is worth to consider decreasing the radius since it allows the phase change to occur at much lower currents.

**Figure 14.** The dependence of the (**a**) opening and (**b**) closing currents on the radius of the VO2 resistor.

During the hysteresis domain, the resulting V-I curve also has hysteric nature (refer to Figure 9b). Maximum difference between the resistance of VO2 resistor during heating and cooling phases is called the maximum switching resistance. It can be observed that the maximum switching resistance varies linearly versus thickness Figure 15a, in contrast, it has a nonlinear dependence on the radius Figure 15b.

**Figure 15.** The dependence of the maximum switching resistance on (**a**) the thickness and (**b**) the radius of the VO2 resistor.

#### **5. Conclusions**

The results of this article offer a good base of sample sizes estimates for further studies in submicron modeling of VO2 resistors. A valid operating model will pave the way for designing thermal-electronic logic circuits. The results provided in our previous work [39] about the lateral structures and results of this article about the vertical structures will straiten the range of samples sizes for producing and studying. This will save both time and effort.

**Author Contributions:** Conceptualization, P.N. and M.D.; methodology, M.D.; software, L.P.; validation, M.D.; technology, J.M.; sample creation and measurement, P.N.; writing—original draft preparation, M.D., L.P., P.N.; writing—review and editing, L.P., M.D.; visualization, M.D.; supervision, L.P.; project administration, J.M.; funding acquisition, P.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research reported in this paper and carried out at the Budapest University of Technology and Economics was supported by the "TKP2020, National Challenges Program" of the National Research Development and Innovation Office (BME NC TKP2020). This research also received fund by the Higher Education Excellence Program of the Ministry of Human Capacities in the frame of Nanotechnology research area of the Budapest University of Technology and Economics (BME FIKP-NANO) and the Science Excellence Program at BME under the grant agreement NKFIH-849-8/2019 of the Hungarian National Research, Development and Innovation Office is also acknowledged. The research reported in this paper was partially supported by the Stipendium Hungaricum Scholarship Programme of the Hungarian Government.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
