The Piston Effect inside a microchannel with Carbon Dioxide near Critical conditions


 Heat transfer near the critical condition of Carbon Dioxide due to thermo-acoustic waves in a 100-µm high microchannel was numerically studied. The temperature at a point farthest away from the heated surface was compared between computational fluid dynamics (CFD) models and a pure conduction model. The comparison revealed that the CFD model predicted a temperature increase furthest from the surface much faster than the time constant required for such increase purely by conduction. It is believed that another heat transfer process, termed the piston effect (PE), which is associated with pressure waves in the fluid, was responsible for this increase. Explicit unsteady methodology in the fluid model indicated that propagation of pressure waves due to a rapid expansion of the boundary layer and the associate change in the fluid density distribution resulted in this temperature raise. It was confirmed that natural convection wasn’t responsible for the temperature increase under quiescent conditions. In addition, it was discovered that the PE is significant for certain forced convection conditions.


Introduction
Introductory textbooks typically cover three established heat transfer modesconduction, convection, and radiation. However, in the early 1970's it was found that around the critical condition of fluids, a fourth heat transfer mode exists, termed the "piston effect (PE)," which is activated by thermo-acoustic waves due to the fluid's compressibility [1,2]. When heat is suddenly generated on a surface, the fluid gradually heats up in the vicinity of the heated surface due to conduction heat transfer. Since fluids near the critical conditions exhibit large density swings with temperature, the density inside the boundary layer rapidly decreases causing rapid expansion of the fluid near the wall, which like a 'piston' compresses the fluid outside the boundary layer. The mechanical energy generated by the high-pressure waves of the PE is then converted back to heat that is dissipated into the bulk fluid. Thus, a two-prong mechanism to transfer heat from the wall to the bulk fluid is established, which bypasses the more common convective heat transfer process.
Initially, the PE was observed in microgravity [3], and later it was also observed in natural convection under terrestrial conditions [4] and even in boiling [5,6,7]. For instance, microgravity experiments in the MIR space station suggested that the PE effect can also be present during boiling in fluids such as CO2 and SF6 [7], and that overheating of vapor bubbles is attributed to an adiabatic compression by the rapid expansion of the fluid near a hot surface. Soboleva [8] simulated a two-dimensional square cavity with a heater at the bottom using the Van der Waals equation of state, and based on the rapid thermalization of the bulk fluid, it was concluded that the PE was significant even in the presence of gravity. Experimental [7], analytical [9], and numerical [10] studies suggested that this process is typically orders of magnitude faster than the conduction and/or advection heat transfer processes. For instance, Zappoli et al. [10] found that the thermo-mechanical couplings are responsible for an equilibrium state with a time scale that is 1% of the diffusion time scale. Similarly, Nitsche and Straub [11] experimentally discovered that under microgravity condition, thermal equilibration of fluids occurred in a much faster relaxation times than for conduction, which was attributed to the PE.
While previous studies revealed the importance of the PE under certain conditions, it was never examined at the micro scale where natural convection is expected to diminish and time scales are much faster than at the conventional scale.
This article attempts to address this shortcoming by reporting on a numerical study preformed to examine the piston effect inside a microchannel near the critical condition of CO2. Three types of heat transfer models were developed and compared to quantify the relative importance of the PE-conduction, natural convection, and forced convection under micro gravity and terrestrial conditions. (While the Rayleigh number at the micro scale is low and natural convection is expected to be suppressed, it was examined to eliminate potential contributing factors to the bulk fluid temperature changes that were detected by the simulation.) The numerical results suggest that the PE is significant even under certain forced convection conditions, and the time scales associated with adiabatic thermalization (caused by thermo-acoustic waves) are 5 to 6400 times faster than the diffusion time scales. Starccm+, a commercially available CFD solver, which was successfully implemented by other researchers [12,13] for the flow of supercritical Carbon Dioxide inside microchannels, was chosen to model PE in this study.

Results and discussion
The 2-D domain used to study the PE at the micro scale is shown in Fig. 1 consisting of a rectangular microchannel with a depth of 100 µm, a length of 7.5 mm, and a heater element with a length of 1.5 mm placed at the center of the top surface.
A total of four cases were analyzed including conduction only, quiescent flow (i.e., no forced convection) in microgravity and terrestrial conditions, and forced convection with micro gravity. These conditions were compared among each other to quantify the PE. Properties, such as pressure, density, and temperature were predicted from the simulation at a location that was 91 microns vertically down from the heated surface as shown in Fig. 1. (a) Case 1: Conduction only To quantify the PE, a baseline case was established that assumes only conduction in the CO2 (i.e., only the energy equation was solved with no flow). This allowed to decouple the PE from the conduction heat transfer process. It should be noted that because the Rayleigh number is zero for zero-g, natural convection is completely suppressed (i.e., no advection due to buoyancy effects), and absent PE, heat transfer occurs only by conduction. Thus, when no flow is allowed in or out of the domain, but advection is allowed by the simulation (i.e., the domain is treated as a fluid with pertinent properties of the CO2), it will be a result of the PE, and any deviation from the conduction simulation can be attributed to thermo-acoustic waves. An unsteady simulation was carried out at a time step of 1 nano second (till the simulation time of 30ms and a time step of 10 nano seconds afterwards) to predict the conduction time scale. The entire domain was initially kept at a temperature of 296.46 K and a constant heat flux of 6.75 W/cm 2 was supplied to the heater. Fig. 2 shows the temperature of the probe w.r.t time. The first temperature change was detected at the probe position after 19 ms. The temperature gradually increased and reached a value of 297.524 K after 0.1 s. Figure 2. Temperature as a function of time (b) Case 2: Quiescent medium in microgravity Explicit unsteady method with a Courant number of 0.5 was used, resulting in a time step as low as 3.0 ns. This time scale was required due to the small dimensions of the channel and the timescale associated with thermo-acoustic waves. Once the acoustic/pressure wave was generated, it propagated in all directions and reflected back from the adiabatic walls. These acoustic waves carry energy that raised the temperature of the CO2 at the probe position. The initial temperature of the fluid was 296.46 K, the pressure was 7.308675 MPa and, the heater heat flux was 6.75 W/cm 2 .
The simulation revealed that the first pressure disturbances felt by the probe was at a time of 0.347 µs. Fig. 3 shows the pressure and transverse velocity distributions inside the microchannel 2.73 µs after the heater was turned on. A pressure wave originated from the heater surface, propagated throughout the channel, and triggered fluid motion (i.e., advection) in the stagnant fluid even in the absence of gravity. The fluid adjacent to the heater had a velocity of 225 μm/s, vertically downwards, while the rest of the fluid was almost stagnant. The fluid layer next to the heater surface with higher velocity had a thickness of approximately 3 microns. As the fluid moves away from the heater (vertically downwards), it hits the adiabatic wall and appears to be deflected sideways, forming a very thin boundary layer. In all contour plots the probe is shown as a red dot. A new pressure wave that was generated by the heater and propagated throughout the micro channel is converted back to heat increasing the temperature of the fluid. The temperature detected by the probe at this time was 1.3 mK higher than the initial temperature and the entire cavity (away from heater) witnessed the same temperature raise (i.e., a raise of 1.3 mK) as the acoustic waves move through the cavity. The acoustic waves and the associated expansion of the fluid near the heater produced fluid motion in the quiescent fluid around 262 μm/s in the direction perpendicular to the heater.  Initially, each unique acoustic wave resulted in a fluid temperature increase of ~0.7 mK. However, with time, the temperature of the fluid rose rapidly as the acoustic waves continued to pressurize the entire cavity (e.g., after 1 ms the temperature increased by 26.52 mK). Initially, each acoustic wave appeared to carry a strength of 560 Pa, and after a time of 5 ms, the temperature of the entire cavity (away from the heater) rose by around 140 mK. (The conduction model predicted a temperature increase only after about 19 ms).
The rapid decrease in the fluid's density due to the temperature increase near the heated wall resulted in rapid compression of the remaining fluid in the region outside the developing thermal boundary layer. Although the drop in density of the fluid near the heater was rapid, the raise in density of the fluid away from the heater was less because the region occupied by the compressed fluid was significantly larger compared to that of expanding fluid. After 14 ms, the density near the heater dropped by 321 kg/m 3 while away from the heater the density increased by around 18 kg/m 3 . Fig. 6 shows the transverse velocity distribution inside the microchannel at different times. Velocity plumes are originating from the heater surface at around 6.85 ms. This observation was consistent with the observations of other studies, such as Soboleva [8] who observed a velocity plume that started from the heater surface in the absence of gravity.
As time progressed, the velocity plumes disappeared, energizing the adjacent fluid. It was also observed that the plumes had a maximum transverse velocity of 529 μm/s around 7.93 ms. At about 8.04 ms, almost all the velocity plumes disappeared, heating up the entire fluid near to the heater surface and supporting the fluid heating due to the PE. It took only 1.19 ms for these velocity plumes to disappear and energize the adjacent fluid layers. Animations (not showed here) revealed that this fluid layer oscillated like an automobile piston, transferring the heat from the heated surface to the cold quiescent fluid.In this period of 1.19 ms the temperature of the bulk fluid (away from the heater) rose by ~30 mK, from 296.65 K to 296.68 K. The body force term with a gravitational acceleration of 1-g was considered while solving the Navier-Stokes equations in terrestrial conditions. All the other boundary conditions, including the heater heat flux, were kept the same as for the zero-gravity case. The transient temperature of the probe position was compared to the zero-gravity condition (i.e., Case 2), and was concluded that the gravitational effects are negligible leading to conclude that the PE is significant even in terrestrial conditions.
(d) Case 4: Forced convection A total of four cases with different mass fluxes were examined for PE under forced convective conditions. Initially a steady state solution was obtained for the cold flow, for each mass flux, before turning on the heater to ensure that the flow is hydrodynamically fully developed. Once the steady state flow condition was achieved, the heater was turned on (with constant surface heat flux of 6.75 W/cm 2 ) and the solution was allowed to run transiently. The temperature raise of the fluid at the probe location was predicted and compared with the temperature raise of the fluid in the quiescent case (Fig. 7) w.r.t to time. Tmax is the maximum temperature the probe detected, and Ti is the initial fluid temperature, i.e., 296.46 K. The probe temperature peaked before falling back to the free stream temperature. The initial increase in the probe temperature attributed to the PE that took place immediately after turning on the heater. The probe temperature increased with decreasing mass flux suggesting that the PE dominance increases with diminishing mass flux and takes more time to reach back to the free stream temperature.
Normalized temperature is plotted against the Fourier number (ratio of heat conducted to heat stored) in Fig.8. Where Ө = [T(t)-T(t=0)]/[Tmax-T(t=0)]. Tmax is the maximum temperature attained by the probe and T(t=0) is the initial temperature of the fluid, i.e., 296.46 K. As the Fourier number is low, the heat transfer due to diffusion is less than the heat stored inside the fluid. At these conditions the heat transfer to the fluid due to PE is significant.

Validation
The results of the simulations in this study were validated using studies conducted by Hegg [3]. The density change with respect to the base density at the probe location as a function of time is shown in Fig.10. The thermo-acoustic waves that were generated by the heater reached the probe location at acoustic time scales. The expanding thermal boundary layer adjacent to the heater surface compressed the bulk fluid, and as a result increased the density of the bulk fluid. This pattern is consistent with the that reported by Hegg [3].

Methods
The computational domain was constructed in Starccm+ with a structured mesh of 3-µm (cell length). The total number of elements in the domain were 85,000. The top and bottom sides, except for the heater, were modeled as adiabatic walls; the left and right sides were modeled as symmetry planes for the quiescent flow simulations and mass flow inlet and outlet for the forced convective flow simulations.
For the Case 1, only energy equation was solved with no flow. Implicit unsteady methodology was employed, and the 2 nd order temporal discretization and 2 nd order spatial discretization schemes were selected to solve the unsteady energy equation. The Gauss-Seidel relaxation scheme(iterative method) was used to solve the governing equation.
In Case 2, an explicit-unsteady methodology(Runge-kutta scheme) was used solve the unsteady flow and energy equations in a coupled fashion with a Courant number of 0.5 to achieve a time step size of 3 ns. The body force term in the Navier-Stokes equations was not considered in the solution. Residual smoothing and preconditioning of equations were disabled to preserve the temporal accuracy of the unsteady solution. Upwind based AUSM+ fluxvector splitting scheme was employed to evaluate the inviscid fluxes. 2 nd order temporal and spatial discretization schemes were employed to solve the unsteady flow and energy equations. In Case 3, the same methodology used as in the Case 2, but the body force term was included in the Navier-Stokes equations. The reference gravity of -9.81 m/s 2 and the reference density of 721.47 kg/m 3 were used to predict the gravitational effects.
In Case 4, first the quasi-steady state solution was obtained using implicit formulation for the cold flow, for each mass flux, before the heater was turned on. The Courant number was ramped linearly from 0.1 to 200 to achieve faster convergence. The Gauss-Seidel relaxation scheme was used to solve the governing equations in the steady state solution, and a second order convergence was achieved. After the heater was turned on, the solution was allowed to run transiently and the same methodology as in the Case 2 was employed, but with a Courant number of 0.95 to achieve faster time marching.
In all the above Case 1 to Case 4, Venkatakrishnan gradient limiters [14], flow boundary diffusion, and secondary gradients were turned on and laminar flow settings were considered. Also, in all cases the initial temperature of the fluid was 296.46 K and the initial pressure was 7.308675 MPa. The CO2 in the simulations was treated as a real gas and the Peng-Robinson equation of state [14] was used to calculate the pressure at different temperatures and densities. The properties of CO2, such as specific heat, dynamic viscosity, and thermal conductivity were supplied to the Starccm+ solver via tables taken from NIST property database miniREFPROP [15].

Conclusions
Using the explicit-unsteady methodology, this study successfully revealed the piston effect, which and can match experimental results inside microchannels. The time scales at which the heat transport occurred due to PE are orders of magnitude faster than the diffusion time scales. The Peng-Robinson equation of state was successfully used to approximate the density of Carbon Dioxide near its critical point. The adiabatic thermalization of the bulk fluid due to the PE is significant under terrestrial conditions and even under some forced convective conditions. The time periods required to reach the maximum temperature indicate that the PE is significant in forced convection.

Data availability
The data generated and analyzed can be available from the corresponding author upon reasonable request.