Numerical investigation of heated gas flow in a thermoacoustic device

Research PaperNumerical investigation of heated gas flow in a thermoacousticdeviceKazuto Kuzuu ⇑,Shinya HasegawaDepartment of Prime Mover Engineering,Tokai University,Hiratsuka,Kanagawa 2591292,Japanh i g h l i g h t sThe self-sustained oscillation in thermoacoustic engine is reproduced by CFD. Flow field in CFD is assigned to thermoacoustic properties from linear theory. Unique flow behavior due to non-uniform temperature gradient is observed in CFD.a r t i c l e i n f o Article history:Received 10March 2016Revised 1July 2016Accepted 16August 2016Available online 18August 2016Keywords:Thermoacoustics Heat exchangerComputational fluid dynamics Oscillatory flowa b s t r a c tSelf-sustained oscillation in a standing wave thermoacoustic device is reproduced via computational fluid dynamics simulations,and heated flow behavior in the device is explored using the results obtained.The straight-type thermoacoustic device is composed of two resonance tubes,two heat exchangers and a stack.In these simulations,both the acoustic characteristics and the temperature field during self-sustained oscillation are considered.Therefore,to reproduce the self-sustained oscillatory flow,an acous-tic signal is injected into the computational domain as a trigger pulse.From the results,oscillatory fluid motion around the engine is investigated,and characteristics,including the acoustic field or work flow,the energy dissipation,the work source and other associated aspects,are estimated.The results agree well with those of linear theory,although high energy dissipation caused by vortex generation is observed near the engine.This result verifies the computational fluid dynamics simulation results.The temperature field around the engine is also investigated.The results show occurrence of asymmetrical temperature oscillations within the heat exchangers.This behavior cannot be predicted using linear the-ory because the non-uniform temperature gradient in the engine unit is transferred in stream-wise by convection.Finally,a modification to conventional linear theory is suggested to reproduce this behavior.Ó2016Elsevier Ltd.All rights reserved.1.IntroductionThe thermoacoustic device reviewed by Swift in 1988combines the advantages of the inherent thermal efficiency properties of the Stirling cycle with the ability to work using a minimum number of moving parts [1].Swift et al.demonstrated the feasibility of the device using an inexpensive prototype [2],and since then,ther-moacoustic devices have attracted considerable attention as appli-cations of renewable heat energy.Many basic and practical studies of these devices have been performed in recent years.Thermoacoustic devices work using thermoacoustic phenom-ena.Theories to explain these phenomena were originally formed based on acoustic theory [3–7],and the linear theories developed by Rott [8,9]and Tijdeman [10]are still commonly applied to developments in these devices.For example,the Design Environ-ment for Low-amplitude Thermoacoustic Engines (DELTAE),which is a numerical analysis code developed by Swift et al.[11],uses the linear theory of Rott [8]and is useful for thermoacoustic device design.Linear theory is thus significant when dealing with thermoa-coustic phenomena that occur in a thermoacoustic device.However,some phenomena cannot currently be explained using linear theory.It is not actually possible to predict non-linear phe-nomena like two or three-dimensional vortex generation.Heat transfer in the heat exchanger is another problem to be solved in this field,and the problem must be studied using methods beyond linear theory.Many studies of heat transfer in thermoacoustic devices have been performed recently,and they can be classified into several different types of approaches.The first approach combines linear theory and numerical calcu-lations.Piccolo et al.introduced a simple energy conservation model coupled with classical linear thermoacoustic theory,and/10.1016/j.applthermaleng.2016.08.0931359-4311/Ó2016Elsevier Ltd.All rights reserved.⇑Corresponding author.E-mail addresses:kuzuu@tokai-u.jp(K.Kuzuu),s.hasegawa@tokai-u.jp(S.Hasegawa).used it to estimate the heat transfer properties of thermoacoustic heat exchangers composed of several parallel plates[12,13].In another method,de Jong et al.proposed a heat transfer model for one-dimensional oscillatoryflow.This model can be applied to parallel-plate thermoacoustic heat exchangers,and was used to investigate their heat transfer properties[14].Experimental approaches can also be used to investigate heat transfer.Wakeland and Keolian et al.investigated the heat transfer of parallel-plate heat exchangers in oscillatoryflow environments [15].They estimated heat exchanger effectiveness based on its temperaturefield,which was measured using thermistor probes placed at the heat exchanger exits and entrances.They also com-pared the results of their study with those obtained using the DEL-TAE code[11].Additionally,Jaworski et al.investigated the heat transfer properties of parallel-plate heat exchangers through acetone-based planar laser-inducedfluorescence(PLIF)measure-ments[16,17],and compared their measured results with numer-ical results[18,19].Thefinal approach is based on computationalfluid dynamics (CFD).This approach is advantageous because the assumptions in linear theory are excluded,and because the heat transfer proper-ties of a thermoacoustic device can be calculated directly from its temperaturefield.Some numerical simulations of oscillatoryflow in heated pipes were performed by Zhao and Cheng[20,21].While their simulations were for oscillatoryflow induced by acoustic wave propagation,their results provide the characteristic features of an oscillatory temperaturefield in the tube.In another numeri-cal simulation,Cao et al.investigated the energyflux density in a thermoacoustic couple under acoustic standing wave conditions [22].In this study,they estimated the effects of the displacement amplitude on heat transfer.The study is significant for heat trans-fer estimation because the displacement amplitude in Rott’s theory [8,9]is assumed to be negligible when compared with the device length.Ishikawa and Mee also studied theflowfields and energy transport near thermoacoustic couples through numerical simula-tions using a2D full Navier–Stokes solver[23].Marx and Blanc-Benon performed numerical simulations of a thermoacoustic refrigerator that consisted of a resonator and a parallel plate stack [24].They compared their results with those predicted by linear theory,and showed that there is a difference in mean temperature between thefluid and the plate.Mohd and Jaworski also investi-gated the oscillatoryflow and heat transfer of parallel heat exchangers,and compared their results with experimental data [25].Using both numerical results and experimental data,they demonstrated the effect of the temperaturefield on oscillatory flow and the dependencies of heat transfer on the Reynolds number.Additionally,self-sustained oscillation is also reproduced and discussed with respect to CFD simulations of thermoacoustic phe-nomena.Hantschk and Vortmeyer simulated self-sustained oscilla-tion in a Rijke tube[26].While the simulated tube includes only heating elements,rather than heat exchangers,they discussed non-linearity in the heat transfer process.Recently,other CFD sim-ulations of self-sustained oscillation in thermoacoustic engines have been performed.Spoelstra et al.simulated a travelling-wave thermoacoustic engine,and showed strong non-linear effects for high-amplitude thermoacoustic systems[27].Zink et al.showed the transition from initial disturbance to self-sustained oscillation in a thermoacoustic engine,and explored the effects of a curved resonator[28].Dai et al.simulated self-sustained oscillation in a 300Hz standing wave thermoacoustic engine[29]by visualizing flowfields at the ends of the stack,and discussed multi-dimensional effects that occurred because of the abrupt change inflow area.As described above,many studies of thermoacoustic devices have been performed using a variety of approaches.CFD simula-tions attract particular attention because of their ability to repro-duce non-linear effects in thermoacoustic phenomena,and the technique has progressed such that the non-linearity offlow behavior in self-sustained oscillation can be discussed.However, while the acousticfields obtained from these simulations are com-parable with those obtained using conventional linear theory,the non-linearity of heatedflow behavior,which must affect the sys-tem heat transfer,has not been discussed adequately to date.In this study,to investigate such heatedflow behavior occurring in a thermoacoustic device in greater detail,unsteady CFD simula-tions of self-sustained oscillatoryflow are performed.For this pur-pose,a full-scale device model was set up and the engine unit was given a suitable temperature profiing this setup,self-sustained thermoacoustic oscillation can be reproducedNomenclatureq densityq m mean(time-averaged)densityu iflow velocity vectorh u i r velocity amplitude at cross-sectionp absolute pressurep m mean(time-averaged)pressurep1pressure amplitude at cross-sectionu phase of acoustic waveu pu phase difference between pressure and velocity T absolute temperatureT m mean(time-averaged)temperatureT1temperature amplitude at cross-sectionT w wall temperatureT Hw,T Cw wall temperatures of HEX and CEX,respectively n i normal vector of surface elementt timeD t period of acoustic wavef frequency of acoustic waveV volume element of control volumeS surface element of control volumeA area of cross-section R gas constants ij shear stress tensorq j heatfluxl viscositym kinetic viscosityk thermal conductivitya thermal diffusivityc specific heat ratior Prandtl numberd ij Kronecker deltaj imaginary unitx angular frequencyv m,v a thermoacoustic functionI workflowW m,W p kinetic and potential energy dissipationW prog,W stand travelling and standing wave components of work source,respectivelye i internal energy/viscous dissipation functionC p specific heat at constant pressure1284K.Kuzuu,S.Hasegawa/Applied Thermal Engineering110(2017)1283–1293numerically.A linear analysis based on conventional thermoacous-tic theory is performed simultaneously,and the results are com-pared.In particular,rather than simply compare the results for the acoustic characteristics,the relationship between the acoustic characteristics when treated using linear theory and theflowfields obtained from CFD is clarified,and the linear and non-linear char-acteristics of oscillatoryflows are discussed.Finally,the tempera-turefield behavior that cannot be explained using linear theory is investigated.2.Calculations2.1.Calculation method in CFDCFD simulations are carried out using the LS-FLOW unstruc-tured compressibleflow solver developed by the Japan Aerospace eXploration Agency(JAXA)[30].The solver is based on three-dimensional unsteady compressible Navier–Stokes equations.The basic equations are as follows.@ @tVQd Vþt SðF eÀF vÞdS¼0ð1ÞQ¼qq u xq u yq u zE2666666437777775;F e¼q Uq u x Uþn x pq u y Uþn y pq u z Uþn z pðEþpÞU2666666437777775;F v¼s xx n xþs yx n xþs zx n xs xy n yþs yy n yþs zy n ys xz n zþs yz n zþs zz n zðu i s ijÀq jÞn j2666666437777775U¼u x n xþu y n yþu z n z;u i¼u xu yu z264375;sij¼À23lðrÁUÞd ijþl@u i@x jþ@u j@x iE¼1cÀ1pþ12q U2;q i¼Àk r T;where q,u i,p,T,l and k are the density,the velocity vector,the pressure,the temperature,the viscosity and the heat conductivity of the gas,respectively,and n i is the vector normal to the volume element surface.The main numerical schemes used in this simula-tion are given in Table1.2.2.CFD calculation modelFig.1shows the calculation model used in this simulation.The engine unit has three components:hot and cold heat exchangers (HEX and CEX),and a stack(STK).Each component contains sixflat plates and seven channels.The computational domain includes the buffer region.The individual part sizes and boundary conditions are shown in Fig.1(a).The model is two-dimensional;this is achieved by providing symmetrical conditions in the z direction. The wall temperatures,excluding those of the sixflat plates of the engine unit,are allfixed at room temperature(298.15K).The temperatures in the engine unit are423.15K for HEX and 298.15K for CEX.A linear distribution ranging between423.15K and298.15K is used for the STK temperature.The working gas is air,the working pressure is p=1.01325Â105Pa,and the temper-ature T=298.15K.In this simulation,to induce self-sustained oscillatoryflow,the pressure disturbance was injected as a trigger pulse from the open end.The simulations actually start from a static initial condition and continue for a short time.Then,one-half cycle of a sinusoidal acoustic wave is injected as a disturbance with^p%283Pa and f=21.2Hz from the open end into the buffer region.The supply of this acoustic wave is then terminated.The mesh configuration is composed of hexahedral cells that were converted from non-uniform Cartesian grids.The concen-trated mesh state of the representative parts is shown in Fig.1 (b).For two-dimensional calculations,the number of divisions in the z-direction is set at one.The total number of cells is approxi-mately300,000,and the minimum mesh size,which corresponds to the distance from the boundary wall to the adjacent mesh,is 0.0237mm.The time step used for this simulation is2.0l s,which corresponds to a Courant–Friedrichs–Lewy(CFL)number of30 based on the sound velocity.2.3.Numerical calculations based on linear theoryUsing Rott’s linear theory[8],the basic equations for one-dimensional thermoacoustic oscillatoryflow within a tube can be developed and can thus be expressed as simultaneous differential equations with respect to velocity and pressure.@p1@x¼Àjxqmð1Àv mÞh u i r;ð2Þ@h u ir@x¼Àjxpm1ÀcÀ1cð1ÀvmÞ!p1þv aÀv mð1Àv mÞð1ÀrÞ1T m@T m@xh u ir;ð3Þwhere v m and v a are the viscous and thermal thermoacoustic func-tions,respectively,and have two-dimensional formulations as follows.v m¼tanhðð1þjÞffiffiffiffiffiffiffiffiffix s mpÞð1þjÞffiffiffiffiffiffiffiffiffix s mp;v a¼tanhðð1þjÞffiffiffiffiffiffiffiffiffix s apÞð1þjÞffiffiffiffiffiffiffiffiffix s apð4Þs m¼r20=2m;s a¼r20=2aHere,r0is one-half of the two-dimensional channel width.These equations are calculated numerically.For example,Ueda et al.proposed a numerical method based on the fourth-order Runge–Kutta method and calculated the critical temperature ratio for self-sustained thermoacoustic oscillation[34].In the linear analysis,the pressure amplitude at the antinode point(closed end)is given as a boundary condition.A value of 495Pa is calculated using CFD.For the temperature conditions, the time and section-averaged temperature obtained from CFD is used.These conditions lead to the properties shown in Table2.3.Results and discussion3.1.Production of self-sustained oscillationTo confirm self-sustained oscillation,time variation of the phys-ical values after termination of the trigger injection is observed. Fig.2shows the variations in the axial velocity and pressure at (x,y)=(1.04,0.0).The acoustic signal injection time is t=2.054–2.078s.As shown in thefigure,the phenomenon is considered toTable1Numerical schemes used in CFD simulations.Calculation Name of schemeTime integration Three points backward stepapproximationImplicit solution LU-SGS[31]Interpolation MUSCL scheme by Green-Gauss methodNumericalflux ConvectivetermSLAU[32]Viscous term Wang’s method[33]K.Kuzuu,S.Hasegawa/Applied Thermal Engineering110(2017)1283–12931285be in a periodic steady state after4s and is then in self-sustained oscillation.The enlarged graph shows the period around the occurrence of the trigger pulse and the periodic steady state.3.2.Verification of CFDOscillatoryflow behavior in the device is investigated,and the CFD results are compared with those from linear theory.Here, we show the velocity profiles at two device cross-sections.The velocity profile of oscillatoryflow varies with time.Fig.3shows 12phases in a single cycle.In Fig.3,the velocity amplitude is the section-averaged value at the midpoint of the STK(x=1.045m), and the corresponding displacement amplitude is also shown.Fig.4shows the variations in the velocity profiles.Each graph corresponds to profiles on the section including the STK and the resonance tube(x=1.045and2.00m).Both profiles agree with the results of linear theory.This implies that the fundamental oscillatoryflow behavior obeys linear theory within the device. 3.3.Acoustic propertiesIn thermoacoustic theory,the acoustic power is expressed as a workflow.From the CFD results,we can calculate this value using Eq.(5).I¼x2pÂZ tþ2p=xtt AðpÀp mÞu x dAf g dtð5ÞIn Eq.(5),the workflow,I,is the time-averaged value over a sin-gle acoustic wave cycle.Fig.5shows a comparison of the CFD simulation results with those from linear analysis.Here,the temperature gradient used in the linear analysis is based on the CFD results.In thefigure, the CFD results almost agree with the power gain obtained in the linear analysis.However,the CFD power gain is approximately 9%smaller than that determined by linear analysis.This may be caused by differences between the CFD simulations and the linear analysis in terms of temperature conditions at each section.In fact,putationalTable2Properties of air within the engine.T(K)q(kg/m3)m(m2/s)s m x4020.8782 2.6299-05 5.617298.15 1.1847 1.5561-059.4931286K.Kuzuu,S.Hasegawa/Applied Thermal Engineering110(2017)1283–1293while the temperature is assumed to be common at each section of the engine channels in the linear analysis,the temperature differs in every channel in the CFD simulations.This difference is described in detail later.The enlarged graph shows the characteristics around the engine unit.3.4.Energy dissipation and work sourceIn this section,we calculate the energy dissipation and work source from the CFD results and compare them with the results from linear theory.First,we consider the energy dissipation and the work source as thermoacoustic properties.To explain the cause of the workflow, Tominaga[35]introduced the concept of energy dissipation and the work source based on linear theory offlow in a tube.In this theory,energy dissipation is classified in terms of the kinetic and potential energy dissipations,or W m and W p,respectively.Addition-ally,the work source is divided into two parts,i.e.,the travelling and standing wave components,W prog and W stand,respectively.These parameters are defined as follows.W m¼ÀA2Rej xq mð1Àv mÞ!jh u irj2ð6ÞW p¼ÀARej xc m1þðcÀ1Þv aÀÁ!j p1j2ð7ÞW prog¼AReðv aÀv mÞð1Àv mÞð1ÀrÞT m!j p1jjh u irj cos u pu@T m@xð8ÞW stand¼AImðv aÀv mÞð1Àv mÞð1ÀrÞT m!j p1jjh u irj sin u pu@T m@xð9ÞHere,we must extract the energy dissipation and the work source from the CFD results.K.Kuzuu,S.Hasegawa/Applied Thermal Engineering110(2017)1283–12931287Thefirst law of thermodynamics is commonly given asd Q¼de iþp d vð10Þd Q is the heat energy added to thefluid element from an exter-nal system,d e i is the increment in the internal energy,and p d v cor-responds to the work in which the element acts for the external system.Given that d Q is the heat energy added to thefluid element by heat conduction,the energy equation that applies to two-dimensionalfluid motion is obtained as follows.q De i¼À@@xÀk@T@xþ@@yÀk@T@y&'þðÀpÞ@u@xþ@v@yþs xx@u@xþs yx@u@yþs xy@v@x þs yy@v@yð11ÞIn the above equation,the third term corresponds to kinetic energy dissipation by viscosity,while the second term can be regarded as a summation of the potential energy dissipation and the work source.By integrating these terms with respect to both time and cross-section,the energy dissipation and the work source properties can be obtained as follows.Kinetic energy(viscous)dissipation:/¼s xx@u@xþs yx@u@yþs xy@v@xþs yy@v@y¼lÂ2@u@x2þ@v@y2!þ@v@xþ@u@y2À23@u@xþ@v@y2!ð12Þh/it¼x2pÂZ tþ2p=xtt A/dAf g dtð13ÞSummation of potential energy dissipation and work source: w¼ÀðpÀp aÞÂ@u@xþ@v@y¼ðpÀpaÞqÂD qð14Þh w it¼x2pÂZ tþ2p=xtt A w dAf g dtð15ÞFig.6shows the energy dissipations and the work source calcu-lated using the above equations.The work source in the resonance tube can be neglected because of the lack of a heat source.From 1288K.Kuzuu,S.Hasegawa/Applied Thermal Engineering110(2017)1283–1293this,it can be said that both the viscous and potential energy dis-sipations almost follow linear theory in the resonance tube.Additionally,we show the kinetic energy dissipation within and near the engine.In CFD,this dissipation is calculated using Eq.(12).The enlarged graph in Fig.6shows a comparison of the results of CFD and linear analysis.The graph shows good agreement,but dif-ferences can be observed near the heat exchanger extremities (x =1.00and 1.09m).This may be the effect of vortex generation around the corners of the engine plates.However,the work flow is not so strongly affected by this difference because the total energy loss is decided by the integration in the x-direction and the peak value is limited to quite a narrow area.Next,the effects of potential energy dissipation and the work source around the engine unit are investigated.In CFD,this effect is described using Eq.(14),and is the summation of the potential energy dissipation and the work source.Therefore,the analytical results for comparison must be expressed using the summation of W p ,W prog and W stand from Eqs.(7)–(9).Fig.7compares theK.Kuzuu,S.Hasegawa /Applied Thermal Engineering 110(2017)1283–12931289results from CFD and linear theory.They show almost perfect agreement.Compared values:h w i t in CFD and W p+W prog+W stand in the analysis.From thermoacoustic theory in a tube,the workflow must be consistent with the integration of the above values,as follows:I¼Z xW mþW pþW progþW standÀÁdxð16ÞFig.8compares the results of the integration in Eq.(16)above with the workflow calculated from Eq.(5),and shows good agree-ment between these results.The work source and the dissipation are thus considered to have been correctly estimated using the CFDflow data.3.5.Flow behavior in a temperaturefieldFinally,the temperaturefield around the engine unit is investigated.In the linear theory of Rott[8],it is assumed that the time-averaged temperature of thefluid in each section coincides with the wall temperature.However,because of heat transfer between thefluid and the solid wall,the averagedfluid temperature may possibly differ from the wall temperature.In fact,while the wall temperature of the HEX(423.15K)is given as shown in Fig.1, thefluid temperature in the CFD simulation is lower than the wall temperature.Fig.9shows the time and section-averaged tempera-ture distributions in the x-direction around the engine.White cir-cles denote the results for each section along the device,and in the engine region(x=1.00–1.09m),the temperature distribution for each channel is also shown in the graph.Each result corresponds to one of the center(black circles),bottom(white triangles)or top(black triangles)channels in the engine.The wall temperatures of the engine plates are plotted using a dotted line.There is a dif-ference between the temperature of thefluid and that of the engine plates;specifically,thefluid temperatures in the bottom and top channels are much lower than that of the engine plates.This is because one side of each of these channels is set at room1290K.Kuzuu,S.Hasegawa/Applied Thermal Engineering110(2017)1283–1293temperature (298.15K).As a result,the effective temperature,when averaged over all engine channel sections (white circles),is further reduced when compared with the engine plate temperatures.On the other hand,a high-temperature region on the outside of the HEX (x =0.95–1.00m)can also be observed.This is caused by leakage of heated air from the HEX;the area is called a thermal buffer region.The rapid increase in temperature at x =1.0m can be considered to be affected by the wall temperature at the engine plate extremity.The behavior of the heated air is visualized in Fig.10.In this figure,the high temperature region of the thermal buffer is asymmetrical around the central line and deflects towards the upper wall.This is because the temperature is averaged over one period of an acoustic wave at most,despite the flow being dis-turbed by vortex generation.Next,to investigate the temperature field behavior in detail,the phase variations of the temperature are observed.The observed area is the central channel of the engine.The CFD results are also compared with the values that were predicted using Eq.(17)based on linear theory.T 1¼1C p q m ð1Àf a Þp 1þf a Àf m 1Àv m Þð1Àr ÞÀ1Àf m 1Àv mh u i r 1j x @T m@x ð17ÞIn Eq.(17),the temperature gradient in the x-direction is taken from the temperature distribution on the center channel (black cir-cles)shown in Fig.9.Fig.11shows the distribution within the central section of STK.The normalized temperature H must now be introduced:H ¼T ÀT wT Hw ÀT Cw;ð18Þwhere T w ,T Hw ,and T Cw are the wall temperatures of the local sec-tion,the HEX,and the CEX,respectively.As the graph shows,the CFD results closely follow those of the linear theory.This is because the temperature gradient in the x-direction is almost uniform in this region.In contrast,at the boundary area between the heat exchanger (HEX and CEX)and the stack (STK),the situation is quite different.As shown in Fig.12,the temperature oscillation from CFD does notagree with the results predicted using Eq.(17).The actual temper-ature oscillation in the CFD results is asymmetrical.In Eq.(17),the temperature gradient is assumed to be constant within the dis-placement amplitude of the local fluid element.In this case,the actual fluid element moves to an extent of 15mm,as shown in Fig.3.When the temperature distribution shown in Fig.9is consid-ered,the above assumption cannot be applied in this area.This means that the temperature gradient varies in the oscillatory cycle.To improve this point,the convective term of Eq.(17)must therefore be modified.Here,it is assumed that the mean tempera-ture gradient can be described using the linear transfer equation@F @t þu x @F @x¼0;F ðx ;t Þ¼@T m@x:ð19ÞThe solution to the above equation is:F ðx ;t Þ¼F ðx Àu x t Þ:ð20ÞHere,u x t corresponds to the fluid element displacement in the oscillatoryflow.Fig.10.Temperature distributions around the engine.-1.5-1-0.50.511.5-0.4-0.3-0.2-0.1 0y (m m )ΘΦ01Φ03Φ05Φ07Φ09Φ11parison of temperature distributions linear analysis).。