## 2 Computational Model
In the first, the walls of heat exchangers are treated as adiabatic to replicate the experiments with no temperature gradient imposed on the heat exchangers (essentially ambient laboratory conditions),  
Figure 2: (**a**) A quarter-wavelength standing wave thermoacoustic rig; (**b**) computational domain of the rig covering the 7.4 m length of the resonator; (**c**) mesh generated in the area of heat exchangers and its neighbourhood; (**d**) designations/locations within the individual channel used for analysis. mainly to make comparisons between the experimental and numerical results for vortex structures. In the second, the temperature distributions were based on experiments using a linear interpolation of experimental data reported by [PERSON] et al. [17], as shown in Figure 3, mainly for validation of the CFD model itself. Further extended studies into the system behaviour under arbitrary conditions assumed \"flat\" temperature profiles (as explained later). The flow is modelled in ANSYS FLUENT [12] using two-dimensional Navier-Stokes equations as described by:  
\[\frac{\partial\rho}{\delta t}+\
abla\cdot(\rho\mathbf{v})=0, \tag{2}\]  
\[\frac{\partial}{\partial t}(\rho\mathbf{v})+\
abla\cdot(\rho\mathbf{v} \mathbf{v})=-\
abla p+\
abla\cdot(\mathbf{\tau})+\rho\mathbf{g}, \tag{3}\]  
\[\frac{\partial}{\partial t}(\rho E)+\
abla\cdot(\mathbf{v}(\rho E+p))=\
abla \cdot(k\
abla T+(\mathbf{\tau}\cdot\mathbf{v})). \tag{4}\]  
The fluid used in this study is nitrogen modelled as an ideal gas with temperature-dependent properties. A power law model suggested by [PERSON] [1] is used for the temperature-dependent viscosity while a seventh order polynomial model suggested by [PERSON] et al. [21] is selected to model the temperature-dependent thermal conductivity:  
\[p=\rho RT\ ;\mu=1.82\times 10^{-5}\left(\frac{T}{T_{0}}\right)^{0.69}; \tag{5}\]  
\[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.011\times 10^{-25}T^{7}.\]  
In Equations (2) to (5), \(\rho\), \(\
u\), \(t\), \(p\), \(\mathbf{g}\), \(E\), \(k\), \(T\), \(T_{0}\), \(\mathbf{\tau}\), \(R\) and \(\mu\) are density, velocity vector, time, pressure, gravity, energy, thermal conductivity, temperature, reference temperature (300 K), stress tensor, gas constant and viscosity, respectively [12]. Viscous heating (the second term on the right-hand-side of Equation (4)) is taken into account in this model because of the presence of heat exchangers with large surface areas. The flow is assumed to be laminar due to the small Reynolds number involved in the case investigated. The Pressure-Implicit with Splitting Operators (PISO) algorithm is selected for pressure-velocity coupling as the algorithm provides better solutions for transient cases [12; 22; 23; 24; 25]. A first order implicit scheme is used for the discretisation of time due to the presence of the moving wall [12]. The momentum and energy equations are discretised using a second order upwind scheme with the convergence set to the absolute values of \(1\times 10^{-4}\) for the continuity and momentum equation, and \(1\times 10^{-7}\) for energy equation. The density is calculated using the second order upwind numerical scheme. Density is related to pressure and temperature by using the equation of state as shown in Equation (5). The transient solution of the flow problem is solved in a segregated way using the pressure-based solver [12; 25]. The time step size was chosen so that solution converged within 15  
Figure 3: Cold (CHX) and hot (HHX) heat exchanger wall temperature. to 18 iterations in every time step. If the time step size was too large, the solution was found not to converge in every time step. If the size was too small, convergence occurred too fast within a certain time step (sometimes only requiring one iteration in one time step). The best time step for convergence was determined and set at 1200 steps per acoustic cycle.