## 2 Numerical Methodology
### Model Assumptions and Governing Equations  
The following assumptions were made to simplify the simulation model. 1. The flow and temperature fields are symmetric with respect to the geometric symmetry planes of the heat sink;
2. The cooling medium is an incompressible fluid, and its thermophysical properties are adopted at the average temperature of the inlet and outlet;
3. As the temperature of the device will not exceed 100 \({}^{\circ}\)C, the thermal radiation of the device to the external environment is negligible;
4. The influence of gravity is negligible;
5. The metallic PCM is assumed to be uniformly embedded between the contact surface of chips and the heat sink;
6. Convection heat transfer of the melted PCM can be assumed to be negligible due to the small thickness. [PERSON] et al. [29] gave a critical \(Re\) number of 980 for the axial jet and 625 for the radial jet, and the critical \(Re\) number is related to the jet geometry. [PERSON] [30] mentioned that a free jet is fully turbulent for jet \(Re\) number above approximately 1000-2000 and laminar for \(Re\) number below about 100. According to these studies, it is considered that all the jet flows are turbulent in this paper. In the current numerical calculation of engineering turbulence, the most widely used method is the Reynolds-averaged Navier-Stokes (RANS) method. Combined with the above assumptions, the following governing equations need to be satisfied for incompressible jet impingement [31; 32]. Continuity equation  
\[\frac{\partial}{\partial x_{i}}(u_{i})=0 \tag{1}\]  
Figure 1: Numerical calculation geometry models: (**a**) 3D; (**b**) 2D axisymmetric. Momentum equation  
\[\frac{\partial}{\partial t}(\rho u_{i})+\frac{\partial}{\partial x_{j}}(\rho u_{i} u_{j})=-\frac{\partial p}{\partial x_{i}}+\frac{\partial}{\partial x_{j}}\bigg{(}u \frac{\partial u_{i}}{\partial x_{j}}-\rho\overline{u_{i}^{\prime}u_{j}^{ \prime}}\bigg{)} \tag{2}\]  
Energy equation  
\[\frac{\partial}{\partial t}(\rho T)+\frac{\partial}{\partial x_{j}}\big{(}\rho \overline{Tu}_{j}\big{)}=\frac{\partial}{\partial x_{j}}\bigg{(}\frac{\mu}{ Pr}\frac{\partial T}{\partial x_{j}}-\rho\overline{T^{\prime}u_{j}^{ \prime}}\bigg{)} \tag{3}\]  
where, \(u\) is the instantaneous velocity, m/s; \(i,j\) represent three coordinate direction identifiers; \(\rho\) is the fluid density, kg/m\({}^{3}\); \(\mu\) is the dynamic viscosity of fluid, kg\(\cdot\)m\({}^{-1}\cdot\)s\({}^{-1}\); \(Pr\) is the Prandtl number of fluid; and \(-\rho\overline{u_{i}^{\prime}u_{j}^{\prime}}\) and \(\rho\overline{T^{\prime}u_{j}^{\prime}}\) are the additional terms containing pulsations generated after time homogenization, called Reynolds stress and specific turbulent heat flux, respectively. To close the above RANS equations, it is necessary to connect the additional terms of turbulence pulsation values with the homogenization equation values, that is, to introduce the turbulence model. Choosing a suitable turbulence model is very important to improve the reliability of numerical calculations. The \(k\)-\(\omega\) turbulence model in COMSOL Multiphysics software is based on the [PERSON] revised \(k\)-\(\omega\) model [33], which takes into account the low Reynolds number, compressibility and predicts the propagation rate of free shear flow. It performs well in internal flow problems and calculates more accurately flows with strong curvature, separated flow, adverse pressure gradient, jet, and other phenomena. Therefore, the \(k\)-\(\omega\) model was selected. For the near-wall surface processing, the automatic wall treatment was adopted, which can switch between low Reynolds number processing and wall function processing according to the grid resolution and obtain more accurate results with strong robustness. The apparent Heat Capacity Method is used to model phase transition processes in metallic PCM, where the PCM heat capacity is given as a temperature-dependent function [34; 35].