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

: 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 e ﬃ ciently derive a low-computing-e ﬀ ort 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-ﬂuidic 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 ﬁrst 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.


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

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

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  Table 1 supplies the thermal properties used to establish the numerical model of this BGA package.

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, 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.

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.  Following JEDEC standards, the chip behavior is characterized by a metric, which is called junction-to-ambient thermal resistance, defined by Equation (1): That metric indicates the flowing capacity of a uniform power ( ) dissipated in the device through all the thermal paths between the chip junction ( ) 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.

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.
Energies 2020, 13, x FOR PEER REVIEW 4 of 17 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.

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.  Following JEDEC standards, the chip behavior is characterized by a metric, which is called junction-to-ambient thermal resistance, defined by Equation (1): That metric indicates the flowing capacity of a uniform power ( ) dissipated in the device through all the thermal paths between the chip junction ( ) 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. Following JEDEC standards, the chip behavior is characterized by a metric, which is called junction-to-ambient thermal resistance, defined by Equation (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 (T J ) and the ambient air. This parameter can be easily calculated with measured temperatures and power. Table 2 gives the reference values of R JA used to validate the numerical models.

Definition of the Mathematical Model
Let Ω be a domain made of two disjoint sub-domains Ω 1 and Ω 2 , as illustrated in Figure 5.

Definition of the Mathematical Model
Let be a domain made of two disjoint sub-domains 1 and 2 , as illustrated in Figure 5. The boundary between each sub-domain is referred to as . An interface thermal resistance accounts for imperfect contact ( ).
Let 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 ℎ gathers convection and radiation phenomena. The thermal conductivity and the thermal capacity are respectively defined as and . Heat exchanges are modeled by the heat equation: The heat flux density is conserved through , but the imperfect contact creates a temperature discontinuity.
where is the temperature of a given sub-domain . The matrix formulation is established using classic spatial discretization by finite elements: Matrix , is a rectangular matrix that ensures the coupling between substructures i and k. U is the vector representing solicitations.

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. The boundary between each sub-domain is referred to as Γ. An interface thermal resistance accounts for imperfect contact (R c ).
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: The heat flux density ϕ is conserved through Γ, but the imperfect contact creates a temperature discontinuity.
k ∇T 1 · n 1 = −k ∇T 2 · n 2 = ϕ on Γ (4) where T i is the temperature of a given sub-domain Ω i . The matrix formulation is established using classic spatial discretization by finite elements: Matrix J i,k is a rectangular matrix that ensures the coupling between substructures i and k. U is the vector representing solicitations.

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. 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. Table 3 gives the adjusted thermal properties of the calibrated numerical model of the board substrate. 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. It occurs that the discrepancy of the numerical model ( ) in comparison to measurements ( ) 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.

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.  Table 3 gives the adjusted thermal properties of the calibrated numerical model of the board substrate. 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 R JA 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 (R N JA ) in comparison to measurements (R M 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.

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.

Modal Formulation
To ensure the coupling between both sub-domains (Ω i ), the temperature is decomposed on a Dirichlet-Steklov base [11]: Dirichlet's modes are defined as follows: where λ i 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.
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.

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 i and are expressed as a linear combination of the original modes V i according to The coefficients α ℵ i,p 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.

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.
where Λ D i and Λ S i are diagonal matrices of eigenvalues such that Energies 2020, 13, 2968 8 of 17

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.  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: 1. One "Junction": Maximum temperature of the chip, 2. One "top inner": Projected chip area on top surface, 3. Two "top outer": The four regrouped corners and four remaining top surfaces, 4. One "Bottom inner": Keep-out ball area 5. Three "Bottom outer" according to ball footprint patterns [8] 6. One "Sides": Regrouped lateral surfaces excluding the balls layer.
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.
Creation of a 3D model  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: 1.
One "Junction": Maximum temperature of the chip, 2.
One "top inner": Projected chip area on top surface, 3.
Two "top outer": The four regrouped corners and four remaining top surfaces, 4.
One "Sides": Regrouped lateral surfaces excluding the balls layer. 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. 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.

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 (2 9 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.
Energies 2020, 13, x FOR PEER REVIEW 9 of 17 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.

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 (2 9 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. 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 ℎ = ℎ = 20 /( 2 . ) and ℎ = 800 /( 2 . ). 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 h TOP = h SIDES = 20 W/(m 2 .K) and h BOTTOM = 800 W/(m 2 .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. 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.

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.

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.

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.

Steady-State Results
Two cases are presented. They highlight the component thermal behavior when: 1. Test 1: A power dissipation of 2.6 watts is applied on zone 2 2. Test 2: Zones 3, 5, and 8 are respectively submitted to a power dissipation of 0.41, 0.675, and 0.0975 watts.
The mathematical calculations assume the boundary conditions on package external surfaces presented in Table 6.

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.

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.

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.

Steady-State Results
Two cases are presented. They highlight the component thermal behavior when:

1.
Test 1: A power dissipation of 2.6 watts is applied on zone 2 2.
The mathematical calculations assume the boundary conditions on package external surfaces presented in Table 6.  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 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. 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. 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 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.
Test 2 1000 40 100 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. 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. 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 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.

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/(m 2 .K)). 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. 3.46·10 9 0.14 0-50

Caseh
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. 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 .
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. 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. 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 .
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. 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. 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. 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.

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

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

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 Energies 2020, 13, 2968 15 of 17 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].