Numerical simulation of the effect of plate spacing on heat transfer characteristics within a parallel-plate heat exchanger in a standing wave thermoacoustic system

Power can be converted with high efficiently between thermal energy and mechanical (acoustic) energy by using thermoacoustic technologies. Thus, the heat transfer characteristics are significant to the understanding of mechanisms, and improvement of efficiency for thermoacoustic devices, notably in heat exchangers. This paper introduces a two-dimensional computational fluid dynamics model of flow across a parallel-plate heat exchanger and investigates the effect of plate spacing on heat transfer characteristics. The open source CFD software OpenFOAM is applied because of the highly customizable capabilities to vary the control parameters. Firstly, the computational model including geometry, boundary conditions, equations, discretization scheme, turbulence and thermophysical properties’ models are presented, and then grid-independence validation is presented to verify the quality of mesh. The simulation results show that plate spacing influences the heat transfer between plates and adjacent area of heat exchanger, and the heat transfer coefficient goes up when the plate spacing decreases. The analysis also indicates that a possible flow transition to turbulence occurred within Re number between 247.2 and 321.4. The results in this work can help the understanding of heat transfer inside thermoacoustic system, and form a basis for future research.


Thermoacoustic systems
In general, there are two types of thermoacoustic systems, engines and refrigerators. Thermoacoustic engines can convert temperature gradient into mechanical (acoustic) energy while thermoacoustic refrigerators are able to convert mechanical (acoustic) energy back to temperature gradient. Compared to traditional engines and heat pumps, there are several advantages for thermoacoustic devices. The conventional engines and heat pumps need moving parts to work, while thermoacoustic devices do not, which means thermoacoustic engines and refrigerators can be highly reliable and less expensive. Another advantage is that the thermoacoustic systems do not use any exotic materials, which means that, in waste heat recovery applications, heat can be easier to recycle. As a result, the systems can be more environmentally friendly at reasonable Carnot's efficiencies or coefficients of performance [1]. For a typical thermoacoustic refrigerator, the acoustic driver at one end of resonator supplies power to gas medium inside the resonator, so that the gas "parcels" can execute a thermodynamic cycle within so called thermoacoustic core (comprising of hot heat exchanger, porous medium and cold heat exchanger). As gas parcels are compressed they move to the hot end of the core with the temperature increases to be slightly higher than the local solid temperature. Consequently, heat is being transferred from the gas parcel to the solid. As the flow reverses, the gas parcel expands and the temperature decreases to slightly lower than the local solid temperature, so the working gas is able to absorb heat from the cold end and the parcel returns to original thermodynamic state. With this cycle, the heat is moved from cold heat exchanger to hot heat exchanger, and the acoustic energy is transformed to thermal energy in the form of temperature difference. Conversely, a high enough temperature difference can also drive the gas parcel, and acoustic power will be self-exited and generate work.

Heat transfer in thermoacoustic heat exchanger
Although the principles of thermoacoustic systems seem straightforward, their application is still limited, partly because of the lack of understanding of flow and heat transfer across the porous structures such as stacks and heat exchangers. A number of flow phenomena has been investigated, which is helpful for a better understanding of the device operation [1][2][3][4]. In addition, there are also some research works focusing on the heat transfer. Paek et al. [5] carried out experimental studies of the heat transfer coefficients for heat exchangers in a thermoacoustic cooler. The authors found that despite the boundary layer model being widely used for thermoacoustic calculations, it did not predict the heat transfer accurately enough. However, by using a steady flow correlation and a modified Reynolds number, an accurate prediction was obtained. In another experiment, Nsofor et al. [6] indicated the relationship between the oscillatory heat transfer coefficient, the mean pressure and frequency of oscillation.
Using numerical research methods, Marx et al. [7] found that there are nonlinear temperature oscillations which are not predicted by linear models. However, after processing CFD results from a two-dimensional thermoacoustic engine, Nijeholt et al. [8] found that the nonlinear phenomena such as vortex formation and streaming flows can be captured by numerical codes. Numerical research [9] on an idealized thermoacoustic refrigerator under low Mach number indicated that the heat transfer peaks at a well-defined heat exchanger length, stack gap and distance between the heat exchangers and the stack plates. Piccolo [10] also pointed out that the fin length and spacing affect the heat exchanger behaviour and have to be optimized to decrease losses.
This work aims to investigate the effect of plate spacing on heat transfer characteristics numerically. The geometry and some parameters are based on the earlier experimental [11] and numerical [12] work.

Computational model 2.1 Physical model and boundary conditions
The geometry of simulation model is based on a thermoacoustic heat exchanger test rig, which has been described in detail by Shi et al. [11], as shown in Fig.1. The red dashed box indicates the computational domain. It is a loudspeaker-driven, standing-wave, 7.4 m-long resonator with a cross-section formed by a square duct of 134 mm × 134 mm. The working gas is nitrogen under atmospheric pressure and the resonator works in the quarter-wavelength mode with the loudspeaker frequency of 13.1 Hz. The drive radio in the rig, pressure amplitude divided by mean pressure, is 0.83%. The details of heat exchangers are shown [11].
Nitrogen driven by a loudspeaker flows across the spaces between plates of heat exchangers, and the general flow direction is parallel to x axis. Since the velocity and heat transfer is negligible along the direction perpendicular to the page, the flow is considered two-dimensional in x-y plane (cf. Fig. 1). Because the pressure and velocity of nitrogen far away from the heat exchangers can be estimated from the linear thermoacoustic theory [Błąd! Nie można odnaleźć źródła odwołania.], a "short model" (limited to the red dashed line) was applied to reduce mesh number. As a result, the boundary conditions need to be set according to the solutions from the linear model. The schematic of boundary condition is shown in Fig. 2. The pressure in left edge of domain and the velocity in right edge are defined using dynamic code so that the boundary conditions are able to change with time.
The periodic velocity on the rightmost end of domain and the periodic pressure on the leftmost end of domain are defined respectively as: where Pa is pressure amplitude at the antinode, kω is wave number, f is frequency, ρ is density and is speed of sound. Phase θ is set to follow the standing-wave criterion where pressure and velocity are 90 degrees out of phase. More details are presented by Mohd Saat & Jaworski [12]. In order to investigate the effect of plate spacing on heat transfer, the plate spacing D in this research is 3 mm, 4.5 mm and 6 mm respectively.

Simulation settings
The open source software OpenFOAM 6.0 was employed in the simulations and the mesh was generated by BlockMesh. There are 3 cases with different heat exchanger plate spacing.
With the mesh number of 140200 and max skewness of 2.3 × 10 -12 , a high quality mesh is used for accurate simulations. The nearest node from the wall, y+, for all cases is less than 0.1 after 20 cycles (all the fields exhibit regular periodic changes after 20 cycles). The simulation used buoyantPimpleFoam, a solver for transient and compressible flow with buoyancy effect. The solutions rely on the RANS equations and shear-stress-transport (SST) k-omega turbulence model, and the nitrogen thermophysical properties including density, viscosity, isobaric specific heat and thermal conductivity refer to REFPROP software based on data from NIST (National Institute of Standards and Technology) [14,15]. The second order implicit discretisation scheme, Crank-Nicolson scheme, with the off-centering coefficient 0.9 is selected for discretisation of time, and the second order upwind scheme is used for turbulent equations and transport equations. The relative tolerances for the equations are set as 10 -6 . The detailed governing equations and boundary conditions are given by Mohd Saat & Jaworski [12].

Results and discussions
For convenience, three cases studied, d=6 mm, 4.5 mm and 3 mm, are named A, B and C, respectively. The numerical simulations are started from a state of rest and are progressed until stationary conditions are reached. Convergence properties of simulations are ensured before analysing the heat transfer. Fig. 3 presents the axial velocity magnitude at point 'm' in phase 6 and 16 versus mesh number for case A (in phase 6 the velocity is a positive value, while in phase 16 it is negative). The velocities under the mesh number of 43200, 88200, 140200 and 200000 are investigated respectively, and the result shows that the model with a total of 140200 cells is sufficient to provide a grid-independence solution. In addition, the value of y+ is also analysed. The results show that the y+ is less than 1 and the maximum value oscillates regularly after t = 2.3 s, which means that the mesh thickness in the first layer next to wall meets the demands of SST k-omega turbulence model. In this study the heat transfer within one steady oscillatory cycle is investigated according to the time frame defined in Fig. 4. The cycle is divided into 20 phases, and the relationship between pressure, velocity and gas displacement for the phase is illustrated.   (1-20) for different plate spacing in hot plates and cold plates respectively. It is found that the variation of average heat flux in case A has a pseudo sinusoidal character as a function of phase while the trend is not so obvious in case B and C. However, according to Fig. 6 which shows the total space-averaged heat flux (where qtotal = qhot + qcold), the sinusoidal character is obvious, and similar results are mentioned by Shi et al. [11]. It can also be seen from the figure that the heat flux goes up when plate spacing decrease. The explanation of the phenomena can be found by analysing heat transfer coefficient (for short HTC). The HTC is defined as follow: where qavg is average wall heat flux, Tw,avg is average wall temperature, Ti is selected as the temperature at point "m" in Fig. 1 because the temperature reflects the interaction between the hot plate and cold plate. Fig. 7 shows the heat transfer coefficient versus phase 1-20 for different plate spacing in hot plates and cold plates. It is obvious that the HTCs in cases with smaller plate spacing are higher than those with larger plate spacing, both for hot plates and cold plates. This is mainly because with the same velocity from the far end, smaller plate spacing means higher flow velocity between plates and higher velocity always leads to higher heat transfer rates. For example, the peak velocity in the far end is about 1.89 m/s for all cases, but the peak average velocities between plates in case A, B, and C are 3.01 m/s, 3.38 m/s, 4.13 m/s, respectively.
It is also found that the HTC of case C is much higher than that of case A and B in hot plates from 0° to 180° and in cold plates from 180° to 360°, and the HTC of case B is only a little higher than that of case A. A possible reason for the phenomenon is that the flow transitions to turbulence in case C because of the high velocity while cases A and B are still laminar. The Reynolds numbers based on the thickness of Stokes' layer (Reδ) is defined as follows: with um, ν and ω representing the axial velocity amplitude, kinematic viscosity and angular velocity, respectively. An oscillatory flow can be considered laminar as long as the Reδ is less than a certain value, and according to Merkli & Thomann [16], the value is 400. The maximum Reδ for case A, B and C are 191.2, 247.2 and 321.4 respectively in this research, which are all less than 400; so the flow should be laminar according to the literature. However, Merkli & Thoman [16] have also mentioned that other researchers recommend different values from 150-800 based on the experiments, which means the value is not accurate for different condition. Fortunately these values have the same order of magnitude [16]. As a result, it can be estimated that the critical Reδ is between 247.2 and 321.4 based on the simulations. More research is needed to determine the exact value.
It can be seen from Fig. 5 that the curve of cases B and C are less likely a sinusoidal function. The reason is shown in Fig. 7. The peak Reynold number in all phases is high enough so that von Karman vortex street is created during last several phases, and then the vortices are "sucked" into the space between plates; as a result, the temperature boundary layer near the wall is broken and the distribution on the surface of plates becomes more non-uniform, which affects the heat flux and heat transfer coefficient. Fig. 8 shows the detailed difference between uniform and non-uniform temperature fields near the wall.

Conclusions
A two-dimensional OpenFOAM CFD model has been set up to research the effect of plate spacing on heat transfer characteristics within a parallel-plate heat exchanger, and a "short model" is applied according to previous research. The flow and heat transfer of working gas near the heat exchanger plates are analysed and the main conclusions of this paper are as follows: 1) the heat flux and heat transfer coefficient goes up with the plate spacing decrease.
2) a possible flow transition to turbulence occurred under Re number within 247.2 and 321.4, which differs from previous research and the transition affects the heat transfer on the surface of plates a lot.
3) the von Karman vortex street influence the temperature field near the plates notably, making the distribution more non-uniform. However, whether the heat transfer is enhanced or deteriorated is still unclear, which needs detailed research in the future.
This work can be of help for better understanding of heat transfer under oscillation condition. However, the presented work is of preliminary character, and further research needs to be continued.