Fast simulation of transient temperature distributions in power modules using multi-parameter model reduction

: In this study, a three-dimensional model with multi-parameter order reduction is applied to the thermal modelling of power electronics modules with complex geometries. Finite element or finite difference method can be used to establish accurate mathematical models for thermal analyses. Unfortunately, the resulting computational complexity hinders the analysis in parametric studies. This study proposes a parametric order reduction technique that can significantly increase simulation efficiency without significant penalty in the prediction accuracy. The method, based on the block Arnoldi method, is illustrated with reference to a multi-chip SiC power module mounted on a forced air-cooled finned heat sink with a variable mass flow rate.


Introduction
The design of modern power electronic systems is increasingly complex, requiring multi-domain optimisation encompassing the interactions among the electrical, thermal and mechanical domains. Nevertheless, due to the increasing demands for higher power density, the thermal analysis and management of power electronics systems is becoming more and more important. Reliability of a power electronics converter is significantly affected by operating temperature and temperature cycling. It is well known that components' lifetime decreases exponentially with temperature and that thermal and thermo-mechanical failure modes in devices and packaging are accelerated by temperature cycling.
The thermal management of power converters is becoming increasingly demanding because of its strong industrial drive for smaller, more efficient power electronics systems [1,2]. Heat exchangers account for a significant portion of power converters' mass. Increased converter efficiency and better thermal design can contribute to significant reduction of heat exchangers' mass and therefore increased power densities. Therefore, accurate modelling tools for thermal analyses can significantly aid the design optimisation of power converters, helping the design engineer to select the optimal system design with the required heat dissipation. However, the extent to which the system design can be optimised for size and weight is limited by the maximum rated components temperature, which cannot be exceeded during normal operation [3].
Additionally, the life expectancy of the power converter is reduced due to thermal cycling, where the damage from multiple heating and cooling events accumulates until the system failure occurs [4]. As a result, temperature monitoring systems to protect products and meet regulatory requirements are increasingly being adopted in safety critical applications. Compact thermal models can be used as an aid for the estimation and monitoring of the temperature of components during real-time operation.
A vast literature on models for the thermal analysis of power electronic systems has been published. The simplest methods use compact thermal models typically based on empirically derived lumped element models such as those based on Foster or Cauer networks. Although computationally efficient, these lumped parameter models typically require experimental calibration and cannot be easily employed in parametric studies where geometries or operating conditions change. More accurate and physically representative methods for thermal modelling of power assemblies including devices, packaging and heat exchangers rely on a number of well-established numerical tools that discretise the distributed partial differential equations (PDEs) that model the heat-transfer problem using finite element method (FEM), finite difference method (FDM) or compact thermal model based on an analytical or empirical lumped parameter model (LPM). In addition, accurate numerical modelling of heat exchangers with natural or forced convection typically requires the use of computational fluid dynamics (CFD) methods. Typically, CFD software tools can simultaneously solve conductive and convective heat transfer problems, providing the most accurate and detailed temperature distribution for power electronic systems. Unfortunately, CFD analyses are extremely demanding in terms of computing resources and calculation time [5]. A number of model order reduction (MOR) techniques have been proposed to alleviate the problems of computational complexity arising from the simulation of complex and distributed dynamical systems. MOR techniques applied to thermal problems use the discretised version of the underlying PDEs generated using either FEM or FDM to produce a reducedorder model that significantly reduces computational complexity while guaranteeing reasonably accurate results [6]. In this paper, an FDM with MOR is selected to establish a mathematical model.
A number of MOR techniques have been proposed for application in the thermal modelling problem. Among the most effective strategies, Guyan reduction [7] and Krylov subspace methods have been proposed. The thermal performance of a power converter depends not only on the layout of components but also on boundary conditions such as the coolant mass flow rate. It is therefore important that the compact models used in thermal analyses conserve the dependency on these design parameters and operating conditions. Unfortunately, once MOR techniques are applied to the original model formulation, the dependency on parameters, e.g. the coolant mass flow rate, disappears. This results in the need to repeat the MOR process for every different operating condition, making parametric studies of system operation in different operating conditions (e.g. different ambient temperature or coolant mass flow rate) extremely tedious.
The paper presents a parametric MOR method that conserves one or more parameters in the reduced-order model, making analyses of the converter in different operating conditions computationally efficient. The method, based on multi-moment matching and block Arnoldi's orthogonalisation on standard Krylov subspaces, is analytically derived. The method is illustrated, and its benefits are demonstrated with reference to a power module mounted on a forced air-cooled finned heat sink. Detailed comparisons with commercial CDF software demonstrate the accuracy and computational efficiency of the proposed method.

Parametric model order reduction
A geometry-based mathematical model is needed to construct the thermal model of the power module and its cooling assembly. A geometry-based method has the advantage that can be used as a tool in the module design process by facilitating the optimisation of components' placement, distances etc., since the topology is directly taken into account.

Conventional model order reduction
In each thermal simulation, the temperature distribution is computed on a discrete grid, and its size can produce millions of ordinary differential equations, depending on the complexity of the geometry [7]. If a complete model is used, system-level simulations quickly become unmanageable. The three-dimensional PDE describing the conductive heat transfer problem can be discretised into a system of ordinary differential equations (ODEs) as where C is the thermal specific heat matrix, K is the thermal conductivity matrix, Q is the heat generation vector and T is the vector of temperatures in all the n points of the discretised domain. F ∈ R n×m and E ∈ R n×p are the input and the output matrices, and m and p denote the number of inputs and outputs, respectively [7][8][9]. As a result, transforming (1) into the frequency domain result in a matrix-valued rational transfer function G: C C p×m given by Arnoldi-based reduction is a well-established MOR tool [7], whose goal is to transform the equation system (1) into a system of lower dimensionality but in the same form [7]: where z ∈ R r is obtained by projecting the original state T of dimension n to a sub-space of dimension r ≪ n verifying The transformation is obtained by a projection process based on the Padé-type approximation where the reduced-order system matrices are obtained as follows [8,9]: and V is an output of the Arnoldi algorithm. Before the block Arnoldi can be employed, the two matrices C and K have to be reduced to a single matrix, denoted by A in the following. This can be done by rewriting (2) as follows: subspace after building block Krylov subspaces. The matrix V is composed from r -dimensional vectors that form a basis for the right Krylov subspace of the dimension r: After building the block Krylov subspaces, Arnoldi's orthogonalisation, which is shown in Table 1, is carried to extend the classical Arnoldi algorithm to block Krylov subspaces.

Multi-parameter model order reduction
In this section, a parameter-independent MOR method is proposed based on multi-series expansion with respect to a set of heat transfer coefficients. As in the non-parametric case, ODEs of the form (1) and (2) are considered. In this case, the convective boundary layer is assumed to have a multi-parameter dependency on air mass flow rate. In an air-cooled system, the temperature variation in the mass of air can be neglected compared with the solid part of the power module substrate and heat sink assembly and the interaction with the cooling medium described by a convective boundary layer. Consequently, the multi-parameter condition only happens in the matrix of conductance. Then, the ODEs can be rewritten as where p i represent the parameters that are required to be kept in the reduced model. The projection matrixV can be used to calculate the reduced-order temperature vector z whose dynamics are described as Similar to the conventional MOR, the transfer function of the system in (5) is formed as which can be written as Many methods for multi-parametric order reduction have been proposed. There are two main strategies based on MOR with or without moment matching, such as [10][11][12][13][14][15][16] or reduction without multi-moment matching [17,18]. In this paper, reduction with multi-moment matching is introduced. The process is based on the Taylor-series expansion of the transfer function H(s) around a certain point s 0 .The moments of the transfer function (10) are the coefficients of its Taylor series expansion. It is worth noting that only when there is a weak correlation between the parameters, the mixing moment can be ignored without affecting the precision [19]. This is typically the case for thermal problems [20,21]. For the problem under investigation, i.e. the thermal analysis of power modules with cooling system, the parameters series p 1 p 2 ... p i and submatrix K 1 K 2 ... K i only appear in the equations of the states on boundary layer of the baseplate. The block Arnoldi's orthogonalisation based on standard Krylov subspaces for multimoment matching needs to be applied [7][8][9]. The next step is to make another expansion but this time in series of each parameter p i , for each moment. For the first moment For the second moment For the jth moment Order reduction with moment matching needs moment m 0 to m j be independent on parameter series p 1 p 2 ... p i . Equation (13) shows that the moments are combination of the matrices which means that each moment lies in the subspace spanned by the columns of the matrices in (10). These matrices are then taken to construct the projection matrix V, which is located on the first r columns of (14). This can be done by rewriting matrices A and B in Section 2 as m columns of the matrix B = [ B 1 B 2 … B m ] are the starting vectors of the so-called block Krylov-subspace. The following calculation is the same with Section 1. The details of this algorithm are shown in Table 2.

Thermal model and simulation results
A simplified power module mounted on a parallel-plate finned heat sink is considered, as shown in Fig. 1. The power module contains six SiC MOSFETs and is assumed to be mounted on a direct copper bonded ceramic substrate and attached to the heat sinks. The boundary condition is complex due to the forced air cooling. A simplified analytical model of the heat transfer coefficient varying axially along the direction of the air flow has been established based on [22,23]. In this paper, a dimensionless fluid dynamic entry length is introduced is the heat sink channel aspect ratio and ϵ = fin thickness/channel space .
An analytical model for the Nusselt number Nu A in [23] is suitable for the heat sink model, as follows: for j = 1, ... , j max − 1 One matrix multiplication where m is the model blending parameter provided in [22] and other parameters of (17) are given in Table 3. The Nusselt number decreases along the thermal entry length [23] and settles to a constant value. It should be noted that 0.1 < Pr < ∞ is valid for most heat exchanger applications. z* is the dimensionless thermal axial position. The friction factor Reynolds product equation (18) and (19) describes the effect of the boundary layer velocity profile on the mass transfer [24]: f Re A, f d = 12 With this and with Nusselt number Nu A , the heat transfer coefficient becomes The resulting heat transfer coefficient as a function of the axial distance from the inlet for the heat sink in Fig. 1 for three different values of air mass flow is shown in Fig. 2.
The proposed MOR is applied to the system. Figs. 3a-c illustrate surface temperature responses of the three MOSFETs identified in Fig. 1 calculated by reduced-and full-order simulation and compared to results obtained with commercial FE software ANSYS mechanical and the CFD tool ICEPAK. The 'ANSYS mechanical' simulations are obtained by solving the conduction heat transfer problem, removing the fins and setting the bottom baseplate surface with a convective boundary condition with heat transfer coefficient as in Fig. 2. As can be seen, the agreement among the finite difference full order, the proposed reduced-order method and the ANSYS FE method is excellent. However, some discrepancies are present when compared with the CFD results. This is due to the approximations resulting from the semianalytical model of the variable heat transfer coefficient. The discretisation employed in the full-order model results in a system with 5711 nodes, while the reduced order has 108 states corresponding to 18 temperature nodes per MOSFET. On the same computer and with the same mesh size, CFD takes over 300 min, the full-order simulation needs 20 min, while the reduced-order simulation only takes about 5 s.

Conclusion
In this paper, a novel multi-parameter order reduction is developed and applied to a power module with forced air-cooled systems. The multi-moment matching technique is used to preserve in the reduced order a number of parameters, making calculations in variable operating conditions significantly more efficient. An example of a power module cooling system with different mass air flow rates is reported.
A high degree of accuracy compared to that of conventional FE and CFD tools is shown. A significant increase in computational efficiency is demonstrated resulting in faster calculation time and memory requirements.

Acknowledgment
This work was supported by the European Commission Horizon 2020 -Mobility for Growth Program under Grant 636170.