Next Article in Journal
Nuclear Power Plant Prestressed Concrete Containment Vessel Structure Monitoring during Integrated Leakage Rate Testing Using Fiber Bragg Grating Sensors
Next Article in Special Issue
Modeling of Heat Transfer and Oscillating Flow in the Regenerator of a Pulse Tube Cryocooler Operating at 50 Hz
Previous Article in Journal
Fusion of Intraoperative 3D B-mode and Contrast-Enhanced Ultrasound Data for Automatic Identification of Residual Brain Tumors
Previous Article in Special Issue
Excitation of Surface Waves Due to Thermocapillary Effects on a Stably Stratified Fluid Layer
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Effect of Temperature Field on Low Amplitude Oscillatory Flow within a Parallel-Plate Heat Exchanger in a Standing Wave Thermoacoustic System

1
Centre for Advanced Research on Energy, Faculty of Mechanical Engineering, Universiti Teknikal Malaysia Melaka, Hang Tuah Jaya, Durian Tunggal 76100, Malaysia
2
Faculty of Engineering, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK
*
Author to whom correspondence should be addressed.
Submission received: 25 January 2017 / Revised: 19 March 2017 / Accepted: 17 April 2017 / Published: 20 April 2017
(This article belongs to the Special Issue Heat Transfer Processes in Oscillatory Flow Conditions)

Abstract

:
Thermoacoustic technologies rely on a direct power conversion between acoustic and thermal energies using well known thermoacoustic effects. The presence of the acoustic field leads to oscillatory heat transfer and fluid flow processes within the components of thermoacoustic devices, notably heat exchangers. This paper outlines a two-dimensional ANSYS FLUENT CFD (computational fluid dynamics) model of flow across a pair of hot and cold heat exchangers that aims to explain the physics of phenomena observed in earlier experimental work. Firstly, the governing equations, boundary conditions and preliminary model validation are explained in detail. The numerical results show that the velocity profiles within heat exchanger plates become distorted in the presence of temperature gradients, which indicates interesting changes in the flow structure. The fluid temperature profiles from the computational model have a similar trend with the experimental results, but with differences in magnitude particularly noticeable in the hot region. Possible reasons for the differences are discussed. Accordingly, the space averaged wall heat flux is discussed for different phases and locations across both the cold and hot heat exchangers. In addition, the effects of gravity and device orientation on the flow and heat transfer are also presented. Viscous dissipation was found to be the highest when the device was set at a horizontal position; its magnitude increases with the increase of temperature differentials. These indicate that possible losses of energy may depend on the device orientation and applied temperature field.

1. Introduction and Literature Review

Thermoacoustic systems are usually divided into engines and refrigerators depending on the direction of energy conversion between the acoustic and thermal energy. The working principle of the thermoacoustic refrigerator is illustrated in Figure 1. The acoustic driver at the end of the resonator supplies acoustic power, W a c , to the working gas inside the resonator. The standing wave acoustic field induces oscillatory motion of “gas parcels” coupled with their cyclic compression and expansion. Internal structures such as heat exchangers and stacks are placed at a location where the oscillating pressure and velocity are non-zero. A fragment of a single stack plate has been magnified at the bottom of Figure 1 together with a gas parcel undergoing an acoustic oscillation. The gas parcel is compressed as it moves to the right and the parcel temperature increases to be slightly higher than the plate temperature. Consequently, heat is being transferred from the parcel to the plate. As the flow reverses, the gas parcel expands and the temperature drops slightly lower than that of the plate. As a result, heat is being absorbed from the plate and the parcel returns to its original thermodynamic state. This completes the thermodynamic cycle. In essence, the acoustic wave provides power that allows the heat to be pumped up the temperature gradient. Conversely, if a high enough temperature gradient is imposed across the stack/regenerator, an acoustic power will be self-excited and useful energy will be produced [1]. This in turn could be extracted by a linear alternator to produce electricity, or used to drive a coupled thermoacoustic cooler.
The challenges in achieving efficient thermoacoustic devices were discussed comprehensively by Swift [1]. The difficulties with complex flows across the porous structures (stack, regenerator, and heat exchangers) that form the core of the thermoacoustic device have been raised as important research issues. Understanding phenomena such as vortex flow at the end of stacks [2], entrance effects [3], joining conditions [1], and streaming [4], is important as these could hold the key to an efficient device operation.
The structure of the vortex pattern, as observed for example by Mao et al. [5], is expected to influence the flow and possibly induce losses that may affect the efficiency of thermoacoustic systems. This has been illustrated through the observations of streamlines and viscous dissipation presented in the numerical investigation of Worlikar and Knio [6]. The viscous dissipation represents mechanical energy losses within the model. The choice of location for the stack/heat exchanger and blockage ratio (gap between plates) were reported as factors that influence mechanical energy losses.
The same group has made several improvements to their model [7,8,9]. These involve adiabatic stacks and heat exchangers, all working in a thermoacoustic standing wave environment. In all models, a simplified periodic configuration was used covering the area between the plates and some open area next to the plates to include the possible contribution from the vortex shedding phenomenon. The natural convection was assumed to be very small compared to the oscillatory flow magnitude, and was hence neglected. The effect of natural convection in an oscillatory flow has been previously tested in an experiment involving a heated wire located in an empty resonator [10]. The natural convection was found to dominate the flow under a small velocity regime. This result, however, is not very conclusive because the temperature difference investigated was not widely varied. Furthermore, it is speculated that the presence of additional structures such as the stack and heat exchanger may alter the influence of the natural convection effect on the flow even when the velocity is relatively small [11].
It is important to highlight that most numerical modelling work presented in this review applies a simplified model whereby a thin plate or a pair of adjacent plates with an implementation of periodic boundary conditions to replicate the array of plates is considered. In modelling, a periodic boundary condition means that the flows at the matching periodic boundary are linked so that flow conditions are shared at this boundary. As a result the flow will be duplicated for all other channels in a periodic manner [12]. While this has advantages in terms of computational cost, it must be carefully applied as it might give misleading information pertaining to the flow structure. This has been demonstrated experimentally by Guillaume and LaRue [13]. Their experimental results were later compared to a numerical model developed using commercial software, Fluent [14], using the same number of plates as in the experiment. Their works show that, for a structure with an array of plates, the flow structures at the end and in the middle of the plates are influenced by the adjacent plates. The Strouhal number (the ratio of a product of a characteristic length and the frequency of the vortex shedding to a reference velocity) is reported to have changed in magnitude when a comparison was made between the two arrangements. It is important to highlight here that the investigations reported in [13,14] were carried out in steady flow and might not be directly applicable to the oscillatory condition, but they do give an idea of possible sources of discrepancy if investigations of an array of plates are carried out using single plate approximation. The observation of asymmetry and non-periodic development of vortex structures at the end of a plate array was also reported by Mao et al. [5]. This suggests that a physical periodic structure (array of plates) does not guarantee that the flow structures are periodic. Hence modelling a single plate may not be adequate to describe the physics of flow around a structure with an array of plates.
Similarly, considering a single plate, the flow structures at two ends may not be identical in practice for reciprocating flows. Most investigations (cf., [7,9,12]) assumed that flows at the two ends are actually symmetrical and therefore the analysis was mainly focused on just one end. This seems reasonable for investigations that do not involve a temperature gradient (adiabatic stack). However, the presence of a temperature gradient on solid boundaries is necessary for the thermoacoustic effect to take place. Temperature effects must be considered to account for additional phenomena such as temperature-driven buoyancy [11,15,16,17,18] and non-linear effects [4] that could affect the flow symmetry. Experimental studies often reveal signs of nonlinearities that are insufficiently addressed by the linear model [1,11]. In most numerical studies (cf. [7,9,18]) the natural convection effect was neglected. This effect however was observed in many experimental works [1,11,17].
Current computational fluid dynamics (CFD) study aims to improve the understanding of physical processes involved in the phenomena reported in the earlier experimental work of Shi et al. [17]. However, bearing in mind more complex operating conditions (e.g., tilted-angle solar powered systems [19,20] or space applications), the scope of CFD is somewhat wider and extends to include gravity effects and effects of the device orientation on the physics of thermo-fluid processes across the heat exchanger plates. The novelty of this study stems from the fact that the model used attempts to take into account a far wider range of physical effects which are commonly neglected in numerical works related to thermoacoustics (i.e., natural convection or the use of a full array of plates to handle asymmetrical flow features). This is undertaken to explain the phenomena reported in [11,17]. In addition, the study of device orientation provides a new perspective of the importance of considering such seemingly minor details in modelling heat transfer phenomena in thermoacoustics.

2. Computational Model

A two-dimensional computational domain based on the experimental setup of [11,17] is shown by the red-dashed box in Figure 2a. It covers a rectangular area of 0.132 m × 7.4 m of the quarter-wavelength resonator. The mesh is designed to be denser within the vicinity of the heat exchanger (300 mm to the left and right from the heat exchangers) to resolve vortex structures issuing from and returning into the heat exchanger assembly. Elsewhere, a much coarser mesh is adopted as the problem consists of solving the governing equations for an empty resonator. The coarse mesh is illustrated in Figure 2b and the fine mesh in Figure 2c. This approach was validated against the theoretical predictions of linear thermoacoustic theory [1]: 300 mm away from the heat exchanger, the pressure and velocity simply follow the linear model. A structured mesh of quadrilateral type was used for this model; the long computational domain is divided into several parts so that the distribution of mesh sizes could be made with acceptable quality. The minimum orthogonal quality was recorded to be 0.74 with maximum skewness of 0.31 (i.e., close to 1 for minimum orthogonality and close to 0 for maximum skewness [15]). The minimum mesh spacing in the x- and y-directions is 0.08 mm and 0.25 mm, respectively.
The plates of the hot and cold heat exchangers (subsequently denoted as HHX for “hot heat exchanger” and CHX for “cold heat exchanger” throughout the paper) have thickness, d , of 3 mm and are arranged in a parallel configuration with the spacing, D , of 6 mm between them. This is shown in Figure 2d as a magnified view of a single channel between a pair of hot and cold plates; point ‘m’ is the location of the joint where the hot and cold plates meet. The HHX and CHX are located next to each other with a joint positioned at 0.17 λ from the pressure antinode. The wavelength, λ , is defined as λ = a / f , where a is the speed of sound and f is the frequency of the flow (13.1 Hz). Location of the pressure antinode for this 1/4-wavelength rig is at the end wall—cf. Figure 2a. The arrangement of the heat exchanger plates (10 pairs) follows the physical setup with five pairs of heated/cooled fins and five pairs of unheated ones, treated as adiabatic plates. The unheated plates provide consistent porosities for the flow (cf., [17]).
The model extends over the whole height of the resonator to resolve the flow asymmetries observed experimentally [17]. The use of symmetry or periodic boundary conditions is inappropriate due to the presence of temperature-driven buoyancy effects. In addition, the experimental data are based on phases determined by comparing the phase of the velocity between the plates to the phase of the pressure at the pressure antinode. As the intention was to replicate as closely as possible the phenomena taking place in the physical apparatus, a full two-dimensional model was chosen.
The boundary condition at the right end of the resonator is defined as a moving wall to replicate the acoustic displacement induced by the loudspeaker. The displacement, δ , is simply modelled as:
δ = [ p a ω ρ m a sin ( k w x s ) ] cos ( ω t ) .
Here, ω , ρ m , a , k w , and t are the angular velocity, mean density, speed of sound, wave number, and time, respectively. The term p a is the oscillating pressure at the location of the pressure antinode and x s is the distance from the joint to the pressure antinode (shown as 4.6 m in Figure 2a). The mean pressure, p m , is 0.1 MPa. The drive ratio, p a / p m , of the low amplitude flow investigated throughout this paper is 0.3%. The resonator wall was treated as adiabatic. The model was solved for two thermal conditions. 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), 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 Shi 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:
ρ δ t + · ( ρ v ) = 0 ,
t ( ρ v ) + · ( ρ v v ) = p + · ( τ ) + ρ g ,
t ( ρ E ) + · ( v ( ρ E + p ) ) = · ( k T + ( τ · v ) ) .
The fluid used in this study is nitrogen modelled as an ideal gas with temperature-dependent properties. A power law model suggested by Swift [1] is used for the temperature-dependent viscosity while a seventh order polynomial model suggested by Abramenko et al. [21] is selected to model the temperature-dependent thermal conductivity:
p = ρ R T   ; μ = 1.82 × 10 5 ( T T 0 ) 0.69 ;   k = 8.147 × 10 4 + 1.161 × 10 4 T 1.136 × 10 7 T 2 + 1.062 × 10 10 T 3 5.406 × 10 14 T 4 + 1.454 × 10 17 T 5 1.942 × 10 21 T 6 + 1.011 × 10 25 T 7 .
In Equations (2) to (5), ρ , v , t , p , g , E , k , T , T 0 , τ , R and μ 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 × 10−4 for the continuity and momentum equation, and 1 × 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 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. The area-weighted-average of pressure at the end wall, known as the pressure antinode, was monitored until a steady state oscillatory flow was obtained. This is defined as a state when pressure, velocity, and temperature do not change much from cycle to cycle. By way of example, Figure 4 shows the history of oscillating pressure and velocity monitored for location “m” as defined in Figure 2. A steady oscillatory flow condition is achieved after six flow cycles. However, as will be discussed later, results presented in this paper are obtained after 70 flow cycles, when both the flow and thermal oscillatory flow conditions have reached a steady oscillatory state. The solutions are also monitored so that they converge in every time step at every flow cycle (through the selection of the time step size).
When the steady oscillatory state condition is reached, the model is validated and used for analysis. In this study, the development of flow within one oscillatory cycle is investigated according to the time frame defined in Figure 5. The cycle is subdivided into twenty phases starting from the first phase, ϕ 1. The relationship between pressure, velocity, and gas displacement, for the 20 phases of a flow cycle is illustrated. Phase ϕ 1 is set for the maximum value of the oscillating pressure at the location of the pressure antinode (rigid wall of the resonator). The period of one cycle is the inverse of the operating frequency—13.1 Hz.
The model was solved using the SPECTRE High Performance Computing facility at the University of Leicester (4 login nodes for solving one model, with 28 cores of 56 threads and 256 GB memory per one login node). With this computing facility, the 2D unsteady calculation for one cycle takes about 30 min. Therefore the solution for 490 cycles (as discussed later) takes about 9.5 days.

2.1. Grid Independency

The grid independency tests are illustrated in Figure 6. In essence, the mesh was progressively refined in the part of the domain surrounding the heat exchangers (there was practically no need to do this in an empty resonator). The description of mesh density is given as cell counts in Figure 6a for five cases. Here, the velocity values at point “m” for phases 6 and 16 are taken as test variables (+ve and −ve values due to the change of direction). Similarly, the resulting axial velocity profiles (from the wall to channel centre), for selected phases in the cycle are further illustrated in Figure 6b. For clarity, velocity profiles are shown only for three cell counts. At low mesh density, the velocity profile is slightly over-predicted at all 6 phases shown. The velocity within the boundary layer appears to be the same for all cases because the mesh is always designed to be denser in that area. It is found that the model with a total of 45,910 cells is sufficient to provide a grid-independent solution.

2.2. CFD Model Validation

The velocity amplitude obtained from the converged model that reaches a steady oscillatory state is compared to the experimental result to validate the model. The velocity amplitudes shown in Figure 7a,b are for the location at the centre of the gap between the plates of the heat exchanger, 1 mm away from the joint—cf., Figure 2c—for all 20 phases in a cycle. In both graphs, the line denoted as “Theory_dT = 0” represents the theoretical/analytical solution for oscillatory velocity calculated from the general thermoacoustic theory as given by [1]:
u 1 = i ω ρ m [ 1 h υ ( y , z ) ] p 1 x ,
where ω , ρ m , and p 1 / x are the angular velocity, mean density, and the pressure gradient of the oscillatory flow, respectively. The term h ν is a viscous shape factor which varies according to the geometry of the internal structure involved in the system. The shape factor for parallel plate geometry is given as [1]:
h κ , ν = cosh [ ( 1   +   i ) y / δ κ , ν ] cosh [ ( 1   +   i ) y 0 / δ κ , ν ] .
The shape factor defines the viscous, ν , or thermal, κ , effects depending on the definition of the penetration depths used. The terms: y , y o , and δ κ , ν are the distance from the centre of the gap to the wall of the plate, the centre of the gap, and the thermal ( δ κ = 2 κ / ω ) or viscous ( δ υ = 2 υ / ω ) penetration depth.
Figure 7a, for the adiabatic case, compares the experimental results of Shi et al. [11] (“Experiment_dT = 0”) with the current numerical results (“CFD_Laminar_dT = 0”). The theoretical solution is calculated for the mean density taken at 300 K. A good match between the experiment, CFD, and analytical solution can be seen. Maximum discrepancies between the three methods appear at phases ϕ 6 and ϕ 16 representing the highest velocity values for the two parts of the flow cycle. Comparison between the theory and experiment shows that the theoretical value is 9.3% lower than the experimental value at ϕ 6, and 2.6% lower at ϕ 16. However, CFD results differ from the experimental results only by 0.3% and 1.8% for ϕ 6 and ϕ 16, respectively.
Figure 7b shows the velocity oscillations obtained from the experiment and CFD for the case with imposed temperature difference (the theoretical adiabatic solution is only plotted for reference). Experimental and CFD results show that, with the imposed wall temperature condition, the flow in the second half of the cycle exhibits a larger amplitude of velocity compared to the first half of the cycle. This is counterintuitive since one would expect that a hotter fluid would travel with higher velocity in the first half of the cycle, while cold fluid would flow with a lower velocity in the second half of the cycle. However it is possible that there is a mean flow along the channels (streaming effect) due to the convective currents as discussed later in Section 3.2.2.
Figure 8 shows the validation of temperature profiles for a location 15 mm from the joint—cf., Figure 2c—above both the cold and hot plates. The maximum relative discrepancies of 6% were found between the numerical and experimental values when using the absolute temperature scale (17 K on the HHX and 18 K on CHX). However these discrepancies need to be seen in the context of high experimental uncertainty of the temperature field when using acetone based planar laser induced fluorescence (PLIF) method (16 K reported in [17]) as well as a large temperature span between the heated and cooled plates—cf. Figure 3. Clearly, the trend of the temperature profiles agrees very well qualitatively with the experimental values. Indeed, the original PLIF experiments were more concerned with temperature gradients than absolute temperature values in order to calculate surface heat fluxes.
Of course, even if the temperature measurement had zero error, there would still inevitably be discrepancies in CFD arising from the inability to model all aspects of the physical process. For example, in reality the walls of the rig are not adiabatic, and in addition they do accumulate heat and tend to induce thermal inertia processes with the scale of hours. The role of the three-dimensional flow effects within the system is also unclear—unfortunately, the experiment could only provide “planar” results close to the resonator centreline.

3. Results and Discussions

Results presented in this paper look at four aspects of CFD modelling: the choice of temperature field initialisation, the effects of wall temperature on the flow field, device orientation relative to gravity force and viscous dissipation.

3.1. Investigation of the Effect of the Initial Fluid Temperature on Flow and Heat Transfer

The current CFD study deals with flow fields with relatively high imposed temperature gradients which results in difficulties to achieve numerically a temperature field similar to that obtained from the experiment. A possible way of bringing these close to each other, within a reasonable computational time, was through setting an appropriate initial fluid temperature in the model. Three approaches were investigated and these are discussed below. Their outcomes are illustrated for phase ϕ 1 in Figure 9; the experimental data shown in Figure 9a are compared against CFD data in Figure 9b–d.
First, the model was initialised with a uniform fluid temperature, T i = 300   K . The model was iterated until it reached 70 flow cycles which is equivalent to 5 s of real time. Figure 9b shows that the resulting temperature field is still dominated by the low temperature of the initial 300 K. On the other hand, the experimental data in Figure 9a show that the average fluid temperature is considerably higher due to the fact that the experimental data was collected after the flow has conveniently settled in the steady oscillatory state. At that stage, the heat has been accumulated in the system giving rise to the mean temperature of the flow. For the computational model with initial temperature of 300 K, a significant amount of computational time and effort would be required if a similar heat accumulation as in the experiment were to be achieved.
Second, the model was initialised with a uniform fluid temperature close to the mean temperature from the experiments, T i = 360   K and was similarly iterated for 70 cycles. Although the resulting temperature field in Figure 9c is closer to the experimental results, the effect of natural convection in the open area next to the CHX is not visible in the computational result. In addition, the temperature of the fluid within the channel of the CHX reduces from the initial high temperature to a lower temperature due to the cooling effect at that location. The reason for the mismatch is not clear, especially as the wall of the heat exchanger is already modelled following the temperature measured by the thermocouples as reported in the experimental work [17]. However, it is possible that additional features need to be introduced into the model such as heat losses (the modelled resonator is adiabatic, in reality there are considerable heat losses), heat accumulation, or heat leakages. The known issues of experimental temperature uncertainties need to be borne in mind too.
The third approach tried to improve the predictions of natural convection observed experimentally at the plate ends by initialising the model with an approximate temperature field estimated from the experimental results ( T i T ( x , y r ) ). The resulting temperature contour after 70 cycles is shown in Figure 9d. The effect of the initialisation temperature seems to dominate the temperature field especially at the open area next to the plates. In the area between the plates of the heat exchanger, the wall temperature seems to control the temperature within the area.
Figure 10 shows the effect of initial temperatures on the velocity amplitudes at a location 1 mm away from the joint above the cold plate. Evidently, a higher mean temperature of 360 K results in a slightly higher amplitude of the velocity, especially during the second half of the cycle. This may suggest that the experimental temperature field is affected by the accumulation of heat within the investigated area since raising the initial temperature from 300 to 360 K brings the solution closer to the experimental results. The same observation can be made for all locations along the heat exchanger plate as illustrated in Figure 10b. The velocity amplitudes at ϕ 6 and ϕ 16 shown in Figure 10b represent the highest velocity amplitude during the first and second parts of the flow cycle. The velocity amplitude for the model initialised at 360 K is always higher than the other two models for all the locations along the heat exchanger plate.
Finally, the effect of initial temperature on the heat flux is shown in Figure 11. Here, the heat flux is spatially averaged over the length of the heat exchanger (the same way as in [8]), and is presented as a function of time (phase in the cycle):
q h , c = 1 l 0 l k d T d y | w a l l d x .
Subscripts h and c refer to HHX and CHX, respectively. The results were taken at cycle 70 after the model reached a steady oscillatory condition.
Overall, the space-averaged wall heat fluxes from CFD are in-phase with the experimental values at both heat exchangers. However the magnitude differs at the HHX, especially when the flow starts slowing down at ϕ 6 and changes direction. The numerical space-time-averaged wall heat flux is calculated to be 30% lower than the experimental values. However, the difference is within the range reported in other numerical investigations [26]. The wall heat flux predicted by the model initialised at 360 K seems closest to the experiment. Also, the predictions for CHX seem to follow the experiment better. The above discussion shows that the choice of initial temperature field is an important factor in the model, and it may lead to difficulties with convergence and accuracy.
Figure 12 shows the temperature profiles plotted along the vertical direction y r , between the bottom and top wall of the resonator, at locations in the open areas: 40 mm from the cold end and 38 mm from the hot end of the heat exchangers, that is where most of natural convection effects are identified experimentally. Clearly, initialising the temperature at 360 K (dotted lines) gives the wrong values of temperature especially at the lower area of the domain compared to the experiments. Initialising the temperature field as T ( x , y r ) (dashed lines), seems to give a closer match in this area.
Finally, it is interesting to look at the development of the temperature field as a function of the number of cycles taken in the solution. Here the model initialised at 300 K is taken as an example; the temperature profiles are plotted in Figure 12 for increasing the total number of cycles as coloured lines between the dash-dotted black line (45 cycles) and solid black line (490 cycles). The temperature profiles seem to tend towards the experimental values. However, the changes become smaller and smaller (especially after about 70 cycles) and it is unlikely that the profiles would ever reach the experimental data. Following this observation, all the models used in the analysis reported in this paper were run for at least 70 cycles so that the steady oscillatory thermal condition is reached from a practical point of view, in addition to the steady oscillatory flow condition reported in Figure 4.
Figure 13 shows the temperature field within and around the heat exchanger assembly for the initial temperature of 300 K after 490 cycles. This illustrates the convective currents at both the hot and cold ends to supplement Figure 12. The heat accumulation phenomena shown in Figure 13 were also seen in [8]. The results illustrate a potential source of streaming in thermoacoustic systems [1] caused by buoyancy driven flows, which may interfere with the main thermoacoustic processes.
Figure 14 shows the temperature profiles plotted from the heat exchanger wall to y = D / 2 (centre of the channel between the plates) at a location 10 mm from the joint above both the hot and cold plates. These are plotted for selected phases in the cycle ( ϕ 6, ϕ 8, ϕ 12, and ϕ 16) for the experiment and three temperature initialisation cases, after 70 cycles. The top four plots are for the cold plate; the bottom four plots are for the hot plate. An initial temperature of 360 K leads to profiles with magnitudes larger than the experiment and it was already shown that it models the natural convection in the open areas poorly. The use of temperature distribution T ( x , y r ) provides solutions almost similar to a model initialised at 300 K. Both offer temperature profiles closer to the experiment. However, from the practical perspective of experimental accuracies, it would be hard to obtain the actual experimental temperature distributions for the sake of CFD, and therefore all subsequent models developed in this study are initialised at 300 K and iterated for 70 cycles.
Figure 15a shows the experimental temperature contours for all twenty phases of the cycle inside the channel formed by the HHX and CHX walls, while Figure 15b shows the CFD solutions for the same “viewing area” as used in the experiment. Qualitatively, there seems to be a fairly good agreement between the two in terms of flow physics. The differences in magnitude, particularly at the right side of the viewing area, are mainly due to the difficulty of achieving numerically the experimentally recorded mean temperature of the flow. However the trends of the temperature distributions are the same for all phases in the cycle. Hot fluid starts flowing into the channel from the left during the first part of the cycle ( ϕ 1– ϕ 10). As the flow reverses ( ϕ 11– ϕ 20), the cold fluid starts flowing from the right. The temperature within the channel is bounded by the experimental temperature profile set at the cold plate (cf. Figure 3). Therefore achieving higher mean temperature within the cold channel, as observed in the experiment, cannot be obtained under the current numerical time integration and modelling. Of course the discrepancies between CFD and the experiments may come from a variety of sources including the already discussed experimental errors as well as incorrect assumptions behind the numerical model compared to reality, e.g., the already discussed adiabatic boundary conditions, heat accumulation (vs. leakage) problems, or unknown (and experimentally unverifiable) three-dimensional effects.

3.2. The Effect of the Heat Exchanger Wall Temperature on the Flow Field

The presence of the temperature field introduces new elements to the flow physics including natural convection, stronger forced convection due to gas thermal expansion, or thermo-viscous interplay in the boundary layers. A few different sets of thermal boundary conditions were applied to heat exchangers, as summarised in Table 1. Here, the temperature profile was assumed to be “flat” as only one case had experimental data available. However, this particular case ( T H = 200 °C, T C = 30 °C) was solved in two ways: first, with the wall temperature from Figure 3, for code validation, and second, with the “flat” temperature profiles to compare between different numerical cases. The drive ratio (defined as the ratio of pressure at the antinode to the mean pressure) of the flow is maintained at 0.3%.

3.2.1. Study of Flow Using the Adiabatic Model

The investigation of flow across the parallel plate structure starts with a model where the walls of the heat exchangers are treated as adiabatic. The resulting velocity profiles plotted at a location of 1 mm from the joint above the cold plate are shown in Figure 16 for comparison to the experimental results of Shi et al. [11]. The agreement with the experimental data is very good. Figure 17 shows the contours of vorticity, ω z , plotted for selected phases of the oscillatory flow, obtained as:
ω z = v x u y ,
where u and v are the velocity components in the x and y directions of the flow. The contours are plotted within the same “viewing area” as in the experiment [11]—cf. also Figure 2d. For phases ϕ 1 to ϕ 10, the fluid is flowing into the channel from the left; for phases ϕ 11 to ϕ 20 from the right. At the beginning of the cycle, a pair of small vortices appear at the entry to the channel followed by a pair of stronger counter-rotating vortices, which appear to “leap-frog” the original pair around phase ϕ 4. However, both vortex structures dissipate by about phase ϕ 8. A similar feature is present at the other end of the channel when the flow changes direction. The laminar model captures the features of the oscillatory flow at 0.3% drive ratio very well.

3.2.2. The Effect of Wall Temperature on the Flow and Heat Transfer

The effect of temperature on the flow and heat transfer is investigated by setting the heat exchanger wall temperatures to the values shown in Table 1. It should be noted that all comparisons to experiments were made using the model that applies the experimental temperature profiles of Figure 3 on the hot and cold plates (of course, except cases marked as “adiabatic” or “dT = 0”). However, for the sake of comparisons between numerical cases, “flat” temperature profiles were set up, as already mentioned in the introduction to Section 3.2. Figure 18 shows the comparison between the experimental and numerical velocity profiles for the flow at the location of 10 mm from the joint above the cold plate for comparable temperature gradients set up along the plates. The slight differences in the magnitude between the experimental and CFD data can be accounted for by the typical errors of particle image velocimetry (PIV), although it is worth mentioning that thermophoresis which affects seeding particles in the presence of temperature gradients [27] can also have an additional impact on the experimental data.
As shown in Figure 9a, buoyancy effects are observed in the experimental results. The heat supplied by the HHX causes reduced density of the fluid, which results in a buoyancy driven flow in the open area next to the HHX superimposed on the forced convection due to the flow being ejected from the channel by acoustic excitation as well as thermal expansion of the fluid travelling in the hot part of the channel. These temperature effects are likely contributors to the flow asymmetries observed in Figure 7b. These effects are also seen through a comparison of velocity profiles between the models with and without the temperature gradient imposed on the heat exchanger walls (cf., Figure 16 vs. Figure 18).
The above asymmetry can be observed in two ways: Firstly, the whole structure of the velocity profile plotted in Figure 18 seems to be “leaning” to the left in comparison to Figure 16. This is clearly illustrated in Figure 19, where the presence of temperature makes the velocity profile of ϕ 9 and its counterpart in reversed flow, ϕ 19, “lean” to the left.
Secondly, the maximum magnitude of velocity fluctuation extends further to the centreline when temperature effects are included in the model as seen in Figure 20a. In the absence of temperature effects (adiabatic wall condition shown in Figure 20 as dT = 0), the velocity profile is similar at both CHX and HHX. When the temperature effect is considered, the velocity profile at the HHX seems shifted up towards the centreline ( y = 3 mm) and the magnitude of the velocity is larger compared to an equidistant location at the CHX. On the other hand, focusing on one location, Figure 20b shows that the velocity profiles at the selected point of the CHX in analogous phases of the two flow directions (i.e., ϕ 9 vs. ϕ 19 here), are also affected by the presence of temperature effects. For ease of comparison, the axial velocity of ϕ 19 in Figure 20b is presented as a negative (−ve) value. The presence of temperature effects results in a lower centreline velocity at ϕ 9 but bigger magnitudes in the reversed flow direction at ϕ 19. The results indicate the changes in the viscous boundary layer due to the presence of temperature.
The asymmetry of the flow structure can be further illustrated through the vorticity contours at the open area next to the heat exchanger for three different wall temperature conditions (Figure 21, Figure 22 and Figure 23). The vorticity is plotted for two phases corresponding to the opposite flow direction ( ϕ 1 and ϕ 11). For each phase, the contour is shown with two views of the two ends of the heat exchanger assembly.
Figure 21 illustrates the case of adiabatic conditions and it is clear that the flow around two ends of the heat exchanger assembly is symmetrical, without any noticeable buoyancy driven convective patterns. Figure 22 shows the structure of the vortices when there is a temperature difference imposed: 200 °C at the hot plate, and 30 °C at the cold plate. The symmetry of the flow is broken: as the viscosity of the gas increases with temperature, the thicknesses of the viscous penetration depth increases, as seen, for example, in different separations of the vorticity contours on the right and left ends of the channel, which is also congruent with the velocity profiles in Figure 19a. Consequently, the vortex structures rolling up on the two sides of the plate assembly differ in size and strength. However, in addition, the temperature-driven flow at the open area at the left introduces some differences in the vortex patterns between the top and bottom side of the plate—note the hot plume observed in Figure 9 and Figure 13 which adds a “cross-flow” (vertically upward) component relative to the hot fluid ejected from the channels. Interestingly Figure 22a, for phase ϕ 11, shows the artefacts of a vortex structure on the top of the hot plate, but not on the bottom. These flow features are strengthened in Figure 23, for temperature conditions of 300 °C at the hot, and 30 °C at the cold plates.
There are two more effects that could be speculated on: First, that the large vortex structures created on the left side of hot plates may provide additional resistance to the flow when it reverses and starts to move from left to right by “bottlenecking” the channel. The same effect may be at play on the right side of the CHX when the flow reverses and starts moving to the left, but such flow would only need to handle much smaller vortices. Second, the temperature driven flow in the open area to the left of the hot plates may lead to setting up a variety of convective “cells” (loops), which could be contained on one side of the plate assembly, but equally induce a local streaming current between the plates (some of these being also adiabatic plates above and below the heated/cooled channels). Both of these effects could be further contributors to breaking the flow symmetry in the acoustic cycle.
Figure 24 summarizes all studies with the varied temperature difference between the hot and cold plates (cf., Table 1) in terms of the centreline velocity amplitude along the length of the heat exchanger channel. Two selected phases, ϕ 9 and ϕ 19, representing the first half (flow to the right) and the second half of the cycle (flow to the left) are chosen, respectively; the legend with thermal conditions is shown at the top. The variation of the plots with varying thermal conditions is indicative of a very complex physics with a number of effects at play. Clearly, the variation of the channel wall temperature and corresponding gas expansion (or contraction) causes gas particles to travel to a different distance back and forth during its reciprocating movement, causing asymmetrical velocity amplitudes on the centreline. However, the exact shape of the velocity profile depends on thermal (and thus viscous) conditions across the oscillatory boundary layers causing velocity “overshoots” and “deficits”. These may correspond to the counterintuitive reduction in the centreline velocity for the cold channel in phase ϕ 9 when the case ( T H = 100 °C, T C = 30 °C) is replaced by cases ( T H = 200 °C, T C = 30 °C) and ( T H = 300 °C, T H = 30 °C) and an increase for case ( T H = 300 °C, T C = 100 °C). Similar counterintuitive phenomena can be seen in other places of Figure 24 for both selected phases. The already mentioned convection “cells” appearing across the overall domain and the varying size and structure of vortices rolling up on the plate ends (and subsequently sucked back into the channels) may also need to be considered as contributors to the appearance of plots in Figure 24.

3.3. The Effect of Gravity and Device Orientation on Flow and Heat Transfer

The orientations of the resonator, relative to the gravity field, investigated in this study are schematically illustrated in Figure 25 (the three pairs of plates are for illustration only—all ten pairs are included in the model as before). The direction of the gravity field is modelled accordingly following the orientation. 0° corresponds to the horizontal layout in the figure, 90° indicates the hot plates above the cold plates, and 270° indicates the hot plates below the cold plates. The case with gravity switched off (g = 0) is also considered for reference. Experimental temperature profiles shown in Figure 3 are applied at the cold and hot plates in the models.
Figure 26 shows the effect of gravity and device orientation (tilt angle 0°, 90°, and 270°) on the velocity profile plotted 15 mm from the joint above the cold and hot plates for selected phases ϕ 9 and ϕ 19. In order to differentiate between the data for the cold and hot channels, the lines for the velocity profiles at the hot channel are assigned a “diamond” symbol at the top end of the lines. The effect of gravity can be seen by comparing cases assigned as “g, 0°” (presented by solid-lines) and “g = 0” (presented by dashed-lines). For phase ϕ 9 (Figure 26b), the comparison shows that gravity has little effect on the flow structure in the horizontal orientation. However, the effect of gravity is seen to be larger at ϕ 19. Comparison between the magnitude of the velocity between ϕ 9 and ϕ 19 for both “g, 0°” and “g = 0” shows that the temperature-driven buoyancy flow seems to cause a stronger asymmetry between the first and second parts of the flow cycle. The temperature-driven flow due to gravity seems to significantly assist the flow at ϕ 19 and slightly resist the flow at ϕ 9.
For an orientation other than horizontal, the gravity field affects the flow depending on the direction of flow and the relative locations of the HHX and CHX. At an orientation of 90° (presented by the dashed-dotted lines), the cold plate is located below the hot plate resulting in a reduced velocity magnitude during ϕ 9 of the flow cycle, because the buoyancy effects are resisting the flow. As the flow reverses (shown as ϕ 19), the buoyancy effect is helping the flow. As a result, the velocity magnitude becomes bigger than that at the first half of the cycle.
Conversely, when the device is at an orientation of 270° (presented by the dotted-lines), the temperature-driven buoyancy effect assists the flow at ϕ 9 (Figure 26b), resulting in a higher velocity amplitude and hence a greater fluid displacement between the plates. At ϕ 19 (Figure 26a), the flow reverses and the temperature-driven buoyancy effect is resisting the flow. Consequently, the velocity amplitude within the channel is lower than that in the first part of the cycle. For all device orientations, the velocity amplitude within the hot channel is always bigger than that at the cold channel. This is a logical consequence of the temperature-dependent density as previously discussed.
The effect of gravity and device orientation on the wall heat transfer is shown in Figure 27. The heat flux is averaged over the length of the heat exchanger and is calculated using Equation (8). The effect of gravity is less pronounced in the CHX compared to the HHX. This is in accordance with the effect of buoyancy driven flow, which is more pronounced at the HHX compared to the CHX. These effects can be looked at from the point of view of temperature profiles shown in Figure 28. It can be seen that the orientation changes the temperature field depending on the location of the HHX, the direction of flow, and the resulting temperature-driven buoyancy effect. When the flow is moving in the positive direction ( ϕ 9), the effect of device orientation on the temperature profiles is more pronounced at the CHX. When the flow reverses ( ϕ 19), the effect appears to be more obvious at the location of the HHX. Natural convection effect at ϕ 19 for HHX is stronger than that at ϕ 9 for CHX because the ‘hot plume’ appears next to the HHX (as seen in Figure 13). This is congruent with the weak gravity effects on the heat flux at the CHX compared to the HHX as seen in Figure 27. Of course, at an orientation of 270° (dashed line), when the HHX is located below the CHX, the temperature between both the HHX and CHX increases due to the temperature-driven buoyancy effect. At an orientation of 90°, presented by a dashed-dotted line, the temperature between the plates and hence the wall heat flux is lower because heat is driven away from the plates of the heat exchanger by the buoyancy effect.

3.4. Viscous Dissipation

Typically, it is assumed that viscous dissipation effects play a small role for low speed flows. However, in the flow situations presented in this paper, the contribution of viscous dissipation due to the existence of internal structures within the investigated domain may be significant and hence it is considered here. Dimensional analysis presented in Appendix A reveals the dimensionless parameters that need to be considered. In particular, it is found that comparison between viscous dissipation in areas of empty resonator and heat exchanger plates should take into consideration porosity, σ , and Eckert number, Ec—cf., Equation (A7).
In a viscous flow, the energy of the fluid motion (kinetic energy) is transformed into the internal energy of the fluid through the existence of viscous dissipation. Viscous dissipation for a two-dimensional model is defined as [28]:
Φ = 2 { ( u x ) 2 + ( v y ) 2 } + ( v x + u y ) 2 2 3 ( u x + v y ) 2 .
The computational domain is divided into four regions defined in Figure 29. The normalised viscous dissipation is averaged over the cycle and space, A , corresponding to each region and calculated as:
Φ * = 1 2 π A 0 A 0 2 π ( μ * Φ * ) d ϕ d A .
The term in the bracket on the right-hand-side of Equation (11) is calculated as follows:
μ * Φ * = 1 μ c ( d c u c σ ) 2 μ [ 2 { ( u x ) 2 + ( v y ) 2 } + ( v x + u y ) 2 2 3 ( u x + v y ) 2 ] .
Following the dimensional analysis, the weighted-average value is further multiplied by the dimensionless parameters identified as porosity, σ , and Eckert number, E c :
Φ = σ 2 E c Φ * .
The reference velocity is taken at 300 mm from the joint of the heat exchangers (labelled as m in Figure 2 and Figure 29) towards the right end of the figure within the region called the “cold end”. The reference location is selected to avoid the influence of temperature on the open area next to the HHX, where a hot plume is observed. A value of porosity, σ , is appropriately used when calculating dissipation for the open area and the area of the heat exchanger following the earlier discussion. For the open area, the porosity is set to one while an appropriate value of porosity (0.65 for the geometry investigated) is used for the heat exchangers.
Figure 30 shows the effect of the heat exchanger wall temperature on the dimensionless viscous dissipation. The viscous dissipation at the open areas at both ends of the heat exchangers are one order of magnitude lower than the dissipation in the area containing the plates. This is consistent with the fact that the working medium interacts with more wall surfaces within the area of the heat exchanger where viscous resistance is significant, hence causing more energy dissipation. When the HHX becomes hotter, the dissipation becomes bigger. This is a consequence of gas viscosity that becomes bigger as the temperature increases. Figure 30b has an expanded y-axis to show dissipation levels at the ends of the heat exchanger assembly. It shows that the dissipation of energy at the open area next to the hot end is bigger than that next to the cold end and that it increases with an increase of temperature of the HHX. This is likely the combined effect of the larger gas viscosity at higher temperature and the more pronounced gravity-driven flow as seen by the development of the hot plume at that end.
Figure 31 shows that viscous dissipation is lower in all areas when the device is set at a vertical orientation. Gravity determines the flow direction of the hot plume. For vertical arrangements (90° and 270°), the temperature-driven flow is in line with the direction of the flow. The temperature-driven flow can either assist or obstruct the flow depending on the device orientation and the direction of flow within a flow cycle. This may affect the magnitude of the velocity as the fluid flows back and forth within the heat exchanger, but the gradients of velocity near the wall are approximately the same, as seen in Figure 26. However, the low dissipation of the device arranged at vertical orientation indicates that the combined effect of temperature-driven flow and the mean flow may help reduce the viscous dissipation that may relate to the re-circulation of the fluid. The horizontal arrangement (0°) results in a temperature-driven flow perpendicular to the direction of the flow. Hence, re-circulation tends to be bigger in comparison to the vertical arrangement.

4. Conclusions

The two-dimensional ANSYS FLUENT CFD model has been set up to investigate the oscillatory flow across the parallel-plate heat exchanger at the low drive ratio of 0.3%. The flow, as represented by the velocity profiles, is well predicted by the model. For example, the velocity amplitudes at the centre of the heat exchanger channel have a maximum discrepancy of 1.8% compared to the experimental data. On the other hand, the numerically calculated temperatures have the maximum discrepancies of 6% compared to the experimental values. The differences noted in the temperature profiles and the resulting heat fluxes are likely to be the result of the combination of factors related to experimental uncertainties and the inability of the model to capture all the complexities of the real experimental setup. Overall, the current model under-predicts the space-time-averaged wall heat transfer by up to 30% in comparison to the experimental values.
The investigation of the effect of initial temperature on the flow and heat transfer helps explain the selection of an appropriate initial temperature that best suits the model. The current approaches of prescribing the initial temperature fields showed that the initialisation procedure influences both the magnitude of heat flux and velocity profiles within the structures. The investigation suggests that the current model is best initialised at 300 K. The model initialised with a temperature field defined more closely to the experimental conditions resulted in a change of the flow character, incongruent with the real flow behaviour. It can be speculated that the reason for this is the type of boundary conditions imposed on the model, which may not fully represent the reality of the experiment with its long time-history (hours), far beyond the capabilities of the current CFD model.
The imposed temperature boundary conditions (on the heat exchanger plates) is shown to cause asymmetry to the flow seemingly related to the nonlinearity caused by the combined effects of the temperature-dependent density and viscosity and the temperature-driven flow (buoyancy). The change in fluid density with respect to the temperature results in lower velocity amplitudes during the first part of the cycle and relatively higher magnitudes of velocity for the second part of the cycle. Here, an “asymmetrical boundary” creates asymmetries in the temperature and velocity fields. In addition, the temperature-driven buoyancy flow creates further “distortions” to the already asymmetrical features of the velocity and temperature fields. The results indicate that there are “net flows” through channels (or in other words “streaming”). These effects are well illustrated through changing the orientation of the device relative to the gravity field.
Temperature-dependent gravity effects influence the magnitudes of velocity profiles and temperature profiles between the plates depending on the direction of the flow and the locations of the CHX/HHX. Consequently, the heat fluxes also change. The presence of the imposed temperature field is shown to influence the velocity profile and gas displacement across the plates, hence breaking the symmetry of the flow. The vortex structures at the end of the plates, and the shear layer within the area bounded by the plates, change with temperature. Device orientation influences the flow and heat transfer due to temperature-driven buoyancy effects.
The viscous dissipation is also affected. The vertical arrangement is shown to provide a lower viscous dissipation but care should also be taken when dealing with this kind of arrangement because the temperature-driven flow changes the temperature and velocity field within the channel.

Acknowledgments

The first author would like to acknowledge the sponsorship of the Ministry of Education, Malaysia and Universiti Teknikal Malaysia Melaka (FRGS/1/2015/TK03/FKM/03/F00274). The second author would like to acknowledge the sponsorship of the Royal Society Industry Fellowship (2012–2015), grant number IF110094 and EPSRC Advanced Research Fellowship, grant GR/T04502.

Author Contributions

Fatimah A.Z. Mohd Saat performed the numerical investigations, analysed the data, and drafted the paper; Artur J. Jaworski planned and supervised the research work, outlined the paper, and subsequently worked with first author on consecutive “iterations” of the paper until the final manuscript was produced.

Conflicts of Interest

The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Appendix A. Dimensional Analysis

The dimensional analysis required in Section 3.4 is based on the analysis of the transport equations that govern the flow. The analysis is carried out using a two-dimensional Navier–Stokes equation for compressible flow. In this study, the normalisation parameters are introduced as:
u * = u σ u c ;   v * = v σ u c ;   x * = x d c ;   y * = y d c ;   t * = t ω ;   p * = p ρ c σ 2 u c 2 ; ρ * = ρ ρ c ;   g * = g g c ;   β * = β β c ;   μ * = μ μ c ;   T * = T T c T H T c ;   c p * = c p c p , c ;   k * = k k c .
The superscript * describes the dimensionless variable. The subscript c refers to a reference value, preferably selected in the open area a distance away from the stack/regenerator where the properties and behaviour of the fluid are not affected much by the presence of the object. The velocity is normalised by the velocity at the reference point multiplied by the porosity, σ . The dimensions x and y are both normalised by the characteristic dimension, d c —here the length of the heat exchanger plate (35 mm) was used. The temperatures T H and T C selected represent the hot and cold plate/fluid temperatures, respectively.
The resulting dimensionless equations for the two-dimensional Navier–Stokes equations are:
Continuity equation
R e ω ρ * t * + σ R e d ( ρ * u * ) x * + σ R e d ( ρ * v * ) y * = 0 ,
x-momentum equation
R e ω t * ( ρ * u * ) + σ R e d x * ( ρ * u * 2 ) + σ R e d y * ( ρ * u * v * ) + σ R e d p * x * = G r σ R e d ρ * g * β * T * + x * [ 2 μ * u * x * 2 3 μ * ( u * x * + v * y * ) ] + y * [ μ * ( v * x * + u * y * ) ] ,
y-momentum equation
R e ω t * ( ρ * v * ) + σ R e d x * ( ρ * u * v * ) + σ R e d y * ( ρ * v * 2 ) + σ R e d p * y * = G r σ R e d ρ * g * β * T * + y * [ 2 μ * v * y * 2 3 μ * ( u * x * + v * y * ) ] + x * [ μ * ( v * x * + u * y * ) ] ,
Energy equation
R e ω θ t * ( ρ * c p * T * ) + σ R e d θ [ x * ( ρ * c p * u * T * ) + y * ( ρ * c p * v * T * ) ] = θ P r [ x * ( k * T * x * ) + y * ( k * T * y * ) ] + σ 2 E c [ R e ω p * t + σ R e d ( x * ( p * u * ) + y * ( p * v * ) ) ] + σ 2 E c ( μ * Φ * ) .
The term Φ * is defined as:
Φ * = 2 { ( u * x * ) 2 + ( v * y * ) 2 } + ( v * x * + u * y * ) 2 2 3 ( u * x * + v * y * ) 2 ,
The dimensionless parameters discovered through this procedure are:
R e ω = ω d c 2 v c   ;   R e d = u c d c 2 v c   ;   G r = g c β c ( T H T c ) d c 3 υ c 2   ;   θ = T H T C T c   ; E c = u c 2 c p , c T c   ;   P r = μ c c p , c k c
The kinetic Reynolds number, R e ω , describes the influence of frequency on the flow. The hydraulic Reynolds number, R e d , describes the flow based on the amplitude of the velocity, u c . The combined effect of frequency and flow amplitude is related to a dimensionless number known as the Keulegen-Carpenter number, K C = R e d / R e ω . The temperature-driven characteristic is represented by the Grashof number, G r . The dimensionless temperature, θ , represents the effect of temperature on the flow. Depending on the definition of the reference temperature used, this dimensionless number may represent the effect of heat accumulation within the investigated area. The effect can be seen in the transient and convection part of the energy equation. This indicates that heat accumulation changes over time and is very much dependent on the amplitude of the flow. The Eckert number, E c , expresses the relationship between the flow kinetic energy and enthalpy, and is used to characterise dissipation. The Prandtl number, P r , describes the property of fluids, which is useful when comparing cases with different fluid media.

References

  1. Swift, G.W. Thermoacoustics: A Unifying Perspective for Some Engines and Refrigerators; Acoustical Society of America Publications: Sewickley, PA, USA, 2002. [Google Scholar]
  2. Mao, X.; Jaworski, A.J. Oscillatory flow at the end of parallel-plate stacks: Phenomenological and similarity analysis. Fluid Dyn. Res. 2010, 42, 055504. [Google Scholar] [CrossRef]
  3. Jaworski, A.J.; Mao, X.; Mao, X.; Yu, Z. Entrance effects in the channels of the parallel plate stack in oscillatory flow conditions. Exp. Therm. Fluid Sci. 2009, 33, 495–502. [Google Scholar] [CrossRef]
  4. Matveev, K.; Backhaus, S.; Swift, G. On Some Nonlinear Effects of Heat Transport in Thermal Buffer Tubes. In Proceedings of the 17th International Symposium on Nonlinear Acoustics, Pennsylvania, PA, USA, 18–22 July 2006. [Google Scholar]
  5. Mao, X.; Yu, Z.; Jaworski, A.J.; Marx, D. PIV studies of coherent structures generated at the end of a stack of parallel plates in a standing wave acoustic field. Exp. Fluids 2008, 45, 833–846. [Google Scholar] [CrossRef]
  6. Worlikar, A.S.; Knio, O.M. Numerical simulation of a thermoacoustic refrigerator I. Unsteady adiabatic flow around the stack. J. Comput. Phys. 1996, 127, 424–451. [Google Scholar] [CrossRef]
  7. Worlikar, A.S.; Knio, O.M.; Klein, R. Numerical simulation of thermoacoustic refrigerator II. Stratified flow around the stack. J. Comput. Phys. 1998, 144, 299–324. [Google Scholar] [CrossRef]
  8. Besnoin, E.; Knio, O.M. Numerical Study of thermoacoustic heat exchangers in the thin plate limit. Numer. Heat Transfer. Part A 2001, 40, 445–471. [Google Scholar]
  9. Besnoin, E.; Knio, O.M. Numerical study of thermoacoustic heat exchangers. Acta Acust. United Acust. 2004, 90, 432–444. [Google Scholar]
  10. Mozurkewich, G. Heat transfer from a cylinder in an acoustic standing wave. J. Acoust. Soc. Am. 1995, 98, 2209–2216. [Google Scholar] [CrossRef]
  11. Shi, L.; Yu, Z.; Jaworski, A.J. Application of laser-based instrumentation for measurement of time-resolved temperature and velocity fields in the thermoacoustic system. Int. J. Therm. Sci. 2010, 49, 1688–1701. [Google Scholar] [CrossRef]
  12. ANSYS Fluent; 13.0; User Manual; ANSYS Inc.: Canonsburg, PA, USA, 2010.
  13. Guillaume, D.W.; LaRue, J.C. Comparison of the vortex shedding behaviour of a single plate and plate array. Exp. Fluids 2001, 30, 22–26. [Google Scholar] [CrossRef]
  14. Guillaume, D.W.; LaRue, J.C. Comparison of the numerical and experimental flow field downstream of a plate array. J. Fluids Eng. 2002, 124, 184–186. [Google Scholar] [CrossRef]
  15. Parang, M.; Salah-Eddine, A. Thermoacoustic Convection Heat-Transfer Phenomenon. AIAA 1984, 22, 1020–1022. [Google Scholar]
  16. Lin, Y.; Farouk, B.; Oran, E.S. Interactions of thermally induced acoustic waves with buoyancy induced flows in rectangular enclosures. Int. J. Heat Mass Transf. 2008, 51, 1665–1674. [Google Scholar] [CrossRef]
  17. Shi, L.; Mao, X.; Jaworski, A.J. Application of planar laser-induced fluorescence measurement techniques to study the heat transfer characteristics of parallel-plate heat exchangers in thermoacoustic devices. Meas. Sci. Technol. 2010, 21, 115405. [Google Scholar] [CrossRef]
  18. Kuzuu, K.; Hasegawa, S. Numerical investigation of heated gas flow in a thermoacoustic device. Appl. Therm. Eng. 2017, 110, 1283–1293. [Google Scholar] [CrossRef]
  19. Shen, C.; He, Y.; Li, Y.; Ke, H.; Zhang, D.; Liu, Y. Performance of solar powered thermoacoustic engine at different tilted angles. Appl. Therm. Eng. 2009, 29, 2745–2756. [Google Scholar] [CrossRef]
  20. Yassen, N. Impact of temperature gradient on thermoacoustic refrigerator. Energy Procedia 2015, 74, 1182–1191. [Google Scholar] [CrossRef]
  21. Abramenko, T.N.; Aleinikova, V.I.; Golovicher, L.E.; Kuz’mina, N.E. Generalization of experimental data on thermal conductivity of nitrogen, oxygen, and air at atmospheric pressure. J. Eng. Phys. Thermophys. 1992, 63, 892–897. [Google Scholar] [CrossRef]
  22. Tasnim, S.H.; Fraser, R.A. Computation of the flow and thermal fields in a thermoacoustic refrigerator. Int. Commun. Heat Mass Transf. 2010, 37, 748–755. [Google Scholar] [CrossRef]
  23. Cha, J.S.; Ghiaasiaan, S.M. Oscillatory flow in microporous media applied in pulse-tube and Stirling-cycle cryocooler regenerators. Exp. Therm. Fluid Sci. 2008, 32, 1264–1278. [Google Scholar] [CrossRef]
  24. Frederix, E.M.A.; Stanic, M.; Kuczaj, A.K.; Nordlund, M.; Geurts, B.J. Extension of the compressible PISO algorithm to single-species aerosol formation and transport. Int. J. Multiph. Flow 2015, 74, 184–194. [Google Scholar] [CrossRef]
  25. Versteeg, H.K.; Malalasekera, W. An Introduction to Computational Fluid Dynamics the Finite Volume Method, 2nd ed.; Prentice Hall: Essex, UK, 2007; pp. 179–197, 263–265. [Google Scholar]
  26. Piccolo, A. Numerical computation for parallel plate thermoacoustic heat exchangers in standing wave oscillatory flow. Int. J. Heat Mass Transf. 2011, 54, 4518–4530. [Google Scholar] [CrossRef]
  27. Talbot, L.; Cheng, R.K.; Schefer, R.W.; Willis, D.R. Thermophoresis of particles in a heated boundary layer. J. Fluid Mech. 1980, 101, 737–758. [Google Scholar] [CrossRef]
  28. Schlichting, H.; Gersten, K. Boundary Layer Theory, 8th revised and enlarged edition; Springer: New Delhi, India, 2001; pp. 76, 416–507. [Google Scholar]
Figure 1. Schematic of a simple thermoacoustic cooler arrangement (top). The acoustically induced compression and expansion of fluid elements causes heat pumping effects along the stack (bottom).
Figure 1. Schematic of a simple thermoacoustic cooler arrangement (top). The acoustically induced compression and expansion of fluid elements causes heat pumping effects along the stack (bottom).
Applsci 07 00417 g001
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.
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.
Applsci 07 00417 g002
Figure 3. Cold (CHX) and hot (HHX) heat exchanger wall temperature.
Figure 3. Cold (CHX) and hot (HHX) heat exchanger wall temperature.
Applsci 07 00417 g003
Figure 4. Oscillating pressure and velocity monitored for identification of the steady oscillatory condition.
Figure 4. Oscillating pressure and velocity monitored for identification of the steady oscillatory condition.
Applsci 07 00417 g004
Figure 5. Pressure, velocity, and displacement over 20 phases of a flow cycle.
Figure 5. Pressure, velocity, and displacement over 20 phases of a flow cycle.
Applsci 07 00417 g005
Figure 6. Grid independency test for (a) axial velocity amplitude at the centre of the channel for phases ϕ 6 and ϕ 16 and (b) velocity profiles near the wall; both taken 15 mm from the joint above the cold plate.
Figure 6. Grid independency test for (a) axial velocity amplitude at the centre of the channel for phases ϕ 6 and ϕ 16 and (b) velocity profiles near the wall; both taken 15 mm from the joint above the cold plate.
Applsci 07 00417 g006
Figure 7. Centreline velocity amplitude (1 mm from the joint, above the cold plate) obtained from the model with heat exchanger walls: (a) adiabatic and (b) with the temperature profile based on the experiment.
Figure 7. Centreline velocity amplitude (1 mm from the joint, above the cold plate) obtained from the model with heat exchanger walls: (a) adiabatic and (b) with the temperature profile based on the experiment.
Applsci 07 00417 g007
Figure 8. Temperature profiles for HHX and CHX—comparison between experimental (EXP) and computational (CFD) results at the location 15 mm from the joint (above the HHX and CHX plates).
Figure 8. Temperature profiles for HHX and CHX—comparison between experimental (EXP) and computational (CFD) results at the location 15 mm from the joint (above the HHX and CHX plates).
Applsci 07 00417 g008
Figure 9. Temperature contours at phase ϕ 1 from the experiment (a) and three numerical models initialised with different temperature fields (bd).
Figure 9. Temperature contours at phase ϕ 1 from the experiment (a) and three numerical models initialised with different temperature fields (bd).
Applsci 07 00417 g009
Figure 10. The effect of initial temperature on (a) the oscillating velocity at location 1 mm from the joint above the cold plate and (b) the velocity amplitude along the heat exchanger’s plate taken at the centre of the channel for ϕ 6 and ϕ 16.
Figure 10. The effect of initial temperature on (a) the oscillating velocity at location 1 mm from the joint above the cold plate and (b) the velocity amplitude along the heat exchanger’s plate taken at the centre of the channel for ϕ 6 and ϕ 16.
Applsci 07 00417 g010
Figure 11. The effect of initial temperature on the wall heat flux of the heat exchangers.
Figure 11. The effect of initial temperature on the wall heat flux of the heat exchangers.
Applsci 07 00417 g011
Figure 12. Temperature profiles at the open area next to the (a) cold and (b) hot ends of the heat exchanger assembly—comparison between the experimental and numerical model with different approaches of temperature initialisation.
Figure 12. Temperature profiles at the open area next to the (a) cold and (b) hot ends of the heat exchanger assembly—comparison between the experimental and numerical model with different approaches of temperature initialisation.
Applsci 07 00417 g012
Figure 13. Illustration of the convective currents occurring at the top area of the heat exchanger assembly.
Figure 13. Illustration of the convective currents occurring at the top area of the heat exchanger assembly.
Applsci 07 00417 g013
Figure 14. The effect of initial temperature on the temperature profiles between the heat exchanger plates at a location 10 mm from the joint above the cold (CHX), and hot (HHX) plates.
Figure 14. The effect of initial temperature on the temperature profiles between the heat exchanger plates at a location 10 mm from the joint above the cold (CHX), and hot (HHX) plates.
Applsci 07 00417 g014
Figure 15. Temperature contours in the area bounded by the heat exchanger walls—comparison between (a) experiment of Shi et al. [17] and (b) simulation.
Figure 15. Temperature contours in the area bounded by the heat exchanger walls—comparison between (a) experiment of Shi et al. [17] and (b) simulation.
Applsci 07 00417 g015
Figure 16. Velocity profiles from computational fluid dynamics (CFD) model (black and grey lines) and the experiment (red and black symbol) for all 20 phases of a flow cycle. The heat exchanger walls are adiabatic.
Figure 16. Velocity profiles from computational fluid dynamics (CFD) model (black and grey lines) and the experiment (red and black symbol) for all 20 phases of a flow cycle. The heat exchanger walls are adiabatic.
Applsci 07 00417 g016
Figure 17. Vorticity contours within the channel with plates treated as adiabatic walls—comparison between the results from (a) the experiment of Shi et al. [11] and (b) CFD.
Figure 17. Vorticity contours within the channel with plates treated as adiabatic walls—comparison between the results from (a) the experiment of Shi et al. [11] and (b) CFD.
Applsci 07 00417 g017
Figure 18. Velocity profiles from CFD (black and grey lines) and the experiment (red and black symbol) for 20 phases of a flow cycle with the influence from heat exchanger wall temperature.
Figure 18. Velocity profiles from CFD (black and grey lines) and the experiment (red and black symbol) for 20 phases of a flow cycle with the influence from heat exchanger wall temperature.
Applsci 07 00417 g018
Figure 19. The effect of the imposed temperature gradient (dT) on the velocity profile taken at 15 mm from the joint above the cold plate (CHX).
Figure 19. The effect of the imposed temperature gradient (dT) on the velocity profile taken at 15 mm from the joint above the cold plate (CHX).
Applsci 07 00417 g019
Figure 20. Velocity profile; (a) for ϕ 9 at CHX and HHX, both at a location 15 mm away from the joint; (b) at CHX only at the different flow direction.
Figure 20. Velocity profile; (a) for ϕ 9 at CHX and HHX, both at a location 15 mm away from the joint; (b) at CHX only at the different flow direction.
Applsci 07 00417 g020
Figure 21. Vorticity contour for the case of the heat exchanger set as the adiabatic wall at the (a) hot and (b) cold plates. (1 bar in nitrogen gas, 0.3% drive ratio, 13.1 Hz).
Figure 21. Vorticity contour for the case of the heat exchanger set as the adiabatic wall at the (a) hot and (b) cold plates. (1 bar in nitrogen gas, 0.3% drive ratio, 13.1 Hz).
Applsci 07 00417 g021
Figure 22. Vorticity contour for the case of the heat exchanger set at 200 °C at the (a) hot and 30 °C at the (b) cold plates. (1 bar in nitrogen gas, 0.3% drive ratio, 13.1 Hz, temperature profiles assumed “flat”).
Figure 22. Vorticity contour for the case of the heat exchanger set at 200 °C at the (a) hot and 30 °C at the (b) cold plates. (1 bar in nitrogen gas, 0.3% drive ratio, 13.1 Hz, temperature profiles assumed “flat”).
Applsci 07 00417 g022
Figure 23. Vorticity contour for the case with the heat exchanger set at 300 °C at the (a) hot and 30 °C at the (b) cold plates. (1 bar in nitrogen gas, 0.3% drive ratio, 13.1 Hz, temperature profiles assumed “flat”).
Figure 23. Vorticity contour for the case with the heat exchanger set at 300 °C at the (a) hot and 30 °C at the (b) cold plates. (1 bar in nitrogen gas, 0.3% drive ratio, 13.1 Hz, temperature profiles assumed “flat”).
Applsci 07 00417 g023
Figure 24. Axial velocity at the middle of the channel along the hot (HHX) and cold (CHX) heat exchangers; temperature profiles assumed “flat”.
Figure 24. Axial velocity at the middle of the channel along the hot (HHX) and cold (CHX) heat exchangers; temperature profiles assumed “flat”.
Applsci 07 00417 g024
Figure 25. Illustration of the thermoacoustic device orientation. Red and blue plates symbolise hot and cold heat exchangers, respectively.
Figure 25. Illustration of the thermoacoustic device orientation. Red and blue plates symbolise hot and cold heat exchangers, respectively.
Applsci 07 00417 g025
Figure 26. The effect of gravity and device orientation on the velocity profiles at a location 15 mm from the joint above the cold (CHX) and hot (CHX) plates for phases (a) ϕ 19 and (b) ϕ 9.
Figure 26. The effect of gravity and device orientation on the velocity profiles at a location 15 mm from the joint above the cold (CHX) and hot (CHX) plates for phases (a) ϕ 19 and (b) ϕ 9.
Applsci 07 00417 g026
Figure 27. The effect of gravity and device orientation on the wall heat flux of the hot (HHX) and cold (CHX) plates.
Figure 27. The effect of gravity and device orientation on the wall heat flux of the hot (HHX) and cold (CHX) plates.
Applsci 07 00417 g027
Figure 28. Temperature profiles for two selected phases, when there is no gravity (g = 0) and with gravity, g, at different tilt angles; 0°, 90°, 270°, at a location 15 mm from the joint above the cold (CHX) and hot (HHX) plates.
Figure 28. Temperature profiles for two selected phases, when there is no gravity (g = 0) and with gravity, g, at different tilt angles; 0°, 90°, 270°, at a location 15 mm from the joint above the cold (CHX) and hot (HHX) plates.
Applsci 07 00417 g028
Figure 29. Illustration of the region defined for the analysis of viscous dissipation.
Figure 29. Illustration of the region defined for the analysis of viscous dissipation.
Applsci 07 00417 g029
Figure 30. (a) The effect of the heat exchanger wall temperature on viscous dissipation and (b) the enlarged view for dissipations at the open area next to the heat exchangers.
Figure 30. (a) The effect of the heat exchanger wall temperature on viscous dissipation and (b) the enlarged view for dissipations at the open area next to the heat exchangers.
Applsci 07 00417 g030
Figure 31. (a) The effect of the device’s orientation on viscous dissipation and (b) the enlarged view for viscous dissipation at the open area next to the heat exchangers.
Figure 31. (a) The effect of the device’s orientation on viscous dissipation and (b) the enlarged view for viscous dissipation at the open area next to the heat exchangers.
Applsci 07 00417 g031
Table 1. Wall temperature condition of the heat exchanger.
Table 1. Wall temperature condition of the heat exchanger.
Drive Ratio (%)Working MediumFrequency (Hz)Mean Pressure (bar)Heat Exchanger Wall Temperature (°C)
Hot, T H Cold, T C
0.3Nitrogen13.11AdiabaticAdiabatic
20030
10030
30030
300100

Share and Cite

MDPI and ACS Style

Mohd Saat, F.A.Z.; Jaworski, A.J. The Effect of Temperature Field on Low Amplitude Oscillatory Flow within a Parallel-Plate Heat Exchanger in a Standing Wave Thermoacoustic System. Appl. Sci. 2017, 7, 417. https://0-doi-org.brum.beds.ac.uk/10.3390/app7040417

AMA Style

Mohd Saat FAZ, Jaworski AJ. The Effect of Temperature Field on Low Amplitude Oscillatory Flow within a Parallel-Plate Heat Exchanger in a Standing Wave Thermoacoustic System. Applied Sciences. 2017; 7(4):417. https://0-doi-org.brum.beds.ac.uk/10.3390/app7040417

Chicago/Turabian Style

Mohd Saat, Fatimah A.Z., and Artur J. Jaworski. 2017. "The Effect of Temperature Field on Low Amplitude Oscillatory Flow within a Parallel-Plate Heat Exchanger in a Standing Wave Thermoacoustic System" Applied Sciences 7, no. 4: 417. https://0-doi-org.brum.beds.ac.uk/10.3390/app7040417

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop