The Generalized Falk Method (GFM) for coordinate transformation, together with two model-reduction strategies based on this method, are presented for efficient coupled field-circuit simulations. Each model-reduction strategy is based on a decision to retain specific linearly-independent vectors, called trial vectors, to construct a vector basis for coordinate transformation. The reduced-order models are guaranteed to be stable and passive since the GFM is a congruence transformation of originally symmetric positive definite systems. We also show that, unlike the Padé-via-Lanczos (PVL) method, the GFM does not generate unstable positive poles while reducing the order of circuit problems. Further, the proposed GFM is also faster when compared to methods of the type Lanczos (or Krylov) that are already widely used in circuit simulations for electrothermal and electromagnetic problems. The concept of response participation factors is introduced for the selection of the trial vectors in the proposed model-reduction methods. Further, we present methods to develop simple equivalent circuit networks for the field component of the overall field-circuit system. The implementation of these equivalent circuit networks in circuit simulators is discussed. With the proposed model-reduction strategies, significant improvement on the efficiency of the generalized Falk method is illustrated for coupled field-circuit problems.
Nonlinear circuitspower electronicsIGBTMOSFETmodel-order reductionmodelingmodule generationpower modeling and estimationelectrothermal simulationIntroduction
Except for very simple systems, computer simulations involving coupled field-circuit problems are hindered by the very large number of degrees of freedom involved in the discretized field problem. The evaluation of the exact natural frequencies and mode shapes for large complex systems often requires a large number of numerical operations. The direct engineering significance of this information, however, may be of limited value. For most analyses, one potentially useful approach to overcome this computational difficulty is model-order reduction, where parts of the original large-order model are replaced by models of substantially lower order, yet capable of capturing the behavior of the original subsystem with sufficient engineering accuracy.
Reduced-order models have been developed in many areas of engineering, such as circuits (e.g., [1–4]), power systems (e.g., [5–7]), electromagnetics (e.g., [8 –10]), fluid mechanics (e.g., [11–13]), nonlinear structural mechanics and earthquake engineering (e.g., [14]), nonlinear hydraulic fracturing problems (e.g., [15]), etc., to cite a few. Deep-learning artificial neural networks have been introduced to build on more traditional model-order reduction methods, such as the Proper Orthogonal Decomposition (POD)^{1}
Proper Generalized Decomposition was used for model-order reduction in nonlinear hydraulic-fracture problems in [15].
, to increase computational efficiency [16–21].
The modal superposition method and the component-mode synthesis method for model reduction in coupled field-circuit simulations require solving a large and expensive eigenvalue problem [22]. Several alternative methods for model reduction without solving an eigenvalue problem are the well-known method in [23] and the WYD (Wilson–Yuan–Dickens) method in [24], both belonging to the general Krylov subspace methods [25, 26]. The basic idea is to find a convenient and inexpensive transformation of coordinates to simplify the equations of motion that are easier to solve. The WYD reduction method—sometimes called the Ritz reduction method [27] or Ritz vector method [28]—is based on the direct superposition of Ritz vectors constructed from the spatial distribution of the specified dynamic excitation; the generated Ritz vectors are thus excitation dependent. Another reduction technique called the AWE (Asymptotic Waveform Evaluation) method, a moment-matching technique based on the Padé approximation in [29], was also applied to coupled field-circuit simulation for integrated circuits (ICs) [30]. In the AWE method, the system behavior is approximated with a lower-order model, obtained from a Taylor series expansion around the Laplace transform variable s = 0. The AWE method is efficient to extract the low-frequency behavior of the circuit.
The WYD method is an efficient alternative to the classical mode superposition as discussed in [22]. In the WYD method, the solution is a direct superposition of a special class of Ritz vectors (approximated eigenvectors) generated from the spatial distribution of the excitation source; it eliminates the requirement for exact evaluation of the natural frequencies and mode shapes. It is important to note here that while the WYD method uses the excitation-dependent trial vectors to form a basis to obtain a reduced eigenvalue problem, whose solution leads to the Ritz vectors (or approximated eigenvectors), our methodology based on the Generalized Falk Method (GFM), introduced in [31], bypasses the solution of the reduced eigenvalue problem and the use of the Ritz vectors, but instead uses directly the trial vectors. To minimize the computational effort while obtaining good accuracy, the WYD method uses only Ritz vectors having large participation factors, i.e., Ritz vectors that carry a significant amount of information about the excitation.
The recursive iterative scheme used in the WYD method was shown in [31] to bear a great mathematical similarity to the Lanczos method. If the same starting vector were used in both the WYD and the Lanczos method, then the generated trial vectors from both method are very similar to each others.
Model reduction techniques based on Krylov-subspace iterations, especially the Lanczos method and the Arnoldi method^{2}
The Lanczos method is applicable to Hermitian matrices, and is a special case of the Arnoldi method, which is applicable to non-Hermitian matrices.
, are popular tools to tackle the large-scale time-invariant linear dynamical systems that arise in the simulation of electronic circuits, e.g., [32–40], and electromagnetics, e.g., [9, 41]. In [42] and [43], the AWE and Padé-via-Lanczos (PVL) model reduction procedures are extended to the use of electromagnetic analysis. Model reduction based on the Padé approximation may, however, produce positive poles; and is thus numerically unstable [36, 40]. The resulting reduced-order models must be post-processed to eliminate the unstable modes; such post-processing method, which must preserve the moment-matching properties, does not guarantee the elimination of unstable modes [39].
The Generalized Falk Method (GFM) proposed in [31] is shown here (Section 2.2) to be able to avoid unstable positive poles (Section 2.2.3), as the GFM can predict actually higher-order poles. We also show that the GFM is more efficient compare to other Lanczos-type methods. In [37], the block Arnoldi method is implemented via a simple congruence transformation, which guarantees the passivity and thus stability of the generated reduced-order model. Unlike the GFM introduced here, eigen decomposition is still needed, however, in the Passive Reduced-Order Interconnect Macro modeling Algorithm (PRIMA). Similar to the block Arnoldi method, the number of moments matched in PRIMA implementation is essentially the same as that for the block Arnoldi method, which is 50% of those moments matched by the block Lanczos method. The GFM applied to coupled field-circuit problems such as electrothermal simulations, on the other hand, has exactly the same order of accuracy as the Lanczos method since the resulting finite element matrices are symmetric and nonnegative. In particular, the capacitance matrix is positive definite, and thus invertible. For general circuit problems where the capacitance and conductance matrices are not symmetric, the GFM has the same order of accuracy as the Arnoldi method.
In this paper, to study the truncation effect in the proposed GFM, we first introduce the concept of participation factors for the selection of the trial vectors (Section 2.3).
To motivate, real-world applications of the electronic devices analyzed in this paper are presented in a top-down fashion, going from the big-picture system level down to the detailed component level (Section 2.5). There are many applications of semiconductors in power electronic devices used in data centers, electric vehicles, renewable energy (solar, wind, wave and tidal power). An important real-world application system is the Tesla Model S^{3}
See, e.g., the video “How does an Electric Car work? Tesla Model S” Youtube, from time 3:56, where the motor inverter and the induction motor correspond to the motor drive and the electric motor/generator in Fig. 1, respectively. See also “In a Tesla Model S, there is no IGBT packaging trick”, Website (Internet archive), where the brand name on the 3-pin IGBT chip in TO-247 package is International Rectifier (IR), with the circle between the letters I and R being the diagram of a diode drawn vertically; see Wikipedia version 09:24, 10 November 2020. International Rectifier was acquired by Infineon (Fig. 5) in 2015. See also p. 17 of the IR 2013 annual presentation to investors Online pdf (Internet archive) and the video “Pioneers in Clean Technology-Marc Tarpenning-Tesla Motors” Youtube at time 1:02:51 when Tarpenning, a Tesla founder, mentioned the new, efficient IGBT by IR, which responded to the need to switch power at 400 V and 1000 Amp for the Tesla.
, for which the wheels are turned by an induction motor driven by a 3-phase AC generated by a motor inverter taking as input the DC from its battery pack; see Fig. 1. At the detailed component level, the design of the Tesla motor inverter relied on IGBTs^{4} [44
IGBT = Insulated-Gate Bipolar Transistor, “a three-terminal power semiconductor device primarily used as an electronic switch which, as it was developed, came to combine high efficiency and fast switching”; see Wikipedia version 11:35, 28 November 2020. See Figs. 5–10 in Section 2.5.
–48].
<italic>Battery electric vehicle schematic</italic> (<xref ref-type="sec" rid="sec1">Section 1</xref>). The motor drive and electric motor/generator correspond to the motor inverter and the induction motor in <xref ref-type="fn" rid="fn1_3">Footnote 3</xref>. The on-board charger is an AC-to-DC converter. The power converter transform DC to DC at different level to charge the battery and to provide DC to the motor drive (inverter). <uri xlink:href="https://ocw.tudelft.nl/course-readings/2-1-2-lecture-notes-parts-of-an-ev-recap/">OpenCourseWare, TU Delft</uri> <uri xlink:href="https://creativecommons.org/licenses/by-nc-sa/4.0/">(CC BY-NC-SA 4.0)</uri>
Remark 1.1 In writing this paper, we are mindful of the general readers, not just experts in power electronics. As a result, sufficient motivational applications and explanation for the electronic devices used here are provided. Whenever possible, permanent links to references and supporting documents on the web are provided. Wikipedia articles are provided with specific versions that are permanently fixed with no further edits, such as Wikipedia version 09:24, 10 November 2020. While a link such as this Online pdf may no longer exist in the future, its Internet archive version will never disappear. See also, e.g., Footnote 3.
Then we introduce two model-reduction strategies to further improve the efficiency of coupled field-circuit problems (Section 2.5), such as those encountered in electrothermal simulation or in the analysis of IC interconnects.
To illustrate the proposed model-reduction methodology, we consider a case study of coupled field-circuit problems involving electrothermal simulations; the same methodology can be applied to other coupled field-circuit problems such as IC interconnect analysis (Section 3). We next describe methods to develop simple equivalent circuit networks for the field component of the overall coupled field-circuit system (Section 3.2). We also discuss how to implement these field circuit networks into general circuit simulators (Sections 3.3, 3.4).
Numerical examples involving a full bridge converter^{5}
A full-bridge converter transform a direct current (DC) at a certain voltage to a DC at a different voltage for use in many electronic devices. See, e.g., the video “Full Bridge Converter-Basics of Switching Power Supplies”, [49], for an explanation of hard switching, soft switching, and phase-shift full-bridge converter.
using IGBTs (Section 4.2) and a Voltage Regulator Module (VRM)^{6}
“A voltage regulator module (VRM), sometimes called processor power module (PPM), is a buck converter that provides a microprocessor the appropriate supply voltage, converting +5 V or +12 V to a much lower voltage required by the CPU, allowing processors with different supply voltage to be mounted on the same motherboard. On personal computer (PC) systems, the VRM is typically made up of power MOSFET devices, see Wikipedia version 22:07, 31 December 2020. For application of VRM in data centers, see [50].
using MOSFETs^{7}
MOSFET = Metal-Oxide-Semiconductor Field-Effect Transistor, is a kind of transistor in which electric field controls current flow, and the device is fabricated by controlled oxidation of silicon semiconductor; see Wikipedia 23:12, 27 November 2020.
(Section 4.3) establish the computational ease of use, accuracy, and efficiency of the present methodology; see also [51, 52]. A comparision of CPU times obtained from MATLAB and from the circuit simulator SABER is provided (Section 4.4).
Generalized Falk Method (GFM) and Model Reduction
In this paper, we present two methods to produce reduced-order models based on a combination of the proposed Generalized Falk Method (GFM) and the traditional WYD methods. These new model-reduction techniques are more stable than the traditional Lanczos method, and maintain a satisfactory accuracy. The selection of the trial vectors that form the basis for the transformation of coordinates of the original system to a reduced-order system can be based on the participation factors of these trial vectors to the dynamic response of the system. Once a reduced-order model is obtained, the dynamic response can be computed by direct numerical integration or by modal superposition method [22].
Case Study for Coupled Field-Circuit Problems
As mentioned above, to illustrate the proposed methodology, we consider electrothermal simulation as a case study of coupled field-circuit problems. The same methodology can be applied to other coupled field-circuit systems such as IC interconnects.
A typical nonlinear electrothermal coupled system can be generally expressed as
u˙=f(u,e,T),P=P(u),
where u is the matrix containing the nodal voltages or currents, e the electrical input, T the temperature, and P the electrical power loss. The nonlinear electrical system is governed by the semiconductor equations (PDE’s) and circuit equations (ODE’s). The electrical power loss is originated from the heat generated inside the semiconductor devices. In practice, the electrical device needs to be treated as a distributed heat source. Here we neglect the distributed heat source, and assume the heat source is a concentrated input load to the thermal domain.
The heat diffusion problem on a domain Ω with boundary ∂Ω=Γ1∪Γ2¯ is governed by the partial differential equation
where T:Ω×R+→R is the temperature (a function of space and time), κ:Ω→R the thermal conductance, ρ:Ω→R the mass density, cp:Ω→R the specific heat, P:Γ1×R+→R the input power from boundary Γ_{1}, A the input-power cross-section area, h the convection coefficient, n the outward normal vector to the boundary, and T_{a} the ambient temperature. In the above, the set of real numbers is denoted by R, with R+:={t∈R∣t≥0}, denoting the set of non-negative numbers. The domain of the heat sink is denoted by Ω⊂Rs, with s = 1, 2, 3 being the space dimension. The boundary ∂ Ω of the heat sink is decomposed into two parts: Γ_{1} for the input-power boundary condition Eq. (3), and Γ_{2} for the convective boundary condition Eq. (4). We assume in the above electrothermal coupled system that it does not include any heat source in the thermal domain Ω.
Employing a Galerkin finite element projection [53, 54] on the PDE Eq. (2), we obtain the following semi-discrete equation for the thermal system:
Cd˙+Kd=f=Ip+Ic,d(0)=d0,
where C∈Rn×n, with n being the number of element nodes (or temperature unknowns/degrees of freedom or dofs), is the capacitance matrix, K∈Rn×n the conductance matrix, f∈Rn×1 the heat supply matrix, I_{p} the input power supply from the electrical component, I_{c} the heat supply coming from the convective boundary condition, and d(t)={di(t)}∈Rn×1 the column matrix of the time-dependent nodal temperatures di(t) at the finite-element (FE) node i, for i = 1,… , n. Let N_{i} be the shape function associated with node i, then the temperature field over the domain of the heat problem is approximated by T(x,t)≅∑i=1nNi(x)di(t). The capacitance matrix C can be realized by capacitor elements, the conductance matrix K by resistor elements, and the heat supply vector f by current sources, so that thermal effects for semiconductor devices (transistors) like an IGBT or MOSFET can be dynamically incorporated into circuit simulators with thermal connection. The global matrices C, K, f are assembled from the elemental matrices c_{e}, k_{e} and f_{e}.
A complex mechanical system can be transformed into a simple one-dimensional mass-spring system by transforming its coordinates (degrees of freedom) using a set of linearly-independent vectors, termed trial vectors, generated by the Falk method [55]. The transformed matrices are diagonal (mass) and tridiagonal (stiffness), which upon scaling represent a one-dimensional lumped system, and have the same eigenvalues as the original system.
Generalized Falk Method (GFM)
The original Falk method suffers from the loss of orthogonality similar to that found in the original Lanczos method [56, 57], since each trial vector in the Falk method is orthogonalized with respect to only the two previously generated trial vectors. By contrast, in the GFM, each trial vector is orthogonalized with respect to all previously generated trial vectors, with the transformed tridiagonal conductance matrix constructed by formal matrix multiplication [31].
Conceptual Algorithm
To ensure that the generated trial vectors^{8}
Even though the terms “vector” and “column matrix” are used here interchangeably here, the latter is preferred when it adds to clarity.
contain important information from the static response, the starting trial vector, denoted by w1∗∈Rn×1 in Line 1 of Algorithm 2.1, with n being the total number of dofs, is selected to be the static response when the spatial distribution of the dynamic excitation is applied on the system.
The traditional WYD method leads to an identity capacitance matrix C∗=WnTCWn=In∈Rn×n, with I_{n} being the identity matrix and Wn=[w1,…,wn]∈Rn×n a matrix whose columns are the trial vectors w_{i}, but a full conductance matrix K∗=WnTKWn∈Rn×n, thus making the construction of an equivalent circuit complicated [31]. The WYD process can be truncated based on certain measure (to be discussed shortly) to obtain reduced-order models.
On the other hand, the GFM leads to an identity capacitance matrix C∗=In∈Rn×n, and a tridiagonal conductance matrix K∗∈Rn×n. The GFM can also be truncated, and the generated trial vectors can be used to construct reduced-order models. There is a trade-off in the solution of the method to use to reduce the order of a large system in a cost-effective manner, while maintaining high accuracy.
Implementable Algorithm
The implementation of the present GFM is presented as a pseudocode in Algorithm 2.2, with particular attention paid to the computational efficiency.
GFM <italic>vs</italic>. Other Methods: No Unstable Positive Poles
It is well known that direct Padé approximation that forms the basis for methods of type Lanczos does not guarantee that the reduced-order model is passive or stable [36, 37, 58]: It is possible to have unstable positive poles. By contrast, no positive pole is generated when using the GFM in circuit simulation since the higher order poles tend to be accurately predicted.
Unlike the Padé-via-Lanczos (PVL) method, Table 1 clearly shows that the GFM does not generate destabilizing positive poles in the system, which was originally stable. The PVL method was shown in [36], using the simple linear RC-circuit in Fig. 2, to possibly produce a model that was unstable even though the system was described by a symmetric, positive definite matrix. The circuit parameters are R1 = R2 = R3 = R4 = R5 = R = 1.3955 and C1 = C2 = C3 = C4 = C5 = 1. The original system capacitance matrix and conductance matrix are
and they are reduced from 5× 5 to 3× 3 in the PVL, Arnoldi, and GFM iterations. The GFM does not suffer from positive poles.
<italic>Poles obtained from different model-reduction methods using Krylov vectors</italic> (<xref ref-type="sec" rid="sec2_2">Section 2.2</xref>): Padé-via-Lanczos (PVL), Arnoldi, GFM. Linear RC circuit in <xref ref-type="fig" rid="fig2">Fig. 2</xref>, with all resistances equal to 1.3955, and all capacitances equal to 1
Krylov methods
Negative poles
Positive poles
1
2
3
4
5
1
Exact
−2.4715
−1.8999
−1.1513
−0.4633
−0.0544
–
PVL
–
–
–
−0.4954
−0.0581
4.6671
Arnoldi
–
–
−1.2865
−0.4972
−0.0581
–
GFM
–
–
−1.6881
−0.6572
−0.0581
–
<italic>Simple linear RC circuit</italic> (<xref ref-type="sec" rid="sec2_2">Section 2.2</xref>). A simple linear circuit with passive elements shows that the Padé-via-Lanczos (PVL) method generates positive poles that are unstable
To detect the positive poles generated from the PVL method and the circuit resistance values R that lead to the positive poles for the circuit in Fig. 2, we swept through a set of values for the resistance R, and plot the poles generated from these three methods. Fig. 3 (left) shows clearly that the PVL method generated unstable poles, whereas the Arnoldi method (Fig. 3, right) and the GFM method (Fig. 4) did not.
<italic>Simple linear circuit</italic> (<xref ref-type="sec" rid="sec2_2">Section 2.2</xref>). Values of the 3rd pole as a function of the resistance <italic>R</italic>, as generated by the Padé-via-Lanczos (PVL) method (left), and by the Arnoldi method (right) for the circuit example in <xref ref-type="fig" rid="fig2">Fig. 2</xref>
<italic>Simple linear RC circuit</italic> (<xref ref-type="sec" rid="sec2_2">Section 2.2</xref>). Values of the 3rd pole as a function of the resistance <italic>R</italic>, as generated by the GFM for the circuit example in <xref ref-type="fig" rid="fig2">Fig. 2</xref>
Participation Factors
To study the truncation effect in our proposed coordinate-transformation methods, we introduce three quantities that will serve as measures to guide the truncation of the basis of trial vectors that are used to construct reduced-order models: The projected loadp, the load participation factorβ, and the response participation factorα.
The superscripts L, W, F, G over these three quantities will be used to designate their association with the Lanczos, WYD, Falk, and GFM methods, respectively. Let yr(t)={y1(t),…,yr(t)}T∈Rr×1, with r ≤ n being a mnemonic for “retained” after truncation of the full basis and n the dimension of the original system (i.e., the size of the full basis, see Line 3 in Algorithm 2.1), denote the state variables in the new coordinates corresponding to the basis of the trial vectors generically written as {w1,…,wr}∈Rn×r, which can be either Lanczos vectors (L), WYD vectors (W), Falk vectors (F) or Generalized Falk vectors (G). The projected force on the basis of trial vectors is represented by pr(t)={p1,…,pr}T∈Rr×1. The load participation factor β_{i} and the response participation factor α_{i} are defined as
βi(t):=∥wipi(t)∥,αi(t):=∥wiyi(t)∥,
respectively.
Lanczos Method
A projection of the matrix ODE in Eq. (6) onto the Lanczos basis vectors yields the following transformed matrix ODE [27] (see Eq. (43) in Appendix A):
Try˙r+Iryr=ZrTCK−1f=:prL,
where Tr∈Rr×r, with r ≤ n and n being the size of the full basis as mentioned above, is a tridiagonal matrix consists of orthogonalization coefficients given in Eq. (39) and Eq. (41), which are similar to the coefficients c_{j} in Algorithm 2.1, Ir∈Rr×r an identity matrix, Zr:=[w1L,…,wrL]=[z1,…,zr]∈Rn×r the matrix containing the r Lanczos trial vectors (see Appendix A), C∈Rn×n the capacitance matrix, K−1∈Rn×n the inverse of the conductance matrix, and the (truncated) projected load prL={p1L,…,prL}T∈Rr×1 with r retained components of pnL∈Rn×1 is as defined on the right-hand side. The circuit network equivalent to Eq. (10) is discussed in Section 3.2.1.
The coefficients of the full projected load matrix pnL={p1L,…,pnL}T∈Rn×1 can be interpreted as the n components of the original excitation matrix f with respect to the basis vectors that are columns of the matrix KZn∈Rn×n as follows. Before truncation, using the full basis, we can express f as
f(t)=KZnpnL(t)=∑j=1nKzjpjL(t)=K∑j=1nzjpjL(t).
It can be verified that prL is indeed the excitation matrix (right-hand side) in Eq. (10) by premultiplying Eq. (11) by ZrTCK−1, with Zr=[w1L,…,wrL]=[z1,…,zr]∈Rn×r, and by using the C-orthonormality conditions for the column matrices in Zn∈Rn×n or in Zr∈Rn×r
Using Eq. (13) with r = n, we can also verify that
f=KZnpnL=KZn(ZnTCK−1f)=K(ZnZnTC)K−1f.
Since ZnTCZn=In⇒ZnTCZnZn−1=Zn−1⇒ZnZnTC=In, thus
f=KK−1f=f.
If we truncate the full basis of the n Lanczos trial vectors, and retain only r of those trial vectors, then:
f≈frL:=K∑j=1rzjpjL.
The error committed in representing f by a truncated series frL is thus:
erL:=∥f−frL∥2.Wilson–Yuan–Dickens (WYD) Method
The reduced-order ODE in Eq. (6) for the WYD method, (Appendix B) [24], takes the form:
Iry˙r+Kr∗yr=WrTf=:prW∈Rr×1,withKr∗:=WrTKWr∈Rr×r
where Wr:=[w1W,…,wrW]=[w1,…,wr]∈Rn×r contains the WYD trial vectors, with the projected load denoted by prW:=WrTf∈Rr×1 on the right-hand side of Eq. (18). With the full basis (no truncation), the excitation f can be represented in terms of the projected load pnW as follows:
f(t)=CWnpnW(t)=C∑j=1nwjpjW(t).
Premultiplying Eq. (19) by WrT∈Rr×n and using the C-orthonormality conditions
A truncated representation of the excitation matrix f using r ≤ n WYD (trial) vectors, and the resulting error on such representation can be expressed as
f≈frW:=C∑j=1rwjpjW;erW:=∥f−frW∥2.
For all four methods (Lanczos, WYD, Falk, Generalized Falk), the response dr(t) in the original (physical) coordinates can be computed from the response yr(t) in the reduced trial-vector basis denoted generically here as W_{r} (for all four methods) as follows:
dr(t)=Wryr(t)=∑i=1rwiyi(t),with r≤n.
In general, it is not necessarily true that the projected load along the first trial vector would have the highest magnitude, and then decreases along subsequent trial vectors (see Fig. 14 for example). The load participation factors βiL:=∥zipiL∥ for the Lanczos method, and βiW:=∥wipiW∥ for the WYD method (see Eq. (9)_{1} for the definition, which comes from Eq. (16) for the Lanczos method, and Eq. (24)_{1} for the WYD method) are not reliable indicators for the selection of the trial vectors to be retained for model reduction. In other words, a trial vector with a low value of load participation factor may actually play a key role in the response dr(t). This reason is why we introduce the response participation factors αiL:=∥ziyi∥ and αiW:=∥wiyi∥ (see Eq. (9)_{2} for the definition, which comes from Eq. (25) for the reconstructed solution), which can be used as indicators to evaluate the importance of the contribution of each trial vector to the response, and to guide the truncation of the trial-vector basis for model reduction.
Falk Method and Stability of Generalized Falk Method (GFM)
The GFM guarantees passivity and stability for efficient simulation of coupled field-circuit problems. We show that the GFM performs better in terms of loss of orthogonality compared to the original Falk method. Even though the Falk trial vectors and the GFM trial vectors are not identical, the final transformed systems have the same mathematical structure, i.e.,
Iry˙r+Kr∗yr=WrTf,
where Wr∈Rn×r, with r ≤ n, denotes the reduced basis formed by either the Falk trial vectors or the GFM trial vectors (Line 6 in Algorithm 2.1). As a result, the transformed stiffness matrices Kn∗ (Line 8 in Algorithm 2.1), even though having the same tridiagonal topology, are also different for different method. Since we normalize the Falk trial vectors and the GFM trial vectors with respect to the capacitance matrix C, i.e., WnTCWn=In (Line 7 in Algorithm 2.1), the projected loads are similar to the case of the WYD method, i.e.,
pnF:=(WnF)Tf,pnG:=(WnG)Tf,
with the excitation matrix f≡ f_{n} decomposed as follows:
f=CWnFpnF=CWnGpnG.
Similarly, the truncated representation of f, and the error committed due to the truncation are given respectively as
frF=C∑j=1rwjFpjF,frG=C∑j=1rwjGpjG.erF:=∥f−frF∥,erG:=∥f−frG∥.Applications of Discrete IGBTs and Model Participation Factors
We begin this section by providing a big picture of how discrete IGBTs are used in real-world engineering systems (Section 2.4.1), followed by an idealized model for electrothermal analysis of a discrete IGBT to deduce the desired participation factors to select modes with high participation factors for model-order reduction (Section 2.4.2).
Real-World Applications of IGBTs
Fig. 5 shows a single-chip power package used in power electronic circuits operating at high frequencies to achieve vastly improved performance in converters, either AC-to-DC or DC-to-DC^{9}
See Fig. 1 in the document ‘TRENCHSTOP^{TM} 5 IGBT (Footnote 4) in a Kelvin Emitter Configuration. Performance Comparison and Design Guidelines’. Application Note Revision 1.0, 2014-10-16, Online pdf, Internet archive. See also ‘Explanation of discrete IGBTs’ datasheets’. Application Note V1.0, 2015-09-18. Online pdf, Internet archive. ‘IGBT TRENCHSTOP 5 technology IGZ75N65H5 650 V high speed series 5th generation-Data sheet’, Online pdf, Internet archive.
. The chip is a discrete IGBT embedded inside a TO-247 package, which provides mechanical support and protection, heat sink for cooling, electrical connection (through the 4 pins in Fig. 5a) and isolation. The IGBT power package with 4 pins in Fig. 5 provides much better energy efficiency during switching compared to earlier models with 3 pins, allowing for higher swithching speed, lower operating temperature, and thus better overall system performance.
<italic>IGBT chip in TO-247 package</italic> (<xref ref-type="sec" rid="sec2_4_1">Section 2.4.1</xref>). (a) Single IGBT (<xref ref-type="fn" rid="fn1_4">Footnote 4</xref>) chip power package with through hole for mounting. Infineon (Brand name) IGZ75N65H5 (product code) device, with the TRENCHSTOP 5 IGBT chip on the Transistor-Outline (TO) 247 package with 4 pins. The dimension of the package, without the pins, is roughly <italic>W</italic> × <italic>L</italic> × <italic>H</italic> = 16 mm× 21 mm× 5 mm. (b) Material layers within the package. For the Tesla Model S electric vehicle, see <xref ref-type="fn" rid="fn1_3">Footnote 3</xref> on the use of IGBTs with 3 pins, and <xref ref-type="fig" rid="fig9">Fig. 9</xref> for the conceptual schematic of a half-bridge inverter. See <xref ref-type="fn" rid="fn1_9">Footnote 9</xref>
<italic>IGBT power packages attached to heat sink</italic> (<xref ref-type="sec" rid="sec2_4_1">Section 2.4.1</xref>). The middle device is a single IGBT chip with 3 pins embedded in a TO-247 package. Each IGBT power package is screwed to the heat sink via its through hole; see <xref ref-type="fig" rid="fig5">Figs. 5</xref> and <xref ref-type="fig" rid="fig7">7</xref>
<italic>General purpose IGBT stack for power converter/inverter</italic> (<xref ref-type="sec" rid="sec2_4_1">Section 2.4.1</xref>). (a) The aluminum heat sink with 8 discrete (3-pin) IGBT power packages (<xref ref-type="fig" rid="fig6">Fig. 6</xref>), which are screwed on via through holes, is cooled by air pushed through by the DC fan (black) on the right end of the heat sink. (b) The assembly in (a) is then mounted on a Power Circuit Board (PCB) for a converter/inverter device; see <xref ref-type="fn" rid="fn1_11">Footnote 11 [59</xref>]. (Figures reproduced with permission of the authors)
<italic>Tesla motor inverter</italic> (<xref ref-type="sec" rid="sec2_4_1">Section 2.4.1</xref>) (a) IGBT array for 2 switches, each with 14 discrete IGBTs in parallel (2 rows of 7 pairs of IGBTs). (b) Motor inverter containing an IGBT array shown in (a) for each of the 3 AC phases. See also <xref ref-type="fn" rid="fn1_3">Footnote 3</xref> for the big picture of how these components are used in the Tesla S electric vehicle, and <xref ref-type="fig" rid="fig9">Fig. 9</xref> for the generic schematic of a half-bridge inverter, having (a) as a specific implementation. <uri xlink:href="http://www.twinkletoesengineering.info/hybrid_car.htm">Twinkletoesengineering.info</uri>
Fig. 6 shows 5 different types of discrete IGBT power packages mounted on a heat sink^{10}
See ‘Extruded heat sink for PCB mounting’ by Astrel Website. Internet archive.
. Fig. 7 is a General Purpose IGBT Stack (GPIS) for use in a power converter or inverter^{11}
A converter transforms AC to DC, whereas an inverter transforms DC to AC.
“rated up to 20 kVA^{12}
The unit “kVA” or 1,000 Volt-Amperes (VA) is the “apparent power”, proportional to real power measured in “watt” (W) by a “power factor”. “The apparent power equals the product of root mean square voltage and root mean square current. In DC circuits, this product is equal to the real power in watts. Volt-amperes are usually used for analyzing AC circuits. The volt-ampere is dimensionally equivalent to the watt. VA rating is most useful in rating wires and switches (and other power handling equipment) for inductive loads”, Wikipedia version 13:20, 06 January 2021.
, and can be used in a variety of power electronic applications such as motor drives, energy storage systems, grid-tied and off-grid renewable energy systems” [59].
Fig. 8 shows a much more advanced motor inverter in the Tesla Model S electric vehicle that uses arrays of IGBTs (Fig. 8a) to convert DC from the battery pack to 3-phase AC, and thus a triangular cross section (Fig. 8b), to drive the induction motor that turns the wheels, Footnote 3.
Fig. 9 is a schematic of a half-bridge inverter in the Tesla patent [45], explaining how the IGBT array in Fig. 8a work to transform DC to the sine wave current for one of the three AC phases.
Fig. 10a shows a different geometry and material layering of the Tesla packaging from the patent application [44], compared to the packaging geometry and material layering in Fig. 5b, but the basic principles were the same.
<italic>Tesla half-bridge inverter</italic> (<xref ref-type="sec" rid="sec2_4_1">Section 2.4.1</xref>). (a) Schematic of a 3-pin IGBT (<xref ref-type="fn" rid="fn1_3">Footnote 3</xref> and <xref ref-type="fn" rid="fn1_4">Footnote 4</xref>, Gate G, Collector C, Emitter E), compared to that of a 4-pin IGBT in <xref ref-type="fig" rid="fig5">Fig. 5a</xref>. (b) Half-bridge inverter <bold>205</bold>, signal from driver <bold>305</bold> goes into the gate G of the upper switch (representing the top two rows, each with 7 IGBTs, in <xref ref-type="fig" rid="fig8">Fig. 8a</xref>) to turn it on or off, likewise for driver <bold>310</bold> for the lower switch (bottom two rows of IGBTs in <xref ref-type="fig" rid="fig8">Fig. 8a</xref>), output sine wave current <italic>I</italic><sub>1</sub> for phase 1 of a 3-phase current, current sensor <bold>350</bold> to feed output current <italic>I</italic><sub>1</sub> back to controller system (not shown). This subfigure is part of Fig. 2 in the Tesla patent [<xref ref-type="bibr" rid="ref45">45</xref>] (rotated by −90<sup>∘</sup> to match the 2 switches in <xref ref-type="fig" rid="fig8">Fig. 8a</xref>)<italic>Tesla IGBT packaging and idealized IGBT in TO-427 package</italic> (<xref ref-type="sec" rid="sec2_4_1">Sections 2.4.1</xref>, <xref ref-type="sec" rid="sec3_3">3.3</xref>, <xref ref-type="sec" rid="sec4_1">4.1</xref>). (a) <italic>Tesla IGBT packaging</italic>. (From top to bottom layer.) An encapsulant <bold>390</bold> encases a portion of the device; A silicon die <bold>310</bold> is sintered to a semiconductor device comprising layers <bold>330</bold>–<bold>350</bold>, below sintering layer <bold>320</bold>. The semiconductor <bold>340</bold> (IGBT or MOSFET) is interspersed between two copper cladding layers <bold>330</bold> and <bold>350</bold>, which lays on top of a direct bonded-copper substrate <bold>355</bold>. The whole structure is connected to the heat sink <bold>370</bold> via a sintering thermally-conducting silver layer <bold>360</bold>. Tesla inverter US patent application [<xref ref-type="bibr" rid="ref44">44</xref>]. <xref ref-type="fig" rid="fig5">Fig. 5b</xref> for TO-247 packaging. (b) <italic>Coupled transistor-thermal model</italic>. IGBT chip on heat sink, compared to <xref ref-type="fig" rid="fig5">Figs. 5b</xref> and (a), used for electrothermal analysis and proof-of-concept for GFM. See also <xref ref-type="fig" rid="fig29">Fig. 29</xref> in <xref ref-type="sec" rid="sec3_3">Section 3.3</xref> on SABER implementationParticipation Factors of IGBT Electrothermal Model
We use here an idealized 2-D thermal model problem shown in Fig. 10b to illustrate the properties of the projected load p_{i}, the load participation factor β_{i}, and the response participation factor α_{i} for each of the four methods discussed in the previous section: Lanczos (Figs. 11–13), WYD (Figs. 14–16), Falk (Figs. 17–19), and GFM (Figs. 20–22). The finite-element (FE) discrete model has 368 triangular elements and 215 nodes.
<italic>Projected load</italic> <inline-formula id="ieqn-58"><mml:math id="mml-ieqn-58"><mml:msubsup><mml:mi>p</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>L</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. Lanczos vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>)<italic>Load participation factor</italic> <inline-formula id="ieqn-59"><mml:math id="mml-ieqn-59"><mml:msubsup><mml:mi>β</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>L</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. Lanczos vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>)<italic>Response participation factor</italic> <inline-formula id="ieqn-60"><mml:math id="mml-ieqn-60"><mml:msubsup><mml:mi>α</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>L</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. Lanczos vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>): Only first 120 trial vectors among 215 trial vectors are shown<italic>Projected load</italic> <inline-formula id="ieqn-61"><mml:math id="mml-ieqn-61"><mml:msubsup><mml:mi>p</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>W</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. WYD vector number</italic> (<xref ref-type="sec" rid="sec2_3_2">Sections 2.3.2</xref>, <xref ref-type="sec" rid="sec2_4_2">2.4.2</xref>): Only the first 120 trial vectors among the 215 trial vectors are shown
The apparent randomness in the projected loads p_{i} and in the load participation factors β_{i} in the Lanczos method (Figs. 11–13) and in the original Falk method (Figs. 17–19), as compared to the smoother projected loads in the WYD method (Figs. 14–16), and in the GFM (Figs. 20–22), is clearly a result of the random starting vectors used in the Lanczos process.
For the WYD method, the projected loads p_{i} and the load participation factors β_{i} start at low values for the static solution, then increase for subsequently generated trial vectors up to some maximum values, then decrease almost monotonically for all subsequent trial vectors. From approximately the 40th trial vector on, the values of p_{i} and of β_{i} in the WYD method are practically zero.
For the GFM, the projected loads p_{i} and the load participation factors β_{i} have significant values concentrated mostly in the first few trial vectors (less than 10).
<italic>Load participation factor</italic> <inline-formula id="ieqn-62"><mml:math id="mml-ieqn-62"><mml:msubsup><mml:mi>β</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>W</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. WYD vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>). Only first 120 trial vectors among 215 trial vectors are shown
The response participation factors α_{i} should be used as indicators to select the participating trial vectors for model reduction. For the Lanczos and WYD methods, the response participation factors α_{i} (Fig. 13 for Lanczos, and Fig. 16 for WYD) have significant values concentrated in the first few trial vectors, as opposed to the more spread-out projected loads p_{i} (Fig. 11 for Lanczos, and Fig. 14 for WYD) and load participation factors β_{i} (Fig. 12 for Lanczos, and Fig. 15 for WYD). On the other hand, the response participation factors α_{i} for both the Falk method (Fig. 19) and the GFM (Fig. 22) are spread out over some 50 trial vectors, while the projected loads p_{i} (Fig. 17 for Falk, and Fig. 20 for GFM) and the load participation factors β_{i} (Fig. 18 for Falk, and Fig. 21 for GFM) are more concentrated in the first few trial vectors. It should also be noted that the response participation factor α_{i} in the GFM decrease monotonically and faster as the trial-vector number increases, compared to those in the Falk method in its original form.
<italic>Response participation factor</italic> <inline-formula id="ieqn-63"><mml:math id="mml-ieqn-63"><mml:msubsup><mml:mi>α</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>W</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. WYD vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>). Only first 50 trial vectors among 215 trial vectors are shown<italic>Projected load</italic> <inline-formula id="ieqn-64"><mml:math id="mml-ieqn-64"><mml:msubsup><mml:mi>p</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>F</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. Falk vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>)
Finally, because the transformed matrices have the tridiagonal structure, the reduced-order models in the Lanczos method and in the GFM require less time to solve than that in the WYD method (Fig. 23). The accuracy of these methods are compared in Fig. 24, where the error norm was plotted in terns of the number of trial vectors retained. The WYD method converges fastest, followed by the GFM, then the Lanczos method, and the Falk method. The GFM method converges faster than the Lanczos method, e.g., the error for the GFM with 10 retained trial vectors is equivalent to that of the Lanczos method with 30 to 40 retained trial vectors, as indicated by the red line in Fig. 24. A comparison of Figs. 23 and 24, together with all previous Fig. 11 to Fig. 22, thus reveals a clear niche for the GFM.
<italic>Load participation factor</italic> <inline-formula id="ieqn-65"><mml:math id="mml-ieqn-65"><mml:msubsup><mml:mi>β</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>F</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. Falk vector number</italic>. (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>)<italic>Response participation factor</italic> <inline-formula id="ieqn-66"><mml:math id="mml-ieqn-66"><mml:msubsup><mml:mi>α</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>F</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. Falk vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>). Only the first 160 trial vectors among the 215 trial vectors are shown<italic>Projected load</italic> <inline-formula id="ieqn-67"><mml:math id="mml-ieqn-67"><mml:msubsup><mml:mi>p</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>G</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. GFM number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>). Only the first 50 trial vectors among the 215 trial vectors are shown<italic>Load participation factor</italic> <inline-formula id="ieqn-68"><mml:math id="mml-ieqn-68"><mml:msubsup><mml:mi>β</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>G</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. GFM vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>): Only first 50 trial vectors among 215 trial vectors are shown<italic>Response participation factor</italic> <inline-formula id="ieqn-69"><mml:math id="mml-ieqn-69"><mml:msubsup><mml:mi>α</mml:mi><mml:mrow><mml:mi>i</mml:mi></mml:mrow><mml:mrow><mml:mi>G</mml:mi></mml:mrow></mml:msubsup></mml:math></inline-formula> <italic>vs. GFM vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>): Only first 60 trial vectors among 215 trial vectors are shown
Also note that for the heat conduction problem the response participation factors α_{i} vary with time, but their magnitude relative to each other remain the same, see Fig. 13 (Lanczos), Fig. 16 (WYD), Fig. 19 (Falk) and Fig. 22 (GFM).
Two Model-Reduction Strategies
Based on the algorithms presented above, several model-reduction strategies using a combination of the above four methods can be introduced. Two strategies are suggested below.
<italic>Simulation time for various reduced models vs. trial-vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>)<italic>Comparison of error norm vs. trial-vector number</italic> (<xref ref-type="sec" rid="sec2_4_2">Section 2.4.2</xref>)
The first strategy consists of (i) using the WYD method to obtain a reduced-order model by truncating the generation of the WYD trial vectors, followed by (ii) the use of the GFM to transform the resulting reduced-order model to a very simple form (identity capacitance, tridiagonal conductance, Fig. 25), as discussed in the previous sections. Actually, we can even obtain a further reduced-order model by truncating the generation of the GFM trial vectors.
Recall that such truncation can be based on the use of the response participation factor α_{i} discussed above. At the end of this two-stage model-reduction strategy, we obtain a reduced-order model with an identity capacitance matrix and a tridiagonal conductance matrix, which is represented by a simple 1-D equivalent circuit, as discussed in the next sections.
<italic>Two-step model reduction strategy</italic> (<xref ref-type="sec" rid="sec2_5">Section 2.5</xref>). <italic>Strategy 1:</italic> The first Step consists of using the WYD method, whereas the second Step relies on the GFM method. <italic>Strategy 2:</italic> The first Step is based on the Lanczos method, and the second Step on the original Falk method
The second model-reduction strategy is based on the Lanczos and the original Falk methods. The first strategy based on the WYD and GFM methods is more stable than the second strategy based on the Lanczos and Falk methods, especially if the eigenspectrum has close eigenvalues.
Coupled Circuit-Thermal Simulation with Transformed SystemsTransistor-Thermal Model
In our case-study of electrothermal simulation of circuits with semiconductor devices (e.g., discrete IGBTs or MOSFETs), the circuit network model, the semiconductor device models, and the equivalent thermal circuit network modeling the heat sinks are all coupled together, and solved in a circuit simulator. The electrothermal semiconductor models use the instantaneous device temperature (temperature at the silicon chip surface d_{node}) to evaluate the temperature-dependent properties of silicon and the temperature-dependent model parameters. These temperature-dependent values are then used by the physics-based semiconductor device models to describe the instantaneous electrical characteristics and the instantaneous dissipated power. In [22 , 54], we proposed a methodology to develop equivalent thermal circuit networks based on a finite-element discretization of the heat-diffusion equation over the domain of the heat sink. Once a thermal network component is developed, and connected to the electrical networks of power electronic systems to provide complete electrothermal models that can be conveniently used in any circuit simulator, it can be used over and over again in many other electrothermal circuit simulation problems.
Equivalent Circuits for Thermal Part
The methodology developed in [15, 22] could lead to complex thermal networks that present some challenge to implement in circuit simulators such as SABER. If we could transform the original discrete field model into a simple equivalent model, having simple circuit representation that can be easily implemented in any circuit simulator, then we would have considerably simplified the building of the models in coupled field-circuit simulation problems. Indeed, as already mentioned above, the coordinate transformation methods developed in the previous sections can transform any complex thermal model into an extremely simple (1-D) model represented by the “1-D” circuit shown in Fig. 26 (Lanczos) and Fig. 27 (GFM).
<italic>Equivalent circuit by Lanczos method</italic> (<xref ref-type="sec" rid="sec3_2_1">Section 3.2.1</xref>)<italic>Equivalent circuit by GFM</italic> (<xref ref-type="sec" rid="sec3_2_2">Section 3.2.2</xref>)
Remark 3.1 As presented in [31], the prescribed initial condition Eq. (7) in the new coordinates is computed as follows:
d(0)=Wy(0)⇒y(0)=WTCd(0),
which can be used as input into the circuit simulator SABER, or other any other circuit simulators, as initial nodal voltages. From the physical standpoint, the initial condition can be obtained as the static solution
Kd(0)=f(0)=Ic
where the ambient temperature T_{a} in the forcing term f(0)=Ic dictates the initial condition. In a circuit simulator, the initial condition as represented by Eq. (32) above is obtained through an initial steady-state solution, prior to the transient solution stage. The modeling of the initial condition in Eq. (5) in our thermal simulation is based on the input initial voltage at each node. From a physical standpoint, before any simulation, SABER solves the steady-state solution first as initial condition. The user then only needs to input the ambient temperature, but does not need to impose the initial voltage at each node. For a general mathematical heat-conduction problem, this is a particular case to bypass the process of reading in an initial condition.
Equivalent Circuit by Lanczos Method
Based on the method for constructing equivalent circuit networks from symmetric matrices in [54], a symmetric global matrix can be decomposed into several elemental matrices, and each elemental matrix represented by a circuit component. A symmetric capacitance matrix C can be realized by a capacitor network, whereas a symmetric conductance matrix K can be realized by a resistor network. After a Lanczos coordinate transformation, the new capacitance matrix T_{r} becomes tridiagonal, and the conductance matrix K^{*} becomes an identity matrix, Eq. (10) and Appendix A. As a result, in the 1-D equivalent circuit, there is one capacitor in parallel with one resistor between each node and the ground. The capacitance C_{ii} = c_{i, i −1} +c_{ii} + c_{i, i+1}, where c_{i, i −1}, c_{i, i −1}, and c_{i, i −1} are the coefficients of T_{r}. The resistor R_{ii} = 1. Between two adjacent nodes i and j, there is a “mutual capacitance” C_{ij}, whose values are the off-diagonal coefficients of T_{r}, as shown in Fig. 26 (Lanczos), where I_{p} corresponds to the input power, and I_{c} the heat-convective boundary conditions, see Eq. (6).
Equivalent Circuit by GFM
Unlike the Lanczos method, the GFM yields an identity capacitance matrix C^{*} and a tridiagonal conductance matrix K^{*}, Eq. (26) and Algorithm 2.1Instead of having the capacitors between two adjacent nodes, we now have the resistors, as shown in Fig. 27. As a result, we have
In addition, there is no capacitor between adjacent nodes since C^{*} is an identity matrix. The resistor between two adjacent nodes is
Rij=−1kij.
Again, in Fig. 27, I_{p} corresponds to the input power, and I_{c} the heat-convective boundary condition, see Eq. (6).
Implementation in Circuit Simulators
We have written a MATLAB code to automatically generate element templates^{13}
A template is the mathematical description of a subsystem and is contained in a text file. The characteristic equations implemented in a template can be any combination of linear or non-linear algebraic or differential equations.
of the 1-D field networks^{14}
The field networks are the thermal networks in the case study of electrothermal simulations.
in Fig. 26 (Lanczos transformation) and Fig. 27 (Generalized Falk transformation) for the SABER circuit simulator using the MAST Hardware Description Language^{15}
The MAST Hardware Description Language is the modeling language, originally developed by Analogy, Inc., that allows users to create device simulation models that can be defined in mathematical terms, in any technology, using the native equations and units, without having to resort to electrical equivalents [60].
. Models of field packages can be built from these element templates.
A 2-D thermal problem employed in [22] was used to test the simplified circuit networks. Fig. 28 shows the flowchart of our MATLAB code to implement a thermal component model. A finite-element code written in MATLAB for 2-D thermal analysis will generate a 2-D mesh, perform finite element analysis and coordinate transformation by the methods we presented above. The transformed system (transformed capacitance matrix C^{*}, transformed conductance matrix K^{*} and transformed load vector f^{*}) and the trial vectors (W) generated by the coordinate transformation methods will be used to construct the equivalent circuit model (SABER netlist file) and the above mentioned element templates.
Since SABER electrothermal model was constructed to have only one thermal terminal, therefore, we will define only one thermal terminal in our thermal component model by defining the nodal temperature (across variable) and the power input (through variable) at that node as system variables. In general, after coordination transformation, the input power supply in the original idealized physical model as shown in Fig. 10b will distribute out to every node in the transformed model. In other words, the transformed excitation (f^{*}) will have nonzero component in all directions. Therefore, CCCS’s (Current Controlled Current Sources) in the equation section^{16}
The equation section contains the terminal (connection point) equations of the model, which are often equivalent to the characteristic equations. The relationships involving the through and across variables must be defined in this section [60].
of our element templates are used for the power loss from the thermal terminal of the electrothermal model (excitation in the physical coordinate) to control excitations in all other nodes (transformed coordinate).
<italic>Implementation flowchart</italic>. (<xref ref-type="sec" rid="sec3_3">Section 3.3</xref>). MATLAB code flowchart to implement a thermal component model. In the SABER circuit simulator, the Current Controlled Current Sources (CCCS) are used in the element template to connect to the thermal terminal of electrothermal models, and the Voltage Controlled Voltage Sources (VCVC) to obtain the nodal temperatures in the transformed coordinate system, for later recovery back into the physical coordinate system
A typical SABER thermal component template is written as follows:
template heat2d_mesh1 tnode
thermal_c tnode
{
thermal_c y1, y2, y3, y4, y5,… # local connection point
number p1 = 3.89146e + 01, h1 = 2.30179e-01,
p2 = 4.03557e + 02, h2 = 1.12042e − 01,
…
w1 = 3.89146e + 00, # trial vector component
w2 = 4.43737e + 01,
…
var i ip # controlling current
equations
{
p(tnode) += ip
ip: tc(tnode) = w1*tc(y1)+ w2*tc(y2)+…# (VCVS)
p(y1) −= p1*ip + h1 # (CCCS)
p(y2) −= p2*ip + h2 # (CCCS)
…
}
c_th.c1 p:y1 m:0 = 1, 6.96555e + 00
g_th.g1 p:y1 m:0 = 3.75182e − 01
…
netlist which describes the transformed 1-D RC network
}
where : p(tnode) is the across variable at the power input node that will be defined as template thermal terminal, which is physically the temperature at that node. The variable : ip is the through variable at the defined terminal. The variable : ip should also be defined as system variable since all power input at the nodes of the transformed coordinates will be controlled by the CCCS’s via the input power : ip (variable P in Eq. (3), Table 2) at the thermal terminal of the physical coordinates as shown in Fig. 29. From Eq. (18) and Eq. (6), the transformed force can be expressed as
WTf=WTIp+WTIc=p+h.
The above expression is given in matrix form as follows:
The controlling input power : ip (“P” in Eqs. (3) and (34), Table 2) is used in the CCCSs as follows:
WTf=[p1p2…pm]P+[h1h2…hm]=p+h.
The nodal temperature of the defined thermal terminal will be recovered from the transformed coordinates to the physical coordinates by a VCVS written in the above equation section of the element template. Basically, it is the coordinate-transformation equation at that node:
dnode=∑j=1mwjyj=w1y1+w2y2+…+wmym
where y={y1,y2,…,ym}T is the temperature at the defined thermal terminal in the transformed coordinates.
Table 2 gives the correspondence between the system variables used in our formulation and those used in our MATLAB or SABER implementation.
<italic>System variables</italic> (<xref ref-type="sec" rid="sec3_3">Section 3.3</xref>). Correspondence between symbols used in formulation and variable names in MATLAB/SABER implementation, with related equation numbers. Nodal temperature in physical coordinates. Modal temperature in transformed coordinates
System variables
Formulation
MATLAB/SABER
Eq. No.
Input power
P
: ip
Eq. (3)
Nodal temperature
d
: tc
Eq. (6)
Trial vectors
W
: w
Algorithm 2.1, Eq. (26)
Modal temperature
y
: y
Eqs. (25), (26)
<italic>Implementation of the 1-D equivalent circuit in SABER</italic> (<xref ref-type="sec" rid="sec3_3">Section 3.3</xref>). See <xref ref-type="fig" rid="fig10">Fig. 10b</xref> for the idealized physical domain, <xref ref-type="fig" rid="fig30">Fig. 30</xref> for the overall schematic of equivalent electrothermal circuit network, and <xref ref-type="table" rid="tab1_2">Table 2</xref>Connection to Other Circuit Components
The input power to the thermal system is supplied from the thermal terminal of the electrothermal semiconductor model via the defined thermal terminal of the thermal component model template. The thermal terminal of the semiconductor electrothermal model is thus connected to only the thermal terminal of the thermal component model. By defining the power flowing into the thermal terminal of the thermal component model as system variable, the calculated nodal temperature is dynamically coupled to the semiconductor model.
In Fig. 30, the original complex equivalent thermal circuit network is first transformed by the trial vector basis as d = Wy, Eq. (25), into a 1-D simple equivalent circuit network by our proposed methodology, and the transformation in SABER is implemented by the CCCS’s and VCVS as explained above. The nodal temperature in the transformed coordinates is solved by the SABER built-in efficient numerical solver. The solved nodal temperatures in y in the transformed coordinates is then recovered back to the physical coordinates d using the same trial vectors. Fig. 29 shows the schematic representation of the implementation of the 1-D equivalent circuit into SABER.
<italic>Schematic of equivalent electrothermal circuit network</italic> (<xref ref-type="sec" rid="sec3_3">Section 3.3</xref>). The thermal part becomes a simple 1-D circuit network after coordinate transformation. See <xref ref-type="fig" rid="fig29">Fig. 29</xref> for implementation in SABER circuit simulatorNumerical Examples with Model Reduction and Transformation
In this section, we apply the model reduction strategies presented above to extract the reduced-order circuit models for the field component in coupled field-circuit problems to illustrate the efficiency of the proposed methodology. Simulation times for both full-order models and for the reduced-order models are presented for comparison.
Coupled Transistor-Thermal Models, Verification
The transistor-thermal models obtained from our model-reduction strategies can be verified by considering the 2-D problem in Fig. 10b, showing an Silicon (Si) chip (IGBT or MOSFET) on an idealized TO-247 package. The finite element mesh of this device contains 368 linear triangular elements and 215 nodes.
Used in the analysis are both (1) the full-order model and (2) the reduced-order model obtained by using a combination of WYD method followed by the GFM or the original Falk method. To generate the reduced-order transformed model, a tolerance of TOL = 1E −3 based on the truncation of projected load p defined in Section 2.4, was prescribed in the WYD algorithm, and automatically resulted in a model with 51 nodes, reduced down from 215 nodes. A power of P = 10 W is inputted at the top of the transistor chip, and is distributed uniformly over an area A = 0.1 cm^{2} (parameter a in Table 4). The circuit simulator SABER and MATLAB were used to implement all developed models, with details of how the full-order model was implemented in SABER described in [54].
Fig. 31 shows a comparison of the transient temperature obtained using the full-order model with m = 215 unknowns, and using a reduced-order model with different numbers r = 11, 31, 41, 51 of retained GFM trial vectors. There is a complete agreement between the two models (full and reduced order) for the temperature at the top of the Si transistor chip.
<italic>Transient temperature in a pn junction</italic> (<xref ref-type="sec" rid="sec4_1">Section 4.1</xref>). GFM reduced-order model <italic>vs</italic>. full-order model (left). GFM method: Temperature error <italic>vs</italic>. number of GFM vectors (right)Full-Bridge Converter, Electrothermal Simulation
Considerable attention has been focused on switched-mode technology to regulate power supply because it is possible to achieve lossless power conversion to meet the continuing and increasing demand of power electronic devices with reduced size and weight, and with increased efficiency. To turn on and off the flow of energy, and thus achieve regulation, duty-cycle control is employed in switching elements. It has the added advantage when applied to off-line applications of giving significant size reduction in the voltage transformer and energy storage elements, The size of voltage transformer and energy-storage elements in off-line applications can be significantly reduced using switched-mode technology. Such size reduction is another advantage of this technology, in addition to regulating energy flow [48].
The full-bridge converter is typically used in switching power supplies at power levels of approximately 750 W and greater. A full-bridge buck converter^{17}
The word “buck” came from the fact that the inductor in the circuit always “bucks” or acts against the input voltage. See ‘Buck DC/DC Converters,’ Power Supply Technology, Mouser Electronics Website.
with four discrete IGBTs and an isolated transformer, shown in Fig. 32, is simulated with the circuit simulator SABER [22]. The voltages that drive the four discrete IGBTs in the full-bridge buck converter in Fig. 32 are shown in Fig. 33, where the on-time is T_{on} = 15 μs, the dead-time Δ = 10 μs when all switches are turned off, and the switching-time (or the impulse period) T_{s} = 50 μs. In the first impulse, lasting for on-time T_{on}, i.e., the interval 0 < t < T_{on}, the two discrete IGBT1 and IGBT4 allow current to pass through, while the voltage Vin is inputted into the primary winding of the transformer. The secondary winding of the transformer is center-tapped^{18}
See the definition of “center-tapped transformer” as a design to allow for two separate voltages with a common connection in Multiple-winding transformers.
, with each half having a voltage of nVin, with n being the ratio of the primary turn over the secondary turn. The positive side of the potential is indicated by a polarity mark (dot). The diode D_{5} in Fig. 32 is placed to allow the clockwise current flow, from positive pole (indicated by the polarity mark) to negative pole, but prevents the counterclockwise current flow, and is thus called forward-biased. Even though the diode D_{6} appears to be in the same direction of diode D_{5}, it only allows the counterclockwise current to flow, and is called reverse-biased. Across the filter, the input voltage is therefore nVin, with the filter inductor current i(t) flowing through the diode D_{5}. During the dead-time in the interval T_{on} < t < T_{on} + δ, all four discrete IGBTs are switched off, and thus the transformer voltage is Vo = 0 (output voltage). During the time interval T_{on} + δ < t < T_{s} − δ, the discrete IGBT2 and IGBT3 and the diode D_{6} allow current to flow through.
<italic>Full-bridge converter</italic> (<xref ref-type="sec" rid="sec4_2">Section 4.2</xref>). As a switching device, a full-bridge converter uses four discrecte IGBTs. See <xref ref-type="fig" rid="fig33">Fig. 33</xref> for the time histories of the driving voltages <italic>Vg</italic>1 to <italic>Vg</italic>4 for <italic>IGBT</italic>1 to <italic>IGBT</italic>4, respectively. See also <xref ref-type="table" rid="tab1_3">Tables 3</xref> and <xref ref-type="table" rid="tab1_4">4</xref> for the parameters used<italic>Full-bridge converter</italic> (<xref ref-type="sec" rid="sec4_2">Section 4.2</xref>). Waveforms of the driving voltages for the full-bridge converter in <xref ref-type="fig" rid="fig32">Fig. 32</xref> as a switching device. The voltage time history <inline-formula id="ieqn-73"><mml:math id="mml-ieqn-73"><mml:mi>V</mml:mi><mml:mi>g</mml:mi><mml:mn>1</mml:mn><mml:mrow><mml:mo>(</mml:mo><mml:mi>t</mml:mi><mml:mo>)</mml:mo></mml:mrow><mml:mo>=</mml:mo><mml:mi>V</mml:mi><mml:mi>g</mml:mi><mml:mn>4</mml:mn><mml:mrow><mml:mo>(</mml:mo><mml:mi>t</mml:mi><mml:mo>)</mml:mo></mml:mrow></mml:math></inline-formula> are for the two discrete IGBT1 and IGBT4, and <inline-formula id="ieqn-74"><mml:math id="mml-ieqn-74"><mml:mi>V</mml:mi><mml:mi>g</mml:mi><mml:mn>2</mml:mn><mml:mrow><mml:mo>(</mml:mo><mml:mi>t</mml:mi><mml:mo>)</mml:mo></mml:mrow><mml:mo>=</mml:mo><mml:mi>V</mml:mi><mml:mi>g</mml:mi><mml:mn>3</mml:mn><mml:mrow><mml:mo>(</mml:mo><mml:mi>t</mml:mi><mml:mo>)</mml:mo></mml:mrow></mml:math></inline-formula> are for IGBT2 and IGBT3
For the full-bridge buck converter in Fig. 32, the parameters were selected as shown in Table 3 so to ensure that a thermal run-away^{19}
When an electronic equipment continues to generate heat at a rate faster than the heat can be dissipated, a phenomenon called “thermal run-away”, the equipment would often fail or a fire would break out.
, which cannot be predicted using a conventional circuit simulator with fixed temperature, would happen.
<italic>Full-bridge converter</italic> (<xref ref-type="sec" rid="sec4_2">Section 4.2</xref>). Parameters. See <xref ref-type="fig" rid="fig32">Fig. 32</xref>
Parameter
Symbol
Value
Input voltage
Vin
100 V DC
Primary to secondary turn ratio
n = n_{1}/n_{2}
2
Gate resistance
R_{gi}, i = 1,… , 4
100 Ω
Low-pass filter inductance
L_{o}
300 μH
Low-pass filter capacitance
C_{o}
100 μF
Output resistance
R_{o}
5 Ω
<italic>SABER IGBT model</italic> (<xref ref-type="sec" rid="sec4_2">Section 4.2</xref>). Parameters. See <xref ref-type="fig" rid="fig32">Fig. 32</xref> and <xref ref-type="fn" rid="fn1_20">Footnote 20</xref>
Parameter
Code variable
Input
Value or unit
High level excess carrier lifetime
tauhl
7.1u
7.1× 10^{−6} s
Temperature exponent for tauhl
tauhlexp
1.5
–
Metallurgical base width
wb
10m
10× 10^{−3} cm
Base doping concentration
nb
2.0e14
cm^{−3}
Device active area
a
0.1
cm^{2}
Gate drain overlap active area
agd
0.05
cm^{2}
Emitter electron saturation current
isne
10f
10× 10^{−15}A
Temperature exponent for isne
isnetexp
0.5
–
MOSFET channel threshold voltage
vt
5.0
V
Temperature coefficient for vt
vttco
−0.009
–
Intrinsic anode series resistance
rs
0.01
Ω
Parameter for transconductance reduction
theta
0.01
V^{−1}
due to transverse electric field in MOSFET
Temperature exponent for theta
thetatexp
0
–
MOSFET empirical factor representing the ratio of
kf
2.0
–
kp in triode region to that in saturation
(=kp, linear/kpt)
Temperature exponent for kf
kftexp
0
–
MOSFET channel transconductance in
kp
0.25
AV^{−2}
saturation region
Temperature exponent for kp
kptexp
1.5
–
Gate to source capacitance
cgs
1n
10^{−9} F
Gate drain oxide capacitance
coxd
2n
2× 10^{−9} F
Gate drain overlap depletion threshold
vtd
0
V
Temperature coefficient for vtd
vtdtco
0
–
Junction curvature factor
bvf
1.0
–
Temperature coefficient for bvf
bvftexp
0.35
–
Avalanche multiplication exponent
bvn
4.0
–
Temperature coefficient for bvn
bvntexp
0
–
Nominal temperature used for parameters
tnom
27
deg. C (300^{∘} K)
Temperature exponent for mobilities
alpha
2.54
–
Minimum slope for MOSFET current
gmin
1p
10^{−12} A
In the circuit simulator SABER, the parameters for the electrothermal IGBT model are shown in Table 4, based on the values suggested in [61 , 62]^{20}
See also the default values for the parameters of the n-type IGBT model in the Mathworks Spice circuit simulator. Website. Internet archive. SABER numerical inputs use the metric prefix as suffix, e.g., “10f” stands for “10 femto = 10×10^{−15}”; see ‘Metric prefix’, Wikipedia version 08:11, 19 February 2021. Among the four references—[61–64]—only the conference paper [61] mentioned the coefficients K_{p} and K_{f} with the corresponding code variables kp and kf in Table 4.
. The simulation results of the short-circuit test of the IGBT model using the parameters in Table 4 are similar to those in [63]. The GFM transformed model used here also established the efficiency in coupled transistor-thermal simulations.
Remark 4.1 Since an IGBT is a cross between a bipolar junction transistor (BJT)^{21}
See ‘Bipolar junction transistor’, Wikipedia version 15:45, 26 February 2021.
and a MOSFET, there are parameters for MOSFET in Table 4, even though this table is for IGBT parameters, [62]. Before the 1970s, when MOSFET was invented, BJT was the only device used in power electronics. BJT requires high base current to turn on, and slow to turn off, and is subjected to thermal runaway. Unlike BJT, for which the switching is current controlled, MOSFET is voltage controlled, and can limit or stop thermal runaway, and thus became the go-to device in power switch design. IGBT, introduced in the 1980s, has the high-current handling of BJT, and the ease of voltage-control of MOSFET. Generally, IGBT is chosen for low-frequency ( < 20 KHz) and high-voltage ( > 1000 V) applications, whereas MOSFET is chosen for high-frequency ( > 200 KHz) and low-voltage ( < 200 V) applications. In between, either device can be used, depending on application requirements, such as cost, size, speed, thermal specifications^{22}
See ‘IGBT or MOSFET: Choose wisely,’ by C. Blake and C. Bull, International Rectifier. Online pdf. Internet archive.
.
The reduced-order model obtained by using a combination of WYD method followed by the GFM^{23}
The same results were obtained from the original Falk method compared to those from the GFM.
, produced a speed-up factor of close to 9 for the transient-simulation time, compared to the original 2-D-equivalent-circuit thermal full model, as shown in Table 5.
<italic>Full-bridge converter</italic> (<xref ref-type="sec" rid="sec4_2">Section 4.2</xref>). Comparison of simulation time. The original model has 215 unknowns (nodal temperatures). The reduced-order model (ROM) has 51 unknowns (modes). See <xref ref-type="table" rid="tab1_6">Table 6</xref> for the Voltage Regulator Module (VRM)
Simulation model
SABER transient analysis
Speed-up factor
Full model
1190 s
GFM/WYD
136 s
1190/136 = 8.75
<italic>Full bridge converter</italic> (<xref ref-type="sec" rid="sec4_2">Section 4.2</xref>). Output power and temperature <italic>vs</italic>. time. Full-order model (left) <italic>vs</italic>. GFM reduced-order model (right)
The temperatures Tj1 and Tj2 of the discrete IGBT1 and IGBT2, with heat sink provided by the TO-247 packages, vs. simulation time, together with the time history of the converter output power (V_{o} * I_{o}), are given in Fig. 34, with and without thermal effects. The result without thermal effects is obtained from the simulation where the Si-chip temperature is fixed to be the ambient temperature. The results clearly show the output power difference between the simulation with thermal effects and that without thermal effects. Such output-power difference would be neglected if we did not incorporate the electro-thermal coupling in the power switching devices into the simulation.
Remark 4.2 In power electronics, the meaning of the term “load” could be confusing. It is best to think of two types of loads: Current load and resistance load. In general, the term “load” is often used to designate “current load”, which is related to output power. Thus “higher load” often means “higher current load”, and thus higher output power. Further, with the relation V = RI, increasing current load corresponds to decreasing resistance load, under constant voltage.
Remark 4.3 A discrete IGBT can be modeled by a resistor, an inductor, and a DC voltage source in series, together with a current switch, controlled by a logical signal, gate on g > 0 or gate off g = 0^{24}
See, e.g., ‘IGBT: Implement insulated gate bipolar transistor’, Simscape Electrical version R2020b, Mathworks, Website.
. In general, for large output voltage, the voltage drop from input to output in IGBT is small. The reason for a small voltage drop across an IGBT device is because this device is roughly equivalent to a voltage source (in series with a resistor and an inductor). Smaller voltage drop in IGBT means larger voltage output which means less power losses. That is why IGBT devices are usually used in higher power applications.
Remark 4.4 Switched-mode semiconductor devices (i.e., IGBT and MOSFET) are employed in high-efficiency power converters. When a semiconductor device operates in the off state, its current is zero and hence its power dissipation is zero. When the semiconductor device operates in the on (saturated) state, its voltage drop is small and hence its power dissipation is also small. In either event, the power dissipated by the semiconductor device is low [48].
Simulation of a Voltage Regulator Module (VRM)
The ability to work at lower voltage and higher current is required of modern microprocessors to meet the demand of faster and more efficient data processing. Moreover, each modern processor is packed with increasingly more devices, and operates at increasingly higher frequencies. Voltage regulator modules (VRMs) are special on-board power supply modules that minimize the effects of parasitics in interconnections, and can provide highly accurate regulation of supply voltage, as a centralized power system cannot realize these goals. In microprocessors capable of operating at lower voltage and higher current, a VRM located near the load on the motherboard is required [65]. The DC-to-DC converter in a distributed power system produces an intermediate voltage appearing at the computer backplane. Locally regulated low voltages are generated by high-density DC-to-DC converters on each card. Because there is not a lot of space on a crammed motherboard, it is important that on-board converter modules, such as the VRM, be of high power density and high efficiency. A serious design challenge is to have the power conversion perform at high switching frequency to provide fast transient response and to satisfy the requirements of high power density and high efficiency. Each Si chip has an operating output-voltage range, with some tolerance outside of which the chip would fail (blow up). Moreover, there is a voltage range, inside the operating-voltage range, within which the Si chip would be most efficient. For these reasons, the regulated voltage is designed to be fixed.
Fig. 35 shows the circuit schematic of a synchronous rectifier buck-converter used as VRM. The load, represented by the resistor R (i.e., the microprocessor) is powered up from a power supply with regulated output voltage Vg (or Vin). The circuit simulator SABER with our own equivalent circuit model for the idealized TO-247 thermal field component shown in Fig. 10b was employed to simulate the VRM in Fig. 35 [66] with the driving voltage waveforms^{25}
The results for the VRM are simulations designed to show the difference between model with thermal effects and that without thermal effects. The parameters selected may not be practical. Too small input voltage did not make a significant difference between the cases with and without thermal effects. We changed these parameters to push the electro-thermal model to temperature limits. The input voltage 60 V is for the buck-converter circuit as a whole, not only for the MOSFETs. For more on VRM, see also [67].
shown in Fig. 36, and with the operating frequency set at 500 kHz. The major VRM components were (Fig. 35 from left to right): Input voltage V_{in} = V_{g} = 60 V, MOSFET 1 main switch SW-2XIRF7811, MOSFET 2 synchronous rectifier SR-4XIRF7811, inductor L = 500 nH, output-filter capacitor C = 100 μF, load resistance R = 15 Ω^{26}
Shown in Fig. 34 were the junction temperature results using the transformed model. These results were the same for the full model.
.
<italic>Voltage Regulator Module (VRM)</italic> (<xref ref-type="sec" rid="sec4_3">Section 4.3</xref>). Buck converter with synchronous rectifier, with MOSFET1 being the main switch, and MOSFET2 the synchronous rectifier. See <xref ref-type="fig" rid="fig36">Fig. 36</xref> for the time histories of the driving voltages for these MOSFETs. The output voltage <italic>Vo</italic> is across the resistance <italic>R</italic> on the right<italic>Voltage Regulator Module (VRM)</italic> (<xref ref-type="sec" rid="sec4_3">Section 4.3</xref>). Buck converter, driving voltage wave-form excitations <italic>V<sub>SW</sub></italic> (Synchronous Switch, MOSFET1, <xref ref-type="fig" rid="fig35">Fig. 35</xref>) and <italic>V<sub>SR</sub></italic> (Synchronous Rectifier, MOSFET2, <xref ref-type="fig" rid="fig35">Fig. 35</xref>)
Using the GFM transformed model of the VRM^{27}
The same results were obtained using the Falk method in its original form.
, the transient simulation time was reduced by 17 times (i.e., the speed-up factor), compared to that of the full two-dimensional thermal model, as documented in Table 6.
<italic>Voltage Regulator Module (VRM)</italic> (<xref ref-type="sec" rid="sec4_3">Section 4.3</xref>). Comparison of simulation time. The original model has 215 unknowns (nodal temperatures). The reduced-order model (ROM) has 51 unknowns (modes). See also <xref ref-type="table" rid="tab1_5">Table 5</xref> for the full-bridge converter
Simulation model
SABER transient analysis
Speed-up factor
Full model
6690 s
GFM/WYD
391 s
6690/391 = 17.11
Remark 4.5 Most published results in the VLSI-CAD literature to demonstrate simulation speed-up have been in the frequency domain, essentially using methods of the Lanczos type based on the moment-matching technique. Many widely cited papers (e.g., [29, 68]) showed speed-up ratios essentially restricted to the linear part of the circuit (interconnects). On the other hand, the way we computed the speed-up ratio is different: We do not restrict the comparison to the reduced field problem (heat sink) alone, but include the complete field-circuit network in the calculation. If we restricted our calculation of the speed-up ratio to the heat sink alone, we would obtain a much higher speed-up ratio (with orders of magnitude in improvement). The speed-up ratio in our circuit examples, which involved highly nonlinear components (IGBT, MOSFET), was based on the reduction of total simulation time, and not on just the reduction of simulation time for the heat sink (field component), as mentioned above.
Comparison of CPU Times for Different Methods
In this section, we provide a comparison of the CPU times spent by the full-order models and by two reduced-order models implemented in the circuit simulator SABER.
For a fair comparison between the full-order model and the 1-D transformed models, we need to clock the transformation time in MATLAB. Table 7 shows the MATLAB CPU time spent to generate the transformed 1-D thermal network used in the above simulations. Once the 1-D model had been generated, it can be reused for other design/redesign applications. Such reuse capability is especially advantageous for large complex coupled field-circuit systems.
<italic>Generation of 1-D transformed models</italic> (<xref ref-type="sec" rid="sec4_4">Section 4.4</xref>). Mesh 1 has 16 nodes. Mesh 2 has 215 nodes. Comparison of CPU times (s)
Simulation model
CPU times (s)
GFM
Original falk
Mesh 1: Coordinate transformation
0.03333
0.03333
Mesh 2: Coordinate transformation
141.5
117.3
Mesh 2: Coordinate transformation and reduction
15.58
14.22
The transformation times for a model of the heat sink (field component) with 1,226 dofs using four different Krylov-based transformation methods considered here, i.e., GFM, original Falk method, WYD method, and Lanczos method, are given in Table 8. The results clearly show that the original Falk method and the GFM are more advantageous than the WYD method and the Lanczos method. The reason is because there is no inversion of the system conductance matrix K in the Falk method and the GFM, in contrast to the WYD method and the Lanczos method.
<italic>Generation of 1-D transformed models</italic> (<xref ref-type="sec" rid="sec4_4">Section 4.4</xref>). Mesh 3 has 1,226 nodes (dofs). Comparison of CPU times (s)
CPU times (s)
Simulation model
GFM
Falk
WYD
Lanczos
Coordinate transformation
613.33 s
609.71 s
1,854.1 s
1,834.4 s
The results in Fig. 37, obtained from FE Mesh 2 for the buck-converter VRM, show that there was a significant increase in the bulk-drain total energy dissipation of the MOSFET with simulation time, for both the full-order model and the GFM model.
The Si-junction temperature rise and the converter regulated output voltage vs. simulation time are given in Fig. 38. Again, the result without thermal effects is obtained from simulations where the Si-junction temperature is fixed at the ambient temperature. The regulated voltage difference increases with time. Such increase also indicates power dissipation due to thermal effects.
The results show that accounting for the field-circuit coupling is essential for accurate simulations. By applying our proposed model reduction strategies based on the GFM, the speed-up ratio can go up to about 9 for the full bridge converter example (Table 5) and 17 for the VRM example (Table 6), which is highly competitive compared to the traditional Lanczos method.
<italic>VRM:</italic> Total energy dissipation in MOSFET <italic>vs</italic>. time, using GFM (Mesh 2)<italic>VRM: Si</italic>-junction temperature and output voltage <italic>vs</italic>. time, using GFM (Mesh 2)Closure
To study truncation effects in the generalized Falk method (GFM) of coordinate transformation, we introduce the concept of response participation factors α_{i} for the selection of trial vectors for model reduction, and show that these response participation factors are more effective than the load participation factors β_{i} traditionally employed. We also compared the simulation times to generate the trial vectors and the error norm committed on the load representation for various coordinate-transformation methods. Unlike the Lanczos method, the GFM does not produce unstable positive poles in circuit simulation, and is more efficient compared to the Lanczos method and the Arnoldi method in the case where a lumped capacitance matrix is used in electrothermal simulation, a case study to illustrate the effectiveness of the proposed methodology for general coupled field-circuit problems, in particular IC interconnects.
At the end of the combined two-stage model-reduction strategy, we obtain a reduced-order model with an identity capacitance matrix and a tridiagonal conductance matrix represented by a simple 1-D equivalent circuit. For the field components in overall coupled field-circuit systems, we develop simple equivalent field circuit networks, which are then implemented in circuit simulators.
Numerical examples of coupled field-circuit problems involving a full-bridge converter with four IGBTs as switching devices and a Voltage Regulator Module with two MOSFETs showed a remarkable efficiency of the proposed coordinate-transformation methods and model-reduction strategies.
Funding Statement: The authors thank the research support of the National Science Foundation.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
ReferencesLombardi, L., Tao, Y., Nouri, B., Ferranti, F., Antonini, G. et al. (2019). Parameterized model order reduction of delayed PEEC circuits. Hiruma, S., Igarashi, H. (2020). Model order reduction for linear time-invariant system with symmetric positive-definite matrices: Synthesis of cauer-equivalent circuit. Hasanpour, S., Baghramian, A., Mojallali, H. (2019). Reduced-order small signal modelling of high-order high step-up converters with clamp circuit and voltage multiplier cell. Cangellaris, A. C., de Zutter, D. (2004). Special issue on model-order reduction methods for computer-aided design of RF/microwave and mixed-signal integrated circuits and systems-guest editorial. Dong, X., Griffo, A., Wang, J. (2020). Multiparameter model order reduction for thermal modeling of power electronics. Acle, Y. G. I., Freitas, F. D., Martins, N., Rommes, J. (2019). Parameter preserving model order reduction of large sparse small-signal electromechanical stability power system models. Wang, Z., Li, Y., Li, Z., Zhao, C., Gao, F. et al. (2019). Reduced-order DC terminal dynamic model for multi-port AC-DC power electronic transformer. Nicolini, J. L., Na, D. Y., Teixeira, F. L. (2019). Model order reduction of electromagnetic particle-in-cell kinetic plasma simulations via proper orthogonal decomposition. Zimmerling, J., Druskin, V., Zaslavsky, M., Remis, R. F. (2018). Model-order reduction of electromagnetic fields in open domains. Zhai, Y., Vu-Quoc, L. (2007). Analysis of power magnetic components with nonlinear static hysteresis: Proper orthogonal decomposition and model reduction. Cordesse, P., Remigi, A., Duret, B., Murrone, A., Ménard, T. et al. (2020). Validation strategy of reduced-order two-fluid flow models based on a hierarchy of direct numerical simulations. Alfaro-Isac, C., Izquierdo-Estallo, S., Sierra-Pallares, J. (2020). Reduced-order modelling of equations of state using tensor decomposition for robust, accurate and efficient property calculation in high-pressure fluid flow simulations. Wang, Y., Ma, H., Cai, W., Zhang, H., Cheng, J. et al. (2020). A POD-Galerkin reduced-order model for two-dimensional Rayleigh–Bénard convection with viscoelastic fluid. Wang, J., Fang, M., Li, H. (2020). An adaptive substructure-based model order reduction method for nonlinear seismic analysis in opensees. Wang, D., Zlotnik, S., Dez, P. (2020). A numerical study on hydraulic fracturing problems via the proper generalized decomposition method. Mohan, A., Gaitonde, D. (2018). A deep learning based approach to reduced order modeling for turbulent flow control using LSTM neural networks. https://arxiv.org/abs/1804.09269.Lee, M. W., Dowell, E. H. (2020). Improving the predictable accuracy of fluid Galerkin reduced-order models using two POD bases. Sugar-Gabor, O. (2020). Parameterised non-intrusive reduced-order model for general unsteady flow problems using artificial neural networks. Zhong, H., Xiong, Q., Yin, L., Zhang, J., Zhu, Y. et al. (2020). CFD-based reduced-order modeling of fluidized-bed biomass fast pyrolysis using artificial neural network. Das, A., Khoury, A., Divo, E., Huayamave, V., Ceballos, A. et al. (2020). Real-time thermomechanical modeling of PV cell fabrication via a pod-trained RBF interpolation network. Vu-Quoc, L., Humer, A. (2022). Deep learning applied to computational mechanics: A comprehensive review, state of the art, and the classics. Hsu, J. T., Vu-Quoc, L. (1996). A rational formulation of thermal circuit models for electrothermal simulation-part II: Model reduction techniques. Lanczos, C. (1950). An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. Wilson, E. L., Yuan, M. W., Dickens, J. M. (1982). Dynamic analysis by direct superposition of Ritz vectors. Remis, R. F. (2004). Low-frequency model-order reduction of electromagnetic fields without matrix factorization. Wu, H., Cangellaris, A. C. (2004). Model-order reduction of finite-element approximations of passive electromagnetic devices including lumped electrical-circuit models. Xu, J. (1989). An improved method for partial eigensolution of large structures. Bertolini, A. F. (1998). Review of Eigen solution procedures for linear dynamic finite element analysis. Pillage, L. T., Rohrer, R. A. (1990). Asymptotic waveform evaluation for timing analysis. Verghese, N. K., Lee, S. S., Allstot, D. J. (1993). A unified approach to simulating thermal and electrical substrate coupling interaction in IC’s. Proceedings of 1993 International Conference on Computer Aided Design, pp. 422–426. IEEE.Vu-Quoc, L., Zhai, Y. H., Ngo, K. D. T. (2004). Efficient simulation of coupled circuit-field problems: Generalized Falk method. Floros, G., Evmorfopoulos, N., Stamoulis, G. (2020). Frequency-limited reduction of regular and singular circuit models via extended Krylov subspace method. Li, M., Chen, R. (2016). Method of moments wide-band simulations of the microstrip circuits with second-order Arnoldi reduced model. Zeng, X. A., Yang, F., Su, Y. F., Cai, W. (2010). NHAR: A non-homogeneous Arnoldi method for fast simulation of RCL circuits with a large number of ports. Li, Z., Shi, C. J. R. (2006). A quasi-Newton preconditioned Newton-Krylov method for robust and efficient time-domain integrated circuit simulation. Silveira, L., Kamon, M., Elfadel, I., White, J. (1999). A coordinate-transformed Arnoldi algorithm for generating guaranteed stable reduced-order models of RLC circuits. Odabasioglu, A., Celik, M., Pileggi, L. T. (1998). PRIMA: Passive reduced-order interconnect macromodeling algorithm. Kerns, K. J., Yang, A. T. (1998). Preservation of passivity during RLC network reduction via split congruence transformations. Bai, Z., Feldmann, P., Freund, R. (1997). Stable and passive reduced-order models based on partial Padé approximation via the Lanczos process. Technical Report 97/3-10, Bell Laboratories, Lucent Technologies, 700 Mountain Avenue, Murray Hill, NJ 07974-0636..Feldmann, P., Freund, R. W. (1995). Efficient linear circuit analysis by Padé approximation via the Lanczos process. Zimmerling, J., Wei, L., Urbach, P., Remis, R. (2016). A Lanczos model-order reduction technique to efficiently simulate electromagnetic wave propagation in dispersive media. Bracken, J. E., Sun, D., Cendes, Z. (1999). Characterization of electromagnetic devices via reduced-order models. Cangellaris, A. C., Zhao, L. (1999). Passive reduced-order modeling of electromagnetic systems. Liu, W., Ramm, R. J., Tepe, A. D., Campbell, C. K., Sasaridis, D. (2018). US patent application No. US20180114740A1.Kroeze, R., Campbell, C., Kalayjian, N. R. (2014). Fast switching for power inverter. US patent No. US8760898B2.Un-Noor, F., Padmanaban, S., Mihet-Popa, L., Mollah, M. N., Hossain, E. (2017). A comprehensive study of key electric vehicle (EV) components, technologies, challenges, impacts, and future direction of development. Reimers, J., Dorn-Gomba, L., Mak, C., Emadi, A. (2019). Automotive traction inverters: Current status and future trends. Erickson, R. W., Maksimovic, D. (2020). Toshiba Electronic Devices & Storage Corporation, Full Bridge Converter-Basics of Switching Power Supplies-Video Gallery. Original website.Ahmed, M. H., Lee, F. C., Li, Q. (2020). Two-stage 48 V VRM with intermediate bus voltage optimization for data centers. Ammous, A., Ghedira, S., Allard, B., Morel, H., Renault, D. (1999). Choosing a thermal model for electrothermal simulation of power semiconductor devices. Ammous, A., Ammous, K., Morel, H., Allard, B., Bergogne, D. et al. (2000). Electrothermal modeling of IGBT’s: Application to short-circuit conditions. Hughes, T. J. R. (1987). Hsu, J. T., Vu-Quoc, L. (1996). A rational formulation of thermal circuit models for electrothermal simulation-Part I: Finite element method. Falk, V. S. (1955). Die abbildung eines allgemeinen Schwingungssystems auf eine einfache Schwingerkette. Paige, C. C. (1972). Computational variants of the Lanczos method for the eigenproblem. Simon, H. D. (1984). The Lanczos algorithm with partial reorthogonalization. Freund, R. W. (1999). Reduced-order modeling techniques based on Krylov subspaces and their use in circuit simulation. Sridevi, N., Venkatramanan, D., Adapa, A. K., John, V. (2016). Analogy, Inc (1997). Hefner, A. R., Diebolt, D. M. (1991). An experimentally verified IGBT model implemented in the Saber circuit simulator. 22nd Annual IEEE Power Electronics Specialists Conference, pp. 10–19. MIT, Cambridge, MA, IEEE. DOI 10.1109/PESC.1991.162645.Hefner, A. R. (1992). Semiconductor measurement technology: INSTANT-IGBT network simulation and transient analysis tool. In: Hefner, A., Blackburn, D. (1993). Simulating the dynamic electrothermal behavior of power electronic circuits and systems. Hefner, A. R., Diebolt, D. M. (1994). An experimentally verified IGBT model implemented in the saber circuit simulator. Zhou, X. W., Wong, P. L., Xu, P., Lee, F. C., Huang, A. Q. (2000). Investigation of candidate VRM topologies for future microprocessors. Panov, Y., Jovanovic, M. M. (2000). Design considerations for 12-V/1.5-V, 50-A voltage regulator modules. IEEE Applied Power Electronics Conference and Exposition, vol. 1, 19000210. Piscataway, NJ, USA, New Orleans, LA, USA: IEEE.Shenai, K., Cavallaro, C., Musumeci, S., Pagano, R., Raciti, A. (2003). Modeling low-voltage power MOSFETs as synchronous rectifiers in buck converter applications. 38th Annual Meeting of the Industry-Applications-Society. IEEE, vol. 1–3, pp. 1794–1801. Salt Lake City, UT.Feldmann, P., Freund, R. (1995). Reduced order modeling of large linear subcircuits via a block Lanczos algorithm. Proceedings of the 32nd annual ACM/IEEE Design Automation Conference, Piscataway, NJ, USA, IEEE.Appendix A. Lanczos method
Lanczos method The Lanczos method [23] was originally proposed as a technique for tridiagonalizing matrices. A sequence of trial vectors are formed by repeatedly multiplying the matrix to be reduced with a starting vector. These trial vectors form a Krylov sequence. Each new trial vector is orthogonalized with respect to the two previous trial vectors. This orthogonalization procedure can be shown to be sufficient to obtain orthogonality with all previously calculated vectors. The coefficients computed from the orthogonalization process are then combined to form a tridiagonal matrix that theoretically has the same eigenvalues as the original matrix, after all n (Lanczos) vectors has been calculated, with n being the order of the system.
(7) Approximated eigenpairs of the original system:
X=ZrY∈Rn×r,ωi=1λi,i=1,…,r,
with Zr=[z1,…,zr]∈Rn×r and Y∈Rr×r.
Appendix B. Wilson–Yuan–Dickens (WYD) method
The WYD method, first introduced in [24], has become popular in structural dynamics, and is based on the direct superposition of a special class of Ritz vectors generated from the spatial distribution f(s)∈Rn×1 of the dynamic load that can be put under the form f(s)g(t), where g(t) is function of time. This method eliminates the requirement for exact evaluation of the free vibration frequencies and mode shapes. The WYD approach does not provide a full uncoupling of the equilibrium equations but has proved to be more efficient than the traditional eigensolution approach. The sequence of WYD vectors is generated by taking into account the spatial distribution of the external excitation, which is important information neglected by direct use of exact mode shapes. The first WYD vector is the steady-state vector obtained from a steady-state analysis using the spatial distribution f of the dynamic excitation as input. The other vectors are generated from a recurrence relationship in which the capacitance matrix C is multiplied by the last WYD vector; the resulting vector is then used as the excitation for the next steady-state solution.
WYD Algorithm
(1) Given C, K∈Rn×n, and f∈Rn×1
(2) Triangularize K = LDL^{T}
(3) Solve for load-dependent starting vector w1∗, and C-normalize w1∗ to obtain w_{1}
Kw1∗=f⇒w1∗=K−1fw1=w1∗w1∗TCw1∗
(4) Solve for additional vectors i = 2,… , r, with r ≤ n:
Kwi∗=Cwi−1⇒wi∗=K−1Cwi−1ci,j=wjTCwi∗wi∗∗=wi∗−∑j=1i−1ci,jwj(C-orthogonalization)bi=wi∗∗TCwi∗∗wi=wi∗∗bi(C-normalization to obtain wi)
In the reduced basis of WYD trial vectors, the capacitance matrix Cr∗:=WrTCWr=Ir∈Rr×r is an indentity matrix, whereas the stiffness matrix Kr∗:=WrTKWr∈Rr×r is full in general.