*Article* **Thermo-Fluidic Characterizations of Multi-Port Compact Thermal Model of Ball-Grid-Array Electronic Package**

#### **Valentin Bissuel 1,\*, Frédéric Joly 2, Eric Monier-Vinard 1, Alain Neveu <sup>2</sup> and Olivier Daniel <sup>1</sup>**


Received: 31 March 2020; Accepted: 3 June 2020; Published: 9 June 2020

**Abstract:** The concept of a single-input/multi-output thermal network was proposed by the Development of Libraries of Physical models for an Integrated design environment (DELPHI) consortium more than twenty years ago. The present work highlights the recent improvements made to efficiently derive a low-computing-effort model from a fully detailed numerical model and to characterize its performances. The temperature predictions of a deduced ball-grid-array (BGA) dynamic compact thermal model are compared to those of a realistic three-dimensional representation, including the large set of internal copper traces, as well as its board structure, which has been validated by experiment. The current study discloses a method for creating an amalgam reduced-order modal model (AROMM) for that electronic component family that allows the preservation of the geometry integrity and shortening scenarios computation. Typically, the AROMM method reduces by a factor of 600 the computation time needed to obtain the solution while keeping the error on the maximum temperature below 2%. Then, a meta-heuristic optimization is run to derive a more practical low-order resistor capacitor model that enables a thermo-fluidic analysis at the board level. Based on the calibrated numerical model, a novel AROMM method was investigated in order to address the chip behavior submitted to multiple heat sources. The first results highlight the capability to enforce a non-uniform power distribution on the upper surface of the silicon chip. Thus, the chip design layout can be analyzed and optimized to prevent thermal and reliability issues.

**Keywords:** BCI-DCTM; ROM; modal approach; BGA; experimental validation

#### **1. Introduction**

The thermal behavior of electronic components can be finely predicted by thermal and fluidic numerical simulations [1]. However, as the geometrical and functional description of the component grows in complexity, the time needed for simulations becomes unbearable for parametric studies if fully detailed numerical representations are used.

To overcome the inherent time-consuming computation of such a detailed model, new methods for creating surrogate models able to properly reproduce its steady-state response, as well as transient ones, in a shorter time, were developed. This statement leads research, at first, toward dynamic compact thermal models (DCTMs) [2], which aim to predict the key thermal characteristics of a component. The Development of Libraries of Physical models for an Integrated design environment (DELPHI) approach promotes the use of a matrix of thermal resistances that link the sub-divided exterior surfaces of a component to its junction, which is the highest temperature of the component. The construction of a DCTM required training data, obtained by numerical simulations or experiment results.

Modal models are an alternative to DTCMs. They can be seen as an extension of the classical Fourier decomposition. The temperature is searched as a sum of known elementary spatial functions, called the modes, weighted by unknown coefficients. Different methods are based on this principle: The most popular is the proper orthogonal decomposition method (POD), which requires knowledge of thermal fields from experimental or numerical data [3,4]. However, the challenge is to find a modal base independent of boundary conditions. In that perspective, Codecasa et al. used a multi-point moment-matching method algorithm [5]. Joly et al. [6] used 'branch modes' to solve problems associated with time-dependent boundary conditions. The originality of the branch modes is that 'the branch eigenvalue problem' uses Steklov boundary conditions. This method has been extended to the amalgam reduced-order modal model (AROMM) method: A modal base is calculated by solving an eigenvalue problem. The reduced model is obtained by reducing the initial base by the amalgam method [7]. These developments gave birth to a new hybrid methodology to build DELPHI-inspired DCTMs by replacing the full detailed models by a reduced-order model based on the modal approach [8,9].

However, more complex configurations demand more sophisticated methods, in which simple reduced models are built and then connected to each other. Following that idea, Grosjean et al. developed a substructuring modal method that allows the reduction of an electronic board with several active components [10,11]. Codecasa et al. also coupled boundary-independent reduced models of components [12]. However, in those studies, the elementary components were simple (Quad Flat No-leads (QFN) package 16 or 32 with a single heat source, or Insulated Gate Bipolar Transistor (IGBT) with heat sources activated together), as the challenge was to couple those independent reduced models efficiently.

A single component with independent multiple heat sources has also been recently considered using the MPMM method [9] or AROMM [13]. Those studies have been limited to a component with a couple of independent sources. The objective of this paper is to present a reduced modal model of a ball grid array package with nine independent heat sources located on the top of the chip.

The paper is organized as follows. The position of the problem, as well as the studied material, is presented first. Then, the experimental setup and measurements are introduced. In the second step, the numerical model is built and experimentally validated for different environments when a uniform power distribution is localized on top of the chip. The creation of the reduced-order model by the modal method is then presented, and two different use cases are outlined. The first use case highlights how reduced-order modal models can be used to replace the detailed model for the creation of the boundary condition-independent dynamic compact thermal model. The second one presents a purely numerical study in which the power distribution is not uniformly applied on the top surface of the chip and demonstrates the relevance of this approach applied to this complex configuration.

#### **2. Position of the Problem**

The ball grid array package with 208 solder balls is a very popular package for integrated circuits (IC) with a large interconnection number. That kind of IC package usually consists of an active centered semiconductor chip that is glued on a two-copper-layer laminate, as described in Figure 1.

The semiconductor chip, as well as its gold wire bonds, is over-molded with a plastic resin. The diameter of the gold wires is 25.4 μm, which highlights the aspect ratio constraints.

The laminate structure routes functional signals from the input/output of the chip to the solder balls and also acts as a radiator for heat spreading. It is made up by two thin signal layers interconnected by vias. Figure 2 shows the copper patterns of the three layers of the studied ball-grid-array (BGA) substrate.

**Figure 1.** Semi-cross-section of the studied BGA208.

**Figure 2.** BGA208's laminate description.

Table 1 supplies the thermal properties used to establish the numerical model of this BGA package.


Dielectric layer FR4 0.38 Solder ball 63Sn37pb 51

**Table 1.** Reference material properties.

#### **3. Experimental Measurements**

#### *3.1. Experimental Setup*

The complex component package cannot be tested independently of a printed circuit board (PCB). Thus, a set of standardized tests [14] was performed for two standardized PCBs, named 2s0p [15] and 2s2p [16], as well as for various stabilized airflow boundary conditions (still air [17] and forced moving air [18]) in order to check the component thermal performances according to JEDEC recommendations.

Figure 3 shows the stack-up of seven layers that alternate between high- (1, 3, 5, 7) and very-low (2, 4, 6)-conductivity layers that are defined for a JEDEC 2s2p thermal test board.

The "*s*" refers to the signal layers and "*p*" to the buried power (or ground plane) layers. The two internal quasi full-covered copper layers act as efficient in-plane heat spreaders. The overall length, width, and thickness of the test board are respectively 102, 112, and 1.6 mm. The typical width of a copper trace is 300 μm.

**Figure 3.** JEDEC 2s2p thermal test boards: Copper traces definition (**a**) and stack-up (**b**).

#### *3.2. Measurements*

Experimental measurements have been performed by the thermal test services of Analysis Tech following all JEDEC *n*◦ 51 requirements. Figure 4 displays the standardized experimental setup used to characterize this BGA208.

**Figure 4.** Standardized experimental setups for natural convection (**a**) and forced convection (**b**).

Following JEDEC standards, the chip behavior is characterized by a metric, which is called junction-to-ambient thermal resistance, defined by Equation (1):

$$R\_{JA}(Q) = \left[T\_I - T\_\infty\right] / Q \tag{1}$$

That metric indicates the flowing capacity of a uniform power (*Q*) dissipated in the device through all the thermal paths between the chip junction (*TJ*) and the ambient air. This parameter can be easily calculated with measured temperatures and power.

Table 2 gives the reference values of *RJA* used to validate the numerical models.

**Table 2.** Experimental *RJA* measurement of BGA208 mounted on the 2s2p board.


#### **4. Detailed Numerical Model**

#### *4.1. Definition of the Mathematical Model*

Let Ω be a domain made of two disjoint sub-domains Ω<sup>1</sup> and Ω2, as illustrated in Figure 5.

**Figure 5.** Sub-decomposition of the domain.

The boundary between each sub-domain is referred to as Γ. An interface thermal resistance accounts for imperfect contact (*Rc*).

Let *T* be the temperature field of the domain Ω. This latter is subjected to an internal power generation, named . The generated heat is exchanged from the outside surfaces (∂Ω) considering a Fourier boundary condition. The heat transfer coefficient *h* gathers convection and radiation phenomena. The thermal conductivity and the thermal capacity are respectively defined as *k* and *C*. Heat exchanges are modeled by the heat equation:

$$\mathbf{C} \cdot \dot{\mathbf{T}} = \underline{\nabla} \cdot \begin{pmatrix} k & \underline{\nabla}T \end{pmatrix} + \phi \text{ on } \Omega \tag{2}$$

$$\begin{array}{ccccc}k & \underline{\nabla}T\_i \cdot \underline{\eta} = h & (T\_{\circ \circ} - T) \quad \text{on} \quad \partial \Omega & & \text{on} \end{array} \tag{3}$$

The heat flux density ϕ is conserved through Γ, but the imperfect contact creates a temperature discontinuity.

$$k\begin{array}{c}\underline{\nabla}T\_1\cdot\underline{n\_1}=-k\ \underline{\nabla}T\_2\cdot\underline{n\_2}=\varphi\ \text{on}\ \Gamma\tag{4}$$

$$T\_2 - T\_1 = \varphi / \mathcal{R}\_{\mathbb{C}} \text{ on } \Gamma \tag{5}$$

where *Ti* is the temperature of a given sub-domain Ω*i*.

The matrix formulation is established using classic spatial discretization by finite elements:

$$\forall \ i, k \in \{1, 2\}, \ i \neq k, \quad \mathbf{C}\_i \ \dot{T}\_i = -[\mathbf{K}\_i + \mathbf{H}\_i] \ T\_i + \mathbf{J}\_{i,k} \ T\_k + \mathbf{U}\_i \tag{6}$$

Matrix *Ji*,*<sup>k</sup>* is a rectangular matrix that ensures the coupling between substructures *i* and *k. U* is the vector representing solicitations.

#### *4.2. Experimental Validation of The Numerical Model*

To be relevant and adequate, the thermo-fluid simulations were made using a full description of every detail of the laminate structure. As seen in Figure 6, with this fine three-dimensional description of the substrate, meshing the BGA requires around 600,000 degrees of freedom (DoF) and, consequently, high computing resources.

Moreover, the JEDEC test board, described in Figure 3, is completely modeled in 3D to minimize the modeling assumptions, and the experimental setups, displayed in Figure 4, are converted to numeric boundary conditions.

**Figure 6.** BGA208's numerical model.

Table 3 gives the adjusted thermal properties of the calibrated numerical model of the board substrate.

**Table 3.** Reference properties for both JEDEC boards.


Numerical simulations were performed using four distinct pieces of computational fluid dynamic software [1,19] and demonstrates, whatever the thermal test board, a very good agreement between the experimental measurements and numerical results on the *RJA* computation, as reported in Table 4, for the 2s2p PCB.


**Table 4.** Fitting of the 2s2p thermal metrics (Icepak®).

It occurs that the discrepancy of the numerical model (*RN JA*) in comparison to measurements (*RM JA*) is lower than 2%. The 3-D numeric model of the BGA has been validated by experimental measurements, so the first step of model reduction has been achieved.

The mesh size of that realistic numerical model is 28.9 million cells. In the steady state, the convergence for each set of boundary conditions applied to that model is reached in 8h00 using 16 cores and 48 GB of RAM. Cleary, a model order reduction is mandatory to act on the design for overpopulated industrial electronic boards.

#### **5. Reduced-Order Modal Model**

Inspired by the classical decomposition in Fourier series, the temperature is searched as a sum of known elementary spatial functions, called the modes, weighted by unknown coefficients. However, the creation of the modal model is more complex and the current study focuses on the substructuring modal method [20]. This latter allows the chip to be handled separately and to reduce it more efficiently.

#### *5.1. Modal Formulation*

To ensure the coupling between both sub-domains (Ω*i*), the temperature is decomposed on a Dirichlet–Steklov base [11]:

$$T\{\underline{M}, t\} = \sum\_{i} \mathbf{x}\_i^D(t) \cdot V\_i^D(\underline{M}) + \sum\_{i} \mathbf{x}\_i^S(t) \cdot V\_i^S(\underline{M}) \tag{7}$$

Dirichlet's modes are defined as follows:

$$
\underline{\nabla} \cdot \begin{pmatrix} k & \underline{\nabla} V\_i^D \end{pmatrix} = \lambda\_i \begin{array}{cc} \mathbb{C} & V\_i^D \end{array} \text{on } \Omega \tag{8}
$$

$$V\_i^D = 0 \text{ on } \partial \Omega \tag{9}$$

where λ*<sup>i</sup>* are the eigenvalues.

Temperature fields that can be rebuilt using Dirichlet modes are null on the boundary. Therefore, these fields belong to a subspace of the admissible thermal fields, but smaller. Thus, it is necessary to add a second subspace so that the union of the eigenbasis of the two subspaces gives the space of the admissible thermal fields. This is the role of the Steklov base, whose modes verify the following eigenvalue problem.

$$\underline{\nabla} \cdot \begin{pmatrix} k \ \underline{\nabla} V\_i^S \end{pmatrix} = 0 \text{ on } \Omega \tag{10}$$

$$k\,\,\underline{\nabla}V\_i^{\rm S}\cdot\underline{\eta} = -\lambda\_i\,\,V\_i^{\rm S}\,\,\,on\,\,\partial\Omega\,\,\tag{11}$$

By construction, the union of the eigenbasis of these two subspaces gives the space of the admissible thermal fields. Thus, temperature fields can be rebuilt on the entire domain.

#### *5.2. Modal Reduction: The Amalgam Method*

The modal formulation only shifts the problem: Instead of computing temperature values at the nodes of a mesh, temporal states (or amplitudes) are searched. The next step consists of reducing the size of the model, i.e., reducing the number of degrees of freedom from *N* to *N*, where *N N*. This is done by the amalgam method, where the most prominent modes are selected and the remaining ones are added to them, weighted by a coefficient [7]. These new amalgamated modes are referred to as *V<sup>i</sup>* and are expressed as a linear combination of the original modes *Vi* according to

$$\overline{V}\_{i} = \sum\_{p} a\_{i,p}^{\Re} V\_{i,p}^{\Re}\text{ where }\Re \in \{D, S\} \tag{12}$$

The coefficients αℵ *<sup>i</sup>*,*<sup>p</sup>* are determined by minimizing, in the modal space, the distance between the modal model and a reference model. The quality of the approximation is, thus, dependent of this reference model.

#### *5.3. The State Equation*

The state equation is obtained by replacing the temperature field in Equation (6) by its modal decomposition (Equation (7)), while the test functions are the eigenmodes. A simplified version is given here, where it is supposed that the conductivity and capacity used in Equations (2) and (8)–(11) are identical, and where orthogonality properties are used to simplify the problem.

$$\dot{\boldsymbol{X}}\_i^D + \boldsymbol{V}\_i^D \mathbf{C}\_i \boldsymbol{V}\_i^S \dot{\boldsymbol{X}}\_i^S = -\boldsymbol{\Lambda}\_i^D \boldsymbol{X}\_i^D + \boldsymbol{V}\_i^D \boldsymbol{\mathcal{U}}\_i \tag{13}$$

$$\forall \forall \quad i, k \in \{1, 2\}, \ i \neq k, V\_i^S \mathbf{C}\_l V\_i^D \dot{X}\_i^D + V\_i^S \mathbf{C}\_l V\_i^S \dot{X}\_i^S = - \left(\mathbf{A}\_l^D + V\_i^S \mathbf{H}\_l V\_i^S\right) \mathbf{X}\_i^S + V\_i^S \mathbf{J} V\_k^S + V\_i^S \mathbf{L}\_l \tag{14}$$

where **Λ***<sup>D</sup> <sup>i</sup>* and **<sup>Λ</sup>***<sup>S</sup> <sup>i</sup>* are diagonal matrices of eigenvalues such that **<sup>Λ</sup>***<sup>D</sup> <sup>i</sup>* (*k*, *<sup>k</sup>*) <sup>=</sup> <sup>λ</sup>*<sup>D</sup> i*,*k* .

#### **6. Utilization of Modal Model to Create a Dynamic Compact Thermal Model**

The proposed global hybrid procedure for the creation of the dynamic compact thermal model (DCTM) is outlined in Figure 7. By coupling the model-order-reduction (MOR) technique based on the modal approach [8] and a meta-heuristic optimization [21], that procedure allows us to reduce both creation and simulation times of a suitable model of a sophisticated BGA package. The most relevant benefit is achieved for transient calculations.

**Figure 7.** Hybrid dynamic compact thermal model (DCTM) creation flow.

*In fine*, the developed reduction process enables an amalgam reduced-order modal model to be generated and then a practical dynamic compact thermal model to be derived, both being highly reliable whatever the environmental conditions. In these two cases, models are boundary conditions-independent (BCI) by construction.

The overall DCTM creation time is reduced by 86% using Reduced Order Model mathematical calculations instead of time-consuming Detailed Thermal Model numerical simulations to generate training data required for Genetic Algorithm optimization, as reported in [8].

#### *Example of DCTM Network Definition*

The derived BGA208 surrogate model is made to handle multiple thermal paths, so the DCTM network is circumscribed, in this case, to nine nodes corresponding to:


The thermal predictions of the deduced DCTM were evaluated for each boundary conditions set, as well as for each thermal board test, and then compared to experimental results, as shown in Table 5.


**Table 5.** Approximation of 2s2p thermal metrics (Icepak®).

The model agreement is good with a discrepancy lower than 4% while the simulation speed is greatly improved. Thus, the good accuracy of the DCTM permits us to integrate this model inside the system/subsystem simulation to quickly identify thermal issues and optimize cooling solutions.

#### **7. Impact of Chip Power Dissipation Layout**

In realistic applications, the functions burnt on the chip are numerous, varied, and dissymmetric, and their activations depend on used or implemented logical functions. Thus, the power distribution of the silicon chip is not uniform and additional thermal analyses need to be carried out to predict the influence of various power dissipation patterns.

As seen in Figure 8, the chip is now partitioned in nine zones (29 possible combinations) to model more accurately the heating due to the individual activation of various logical functions. The largest source and the smallest one represent respectively 22.4% and 1.8% of the chip volume.

**Figure 8.** Chip partitioning in nine zones.

In this fictive case, the conventional method used to create dynamic compact thermal models can hardly be applied [22]. First, based on the superposition principle, the number of mandatory simulations (JEDEC 38-set scenarios) to correctly identify the resistances network must be multiplied by ten. Each zone is separately activated, and then all of them.

Second, the meaning of the junction temperature for a nine-heat-sources network is not trivial. For instance, Figure 9 highlights two temperature mappings (only the chip and the copper traces are displayed) when a power dissipation of 2.6 W is uniformly applied on the upper surface of the chip (Figure 9a), or, at the opposite, concentrated on a peculiar zone (Figure 9b), named zone 7 (refer to Figure 8). Both numerical simulations assume similar boundary conditions on package external surfaces, such as *hTOP* = *hSIDES* = 20 W/(m2.K) and *hBOTTOM* = 800 W/(m2.K).

Obviously, for a smaller surface dissipation, the maximum temperature reached by the chip rises significantly (24 ◦C) and its location is not centered anymore. This phenomenon will be especially exacerbated for dynamic simulations. Indeed, the location of the maximum temperature moves at each time step following the transient power profile applied on each zone. Thus, applying our previous DCTM creation flow seems difficult, and a modal approach is chosen.

**Figure 9.** Sensitive temperatures encountered by the chip with uniform power dissipation (**a**) and localized on a peculiar area (**b**).

#### **8. Reduced-Order Modal Model for Multiple Heat Sources**

#### *8.1. Multi-Source Numerical Model*

As stated in the introduction, this part is not validated experimentally, as the experimental setup is still under development. The reduced-order model is compared to the finite elements model. However, this latter has been validated experimentally in Section 4.2 for a uniform power distribution.

#### *8.2. Computation of the Dirichlet–Steklov Base*

The BGA 208 package is split in two substructures: The chip is modeled separately by 8000 DoF. The substructuring technique allows its complete modal base to be deduced in 2.5 min. For the rest of the components (resin, balls, copper tracks), the complete base computation is not feasible, so only reduced percentages of Dirichlet modes and Steklov modes are selected, as commented in [13]. The computation of 17,800 modes (11,200 Dirichlet + 6600 Steklov) is made in 6.5 h.

#### *8.3. Reduction of the Dirichlet–Steklov Base*

In the perspective of an industrial application, reference simulations needed by the amalgam procedure should be carried out at a low computational cost and should be easy to conceive. Their objective is not to provide precise temperature fields but to trigger the relevant modes for the amalgam procedure. According to the heat sources number, ten cases are simulated and then concatenated: Each single zone is successively active and, finally, all of them. These reference simulations are obtained via a first-order Euler scheme with constant time steps. The whole process, corresponding to reference simulations and the amalgam procedure, can be performed in 1.5 h. *In fine*, 50 modes are retained for the chip, and 250 for the rest of the package, leading to a reduced modal model of order 300. Consequently, the number of DoF has been reduced by a factor of 2000.

#### *8.4. Steady-State Results*

Two cases are presented. They highlight the component thermal behavior when:


The mathematical calculations assume the boundary conditions on package external surfaces presented in Table 6.


**Table 6.** Heat transfer coefficients definition (W/(m2.K)).

Figure 10 presents the temperature field calculated by the reduced model for the whole package, as well as for the chip on the BGA substrate. The computation time of the derived ROM is lower than 0.5 s for a temporal simulation of 60 s (the time required to reach steady state). The main interest of the modal method lies in its ability to compute, at low cost, the whole temperature field, even for complex geometries such as the BGA packages family. Then, 4.5 s are needed to rebuild the whole temperature field for 101 snapshots. The hot-spot location on the chip, our concern, is properly identified, but fine details are also recovered as the temperature elevation on the copper tracks. The error between the reduced and the finite elements model is also displayed in Figure 10c. In most of the chip, and copper track, the error is below 1 ◦C (0.8%), which is a very interesting result for this preliminary investigation.

**Figure 10.** Temperature field with molded resin (**a**), without (**b**), and error field (**c**) at steady state for Test 1.

The temperature distribution at steady state for test 2 is presented in Figure 11 As the boundary conditions differ significantly from test 1 and the dissipated power is reduced, the maximum temperature reached by the chip is much lower and is predicted by the modal model with a maximum error of 0.36 ◦C (1.6%). Obviously, the hot-spot location moves as the different zones are activated, which is well predicted.

**Figure 11.** Temperature field for three active heat sources.

Finally, as the temperature field on the chip is different, it substantially affects the heat spreading on the tracks and, thus, the heat distribution on the ball array. The knowledge of the whole temperature field enables the computation of temperature gradients, and opens the way to thermomechanical consideration.

#### *8.5. Transient Results*

Dynamic simulations are conducted to compare the thermal prediction of the computed ROM with the FEM simulation (assumed to be the reference) on two test cases. FEM simulations need roughly 47 min to perform a 50 s transient simulation with multi-activations.

Two sets of boundary conditions (different from those presented in Table 6) were chosen and are summarized in Table 7.


**Table 7.** Heat transfer coefficients definition (W/(m2.K)).

Boundary conditions of Case 1 correspond to a component mounted on a PCB in vertical natural convection plus radiation. Case 2 corresponds to a component sandwiched between the PCB and thermal drain reported on top of the component (*hTOP* integrates all thermal paths: Contact resistances, thermal interface material, and aluminum drain).

Two power transient scenarios are defined: One in which different zones are successively activated, as presented in Table 8, and a second one in which different zones are simultaneously activated, as presented in Table 9. This latter describes a realistic operating case of a BGA. Indeed, power Input-Outputs and firmware are always on and the other functional areas have dynamic activation imposed by software operations.

**Table 8.** Activation of the different zones of the chip for transient simulation number 1.


**Table 9.** Activation of the different zones of the chip for transient simulation number 2.


The computation time required by the reduced model is 4.5 s, which is 600 times faster than that of the finite elements model.

The maximum temperature reached by the chip computed by the amalgamated reduced-order modal model (AROMM) and by the finite elements model has been compared for both cases. Figure 12a,b present the comparison for boundary conditions case 1 with activation profile 1 (Table 8) and boundary conditions case 2 with activation profile 2 (Table 9), respectively.

**Figure 12.** Temporal evolution of the maximum temperature reached by the ball-grid-array (BGA)–difference on this parameter between the reduced and the finite element model for test case 1 (**a**) and 2 (**b**).

The agreement on this critical parameter is very good, as it never outreaches 0.25 ◦C, i.e., a relative difference less than 1.6% for the first case and 1% for the second. A sudden rise in the difference between models is noticed when the power changes. This effect is induced by modal reduction as modes with a high time constant have been discarded.

This very good accuracy on the maximal temperature is accompanied by a satisfying precision on the entire chip. Figure 13 presents the temperature field computed by the reduced model at *t* = 37 s, i.e., at the time where the error is the most important. The maximal temperature difference on the chip between the reduced model and the FE one is less than 0.5 ◦C. On most of the chip (and the copper etches), the error is below 0.2 ◦C, yielding an average error (in time and space) of 0.08 ◦C.

**Figure 13.** Temperature field (**a**) computed by the reduced model of order 300 at *t* = 37 s. Error (**b**) with the finite elements simulation at the same time.

The error field is erratic, which is characteristic of modal reduction. Thus, the location of the maximum error cannot be known a priori with this method.

Thus, both test cases using different boundary conditions and power profiles confirm the very good accuracy of the ROM, as the error is below 2% for each time step on the chip. Those results validate the new substructuring model order reduction approach to create an AROMM of complex components.

Moreover, as modal methods compute the temperature field in its integrality, there is no need for an a priori definition of the outputs. Indeed, the localizations of the hottest spot of the chip during the transient simulation (depicted by circles) are highlighted in Figure 14.

**Figure 14.** Localization of the maximal temperature reached by the chip during the transient simulation.

Obviously, the hot spot moves as the different zones are activated. This simple fact questions the notion of junction temperature. This is confirmed by Figure 15, which compares the maximum temperature reached by the chip to the temperature at the center of the chip: An output at the center of the chip would underestimate the temperature by up to 4 ◦C, i.e., a relative error of 25%, which is by far greater than the temperature prediction error of the reduced model.

**Figure 15.** Temporal evolution of the maximum temperature reached by the BGA and the temperature at the center of the chip—Difference between those two quantities.

Another main interest of this substructuring modal approach is to have two distinct models, one for the chip and one for all the other parts of the BGA. If one of them is modified, only the latter must be regenerated. In the case of the component, all constituting parts except the die are imposed by the manufacturer, so this model is realized only once. Then, the model of the chip can be easily regenerated to take into account the new spatial power profile or correction of semiconductor thermal properties.

The substructuring modal approach offers a solution to integrate the real spatial power distribution of the component without additional creation and simulation time. Indeed, this power distribution evolves during the development cycle from the uniform power distribution to the real profile based on electric simulation.

#### **9. Conclusions**

This study presents a procedure to validate a numerical thermo-fluid model of a complex electronic component, in this case, a ball grid array package of 208 balls. Then, this detailed thermal model is used to derive a dynamic compact thermal model, inspired by the DELPHI methodology, which can be substituted by the DTM to perform a set of thermo-fluidic simulations while preserving the high level of accuracy. To quicken the reduction process, an amalgamated reduced-order modal model is coupled to meta-heuristic optimizations using genetic algorithms. As a main benefit, the overall DCTM creation time is reduced by 86%.

Further, the AROMM method coupled to a substructuring modal method is applied to a BGA208 with, this time, multiple internal heat sources. This novel method allows the building of reduced models independent of boundary conditions (BCI-AROMM). A reduced-order model of only 300 modes was built and will be improved in future works. However, the time needed to create the model remains important, as 8 h of computation were used. Nevertheless, the first deduced model offers very satisfying results as the error on the maximum temperature never outreaches 2%, as well as for steady-state and transient simulations, for a reduction factor of 600 in computation time. Moreover, this model permits us to study, quickly and accurately, all 3-D thermal phenomena involved by the complex structure of the real component. These numerical results should now be confirmed by experimental data, and an experimental setup is being conceived.

Further, a transient characterization is under investigation. The definition of the adjusted heat capacity parameters is based on one-dimensional network identification using stochastic Bayesian deconvolution [23].

**Author Contributions:** V.B. proposed this approach, performed the numerical models simulations and contributed to the paper redaction. F.J. ran the computations on the substructured model and contributed to the paper redaction. E.M.-V. designed the test vehicle, provided experimental measurements and contributed to the paper redaction. A.N. elaborated the modal substructuration theory. O.D. implemented and optimized the Genetic Algorithms. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank Vincent Fox for his contribution to this work.

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

#### **Nomenclature**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
