## 2 Computational Model  
The computational domain used in this paper is based on the experimental setup of [PERSON] et al. [16; 17], sketched at the top of Figure 1. It is a loudspeaker-driven, standing-wave, 8.3 m-long resonator (of which 7.4 m is a constant cross-section square duct of 134 mm \(\times\) 134 mm). It is filled with nitrogen at atmospheric pressure and operates in the quarter-wavelength mode, with the frequency of 13.1 Hz. The pressure antinode is at the end-wall. The resonator contains a parallel-plate hot heat exchanger (HHX) coupled with a cold heat exchanger (CHX). Overall, there are 10 pairs of individual plates: five pairs are heated/cooled, and five pairs (at the top/bottom) provide a uniform porosity (flow resistance) across the cross-section. The dimensions of the plates and flow channels are shown at the bottom of Figure 1. Point'm' denotes the \"joint\" where the hot and cold plates meet. It is at 0.17\(\lambda\) from the pressure antinode. The electrical heating elements in HHX and cooling water flow in CHX aim to provide flat temperature profiles over the fins: 200 \({}^{\circ}\)C and 30 \({}^{\circ}\)C, respectively; but in practice, the heat leaks produce a slightly distorted temperature distribution as reported in [17; 21]. Particle image velocimetry (PIV) and planar laser-induced fluorescence (PLIF) techniques provide the oscillatory velocity and temperature field information around the heat exchangers. Preliminary numerical studies [22] have shown that the pressure and velocity in the flow far away from the heat exchangers can be estimated from the linear thermoacoustic theory. Thus, a \"short model\" (excluding the rest of the resonator) is acceptable if the incoming/outgoing oscillatory flow does not interfere with the plate structure [11]. The computational domain consists of two regions: the area within the plate structure and two areas outside (i.e., the open areas on the left and right from the plate structure) bounded by the resonator walls. Based on [PERSON] and [PERSON]'s correlation [13], the ratio of length to diameter (\(L/D\)) for the case investigated here should not be larger than 11.4 for flow in the channel and less than 0.52 for the flow outside. The computational domain developed based on the length of 270 mm either way from the joint between HHX and CHX (cf. the middle of Figure 1)  
Figure 1: A “short model” shown as a red-dashed box on the schematic diagram (top); computational domain (middle); and the area of the parallel-plate channel investigated (bottom). complies with these requirements. Due to flow asymmetry caused by the natural convection observed in the experiment, the full height covering all 10 pairs of plates is covered by the model. A quadrilateral structured mesh (cf. the middle of Figure 1) was used. It was refined in both the \(x\)- and \(y\)-directions to check for grid independency. The model was tested for cell counts of 47,360, 58,510, 71,460 and 85,560. It was found that the cell count of 71,460 provides a solution that is independent of the grid size and hence used throughout the current study. The minimum orthogonal quality is 0.86 with the maximum skewness of 0.13, which show that the mesh quality is very close to a good quality mesh (close to one for minimum orthogonality and close to zero for maximum skewness) [23]. The nearest node from the wall, \(y+\), for all cases is between 0.1 and 1.24. This value differs depending on the drive ratio and flow amplitude during each phase of the flow cycle; and serves as a reference for the resolution of the mesh near the wall. The models were solved using a high performance computer (4 login nodes, with 28 cores of 56 threads and 256 GB memory per one login node). The calculation for one cycle takes about 30 min. The transient models were calculated until a steady oscillatory condition is achieved, i.e., where pressure, velocity and temperature do not change with time.