## 2 Computational Model
The boundary conditions are calculated from the lossless equation introduced by [PERSON] [19]:  
\[P_{1}=P_{a}cos(k_{w}x_{1})cos(2\pi ft), \tag{6}\]  
\[m_{2}^{\prime}=\frac{P_{a}}{a}sin(k_{w}x_{2})cos(2\pi ft+\theta), \tag{7}\]  
Figure 1.0  
Oscillating pressure, \(P_{1}\), and mass flux, \(m_{2}^{\prime}\), are assigned at locations \(x_{1}\) and \(x_{2}\), respectively (cf. Figure 1), far enough from the heat exchangers. Condition (8) is to ensure that when the flow reverses, its temperature is equal to that of the cells next to the boundary. The wave number, \(k_{w}=2\pi f/a\), is constant because the frequency, \(f\), is constant. The terms \(a\) and \(P_{a}\) refer to the speed of sound and the pressure amplitude at the antinode, respectively. The phase, \(\theta\), is set to follow the standing-wave criterion where pressure and velocity are \(90^{\circ}\) out of phase. The adiabatic non-slip wall was assigned at the resonator walls. The heat exchanger plates were assigned temperature profiles based on experiments [17]. Figure 2 shows the interpolated profiles used in the current work. Nitrogen inside the resonator is modelled as an ideal gas. The viscosity, \(\mu\), and thermal conductivity, \(k\), of the gas change with temperature, \(T\), and the equations are given after [PERSON] et al. [24] as:  
\[ \begin{array}{c}\mu=1.82\times 10^{-5}\left(\frac{T}{T_{0}}\right)^{0.69};\\ k=-8.147\times 10^{-4}+1.161\times 10^{-4}T-1.136\times 10^{-7}T^{2}+1.062 \times 10^{-10}T^{3}-\\ 5.406\times 10^{-14}T^{4}+1.454\times 10^{-17}T^{5}-1.942\times 10^{-21}T^{6}+1.0 11\times 10^{-25}T^{7}.\end{array}  \tag{9}\]  
The influence of gravity is modelled with gravitational acceleration set to 9.81 m/s\({}^{2}\). Figure 2: Temperature profiles applied to cold and hot heat exchangers based on experimental points. CHX—cold heat exchanger; HHX—hot heat exchanger; CFD—computational fluid dynamics. Additional equations, which depend on the type of turbulence model solved, are used to solve Reynolds stresses appearing on the right-hand side of the RANS Equation (2). For all turbulence models tested, additional boundary conditions for the turbulence model are set at the boundaries \(x_{1}\) and \(x_{2}\) by assigning the value of turbulence intensity, \(I\), and turbulence length scale, \(l_{t}\), defined as:  
\[I\equiv\frac{u^{\prime}}{u_{average}}=0.16(Re)^{-\frac{1}{8}}, \tag{10}\]  
\[l_{t}=0.07D. \tag{11}\]  
The Reynolds number, \(Re=\rho u_{m}D/\mu\), is calculated using the velocity amplitude at location \(m\), \(u_{m}\), obtained from [PERSON] et al. [16; 17] and the gap between plates, \(D=6\) mm, as shown in Figure 1. The density, \(\rho\), and viscosity, \(\mu\), used in calculating Equation (10) are taken at 300 K. ANSYS FLUENT 13.0 [23] uses the inputs defined in Equations (10) and (11) to estimate the inlet distributions of turbulence kinetic energy, \(k_{e}\), and turbulence dissipation, \(\varepsilon\), using the following equations [20; 22]:  
\[k_{e}=\frac{2}{3}(u_{avg}I)^{2}, \tag{12}\]  
\[\varepsilon=C_{\mu}^{3/4}\frac{k_{e}^{3/2}}{l_{t}}, \tag{13}\]  
where \(u_{avg}\) is the characteristic flow velocity and \(C_{\mu}\) is an empirical constant that varies depending on the turbulence model used (approximately 0.09). The sensitivity of the results to these inlet conditions is tested by changing the value of turbulence intensity, \(I\), to be the value calculated using the velocity amplitude at the boundaries obtained from the equations of lossless theory. It is found that the solution is insensitive to the turbulence boundary conditions. This shows that the domain is sufficiently long to avoid the interaction between the flow around the plates and the inlet/outlet flows into/out of the domain. Default values are retained for all constants of all models used [23]. A pressure-based solver is used for all models with the application of the pressure-implicit with splitting operators (PISO) algorithm for the pressure-velocity coupling.