Circuit-Based Electrothermal Simulation of Multicellular SiC Power MOSFETs Using FANTASTIC

Abstract

This paper discusses the benefits of an advanced highly-efficient approach to static and dynamic electrothermal simulations of multicellular silicon carbide (SiC) power MOSFETs. The strategy is based on a fully circuital representation of the device, which is discretized into an assigned number of individual cells, high enough to analyze temperature and current nonuniformities over the active area. The cells are described with subcircuits implementing a simple transistor model that accounts for the utmost influence of the traps at the SiC/SiO₂ interface. The power-temperature feedback is emulated with an equivalent network corresponding to a compact thermal model automatically generated by the FANTASTIC tool from an accurate 3D mesh of the component under test. The resulting macrocircuit can be solved by any SPICE-like simulation program with low computational burden and rare occurrence of convergence issues.

1. Introduction

Silicon carbide (SiC) power devices are promising candidates for energy distribution, as well as for automotive, aircraft, and spacecraft applications, by virtue of their inherent features like high breakdown voltage, low on-state resistance, and excellent high-temperature capability [1].

Unfortunately, such devices often operate under critical conditions with a large amount of heat generation, which may lead to reliability degradation or even to an irreversible device failure in harsh cases. As a consequence, reliable simulation tools accounting for electrothermal (ET) effects are highly desired to define the thermal dissipation constraints and optimize the design of the transistor layout and/or of the cooling system. Such tools must be suited to describe temperature and current nonuniformities, which are often responsible for the safe operating area shrinking of transistors with a multicellular pattern. However, conceiving and developing a viable simulation strategy are challenging tasks due to multiple reasons. Fully numerical 3D ET analyses with device simulators concurrently solving the semiconductor and heat transfer equations are computationally unfeasible. Commonly-adopted approaches rely (i) on the interaction between a circuit simulation program and a 3D thermal-only numerical solver in a relaxation procedure [2,3], or (ii) on the extension of a finite-volume/-element software package to account for the electrical behavior of the transistor with simplified models [4,5]. However, for the specific case of SiC power devices, results can be trustworthy only by using models that accurately describe the key physical parameters and their non-intuitive temperature dependences, which are rather different compared to the traditional silicon (Si) counterparts. In addition, regardless of the technology, the pre-processing geometry/mesh construction within the environment of the thermal solver is onerous and troublesome. Lastly, these approaches are very resource-hungry and prone to convergence failures, especially if dynamic simulations under critical conditions have to be performed.

In this paper, an innovative circuit-based ET simulation approach is proposed for multicellular SiC power MOSFETs, with the ambition of optimizing the trade-off between computational efficiency and accuracy. The strategy is articulated as follows: (i) the device is discretized into an assigned number N of individual cells (each associated to an independent heat source) described with a simple, yet accurate model accounting for the relevant influence of SiC/SiO₂ interface traps; (ii) the cell model is implemented with a SPICE-compatible subcircuit; (iii) an exceptionally accurate 3D finite-element method (FEM) description of the device is effortlessly obtained through a commercial solver aided by an in-house routine; (iv) the FANTASTIC code [6] is invoked, which automatically generates a dynamic compact thermal model (DCTM) of the device from the FEM representation, and thus builds an electrical network emulating the power-temperature feedback; (v) such a network is enriched with suitable voltage sources to account for nonlinear thermal effects, and the resulting circuit is referred to as thermal feedback block; (vi) a purely-electrical macrocircuit describing the whole ET behavior of the power device is constructed by connecting the N cell subcircuits to the thermal feedback block; (vii) the macrocircuit can be solved by any commercial circuit simulator in very short times and with unlikely occurrence of convergence problems. A multicellular 4H-SiC power MOSFET soldered on a direct bonded copper (DBC) substrate and operated under dc, short-circuit (SC), and unclamped inductive switching (UIS) conditions is considered as a case study.

This work extends the preliminary contribution [7], where a fully circuital representation of the SiC power MOSFET was obtained and solved with a similar approach. However, the generation of the thermal feedback block, based on Foster networks, was carried out in a long pre-processing stage involving N 3D FEM transient simulations, as well as a procedure to determine the proper number of RC pairs of each network and identify their values; in addition, the thermal interactions among horizontally-far cells were either coarsely described or even disregarded.

The remainder of the paper is articulated as follows. In Section 2, the device selected to test the approach and the experimental setup are described. Section 3 offers details concerning the transistor model used for the elementary cells. Section 4 probes into the ET simulation approach. Results are reported and discussed in Section 5. Conclusions are then drawn in Section 6.

2. Device under Test and Experimental Setup

The whole analysis was performed on the CREE 4H-SiC vertical double-diffused MOSFET (VDMOS) denoted as CPMF-1200-S080B, rated 1200 V, 50 A, 80 mΩ, and targeted at solar inverters, high-voltage dc–dc converters, and motor drives. A top-view picture of the device under test (DUT) is given in Figure 1; indicated are the gate pad, the two source pads, and the gate interconnect tracks. The DUT presents a 4.08 × 4.08 mm² die area and a 3.46 × 3.46 mm² active area (there is a peripheral inactive region), the effective portion of which contributing to the current capability amounts to about 10 mm². The pattern is multicellular, with many thousands of body islands located in the N-drift region, within which there are source-body contacts surrounded by the polysilicon gate that lies beneath the source metal. The Ni/Ag drain contact is on the backside.

Isothermal measurements of I–V transfer and output characteristics of the bare die were performed by means of an in-house 250 A-rated curve tracer suited to apply down to 1 µs-wide current pulses, the device baseplate being set to assigned temperatures T_B through a thermochuck with 1 °C resolution. UIS experiments aimed at the evaluation of the breakdown voltage were carried out by a non-destructive custom tester [8]. The switching behavior was investigated by using a half-bridge converter board configured to perform a standard inductive load switching (ILS) test [9].

3. Transistor Model and Parameter Extraction Methodology

In the last decades, a noticeable effort has been made to develop models for SiC MOSFETs, a review of them being offered in [10]. Here we propose a slightly modified version of the behavioral model presented in [7,9] and used for ET simulations in [7,9,11]. Such a model, unlike those hitherto described in the literature, concurrently enjoys the following benefits: (i) it is simple, yet accurate enough, with a few parameters easy to extract; (ii) it can be implemented with a subcircuit compatible with any SPICE-like program; (iii) it includes all the key physical parameters and their specific temperature dependence up to very high temperatures; (iv) it also accounts for avalanche effects due to impact ionization (II); (v) the nonlinear nature of the intrinsic capacitances can be activated.

3.1. Transistor Model

The model is an enriched variant of the classic SPICE Level 1. In particular, the transistor is represented as the series of an “intrinsic” conventional MOSFET describing the channel behavior, and a resistance for the lowly-doped N-type epitaxial region, as depicted in Figure 2.

Let us consider the following nomenclature:

• T [K] is the transistor temperature, assumed uniform in the regions impacting the device behavior.
• T₀ = 300 K is the reference temperature.
• ΔT = T − T₀ [K] is the temperature rise over T₀.
• T_B [K] the baseplate (and thermochuck) temperature.
• V_GS [V] is the gate-source voltage.
• V_DS [V] is the drain-source voltage.
• I_D [A] is the drain current.
• R_drift [Ω] is the bias- and temperature-dependent resistance of the N-type epitaxial drift region.
• V_DSch = V_DS − V_drift [V] is the V_DS portion falling over the channel (V_drift = R_drift·I_D).
• V_TH [V] is the temperature-dependent threshold voltage.
• K [A/V²] is the temperature-dependent current factor.

The channel region is described with the Level 1 model. If V_DSch is lower than the overdrive voltage V_GS − V_TH, the DUT operates in triode mode, and the II-free drain current I_DnoII is expressed as

$$I_{DnoII}=K\cdot \left[2\cdot \left(V_{GS}-V_{TH}\right)\cdot V_{DSch}-V_{DSch}^2\right]$$

(1)

Conversely, if V_DSch ≥ V_GS − V_TH, then the DUT is driven into pinch-off, and

$$I_{DnoII}=K\cdot {\left(V_{GS}-V_{TH}\right)}^2$$

(2)

The negative temperature coefficient of V_TH is described through the following law:

$$V_{TH}\left(T\right)=\left[V_{TH}\left(T_0\right)-V_{TH\infty }\right]\cdot \exp \left(-a_{V_{TH}}\cdot \mathsf{\Delta }T\right)+V_{TH\infty }$$

(3)

such an exponential model being an improvement with respect to the simple linear relationship used in [7].

The current factor K depends upon T since the electron mobility in the channel is temperature-sensitive; similar to [7,9], such a dependence is taken into account through the power relationship

$$K\left(T\right)=K\left(T_0\right)\cdot {\left(\frac{T}{T_0}\right)}^{-m\left(T\right)}$$

(4)

where the exponent m(T) is

$$m\left(T\right)=-a_m+\left(a_m+b_m\right)\cdot \left[1-c_m\cdot \exp \left(-d_m\cdot \frac{T}{T_0}\right)\right]$$

(5)

II effects are accounted for as follows [7,9]. The bias- and temperature-sensitive avalanche multiplication factor M (≥1) is given by [12]

$$M\left(V_{DS},I_D,T\right)=1+m_{II}\cdot \tan \left\{f_I\left(I_D\right)\cdot \frac{\pi }{2}\cdot {\left[\frac{V_{DS}-R_{II}\cdot I_D}{B{V}_{DS}\left(T\right)}\right]}^{n_{II}}\right\}$$

(6)

where BV_DS(T) is the temperature-dependent drain-source breakdown voltage, expressed as

$$B{V}_{DS}\left(T\right)=B{V}_{DS}\left(T_0\right)\cdot \exp \left({\alpha }_{II}\cdot \mathsf{\Delta }T\right)$$

(7)

and f_I(I_D) is a nondimensional correction term to describe a potential II dependence on current (i.e., on biasing conditions), given by

$$f_I\left(I_D\right)=\exp \left({\beta }_{II}\cdot I_D\right)$$

(8)

Let us introduce the avalanche coefficient ξ = M − 1 (≥0). The II-affected drain current I_D is evaluated as

$$I_D=I_{DnoII}+I_{DII}=I_{DnoII}+\xi \cdot \left(I_{leak}+I_{DnoII}\right)$$

(9)

where I_DII is the additional current component only dictated by II, and I_leak is a small leakage current.

The resistance R_drift is expressed as the sum of (i) a bias- and temperature-dependent resistance R_JFET to model the path composed by the accumulation and JFET regions, and (ii) a temperature-dependent resistance R_epi for the epitaxial region beneath the JFET one [7,9]:

$$R_{drift}\left(V_{GS},V_{drift},T\right)=R_{JFET}\left(V_{GS},V_{drift},T\right)+R_{epi}\left(T\right)$$

(10)

where

$$\begin{array}{l}{R}_{JFET}\left(V_{GS},V_{drift},T\right)=R_{JFET}\left(T_0\right)\cdot {\left(\frac{T}{T_0}\right)}^{m_{R_{JFET}}}\cdot \left(\frac{1}{1+V_1/V_{drift}}\right)\cdot {\left(\frac{V_{GS}}{V_2}\right)}^{-\eta }\\ R_{epi}\left(T\right)=R_{epi}\left(T_0\right)\cdot {\left(\frac{T}{T_0}\right)}^{m_{R_{epi}}}\end{array}$$

(11)

R_JFET(T₀) being the JFET resistance at T = T₀, V_drift » V₁, and V_GS = V₂ (V₁ and V₂ are fitting parameters). This formulation improves the one reported in [7] in the high-current triode region and is derived on the basis of simple arguments. First, the resistance of the accumulation region reduces with gate voltage due to the increased concentration of the attracted electrons; second, under high V_drift values, the high electric field occurring in the JFET region tends to saturate the electron velocity, thus degrading the mobility.

The dynamic transistor behavior is described by improving the Level 1 capacitance models; the nonlinear nature of C_GD and C_DS = C_DB is accounted for with the following expressions [9]:

$$C_{GD}\left(V_{GD}\right)=\left(C_{GD0}-C_{GDMIN}\right)\cdot \left[1+\frac{2}{\pi }\arctan \left(\frac{V_{GD}}{V^{*}}\right)\right]$$

(12)

and

$$C_{DS}\left(V_{DS}\right)=\frac{2}{\pi }\cdot C_{DS0}\cdot \left[\frac{\pi }{2}+\arctan \left(-\frac{V_{DS}}{V^{**}}\right)\right]+C_{DSMIN}$$

(13)

while C_GS was not modified [13].

3.2. Parameter Extraction Procedure

The threshold voltage V_TH and the current factor K were extracted in a very wide T_B range spanning from 300 to 500 K by resorting to the traditional quadratic extrapolation method [14,15] applied to I_D–V_GS transfer characteristics measured under isothermal (pulsed) conditions at V_DS = 20 V with the curve tracer mentioned in Section 2. Then the parameters in (3), (4), and (5) were calibrated so as to ensure the best matching between experimental data and the following relations:

$$V_{TH}\left(T_B\right)=\left[V_{TH}\left(T_0\right)-V_{TH\infty }\right]\cdot \exp \left[-a_{V_{TH}}\cdot \left(T_B-T_0\right)\right]+V_{TH\infty }$$

(14)

$$K\left(T_B\right)=K\left(T_0\right)\cdot {\left(\frac{T_B}{T_0}\right)}^{-m\left(T_B\right)}$$

(15)

with

$$m\left(T_B\right)=-a_m+\left(a_m+b_m\right)\cdot \left[1-c_m\cdot \exp \left(-d_m\cdot \frac{T_B}{T_0}\right)\right]$$

(16)

The comparison between the measured V_TH and the optimized (14) is shown in Figure 3. It can be inferred that the DUT exhibits (i) a high V_TH(T₀) (≈6.4 V) and (ii) a high negative temperature coefficient of V_TH(T) at low/medium T_B compared to similarly-rated Si power MOSFETs. Both findings were attributed to the high density of SiC/SiO₂ interface traps (quantum states originating from the thermal oxidation of the SiC surface); more specifically, (i) electrons are captured by traps and do not contribute to channel formation, thereby leading to a high threshold voltage V_TH(T₀), and (ii) V_TH markedly reduces with temperature due to the concurrent effect of more broken bounds that release electrons, and of the emission of inversion electrons from the traps, the latter effect being almost absent in the Si counterparts [16,17,18]. It must be underlined that the strong negative temperature coefficient in turn contributes to a significant positive temperature coefficient of I_D, which may exacerbate the ET feedback [19,20].

The comparison between the measured K and model (15) and (16) with tuned parameters is shown in Figure 4. The temperature sensitivity of K is only related to that of the channel electron mobility μ_n, which is due to the interplay between (i) the Coulomb scattering with the filled (charged) interface traps, leading to a positive temperature coefficient induced by the trap discharging (release of electrons) with increasing temperature, and (ii) the acoustic-phonon scattering favoring a negative coefficient, where (i) prevails at low temperatures and (ii) dominates at high temperatures [18,21,22,23]. This behavior is accurately described by the m model (16).

The accuracy of the parameter calibration for the V_TH and K models is witnessed by the comparison reported in Figure 5 between (2) and the I_D–V_GS transfer characteristics measured under isothermal conditions at V_DS = 20 V and various T_B values in a current range wherein the DUT operates in pinch-off (note that I_D ≈ I_DnoII). An inspection of the curves reveals the considerable positive temperature coefficient of I_D within a wide range of currents.

The parameters of the drift resistance R_drift in (10), (11) were optimized by comparing the I_D–V_DS output characteristics at various V_GS and different T_B values in triode region with the model (1) applied to the scheme in Figure 2. Figure 6 reveals the accuracy of the extraction at T_B = 303 K.

The parameters of the II model given by (6) with (7) and (8) were tailored on the basis of experimental data determined under UIS conditions, and 2D numerical simulations of the T_B-dependent I_D–V_DS avalanche curve at V_GS = 0 V of an individual semi-cell of the DUT carried out with TCAD Sentaurus [24] (it must be remarked that the term “cell” here does not correspond to the discretization adopted in this work and identified by N, but it is one of the tens of thousands of identical blocks in which the effective active area of the MOSFET can be partitioned).

The parameters used in the expression of the capacitances C_GD (12) and C_DS (13) were calibrated to match the experimental gate and drain waveforms during an ILS turn-off and turn-on transient at different supply voltages.

All the optimized parameter values are reported in

Table 1. Optimized parameter values used in the transistor model.

Parameter	Definition	Value
VTH(T0)	threshold voltage at reference temperature T0	6.398 V
VTH∞	parameter of the threshold voltage model	2.05 V
aVTH	temperature coefficient of the threshold voltage	6 mK−1
K(T0)	current factor at reference temperature T0	0.422 A/V2
am	parameter of the exponent m accounting for the mobility temperature dependence	0.24
bm	parameter of the exponent m accounting for the mobility temperature dependence	2
cm	parameter of the exponent m accounting for the mobility temperature dependence	1.02
dm	parameter of the exponent m accounting for the mobility temperature dependence	0.09
RJFET(T0)	resistance of the accumulation and JFET regions at reference temperature T0, Vdrift » V1, and VGS = V2	0.235 Ω
mRJFET	exponent for the temperature dependence of the resistance RJFET	−1.3
V1	parameter to account for the RJFET dependence on Vdrift	13 V
V2	parameter to account for the RJFET dependence on VGS	20 V
η	exponent to account for the RJFET dependence on VGS	3.45
Repi(T0)	resistance of the deep epilayer region at reference temperature T0	10 mΩ
mRepi	exponent to account for the potential temperature dependence of the resistance Repi	0
BVDS(T0)	drain-source breakdown voltage at reference temperature T0	1750 V
mII	parameter of the avalanche multiplication factor M	1.8
nII	parameter of the avalanche multiplication factor M	2.9
αII	temperature coefficient of the breakdown voltage BV	0.18 mK-1
βII	coefficient to account for the ID dependence of factor M	0 A−1
RII	resistance to account for the ID dependence of factor M	10 Ω
CGD0	zero-bias gate-drain capacitance	0.85 nF
CGDMIN	minimum gate-drain reverse-biased capacitance	0.01 nF
V *	gate-drain capacitance parameter	2 V
CDS0	zero-bias drain-source capacitance	2.8 nF
CDSMIN	minimum drain-source reverse-biased capacitance	0.06 nF
V **	drain-source capacitance parameter	10 V

.

3.3. Impact of Interface Traps

The impact of the interface traps was investigated by resorting to the following strategy. The traps were virtually removed by (i) reducing V_TH(T₀) to 4 V, (ii) decreasing a_VTH to 2 mK⁻¹, (iii) multiplying the current factor K(T₀) by 50 to annihilate the degradation of mobility μ_n(T₀), and (iv) setting c_m = 0 to eliminate the influence of Coulomb scattering on μ_n, which will thus exhibit a negative temperature coefficient over the whole temperature range. Figure 7 shows the comparison between the real 4H-SiC DUT and the ideal traps-free counterpart in terms of transfer characteristics at V_DS = 20 V and various T_B values. It can be inferred that:

• differently from the DUT, which suffers from a marked positive temperature coefficient over the entire current range, the traps-free device is subject to a slight positive coefficient only at low drain current I_D (<20 A), where the negative temperature coefficient of V_TH prevails over the negative coefficient of μ_n, while beyond a zero-temperature coefficient region (also referred to as compensation region), the negative coefficient of μ_n dominates, and I_D reduces with temperature.
• the traps-free transistor benefits from a much higher current capability due to the lower V_TH and the higher μ_n.

To further corroborate the above findings, Figure 8 reports the temperature coefficient of the drain current I_D, given by [19,20,25]

$${\alpha }_T={\frac{\partial I_D}{\partial T}|}_{V_{DS}}$$

(17)

for the real DUT and its traps-free variant at reference temperature T₀ and drain-source voltage V_DS = 20 V; it is again witnessed that the traps at the SiC/SiO₂ interface lead to a detrimental highly-positive α_T within a broad range of currents.

4. Electrothermal Simulation Approach

We chose to resort to a circuit-based approach [7,11,26,27], which is a good trade-off between computational burden and accuracy. Such an approach makes use of the thermal equivalent of the Ohm’s law (TEOL) and can be summarized as follows.

• The whole DUT is subdivided into an assigned number N of elementary cells, high enough to identify potentially-dangerous temperature gradients over the transistor active area, but not too high to prevent intolerably long CPU times or impossible memory storage. The individual cells are described with the model explained in Section 3, where the area-dependent parameters are properly scaled. For the simulation campaign reported in Section 5, all cells are assumed identical, although it is in principle possible to assign different parameter values to different cells to allow a statistical analysis of the influence of parameter fluctuations on the ET behavior of the component.
• Each cell is represented with a SPICE-compatible subcircuit implementing the above model through a macromodeling technique. The subcircuit makes use of (i) a standard MOSFET at reference temperature T₀ = 300 K as a “main” component to describe the channel region, and (ii) linear and nonlinear controlled sources to include all the model features that cannot be accounted for with the basic MOSFET, i.e., the temperature dependence of the threshold voltage V_TH and of the current factor K, the bias- and temperature-dependent drift resistance R_drift, as well as the bias- and temperature-dependent II mechanism. The TEOL is adopted, namely, the temperature rise over ambient ΔT = T − T₀ is actually a voltage, while the dissipated power P_D is treated as a current; this allows (i) enabling the temperature sensitivity of the key physical parameters, and (ii) describing the power-temperature feedback with an electrical network. Besides the standard electrical terminals (gate, drain, source/body), the cell subcircuit is also equipped with an input node carrying the “voltage” ΔT (provided by the thermal feedback block introduced below) and with an output node offering the “current” P_D (to be fed to the thermal feedback block).
• The power-temperature feedback (i.e., the dynamic heat propagation within the structure) is described with a TEOL-based SPICE-compatible thermal feedback block, which is composed by resistances, capacitances, and controlled sources. The inputs of this block are the powers P_D dissipated by the transistor cells (represented with currents), and the outcomes are the individual (nonlinear) temperature rises ΔT (emulated with voltages).
• The thermal feedback block contains an equivalent thermal network (i.e., a purely-electrical circuit relying on the TEOL). The main contribution of the proposed approach is that this network is automatically constructed in the pre-processing stage by invoking the FANTASTIC tool [6] (Section 4.3). FANTASTIC receives as an input an accurate 3D FEM representation of the domain, i.e., a mesh with information about (i) the discretization into elementary cells (each corresponding to an individual heat source), (ii) material parameters, and (iii) boundary conditions, and then extracts a reduced-order model and the associated network without the need of user’s expertise and COMSOL simulations; only 16 min were needed for the case study. The equivalent network accurately accounts for the self-heating of each cell and for the mutual interactions among all cells and describes the linear thermal problem. However, nonlinear thermal effects can be significant if the DUT is simulated under harsh conditions entailing high temperatures. In order to tackle this issue, a properly-tuned Kirchhoff’s transformation [28] is used, which converts the linear temperature rises (ΔT_lin) into their nonlinear counterparts (ΔT) through [29]

  $$\mathsf{\Delta }T=T_0\cdot {\left[m_k+\left(1-m_k\right)\cdot \frac{\mathsf{\Delta }{T}_{lin}+T_0}{T_0}\right]}^{\frac{1}{1-m_k}}-T_0$$

  (18)
  The calibration of the (positive) parameter m_k will be detailed in Section 4.2. Nonlinear voltage-controlled voltage sources emulating (18) are applied to the N temperature rise nodes of the equivalent thermal network; the resulting circuit is referred to as a thermal feedback block. It is worth noting that the FANTASTIC-based approach improves the strategy exploited in [7,26], where the thermal feedback block was based on Foster networks extracted in a rather long pre-processing stage from N onerous transient COMSOL simulations of the DUT (performed by activating one heat source at a time), and the thermal coupling between horizontally-far heat sources was roughly described or even neglected.
• The cell subcircuits are then connected to the thermal feedback block in a commercial circuit simulation tool (like PSPICE, LTSPICE, Eldo, ADS, SIMetrix); as a result, the whole domain under test, composed by the DUT soldered on a DBC substrate, is transformed into a purely-electrical macrocircuit, which suitably accounts for ET effects: the temperature, and thus the temperature-sensitive parameters, are allowed to vary during the simulation run. The solution of this macrocircuit under both static and dynamic conditions is demanded to the powerful and robust engine of the circuit simulation tool, with very low computational effort and minimized occurrence of convergence issues compared to other numerical methods.

4.1. SPICE Subcircuit and DUT Discretization into Cells

A sketch of the SPICE-compatible subcircuit for the transistor cell is represented in Figure 9, where the standard MOSFET instance at reference temperature T₀ is indicated with M_n. The additional input node feeding the “voltage” ΔT and the output node providing the “current” P_D are also represented.

The negative temperature coefficient of the threshold voltage V_TH(T) described by (3) is accounted for by using the following approach. The voltage source A in series with the gate adds V_TH(T₀) − V_TH(T) to V_G, V_TH(T) being given by (3), thus biasing M_n with V_G’ = V_G + [V_TH(T₀) − V_TH(T)]; as a result, the effective overdrive voltage becomes

$$V_{G’S}-V_{TH}\left(T_0\right)=V_{GS}+\left[V_{TH}\left(T_0\right)-V_{TH}\left(T\right)\right]-V_{TH}\left(T_0\right)=V_{GS}-V_{TH}\left(T\right)$$

(19)

and M_n conducts the current I_D(M_n) given by

$$\begin{array}{l}{I}_D\left(M_n\right)=K\left(T_0\right)\cdot \left\{2\left[V_{GS}-V_{TH}\left(T\right)\right]V_{DSch}-V_{DSch}^2\right\} V_{DSch}<V_{GS}-V_{TH}\left(T\right)\\ I_D\left(M_n\right)=K\left(T_0\right)\cdot {\left[V_{GS}-V_{TH}\left(T\right)\right]}^{2 } V_{DSch}\ge V_{GS}-V_{TH}\left(T\right)\end{array}$$

(20)

The temperature dependence of the current factor K(T) expressed by (4) with (5) is emulated by using the current source B to derive the current

$$I_{\mu }=I_D\left(M_n\right)\cdot \left[{\left(\frac{T}{T_0}\right)}^{-m\left(T\right)}-1\right]$$

(21)

Consequently, the II-unaffected current I_DnoII is obtained as

$$I_{DnoII}=I_D\left(M_n\right)+I_{\mu }=I_D\left(M_n\right)+I_D\left(M_n\right)\cdot \left[{\left(\frac{T}{T_0}\right)}^{-m\left(T\right)}-1\right]=I_D\left(M_n\right)\cdot {\left(\frac{T}{T_0}\right)}^{-m\left(T\right)}$$

(22)

where m(T) is given by (5).

II effects are activated by adding an avalanche current I_DII = ξ·(I_leak + I_DnoII) to I_DnoII through the current source C, ξ = M − 1 being bias- and temperature-dependent according to (6). Hence, the drain current I_D flowing through the drift resistance R_drift is given by (9).

The bias- and temperature-dependent R_drift described by (10) with (11) is taken into account by making use of source D, which imposes the voltage drop

$$V_{drift}=I_D\cdot R_{drift}\left(V_{GS},V_{drift},T\right).$$

(23)

Only half of the DUT was represented and simulated by exploiting its inherent symmetry. The effective active region (≈5 mm²) was partitioned into N = 79 cells, each with a 250 × 250 µm² area (the procedure leading to the choice of the N value will be detailed in Section 4.4). As a consequence, compared to those reported in Table 1, the values of the area-dependent parameters were scaled by dividing the current factor and the capacitances by 2 × N, and multiplying the resistances by 2 × N. The popular commercially-available OrCAD PSPICE [30] was chosen to perform circuit simulations. The main schematic (resembling the actual layout) is represented in Figure 10, along with the corresponding cell numbering and a detail of the connections between adjacent cells. The drain current of the resulting macrocircuit (an outcome of the ET simulation) has to be multiplied by 2.

4.2. FEM Representation of the Component under Test

Exhaustive details on the geometry and materials of the DUT (the CPMF-1200-S080B VDMOS) were taken from the datasheet and from reports available on a reverse-engineering website. The DUT was assumed to be soldered on a DBC substrate by means of a 50 µm-thick tin-platinum alloy (SnPt) layer ensuring both mechanical joint and electrical connection with the drain [31]. Figure 11 and Figure 12 show the top view and cross-section of the assembly, which enjoys the same symmetry of the DUT. The soldering process can be automated by means of thin SnPt films, which are preformed and match the size of the die; in addition, alignment masks are typically exploited to ease their positioning. The DBC is composed by two copper (Cu) sheets (namely, bottom and top plates) with a ceramic layer placed in between them. Such a layer provides dielectric insulation to the assembly and is realized with alumina (Al₂O₃) for the proposed case study; however, alternative materials like silicon nitride (Si₃N₄) [32] and aluminum nitride (AlN) [33] are also adopted in the DBC manufacturing process. The drain contact of the DUT is soldered on the top plate, while the bottom plate ensures a good thermal interface with the 3 mm-thick Cu baseplate.

The 3D domain was represented in the environment of the commercial FEM-based COMSOL Multiphysics software package [34] by means of the in-house routine detailed in [35], which allows automatically building an extremely accurate geometry (Figure 13) and performing a smart selective optimization of the tetrahedral mesh (Figure 14); this process conveniently avoids a painstakingly long and prone-to-errors manual procedure for geometry/mesh construction. The number of elements (tetrahedra) and degrees of freedom (DoFs) of the resulting grid are 3.8 × 10⁵ and 5.2 × 10⁵, respectively. Figure 14 plainly illustrates that the mesh is highly fine over the die, while becoming gradually coarser by moving far away from the active region.

As previously mentioned, all transistor cells (Figure 10) were individually associated to heat sources located over the top surface of the 4H-SiC DUT. In this representation, the heat is assumed to be dissipated over the whole effective active area, and not over the tens of thousands of individual channels; conveniently, this unavoidable approximation was demonstrated to be reasonable under dc and transient conditions (at least for not-too-short times) through a simulation analysis performed with TCAD Sentaurus and COMSOL. An isothermal boundary condition at T_B = T₀ was applied at the bottom of the assembly baseplate, whereas all other surfaces were assumed adiabatic (i.e., with zero outgoing heat flux). The material parameters adopted for the thermal simulations are listed in

Table 2. Properties of the materials composing the assembly.

Material	k(T0) (W/mK)	cp (J/KgK)	ρ (Kg/m3)	α	β (W/mK2)
4H-SiC	370 [36,37]	690 [37]	3211 [38,39] (value common to 4H-SiC and 6H-SiC)	1.29 [40]
Al	240 [41]	905 [41]	2707 [41]		0.04 [41]
SnPt	68.8 [41]	228 [41]	7310 [41]		0.02 [41]
Ni	89.5 [41]	445 [41]	8906 [41]		0.08 [41]
Ag	427 [41]	236 [41]	10,524 [41]		0.07 [41]
poly-Si	40 [39]	920 [39]	2330 [39]
SiO2	1.38 [39]	709 [39]	2203 [39]	−0.33 [39]
Al2O3	28 [39]	796 [39]	3900 [39]	1 [39]
Cu	396.8 [41]	384 [41]	8954 [41]		0.05 [41]
Si3N4	18.5 [39]	787 [39]	3100 [39]	−0.33 [39]

. Nonlinear thermal effects were accounted for by including the temperature dependences of the thermal conductivities described by

$$k(T)=k(T_0)\cdot {\left(\frac{T}{T_0}\right)}^{-\alpha }$$

(24)

$$k(T)=k(T_0)-\beta \cdot (T-T_0)$$

(25)

It must be noted that (24) applies to semiconductors and insulators, while (25) is followed by some metals.

The pre-processing calibration of parameter m_k to be used for the Kirchhoff’s transformation (18) was carried out with the following procedure. The domain under test was thermally simulated with COMSOL by varying the power P_D dissipated by the whole effective active area (described with only one heat source) over a wide range (0.1 to 573.5 W for half device, the latter value leading to a peak temperature of 1400 K); the thermal conductivity dependences upon temperature were either activated (nonlinear conditions) or deactivated (linear conditions). The average temperature rise over that area was determined for both cases. Then the Kirchhoff’s transformation was applied to the linear temperature rise, and m_k was adjusted so as to obtain a good agreement between the nonlinear temperature rise computed by the transformation and the realistic one evaluated by COMSOL; it was obtained that m_k = 0.785. The whole process lasted less than 1 h, as each nonlinear simulation was run in about 2 mins and 24 P_D values were applied.

4.3. FANTASTIC-Based Derivation of the Equivalent Thermal Network

The FANTASTIC tool, originally introduced by some of the authors in [6], where it was applied to merely-thermal 3D simulations of a state-of-the-art gallium-arsenide (GaAs) HBT, is based on the truncated balance (TRB)-based moment matching (MM) approach to model-order reduction (MOR) developed by one of the authors in [42].

In this tool, the dynamic heat conduction problem in a semiconductor device, assumed to be linear, is imported from either commercial (e.g., COMSOL) or open-source codes (e.g., SALOME SMESH), comprising: the mesh discretizing the geometry, as well as the definitions of heat sources, materials, and boundary conditions. Both hexahedral and tetrahedral meshes can be used. It is possible to define arbitrary heat capacity and tensorial thermal conductivity distributions. Neumann’s, Dirichlet’s, or Robin’s boundary conditions can be applied. Superficial (i.e., indefinitely thin) and volumetric heat sources can be taken into account.

The FEM model of the thermal problem is then assembled by FANTASTIC. In particular, the mass matrix M and the stiffness matrix K are constructed. High-order basis functions can be used: most commonly, as a trade-off between efficiency and accuracy, tetrahedral meshes and second-order basis functions are considered. The M DoFs of the temperature rise distribution, forming the M-row vector $\mathbf{\vartheta }(t)$, are solutions of the discretized linear dynamic heat conduction problem

$$M\frac{d\mathbf{\vartheta }}{dt}(t)+K\mathbf{\vartheta }(t)=q(t)$$

(26)

in which the power density distribution vector q(t) takes the form

$$q(t)=QP(t)$$

(27)

where P(t) is an N-row vector with the powers dissipated by N independent heat sources, and Q is an M × N matrix, the n-th column of which is the power density distribution vector of the n-th source, with n = 1, …, N. The port temperature rises of the N sources form the N-row column vector ΔT(t) given by [42]

$$\mathsf{\Delta }T=Q^T\mathbf{\vartheta }(t)$$

(28)

As typical for MOR approaches, an $M\times \widehat{M}$ matrix V with $\widehat{M}\ll M$ is defined, which allows approximating $\mathbf{\vartheta }(t)$ by means of a reduced number $\widehat{M}$ of DoFs forming the $\widehat{M}$-vector $\widehat{\mathbf{\vartheta }}(t)$, so that

$$\mathbf{\vartheta }(t)=V\widehat{\mathbf{\vartheta }}(t)$$

(29)

The V matrix is used to project the heat conduction discretized problem (26)–(28) with the Galerkin’s method, thus deriving a DCTM in the form

$$\widehat{M}\frac{d\widehat{\mathbf{\vartheta }}}{dt}(t)+\widehat{K}\widehat{\mathbf{\vartheta }}(t)=\widehat{\mathbf{q}}(t)$$

(30)

being

$$\widehat{\mathbf{q}}(t)=\widehat{\mathbf{G}}P(t)$$

(31)

$$\mathsf{\Delta }T(t)={\widehat{G}}^T\widehat{\mathbf{\vartheta }}(t)$$

(32)

in which

$$\widehat{M}=V^{T}MV$$

(33)

$$\widehat{K}=V^{T}KV$$

(34)

are $\widehat{M}$-order matrices, and

$$\widehat{G}=V^{T}Q$$

(35)

is an $\widehat{M}\times N$ matrix. The V matrix is determined by the algorithm reported below.

At line 1 of the algorithm, the values σ_p, with p = 1, …, P_n, are automatically determined as a function of the desired relative error parameter ε as follows. First, the real positive quantities λ_n < Λ_n are properly estimated for the heat conduction problem with respect to the power impulse thermal response of the current p-th heat source. Next, the value P_n is determined as the smallest integer such that

$$4\exp \left(-P_n{\pi }^2/\log \left(4/k^{\prime }\right)\right)\le \epsilon$$

(36)

being k’ = λ_n/Λ_n. It is then set

$${\sigma }_p={\mathsf{\Lambda }}_n\mathrm{dn}\left(\frac{2p-1}{2P_n}K,k\right)$$

with p = 1, …, P_n, in which K is the complete elliptic integral of the first kind of modulus k, and dn is the homonymous elliptic function of modulus $k=\sqrt{1-{k^{\prime }}^2}$.

At line 3, the temperature response to the n-th heat source is solved in the complex frequency domain at the frequency value σ_p. Thus, equation

$$({\sigma }_{p}M+K){\Theta }_p=Q{e}_n$$

(37)

is solved for ${\Theta }_p$, e_n being an N-row vector having all zero elements but one at the n-th row. Since the complex frequency values are real positive, the coefficient matrices of the resulting linear systems are symmetric positive definite, and the most efficient multigrid iterative solvers can be used for their solution.

At line 4, the V_p matrix is achieved by appending to the columns of V_p−1 a vector derived orthogonalizing ${\Theta }_p$ with respect to the columns of V_p−1.

At line 5, the DCTM $𝒞$_p is determined proceeding as in (30)–(32), matrix V in (33)–(35) being substituted by matrix V_p.

At line 2, the temperature response ${\Theta }_p$ in the complex frequency domain due to the n-th independent heat source is estimated at σ_p using the DCTM $𝒞$_p−1. To this aim, equation

$$({\sigma }_p{\widehat{M}}_{p-1}+{\widehat{K}}_{p-1}){\widehat{\Theta }}_p={\widehat{Q}}_{p-1}{e}_n$$

(38)

is solved for ${\widehat{\Theta }}_p$. This vector is used to determine the approximation ${\tilde{\Theta }}_p=V_{p-1}{\widehat{\Theta }}_p$. At line 3 of the algorithm, such estimation of ${\Theta }_p$ is used to speed up the solution of (37) by setting the initial estimation of the solution.

It is noted that for the n-th heat source, with n = 1, …, N, the complex frequency values σ_p, with p = 1, …, P_n, are sorted in decreasing order. In this way, the linear systems introduced at line 3 of the algorithm are solved from the least to the most onerous ones in terms of CPU time. The computational burden is drastically reduced by this proper ordering, since the iterative solver can start from an estimate of the solution found in the previous step, which becomes increasingly accurate.

The DCTM achieved by this algorithm has dimension $\widehat{M}=P_1+P_2+ \dots +P_N$. Let Z(t) and $\widehat{Z}(t)$ be the N-th order power impulse thermal response [K/W] matrices of the discretized heat conduction problem (26)–(28) and of the DCTM, respectively. As shown in [43], it results in

$${\Vert Z(t)-\widehat{Z}(t)\Vert }_{{\mathscr{H}}_2}\le 2\epsilon {\Vert Z(t)\Vert }_{{\mathscr{H}}_2}$$

being

$$\parallel Z(t){\parallel }_{{\mathscr{H}}_2}=\sqrt{{\displaystyle {\int }_0^{+\infty }\mathrm{tr}}\left(Z^T(t)Z(t)\right)dt}$$

the Z(t) Hankel norm. Similar results can be extended from the time to the frequency domain and can be stated for the whole spatial-temporal temperature distributions due to power impulses [43]. All these results provide a strong theoretical guarantee for the convergence of the method. As a consequence of (36), this convergence is very fast, being exponential with respect to the numbers P of σ_p, with p = 1, …, P_n. As a result, in practice accurate DCTMs can always be obtained with small space-state dimensions, $\widehat{M}$.

It is now observed that by solving at limited cost the generalized eigenvalue problem

$${\widehat{U}}^T\widehat{M}\widehat{U}=I_{ \widehat{M}}$$

$${\widehat{U}}^T\widehat{K}\widehat{U}=\widehat{\Lambda }$$

having as unknown the $\widehat{M}$-order matrix $\widehat{U}$, in which $\widehat{\Lambda }$ is an $\widehat{M}$-order diagonal matrix, and introducing the change of variables

$$\widehat{\mathbf{\vartheta }}(t)=\widehat{U}\widehat{\xi }(t)$$

the DCTM equations (30)–(32) are transformed into

$$\frac{d\widehat{\xi }}{dt}(t)+\widehat{\Lambda }\widehat{\xi }(t)=\widehat{\Gamma }P(t)$$

(39)

$$\mathsf{\Delta }T(t)={\widehat{\Gamma }}^T\widehat{\xi }(t)$$

(40)

where $\widehat{\Gamma }$ is the $\widehat{M}\times N$ matrix ${\widehat{V}}^T\widehat{G}$. The temperature rise distribution is then reconstructed as

$$\mathbf{\vartheta }(t)=\Xi \widehat{\xi }(t)$$

(41)

$\Xi =V\widehat{U}$ being an $M\times \widehat{M}$ matrix like V.

Equations (39) and (40) can be interpreted as the equations ruling the equivalent network sketched in Figure 15.

Such a network, which can in principle describe the thermal behavior of any electronic component, is particularly suited to be solved by means of nodal analysis in SPICE-like simulators, since all elements are voltage-controlled current sources, and thus the number of variables is limited to $\widehat{M}.$ The topology is general and can be implemented into any circuit simulator.

4.4. Pre-Processing Evaluation of N and ε

During the pre-processing stage, two key parameters have to be chosen, namely, the number of elementary cells N (=79, as mentioned in Section 4.1) and the relative error parameter ε, the latter needed for the generation of the equivalent thermal network with FANTASTIC and set to 10⁻³. For the particular component under test, this discretization and this error were found to be a good trade-off between accuracy and CPU time needed for the ET simulation, as evaluated in an additional pre-processing analysis where some N and ε values were tested; considering much higher N and/or much lower ε leads to a marginal accuracy improvement paid with a significant increase in CPU time.

4.5. Construction of the Macrocircuit

The macrocircuit representing the ET behavior of the packaged DUT was constructed in the PSPICE environment as follows. The linear equivalent thermal network, provided in the form of a netlist, was enriched with N nonlinear voltage-controlled voltage sources to account for the Kirchhoff’s transformation (18), and the thermal feedback block was thus obtained. It must be remarked that nonlinear thermal effects could in principle be accurately described by extracting a fully nonlinear equivalent thermal network with a variant of the FANTASTIC tool [11]. However, the complexity of the nonlinear network grows with the discretization N much more rapidly than that of the linear counterpart (used in this paper); for N = 79, such a network would be composed by about 32 × 10⁶ elements, with insurmountable memory-storage problems. As a consequence, the adoption of a calibrated Kirchhoff’s transformation represents the only viable strategy to get accurate enough results. The ΔT and P_D nodes of the subcircuits (the individual transistor cells) were connected to the thermal feedback block. As shown in Figure 10c, all the gate terminals of the subcircuits were shorted together, as well as the drain and source ones, and it is possible to activate an electrical network to include the de-biasing over the source pad. A simplified scheme of the adopted strategy is shown in Figure 16.

After the circuit simulation run, the whole spatial-temporal temperature rise distribution can be reconstructed in a post-processing stage at negligible computational cost and memory storage using (41).

5. Static and Dynamic Electrothermal Simulations

The whole pre-processing activity, involving isothermal measurements, optimization of the model parameters, geometry/mesh construction in COMSOL, proper domain discretization, equivalent network extraction with FANTASTIC, calibration of the Kirchhoff’s transformation, and macrocircuit generation, can last a few days, after which any analysis can be performed in short times on the device of interest.

In particular, the macrocircuit was adopted to perform many PSPICE simulations of the DUT with DBC substrate in both static (dc) and dynamic conditions on a PC with an Intel Core i7-7700 (3.60 GHz) CPU and equipped with a 16 GB RAM. Unfortunately, experimental data for accuracy validation were not available for the examined component. On the other hand, the approach benefits from (i) a careful calibration of the transistor model from experimental data; (ii) a very accurate domain representation in COMSOL aided by the in-house routine mentioned in Section 4.2; (iii) an automated (and unaffected by user’s errors) extraction of the linear equivalent network through FANTASTIC; (iv) a smart pre-processing calibration of the Kirchhoff’s transformation. Hence, the only source of error can be induced by the degree of uncertainty concerning the thermal conductivities of the involved materials, which however affects any ET simulation approach.

First, the I_D–V_DS output characteristics were determined with a V_DS step amounting to 0.1 V under isothermal (at T₀) and ET conditions. Isothermal conditions were obtained by deactivating the thermal feedback block. The CPU time needed to simulate a single ET characteristic was nearly 100 s. Results are reported in Figure 17, which also shows the temperature rise above T₀ averaged over the effective active area (ΔT_av). It can be inferred that the simulation runs were stopped as ΔT_av reached 500 K.

Afterward, SC tests were simulated. Such tests involve large power dissipation, and are typically used to quantify the device robustness under harsh and abnormal events (see e.g., [7,23,25,44,45,46,47], all focused on SiC MOSFETs). In the SC experiment, the DUT is first biased in the OFF state with a given supply voltage applied to the drain, and then turned on with a single gate pulse (a gate resistance R_GATE of 50 Ω was considered). The knowledge of the whole temperature distribution over the active area is important, since the value and position of the temperature peak are needed for reliability considerations. The effects of many combinations of gate and supply voltages were examined. Simulation runs were very fast: a single test required about 300 s with a fine time discretization. The total drain current I_D conducted by the DUT vs. time is shown in Figure 18a for all the analyzed cases, while Figure 18b illustrates the corresponding temperature rises ΔT_av. The first figure reveals that I_D first grows due to the strong positive temperature coefficient induced by (i) the reduction in threshold voltage and (ii) the mobility increase (for the lowered Coulomb scattering), and then drops since the negative temperature coefficient triggered by the acoustic-phonon scattering dominates (thermally stable behavior) as all the traps have released electrons [7,25].

The temperature rises ΔT over cells #01, #23, #47, #61 (identified in the layout of Figure 10) are reported in Figure 18c for two cases, namely, V_GS = 10 V/V_DD = 200 V and V_GS = 20 V/V_DD = 200 V. From the inspection of the waveforms, it is found that a pronounced temperature nonuniformity takes place in the first case (the inner cell #47 suffers from a ΔT two times higher than that of the top-corner cell #01), which can be explained as follows. The milder bias conditions (V_GS = 10 V) allowed the device to safely undergo the test for a longer period, within which the heat had enough time to significantly spread, thereby favoring a stronger impact of the mutual thermal interactions and the consequent exacerbation of temperature gradients.

The whole temperature field in the domain was determined at chosen time instants from the DCTM generated with FANTASTIC in a post-processing step. More specifically, points A (t = 25.5 µs, V_GS = 20 V/V_DD = 200 V), B (t = 338 µs, V_GS = 20 V/V_DD = 50 V), and C (t = 14 ms, V_GS = 10 V/V_DD = 50 V) identified in Figure 18a,b were selected, which approximately share the same ΔT_av value, i.e., 500 K. The computed temperature rise maps are shown in Figure 18d (top view) and Figure 18e (side view). As can be seen, despite the same ΔT_av,

• point A, which falls at a time instant close to the beginning of the test, is endowed with a uniform and superficial temperature field, since the heat is still confined in the top active region close to the generation area;
• the temperature distributions in B and C are increasingly uneven, which is again ascribable to the much longer stress times. In particular, C suffers from a severe temperature focusing over the innermost DUT area, and—as witnessed by the side view—the downward heat had enough time to reach and hit the DBC.

Lastly, the macrocircuit was used to simulate two UIS tests, which are commonly adopted to evaluate the maximum amount of avalanche energy sustainable by the device [8,44]. In the UIS experiment, a load inductor L tied to the drain is first ramped up to a desired current by keeping the DUT in linear mode for a time t_ON; then the DUT is turned off and brought into avalanche by the inductor, which preserves the current continuity. The tests will be hereinafter denoted as case #1 and #2, the specifics of which are as follows:

• case #1: V_DD = 300 V, L = 4.6 mH, R_GATE = 15 Ω
• case #2: V_DD = 600 V, L = 12 mH, R_GATE = 15 Ω

In both cases, a gate voltage equal to 20 V was applied for t_ON = 200 µs, and then lowered to 0. Again, the time elapsed by PSPICE for a single test was about 300–400 s. Concerning case #1, Figure 19a shows the drain current I_D and the drain-source voltage V_DS vs. time in the span 190 to 280 µs, while Figure 19b illustrates the temperature rises of cells #01, #23, #47, #61. An almost uniform temperature distribution was found, as witnessed by the map taken at time instant t* = 220 µs, where the dynamic temperature averaged over the active area peaks (Figure 19c).

In case #2, despite the lower drain current I_D ≈ V_DD·t_ON/L before turn-off (10 A instead of 13 A obtained for case #1), the average temperature reaches higher values due the longer discharging transient, in turn dictated by the higher L, since

$$\frac{d{I}_D}{dt}=-\frac{B{V}_{DS}\left(T\right)-V_{DD}}{L}$$

(42)

Results are shown in Figure 20. As can be inferred from Figure 20b, which reports the temperature rises of the selected cells, and from Figure 20c, which depicts the post-processing temperature maps at t* = 240 µs, a slightly more nonuniform temperature field occurs with respect to case #1.

6. Conclusions

In this paper, an advanced circuit-based approach for the static and dynamic electrothermal simulation of multicellular SiC power MOSFETs has been proposed, which—differently from the strategies encountered in the literature—seems to represent a good trade-off between accuracy and efficiency. The device is discretized into a chosen number of elementary cells (heat sources) and turned into a purely-electrical macrocircuit, where (i) the cells are described with subcircuits accounting for the key and harmful influence of SiC/SiO₂ interface traps, and (ii) the power-temperature feedback is modeled with an equivalent thermal network. This network is obtained through a fully automated process: first, a 3D mesh representing the device is generated in COMSOL with the support of an in-house routine that elaborates the layout files and makes use of further information about thickness of layers, position/shape of the heat sources, material parameters, and boundary conditions; then, such a mesh is provided as an input to the FANTASTIC tool, which derives a dynamic compact thermal model of the device and the related equivalent network. The macrocircuit can be solved in the environment of any SPICE-like simulator with short CPU time and unlikely occurrence of convergence issues. The effectiveness and efficiency of the proposed strategy have been verified on a 1200 V, 50 A multicellular 4H-SiC VDMOS soldered on a DBC package. It has been shown that the simulation of short-circuit and unclamped inductive switching tests requires only 300–400 s on a normal PC, despite the critical conditions and the fine time discretization, and that potentially-dangerous temperature gradients/hogging can be easily identified. It can be concluded that the proposed approach can be helpful for industry engineers who are in charge of optimizing the thermal design of multicellular power devices in any technology.