Study of Miniaturization of a Silicon Vapor Chamber for Compact 3D Microelectronics, via a Hybrid Analytical and Finite Element Method

The interest in silicon vapor chambers (SVCs) has increased in the recent years as they have been identified as efficient cooling systems for microelectronics. They present the advantage of higher thermal conductivity and smaller form factor compared to conventional heat spreaders. This work aims to investigate the potential miniaturization of these devices, preliminary to integration on the backside of mobile device chips, located as close as possible to hotspots. While detailed numerical models of vapor chamber operation are developed, an easy modeling with low computational cost is needed for an effective parametric study. Based on the study of the operating limits, this paper shows the thinning potential of a water filled micropillar for a device operating below 10 W and identify the corresponding vapour core height, and wick thickness.


Introduction
The generation of heat in the circuits became, in a few decades, the factor more limiting in the continuation of the improvement of the performances of the integrated devices for mobile telephony and the wandering applications. The heat, generated in the active zones of the chips, on the transistor level, is all the more concentrated that the latter are miniaturized. The heat flows generated in the last generations of silicon processors reach orders of magnitude comparable with those of the existing energy systems. The question of thermal dissipation for the electronic systems arises since their emergence. It drained an important research task through five essential technological eras: 2 Thermal management in the systems strongly integrated/wandering

Thermal resistance Model
Because of high compactness of the wandering systems such as the tactile telephones or tablets, the capacity to dissipate heat is limited to a few Watts only; the integration of heat-sinks or ventilators being incompatible with the requirement of the thickness reduction. So, this limitation results concretely in two problems of a thermal nature: the high average temperature of the system, and the presence of hot points associated to a dense logical building blocks.
The variations in temperature related to the hot points induce important variations in the performances of the devices between two zones of the same chip [11][12][13][14][15][16][17][18][19]. If it is about a sensor of images CMOS for example, the non-uniformity of temperature within the matrix of pixels can have important consequences on the final quality of the image. As example, the dark current of the pixels doubles for a rise in 6 to 8 degrees. In addition, the engineers must from now on circumvent the risks induced by these hot points by limiting the intense working lives of calculation. In the case of processors multi-hearts, for example, the activity of calculation is sometimes limited to a few tens of millisecond, time to the end of which the acceptable maximum temperature is reached. It is then necessary to stop the activity of the heart so that the system finds an acceptable temperature. The heat generated by a chip can also impact the performances of the surrounding chips in the system . Finally, the lifespan of the chips decreases because of exacerbation of the thermically activated modes of failure.
Before detailing the analytical-numerical structure of our model, we first study the relationship between thermal resistances, both for understanding the model itself and for interpreting the results. As shown in Fig. 1, heat passes through various thermal resistances when it flows from the heat source to the cooling fluid and essentially corresponds to Garimella et al. model [20][21][22]  Due to the primary heat transfer processes in each part of the heat pipe, the temperature drop can be represented using the simple effective heat resistance network as shown in Figure 1. Thermal conduction through the wall , inducing thermal resistances, in the evaporator ( , ) and condenser ( , ); implies that sections are negligible (cf. the thin thickness of wall), which is targeted at ultrathin heat pipes constructed from high-conductivity materials. Lateral conduction along the heat pipes wall is shown as , and , ; These resistances can be calculated using an effective device length, which reflects the varying heat load throughout the vaporizer and capacitor section and gives:: The resistances through the wick in the evaporator (Rw,e ) and condenser ( Rw,c) and laterally through the wick (RW, lateral), are considered as convective only. Convective heat treatment in the porous medium is negligible because of the wall in the 1-D; an efficient saturated wick thermal conductivity is used. The low-temperature, high-conductivity wick's geometry, has resulted in a high strength lateral heat flow. The interfacial phase-change thermal resistances at the wick-vapor interface in the evaporator and condenser sections, indicated by R i e and R i, c, is generally small, then neglected.
The effective thermal resistance of the Rv steam core, calculated by representing the vapor flow field as non-compressible, laminar, full flow between parallel plates, and controlled by continuum physics, to estimate the pressure fall on the core steam over effective length. This decrease is associated with a decrease in saturation temperature by applying the Clapeyron equation and the ideal gas law. According to this approach, effective thermal conductivity is defined for a 1-D lateral resistance along the effective length of vapor: where the thermo -physical properties are evaluated for saturated vapor at the local steam temperature, which changes along the device. To achieve heat loss steam properties, the 1-D steam thermal resistance is discretized and each resistor element is calculated from condenser-side temperature and is repeated in the heat flowing through the steam core. This approach does not take into account convection in the vapor core.
We then obtain an effective heat pipe resistance.

Performance-limiting Conditions
Comparing the performance between the different input heat streams and different geometries, it is possible to identify the restrictive conditions for weak thickness and low input power. Figure 2 shows the contour map of the ratio between the heat pipe and copper heat spreader thermal resistances with the adiabatic length of heat dissipation and the total thickness. The colored area shows the beneficial heat pipe performance and is whitened as the ratio transitions in favor of copper heat spreader. Performance thresholds and capillary limits are also defined. The pressure of 2250 Pa is then considered the maximum pressure, for capillarity. In Figure 2 (a) (c), show the limiting conditions for the case of a much lower input heat flux of 0,2 W / cm². At low input heat fluxes, there is sufficient capillary pressure due to the reduced fluid velocities, and the capillary limit dramatically shifts to lower thicknesses and higher working lengths. For weak thicknesses and with a relatively high temperature flux, the minimum thickness for which RHP/RHS is greater than unity, is governed by the limit of capillary; The advantage of heat pipe performance lies in the small resistances caused by the phase change. This benefit is lost for short, thick heatpipes that reduce effective length and solid heat pipe spreaders provide with a more direct thermal conductivity path than the heat pipe wick.
Except for the thickness-to-length ratio and the capillary limit, there is another limit, which is called a steam-resistance threshold. This threshold is due to the increasing pressure loss (and corresponding drop in saturation temperature) along the steam core. By replacing the different heat resistance components and representing the thicknesses of the cores, wick vapor and walls as the whole thickness fraction, ie tv = rvt, tw = rwt and twl = rwlt, the threshold conditions can be written as [21]: The solution of the equation translates into two positive roots that represent the thickness in both threshold states (vapor resistance and high ratio between thickness and length). The minimum limit thickness set by the vapor resistance threshold is: where a = 1 / (1 -rwl) and Mv is a constant that represents the properties of steam. This limit thickness is independent of the heat inlet; however, the capillary limit should be assessed simultaneously with this threshold, which should prevail over all moderate heat inputs. The properties of the steam that dictate the vapor resistance threshold are represented as a single factor that can be used as a number of merit for the choice of liquid in the design of ultra-thin heat pipes that operate at low temperature inputs: which is not dependent of geometrical parameters.

heat pipe configurations and associated equations
The standard configurations of flat heat pipes, which include wall, porous wick saturated with operating fluid, and vapor space are shown in figure 3 a, b: simple heat, multiple heat sources, underneath the heat source, and the radiator at the top of the heat pipe, multiple heat sources with symmetrical limit conditions. The basic assumptions are two-dimensional, constant material properties, constant saturation temperature and linear temperature profile across the structure of the wick.  -Conduction in wall [9] The non-dimensional temperature, heat flux and coordinates can be defined as: The two-dimensional dimensionless steady-regime heat conduction equation in the wall with constant thermal conductivity is : For the condition on boundaries (X=0 and X=1) are At the wick/wall interface (Y=0), the thermal boundary condition can be obtained assuming a linear temperature profile across the thin wick and a constant saturation temperature at the liquid-vapor interface. Then, the boundary condition are: The boundary condition at the outer wall (Y=0) are: constant heat fluxes in the evaporators and the condenser, with different locations on the heat pipe, and a heat flux null in the adiabatic sectionsb at the surface rest . That is： The non-dimensional temperature for configurations I and II can be considered as an infinite Fourier series expansion: The non-dimensional input heat flux, location-dependant of heat sources and sinks, can be developed as follows: in which : , , , , , , , , express the heat sources and heat sink locations.

-Vapor flow
A parabolic velocity profile is used for vapor flow within the heat pipe. The velocity distribution is represented by a functional product in the x-and y-directions: ( , ) = −6 ( )[ where ( ) is the local mean velocity along the x-axis.
The continuity equation, for the two-dimensional incompressible vapor flow; can be integrated with respect to y, to determine ( ) where , is the vapor interfacial velocity. The heat flux, normal to the liquid-vapor interface; , can be calculated from the conduction model, with the assumption of a linear temperature profile across the wick, getting : Integrating Eq. (20), with respect to x, and using the boundary condition at the beginning of the evaporator -= 0 −, results in the following expression for ( ): Therefore, the quasi-two-dimensional velocity profile can be written: The boundary layer of the x-momentum equation, for steady-state incompressible laminar vapor flow, can be integrated , obtaining the vapor pressure distribution: Assuming a constant vapor pressure at each cross-section, and considering the boundary conditions at the front of the evaporator, = 0 and = , yields to: where is obtained from (21).
The temperature drop in the vapor can be linked to the pressure drop in the vapor area , considering the Clausius-Clapeyron equation, and using the ideal gas law: where R is the ideal gas constant, ℎ is the latent heat, and is the pressure drop obtained from (24).

-Liquid flow
The continuity equation, for an incompressible liquid flow, can be integrated with respect to y, thus obtainingthe liquid axial velocity ( ): where , represents the interfacial velocity (normal to the liquid vapor interface) for the liquid and is related to the vapor interfacial velocity by: The axial liquid velocity can be calculated by integrating (26) with respect to x, using the velocity boundary condition at the beginning of the evaporator (x = 0), = 0: The one-dimensional steady-regime momentum conservation, for incompressible liquid flow in the wick, can be expressed by a Darcy's law, assuming negligible the inertial effects compared to viscous losses, as: in which the permeability K is calculated for a mesh screen as: The liquid pressure is derived by integrating (30) and utilizing the pressure boundary condition at the end of the condenser (x=l), = : We give also capillary pressure

FEM verification results
The schematic diagram of the simulated geometries is shown in the figure 4.
The properties of the working fluid, copper and porous wick materials, used in the vapor core chamber simulations, are shown in table1. Vapor properties are shown for a 40degC temperature. Therefore, the temperature fields predicted by the modeling approaches are shown in Figure 4. (a) shows the temperature profile comparisons along the x directions, of 50W heat source, between COMSOL and publication Maziar Agnvami [14] ; the error between them, while Figure 4 maximum temperature, at the evaporation side, in the top surface and a minimum temperature at the condenser side. Plots reveal a very good accuracy match between COMSOL and the theoretical solution.
The maximum error is about 6% for configuration I and 0.075% for configuration II. From Figure 4 (b) & (d), we observe that, by adding 10W at the heat source, the temperature gains about 3K for configuration I and 4K for configuration II. The Figure 5 compares the velocities at the center of vapor chamber.
As we can see in the figure 5, the vapor has a maximum velocity at the vapor core center and reduces to 0 when it attaches the wall. And for the COMSOL model, the maximum velocity is about 7m/s, that is almost the same for the velocities of paper [14]. For configuration I, there is about 20% difference between analytical solution and COMSOL; especially for the analytical solution, the velocity of configuration I has a non-symmetry form. But for configuration II, the difference in velocity is less than 4%. And by increasing 10W of the heat source, the velocity gains about 1.5m/s for configuration I and 6m/s for configuration II.
Finally, the most important part of the heat pipe is the capillary pressure at the vapor/wick interface. Fig  6 shows the capillary pressure at the wick/vapor interface On the steam, the water is evaporated and the steam has a maximum pressure on the surface; but next to the condenser , condensates water and creates fluid in the water, so the opposite process along with condensation and mass flow in relation to the evaporator, so that the water evaporates to the ground at condensation. It is much slower than the steam core. Figure 6 shows the COMSOL simulation and publication configuration [14], with a capillary difference of up to 20% pressure. But in the configuration II capillary pressure is almost identical. And it's less than the defect 0.5%. In the vapor chamber operation, a steady state model is developed, which allows for multiple, arbitrarily shaped, heat inputs on the evaporator-side face; model offers 2D temperature, pressure and speed ranges in the steam chamber, depending on the size of the camera, liquid and steam flow and material properties hypotheses. The governing mass, power, momentum and energy equations are set on the wall, wick and steam core domains. The model is compared to the numerical model [14].

Standalone chip
The first step consists to create a numerical model representing a standalone system [22]as shown in Figure 7. These ICs do not constitute a set of simple layers, but many small components (cf., electrical communication) that would be not easy to discretize. For the sake of practicability, these components are homogenized by substances in their environment, using the method proposed by Crécy [23]. Flat and simple heat conductors are independently referred to a representing thermal behavior in these domains. Micro heat pipe are in the Si die to the Si molding junction. The geometry of the model containing homogeneous regions is shown in Figure 8. Table 2 summarizes the thermal transitions and measurements used. The thermal conductivity is estimated by force law on the basis of empirical results [24].  The boundary conditions surrounding the chip are natural convection, with a convection coefficient estimated to hsurr=10 W/m²K [22]. The bottom of the model is normally linked to a Printed Circuit Board (PCB), which allows important lateral heat spreading, with an equivalent convection coefficient estimated to hbot=570 W/m²K [22]

Heat management
The heating power is distributed into the 8-core of 360X490 mm² to represent the test chip, as seen in Fig. 9. With a 2W per core, a total of 16 W needs to be removed. For correct operation, the chip temperature, , should be below 125 degrees. Ambient temperature and inlet temperature are set at 40 ° C, which is recommended for indoor air.

Uncooled chip
The reference case is simulated without micro-heating pipes. In this case, chip is exclusively cooled by expansion during the BGA and normal convection. In Figure 10, the reference case is simulated without micro-heating pipes. In our case, the piece is exclusively cooled by expansion during the BGA and normal convection. I expect 10 maximum temperatures above the maximum allowed temperature of 125 ° C, which indicates the need for cooling. The maximum temperature is above the maximum allowed one: 125 ° C, which clearly shows the need for cooling solution.

Conclusion
Analytical and numerical modeling, and experimental simulations of the heat pipes significantly facilitated a much better understanding of the various physical phenomena found in the heat pipes and the development of computer and experimental methods.
FEM models (COMSOL Multiphysics) are very accurate in the 2D and 3D simulation of the heat pipe. Exact limit conditions can be found in the previous section; consider the pressure factor in the calculation and specifically extract the capillary pressure from wick layer's capillary. In the course of future work, the streamlining of the evaporated material needs to be changed to temperature-dependent parameters and due to the simplification requirements; the quantitative analysis of the model error is performed in a range of operating conditions.
This model changes to a 3D model and uses multiple thermal input boundary conditions to simulate the behavior of the steam chamber to demonstrate its ability to solve the 3D thermal response to complex boundary conditions expected in real applications. Before some true applications are introduced, some transfer mechanisms will be included in the model. These include the effect of external boundary conditions on the interface mass flow profile, the vapor convection, and the flux of the interface mass in the adiabatic phase with a 2-D drive. However, these additional mechanisms should be taken into account in detailed numerical models for the application-oriented design and optimization of ultrathin heat pipes. This study highlights the relative importance of the three types of resistances in the microfluidic chip cooling and that conventional thinking is challenging if we take real module configurations into account.