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

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/SiO2 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.


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 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 TB 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].

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.

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.

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.

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. 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 TB 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].

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.

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.
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 then the DUT is driven into pinch-off, and The negative temperature coefficient of V TH is described through the following law: 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 where the exponent m(T) is II effects are accounted for as follows [7,9]. The bias-and temperature-sensitive avalanche multiplication factor M (≥1) is given by [12] M(V DS , I D , where BV DS (T) is the temperature-dependent drain-source breakdown voltage, expressed as 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 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 + ξ · (I leak + I DnoII ) (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]: Energies 2020, 13, 4563 (11) R JFET (T 0 ) being the JFET resistance at T = T 0 , V drift » V 1 , and V GS = V 2 (V 1 and V 2 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]: and while C GS was not modified [13].

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: with 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 0 ) (≈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 2 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 0 ), 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].
Energies 2020, 13, 4563 6 of 27 contribute to channel formation, thereby leading to a high threshold voltage VTH(T0), and (ii) VTH 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 ID, 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 VTH and K models is witnessed by the comparison reported in Figure 5 between (2) and the ID-VGS transfer characteristics measured under isothermal conditions at VDS = 20 V and various TB values in a current range wherein the DUT operates in pinch-off (note that ID ≈ IDnoII). An inspection of the curves reveals the considerable positive temperature coefficient of ID within a wide range of currents. 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).
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 ID, 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 VTH and K models is witnessed by the comparison reported in Figure 5 between (2) and the ID-VGS transfer characteristics measured under isothermal conditions at VDS = 20 V and various TB values in a current range wherein the DUT operates in pinch-off (note that ID ≈ IDnoII). An inspection of the curves reveals the considerable positive temperature coefficient of ID within a wide range of currents. 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 Rdrift in (10), (11) were optimized by comparing the ID-VDS output characteristics at various VGS and different TB values in triode region with the model (1) applied to the scheme in Figure 2. Figure 6 reveals the accuracy of the extraction at TB = 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 TBdependent ID-VDS avalanche curve at VGS = 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 CGD (12) and CDS (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.  experimental model 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 TBdependent ID-VDS avalanche curve at VGS = 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 CGD (12) and CDS (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.  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.

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 0 ) to 4 V, (ii) decreasing a VTH to 2 mK −1 , (iii) multiplying the current factor K(T 0 ) by 50 to annihilate the degradation of mobility µ n (T 0 ), 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] for the real DUT and its traps-free variant at reference temperature T 0 and drain-source voltage V DS = 20 V; it is again witnessed that the traps at the SiC/SiO 2 interface lead to a detrimental highly-positive α T within a broad range of currents.
Energies 2020, 13, 4563 9 of 27 drain current ID (<20 A), where the negative temperature coefficient of VTH 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 ID reduces with temperature.  the traps-free transistor benefits from a much higher current capability due to the lower VTH and the higher μn. To further corroborate the above findings, Figure 8 reports the temperature coefficient of the drain current ID, given by [19,20,25] for the real DUT and its traps-free variant at reference temperature T0 and drain-source voltage VDS = 20 V; it is again witnessed that the traps at the SiC/SiO2 interface lead to a detrimental highly-positive αT within a broad range of currents.

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 T0 = 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 VTH and of the current factor K, the bias-and temperature-dependent drift resistance Rdrift, as well as the bias-and temperature-dependent II mechanism. The TEOL is adopted, namely, the temperature rise over ambient ΔT = T − T0 is actually a voltage, while the dissipated power PD is treated as a current; this allows (i) enabling the temperature sensitivity of the key physical

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 0 = 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 0 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] 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.

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 0 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.
Consequently, the II-unaffected current IDnoII is obtained as where m(T) is given by (5). II effects are activated by adding an avalanche current IDII = ξ‧(Ileak + IDnoII) to IDnoII through the current source C, ξ = M − 1 being bias-and temperature-dependent according to (6). Hence, the drain current ID flowing through the drift resistance Rdrift is given by (9).
The bias-and temperature-dependent Rdrift described by (10) with (11) is taken into account by making use of source D, which imposes the voltage drop Figure 9. Sketch of the SPICE-compatible subcircuit for the transistor cell.
Only half of the DUT was represented and simulated by exploiting its inherent symmetry. The effective active region (≈5 mm 2 ) was partitioned into N = 79 cells, each with a 250 × 250 µm 2 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.
and M n conducts the current I D (M n ) given by 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 Consequently, the II-unaffected current I DnoII is obtained as 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 Only half of the DUT was represented and simulated by exploiting its inherent symmetry. The effective active region (≈5 mm 2 ) was partitioned into N = 79 cells, each with a 250 × 250 µm 2 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.

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]. Figures 11 and 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 (Al2O3) for the proposed case study; however, alternative materials like silicon nitride (Si3N4) [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.

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]. Figures 11 and 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 2 O 3 ) for the proposed case study; however, alternative materials like silicon nitride (Si 3 N 4 ) [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 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 (Al2O3) for the proposed case study; however, alternative materials like silicon nitride (Si3N4) [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 5 and 5.2 × 10 5 , 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.  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 5 and 5.2 × 10 5 , 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.  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 5 and 5.2 × 10 5 , 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 TB = T0 was applied at the 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 0 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. Nonlinear thermal effects were accounted for by including the temperature dependences of the thermal conductivities described by 18.5 [39] 787 [39] 3100 [39] −0.33 [39] 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.

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 ϑ(t), are solutions of the discretized linear dynamic heat conduction problem in which the power density distribution vector q(t) takes the form 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] ∆T = Q T ϑ(t) As typical for MOR approaches, an M ×M matrix V withM M is defined, which allows approximating ϑ(t) by means of a reduced numberM of DoFs forming theM-vectorθ(t), so that 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 ∆T(t) =Ĝ Tθ (t) (32) in whichM is solved forΘ p . This vector is used to determine the approximation Θ p = V p−1Θp . At line 3 of the algorithm, such estimation of Θ 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 dimensionM = P 1 + P 2 + . . . + P N . Let Z(t) andẐ(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 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,M.
It is now observed that by solving at limited cost the generalized eigenvalue problem having as unknown theM-order matrixÛ, in whichΛ is anM-order diagonal matrix, and introducing the change of variablesθ (t) =Ûξ(t) the DCTM equations (30)-(32) are transformed into whereΓ is theM × N matrixV TĜ . The temperature rise distribution is then reconstructed as 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 toM. The topology is general and can be implemented into any circuit simulator. 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  .
M The topology is general and can be implemented into any circuit simulator.

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 -3 . For the particular component under test, this discretization and this error were found to be a good tradeoff between accuracy and CPU time needed for the ET simulation, as evaluated in an additional preprocessing 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.

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 6 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 PD 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.

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 −3 . 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.

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 6 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). 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).

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 0 ) 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 0 averaged over the effective active area (∆T av ). It can be inferred that the simulation runs were stopped as ∆T av reached 500 K. re 16. Schematic representation of the proposed strategy to perform a fully coupled ET analysis circuit simulation tool: the feedback loop between the electrical circuit (left) and the thermal back block (TFB, right) relying on the DCTM-based equivalent network is highlighted. er the circuit simulation run, the whole spatial-temporal temperature rise distribution can be ucted in a post-processing stage at negligible computational cost and memory storage using and Dynamic Electrothermal Simulations whole pre-processing activity, involving isothermal measurements, optimization of the parameters, geometry/mesh construction in COMSOL, proper domain discretization, nt network extraction with FANTASTIC, calibration of the Kirchhoff's transformation, and rcuit generation, can last a few days, after which any analysis can be performed in short times evice of interest. articular, the macrocircuit was adopted to perform many PSPICE simulations of the DUT C substrate in both static (dc) and dynamic conditions on a PC with an Intel Core i7-7700 z) CPU and equipped with a 16 GB RAM. Unfortunately, experimental data for accuracy n were not available for the examined component. On the other hand, the approach benefits a careful calibration of the transistor model from experimental data; (ii) a very accurate representation in COMSOL aided by the in-house routine mentioned in Section 4.2; (iii) an ed (and unaffected by user's errors) extraction of the linear equivalent network through STIC; (iv) a smart pre-processing calibration of the Kirchhoff's transformation. Hence, the urce of error can be induced by the degree of uncertainty concerning the thermal ivities of the involved materials, which however affects any ET simulation approach. t, the ID-VDS output characteristics were determined with a VDS step amounting to 0.1 V under al (at T0) and ET conditions. Isothermal conditions were obtained by deactivating the feedback block. The CPU time needed to simulate a single ET characteristic was nearly 100 s are reported in Figure 17, which also shows the temperature rise above T0 averaged over tive active area (ΔTav). It can be inferred that the simulation runs were stopped as ΔTav 500 K. Afterward, SC tests were simulated. Such tests involve large power dissipation, and are ty used to quantify the device robustness under harsh and abnormal events (see e.g., [7,23,25,44-focused on SiC MOSFETs). In the SC experiment, the DUT is first biased in the OFF state with a supply voltage applied to the drain, and then turned on with a single gate pulse (a gate res RGATE of 50 Ω was considered). The knowledge of the whole temperature distribution over the area is important, since the value and position of the temperature peak are needed for rel considerations. The effects of many combinations of gate and supply voltages were exa Simulation runs were very fast: a single test required about 300 s with a fine time discretizatio total drain current ID conducted by the DUT vs. time is shown in Figure 18a for all the analyzed while Figure 18b illustrates the corresponding temperature rises ΔTav. The first figure reveals first grows due to the strong positive temperature coefficient induced by (i) the reduction in thr voltage and (ii) the mobility increase (for the lowered Coulomb scattering), and then drops sin 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. 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: 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 both cases, a gate voltage equal to 20 V was applied for tON = 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 ID and the drain-source voltage VDS 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 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.
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.

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 2 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 Energies 2020, 13, 4563 24 of 27 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.

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