- Email: [email protected]
A 3D investigation of thermoacoustic fields in a square stack
International Journal of Heat and Mass Transfer 108 (2017) 292–300
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A 3D investigation of thermoacoustic fields in a square stack Ahmed I. Abd El-Rahman ⇑, Waleed A. Abdelfattah, Mahmoud A. Fouad Department of Mechanical Power Engineering, Cairo University, Cairo University Rd., Giza 12613, Egypt
a r t i c l e
i n f o
Article history: Received 10 April 2016 Received in revised form 2 December 2016 Accepted 5 December 2016
Keywords: Thermoacoustic refrigerator Square stack 3D simulation Energy flux density Flow streaming
a b s t r a c t The thermoviscous behavior of the oscillating gas within the porous medium of a thermoacoustic refrigerator enables the conversion of sound into heat in the process of typical standing-wave thermoacoustic refrigeration systems. Several nonlinear mechanisms, such as harmonics generation, self-induced streaming and possible turbulence, cause complicated flow behavior and influence the performance of such devices. Few analytical and numerical approximations (Cao et al., 1996; Worlikar and Knio, 1996; Worlikar et al., 1998) describe the flow field and the energy flux density in standing devices comprising stacks of parallel plates, but almost no 3-D simulation has been developed that models the largeamplitude excitations and enables the prediction of the oscillating-flow behavior within a porous medium made of square channels. Here, we build on existing effort (Abd El-Rahman and Abdel-Rahman, 2014) and report a computational fluid dynamics analysis of a Helium-filled half-wavelength thermoacoustic refrigerator. The finite volume method is used, and the solid and gas domains are represented by large numbers of hexahedral elements. The calculations assume a 3-D periodic cell structure to reduce the computational cost and apply the dynamic mesh technique to account for the adiabatically oscillating wall boundaries. The simulation uses an implicit time integration of the full unsteady compressible flow equations along with the conjugate heat transfer algorithm. For the sake of simulation, the geometry and operating conditions closely follow the experimental setup of Lotton et al. (2009) [5]. A typical run involves about 200,000 computational cells, whereas the working frequency and the applied drive ratio are 592.5 Hz and 2.0%, respectively. The stationary system response is predicted and the numerical values are compared to theory. The results show good agreement with theoretical behavior in the linear regime of drive ratios up to 2.0% while a maximum cooling effect of 10 degrees is captured between the stack ends. This simulation provides an interesting tool for understanding the transient flow behavior, mean energy-flux density and flow streaming, and helps building 3-D computational fluid dynamics models for thermoacoustic devices. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction A typical standing-wave thermoacoustic refrigerator (TAR) consists of a long resonator that is usually filled with an inert gas, a porous medium called ‘Stack’ and two heat exchangers. The stack of a standing-wave refrigerator is usually short and made of ceramic material, whereas the resonator is closed at one end and bounded by a loudspeaker at the other end, thus, allowing the propagation of a half-wavelength acoustic excitation. The hot and cold heat exchangers are made of copper to provide efficient heat transfer between the working inert gas and the external heat source and sink respectively. TARs are interesting because they have no moving parts, unlike conventional refrigerators, and almost no environmental impact exists as they rely on clean ⇑ Corresponding author. E-mail address: [email protected] (A.I. Abd El-Rahman). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.12.015 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
conversion of acoustic and heat energies. Their fabrication process is rather simpler and sizes span wide variety of length scales. In typical low-pressure excitation systems, the thermoacoustic behavior is well described by the 1-D linear theory. However, practical systems fall into the high-pressure regime of acoustic excitation, where non-linear behavior is dominant. The following recent studies showed substantial deviation of most cooling curves from the linear theory. They rather focused on modeling thermoacoustic devices to help predict the general performance, capture the thermal and viscous flow losses in the device, and define key variables affecting the design considerations and system optimization, such as stack porosity, size and location, working gas and appropriate operating mean pressure. Cao et al. [1] was the first to construct a first-order temporal and spatial discretization scheme to solve for the time-averaged energy flux density in the gas. In their model, the stack was represented by short (compared with the acoustic wavelength)
293
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
isothermal parallel plates. Cao’s attempt was followed by the significant effort of Worlikar et al. [2,3] who considered the lowMach number approximation of the full Navier-Stokes equations and solved the streamfunction-vorticity equation. They were able to successfully reduce the full computational domain into a single periodic thermoacoustic couple. Although limited to twodimensional analysis, short stack lengths and constant gas properties, their approach was then traced in similar fashion by other recent workers, such as Blanc-Benon et al. [6]. Blanc-Benon performed a particle image velocimetry (PIV) measurements of the flow field and reported close agreement with the low-Mach number simulation results. Ishikawa and Mee [7] extended Cao analysis to longer plates and investigated the effects of plate spacing using the PHOENICS finite-volume method, however, their single-precision calculations caused significant energy imbalances and restricted the overall accuracy in their heat fluxes results. Latter, Marx and BlancBenon [8,9] considered the time-averaged velocity field above the stack plate and observed the acoustic flow streaming above the plate. They further identified the minimum plate length, as compared to the particle displacement, necessary for the generated temperature harmonics to get displayed above the entire plate surface. Although accurate, this model is too expensive computationally to limit its applicability to devices operating at high-frequency of 20 kHz, which is useful for miniaturization purposes. Zoontjens et al. [10] was the first to implement a FLUENT model to simulate the flow and energy fields in a TAR couple of non-zero thickness at a wide range of drive ratio (1.7–8.5%) and 200 Hz operating frequency. They assumed symmetric computational domain and investigated the time-averaged heat transfer through the stack material. Their model however used a first-order discretization scheme which in turn affected the accuracy of their numerical results. This was recently followed by the work of Tasnim and Fraser [11] who introduced a two-dimensional model using the STARCD finite volume solver. He assumed a periodic thermoacoustic couple and described the laminar-flow and thermal fields at low drive ratio of 0.7%. To our knowledge, no 3-D simulation has been developed that directly models the large-amplitude excitations of standing-wave TARs. Building on existing effort [4], the objective of the present work is to extend the two-dimensional analysis of standing-wave thermoacoustic couples into three-dimensional simulation of the oscillating-flow behavior within a porous medium made of square channels. The next objective includes the investigation of the thermoviscous interactions along the stack and at the stack ends, and to compare the axial velocity profiles with the corresponding theoretical behavior at low and moderate drive ratios. Furthermore, the energy flux density is evaluated and plotted at two different longitudinal sections to illustrate the three-dimensional flow behavior in square pores.
2. Reference model 2.1. Physical model Inspired by the setup of Lotton et al. [5], the present thermoacoustic refrigerator consists of an 86-cm-long Helium-filled resonator that has a square cross section with an edge length of eight centimeters, as schematically drawn in Fig. 1. The halfwavelength resonator has a loudspeaker placed at one end and is closed at the other end. The porous media ‘stack’ is made of cordierite (ceramic material) and consists of a large number of small square channels (600 CPSI), aligned with the direction of wave propagation. Here, each channel is 10-mm-in-length and has a transversal width h ¼ H t s of 0.92 mm and a wall thickness ts
Fig. 1. Schematic of the thermoacoustic refrigerator is shown on top. Bottom: The representative element volume (REV) is drawn with a couple of periodic surfaces along with an acoustically-matching adiabatic boundary conditions imposed at both ends (the edge length of the square pore is h ¼ H t s ).
of 0.12 mm, thus resulting in a stack porosity / of 0.787. The model is constructed such that the stack center is positioned at a distance 14 cm ( k=12:3) from the resonator closed end. According to Poignand et al. [12], such stack position corresponds to the maximum cooling effect when using air as the working gas. Heat exchangers are not included in the present model setup for simplicity. The system is excited with an acoustic frequency of 592.5 Hz and acoustic pressure amplitude in the range from 250 to 2000 Pa while operating at atmospheric pressure and ambient temperature. Details about flow conditions and material thermophysical properties are remarked in Table 1. 2.2. Numerical model It is unreasonable to model the entire domain of the thermoacoustic resonator together with the square stack because of their excessively dispersed length scales. Therefore, an element volume that includes one gas channel is chosen to represent the oscillating-flow behavior within the stack and at the stack ends, as shown in Fig. 1. Such representative element volume (REV) can be thought as a unit cell that is surrounded by large number of similarly repeated cells (Or, images). To account for those image cells in both transversal and spanwise directions, the present REV is considered bounded by periodic surfaces in both directions. In the axial direction, however, the 3-D computational domain is then allowed to extend with a total length of 4.64 times the stack length l. Acoustically matching boundary conditions are thus imposed at both left and right ends of the chosen computational domain,
Table 1 Flow conditions and material thermophysical properties. Property
Value
Units
Flow Conditions Operating frequency, f Ambient Temperature, T m Mean Pressure, pm
592:5 300 101:325
Hz K kPa
0:1625
kg m3
Gas properties: Helium Density, qf Prandtl number, Pr Thermal Conductivity, kf Heat capacity, cp;f Dynamic viscosity, lf
0:68 0:152 5193
Ratio of specific heats, c
1:99 105 1:67
Stack properties: Cordierite Thermal Conductivity, ks Heat capacity, cp;s Density, qs
0:19 1000 1400
W m1 K1 J kg1 K1 kg m1 s1
W m1 K1 J kg1 K1 kg m3
294
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
assuming an ideal acoustic wave propagation. Although this approximation does not allow for accurate replication of the physical process at these boundaries, it can still suggest appropriate modeling at the far field of the stack plates provided that the distance from the two vertical surfaces to the stack ends is big enough, as argued by Cao et al. [1]. For instance, the length of the computational domain to the stack length is chosen to ensure that the flow at the domain left and right boundaries can be accurately modeled as a one-dimensional planar wave. Fig. 2 presents the domain meshing in the transversal and longitudinal directions. The domain is carefully discretized to properly capture the thermoviscous interactions of the oscillating-flow with the stack walls and to maintain grid-independent system response. The flow is assumed laminar. Such flow approximation will be elaborated in the results section through the specification of two key parameters [13,14]: the dimensionless frequency parameter
Table 2 Main meshing details, as illustrated in Fig. 2. Parameter
Value
Units
Fluid Domain - Inside the Stack Boundary cell size (MinjMax)
ð0:03Þ2 0:024jð0:03Þ2 0:12
mm3
ð0:06Þ 0:024jð0:06Þ 0:12 21 21 84
mm3
Number of cells Fluid Domain - Outside the Stack Boundary cell size (MinjMax)
ð0:03Þ2 0:025jð0:03Þ2 0:15
mm3
Centerline cell size (MinjMax)
Centerline cell size (MinjMax) Number of cells Total number of cells Solid Domain Corner cell size (MinjMax) Middle cell size (MinjMax) Total number of cells
2
2
2
2
ð0:06Þ 0:025jð0:06Þ 0:15 25 25 244 189,544
–
mm3 – –
ð0:03Þ2 0:024jð0:03Þ2 0:12
mm3
ð0:03Þ2 0:06j0:03 0:06 0:12 14,784
mm3 –
2
h x=m (Also, referred as the Womeresley number for circular pipe)
^ m =m, based on the amplitude and the flow Reynolds number Re ¼ Ud ^ and the viscous of the cross-sectional mean-velocity variation U pffiffiffiffiffiffiffiffiffiffiffiffi diffusion thickness dm ¼ 2m=x. Here, m and x are the gas viscous diffusivity and the angular frequency, respectively. In the transversal direction, the mesh is distinguished into two concentric regimes, depth of each corresponds to the thermal difpffiffiffiffiffiffiffiffiffiffiffiffiffi fusion thickness dj ¼ 2a=x, where a is the gas thermal diffusivity. In the first regime, the flow experiences steep velocity and temperature gradients, therefore, eight computational cells are enforced to capture the flow behavior very near to the wall. Less steep gradients are expected in the second regime, thus, five larger cells are assumed sufficient since this regime lies in the stack core, where poor transversal fluctuation in the flow properties generally evolves. In the longitudinal direction, a structured hexahedral mesh is maintained such that the cell height varies gradually and smoothly from one to five times the cell width throughout the entire domain, as indicated in Fig. 2(Left). For instance, the resulting total number of cells is around 2 105 . The main meshing details are further presented in Table 2. In the present simulation, the adiabatic thermal boundary condition is applied at the matching boundaries whereas the viscous interaction with the stack inner walls satisfies the no-slip boundary condition. At the solid-fluid interface, thermally coupled boundary conditions are imposed to simulate the thermal interaction between the solid stack walls and the surrounding gas
particles through the conjugate heat transfer algorithm. In such case, the solver requires that the temperature T s and the heat transfer to the wall boundary from a solid element match the gas side temperature T g and heat transfer, such that T s ¼ T g at the stack inner walls and:
Kg
@T g @T s ¼ Ks @n @n
ð1Þ
In Eq. (1), the heat transfer is expressed in terms of the gas and solid thermal conductivities K g and K s , whereas n indicates to the local coordinate normal to the wall. Next, in order to generate the acoustically matching boundary condition, the layering dynamic meshing technique is applied, as detailed in the ANSYS FLUENT user documentation [15]. A preprocessing program is developed to define the instantaneous axial motions of the matching boundaries U MB , consistent with the local ideal standing wave propagation:
U MB ¼ U 0 sinðkxÞ cosðxtÞ
ð2Þ
Here, U 0 represents the maximum amplitude of the acoustic velocity, while k and x refer to the wave number and the axial position of the matching boundaries, measured from the resonator closed end, respectively. x is the angular velocity and t indicates the simulation time. Further, U 0 can be expressed in terms of the pressure amplitude at the pressure antinodes of the standing wave (i.e. resonator
Fig. 2. The left figure presents the 3-D meshing of the computational gas domain, as spanned from inside the stack to the outside of it. The solid mesh is suppressed for clarity. However, the right plot shows a side view -transversal cut- of the whole computational domain (Inset: dj 0:31 mm). Here, the yellow frame represents the solid stack, while the gray structured mesh corresponds to the gas domain with a minimum cell size of 0.03 mm and a maximum of 0.06 mm. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
closed end pL ), the mean density q, and the speed of sound a, such that U 0 ¼ pL =qa. The present analysis uses the ANSYS FLUENT second-order upwind implicit-time algorithm pressure-based CFD solver with the conjugate heat transfer algorithm along with a second-order formulation for the temporal discretization. The CFD solver is based on the finite volume method that solves the full compressible Navier-Stokes equations [15]. The ideal gas approximation is considered. The high-accuracy pressure staggering scheme (PRESTO) is enabled as the pressure interpolation scheme at the computational element faces, while the pressure-implicit with splitting of operators (PISO) pressure-velocity coupling scheme is applied. Double-precision parallel computation is essential to this analysis, therefore, the computational domain is divided into eight partitions and each simulation run is performed on an eight-core workstation and consumes a CPU time of about a month to reach its quasi-stationary response. A fixed time step is chosen 1=120 of the acoustic period and convergence is claimed when the scaled imbalance in each of the flow conservation equations summed 4
over all computational cells (residuals) becomes less than 10 in any two successive iterations, except for the energy equation, in which convergence is satisfied when residual becomes less than 106 . A typical transient behavior of the present computational model is shown in Fig. 3 at drive ratio of 2:0%. As the simulation proceeds, heat is being extracted from the stack cold end, carried by the introduced acoustic oscillations through the working gas and pumped into the stack hot end. This causes an increase in the temperature of the stack hot end (near the resonator rigid end or the velocity node), while the temperature decreases at the stack cold end. The variation with time of the temperature difference DT along with its rate of change are illustrated in Fig. 4. In support to the work of Worlikar et al. [3], the present temperature difference increases rapidly in the few early time steps, then its time rate of change sharply decreases with a noticeable exponential decay to asymptotically approach its steady-state value at sufficiently large simulation time. In the present simulation runs, the calculations are allowed to proceed till the time rate of change drops to less than 2:0% of its initial value while operating at the maximum drive ratio (2:0%). Numerical values of dðDTÞ=dt at large times are next extrapolated through an exponentially-decaying fitting process to yield the exact steady-state value of the ANSYS’s predicted temperature difference, as also shown in Fig. 4, according to:
t tf dðDTÞ ffi C 1 exp dt C2
ð3Þ
Here, C 1 and C 2 are the fitting constants, t f is the total simulation time and DT f is the corresponding temperature difference. For the particular case of 2:0% drive ratio, the fitting constants C 1 and C 2
295
are found to be 0:7034 and 0:9436, respectively. The steady-state value is then given by:
DT ss ffi DT f þ C 1 C 2
ð4Þ
3. Model validation The theoretical framework of both heat and work flows within a thermoacoustic couple in the short stack approximation has been extensively reviewed in literature [16,17], whereas fewer works were reported for stacks having arbitrarily shaped pore cross sections. Arnott [18] derived the time-averaged thermoacoustic heat flow equation due to hydrodynamic transport for square pores stack geometry, as follows:
AP 2 sinð2kxÞ F ðkj Þ I½Fðkj Þ þ rF ðkm Þ _ I 1C QðxÞ ¼ L 4qað1 þ rÞ Fðkm Þ I½Fðkj ÞFðkm ÞÞð1 þ rÞ
ð5Þ
where
C¼
dT m C p tanð2kxÞ dx xXa
ð6Þ
Here, the overbar denotes the time-average component of the second-order heat flow Q_ , whereas I refers to the imaginary part of its argument. A and X refer to the resonator cross-sectional area and the stack porosity, while x indicates the axial position, measured from the resonator closed end. The properties C p and r represent the specific heat and Prandtl number of the working gas, respectively. dT m =dx is the mean temperature gradient along the stack. Eq. (5) describes the heat flow in terms of the geometrical thermoviscous functions Fðkm Þ and Fðkj Þ and their complex conjugates F ðkm Þ and F ðkj Þ. According to Arnott [18], they are defined as:
Fðkm;j Þ ¼
64 X
1
p4 m;n odd m2 n2 Y mn ðkm;j Þ
ð7Þ
Y mn ðkm;j Þ ¼ 1 þ ıp2 ðm2 þ n2 Þ=ð4k2m;j Þ
ð8Þ
pffiffiffi km;j ¼ 2 2 r h =dm;j
ð9Þ
The dimensionless shear wave number km is defined as the ratio of the pore hydraulic radius r h to the gas viscous diffusion thickness dm . Similarly, kj corresponds to the thermal diffusion thickness dj . In the present work, the first 16 terms of the summation series were found sufficient to solve for the real and imaginary components of the thermoviscous functions. To further proceed with the theoretical estimation of the developed temperature gradient along the stack at no load, the timeaveraged heat transferred in the oscillating-gas due to the thermoacoustic heat-pumping effect from the cold to the hot stack
Fig. 3. Temperature history at the stack ends (Left) and variation of the temperature difference DT ¼ T H T C at 2.0% drive ratio (Right).
296
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
Fig. 4. The time rate of change of gas temperature at DR ¼ 2:0% is plotted in the left figure along with the exponential decay fit, while a comparison of the resulting DT m (closed circles) with theory (solid line), following Eq. (11), is indicated in the right plot. The corresponding DeltaEC predictions (open squares) are also included.
ends must be balanced, at dynamic equilibrium, with the diffusive conduction down the yet-generated temperature gradient through both the gas and solid domains. In a similar process to that followed by Wheatley et al. [17], the theoretical amount of heat transfer, presented in Eq. (5), can be expressed as a function of the gas and solid thermal conductivities K g and K s , respectively, as follows:
dT m dT m ¼ ðXK g þ ð1 XÞK s ÞA Q_ ðxÞ ¼ K eff A dx dx
ð10Þ
Substituting Eq. (5) into Eq. (10), solving for dT m =dx, and then integrating with respect to the axial direction over the stack length l yields the temperature difference DT m -Noting that in the shortstack approximation, the spatially-averaged components can be fairly replaced by their corresponding values at the stack center x
DT m
AP2L l sinð2kxÞ 4qað1þrÞ
K eff A þ
AP 2L sinð2kxÞ 4qað1þrÞ
I FFðkðkmjÞÞ
I FFðkðkmjÞÞ
I½Fðkj ÞþrF ðkm Þ C p tanð2kxÞ xXa I½Fðkj ÞFðkm Þð1þrÞ
ð11Þ
To validate the present numerical model, the theoretical values (Eq. (11)) are plotted in Fig. 4 along with the predicted quasisteady state temperature difference across the stack. For consistency and to gain confidence with the above theoretical analysis, the corresponding DeltaEC predictions are also included and found in perfect agreement with the linear theory. The numerical results, as predicted in the range of drive ratios from 0.25% to 2.0%, are qualitatively in a fair agreement with the linear theoretical behavior. Interestingly, the agreement with theory extends fairly well up to 2.0% drive ratio in this particular three-dimensional simulation of square pores, in contrary with corresponding two-dimensional analysis, such as that reported by Worlikar et al. [3] and Abd ElRahman and Abdel-Rahman [4], that declares divergence from linear theory almost at DR 1:0% for parallel-plate stack configuration. It is therefore concluded that the linear theory still holds up to the 2.0%-drive-ratio limit considering this particular stack geometry. It is worth mentioning that discrepancies observed at low drive ratios is partially attributed to less-accurate representation of the steady-state values DT m using the exponentially-decaying fit, in addition to typical numerical dissipation associated with the numerical upwind discretization scheme [19], employed by the ANSYS solver in the present simulation. 4. Results 4.1. Axial velocity profiles Through similar mathematical manipulation, the theoretical axial profiles of the oscillating flow velocity for 45 -increments
over one period are shown in Fig. 5 in comparison with the present numerical values for drive ratios 0.25% and 1.5%, respectively, and at two different longitudinal sections to illustrate the threedimensional flow behavior within the square pores. Nodal values of the axial velocity are first extracted from ANSYS solver and then passed to a specifically-developed MATLAB program for postprocessing. The profiles are plotted at the stack center (l/2-apart from the stack end) with a peak axial velocity of 12 m/s at the stack centerline while opearting at DR ¼ 1:5%, whereas no sign of flow reversal is captured in the present prediction. Fig. 5 shows excellent agreement between the present numerical results and the theoretical axial profiles of the oscillating flow velocity for both low and moderate drive ratios, as predicted at the stack center. This demonstrates the existence of a onedimensional flow behavior within the stack and thus supports the validity of the 1-D linear theory up to this limit (DR ¼ 2:0%). It also suggests the suppression of DC-streaming within the square pores; because otherwise, any induced 3-D motion would influence the axial velocity distribution and cause departure between the theory and the computations. This is in contrast with the conclusion of Worlikar et al. [3], in their 2-D computations considering a stack consisting of parallel plates, which indicates the occurrence of a relatively important 2-D motion within the gap that leads to the departure between the theory and the computations. To our knowledge, this has not been reported before. Furthermore, typical laminar flow oscillations are noticed at both drive ratios, and the system is clearly linear in the present simulation runs. To further demonstrate the propagation of laminar flow, both the dimensionless frequency parameter and the Reynolds numbers are estimated. For instance, at 1.5% drive ratio, the amplitude of the cross-sectional mean-velocity variation ^ 8 m/s. Thus, for a viscous penetration depth of dm 0:26 mm, U 2
the frequency parameter and Reynolds number are h x=m 25:4 ^ m =m 16:74, respectively. According to Hino et al. [14], these and Ud values are fairly low for a transition to turbulence to take place. 4.2. Energy flux density A better description of the heat transport process from the stack cold end to the stack hot end is here introduced through the calculation of the time-averaged energy flux density over the last acoustic period. A specific user defined function is developed to record the instantaneous nodal values of the gas temperature, velocity and temperature gradients and thus, enable the calculation of the time-averaged (mean) components of the energy flux density, according to:
H ¼ C p qðT T m Þu K g rT
ð12Þ
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
297
Fig. 5. Velocity Profiles at drive ratios 0.25% (Top) and 1.5% (Bottom), as compared to the Arnott’s extended linear theory (Solid lines) for oscillating flow in square pores at an axial distance of l=2 from the stack cold end. The profiles are plotted at two different longitudinal sections of the gas square channel, corresponding to 0:25h and 0:5h, respectively.
Further post-processing is achieved using a specific MATLAB program. Fig. 6 shows the streamlines -superimposed by the vector plots- of the mean energy flux density for two different drive ratios at two different longitudinal planes, as predicted in the present runs. The mean energy streamlines are plotted, such that their starting points are positioned at the inner surfaces of the stack cold end, for the sake of visualization. The vectors lengths are automatically scaled for visualization, whereas, the prediction reveals a maximum vector length of 123.89 and 672.82 W/m2 at 0.75% and 2.0%, respectively. It is clearly shown in figure that, because of the thermoacoustic heat pumping effect, the thermal energy is first extracted from the stack inner walls at the cold end, swirls in an eddy-like motion till finds its straight-path through the gas along the stack to finally get rejected into the hot end. In Fig. 6(B) and (D), the predicted energy flux densities are plotted over a longitudinal section located quarter-way the channel width h, while Fig. 6(A) and (C) correspond to those calculated at the longitudinal mid-plane, as shown in the inset figure. The figures illustrate the existence of large-areas of recirculating mean energy fluxes, intensity of which increases as the drive ratio grows from 0.75% to 2.0%, as first noticed by comparing Fig. 6(A) and (C). It is also observed that the axial extension of the recirculating zones closely corresponds to the cross-sectional mean particle dis^ x. Remarkably, such recirculating zones become placement U= more intense while moving from the mid-plane behavior to the 0:25h-corresponding section, such that longer recirculating pathlines are exhibited before getting back into the stack. In addition to that, limited asymmetrical behavior in the energy flux density is
captured at higher drive ratio at both longitudinal sections with more intensive appearance towards the center. Furthermore, the energy fluxes predicted at the channel midplane shown in Fig. 6(A) and (C) are found to be more concentrated rather than detached, with respect to those illustrated in Fig. 6(B) and (D). This is in part consistent with -and is possibly caused by- the induced non-zero DC-streaming and its associated vortical motion at the stack extremities with similar dependence on the imposed drive ratio, as illustrated in the next section. This might in turn assume a significant contribution to the overall thermal energy losses. Interestingly, this is not supported by the above agreement of the mean temperature difference and axial velocity profiles with theoretical predictions at DR ¼ 2:0%. Existing 2-D simulations [3] on parallel-plate stacks reported that the departure between the theory and the computations, in their evaluation for the steady-state stack temperature difference at DR 1:0%, results from the growing mean energy paths and the generated 2-D motion within the gap. In contrary with that, the present 3-D simulation clearly indicates the suppression of the multi-dimensional flow behavior within the square pores as demonstrated through the agreement of the axial velocity profiles with theory at the stack center. The other factor, although evident in the present work, its contribution does not strongly impact the bulk oscillating-flow behavior and thus does not cause significant departure from the linear theory in the studied range of drive ratio. It is therefore concluded that the effect of the eliminated multi-dimensional flow behavior within the stack square pores on the adherence to the 1-D linear theory is more important than
298
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
Fig. 6. Vector plots superimposing streamlines (Solid red lines) of the mean energy flux density at drive ratios of 0.75% (A and B) and 2.0% (C and D). These plots are presented for two different longitudinal sections parallel to (and 0:5h- and 0:25h-apart from) the channel walls, as illustrated for instance in the insets. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
alleviating the increasingly-diffusing large-areas of recirculating mean energy fluxes outside the stack pores. 4.3. Flow streaming In addition to the estimation of the mean energy-flux density, the ANSYS solver enables us to calculate the mean velocity components and flow vorticity to gain insight into the influence of the drive ratio on the induced non-linear flow streaming particularly at the stack extremities. Fig. 7 shows the contour plots of the calculated vorticity field superimposed by vector plots of its respective flow streaming at drive ratios 0.25% and 2.0%, as similarly predicted at the two longitudinal sections located parallel to (and 0:5h- and 0:25h-apart from) the channel walls. The vectors lengths are automatically scaled in the MatLab program for better visualizing the mean flow fields, whereas for instance, the results shown in Fig. 7(A) and (C) indicate a maximum vector length of 0.11 and 2.46 m/s at drive ratios 0.25% and 2.0%, respectively. The right plots of Fig. 7 also shows the calculation of the instantaneous velocity profiles at the stack left extremity to help understanding the additional entrance effect. Flatter velocity profiles are captured during the second half of the acoustic cycle while the gas particles are moving back into the stack channels because of the
non-fully-developed flow behavior at the stack entrance. Such effect is noticeably stronger at the larger drive ratio 2.0%. The streaming results, as illustrated in Fig. 7, also reveal the structure of the evolving three-dimensional vortical motions and their associated flow rotationalities at the corners of each square pore because of the sudden expansion and contraction the flow exhibits as getting in and out of the stack square pores. Concentrated vortices are found when operating at the lower drive ratio, whereas intense vortex shedding is shown at the larger drive ratio. Similar to the energy flux density, the predicted vortical pathlines are shown elongated and detached from the stack extremities at the larger drive ratio. Therefore, in the present particular 3D simulation of the thermoacoustic fields within a square stack, it is concluded that stronger currents of mean swirling motions, and thus vorticity, are captured at larger drive ratios, however, such non-linear effects are limited in terms of vorticity transport- and remarkably localized nearby the stack extremities. No evidence of significant streaming activity has been found at the stack middle section at both drive ratios that could drastically impact the one-dimensionality oscillating-flow condition within the stack; and thus, help reduce the non-linear flow losses and retain efficient transfer of heat within the square pores from the stack cold to hot ends.
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
299
Fig. 7. The left figures show vector plots of the flow streaming superimposing contours of the ‘mean’ vorticity field at drive ratios of 0.25% (A and B) and 2.0% (C and D). These plots are presented for two different longitudinal sections parallel to (and 0:5h- and 0:25h-apart from) the channel walls, following exactly the same arrangement shown in Fig. 6. The arrows lengths are adjusted automatically within the MatLab program for the sake of visualization. The vorticity is expressed in the unit of s1. The right figures plot the velocity profiles at eight different instances over a cross-sectional plane located exactly at the stack cold end.
5. Summary A new three-dimensional finite-volume model is presented that solves the full non-linear thermo-hydrodynamic laminar flow equations in three-dimensions and examines the heat transport between the thermoacoustically oscillating gas and the surrounding surfaces in a thermoacoustic stack made of square pores. The present model enables us to investigate the temperature variation with time at the stack ends along with the evolving axial velocity profiles at the stack center and the mean (time-averaged) energy flux density nearby the stack ends. The flow streaming and mean vorticity field are also calculated and presented in this study to gain insight into the role of the developed acoustic non-linearties
in determining the departure of the stack temperature difference from the linear theory. Here, the theoretical mean temperature difference is deduced based on the Arnott’s extended linear theory for stacks having arbitrarily shaped pore cross sections. In the present simulation runs, the captured instantaneous velocity profiles and the variation with the drive ratio of the estimated mean temperature difference DT m between the stack ends show good agreement with the Arnott’s extended linear theory [18] up to a drive ratio of 2.0%. In contrary with existing 2-D simulations carried out on parallel-plate stack, the present results indicates the prevalence of the 1-D flow within the stack pores and limits the multi-dimensional flow behavior to the stack extremities, such that the departure of the present 3-D computations from the
300
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
theory has been effectively shifted away. Although the swirling pathlines of the mean energy flux densities are clearly evident, and even become more intense with the increase of the drive ratio, in the neighborhood of the stack ends, it is found in the present work that its influence is weak enough to promote departure between the theory and the computations. Therefore, the present work suggests that more useful transfer of heat on and off the stack is expected with the square-pore geometry than the corresponding parallel-plate structure. This study will be extended in our future work to include higher drive ratios, in which the potential of flow transition from laminar to turbulence is no longer trivial.
References [1] N. Cao, J.R. Olson, G.W. Swift, S. Chen, Energy flux density in a thermoacoustic couple, J. Acoust. Soc. Am. 99 (6) (1996) 3456–3464. [2] A.S. Worlikar, O.M. Knio, Numerical simulation of a thermoacoustic refrigerator: I. Unsteady adiabatic flow around the stack, J. Comput. Phys. 127 (2) (1996) 424–451. [3] A.S. Worlikar, O.M. Knio, R. Klein, Numerical simulation of a thermoacoustic refrigerator: II. Stratified flow around the stack, J. Comput. Phys. 144 (2) (1998) 299–324. [4] A.I. Abd El-Rahman, E. Abdel-Rahman, Computational fluid dynamics simulation of a thermoacoustic refrigerator, J. Thermophys. Heat Transf. 28 (1) (2014) 78–86. [5] P. Lotton, P. Blanc-Benon, M. Bruneau, V. Gusev, S. Duffourd, M. Mironov, G. Poignand, Transient temperature profile inside thermoacoustic refrigerators, Int. J. Heat Mass Transf. 52 (21–22) (2009) 4986–4996.
[6] P. Blanc-Benon, E. Besnoin, O. Knio, Experimental and computational visualization of the flow field in a thermoacoustic stack, CR Mec. 331 (1) (2003) 17–24. [7] H. Ishikawa, D.J. Mee, Numerical investigations of flow and energy fields near a thermoacoustic couple, J. Acoust. Soc. Am. 111 (2) (2002) 831–839. [8] D. Marx, P. Blanc-Benon, Computation of the mean velocity field above a stack plate in a thermoacoustic refrigerator, CR Mec. 332 (11) (2004) 867–874. [9] D. Marx, P. Blanc-Benon, Computation of the temperature distortion in the stack of a standing-wave thermoacoustic refrigerator, J. Acoust. Soc. Am. 118 (5) (2005) 2993–2999. [10] L. Zoontjens, C. Howard, A. Zander, B. Cazzolato, Numerical study of flow and energy fields in thermoacoustic couples of non-zero thickness, Int. J. Therm. Sci. 48 (4) (2009) 733–746. [11] S.H. Tasnim, R.A. Fraser, Computation of the flow and thermal fields in a thermoacoustic refrigerator, Int. Commun. Heat Mass Transf. 37 (7) (2010) 748–755. [12] G. Poignand, B. Lihoreau, P. Lotton, E. Gaviot, M. Bruneau, V. Gusev, Optimal acoustic fields in compact thermoacoustic refrigerators, Appl. Acoust. 68 (6) (2007) 642–659. [13] C. Fan, B.-T. Chao, Unsteady, laminar, incompressible flow through rectangular ducts, Z. Angew. Math. Phys. ZAMP 16 (3) (1965) 351–360. [14] M. Hino, M. Kashiwayanagi, A. Nakayama, T. Hara, Experiments on the turbulence statistics and the structure of a reciprocating oscillatory flow, J. Fluid Mech. 131 (1983) 363–400. [15] ANSYS Academic Research, Release 15.0, Help System, User Guide, ANSYS, Inc. [16] G.W. Swift, Thermoacoustic engines, J. Acoust. Soc. Am. 84 (4) (1988) 1145– 1180. [17] J. Wheatley, T. Hofler, G.W. Swift, A. Migliori, An intrinsically irreversible thermoacoustic heat engine, J. Acoust. Soc. Am. 74 (1) (1983) 153–170. [18] W.P. Arnott, H.E. Bass, R. Raspet, General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross sections, J. Acoust. Soc. Am. 90 (6) (1991) 3228–3237. [19] O. Zikanov, Essential Computational Fluid Dynamics, Wiley, 2010.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A 3D investigation of thermoacoustic fields in a square stack Ahmed I. Abd El-Rahman ⇑, Waleed A. Abdelfattah, Mahmoud A. Fouad Department of Mechanical Power Engineering, Cairo University, Cairo University Rd., Giza 12613, Egypt
a r t i c l e
i n f o
Article history: Received 10 April 2016 Received in revised form 2 December 2016 Accepted 5 December 2016
Keywords: Thermoacoustic refrigerator Square stack 3D simulation Energy flux density Flow streaming
a b s t r a c t The thermoviscous behavior of the oscillating gas within the porous medium of a thermoacoustic refrigerator enables the conversion of sound into heat in the process of typical standing-wave thermoacoustic refrigeration systems. Several nonlinear mechanisms, such as harmonics generation, self-induced streaming and possible turbulence, cause complicated flow behavior and influence the performance of such devices. Few analytical and numerical approximations (Cao et al., 1996; Worlikar and Knio, 1996; Worlikar et al., 1998) describe the flow field and the energy flux density in standing devices comprising stacks of parallel plates, but almost no 3-D simulation has been developed that models the largeamplitude excitations and enables the prediction of the oscillating-flow behavior within a porous medium made of square channels. Here, we build on existing effort (Abd El-Rahman and Abdel-Rahman, 2014) and report a computational fluid dynamics analysis of a Helium-filled half-wavelength thermoacoustic refrigerator. The finite volume method is used, and the solid and gas domains are represented by large numbers of hexahedral elements. The calculations assume a 3-D periodic cell structure to reduce the computational cost and apply the dynamic mesh technique to account for the adiabatically oscillating wall boundaries. The simulation uses an implicit time integration of the full unsteady compressible flow equations along with the conjugate heat transfer algorithm. For the sake of simulation, the geometry and operating conditions closely follow the experimental setup of Lotton et al. (2009) [5]. A typical run involves about 200,000 computational cells, whereas the working frequency and the applied drive ratio are 592.5 Hz and 2.0%, respectively. The stationary system response is predicted and the numerical values are compared to theory. The results show good agreement with theoretical behavior in the linear regime of drive ratios up to 2.0% while a maximum cooling effect of 10 degrees is captured between the stack ends. This simulation provides an interesting tool for understanding the transient flow behavior, mean energy-flux density and flow streaming, and helps building 3-D computational fluid dynamics models for thermoacoustic devices. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction A typical standing-wave thermoacoustic refrigerator (TAR) consists of a long resonator that is usually filled with an inert gas, a porous medium called ‘Stack’ and two heat exchangers. The stack of a standing-wave refrigerator is usually short and made of ceramic material, whereas the resonator is closed at one end and bounded by a loudspeaker at the other end, thus, allowing the propagation of a half-wavelength acoustic excitation. The hot and cold heat exchangers are made of copper to provide efficient heat transfer between the working inert gas and the external heat source and sink respectively. TARs are interesting because they have no moving parts, unlike conventional refrigerators, and almost no environmental impact exists as they rely on clean ⇑ Corresponding author. E-mail address: [email protected] (A.I. Abd El-Rahman). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.12.015 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.
conversion of acoustic and heat energies. Their fabrication process is rather simpler and sizes span wide variety of length scales. In typical low-pressure excitation systems, the thermoacoustic behavior is well described by the 1-D linear theory. However, practical systems fall into the high-pressure regime of acoustic excitation, where non-linear behavior is dominant. The following recent studies showed substantial deviation of most cooling curves from the linear theory. They rather focused on modeling thermoacoustic devices to help predict the general performance, capture the thermal and viscous flow losses in the device, and define key variables affecting the design considerations and system optimization, such as stack porosity, size and location, working gas and appropriate operating mean pressure. Cao et al. [1] was the first to construct a first-order temporal and spatial discretization scheme to solve for the time-averaged energy flux density in the gas. In their model, the stack was represented by short (compared with the acoustic wavelength)
293
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
isothermal parallel plates. Cao’s attempt was followed by the significant effort of Worlikar et al. [2,3] who considered the lowMach number approximation of the full Navier-Stokes equations and solved the streamfunction-vorticity equation. They were able to successfully reduce the full computational domain into a single periodic thermoacoustic couple. Although limited to twodimensional analysis, short stack lengths and constant gas properties, their approach was then traced in similar fashion by other recent workers, such as Blanc-Benon et al. [6]. Blanc-Benon performed a particle image velocimetry (PIV) measurements of the flow field and reported close agreement with the low-Mach number simulation results. Ishikawa and Mee [7] extended Cao analysis to longer plates and investigated the effects of plate spacing using the PHOENICS finite-volume method, however, their single-precision calculations caused significant energy imbalances and restricted the overall accuracy in their heat fluxes results. Latter, Marx and BlancBenon [8,9] considered the time-averaged velocity field above the stack plate and observed the acoustic flow streaming above the plate. They further identified the minimum plate length, as compared to the particle displacement, necessary for the generated temperature harmonics to get displayed above the entire plate surface. Although accurate, this model is too expensive computationally to limit its applicability to devices operating at high-frequency of 20 kHz, which is useful for miniaturization purposes. Zoontjens et al. [10] was the first to implement a FLUENT model to simulate the flow and energy fields in a TAR couple of non-zero thickness at a wide range of drive ratio (1.7–8.5%) and 200 Hz operating frequency. They assumed symmetric computational domain and investigated the time-averaged heat transfer through the stack material. Their model however used a first-order discretization scheme which in turn affected the accuracy of their numerical results. This was recently followed by the work of Tasnim and Fraser [11] who introduced a two-dimensional model using the STARCD finite volume solver. He assumed a periodic thermoacoustic couple and described the laminar-flow and thermal fields at low drive ratio of 0.7%. To our knowledge, no 3-D simulation has been developed that directly models the large-amplitude excitations of standing-wave TARs. Building on existing effort [4], the objective of the present work is to extend the two-dimensional analysis of standing-wave thermoacoustic couples into three-dimensional simulation of the oscillating-flow behavior within a porous medium made of square channels. The next objective includes the investigation of the thermoviscous interactions along the stack and at the stack ends, and to compare the axial velocity profiles with the corresponding theoretical behavior at low and moderate drive ratios. Furthermore, the energy flux density is evaluated and plotted at two different longitudinal sections to illustrate the three-dimensional flow behavior in square pores.
2. Reference model 2.1. Physical model Inspired by the setup of Lotton et al. [5], the present thermoacoustic refrigerator consists of an 86-cm-long Helium-filled resonator that has a square cross section with an edge length of eight centimeters, as schematically drawn in Fig. 1. The halfwavelength resonator has a loudspeaker placed at one end and is closed at the other end. The porous media ‘stack’ is made of cordierite (ceramic material) and consists of a large number of small square channels (600 CPSI), aligned with the direction of wave propagation. Here, each channel is 10-mm-in-length and has a transversal width h ¼ H t s of 0.92 mm and a wall thickness ts
Fig. 1. Schematic of the thermoacoustic refrigerator is shown on top. Bottom: The representative element volume (REV) is drawn with a couple of periodic surfaces along with an acoustically-matching adiabatic boundary conditions imposed at both ends (the edge length of the square pore is h ¼ H t s ).
of 0.12 mm, thus resulting in a stack porosity / of 0.787. The model is constructed such that the stack center is positioned at a distance 14 cm ( k=12:3) from the resonator closed end. According to Poignand et al. [12], such stack position corresponds to the maximum cooling effect when using air as the working gas. Heat exchangers are not included in the present model setup for simplicity. The system is excited with an acoustic frequency of 592.5 Hz and acoustic pressure amplitude in the range from 250 to 2000 Pa while operating at atmospheric pressure and ambient temperature. Details about flow conditions and material thermophysical properties are remarked in Table 1. 2.2. Numerical model It is unreasonable to model the entire domain of the thermoacoustic resonator together with the square stack because of their excessively dispersed length scales. Therefore, an element volume that includes one gas channel is chosen to represent the oscillating-flow behavior within the stack and at the stack ends, as shown in Fig. 1. Such representative element volume (REV) can be thought as a unit cell that is surrounded by large number of similarly repeated cells (Or, images). To account for those image cells in both transversal and spanwise directions, the present REV is considered bounded by periodic surfaces in both directions. In the axial direction, however, the 3-D computational domain is then allowed to extend with a total length of 4.64 times the stack length l. Acoustically matching boundary conditions are thus imposed at both left and right ends of the chosen computational domain,
Table 1 Flow conditions and material thermophysical properties. Property
Value
Units
Flow Conditions Operating frequency, f Ambient Temperature, T m Mean Pressure, pm
592:5 300 101:325
Hz K kPa
0:1625
kg m3
Gas properties: Helium Density, qf Prandtl number, Pr Thermal Conductivity, kf Heat capacity, cp;f Dynamic viscosity, lf
0:68 0:152 5193
Ratio of specific heats, c
1:99 105 1:67
Stack properties: Cordierite Thermal Conductivity, ks Heat capacity, cp;s Density, qs
0:19 1000 1400
W m1 K1 J kg1 K1 kg m1 s1
W m1 K1 J kg1 K1 kg m3
294
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
assuming an ideal acoustic wave propagation. Although this approximation does not allow for accurate replication of the physical process at these boundaries, it can still suggest appropriate modeling at the far field of the stack plates provided that the distance from the two vertical surfaces to the stack ends is big enough, as argued by Cao et al. [1]. For instance, the length of the computational domain to the stack length is chosen to ensure that the flow at the domain left and right boundaries can be accurately modeled as a one-dimensional planar wave. Fig. 2 presents the domain meshing in the transversal and longitudinal directions. The domain is carefully discretized to properly capture the thermoviscous interactions of the oscillating-flow with the stack walls and to maintain grid-independent system response. The flow is assumed laminar. Such flow approximation will be elaborated in the results section through the specification of two key parameters [13,14]: the dimensionless frequency parameter
Table 2 Main meshing details, as illustrated in Fig. 2. Parameter
Value
Units
Fluid Domain - Inside the Stack Boundary cell size (MinjMax)
ð0:03Þ2 0:024jð0:03Þ2 0:12
mm3
ð0:06Þ 0:024jð0:06Þ 0:12 21 21 84
mm3
Number of cells Fluid Domain - Outside the Stack Boundary cell size (MinjMax)
ð0:03Þ2 0:025jð0:03Þ2 0:15
mm3
Centerline cell size (MinjMax)
Centerline cell size (MinjMax) Number of cells Total number of cells Solid Domain Corner cell size (MinjMax) Middle cell size (MinjMax) Total number of cells
2
2
2
2
ð0:06Þ 0:025jð0:06Þ 0:15 25 25 244 189,544
–
mm3 – –
ð0:03Þ2 0:024jð0:03Þ2 0:12
mm3
ð0:03Þ2 0:06j0:03 0:06 0:12 14,784
mm3 –
2
h x=m (Also, referred as the Womeresley number for circular pipe)
^ m =m, based on the amplitude and the flow Reynolds number Re ¼ Ud ^ and the viscous of the cross-sectional mean-velocity variation U pffiffiffiffiffiffiffiffiffiffiffiffi diffusion thickness dm ¼ 2m=x. Here, m and x are the gas viscous diffusivity and the angular frequency, respectively. In the transversal direction, the mesh is distinguished into two concentric regimes, depth of each corresponds to the thermal difpffiffiffiffiffiffiffiffiffiffiffiffiffi fusion thickness dj ¼ 2a=x, where a is the gas thermal diffusivity. In the first regime, the flow experiences steep velocity and temperature gradients, therefore, eight computational cells are enforced to capture the flow behavior very near to the wall. Less steep gradients are expected in the second regime, thus, five larger cells are assumed sufficient since this regime lies in the stack core, where poor transversal fluctuation in the flow properties generally evolves. In the longitudinal direction, a structured hexahedral mesh is maintained such that the cell height varies gradually and smoothly from one to five times the cell width throughout the entire domain, as indicated in Fig. 2(Left). For instance, the resulting total number of cells is around 2 105 . The main meshing details are further presented in Table 2. In the present simulation, the adiabatic thermal boundary condition is applied at the matching boundaries whereas the viscous interaction with the stack inner walls satisfies the no-slip boundary condition. At the solid-fluid interface, thermally coupled boundary conditions are imposed to simulate the thermal interaction between the solid stack walls and the surrounding gas
particles through the conjugate heat transfer algorithm. In such case, the solver requires that the temperature T s and the heat transfer to the wall boundary from a solid element match the gas side temperature T g and heat transfer, such that T s ¼ T g at the stack inner walls and:
Kg
@T g @T s ¼ Ks @n @n
ð1Þ
In Eq. (1), the heat transfer is expressed in terms of the gas and solid thermal conductivities K g and K s , whereas n indicates to the local coordinate normal to the wall. Next, in order to generate the acoustically matching boundary condition, the layering dynamic meshing technique is applied, as detailed in the ANSYS FLUENT user documentation [15]. A preprocessing program is developed to define the instantaneous axial motions of the matching boundaries U MB , consistent with the local ideal standing wave propagation:
U MB ¼ U 0 sinðkxÞ cosðxtÞ
ð2Þ
Here, U 0 represents the maximum amplitude of the acoustic velocity, while k and x refer to the wave number and the axial position of the matching boundaries, measured from the resonator closed end, respectively. x is the angular velocity and t indicates the simulation time. Further, U 0 can be expressed in terms of the pressure amplitude at the pressure antinodes of the standing wave (i.e. resonator
Fig. 2. The left figure presents the 3-D meshing of the computational gas domain, as spanned from inside the stack to the outside of it. The solid mesh is suppressed for clarity. However, the right plot shows a side view -transversal cut- of the whole computational domain (Inset: dj 0:31 mm). Here, the yellow frame represents the solid stack, while the gray structured mesh corresponds to the gas domain with a minimum cell size of 0.03 mm and a maximum of 0.06 mm. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
closed end pL ), the mean density q, and the speed of sound a, such that U 0 ¼ pL =qa. The present analysis uses the ANSYS FLUENT second-order upwind implicit-time algorithm pressure-based CFD solver with the conjugate heat transfer algorithm along with a second-order formulation for the temporal discretization. The CFD solver is based on the finite volume method that solves the full compressible Navier-Stokes equations [15]. The ideal gas approximation is considered. The high-accuracy pressure staggering scheme (PRESTO) is enabled as the pressure interpolation scheme at the computational element faces, while the pressure-implicit with splitting of operators (PISO) pressure-velocity coupling scheme is applied. Double-precision parallel computation is essential to this analysis, therefore, the computational domain is divided into eight partitions and each simulation run is performed on an eight-core workstation and consumes a CPU time of about a month to reach its quasi-stationary response. A fixed time step is chosen 1=120 of the acoustic period and convergence is claimed when the scaled imbalance in each of the flow conservation equations summed 4
over all computational cells (residuals) becomes less than 10 in any two successive iterations, except for the energy equation, in which convergence is satisfied when residual becomes less than 106 . A typical transient behavior of the present computational model is shown in Fig. 3 at drive ratio of 2:0%. As the simulation proceeds, heat is being extracted from the stack cold end, carried by the introduced acoustic oscillations through the working gas and pumped into the stack hot end. This causes an increase in the temperature of the stack hot end (near the resonator rigid end or the velocity node), while the temperature decreases at the stack cold end. The variation with time of the temperature difference DT along with its rate of change are illustrated in Fig. 4. In support to the work of Worlikar et al. [3], the present temperature difference increases rapidly in the few early time steps, then its time rate of change sharply decreases with a noticeable exponential decay to asymptotically approach its steady-state value at sufficiently large simulation time. In the present simulation runs, the calculations are allowed to proceed till the time rate of change drops to less than 2:0% of its initial value while operating at the maximum drive ratio (2:0%). Numerical values of dðDTÞ=dt at large times are next extrapolated through an exponentially-decaying fitting process to yield the exact steady-state value of the ANSYS’s predicted temperature difference, as also shown in Fig. 4, according to:
t tf dðDTÞ ffi C 1 exp dt C2
ð3Þ
Here, C 1 and C 2 are the fitting constants, t f is the total simulation time and DT f is the corresponding temperature difference. For the particular case of 2:0% drive ratio, the fitting constants C 1 and C 2
295
are found to be 0:7034 and 0:9436, respectively. The steady-state value is then given by:
DT ss ffi DT f þ C 1 C 2
ð4Þ
3. Model validation The theoretical framework of both heat and work flows within a thermoacoustic couple in the short stack approximation has been extensively reviewed in literature [16,17], whereas fewer works were reported for stacks having arbitrarily shaped pore cross sections. Arnott [18] derived the time-averaged thermoacoustic heat flow equation due to hydrodynamic transport for square pores stack geometry, as follows:
AP 2 sinð2kxÞ F ðkj Þ I½Fðkj Þ þ rF ðkm Þ _ I 1C QðxÞ ¼ L 4qað1 þ rÞ Fðkm Þ I½Fðkj ÞFðkm ÞÞð1 þ rÞ
ð5Þ
where
C¼
dT m C p tanð2kxÞ dx xXa
ð6Þ
Here, the overbar denotes the time-average component of the second-order heat flow Q_ , whereas I refers to the imaginary part of its argument. A and X refer to the resonator cross-sectional area and the stack porosity, while x indicates the axial position, measured from the resonator closed end. The properties C p and r represent the specific heat and Prandtl number of the working gas, respectively. dT m =dx is the mean temperature gradient along the stack. Eq. (5) describes the heat flow in terms of the geometrical thermoviscous functions Fðkm Þ and Fðkj Þ and their complex conjugates F ðkm Þ and F ðkj Þ. According to Arnott [18], they are defined as:
Fðkm;j Þ ¼
64 X
1
p4 m;n odd m2 n2 Y mn ðkm;j Þ
ð7Þ
Y mn ðkm;j Þ ¼ 1 þ ıp2 ðm2 þ n2 Þ=ð4k2m;j Þ
ð8Þ
pffiffiffi km;j ¼ 2 2 r h =dm;j
ð9Þ
The dimensionless shear wave number km is defined as the ratio of the pore hydraulic radius r h to the gas viscous diffusion thickness dm . Similarly, kj corresponds to the thermal diffusion thickness dj . In the present work, the first 16 terms of the summation series were found sufficient to solve for the real and imaginary components of the thermoviscous functions. To further proceed with the theoretical estimation of the developed temperature gradient along the stack at no load, the timeaveraged heat transferred in the oscillating-gas due to the thermoacoustic heat-pumping effect from the cold to the hot stack
Fig. 3. Temperature history at the stack ends (Left) and variation of the temperature difference DT ¼ T H T C at 2.0% drive ratio (Right).
296
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
Fig. 4. The time rate of change of gas temperature at DR ¼ 2:0% is plotted in the left figure along with the exponential decay fit, while a comparison of the resulting DT m (closed circles) with theory (solid line), following Eq. (11), is indicated in the right plot. The corresponding DeltaEC predictions (open squares) are also included.
ends must be balanced, at dynamic equilibrium, with the diffusive conduction down the yet-generated temperature gradient through both the gas and solid domains. In a similar process to that followed by Wheatley et al. [17], the theoretical amount of heat transfer, presented in Eq. (5), can be expressed as a function of the gas and solid thermal conductivities K g and K s , respectively, as follows:
dT m dT m ¼ ðXK g þ ð1 XÞK s ÞA Q_ ðxÞ ¼ K eff A dx dx
ð10Þ
Substituting Eq. (5) into Eq. (10), solving for dT m =dx, and then integrating with respect to the axial direction over the stack length l yields the temperature difference DT m -Noting that in the shortstack approximation, the spatially-averaged components can be fairly replaced by their corresponding values at the stack center x
DT m
AP2L l sinð2kxÞ 4qað1þrÞ
K eff A þ
AP 2L sinð2kxÞ 4qað1þrÞ
I FFðkðkmjÞÞ
I FFðkðkmjÞÞ
I½Fðkj ÞþrF ðkm Þ C p tanð2kxÞ xXa I½Fðkj ÞFðkm Þð1þrÞ
ð11Þ
To validate the present numerical model, the theoretical values (Eq. (11)) are plotted in Fig. 4 along with the predicted quasisteady state temperature difference across the stack. For consistency and to gain confidence with the above theoretical analysis, the corresponding DeltaEC predictions are also included and found in perfect agreement with the linear theory. The numerical results, as predicted in the range of drive ratios from 0.25% to 2.0%, are qualitatively in a fair agreement with the linear theoretical behavior. Interestingly, the agreement with theory extends fairly well up to 2.0% drive ratio in this particular three-dimensional simulation of square pores, in contrary with corresponding two-dimensional analysis, such as that reported by Worlikar et al. [3] and Abd ElRahman and Abdel-Rahman [4], that declares divergence from linear theory almost at DR 1:0% for parallel-plate stack configuration. It is therefore concluded that the linear theory still holds up to the 2.0%-drive-ratio limit considering this particular stack geometry. It is worth mentioning that discrepancies observed at low drive ratios is partially attributed to less-accurate representation of the steady-state values DT m using the exponentially-decaying fit, in addition to typical numerical dissipation associated with the numerical upwind discretization scheme [19], employed by the ANSYS solver in the present simulation. 4. Results 4.1. Axial velocity profiles Through similar mathematical manipulation, the theoretical axial profiles of the oscillating flow velocity for 45 -increments
over one period are shown in Fig. 5 in comparison with the present numerical values for drive ratios 0.25% and 1.5%, respectively, and at two different longitudinal sections to illustrate the threedimensional flow behavior within the square pores. Nodal values of the axial velocity are first extracted from ANSYS solver and then passed to a specifically-developed MATLAB program for postprocessing. The profiles are plotted at the stack center (l/2-apart from the stack end) with a peak axial velocity of 12 m/s at the stack centerline while opearting at DR ¼ 1:5%, whereas no sign of flow reversal is captured in the present prediction. Fig. 5 shows excellent agreement between the present numerical results and the theoretical axial profiles of the oscillating flow velocity for both low and moderate drive ratios, as predicted at the stack center. This demonstrates the existence of a onedimensional flow behavior within the stack and thus supports the validity of the 1-D linear theory up to this limit (DR ¼ 2:0%). It also suggests the suppression of DC-streaming within the square pores; because otherwise, any induced 3-D motion would influence the axial velocity distribution and cause departure between the theory and the computations. This is in contrast with the conclusion of Worlikar et al. [3], in their 2-D computations considering a stack consisting of parallel plates, which indicates the occurrence of a relatively important 2-D motion within the gap that leads to the departure between the theory and the computations. To our knowledge, this has not been reported before. Furthermore, typical laminar flow oscillations are noticed at both drive ratios, and the system is clearly linear in the present simulation runs. To further demonstrate the propagation of laminar flow, both the dimensionless frequency parameter and the Reynolds numbers are estimated. For instance, at 1.5% drive ratio, the amplitude of the cross-sectional mean-velocity variation ^ 8 m/s. Thus, for a viscous penetration depth of dm 0:26 mm, U 2
the frequency parameter and Reynolds number are h x=m 25:4 ^ m =m 16:74, respectively. According to Hino et al. [14], these and Ud values are fairly low for a transition to turbulence to take place. 4.2. Energy flux density A better description of the heat transport process from the stack cold end to the stack hot end is here introduced through the calculation of the time-averaged energy flux density over the last acoustic period. A specific user defined function is developed to record the instantaneous nodal values of the gas temperature, velocity and temperature gradients and thus, enable the calculation of the time-averaged (mean) components of the energy flux density, according to:
H ¼ C p qðT T m Þu K g rT
ð12Þ
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
297
Fig. 5. Velocity Profiles at drive ratios 0.25% (Top) and 1.5% (Bottom), as compared to the Arnott’s extended linear theory (Solid lines) for oscillating flow in square pores at an axial distance of l=2 from the stack cold end. The profiles are plotted at two different longitudinal sections of the gas square channel, corresponding to 0:25h and 0:5h, respectively.
Further post-processing is achieved using a specific MATLAB program. Fig. 6 shows the streamlines -superimposed by the vector plots- of the mean energy flux density for two different drive ratios at two different longitudinal planes, as predicted in the present runs. The mean energy streamlines are plotted, such that their starting points are positioned at the inner surfaces of the stack cold end, for the sake of visualization. The vectors lengths are automatically scaled for visualization, whereas, the prediction reveals a maximum vector length of 123.89 and 672.82 W/m2 at 0.75% and 2.0%, respectively. It is clearly shown in figure that, because of the thermoacoustic heat pumping effect, the thermal energy is first extracted from the stack inner walls at the cold end, swirls in an eddy-like motion till finds its straight-path through the gas along the stack to finally get rejected into the hot end. In Fig. 6(B) and (D), the predicted energy flux densities are plotted over a longitudinal section located quarter-way the channel width h, while Fig. 6(A) and (C) correspond to those calculated at the longitudinal mid-plane, as shown in the inset figure. The figures illustrate the existence of large-areas of recirculating mean energy fluxes, intensity of which increases as the drive ratio grows from 0.75% to 2.0%, as first noticed by comparing Fig. 6(A) and (C). It is also observed that the axial extension of the recirculating zones closely corresponds to the cross-sectional mean particle dis^ x. Remarkably, such recirculating zones become placement U= more intense while moving from the mid-plane behavior to the 0:25h-corresponding section, such that longer recirculating pathlines are exhibited before getting back into the stack. In addition to that, limited asymmetrical behavior in the energy flux density is
captured at higher drive ratio at both longitudinal sections with more intensive appearance towards the center. Furthermore, the energy fluxes predicted at the channel midplane shown in Fig. 6(A) and (C) are found to be more concentrated rather than detached, with respect to those illustrated in Fig. 6(B) and (D). This is in part consistent with -and is possibly caused by- the induced non-zero DC-streaming and its associated vortical motion at the stack extremities with similar dependence on the imposed drive ratio, as illustrated in the next section. This might in turn assume a significant contribution to the overall thermal energy losses. Interestingly, this is not supported by the above agreement of the mean temperature difference and axial velocity profiles with theoretical predictions at DR ¼ 2:0%. Existing 2-D simulations [3] on parallel-plate stacks reported that the departure between the theory and the computations, in their evaluation for the steady-state stack temperature difference at DR 1:0%, results from the growing mean energy paths and the generated 2-D motion within the gap. In contrary with that, the present 3-D simulation clearly indicates the suppression of the multi-dimensional flow behavior within the square pores as demonstrated through the agreement of the axial velocity profiles with theory at the stack center. The other factor, although evident in the present work, its contribution does not strongly impact the bulk oscillating-flow behavior and thus does not cause significant departure from the linear theory in the studied range of drive ratio. It is therefore concluded that the effect of the eliminated multi-dimensional flow behavior within the stack square pores on the adherence to the 1-D linear theory is more important than
298
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
Fig. 6. Vector plots superimposing streamlines (Solid red lines) of the mean energy flux density at drive ratios of 0.75% (A and B) and 2.0% (C and D). These plots are presented for two different longitudinal sections parallel to (and 0:5h- and 0:25h-apart from) the channel walls, as illustrated for instance in the insets. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
alleviating the increasingly-diffusing large-areas of recirculating mean energy fluxes outside the stack pores. 4.3. Flow streaming In addition to the estimation of the mean energy-flux density, the ANSYS solver enables us to calculate the mean velocity components and flow vorticity to gain insight into the influence of the drive ratio on the induced non-linear flow streaming particularly at the stack extremities. Fig. 7 shows the contour plots of the calculated vorticity field superimposed by vector plots of its respective flow streaming at drive ratios 0.25% and 2.0%, as similarly predicted at the two longitudinal sections located parallel to (and 0:5h- and 0:25h-apart from) the channel walls. The vectors lengths are automatically scaled in the MatLab program for better visualizing the mean flow fields, whereas for instance, the results shown in Fig. 7(A) and (C) indicate a maximum vector length of 0.11 and 2.46 m/s at drive ratios 0.25% and 2.0%, respectively. The right plots of Fig. 7 also shows the calculation of the instantaneous velocity profiles at the stack left extremity to help understanding the additional entrance effect. Flatter velocity profiles are captured during the second half of the acoustic cycle while the gas particles are moving back into the stack channels because of the
non-fully-developed flow behavior at the stack entrance. Such effect is noticeably stronger at the larger drive ratio 2.0%. The streaming results, as illustrated in Fig. 7, also reveal the structure of the evolving three-dimensional vortical motions and their associated flow rotationalities at the corners of each square pore because of the sudden expansion and contraction the flow exhibits as getting in and out of the stack square pores. Concentrated vortices are found when operating at the lower drive ratio, whereas intense vortex shedding is shown at the larger drive ratio. Similar to the energy flux density, the predicted vortical pathlines are shown elongated and detached from the stack extremities at the larger drive ratio. Therefore, in the present particular 3D simulation of the thermoacoustic fields within a square stack, it is concluded that stronger currents of mean swirling motions, and thus vorticity, are captured at larger drive ratios, however, such non-linear effects are limited in terms of vorticity transport- and remarkably localized nearby the stack extremities. No evidence of significant streaming activity has been found at the stack middle section at both drive ratios that could drastically impact the one-dimensionality oscillating-flow condition within the stack; and thus, help reduce the non-linear flow losses and retain efficient transfer of heat within the square pores from the stack cold to hot ends.
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
299
Fig. 7. The left figures show vector plots of the flow streaming superimposing contours of the ‘mean’ vorticity field at drive ratios of 0.25% (A and B) and 2.0% (C and D). These plots are presented for two different longitudinal sections parallel to (and 0:5h- and 0:25h-apart from) the channel walls, following exactly the same arrangement shown in Fig. 6. The arrows lengths are adjusted automatically within the MatLab program for the sake of visualization. The vorticity is expressed in the unit of s1. The right figures plot the velocity profiles at eight different instances over a cross-sectional plane located exactly at the stack cold end.
5. Summary A new three-dimensional finite-volume model is presented that solves the full non-linear thermo-hydrodynamic laminar flow equations in three-dimensions and examines the heat transport between the thermoacoustically oscillating gas and the surrounding surfaces in a thermoacoustic stack made of square pores. The present model enables us to investigate the temperature variation with time at the stack ends along with the evolving axial velocity profiles at the stack center and the mean (time-averaged) energy flux density nearby the stack ends. The flow streaming and mean vorticity field are also calculated and presented in this study to gain insight into the role of the developed acoustic non-linearties
in determining the departure of the stack temperature difference from the linear theory. Here, the theoretical mean temperature difference is deduced based on the Arnott’s extended linear theory for stacks having arbitrarily shaped pore cross sections. In the present simulation runs, the captured instantaneous velocity profiles and the variation with the drive ratio of the estimated mean temperature difference DT m between the stack ends show good agreement with the Arnott’s extended linear theory [18] up to a drive ratio of 2.0%. In contrary with existing 2-D simulations carried out on parallel-plate stack, the present results indicates the prevalence of the 1-D flow within the stack pores and limits the multi-dimensional flow behavior to the stack extremities, such that the departure of the present 3-D computations from the
300
A.I. Abd El-Rahman et al. / International Journal of Heat and Mass Transfer 108 (2017) 292–300
theory has been effectively shifted away. Although the swirling pathlines of the mean energy flux densities are clearly evident, and even become more intense with the increase of the drive ratio, in the neighborhood of the stack ends, it is found in the present work that its influence is weak enough to promote departure between the theory and the computations. Therefore, the present work suggests that more useful transfer of heat on and off the stack is expected with the square-pore geometry than the corresponding parallel-plate structure. This study will be extended in our future work to include higher drive ratios, in which the potential of flow transition from laminar to turbulence is no longer trivial.
References [1] N. Cao, J.R. Olson, G.W. Swift, S. Chen, Energy flux density in a thermoacoustic couple, J. Acoust. Soc. Am. 99 (6) (1996) 3456–3464. [2] A.S. Worlikar, O.M. Knio, Numerical simulation of a thermoacoustic refrigerator: I. Unsteady adiabatic flow around the stack, J. Comput. Phys. 127 (2) (1996) 424–451. [3] A.S. Worlikar, O.M. Knio, R. Klein, Numerical simulation of a thermoacoustic refrigerator: II. Stratified flow around the stack, J. Comput. Phys. 144 (2) (1998) 299–324. [4] A.I. Abd El-Rahman, E. Abdel-Rahman, Computational fluid dynamics simulation of a thermoacoustic refrigerator, J. Thermophys. Heat Transf. 28 (1) (2014) 78–86. [5] P. Lotton, P. Blanc-Benon, M. Bruneau, V. Gusev, S. Duffourd, M. Mironov, G. Poignand, Transient temperature profile inside thermoacoustic refrigerators, Int. J. Heat Mass Transf. 52 (21–22) (2009) 4986–4996.
[6] P. Blanc-Benon, E. Besnoin, O. Knio, Experimental and computational visualization of the flow field in a thermoacoustic stack, CR Mec. 331 (1) (2003) 17–24. [7] H. Ishikawa, D.J. Mee, Numerical investigations of flow and energy fields near a thermoacoustic couple, J. Acoust. Soc. Am. 111 (2) (2002) 831–839. [8] D. Marx, P. Blanc-Benon, Computation of the mean velocity field above a stack plate in a thermoacoustic refrigerator, CR Mec. 332 (11) (2004) 867–874. [9] D. Marx, P. Blanc-Benon, Computation of the temperature distortion in the stack of a standing-wave thermoacoustic refrigerator, J. Acoust. Soc. Am. 118 (5) (2005) 2993–2999. [10] L. Zoontjens, C. Howard, A. Zander, B. Cazzolato, Numerical study of flow and energy fields in thermoacoustic couples of non-zero thickness, Int. J. Therm. Sci. 48 (4) (2009) 733–746. [11] S.H. Tasnim, R.A. Fraser, Computation of the flow and thermal fields in a thermoacoustic refrigerator, Int. Commun. Heat Mass Transf. 37 (7) (2010) 748–755. [12] G. Poignand, B. Lihoreau, P. Lotton, E. Gaviot, M. Bruneau, V. Gusev, Optimal acoustic fields in compact thermoacoustic refrigerators, Appl. Acoust. 68 (6) (2007) 642–659. [13] C. Fan, B.-T. Chao, Unsteady, laminar, incompressible flow through rectangular ducts, Z. Angew. Math. Phys. ZAMP 16 (3) (1965) 351–360. [14] M. Hino, M. Kashiwayanagi, A. Nakayama, T. Hara, Experiments on the turbulence statistics and the structure of a reciprocating oscillatory flow, J. Fluid Mech. 131 (1983) 363–400. [15] ANSYS Academic Research, Release 15.0, Help System, User Guide, ANSYS, Inc. [16] G.W. Swift, Thermoacoustic engines, J. Acoust. Soc. Am. 84 (4) (1988) 1145– 1180. [17] J. Wheatley, T. Hofler, G.W. Swift, A. Migliori, An intrinsically irreversible thermoacoustic heat engine, J. Acoust. Soc. Am. 74 (1) (1983) 153–170. [18] W.P. Arnott, H.E. Bass, R. Raspet, General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross sections, J. Acoust. Soc. Am. 90 (6) (1991) 3228–3237. [19] O. Zikanov, Essential Computational Fluid Dynamics, Wiley, 2010.