Thermomechanical Optimization of Three-Dimensional Low Heat Generation Microelectronic Packaging Using the Boundary Element Method

: This paper presents a simulation based on the boundary element method for the optimization of the thermomechanical behavior of three-dimensional microchip-dissipator packaging when the heat generation produced is medium-low. Starting from a basic architecture studied in the literature, different modiﬁcations affecting both elastic boundary conditions and the contact interface between the microprocessor and the heatsink are included and studied in order to improve heat dissipation. A nonlinear interface material is included at the interface of both solids. Thus, a thermal contact conductance as a function of the normal contact traction is simulated. Finally, all these improvements in both contact interface and boundary conditions are applied to study the maximum heat generation that this kind of architecture can efﬁciently dissipate, so that the microchip will not be damaged due to thermal deformations.


Introduction
When the heat generation in the CPU microelectronic packaging is low (<20 W), the most commonly used packaging architecture is the one shown in Figure 1 where the heatsink simply rests on the microprocessor. This situation was modeled in [1], in which the behavior of the complete heatsink-microprocessor assembly was simulated, obtaining the normal traction distribution in the contact zone and the maximum temperatures on the microprocessor. Ref. [1] concluded that above certain constant clamping pressure in the top face of the heatsink, due to the fact that the heatsink is larger than the microprocessor, very high normal contact tractions appeared at the edges of the microprocessor due to a cantilever effect. In addition, this effect is accompanied by a decrease of the normal tractions at the center which could lead to detachments at this area. Moreover, the increment of heat coming from the CPU increased these effects due to the higher thermal deformations. Thus, the microelectronic packaging architecture shown in Figure 1 was limited to low heat generation.
In order to reduce this double effect and to increase the cooling capacity of the microcelectronic packaging, in this work are studied different modifications to the architecture shown in Figure 1. Firstly, the heatsink is going to be supported on its edges in order to minimize the cantilever effect.This modification is going to constrain even more the deformations of the heatsink and the supports are going to support part of the applied clamping pressure on the heatsink, leading to a more uniform normal traction distribution at the contact zone between the heatsink and the microprocessor. Secondly, instead of applying a constant clamping pressure on the upper face of the heatsink, a retention mechanism is Mathematics 2022, 10,1913 2 of 30 simulated by means of Winkler elastic supports [2]. Thus, the clamping pressure is going to depend on the deformation of the face as a function of the elastic supports stiffness (K).
This improvement leads to a better normal traction distribution at the contact zone between both solids. Moreover, it is well known that in thermomechanical contact, the resulting interface between both solids in contact consists of numerous microcontacts and microgaps. Thus, a thermal contact resistance [3][4][5] is created since the heat flow has to go through these microgaps from one surface to the other. This thermal contact resistance increases the nonlinearity inherent in the problem as the temperature values along the contact zone becomes dependent on the traction values. In [1], vacuum conditions were assumed at the interface between both solids, so the heat flow transmission was very inefficient since it occurred only through the microcontacts. Moreover, the contact zone was very dependent on the normal contact traction which is limited by the microprocessor mechanical strength. For this reason, this work also includes a thermal interface material (TIM) fulfilling the microgaps in order to increase the thermal contact conductance as it has been demonstrated in [4,6].
There are several works focused on the study of the thermal impedance [7,8], the cooling performance [9,10] and different thermal resistances models [11,12] on microelectronic packaging. To the best of the authors' knowledge, only [1] has analyzed the thermoelastic contact in microelectronic packaging and so, this work is the first one in which enhancements of three-dimensional low heat generation microelectronic packaging are studied in order to improve its cooling capacity.
In this context, this work is focused in the study of improvements to get a better cooling capacity of low heat generation microelectronic packaging. For this reason, different modifications in both elastic boundary conditions and at the contact zone are studied. As it was stated above, in order to save the nonlinearity inherent in the problem, the BEM has been used in order to get thermomechanical variables in low heat generation microelectronic packaging. This paper is structured as follows: Section 2 presents the governing equations of the problem. The nonlinear thermomechanical contact conditions are presented in Section 3. The boundary element thermomechanical discrete equations have been described in Section 4. Section 5 is focused in the simulation of the improvements performed in the low heat generation microelectronic packaging in order to obtain a better cooling capacity. Moreover, once these improvements have been simulated, the heat generated in the microchip will be progressively increased in order to approximately determinate when the architecture herein presented stops being efficient in heat dissipation and it is necessary to find other type of solutions. Finally, the paper is finished with the summary and the concluding remarks.

Governing Equations and Boundary Conditions
The contact problem between two 3D thermoelastic isotropic bodies Ω l ⊂ R 3 (l = A, B) under quasistatic thermoelastic contact conditions with a piecewise smooth boundary Γ l , in a Cartesian coordinate system (xyz) (see Figure 2). The mechanical and thermal boundary conditions are defined using two partitions of Γ l . The thermal partition is divided into four disjoint parts: Γ l θ , on which the temperature θ is imposed, Γ l q with heat flux q prescribed, Γ l FC on which forced convection conditions θ f , h f are imposed and Γ l c , which is defined as the potential contact zone. The mechanical partition is: being the traction t i assumed on Γ l t , the elastic displacements u i prescribed on Γ l u , elastic retention mechanism modeled by means of Winkler elastic supports prescribed on Γ l K with K known and Γ l c , which is defined as the potential contact zone. The thermomechanical equilibrium equations of this problem under free volume heat sources and in the absence of body forces are: The thermal and elastic fields are coupled through the linear constitutive equation in Ω l for isotropic thermoelasticity where ε ij is defined by Mathematics 2022, 10, 1913 4 of 30 Mechanical and thermal boundary conditions are prescribed on Γ l . The Dirichlet boundary conditions are and the Neumann boundary conditions are given by with n 1 being the outward unit normal to the boundary (see Figure 2) and the heat flux q is defined as q = q i n 1i , being q i = −k t ∂ i θ.

Thermal Contact Conditions
Nonlinear thermal boundary conditions on Γ c can be defined, according to [6,16,18,25,26], for each solid as: where R TC is defined in Equation (9). The expressions above were derived from the energy balance for the contact interface Γ c (i.e., q A + q B = 0), so the contact heat flux can be defined as: When two rough surfaces come into contact under the application of certain mechanical loads, a contact joint is created, formed by several microscopic contacts and microgaps, so that the real contact area is a small fraction of the apparent area, typically a few percent [34], heat flow from one surface to another overcomes a thermal contact resistance. According to [3,4,[35][36][37][38], the thermal contact conductance is a function of several geometric, mechanical and thermal parameters such as: the surfaces roughness, the flatness of the surfaces, the surface microhardness, the thermal conductivities of both solids, the range of the contact tractions or the properties of the intermediate material.
Based on the works of Cooper [3] and works of Song and Yovanovich et al. [4,[35][36][37][38], the general form of the contact conductance can be written as: where Once thermal contact conductance is obtained from Equation (8), the thermal contact resistance is modeled in Figure 3a and calculated as: In the expressions above, it is assumed that the heat flows from one body to another only by means of conduction through the microcontacts (i.e., no interstitial conditions are considered) and no thermal interface material (TIM) is considered. Thus, to take into account both phenomena, the voids are filled with a TIM. The heat goes from one solid to the other by means of conduction through the microcontacts and throughout the microgaps (see Figure 4). According to [4,13,37,39], the gap conductance is modeled as: where the dependence of Y on the normal contact pressure is analyzed in [4,13,37,39] and it can be written as: M is a parameter that takes into account the rarefaction effects of a gas at high temperatures and low pressures. Therefore, according to [37], M = 0 if TIM is not a gas, and otherwise: where θ g = (θ A + θ B )/2. Thus, the overall thermal contact conductance considering an interface material is calculated as: and consequently, R TC = 1 h TC as it is modeled in Figure 3b. On the other hand, in elements belonging to Γ c that stay detached in the final solution, the following thermal conditions are imposed [23,27,28,40]: • Convective conditions are divided in two: -Natural convection: A constant pseudo-resistance at these regions whose value is R TC = 2/h f is considered. So the natural convection conditions are given by: • Conduction: The thermal jump depends on the separation between the pair of nodes (g n > 0). In this case, it is equivalent to consider the influence of a variable thermal resistance whose value is Finally, the boundary conditions on forced convection boundary, Γ FC , are also considered as:

Mechanical Contact Conditions
Under small displacement assumption, a common normal unit normal vector Figure 5. Therefore, the nonlinear boundary conditions on the contact zone Γ c are: Under frictionless contact restrictions, the contact traction t i on Equation (20) satisfies the Signorini's conditions on Γ c : being t n = t · n c , and g n = (g o + u n ), g o = (x A − x B ) · n c (see Figure 5) and

Boundary Element Equations
The boundary integral equations for the 3D steady-state thermoelastic problem can be written for collocation points ξ and boundary points γ on Γ l , according to Aliabadi [22], as: Equations (22) and (23) are discretized on Γ l using triangular planar elements with the node located in the barycenter (see Figure 5), leading to a system of equations such as where vectors θ l , q l , u l and t l contain the values of all nodal temperatures, heat fluxes, displacements and tractions, respectively, of solid Γ l . Matrices Q l , Θ l , H l , G l ,T l ,Q l are the coefficients matrices containing the contributions from q * , θ * , U * ij , T * ij , P * i , Q * i respectively, and are calculated as: being δ ij the Kronecker's delta, r ,i the vector r(ξ, γ) that links the application point of the load ξ with the considered point γ, n i is the component of the normal vector to the boundary, A k is the area of the field element and i, j = n 1 , n 2 , n 3 are local axis directions (see Figure 2).

Thermoelastic Contact Discrete Formulation
The interface discretization on Γ c , considering node-to-node contact, is performed for contact problems.Thus, each node on Γ l c creates a contact pair (I). Therefore, imposing  (24) and (25), the final discrete boundary equations passing the unknowns to the left-hand side, the final coupled systems for thermoelastic contact problems result in: where x l θ collects the thermal nodal external unknows (i.e., the nodal unknows that are outside Γ c . θ l c and q l c collect the thermal nodal unknowns at Γ l c . Finally, b l θ collects the known nodal thermal boundary conditions. being x l e collect the elastic nodal external unknowns outside Γ c . u l c and t l c are the mechanical nodal unknowns at Γ c . Finally, b l e collects the nodal mechanical boundary conditions known and b l θ collects the normal thermal results obtained from Equation (32). The previous systems of equations (Equations (24) and (25)) are solved following the double iterative procedure presented in [1,25].

Case Studies
The proposed formulation and the algorithms have been applied to study the thermomechanical behavior of low heat generation microelectronic packaging. The influence of simulating a retention mechanism at the top of the heatsink is analyzed. Moreover, the influence of the introduction of a thermal interface material at the contact zone between both solids on the thermomechanical variables of the packaging is also studied. Finally, the heat generation at the base of the microprocessor is increased in order to obtain the most efficient dissipation capacity for the improved architecture analyzed in this work.
The case studies presented in this work are organized as follows, starting from an initial geometry common to the three examples, successive modifications affecting the boundary conditions and/or the conditions at the interface between the microprocessor and the heatsink are proposed.
The resulting mechanical, thermomechanical and thermal results of the proposed modifications are, respectively, graphed and compared with the initial situation.
• Mechanical behavior -The differences obtained in the normal contact tractions are analyzed: tractions at the nodes located in the diagonal, in the whole contact zone and averages of normal contact tractions.
• Thermomechanical behavior -The differences obtained in the thermal contact resistance values are analyzed: tractions at the nodes located in the diagonal, in the whole contact zone and averages in thermal contact tractions.
• Thermal behavior -The maximum temperatures produced at the base of the microelectronic package are compared: temperatures at the nodes located in the diagonal, in the entire base of the microprocessor and the absolute maximum temperature.
A brief summary of what is going to be analyzed in the case studies herein presented is shown in Table 1. Moreover, additional figures are included in Appendix A. The heatsink is simply supported on the microprocessor.
• The load is applied to the top face of the heatsink by means of an elastic retaining mechanism (variable clamping pressure). • Heatsink edges are supported.

•
Vacuum conditions are assumed at the interface of the packaging.
• Thermal interface material is introduced at the interface of the packaging.

•
The heat generated at the base of the package is equivalent to a heat generation of 18.8 W. • The convective coefficient is equal to 1000 W/ • Cmm 2 .
• The heat generation at the base is progressively increased. • The convective coefficient of the cooling fluid is modified.

Effect of Retention MECHANISM at the Heatsink
In [1], the thermoelastic problem of a heatsink simply supported on the processor package was solved (see Figure 6a), obtaining the traction distribution in the contact zone and the maximum temperatures that were produced at the base of the microprocessor, when a uniform pressure was applied on the top face of the heatsink. It was observed that as a result of this arrangement, the normal traction distribution in the contact zone was not very uniform, leading to different tendencies: when the clamping pressures were low, the thermal effect predominated over the elastic one, so that these tractions were decreasing towards the edges of the encapsulation. For higher clamping pressures, the opposite effect occurred, with the normal contact tractions increasing as we moved away from the center. These non-uniform normal traction distributions affect significantly the thermal contact resistance at the interface and, in consequence, the final temperatures values at the contact zone. In addition, detachments are produced in the contact zone between the heat sink and the microchip. In this example, the problem shown in Figure 6b is solved. The aim is to propose different improvements in the elastic boundary conditions to the heatsink and package assembly in order to improve the thermomechanical behavior of the packaging.
On the one hand, the clamping pressure is not kept uniform on the top face of the heatsink so the retention mechanism of the assembly is simulated by means of Winkler elastic support with a certain stiffness. Thus, a variable traction distribution is going to be produced depending on the deformation of the upper face of the heatsink. In order to perform the study, the different stiffness values of the Winkler elastic support were obtained from the average clamping pressures applied in each case (see Table 2). The results are going to be plotted and compered with the configuration shown in Figure 6a, where equivalent uniform clamping pressures are obtained. On the other hand, the heatsink is supported on its boundary (see Figure 6a) in such a way as to limit the vertical movements of the heatsink edges in order to obtain a more uniform traction distribution in the contact area. This normal traction distribution is going to reduce the average thermal contact resistance in the contact zone, which is going to lead to an improvement in the heat dissipation.
The boundary conditions applied in the problem solved are described as follows and employed materials parameters are gathered in Table 3. The perpendicular displacements are avoided (u n 1 = 0) at the lower face of the microprocessor (z = 0 mm); -On the rest of the faces, t n 1 = t n 2 = t n 3 = 0. • Thermal boundary conditions:

-
The top surface of the aluminum heatsink (z = 19 mm) is assumed to be cooled by forced convection; The heat generated during operation of the processor is accounted for by treating the bottom surface (z = 0 mm) as a heat source (i.e., thermal gradient) with a total power of 18.8 W (gradθ = 1 • C/mm); -On the rest of the faces, gradθ = 0.
In order to apply the BEM, different meshes have been tested to obtain satisfactory numerical solutions. The influence of different meshes in the contact thermomechical variables is shown in Figure 7. The final discretization chosen is shown in Figure 8. The potential contact zone was discretized by 128 elements and the processor and the heatsink were meshed with 384 and 752 elements, respectively. Faces corresponding to planes YZ and XZ were not meshed since symmetry has been applied to those planes.  Thermal conductance has been calculated using Equation (8) and parameters gathered in Table 3, obtaining h c = 7.650 × 10 −3 |t n 1 | 0.95 W/ • Cmm 2 . Figure 9 shows the normal traction distribution for the different uniform clamping pressures and for the different Winkler elastic support stiffness. The new elastic boundary conditions included in this work improve significantly the uniformity of the traction distribution in the contact for all the cases solved. For low clamping pressures, the mean normal traction increases in the contact zone (always less than 0.7 MPa) and the thermal contact resistance decreases (see Figure 10. Nevertheless, the average normal traction decreases in the contact zone for higher clamping pressures, as the elastic traction predominates over the thermal traction. Thus, the supports located on the edges of the heatsink help to achieve a better distribution of the clamping pressure, relieving the microelectronic packaging. The stiffness of the Winkler elastic support should be less than K = 150 MPa/mm in order to avoid high tractions (greater than 0.7 MPa) that could damage the microelectronic packaging. For uniform clamping pressures lower than 0.07 MPa at the top of the heatsink, debonding occurred at the edges of the contact zone. In contrast, as the heatsink becomes longer than the processor, debonding ocurred in the center of the contact zone when very high clamping pressures are applied. However, no debonding is produced in the contact zone when Winkler elastic supports are introduced since the traction distribution becomes more uniform. It has to be mentioned that due to the singularity at the edges, high tractions would appear at these locations whose exact values cannot be collected with the method used in this work [41,42]. These large traction values at the edges are going to produce small thermal contact resistance values and do not significantly affect the solution of the problem. Thus, it has been decided, for clarity, to not include them in Figure 9.  The complete normal traction distribution of the entire contact zone is shown in Figure A1 of Appendix A.1, where the abovementioned phenomenon of singularity at the edges can be observed. Moreover, Table 4 gathers the median contact normal traction comparison for different clamping pressures and different Winkler elastic support stiffness values.  Figure 10 shows the thermal contact resistance values at the nodes located at the diagonal of the contact zone for different clamping pressures and Winkler elastic support stiffness. On the one hand, Figure A2 of Appendix A.1 shows the complete thermal contact resistance distribution comparison in the potential contact zone for the cases solved. A decrease and greater uniformity in thermal resistances is observed, resulting in better heat dissipation efficiency.

Mechanical Behavior
The median thermal contact resistance for each uniform clamping pressure and Winkler elastic support is gathered in Table 5.

Thermal Behavior
This improvement leads to a smooth descent of the maximum temperatures at the base of the microprocessor and in a more uniform distribution, as shown in Figure 11 and in Table 6.  On the one hand, Figure A3 of Appendix A.1 shows the complete temperature distribution of the base of the microprocessor. The descent and the more uniform distribution of the temperatures as the Winkler elastic support stiffness increases can be clearly seen.
Finally, the complete temperature distribution at the boundary of the microelectronic packaging for the uniform clamping pressure (P = 0.15 MPa) is shown in Figure A4 of Appendix A.1.

Influence of the Inclusion of a Thermal Interface Material on the Contact Thermomechanical Variables
In the previous example in Section 5.1 and in [1], a theoretical problem was simulated in which vacuum conditions were assumed at the interface between the microprocessor and the heatsink, so the heat transfer occurred only through the microcontacts. In reality, the microgaps are filled with an intermediate thermal material (TIM), so the thermal contact resistance decreases significantly compared to the thermal contact resistance when vacuum conditions are assumed. Thus, TIMs are used to provide an effective heat conduction path between both surfaces since the heat transfer is produced through both microcontacts and TIM. In this example, the influence of the introduction of a TIM in the thermomechanical variables in low heat generation microelectronic packaging is analyzed. Geometry used and new boundary conditions are shown in Figure 12. The heatsink is considered to be under the influence of retention Winkler elastic support with stiffness K = 150 MPa/mm and it is also considered to be supported at the edges. Moreover, mesh and material properties are the same as in the previous example. The thermal contact resistance (R TC ) takes into account the existence of different TIMs filling the microgaps by means of Equations (8), (10) and (13) using the parameters gathered in Table 7.   The clamping pressure is barely important to improve heat dissipation of the packaging when grease is introduced at the interface. This pressure should have minimum values to ensure contact between the microprocessor and the heatsink in addition to maximum values to avoid large stresses in the microelectronic packaging. Finally, Figure A5 of Appendix A.2 presents the complete normal contact traction distribution in the potential contact zone for all the different TIMs analyzed.

Thermomechanical Behavior
The thermal contact resistance values is the most thermomechanical parameter affected as it is shown in Figure 14. The thermal contact resistance values are significantly improved due to the reduction of microgaps in the interstitial zone between both solids, which are fulfilled by the TIM. Opposite to the normal contact traction values, thermal contact resistance values decrease as the conductivity of the TIM increases. In addition, Figure A6 of Appendix A.2 shows how thermal contact resistances become more uniform as the conductivity of the TIM increases, with a minor dependency with the normal contact traction values.
For the case of the microgaps filled by air, the values of thermal conductance in the microcontacts are of the same order as the values of the thermal conductance of a gas. Conversely, the thermal conductance of the microcontacts is negligible compared to the thermal conductance of the TIM when the microgaps are fulfilled by grease, as it is shown in Figure 15 and in Table 8.  Figure 15. Microcontacts thermal contact conductance (h c ) and TIM thermal conductance (h g ) % comparison for all TIMs studied at the potential contact zone.

Thermal Behavior
In relation to the thermal behavior, perfect contact is also assumed in order to analyze the efficiency of filling the microgaps with TIM (see Table 9). As a consequence of the decrease of the thermal contact resistance as the conductivity of the TIM increases, the maximum temperatures in the packaging decrease, as it is shown in Figure 16a. When a thermal grease is included in the interface of both solids, the temperature distribution at the base of the microprocessor becomes very close to the ideal situation that no thermal contact resistance is produced in the contact zone (perfect contact conditions), as it is also shown in Figure 16b. Moreover, Figure A7 of Appendix A.2 shows the complete temperature distribution of the base on the microprocessor of the microelectronic packaging for all the TIM types studied.
Defining the efficiency (E i ) as: where θ vacumm is the maximum interface temperature for vacuum conditions, θ i is the maximum temperature for each TIM analyzed and θ 0 is the maximum interface temperature for perfect contact conditions. Finally, Figure A8 of Appendix A.2 illustrates the temperature distribution on the boundary of both microprocessor and heatsink when the interface between both solids is filled by a thermal grease.

Microelectronic Packaging Behavior under Higher Heat Generation at the Microprocessor
In this example, the behavior of the heatsink-microprocessor assembly is simulated by increasing the heat generation at the base of the microprocessor in order to obtain the maximum efficient dissipation capacity for the improved architecture analyzed in previous example in Section 5.2.
As the heat generated in the CPU core increases, the normal traction in the contact zone also increases. In order to avoid high stresses on the microprocessor, the stiffness of the retention mechanism will have to be varied in order to avoid these stresses and also to ensure contact across the interface.
Geometry, mesh, boundary conditions and material properties are the same as in the previous example in Section 5.2 including a thermal grease at the interface of the microelectronic packaging, as shown in Figure 17. Firstly, for each heat generation value, the maximum stiffness of the retention mechanism will be obtained so that the normal contact traction on the microchip are close to 0.7 MPa. The cases studied are gathered in Table 10. where the heat generated at the base of the microprocessor is calculated as follows: q = k · A · gradθ. Being k and A the thermal conductivity and the area of the microprocessor, respectively. Figure 18 shows how the normal contact traction distribution changes as heat generated increases and the stiffness of the retaining mechanism decreases (see Table 10). This occurs because thermal deformations have a greater predominance than the elastic deformations. Normal contact traction distributions are initially increasing functions with a maximum at the edge and then, they reverse this tendency and tend to detachments at locations further from the contact zone center.

Mechanical Behavior
With the method herein used, it is possible to quantify the maximum stiffness that the mechanism can have so that the normal contact tractions do not exceed 0.7 MPa, as well as the limit of the heat generation that can be dissipated efficiently (without debondings at the contact zone) for this architecture, which is approximately 40 W. For higher heat generation, it is not possible to obtain a normal contact traction distribution that does not exceed 0.7 MPa in areas close to the center of the microchip.
It is observed that when medium heat generation (50 W) is simulated, the normal traction distribution is not affected by the stiffness of the system (which is very small). Therefore, it is not possible to design a retention mechanism that ensures both tractions lower than 0.7 MPa and detachments between the heatsink and the microprocessor. For these cases, it is necessary to use different configurations that dissipate better the heat generated at the microchip.

Thermal Behavior
On the other hand, Figure 19 shows the temperature distribution at the nodes located at the diagonal. It is shown that the final temperatures barely depend on the stiffness of the retention mechanism since the heat transfer is mainly produced by the thermal grease. In addition, they remain practically constant, varying a little more as heat dissipation becomes less effective. Figure A9 of Appendix A.3 presents the complete temperature distribution at the base of the microprocessor for all the heat generation values studied.  When a thermal grease is introduced at the interface between both solids, the contact is almost perfect. Therefore, the maximum temperatures could be obtained as a linear function of the microprocessor power.
where C is the coefficient that depends on the convective coefficient (h f [W/mm 2

Conclusions
In this work, a simulation based on the BEM in order to optimize the thermomechanical behavior of the microelectronic packaging for medium-low heat generation is presented. The cooling capacity of the microelectronic packaging was improved, supporting the edges of the heatsink and also replacing the constant clamping pressure applied on the top face of the heatsink by a retention mechanism simulated by means of Winkler elastic supports. In addition, a thermal interface material is introduced at the interface between both solids of the microelectronic packaging in order to improve the thermal contact conductance. The nonlinearity inherent in the problem due to the thermal contact conductance is solved by means of a double iterative procedure. The thermomechanical variables are obtained in a few iteration with the formulation proposed in this work. Following these studies, the main conclusions and remarks of this work are: • Elastic boundary conditions modifications: -Supporting the edges of the heatsink minimizes the cantiliver effects produced. Therefore, the planarity is improved and a more uniform normal contact traction distribution in the contact zone is produced; -Replacing the heatsink clamping pressure by a retention mechanism also improves the uniformity of the normal contact traction distribution in the contact zone and possible detachments at the contact zone are avoided. Moreover, as shown in Figure 10 and in Table 5, the retention mechanism reduced the thermal contact resistance for the equivalent constant clamping pressure. Thus, the cooling capacity of the complete packaging increases; • Modifications at the contact zone:

-
The thermal contact conductance is greatly improved by the inclusion of a TIM. The thermal contact resistance is reduced as the thermal conductivity of the TIM, so the maximum temperatures of the packaging decreases; -Thermal grease filling the microgaps produced a thermal efficiency of the microelectronic packaging very close to the ideal situation in which no thermal contact resistance is produced (see Table 9); -When a thermal grease is introduced at the interface between both solids, the contact is almost perfect.
Finally, the proposed method simulates satisfactorily the thermomechanical behavior of the heatsink-microprocessor assembly for the proposed improved architecture when heat generations are medium-low. However, this formulation cannot be also applied to analyze the thermomechanical contact problem for medium-high heat generation microelectronic packaging. In order to obtain a greater heat dissipation efficiency in this type of packaging, the contact area has to increase, and therefore, so does the heat dissipation area. This would be the result of incorporating an intermediate element that encapsulates the microprocessor, called Integrated Head Spreader (IHS), whose function is to capture the heat generated by the CPU core and distribute it over a wider area and then transfer it to the heatsink (see Figure 23). Therefore, it is going to exit two different contact zones. Moreover, the simulation of heatpipes running through the heatsink in order to capture the CPU heat and distribute it better throughout the block could also improve the heat dissipation efficiency.

Heatsink
Ingegrated heat spreader Thermal inter face material (TIM I )

Microprocessor
Thermal inter face material (TIM II )

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article. Furthermore, if any additional data be required, they are available upon reasonable request to the corresponding author.

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