## 2 Methodology
### Computational Details  
This study solves the governing Equations (1)-(4) using the ANSYS/FLUENT CFD software. The finite-volume method is utilized, and the SIMPLE algorithm is employed to ensure accurate pressure-velocity coupling. A second-order upwind scheme is utilized for the spatial discretization of momentum and energy equations. The transient formulation was obtained using a second-order implicit scheme. To achieve convergence of the iterative solution, continuity, momentum, and energy convergence criteria are predefined as \(10^{-6}\).  
ANSYS/FLUENT software is unable to simulate the vibration and hence requires the development of User-Defined Function (UDF) codes, which consist of displacement  
Figure 1: Computational geometries. (**a**) conventional heat sink (**b**) thin plate-fin heat sink.  
equations of the heat sinks, to enable the software to simulate the vibration of the heat sinks and their interaction with the air domain. The horizontal sinusoidal vibration is applied to the conventional and thin plate-fin heat sink in the direction of \(Y\) (Figure 1). The conventional heat sink and the base of the thin plate-fin heat sink vibrate as a rigid body with displacement Equation (5), while the thin fins flap with the displacement Equation (6) [47; 48].  
\[Y(t)=A_{b}sin(2\pi ft) \tag{5}\]  
\[\begin{array}{c}W_{rel}(Z,t)=0.2345A_{b}\Bigg{(}\begin{array}{c}\frac{sin(2 \pi ft)-\frac{f}{\tau_{B}}sin(2\pi ft)}{\left(\frac{\tau_{B}}{\tau}\right)^{2}- 1}\\ -cosh\left(\frac{\lambda}{H}Z\right)+0.734\left(sin\left(\frac{\lambda}{H}Z \right)-sinh\left(\frac{\lambda}{H}Z\right)\right)\end{array} \tag{6}\]  
Here, \(Y(t)\) represents the displacement at a given time, \(A_{b}\) represents the base amplitude, \(t\) represents the time, \(H\) represents the height of the fin, \(Z\) represents the height between any point on the fin to the base of the fin, and \(f\) and \(f\)n represent the supplied frequency and natural frequency of the fin, respectively. \(\lambda\) represents the solutions of the equation \(\cos\lambda\)\(\cosh\)\(\lambda\) + 1 = 0. The induced flapping motion of the thin fins is based on the first mode of the vibration; thus, \(\lambda\) = 1.8751.  
A mesh sensitivity analysis is performed to ensure that the numerical results are independent of the number of cells that are utilized in each simulation model. Figure 2 illustrates the variation in Nusselt number of different grid sizes in each model. The Nusselt number initially decreases significantly as the number of cells increases. However, beyond the number of cells of 522,764 and 558,771, a further increase in the number of cells results in a deviation in the Nusselt number of less than 1% in conventional and thin plate-fin heat sink models, respectively. Therefore, the meshes with 522,764 and 558,771 cells are chosen as the optimal mesh for all simulation scenarios in the conventional and thin plate-fin heat sink models, respectively.  
Time-step sensitivity analysis is carried out using three different time steps--\(\Delta t=\tau/8\), \(\tau/16\) and \(\tau/32\)--at a frequency of 100 Hz for both simulation models. Figure 3 shows the area average Nusselt number variations in each heat sink model with the cycle. It can be seen that reducing the time step (\(\Delta t\)) from \(\tau/16\) to \(\tau/32\) s leads to a peak-to-peak variation in the Nusselt number of less than 1% in each simulation model. Here, \(\tau\) represents the cycle time of vibration. Hence, the time step of \(\tau/16\) s is utilized as a time step for both simulation models.  
Figure 2: Grid independence test results. (**a**) conventional heat sink model (**b**) thin plate-fin heat sink model.