## 2 Simulation Model and Methods
### Computational Details  
The governing equations solved by ANSYS/FLUENT CFD software, discretized by the finite-volume method and the SIMPLE algorithm, is used to obtain the pressure-velocity coupling. The spatial discretization of momentum and energy equations is obtained by using a second-order upwind scheme. The transient formulation is obtained by using a second-order implicit scheme. The criteria for convergence of the iterative solution are specified as \(10^{-6}\) for the residual of continuity, momentum, and energy.  
Figure 3 shows the computational domain of the current study. While the heat sink is moving with sinusoidal and square wave-shaped vibration, the mesh of the fluid domain must be rearranged. A dynamic mesh setup is used to obtain the mesh motion, and a layering scheme is used to regenerate the mesh.  
A mesh independency analysis is performed to ensure that the numerical results are independent of the number of cells. As shown in Figure 4, eight grids are considered in the current study with numbers of cells of 234,556, 324,568, 435,892, 498,460, 522,764, 580,356, 784,322, and 985,662. When the number of cells is increasing from 234,556 to 522,764, the Nusselt number is found to significantly decrease. However, a further increment in the number of cells leads to less than 1% of deviation in the Nusselt number. Thus, the grid with 522764 cells is used as the optimum mesh for all developed models.  
The quality of a mesh used in the simulation is indicated by skewness and aspect ratio values. The maximum skewness and aspect ratio of the mesh for the present study are 0.62 and 6.12, respectively, meeting the requirements of acceptable mesh quality.  
In addition, a time step independency investigation is carried out using three different time steps \(\Delta t=\tau/8\), \(\tau/16\), and \(\tau/32\) at a frequency of 100 Hz. The variations in the area-averaged Nusselt number of the heat sink with cycle time are shown in Figure 5. Reducing the time step (\(\Delta t\)) from \(\tau/16\) to \(\tau/32\) s gives less than 1% of variation in the peak-to-peak variation in the area-averaged Nusselt number. Here, \(\tau\) is the cycle time of each vibration case. Thus, the \(\tau/16\) s time step is used in all simulated cases.  
Validation of the CFD modeling is performed in terms of the sinusoidal vibrating heat sink, which is simulated under similar experimental conditions to [PERSON], et al. [37]. Numerical computations are performed under a range of vibrational frequencies, while the amplitude, supplied power, and air velocity are kept constant at 0.005 m, 17 W, and 0.6 m/s, respectively. As shown in Figure 6, a reasonable agreement is achieved in the validation against [PERSON], et al.'s [37] experimental measurements with a deviation of up to 8.3%. Therefore, CFD modeling in this study can be utilized with confidence to numerically analyze the heat transfer enhancement of the vibrating heat sink.  
Figure 3: Computational domain. (**a**) Heat sink, (**b**) section view of heat sink and air domain, and (**c**) air domain.  
Figure 4: Grid convergence.  
Figure 5: The variations in the area-averaged Nusselt numbers for different time steps.  
However, it can be seen that the average Nusselt number of the validation study is consistently lower than that of the experimental finding for the entire range of vibrational frequency. [PERSON], et al. [36] performed numerical and experimental investigations and validated their numerical results against the experiments they performed. Their validation results were presented in terms of temperature of the fin, and they showed that their experimental temperature graph was consistently lower than the numerically obtained temperature. In other words, when the temperature graph corresponds to a higher Nusselt number, it can be concluded that the results obtained by the [36] is in agreement in terms of the trend with our validation results.