DNS of turbulent channel flow subject to oscillatory heat flux

In this paper we study the heat transfer in a turbulent channel flow, which is periodically heated through its walls. We consider the flow of air and water vapor using direct numerical simulation. We consider the fluid as a compressible Newtonian gas. We focus on the heat transfer properties of the system, e.g., the temperature difference between the walls and the Nusselt number. We consider the dependence of these quantities on the frequency of the applied heat flux. We observe that the mean temperature difference is quite insensitive to the frequency and that the amplitude of its oscillations is such that its value multiplied by the square root of frequency is approximately constant. Next we add droplets to the channel, which can undergo phase transitions. The heat transfer properties of the channel in the case with droplets are found to increase by more than a factor of two, compared to the situation without droplets.


Introduction
Turbulent flows laden with a large number of small heavy droplets play an important role in many applications such as pollutant dispersion in the atmosphere and heat transfer in power stations.It is important to investigate the heat properties of the flow, which are influenced by the droplet behaviour.In the past fifteen years several studies were conducted in order to investigate droplet-laden turbulent flows.In 1997 Mashayek [1] conducted an Euler-Lagrange simulation study of droplet-laden homogeneous turbulence with twoway coupling in momentum, mass and energy.A study of the mixing layer with embedded evaporating droplets was reported by Miller and Bellan [2].In 2010 Masi et al. [3] investigated the interaction of a non-isothermal dropletladen turbulent planar jet with a cloud of inertial evaporating droplets.In this study we focus on wall-bounded turbulence by considering turbulent channel flow with a dispersed droplet phase undergoing phase transition.As a point of reference we consider the flow of water droplets in air, in which the presence of water vapor is accounted for.Throughout the paper under the carrier phase or carrier gas a mixture of air and water vapor is meant while water droplets will be referred to as the dispersed phase.The top wall of the channel is heated periodically and the bottom wall is cooled in such a way that the total energy of the system is conserved.The key difference between the previously mentioned studies and the current is that we consider conditions in which not only evaporation of droplets but also their growth by condensation are important.Studies [4] and [5] show that a time-periodic external agitation changes the properties of the turbulent flow.In this study we do not change the flow itself but impose a time-periodic heat flux through the walls of the channel which has the form: ܳ ሶ ൌ ܳ ሺͳ ߳ ߱‫ݐ‬ሻ where ߱ will be referred to as the frequency of the heat flux.In [6] the heat transfer properties of the system were investigated in case of a constant heat flux through the walls of the droplet-laden turbulent channel flow.In this study we focus on the resulting response of the effective heat transfer in turbulent channel flow without water droplets in case of an oscillating heat flux.We investigate the dependence of the temperature difference between the two walls on the frequency of the heat flux through the walls.We vary this frequency and determine the mean temperature difference, the amplitude of its oscillations and the phase lag which shows how much the resulting signal of the temperature difference lags behind the applied heat flux.In addition, we perform a simulation with droplets and show that the Nusselt number which measures the effectiveness of the heat transfer between the walls is increased considerably by their presence.The carrier gas is treated using the compressible formulation since the applied heat flux causes a temperature gradient in the wall-normal direction of the channel and, consequently, a gradient in mass density of the carrier gas, which is achieved by a mean motion of gas from the hot to the cold wall .This mean gas flow in the wall-normal direction is quantified by the wall-normal momentum component averaged in the homogeneous directions and in time which is in the incompressible formulation is equal to zero.The organization of the paper is as follows.In section 2 we present the governing equations for the two phases.In section 3 the numerical algorithm for the treatment of the low Mach number turbulent droplet-laden channel flow with phase transitions will be presented.In section 4 we describe the initial conditions and discuss the results.Finally, the concluding remarks are collected in section 5.

Governing equations for the gasdroplets system
In this section we first discuss the flow set-up.Subsequently, the set of partial differential equations for the carrier phase, i.e., the gas consisting of air and water vapor, the system of ordinary differential equations for the dispersed phase and the source terms describing the coupling between the two phases will be presented in separate subsections.

Description of the flow domain
We consider a water-air system in a channel, bounded by two parallel horizontal plates.We focus on a two-phase system, consisting of a carrier phase and a dispersed phase, consisting of liquid water droplets.The carrier phase is represented using the Eulerian approach while the dispersed phase is treated in the Lagrangian manner.In Figure 1 a sketch of the flow domain is presented.The domain has a size of Ͷߨ‫ܪ‬ in the streamwise direction, which is denoted by ‫,ݔ‬ and ‫ܪߨʹ‬ in the spanwise direction, ‫,ݖ‬ where ‫ܪ‬ is half the channel height.In addition, ‫ݕ‬ is the coordinate in the wall-normal direction.The total volume of the domain is defined by ܸ.The top wall of the channel is denoted by ‫ݕ‬ ൌ ‫ܪ‬ and the bottom wall by ‫ݕ‬ ൌ െ‫.ܪ‬A heat flux of the following form is applied through the walls: ܳ ሶ ൌ ܳ ሺͳ ߳ ߱‫ݐ‬ሻ where ߱ is the frequency of the heat flux, ܳ ǡ ߳ are the mean heat flux and the amplitude of its oscillations.The flux at both walls is equal in size but opposite in sign in order to conserve the total energy of the system.Studies done by Kim et al. [7] motivate the use of periodic boundary conditions in ‫ݔ‬ and ‫ݖ‬ directions.In addition, no-slip conditions are enforced for the gas at the walls.

The carrier gas
We solve the following system of equations for (ߩǡ ‫ݑ‬ ǡ ܶǡ ܻ ௩ ሻ: Figure 1 The droplet-laden channel geometry where ߩǡ ‫ݑ‬ ǡ ܶǡ ܻ ௩ stand for mass density, velocity components, temperature and vapor mass fraction, respectively.The system (1)-( 4) was obtained following the asymptotic analysis in powers of the Mach number by Mሷ ller [8].An important result of this analysis is the decomposition of the pressure into two parts: ‫‬ ௧௧ ൌ ‫‬ ሺ‫ݐ‬ሻ ଶ ‫‬ ଶ ሺǡ ‫ݐ‬ሻ, where ‫‬ depends only on time and ‫‬ ଶ depends on time and spatial coordinates.All variables in (1)-( 4) except the pressure are the leading terms in the expressions for the power series of the Mach number.In (2) ߤ denotes the gas dynamic viscosity and ܵ is the compressible form of the rate-of-strain tensor.Moreover, in (4) ௩ stands for the diffusive mass flux of water vapor.The right-hand side of the temperature equation contains ܵ ் which is a complicated expression and reflects several mechanisms of temperature change: heat transfer because of conduction and because of diffusion and viscous heating.In order to close the system (1)-( 4) we use the ideal gas law for the air-water vapor mixture : where ܴ ǡ ܴ ௪௧ stand for the specific gas constants of air and water vapor, respectively.The terms ܳ ǡ ܳ ௩ ǡ ܳ ௧ denote the two-way coupling terms and will be discussed in the next subsection.
The system (1)-( 4) is made non-dimensional using a set of reference scales of the system.The reference temperature, ܶ , is the initial mean temperature, while the reference mass density, ߩ is the initial mean mass density of the fluid.The reference length ‫ܮ‬ is chosen as half the channel height ‫;ܪ‬ specifically, we consider a water vapor-air system flowing between a channel with ‫ܮ‬ ൌ ʹ .The velocity scale ‫ݑ‬ is taken as the initial bulk velocity ‫ݑ‬ ୠ of the gas.As a result of making the governing equations non-dimensional, the final system of equations contains the Prandtl number Pr, the Mach number Ma and the Schmidt number Sc.The initial conditions which will be given in Section 3 define the non-dimensional parameters of the flow.
The system considered in this study is turbulent channel flow at a Reynolds number ఛ ൌ ͳͷͲ which has often been studied in literature [9].From this value and the initial conditions we find the reference velocity and the reference speed of sound and, consequently, the Mach number which is on the order of 0.005.The considered flow belongs to low Mach number compressible flows.

The dispersed phase
The dispersed phase of the system consists of water droplets.In the current study the droplet volume fraction is chosen to fluctuate around ͳͲ ିସ and therefore we consider two-way coupling according to the classification proposed by Elghobashi [10].The point-particle approach is justified by the maximum ratio of the droplet diameter and the Kolmogorov scale which is equal to approximately 0.3 [9].The detailed model for the droplets can be found in [6].
For each droplet we solve a system of ordinary differential equations following the model in [2].The location of the droplet ‫ܠ‬ is governed by the kinematic condition where ‫ܞ‬ stands for the droplet velocity: In the water-air system droplets have a much higher mass density than the carrier phase and, therefore, the Stokes drag force acting on them is dominant.We solve the following equation of motion with the standard Schiller-Naumann drag correlation [11]: where ‫ܠ‪ሺ‬ܝ‬ ǡ ‫ݐ‬ሻ is the velocity of the carrier gas at the droplet position.Moreover, ߬ ௗǡ defines the droplet relaxation time and ௗǡ is the droplet Reynolds number: where ݀ stands for the droplet diameter.We also consider the change in time of the droplet energy ‫ܧ‬ ൌ ܿ ݉ ܶ where ܿ is the specific heat of liquid water and ݉ ǡ ܶ are droplet mass and temperature, respectively.The droplet energy changes because of evaporation and condensation and heat exchange at the droplet surface: where ݄ ௩ ൌ ߣ ܿ ௩ ܶ denotes the specific enthalpy of water vapor with ߣ and ܿ ௩ the latent heat at 0 K and the specific heat at constant pressure of water vapor, respectively.The correlation for forced convection around a sphere is used for the heat transfer coefficient ݄ [12].Finally, ‫ܣ‬ denotes the surface area of the droplet and ܶሺ ǡ ‫ݐ‬ሻ the carrier phase temperature at the droplet position.
In order to close the system of equations we require an expression for ௗ ௗ௧ , the rate of evaporation and condensation of a droplet.We use the following equation [12]: where ܻ ௩ǡఋǡ and ܻ ௩ǡ௦ǡ are the vapor mass fractions at a distance ߜ from the surface of the droplet and on the surface of the droplet, respectively.Moreover, ߜ describes the thickness of the vapor film around the droplet.For the Sherwood number Sh we use the correlation for a sphere [12].In addition, ܻ ௩ǡ௦ is found using Antoine's relation [13] and, subsequently, the ideal gas law to find the water vapor fraction.In order to find ܻ ௩ǡఋ we recall the point-particle approach and treat it as the vapor mass fraction of the carrier phase at the position of the droplet, i.e., ܻ ௩ǡఋ ൌ ܻ ௩ ሺ‫ܠ‬ ǡ ‫ݐ‬ሻǤ The two-way coupling terms present in ( 1)-( 4) ܳ ǡ ܳ ௧ǡ ܳ ௩ are found from the conservation of the total mass, energy and momentum of the system.For example, the two-way coupling term in ( 1) and ( 4) is given by: where the sum is taken over all the droplets in the domain.The delta-functions express that the coupling terms act only at the positions of the droplets, consistent with the underlying point-particle assumption.This completes the description of the governing equations for the two phases.In the next section we present the numerical algorithm for the mathematical model presented in this section.

The low Mach number time-stepping method
The decomposition of the pressure into two parts and the presence of its space-dependent part in equation ( 2) only motivate us to use the pressure-correction method for low Mach number compressible flow.We extend the algorithm proposed by Bell et al. [14] to two-phase turbulent channel flow with phase transitions.
The system of governing equations for the carrier phase (1)-( 4) with the vector of unknowns ‫ܟ‬ ൌ ሾߩǡ ‫ݑ‬ ǡ ܶǡ ܻ ௩ ሿ is written as: where ‫ܮ‬ and ܰ are linear and nonlinear operators, respectively.The linear operator consists only of the term െ ଵ ఘ ߲ ‫‬ ଶ , all the other terms in equations ( 1)-( 4) belong to the nonlinear part.The right-hand sides of the ordinary differential equations for the dispersed phase ( 6), ( 7), ( 9) and ( 10) are also considered as the nonlinear part.For the time integration a three-stage low-storage hybrid implicit-explicit Runge-Kutta method is used [15].In this method only the values of the unknowns from the previous stages are used in the nonlinear operators, while the linear operator also requires the values at the current stage.We do not solve all quantities simultaneously.We do it in the following order.First, we solve the system of equations for the dispersed phase, ( 6), ( 7), ( 9) and ( 10) explicitly using droplets and gas values at the previous stage.The right-hand sides of the droplet equations contain gas properties at the droplet location.In order to determine these, tri-linear interpolation is applied [9].After finding the droplets properties at the new stage, we also obtain the two-way coupling terms at this stage.Next, we solve the system of equations for the carrier phase: first, we find the temperature and vapor mass fraction at the new stage, next, the mass density and finally, the velocity.Such an order is chosen because the values of the temperature, pressure, mass density and vapor mass fraction at the new stage are required for the velocity which is solved partially implicitly.
The temperature and the vapor mass fraction equations, (3) and ( 4), are solved explicitly.These quantities are transported with the speed of the flow and, consequently, the restriction on the time step from the stability condition does not lead to very small time steps.A special procedure is applied to find the mass density of the carrier phase.The mass density can be calculated from equation of state (5) if apart from temperature and vapor mass fraction also ‫‬ is known.In order to find ‫‬ we integrate (1) over the total computational domain ܸ and obtain the total mass at the new stage, ݉.We express the mass density from equation of state (5) and integrate this expression over the computational domain in order to find ݉ using the dependence of ‫‬ on time only.As a result we obtain the following expression in which ݉ǡ ܶ and ܻ ௩ are known: 13) we find ‫‬ and the mass density is obtained from (5).For the carrier gas velocity we apply an extension of the pressure-correction procedure usually applied to incompressible flows, [16].It consists of three steps.First, a provisional velocity from the Navier-Stokes equation, based on velocity, density, temperature and pressure from the previous time is calculated.This velocity does not satisfy the constraint on the divergence of velocity which will be derived momentarily and which in case of incompressible flow is ߘ ‫ڄ‬ ‫ܝ‬ ൌ Ͳ.In order to correct the provisional velocity a Poisson equation for pressure is derived from this constraint and it is solved during the second step.During the third step we correct the provisional velocity with this pressure.The expression for the divergence of the velocity is obtained following a few steps.To find the constraint on velocity for compressible flow we start with the expression for the material derivative of ‫‬ following Bell et al. [14], using equation of state (5): Next we substitute the three material derivatives in (14) from the governing equations ( 1)-( 4).All partial derivatives in ( 14) follow from equation of state (5).
After these substitutions we obtain an expression for ߘ ‫ڄ‬ ‫ܝ‬ in which the only unknown term is బ ௧ . In order to find this missing term we integrate the expression for ߘ ‫ڄ‬ ‫ܝ‬ over the computational domain applying the boundary conditions for the velocity.Using the independence of ‫‬ of the spatial coordinates we obtain the desired expression for బ ௧ and, consequently, we find the expression for the divergence of the velocity.Knowing this expression for ߘ ‫ڄ‬ ‫,ܝ‬ we can apply the pressure-correction procedure.First, we find a provisional velocity ‫ܝ‬ ‫כ‬ in an explicit way from (2).This velocity does not satisfy the requirement for the divergence of the velocity.During the second step we find ‫‬ ଶ at the new stage from the Poisson equation: where ݅ ͳ defines the new stage and ݅ is equal to 0, 1 or 2, and ߚ is one of the coefficients in the time-integration algorithm [15].The expression for ߘ ‫ڄ‬ ‫ܝ‬ ା is known.It contains ‫‬ , the temperature, the vapor fraction and the mass density at the new stage.All these values were found before during the current stage and it is possible to calculate the right-hand side of ( 15) and, consequently, to find ‫‬ ଶ ାଵ .During the third step of the pressure-correction procedure we apply the correction: The final velocity ‫ܝ‬ ାଵ satisfies the expression for the divergence of the velocity referred to previously.This finishes the description of one stage of one time-step.

Channel heat transfer characteristics 4.1 Initial conditions
The simulations are started from a turbulent velocity field obtained from a simulation without droplets and with adiabatic boundary conditions.This simulation was done until the statistically steady state was satisfactorily developed.We use the following initial conditions: atmospheric pressure, room temperature and relative humidity equal to 100%.These conditions correspond to the case examined in [4].The non-dimensional parameters of the simulations along with the thermodynamic parameters of the flow for the reference temperature of 293.15K can be found in [4].The initial condition for the simulations is a snapshot of the solution of the system (1)-( 4) in the absence of droplets and with adiabatic boundary conditions.Consequently, we start from a constant mass density and a constant temperature; the initial vapor mass fraction is found from the initial condition for the relative humidity.We randomly distribute 2,000,000 droplets over the volume of the channel.The initial droplet diameter is given by ௗ ு ൌ ͵ǤͲͻ ൈ ͳͲ ିଷ .These initial conditions permit us to consider the two-way coupling regime but not four-way coupling yet and to treat droplets using the point-particle approach, [9], [10].In this study we consider a periodic heat flux through the walls of the form: ܳ ሶ ൌ ܳ ሺͳ ߳ ߱‫ݐ‬ሻ where ܳ is equal to ͵ʹ ୫ మ .A constant heat flux of this value was applied in [6] to droplet-laden turbulent channel flow and caused a temperature difference approximately equal to 3 K between the two walls in the steady state.We choose ߳ to be equal to ½ to have significant heat transfer modulation caused by the periodic heat flux through the walls.In [6] we found that the transient stage of the system lasts for approximately 0.1 s.This transient arises from the development of a mean temperature gradient.This motivates us to consider the frequency ߱ ൌ ʹͲߨ ିଵ which corresponds to a period equal to the typical time scale of the transient.Along with ߱ ൌ ʹͲߨ ିଵ we also perform simulations of the turbulent channel flow without droplets at ߱ ൌ Ͳǡ ʹߨǡ ʹͲͲߨ ିଵ to investigate the channel heat transfer properties.Next we add droplets to the channel and perform a simulation at ߱ ൌ ʹͲߨ ିଵ .

Motivation for compressible formulation
In this study compressibility refers to the changes in the mass density due to the development of the temperature gradient in the wall-normal direction caused by the applied heat flux through the walls.Consequently, we get a net transfer of gas from the upper hot wall to the bottom cold wall.In order to quantify this, we considered the gas wall-normal momentum component averaged over the homogeneous directions and time during the initial transient stage and in the statistically steady state for case of ߱ ൌ Ͳ, Figure 2.During the initial stage it is everywhere negative which is consistent with the gas movement due to the temperature gradient between the walls.In the statistically steady state this momentum component is close to zero and this result is confirmed by averaging continuity equation (1) over the periodic directions and time.In the incompressible formulation this averaged quantity is always equal to zero [17].The net wall-normal gas movement during the initial stage motivates us to use the compressible formulation.

Turbulent channel flow
In this subsection we analyse the heat transfer characteristics in terms of the behaviour of the temperature difference between the two walls ߂ܶ and the Nusselt number.First, we consider ߂ܶ for ߱ ൌ Ͳ, Figure 3.The temperature difference reaches the steady state at a time approximately equal to 2 s.Starting from this time, we average ߂ܶ in time and obtain a mean value denoted as ߂ܶ equal to approximately 8.47 K, Table 1.

Figure 3
The temperature difference between the walls for ߱ ൌ Ͳ For the case of an oscillating heat flux the temperature difference between the walls ߂ܶ is written as: ߂ܶ ൌ ߂ܶ ሺͳ ߙ ሺ߱‫ݐ‬ ߶ሻሻ ߂ܶ ᇱ ǡ ሺͳሻ where ߂ܶ ǡ ߙǡ ߶ stand for the mean temperature difference, the amplitude of its oscillations and the phase lag.Moreover, ߂ܶ ᇱ denotes the fluctuating contribution due to turbulence.In Figure 4 we show ߂ܶ as a function of time for ߱ ൌ ʹͲߨ ିଵ .

Figure 4
The temperature difference between the walls for ߱ ൌ ʹͲߨ ିଵ .
We perform phase averaging of the resulting temperature difference for all non-zero values of frequency considered in this paper.This means that we collect full periods of a resulting signal starting from the moment of time in which the signal of ߂ܶ is in the statistically steady state.The period is equal to ଶగ ఠ .We divide each period in 700 points.Next we perform averaging over all periods of the signal, starting from the statistically steady state.The phase-averaged signal for ߱ ൌ ʹͲߨ ିଵ is shown in Figure 5. From the phase-averaged signal we find ߂ܶ ǡ ߙǡ ߶ and the turbulent contribution.In particular, ߂ܶ is obtained as the mean temperature of the phaseaveraged signal while ߙ is the amplitude.In addition, the phase lag shows how much the phase-averaged signal lags behind the applied heat flux.We calculate the statistical error caused by the choice of the time interval used for averaging.We perform averaging over a certain time interval: for ߱ ൌ ʹͲߨ ିଵ it is [2 s; 5 s], Figure 4.In order to find the error we divide this interval into N equal parts and we calculate the value of all characteristics of the signal by averaging over these intervals separately.At each interval i we find ܽ where a denotes ߂ܶ ǡ ߙǡ ߶ or RMS of ߂ܶ ᇱ .Then we calculate the averaged value of a which is denoted by ܽ ത.Next the error is computed as: The mean values along with the statistical errors are shown in Table 1.Table 1 The characteristics of the resulting signal for the turbulent channel flow.
We conclude that the mean temperature does not depend on the frequency of the heat flux through the walls.At the same time we find that with increasing frequency the phase lag also increases, Table 1.The RMS of the turbulent contribution does depend on the frequency of the heat flux, but without clear trend .The calculated values of ߙ permit to find its dependence on frequency which is ߙξ߱ ൌ ‫ݐݏ݊ܿ‬ in good approximation, figure 6.

Droplet-laden turbulent channel flow
Also a simulation with droplets at ߱ ൌ ʹͲߨ ିଵ was performed.We also performed phase averaging of the temperature difference between the walls for this case.The resulting ߂ܶ in this case is equal to 3.1 K.The values of ߙ and ߶ are twice smaller than without droplets at the same frequency.In [6] a simulation with droplets at ߱ ൌ Ͳ was performed.From its results we find that ߂ܶ in this case is also approximately equal to 3.1 K. Consequently, the mean temperature difference also does not depend on the frequency of the heat flux also in the presence of water droplets.
To express the efficiency of heat transfer between the walls we consider the Nusselt number, which is defined in the following way: where ݇ ௪ stands for the thermal conductivity of the gas at the wall.The Nusselt number is calculated according to (18) in cases without and with droplets at ߱ ൌ ʹͲߨ ିଵ .The difference in ߂ܶ leads to a Nusselt number which is larger by a factor more than two when droplets are present.

Conclusions
Droplet-laden turbulent channel flow with phase transitions was simulated with a special time integration method developed for two-phase turbulent low Mach number flows undergoing phase transitions.A periodic heat flux was applied through the walls.We varied the frequency of this flux between 0 and ʹͲͲߨ ିଵ Ǥ First, we performed simulations of the turbulent channel flow in the absence of droplets.We found that the mean temperature difference between the walls of the channel is not dependent on the frequency of the applied heat flux.At the same time, the amplitude of oscillations of the temperature difference decreases with increasing frequency.The dependence of this amplitude on the frequency was found to be inversely proportional to the square root of the frequency in good approximation.The phase lag was found to increase if we increase the frequency of oscillations.In the future we plan to investigate this dependence better performing more simulations with intermediate values of frequency.
The presence of droplets which can undergo phase transitions was seen to increase the heat transfer within the channel by more than a factor of 2. We observe that the mean temperature difference is independent of the frequency of the applied heat flux.In addition, the amplitude of the temperature difference oscillations and its phase lag were found to decrease by approximately a factor of 2 in comparison with the case without droplets.In the future, we plan to investigate resonant turbulence further by modulating the turbulent flow and using a periodic time dependent pressure-gradient.

Figure 2
Figure 2 The wall-normal momentum component averaged over the homogeneous directions and over time as a function of the wall-normal coordinate.Time averaging is performed for dashed: [0s; 0.1s], solid: [2s; 4s].

Figure 6
Figure 6 Amplitude ߙ of phase-averaged temperature difference.The solid line shows the dependence ߙξ߱ ൌ ‫;ݐݏ݊ܿ‬ the symbols represent the simulation results.