Thermal Resistance Matrix Extraction from Finite-Element Analysis for High-Frequency Magnetic Components

The thermal management of magnetic components for power electronics is crucial to ensure their reliability. However, conventional thermal models for magnetic components are known to have either poor accuracy or excessive complexity. Contrary to these models, the use of Thermal Resistance Matrices is proposed in this paper instead, which combine both accuracy and simplicity. They are usually used to characterize semiconductor devices, but not for magnetic components. The guidelines to apply Thermal Resistance Matrices for magnetic components are discussed in detail. The accuracy of this model is validated by 3D FEA simulations and experimental results, showing an absolute error lower than 5 ∘C and a relative error between −6.4% and 3.9%, which is outstanding compared to the carried-out literature review.


Introduction
The constantly increasing interest in miniaturising power components makes crucial a proper characterisation of the thermal behaviour of magnetic components, which are one of the limiting factors [1,2].
The simplest and most widely used approach is the empirical Equation (1) from [3], or other simple equations like the one proposed in [4]. Nevertheless, their accuracy is considerably limited, while it serves as a first order approximation.
Another recurrent approach is to use thermal networks for magnetic components [5][6][7] (or for semiconductors [8]), which shows better results than Equation (1). Their accuracy is associated with their granularity, hence to their complexity, which is still limited since a constant film coefficient is commonly assumed. To compensate this fact, some authors model the fluid dynamics by means of non-linear resistances [9][10][11][12], which add extra complexity to the model. Different approaches consist of building a model based on experimental results, explained in [13,14], or using an inverse model in which the correction of some analytical equations is made by means of prototype results [15]. However, they both have the inherent drawback of requiring a prototype and their application is limited to the boundary conditions of the experiments [16].
In this article, the use of thermal resistance matrices for magnetic components is proposed. Some applications of this model to semiconductors can be found in [17][18][19]. This approach offers several advantages compared to the previous ones. This model consists of a simple coefficients matrix of dimension NxN, where N is the number of characterised objects or parts (for example, the core, the windings and the bobbin). This matrix characterizes the temperature at each part, accounting for the self heating effects as well as the mutual heating between each part. Additionally, regardless of the heterogeneity of the device or the complexity of the boundary conditions, the thermal resistance matrix captures the thermal behaviour of the device very accurately.
As an example, the temperature rise of an E25/13/7 transformer predicted by 3D Finite Element Analysis (FEA) simulations for a range of winding power losses with constant core losses. This is illustrated in Figure 1 along with the results obtained by using Equation (1) and the proposed thermal resistance matrix. It can be seen that the last one shows considerably better results.

Proposed Model
The accuracy of 3D FEA simulations, as well as their capability to analyze complex systems, allows the obtaining of a static FEA-based thermal model to estimate the temperature of any object within a system under any operating point [20]. Furthermore, any boundary condition can be recreated by simulation, regardless of their nature or complexity [21,22].
While the concept of a thermal resistance matrix is generic, some special concerns need to be considered when applying it to magnetic components, which are discussed next.

Modeling Transformers and Inductors
A step-by-step guide to obtain a thermal resistance matrix to characterize any magnetic component is explained along this section. A flow diagram of the corresponding procedure steps is depicted in Figure 2, as a summary of the proposed methodology.

Representative Objects in a Magnetic Component
The first step consists of splitting the magnetic component in the least amount of representative objects. In other words, since a thermal resistance matrix of dimension NxN characterizes the interaction between N objects, there is an accuracy-complexity compromise in the model. The main 'objects' to be characterised are described next: • Magnetic core: its temperature rise is crucial to analyze its magnetic behaviour (permeability and saturation flux).
• Main windings: windings with considerably different physical characteristics (geometry, materials, power losses) can be treated as separate objects for higher accuracy or together to reduce N for simplicity. • Auxiliary windings: whereas their power losses are negligible, they can be treated as an additional object if their temperature is relevant for the analysis. • Passive parts: a coil former or the pin connectors to a PCB do not have power losses, but their temperature depends on the heat flux coming from surrounding objects, which could be critical if the maximum ratings are exceeded.

Linearization
In order to build a thermal resistance matrix, the superposition principle must be applied, which is only possible in linear systems. However, only the thermal conduction is linear while convection and radiation are non-linear [23,24]. This effect is shown in Figure 1, where the temperature rise of a certain magnetic core, exposed to natural convection and radiation, is depicted for different power losses. Contrary to semiconductors, where the conductive heat transfer is dominant due to their reduced size and encapsulation, the evolution of this temperature rise is clearly non-linear because the magnetic components are more exposed to convection and radiation. Therefore, the system must be linearized and a proper thermal operating point must be selected since the slope of the linearized 'resistance' will change accordingly. To achieve realistic results, the thermal operating point must be equal (or close) to the limit temperature rise of the corresponding object, ∆T s lim (see Figure 1). Each object can have a different ∆T s lim , depending on the materials and the specifications of the project.
By chosing this ∆T s lim for each object, a value of power losses per object will be calculated so that a thermal resistance per object can be obtained to characterise its temperature rise from T amb to ∆T s lim -assuming a linear behaviour-, as explained in the following steps.
As a summary, this steps consists only in choosing the maximum allowed (or limit) temperature rise, ∆T s lim , for each object.

P test Calculation for the Magnetic Core, the Windings and the Passive Objects
The thermal resistance is equal to the temperature rise divided by the power losses (or heat) flowing through an object if the system is linear, according to [23,24]. To obtain that thermal resistance, certain test power losses (P test , in Watts) are injected and the subsequent temperature rise is measured. For semiconductors, a common choice consists of making P test equal to the maximum rated power of the device, given by the manufacturers.
However, selecting P test for each part of a magnetic component is not obvious due to the non-linear behaviour explained in the previous step. As the temperature rise of each object must be equal (or close) to its ∆T s lim , the value of P test must be selected accordingly. Two methods can be used to obtain P test : by trial and error in simulation, which is timeconsuming and therefore not recommended, or using basic heat transfer equations and some geometrical approximations, which is explained next.
For the magnetic core, P test is calculated in Equation (2) using the convection and radiation expressions from the external surfaces of the core to the ambient, which are defined in Equations (3) and (4), respectively: where Q conv and Q rad are the heat transferred by convection and radiation (in Watts), respectively; A s is the external heat exchange surface exposed to the ambient, in squared meters; is the relative emissivity of the object; σ is the Stefan-Boltzmann constant (5.670373 × 10 −8 W m 2 ·K 4 ); and h s is the film factor coefficient corresponding to each surface, in W/m 2 ·K. The ambient temperature (T amb ) and T s lim must be expressed in Kelvin.
The typical film factors for natural convection are shown in Table 1. The surfaces in contact with a PCB, which are considered adiabatic, are neglected in this step. However, if the cooling capability of heatsinks and/or forced convection are considered, a simplified film factor must be calculated for each surface according to [23][24][25].
On the other hand, the windings of the magnetic components are treated as a group since only the wires in the outer layer are in direct contact with ambient. Then, the whole winding is treated as a vertical cylinder (Table 1), but the external surface of each wire must also be considered in the total heat exchange area, as highlighted in Figure 3. If the external surface of the winding is assumed to be equal to the window height, as expressed in Equation (5), the results would not be accurate. Therefore, the surface of each wire must be considered, leading to the proper external surface value in squared meters, calculated in Equation (6). Table 1. Film coefficient for natural convection for typical geometries [23,24].

Object Film Coefficient
Vertical plate or cylinder A ext = (2 · r · n wires ) · perim wind (5) where n wires is the number of wires in the external layer of conductors, r is their radius and perim wind is the perimeter of that layer. In case the winding is made of litz bundles, homogenization techniques can be used to simplify them [26,27].
Since only the external layer of conductors is directly exposed to the ambient, all the windings must be treated as a block in order to calculate the required P test , according to Equations (2)-(4). In order to simplify the calculations, it is assumed that the same P test is applied to all the windings: Finally, the power losses in the auxiliary windings are usually negligible compared to the main windings. Therefore, they can be considered as passive objects, for which P test is not calculated. This is also applicable for the bobbin or the pin connectors to the PCB.

Superposition
Once ∆T s lim and P test are calculated for every object, the next step is to identify which simulations are required to characterize the complete system.
The generic transformer from Figure 4 is considered as an example. It is constituted by five objects (N = 5): the magnetic core, a plastic bobbin and three windings (primary, secondary and auxiliary). Therefore, the corresponding thermal resistance matrix (in • C/W) is defined next: [R th max ] is built by applying the superposition principle because the system has been previously linearized. Each column contains the maximum temperature rise of each object, normalized by the corresponding P test . On the other hand, the auxiliary winding and the plastic bobbin have no power losses associated, so their corresponding columns (4th and 5th) are filled with zeros. However, these columns must be included in the matrix in order to calculate the maximum temperature of all the objects in the last step of the design process. The diagonal terms of this matrix (i = j, being i and j the indexes for the rows and the columns, respectively) represent the self-heating behaviour of each object due to the internal heat generation, while the remaining terms (i = j) express the mutual heat flux between each object.
As a conclusion, only one simulation per active object (objects with power losses) is required to obtain the temperature of every object for the calculated P test . In this example, three simulations must be performed, as shown in Figure 4. The calculated P test in the previous step are assigned to a single object in each simulation. The subsequent maximum temperature rise per object from that simulation is then stored to obtain the corresponding column of [R th max ] until the matrix is complete.

Using the Model to Estimate the Maximum Temperature Per Object
The combination of the power losses (in Watts) applied in each object, [P loss ], is defined as a vector in Equation (9). The maximum temperature of each part of the magnetic component (in • C) is calculated in Equation (10) by multiplying the thermal coefficients matrix, [R th max ], by any combination of power losses applied in each object, [P loss ].
[P loss ] = [P core , P pri , P sec , P aux , P bobbin ] T (9) It is important to notice that once [Rth max ] is built, the only input required to calculate the temperature at any region of the magnetic component is [P loss ]. This leads to a versatile model that is valid to calculate the temperature of the component for any [P loss ] input, which means for any electrical operating point.
A generic multi-winding transformer is used as an example in this section, but the model can also be applied to inductors. They are just a particular case with just one winding, so the process to obtain [Rth max ] is identical.

Impact and Limitations
The main advantages and the scope of the developed model can be summarized next: • It can be extended to N objects, whose behaviour will be represented in the thermal coefficients matrix. • The self-heating effects, as well as the mutual heat flux between objects, are modelled by this simple matrix, avoiding the unfeasible complexity of the three-dimensional heat transfer equations for such heterogeneous systems as magnetic components [23,24]. • The maximum temperature can be calculated by the proposed matrix for any combination of power losses of the magnetic component. • It is valid for any boundary conditions as long as they can be reproduced in the required simulations. • The accuracy is ensured for temperatures around T s lim , which will be proven in Section 3. As a consequence, it allows us to properly identify the actual thermal limits of the device under test (DUT). • The accuracy is not limited by (nor dependent on) the location of the maximum temperature per component. This is because the matrix is obtained using the maximum temperature per component, not the temperature at certain coordinates from the geometry.
Despite the wide scope and applicability of the proposed approach, it has some limitations: • Since the behaviour of the temperature rise is assumed linear from T amb to T s lim , the accuracy is very high around those temperatures, but it is lower for temperatures in the middle of that range. • The developed model is only valid for the boundary conditions (convection and radiation) established in the simulations. Therefore, a new thermal coefficients matrix must be calculated for different boundary conditions, but it can be easily obtained by changing the parameters in the simulations. • The accuracy for temperatures higher than T s lim is not guaranteed since the linearization of the thermal coefficients matrix is performed from T amb to T s lim .

Experimental Validation
The proposed FEA-based thermal coefficients matrix is validated along this section. First, the electrical setup and the DUT are described. Then, some aspects concerning the required simulations and the generated models are discussed. Finally, the estimated temperature rise of the tested magnetics calculated by means of this model is compared with the thermal measurements and the 3D FEA simulation results. The used simulator is 'Icepak', from 'ANSYS Electronics Desktop 2019 R2'.

Experimental Setup and Tested Devices
The electrical setup is depicted in Figure 5, in which an AC current is injected into the winding of the prototypes to generate core and winding power losses, causing certain temperature rise as a consequence. First, a sinusoidal voltage waveform is generated at certain frequency by a waveform generator (GW Instek, model AFG-2005) to supply the circuit, and it is amplified by means of an RF amplifier (AR, model 150A 100B). The value of the capacitor C res is selected so that the resonant frequency between itself and the magnetic component under test corresponds to the input signal frequency. This capacitor is placed in series with the inductor, but in parallel with the primary winding of the transformer. This way, the impedance of the resonant tank is low enough to achieve a relatively high input current with a relatively low input voltage in both cases, avoiding the saturation of the RF amplifier due to its limited voltage gain. Since the voltage is imposed by the RF amplifier, a variable resistor (R damp in series for the inductor; R load in the secondary of the transformer) is placed in order to adjust the demanded current by the circuit. This way, the voltage and current amplitudes can be selected as desired. Since the power losses in the auxiliary winding of the transformer are assumed negligible, no load is connected to that winding.  A generic inductor and a generic transformer with three windings (primary, secondary and auxiliary) with different core shapes, wire diameter and number of turns are chosen as representative samples in order to cover the scope of the proposed model. They are shown in Figure 6 along with the required resonant capacitors for the electrical setup. The thermal conductivity assumed for each material is also summarized in Table 4. The emissivity have been assumed equal to 0.8 for all surfaces.  Table 4. Thermal conductivity of each material assumed in this study.

Thermal Conductivity
Ferrite core [28] 4 Regarding the thermal measurements, an OEM-PLUS Series optics fibre thermometer [30] is used to measure temperatures approximately located in the hot-spots, according to the cross-section representation of the prototypes from Figure 7. Then, these temperatures are used as a reference for the maximum temperatures.

Model Extraction from FEA Simulations
At this point, the procedure to obtain [R th max ] for both prototypes, described in Section 2.1, is carried out. The ambient temperature of the simulations is set equal to the experiments (T amb = 26 • C). The limit surface temperature is established as T s lim = 100 • C, resulting in a temperature rise of ∆T s lim = 74 • C, which is used to calculate P test for each object. The values of P test for both prototypes are summarized in Table 5, considering natural convection and radiation, as well as including the PCB and the wooden table below the core of the magnetics. The fluid dynamics of the surrounding air, as well as the radiation effects and the surrounding objects (the PCB and the table, whose material properties are defined in Table 4), are simulated in order to recreate the same conditions in the simulations as in the experimental setup. Therefore, the external faces of the boundary box must be defined as 'openings' in ANSYS, and both convection and radiation must be enabled.
The obtained thermal coefficients matrix (in • C/W) for the inductor (object A) is shown next, being core and winding the objects 1 and 2, respectively.
On the other hand, the thermal coefficients matrix (in • C/W) for the transformer (object B) is also shown next, being core, primary, secondary and auxiliary winding the objects 1, 2, 3, and 4, respectively.

Comparison with Measurements and 3D FEA
In this subsection, the thermal coefficients matrix is used to calculate the temperature rise of the prototypes for the same power losses than in the experimental setup. This temperature estimation is compared to the thermal measurements of the prototypes and 3D FEA simulations.
First of all, the Steinmetz equation is used to estimate the core losses [31,32]. Then, the AC resistances are obtained by means of the ANSYS PEmag software [33] to calculate the winding losses [34]. As a result, the calculated losses for the inductor and the transformer experiments are summarized in Table 6.  (11) and (12) for the inductor and the transformer, respectively.
Finally, the comparison between the temperature measurements and the temperature estimation is shown in Figure 8, in which the maximum temperature rise at each region of the inductor and the transformer is depicted. Each region of the magnetics has three different temperatures associated with the bar graphs: the purple one represents the experimental measurements, the orange one is the temperature obtained by 3D FEA simulations and the green one is the temperature estimated by using the proposed [R th max ] matrix. The error margin is depicted in the experimental results (the purple bars) due to the inherent error associated to the measurements (consequence of the resolution of the thermocouples and the deviation of the location of the hotspot). The deviation between the [R th max ] model and the 3D FEA simulations is defined as a percentage over the green bars, which is calculated in Equation (13): This error between [R th max ] and the 3D simulations is a direct measure of the inherent error of the model, avoiding the uncertainties related to measurements. Experimental results are shown as proof and validation.
In addition, thermal images (with a FLUKE-Ti400 thermal camera [35]) are performed to check that the surface temperatures are smaller than the maximum temperatures measured with the thermocouples [30]. The corresponding images are illustrated in Figure 9. The 3D FEA thermal simulations corresponding to the prototype measurements for the inductor and the transformer are shown in Figures 10 and 11, respectively, establishing the same operating conditions than the experiments. The temperature scales are set equal than the thermal camera measurements (Figure 9) and the similarity between the simulations and the experiments can be observed.  It is important to highlight that the model is developed assuming linearization between T amb and T s lim . As a result, the accuracy is higher for temperatures close to these values, but it is compromised in the middle of that range. In these particular experiments, ∆T s lim = 74 • C is used to calculate the P test value for each object. Even though the maximum temperature rise illustrated in Figure 8 ranges between 35 • C and 60 • C, the absolute error between the proposed model and the experimental results is lower than 5 • C, and the relative error between the proposed model and the 3D FEA simulations is between −6.4% and 3.9%. These deviations are quite lower than what other authors consider acceptable [20,36,37].

Conclusions
A methodology to extract thermal resistance matrices from 3D Finite Element Analysis thermal simulations is developed in this paper, reducing the computational requirements while ensuring high accuracy. It can be generalized to any power electronics device, but the corresponding particularities to model inductors and transformers are explained in detail. The main characteristics of this FEA-based thermal model are the accuracy for any combination of power losses in the magnetic component and the simplicity of the thermal coefficients matrix. As a conclusion, the maximum temperature can be predicted by this model at any region of a magnetic component, while some approaches from the literature consider the component just as a single node.
The proposed model is proven valid independently of the core shape, size and winding arrangement by testing a generic inductor and a three-winding transformer. The temperature estimation using this model is compared with the experimental measurements and 3D FEA simulations to validate it.
The absolute error between the proposed model and the measurements is lower than 5 • C and the relative error ranges between −6.4% and 3.9%. This error band is considerably lower than what is considered acceptable in other research [20,36,37]. Therefore, the proposed model can be used to accurately characterize the temperature rise of any magnetic component for any boundary condition.

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