Heat and mass transfer in turbulent multiphase channel flow

(1)

HEAT AND MASS TRANSFER

IN TURBULENT MULTIPHASE

CHANNEL FLOW

HEA

T AND MA

SS TR

ANS

FER IN TURB

ULENT MUL

TIP

H

A

SE CH

ANNEL FLO

W

AN

A

ST

A

SIA B

UKHV

OS

TO

VA

ANASTASIA BUKHVOSTOVA

ISBN: 978-90-365-3881-7

(2)

H

EAT AND MASS TRANSFER

IN TURBULENT MULTIPHASE

CHANNEL FLOW

(3)

De promotiecommissie:

Voorzitter en secretaris:

Prof. dr. ir. A.J. Mouthaan Universiteit Twente

Promotoren:

Prof. dr. ir. B.J. Geurts Universiteit Twente Prof. dr. J.G.M. Kuerten Universiteit Twente

Leden:

Prof. dr. ir. J. Derksen Technische Universiteit Delft Prof. dr. J. Harting Universiteit Twente

Prof. dr. ir. G.J.F. van Heijst Technische Universiteit Eindhoven Prof. dr. D. Lohse Universiteit Twente

Prof. dr. B. M¨uller Norwegian University of Science and Technology

The research presented in this thesis was done in the group Multiscale Modeling and Simu-lation (Dept. of Applied Mathematics), Faculty EEMCS, University of Twente, The Nether-lands as part of the research project of the Foundation for Fundamental Research on Matter (FOM) 08DROP02-2 which is part of the Netherlands Organization for Scientific Research (NWO).

Computing resources were granted by the National Computing Facilities Foundation (NCF), with financial support from the Dutch Organization for Scientific Research (NWO).

Heat and mass transfer in turbulent multiphase channel flow. Ph.D. Thesis, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands. c Anastasia Bukhvostova, En-schede, 2015.

ISBN : 978-90-365-3881-7 DOI : 10.3990/1.9789036538817 Printed by Gildeprint, Enschede

(4)

H

EAT AND MASS TRANSFER

IN TURBULENT MULTIPHASE

CHANNEL FLOW

PROEFSCHRIFT

ter verkrijging van

de graad van doctor aan de Universiteit Twente, op gezag van de rector magnificus,

prof. dr. H. Brinksma,

volgens besluit van het College voor Promoties in het openbaar te verdedigen

op vrijdag 19 juni 2015 om 14:45 uur

door

Anastasia Bukhvostova geboren op 25 april 1989

(5)

Dit proefschrift is goedgekeurd door de beide promotoren:

(6)

Chapter 1 Introduction

1.1 Droplets in turbulent flow

Turbulent multiphase flows with a large number of dispersed droplets in a gas are present in many processes both in nature and in industry. The possible enhancement of heat transfer by the presence of particles is an important aspect in the design of more efficient devices. On larger scales the investigating of droplet dynamics in turbulent atmospheric clouds can help to describe and understand their influence on our climate.

Turbulent flow itself is considered a major problem of classical physics. It can be modeled in various ways and moreover, different numerical and experimental methods were developed during the last century in order to increase our understanding of turbulence. The presence of a large number of particles and droplets in a turbulent flow makes the problem even more chal-lenging from both modeling and computational points of view. Several studies were dedicated to the study of fundamental issues of turbulent droplet dynamics by means of the Lagrangian approach (Shraiman and Siggia 2000, Procaccia 2001). In particular, on the smallest turbulent scales droplets which are heavier than the fluid are propelled away from vortices and gather in regions of strain (Cencini et al. 2006). Moreover, if their volume and mass fraction are suf-ficiently large, they influence the turbulent flow itself. Particles with different Stokes number, which characterizes how fast the particle adopts to the flow, were investigated in Calzavarini et al. (2008). In Figure 1.1, left we observe that heavy particles are subject to preferential clustering.

The current PhD study is part of a research project called DROP, funded by FOM, the aim of which is a better quantitative and fundamental understanding of the dynamics of droplets in turbulent flows. The project is carried out by research groups in the Netherlands from Delft, Eindhoven and Twente, and includes four experimental and four numerical studies. In par-ticular, the study presented in this PhD thesis concentrates on the mass and heat transfer in droplet-laden turbulent channel flow in which droplets undergo phase transition. The study conducted by Kuerten et al. (2011) showed an enhancement of the heat transfer in turbulent channel flow by more than a factor of two by introducing heavy inertial particles. We take the next step in the physical modeling by introducing phase change between dispersed droplets and carrier gas and investigate its effect on the global heat transfer properties of channel flow and on the droplet behavior. The modeling challenge is to take into account not only the exchange of momentum between the two phases but also of mass and energy, taking the dispersed nature of the liquid droplet phase in the continuous gas flow into account using an

(11)

Fig. 1.1: Snapshot of the distribution of heavy particles in a turbulent flow field (right) com-pared to the case of tracers (left). From Calzavarini et al. (2008) .

Euler-Lagrange formulation (Mashayek 1998, Miller and Bellan 1998).

Before presenting the contents of the thesis, we will briefly discuss the choices which we made in physical modeling and numerical treatment of the problem. Next, the key points of each chapter will be presented along with the general outline of this thesis.

1.2 Models and methods in the context of modern computational

capabilities

A turbulent flow can be modeled in various ways. In this thesis we make use of direct numer-ical simulation in which all scales of the flow, except the detailed flow around each droplet, are resolved. Droplets are solved using the point-particle approach which is common practice when the size of the droplets is smaller than the Kolmogorov scale (Marchioli et al. 2006). This permits to maintain the computational costs within acceptable levels compared to the turbulent flow case but without droplets.

We adopt the Euler-Lagrangian formulation which imposes conservation laws for mass, mo-mentum, energy and species of the carrier phase, while using empirical correlations for the interaction between the two phases. Two-way coupling terms are derived reflecting the inter-action of droplets with the carrier gas conserving mass, momentum and energy. The model equations used for the two phases resemble the approach used by Miller and Bellan (1998), Masi et al. (2010).

One of the issues while modeling the carrier gas is whether to use a compressible or an incompressible formulation. The Mach number, which is a non-dimensional parameter char-acterizing the compressibility of the flow, is very low for the conditions which we simulate. It is approximately equal to 0.005. It is common practice to treat the flow as incompressible for Mach numbers smaller than 0.3 (Anderson 2007). In the present study the carrier gas consists of water vapor and dry air. If the carrier gas is assumed to be strictly incompressible, then all

(12)

instantaneous changes in the mass density of vapor and air should cancel each other precisely throughout the domain. A full simulation model was developed for such an incompressible carrier phase by Russo et al. (2014a). Here, we confront the incompressible model with a fully compressible description, which allows to quantify the consequences of non-constant mass density of the carrier phase and indicates in which respects and under what conditions the full compressible formulation becomes essential. While the kinematics of the flow are not very sensitive to the Mach number in the low-Mach regime, more subtle influences, deriving from thermodynamics associated with the phase changes within the flow, display a stronger sensitivity to compressibility.

Another challenge which we face in this study is the choice of a suitable time integration method. We perform the first part of the work using a density-based fully explicit time inte-gration method. In order to avoid excessive computational requirements caused by the sta-bility conditions of this time integration method, this method is applied to low values of the Mach number that are, however, much higher than the actual value of 0.005. In order to enable efficient simulation at the real Mach number we develop a different time integration method which belongs to the class of pressure-based methods for compressible flows. We extend ear-lier work by Bell et al. (2004) to two-phase flows with phase transitions and validate it using the results of the fully explicit solver.

The study of Russo et al. (2014a) and the first chapters of this thesis consider flows without gravity. If gravity in the wall-normal direction is taken into account, droplets migrate contin-uously towards the bottom wall and will typically create a film of water there. This motivates us to investigate the extended setting in which initially a thin film of water covers the bottom wall while a large number of droplets is dispersed in the domain. In this work the film is con-sidered to have a constant height and its surface does not move. Moreover, we keep the total mass of water in the system constant by draining water from the bottom wall and inserting, on average, the same amount in the form of new droplets released at the top wall. We investigate in detail the behavior of the droplets along with the heat transfer properties of the channel in the presence of gravity.

1.3 Structure of the thesis

This thesis comprises three main developments in the field of turbulent multiphase flows with phase transitions. These will be discussed briefly next.

1.3.1 Comparison of compressible and incompressible formulations

The second chapter of the thesis is dedicated to the description of the problem and the mod-eling equations for the two phases using either a compressible or an incompressible formu-lation for the carrier phase. We investigate the sensitivity of the system to a low level of compressibility by confronting the results from these two formulations. We choose ‘mild’ initial conditions of room temperature, atmospheric pressure and 100% relative humidity in order to guarantee ‘modest’ coupling between the phases expressed by the mass transfer. Consequently, only small differences between the results of the two formulations should be

(13)

expected, serving as a point of reference. We concentrate on heat and mass transfer quantities of the channel. The results give evidence of the need for the compressible formulation even for the lowest value of the Mach number of 0.05 investigated. Changes in, e.g., the Nusselt number of up to 15% were observed when the mean system temperature is elevated to about 50 degrees C, while overall the agreement between the two formulations is quite close. In the third chapter we present the results of a study into compressibility effect due to phase transitions and choose significantly lower values of the initial relative humidity. This has the effect of emphasizing evaporation during initial stages which has a marked influence on droplet sizes and associated heat transfer. We concentrate on the Nusselt number (Kuerten et al. 2011) and analyze all contributions to it. The main finding of this study is a 1.5 larger difference in the Nusselt number between the two formulations in case of initial relative hu-midity of 50% than in the results for the mild case of a saturated initial state.

1.3.2 Low Mach number time integration algorithm

The fourth chapter of the thesis is dedicated to the new time integration algorithm for low Mach droplet-laden turbulent channel flow with phase transitions. There are two main nu-merical techniques for the treatment of compressible low Mach number flows: density-based methods and pressure-based methods. Two types of time-marching procedures are applied in density-based methods: explicit and implicit algorithms. The severe restriction on the time step for low values of the Mach number in explicit solvers is dictated by a stability restriction, called Courant-Friedrichs-Lewy (CFL) condition (Morton and Mayers 2005). We faced this problem in the study described in chapters two and three where the fully explicit time inte-gration method applied in the compressible formulation required long computational times. In implicit density-based methods the system of governing equations for compressible tur-bulent flow is ill-conditioned, making iterative solution methods excessively time consuming (Roller and Munz 2000).

To alleviate these problems for unsteady low-Mach flows we focus on a pressure based method which is very similar to pressure correction methods used for incompressible flow. We apply this method to a system of equations obtained by asymptotic analysis at low Mach numbers. The main virtue of this approach is the absence of large eigenvalues which allows to treat most of the terms explicitly. We closely follow the work of Bell et al. (2004) and extend it to multiphase problems. This novel method is validated by comparing results obtained with those of the fully explicit solver from the second and third chapters of this thesis.

1.3.3 Introduction of water film on bottom wall

In the fifth chapter of the thesis we describe a more realistic situation when gravity acts in the wall-normal direction. This force results in the movement of droplets towards the bottom wall where they merge and form a film of water. We model the film as a sufficiently thin layer neglecting internal motion and keeping it to a constant height. In order to maintain the height we add or drain water near the bottom wall. Simultaneously, we keep the total mass of water in the system constant on average by distributing new water droplets into the channel from the

(14)

top wall. We express the film temperature approximating it by a second order polynomial in the wall-normal coordinate (Russo et al. 2014b) and consider the film energy per unit area as dependent variable. The film is coupled to the carrier phase by evaporation and condensation of water vapor at its surface. The coupling between the film and the dispersed phase is through droplets merging with the film as they migrate downward due to gravity. In this chapter we analyze the behavior of droplets and the influence on the heat transfer due to the inclusion of the thin film.

(15)

(16)

References

Shraiman, B.I. and Siggia, E.D. 2000. Scalar turbulence. Nature 405, 639-646. Procaccia, I., 2001. Go with the flow. Nature 409, 993-995.

Cencini., M, Bec, J., Biferale, L., Boffetta, A., Celani, A., Lanotte, A.S., Musacchio, S. and Toschi, F., 2006. Dynamics and statistics of heavy particles in turbulent flows. J. Turbul. 7, 36.

Calzavarini, E., Kerscher, M., Lohse, D. and Toschi, F. 2008. Dimensionality and morphology of particle and bubble clusters in turbulent flow. J. Fluid Mach. 607, 13-24.

Kuerten, J.G.M., van der Geld, C. W. M. and Geurts, B.J., 2011. Turbulence modification and heat transfer enhancement by inertial particles in turbulent channel flow. Phys. Fluids 23, 123301.

Mashayek, F., 1998. Droplet-turbulence interactions in low-Mach-number homogeneous shear two-phase flows. J. Fluid Mech. 367, 163-203.

Miller, R.S. and Bellan, J., 1998. Direct numerical simulation of a confined three-dimensional gas mixing layer with one evaporating hydrocarbon-droplet-laden stream. J. Fluid Mech. 384, 293-338.

Marchioli, C., Picciotto, M., Soldati, A., 2006. Particle dispersion and wall-dependent turbu-lent flow scales: implications for local equilibrium models. J. Turbul. 7, 60.

Masi, E., Simonin, O., B´edat, B., 2010. The mesoscopic Eulerian approach for evaporating droplets interacting with turbulent flows. Flow Turbul. Combust. 86, 563-583.

Anderson, J.D., 2007. Fundamentals of Aerodynamics. McGraw-Hill series.

Russo, E., Kuerten, J.G.M., van der Geld, C.W.M. and Geurts, B.J., 2014. Water droplet condensation and evaporation in turbulent channel flow. J. Fluid Mech. 749, 666-700 . Bell, J.B., Day, M.S., Rendleman, C.A., Woosley, S.E. and Zingale, M.A., 2004. Adaptive

low Mach number simulations of nuclear flame microphysics, J. Comput. Phys. 195, 677-694.

Morton, K.W. and Mayers, D., 2005. Numerical solution of partial differential equations, second ed., Cambridge University Press.

Roller, S. and Munz, C.D., 2000. A low Mach number scheme based on multi-scale asymp-totics, Computing and Visualization in Science 3, 1/2, 85-91.

Russo, E., Kuerten, J.G.M. and Geurts, B.J., 2014. Delay of biomass pyrolysis by gas-particle interaction. J. Anal. Appl. Pyrol 110, 88-99.

(17)

(18)

Chapter 2 Comparison of DNS of compressible and

incompressible turbulent droplet-laden heated

channel flow with phase transition

Abstract

Direct numerical simulation is used to assess the importance of compressibility in turbulent channel flow of a mixture of air and water vapor with dispersed water droplets. The dispersed phase is allowed to undergo phase transition, which leads to heat and mass transfer between the phases. We compare simulation results obtained with an incompressible formulation with those obtained for compressible flow at various low values of Mach number. We discuss differences in fluid flow, heat- and mass transfer and dispersed droplet properties. Results for flow properties such as mean velocity obtained with the compressible model converge quickly to the incompressible results in case the Mach number is reduced. In contrast, thermal properties such as the heat transfer, characterized by the Nusselt number, display a systematic difference between the two formulations on the order of 15%, even in the low-Mach limit. This shows the necessity of the use of a compressible formulation for accurate prediction of heat transfer, even in case of an initial relative humidity of 100%. Mass transfer properties display a difference between the models on the order of 5%, for example in the prediction of the droplet mean diameter near the walls.

2.1 Introduction

Multiphase flows with a large number of droplets dispersed into a gas play an important role in a variety of technological and environmental applications. Examples include thermal processing in food manufacturing, air pollution control and heat transfer in power stations (Barigou et al. 1998). In this paper we investigate a coupled Euler-Lagrange model to sim-ulate droplet-laden turbulent channel flow in which phase transition plays an important role. We study in particular the effect of compressibility at low Mach numbers on the vapor and temperature fields and on the size distribution of the droplets. While the kinematics of the flow are found to be rather insensitive to compressibility at low Mach numbers, we observe a range of systematic differences in heat and mass transfer characteristics.

Not many studies focus on mass and heat transfer in droplet-laden turbulent flow. The first study was done by Mashayek in 1998. He conducted an Euler-Lagrange simulation study, investigating homogeneous turbulence with two-way coupling in momentum, mass and

(19)

ergy of the system. Later an overview of suitable techniques for the treatment of such type of problems was presented in Mashayek and Pandya (2003), and in Wang and Rutland (2006). Later, a study dedicated to the mixing layer with embedded evaporating droplets was con-ducted by Miller and Bellan (1998). In Masi et al. (2010) a cloud of inertial evaporating droplets, interacting with a non-isothermal droplet-laden turbulent planar jet, is considered. In the present paper we extend the work of Mashayek (1998) to wall-bounded turbulence by investigating turbulent channel flow with a dispersed droplet phase undergoing phase tran-sition. All the studies mentioned above focus on evaporating droplets while in this study droplets may also grow by condensation of the vapor phase. To enable comparison between a strictly incompressible and a more general compressible formulation, we focus the attention on a case of initially 100% relative humidity in the channel. The models adopted here are identical to the models used in Miller and Bellan (1998) and Masi et al. (2010).

We study droplet-laden channel flow in which the top wall is heated uniformly and the bottom wall is cooled such that the total energy of the system is conserved. This leads to a temper-ature gradient in the wall-normal direction. As a point of reference we consider the flow of droplets of water in air, in which the presence of water vapor is accounted for. The mixture of air and water vapor will be referred to as the carrier phase or the carrier gas and liquid droplets as the dispersed phase. The system investigated in this study contains a rather small amount of water vapor, which motivates to approximate several transport characteristics of the carrier phase as those of clean air and investigate trends in flow and in heat and mass transfer properties due to the presence of vapor and liquid water droplets.

Incorporation of evaporation and condensation of the dispersed phase raises the question whether or not to include explicit compressibility of the carrier phase. If the carrier gas is assumed to be strictly incompressible then the inclusion of evaporation and condensation is subject to the condition that all instantaneous changes in the local mass density of air and water vapor cancel each other precisely throughout the domain. A full simulation model can be developed for such an incompressible carrier phase (Russo et al. 2014). Here, we confront the incompressible model with a fully compressible description, which allows to quantify the consequences of non-constant mass density of the carrier phase and indicates in which respects and under what conditions the full compressible formulation becomes essential. We investigate the sensitivity of flow and transport characteristics arising from low-level compressibility effects in multiphase systems. For that purpose we analyze the results from simulations using two different approaches, i.e., a fully incompressible flow model and a compressible formulation. We compare thermal quantities of the system, such as the Nusselt number and mean temperature profiles, along with properties of the dispersed phase, such as accumulation of droplets near the walls and droplet size distribution history, in order to conclude for which quantities and flow conditions the compressible formulation should be used and for which other quantities the approximate incompressible model is adequate. Dif-ferences are deliberately kept small by focussing on an initial relative humidity of 100%; this helps to identify the sensitivity of resulting quantities to compressibility effects at low Mach numbers.

The organization of this paper is as follows. In Section 2.2 the mathematical model is formu-lated for the coupled droplet-carrier gas system. The numerical methods for both formulations are presented in Section 2.3. The initial conditions for the two models are described in Sec-tion 2.4. Finally, the results of the incompressible and compressible formulaSec-tions regarding flow properties, heat and mass transfer and dispersed phase characteristics are presented in Section 2.5 and the concluding remarks are collected in Section 2.6.

(20)

2.2 The governing equations

This section is dedicated to the mathematical model. First, the geometry of the problem will be described along with the applied boundary conditions. Subsequently, the set of partial differential equations for the carrier phase, 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. At the end of each subsection the differences between the equations used in the incompressible model (Russo et al. 2014) and in the current study will be emphasized.

2.2.1 The carrier phase

We consider a water-air system in a channel, bounded by two parallel horizontal plates. In Fig.2.1 a sketch of the domain is presented. The domain has a size of 4πH in the stream-wise direction, which is denoted by x, and 2πH in the spanstream-wise direction, z, where H is half the channel height. In addition, y is the coordinate in the wall-normal direction. The total volume of the domain is defined by V . The top wall of the channel is denoted by y = H and the bottom wall by y = −H. Studies done by Kim et al. (1987) motivate us to use periodic boundary conditions in the streamwise and spanwise directions. In addition, no-slip condi-tions are enforced at the walls. A constant heat flux ˙Qis applied through the walls: this heat flux is positive through the top wall and negative trough the bottom wall. The flux at both walls is equal in size in order to conserve the total energy of the system. As a consequence, the gas temperature is higher near the top wall and lower near the bottom, and a gradient of temperature develops across the channel.

2.2.1.1 The compressible formulation

The carrier phase is treated as a compressible Newtonian fluid. We impose conservation of mass, momentum, total energy density and water vapor. The equations can be written as Mashayek (1998):

∂tρ + ∂j(ρuj) = Qm (2.1)

∂t(ρui) + ∂j(ρuiuj) = −∂ip+ ∂j(µ(T )Si j) + Fi+ Qmom,i (2.2)

∂te+ ∂j((e + p)uj) = ∂j(K(T )∂jT) + ∂j(uiµ (T )Si j) + (2.3)

∂j(ρD∂jYv((cpv− cpa)T + λ )) + Qe

∂t(ρYv) + ∂j(ujρYv) = ∂j(ρD∂jYv) + Qm (2.4)

This system is solved for mass density ρ, velocity components ui, total energy density e and

vapor mass fraction Yv. The source terms Qm, Qmom,i, Qeexpress the interaction between the

phases and will be described in detail in subsection 2.2.3. The term Fiis an external force

(21)

a choice for the calculation of this forcing term in order to maintain the flow in the channel. In the current study it is chosen in such a way that the total momentum in the streamwise direction is constant in time.

In our model ρ and ρYv denote the mass density of the carrier gas and water vapor in the

considered volume, respectively. At the same time the carrier phase is composed of pure air and water vapor:

ρ = ρair+ ρvapor (2.5)

In equations (2.1) and (2.4) the same mass source term Qmappears. This term reflects the

change in mass of the carrier gas and water vapor because of evaporation and condensation. Since only vapor can contribute to the change of mass density of the carrier phase through phase transitions, equations (2.1) and (2.4) have the same coupling term.

The pressure p is defined through the ideal gas law for a two-component mixture in the following way (Miller and Bellan 1998):

p= ρT Rair(1 −Yv) + ρT RwaterYv (2.6)

where T denotes the temperature of the carrier gas. In addition, Rairand Rwaterstand for the

specific gas constants for air and water vapor, respectively. They are equal to the universal gas constant, divided by the corresponding molar mass. In equation (2.2) Si jis the compressible

extension of the rate-of-strain tensor defined as Si j= ∂iuj+ ∂jui−2₃∆ δi jwhere ∆ = ∂kuk

de-notes the divergence of the velocity. Here and elsewhere we adopt the summation convention implying summation over repeated indices.

The dynamic viscosity of the carrier gas µ depends on temperature through Sutherland’s law (Sutherland 1893), i.e.:

µ = µre f _T Tre f 3₂ _T re f+ S T+ S (2.7)

where µre f is the dynamic viscosity of air at the reference temperature Tre f. In addition, S

is the Sutherland temperature and for the case of air S = 110.4 K. In (2.7) the presence of water vapor in the carrier phase is not taken into account since the mixture air - water vapor considered here is assumed to have a sufficiently small vapor mass fraction.

In equation (2.3) K stands for the thermal conductivity of the carrier gas which depends on temperature such that ratio of K and µ is constant (Bird et al. 1960). In addition, D in (2.4) denotes the diffusion coefficient of water vapor in air, which also depends on temperature such that ratio ρD/µ is constant (Bird et al. 1960).

The total energy density of the fluid e is the sum of the kinetic energy density and the internal energy density, i.e., e = eint+ ekin, where (Miller and Bellan 1998)

ekin=

1

2ρ ujuj (2.8)

eint = ρcvaT(1 −Yv) + ρcvvTYv+ ρλYv (2.9)

Here cva, cvv are the values of the specific heats at constant volume of air and water vapor

and λ is the latent heat. Likewise, cpvand cpain equation (2.3) denote the specific heats at

constant pressure of water vapor and air, respectively. The first term on the right-hand side of equation (2.9) is the internal energy density of air while the two other terms denote the

(22)

internal energy density of the vapor contained in the carrier phase. Gas temperature T is cal-culated from the total energy density e using (2.8) and (2.9).

The system of governing equations (2.1)-(2.4) is made non-dimensional using a set of refer-ence scales of the system. The referrefer-ence temperature, Tre f, is the initial mean temperature of

each test case, while the reference mass density, ρre f, is the initial mean mass density of the

carrier gas. The reference length Lre f is chosen as half the channel height H; specifically, we

consider a water-air system flowing between a channel with Lre f= 2 cm. The velocity scale

ure f is taken as the bulk velocity ubof the carrier gas. The Reynolds number based on the

reference scales is defined in the following way:

Reb=

ρrefubLref

µ (Tref)

(2.10)

where µ(Tre f) is the dynamic viscosity of air at Tre f. As a result of making the governing

equations non-dimensional, the final system of equations also contains the Prandtl number Pr, the Mach number Ma and the Schmidt number Sc. Pr is the ratio of kinematic viscosity to thermal diffusivity:

Pr =cpaµ (Tref) K(Tref)

(2.11)

The Mach number expresses the ratio of the reference velocity to the speed of sound at Tre f;

Ma = ure f/c(Tre f) where c(Tre f) is the reference speed of sound. The speed of sound is

calculated using the assumption of a sufficiently small amount of water vapor in the carrier gas, and the ideal gas law:

c(T ) = r

γairRT

Mair

(2.12)

where γairis defined as the ratio between the specific heats at constant pressure and constant

volume of air, equal to 1.4. In addition, Mairis the molar mass of air, Mair= 0.0289645 kg/mol.

The Schmidt number Sc denotes the ratio of kinematic viscosity to mass diffusivity:

Sc = µ (Tre f) ρre fD(Tre f)

(2.13)

We take into account the dependence of viscosity on temperature (2.7) and assume specific heats to be constant. The temperature dependence of the thermal conductivity K and diffusion coefficient D can hence be inferred. The values of the main thermal quantities are shown for different Tre f used in this paper in Table 2.1. The Reynolds number is determined by

spec-ifying the mass flow rate - the selected value of Reb corresponds to values of the friction

Reynolds number Reτabout 150, a well-documented case of turbulent channel flow in

liter-ature (Marchioli et al. (2008)). Throughout, we will fully adhere to the values of Reb, Sc and

Pr. The value of Ma is too low to allow full simulations, including the gathering of converged statistics, with the current numerical method. Instead, we will relax the requirement on Ma and choose low, but considerably higher values of 0.05, 0.1 and 0.2 to investigate the trends. Simulation times of 8-10 months are required on 32 CPU’s of a modern supercomputer at

(23)

Ma = 0.05. Infeasible requirements on available computing budgets would be required for Ma = 0.005. Conversely, as this study will show, much, if not all, of the flow physics and heat and mass transfer characteristics can accurately be approximated using Ma in the se-lected range. Current research into tailored low-Mach time-stepping algorithms that would allow simulation at the physical Ma is ongoing and will be presented elsewhere.

Fig. 2.1: : The computational domain

Tre f(K) 293.15 323.15 K(W/m · K) 2.5519 × 10−2 2.7889 × 10−2 µ (kg/m · s) 1.86766 × 10−52.012 × 10−5 D(m2_/s) _{2.0917 × 10}−5 _{2.4244 × 10}−5 ρre f(kg/m3) 1.194 1.0426 ub(m/s) 1.8246 2.25104 c(m/s) 343.111 360.2403 cpa(J/kgK) 1010 1010 Reb 2333 2333 Pr 0.7478 0.7959 Ma 0.00532 0.00625 Sc 0.6283 0.6505

(24)

2.2.1.2 The incompressible formulation

It is known that for values of the Mach number up to 0.3 the fluid mechanics of turbulent channel flow corresponds well to that of incompressible flow (Anderson 2007). There are three main differences between compressible and incompressible flows: a constant mass den-sity, a zero divergence of the velocity and an infinite speed of sound in incompressible flow. The incompressible formulation for the carrier phase adopted by Russo et al. (2014) implies that evaporation and condensation are subject to the condition that all instantaneous changes in the local mass density of air and water vapor cancel each other precisely throughout the domain. This is required to maintain a constant mass density of the carrier phase. In fact, ρ = ρair(x,t) + ρvapor(x,t) where mass densities are functions of position and time.

Conse-quently, we have the following:

∇ρair= −∇ρvapor, ∂tρair= −∂tρvapor (2.14)

which is an unphysical constraint, inherent to the incompressible formulation.

The incompressible treatment of the carrier phase still permits to incorporate changes in mass density into the model using the equation of state and the Boussinesq approximation, which implies the dependence of mass density on temperature and vapor mass fraction. In our prob-lem the divergence of velocity is not equal to zero because of phase transitions. Apart from this, the heat flux in the considered problem leads to a net transfer of gas from the hot to the cold wall during the initial stages during which the mean temperature gradient in the wall-normal direction develops. This will be discussed in more detail in section 2.5.

Apart from the assumption that the flow is incompressible, the treatment of pressure makes the two formulations different. While an explicit equation of state is used to connect pressure and temperature in the compressible case, the incompressible formulation solves the pressure from a Poisson equation derived from the condition of a divergence-free velocity field. Another difference in the governing equations for the carrier phase is due to the viscous dis-sipation term. While it is present through the term ∂j(uiµ (T )Si j) in the total energy density

equation in the compressible formulation, it is absent in the temperature equation in the in-compressible formulation. This term is expected to have only a small influence on the results as will be discussed in Section 2.5.

In addition, there are also some differences in the physical modeling used in the two ap-proaches. In the compressible formulation we incorporate the dependence of the dynamic viscosity µ on the temperature of the carrier phase through Sutherland’s law. In the incom-pressible formulation µ is kept constant. Moreover, in the comincom-pressible approach K depends on temperature such that the ratio of K and µ is constant. In this formulation we do not con-sider the thermal conductivity of a mixture since vapor mass fractions are typically small. In the incompressible approach the thermal conductivity of a mixture is incorporated, however, thermal conductivities of water vapor and pure air are considered constant and not depen-dent on temperature. In addition, the thermal diffusivity D depends on temperature in the compressible formulation such that the ratio ρ D

µ is constant, while in the incompressible

for-mulation D is considered constant. All these are small differences: the typical mass fraction of water vapor contributes on the order of 1% to the mass density of the carrier phase, and the variations in temperature throughout the domain are on the order of 1 − 2% in the test cases considered. We performed a simulation with exactly the same model for the carrier gas in the compressible formulation as used in the incompressible. The found differences in gas and

(25)

droplet quantities are on the order of 0.001% of the mean values. As a consequence, these small differences between the models do not appreciably influence the level of agreement between the two formulations.

In this study we will compare results from the incompressible carrier phase with those ob-tained using a varying mass-density of the carrier gas in the compressible formulation. As a result we may then also appreciate under which physical circumstances the compressible formulation deviates significantly from the incompressible model.

2.2.2 Equations for droplets

The dispersed phase of the system consists of droplets which are distributed over the total volume of the channel. They are assumed to collide elastically with the solid walls. In these collisions no heat is transferred. Periodic boundary conditions are applied to the droplets: if a droplet leaves the domain, it is reinserted at the opposite boundary with the same properties. Referring to the classification proposed by Elghobashi (1994), we focus on droplet volume fractions high enough so that two-way coupling is required, while the volume fraction is low enough for the effects of collisions between droplets to be negligible. This regime corre-sponds to droplet volume fractions from 10−6to 10−3. In the current study the droplet volume fraction fluctuates around 10−4.

For both the compressible and the incompressible model the same system of ordinary differ-ential equations for the dispersed phase is adopted and this will be discussed in detail in this subsection. The model for the dispersed phase uses a point-particle approach. We assume that the droplets are spherical and include evaporation and condensation to describe the response of the droplets to the thermodynamic conditions found in the carrier gas along their trajecto-ries.

The droplets are tracked individually in a Lagrangian manner. We solve a system of ordi-nary differential equations for position, velocity, temperature and mass of each droplet. The location of a droplet is governed by the kinematic condition:

dxi(t)

dt = vi (2.15)

where xiis the location and viis the velocity of droplet i.

For the water-air system droplets have a much higher mass density than the carrier fluid. The Stokes drag force is the dominant force in these circumstances (Armenio and Fiorotto 2001). We do not consider gravity acting on the droplets in the current study and solve the equation of motion for each droplet:

dmivi dt = mi u (xi,t) − vi τd,i 1 + 0.15Re0.687_d,i + v (xi,t) dmi dt (2.16) where miis its mass and u(xi,t) is the velocity of the carrier gas at the droplet position. In

addition, we define Red,i as the Reynolds number based on the diameter di of the droplet

and on the relative velocity, Red,i= diρ |u(xi, t) − v|/µ. Moreover, τd,i= ρld2i/(18µ) defines

(26)

flow. The first term on the right-hand side in (2.16) is the drag force for which the Schiller-Naumann correlation is used (Schiller and Schiller-Naumann 1933). This correlation was reported accurate for Red,i between 0 and 1000, which range of conditions is fully adequate for our

simulations. The second term on the right-hand side of (2.16) reflects the contribution from the momentum of the vapor in case of evaporation or condensation of water vapor. The vapor has the same velocity as the droplet at the moment it condenses or it has just been evaporated at the droplet surface.

In this study we do not consider nucleation of droplets. In fact, homogeneous nucleation only occurs when the relative humidity is significantly larger than 100%, which in the current study does not occur.

In the current study the assumption of a uniform temperature inside a droplet is used. This can be motivated by the Biot number which is the ratio of the heat transfer resistance inside the droplet to the heat transfer resistance at the surface of the droplet, Bi = hmdi/kdwhere hmand

kddenote the convective heat transfer coefficient from gas to droplet and the droplet thermal

conductivity, respectively. This dimensionless number determines whether the temperature varies significantly within the droplet. One may expect an almost uniform temperature to good approximation if Bi 1 (DeWitt et al. 2007). In the current study we use the following correlation for hmfor small droplets, valid for low values of Red,i(Bird et al. 1960):

hmdi

K = 2 + 0.6 Re

1/2

d,i Pr1/3 (2.17)

In actual simulations we found that the pdf of Red,icontains most of its statistical mass in

the range Red,i≤ 1, with some low-probability events up to Red,i= 5. The values of Red,i

are contained in such narrow range in these simulations since gravity was not included. The resulting Biot number is therefore only slightly higher than the ratio of the thermal conduc-tivity of the carrier gas and the droplet, Bi = K/kd= 0.046; this motivates the assumption of

a uniform temperature, Ti, inside a droplet.

The total energy of droplet i is given by Ei=12mi|vi|2+ clmiTi, where cl is the specific heat

of liquid water. This energy may change because of two mechanisms: (i) condensation and evaporation and (ii) the heat exchange by convection at the droplet surface which arises due to the temperature difference between the surroundings and the droplet surface:

d dtEi= hv

dmi

dt + hmAi(T (xi,t) − Ti) (2.18) For the convective contribution, we will use the correlation for forced convection around a sphere (Bird et al. 1960) under the assumptions of uniform surface temperature and small net mass-transfer rate. Under these circumstances we can adopt the expression for the heat transfer coefficient hmas in (2.17). In addition, Aiin the convective term denotes the surface

area of the droplet.

The first contribution on the right-hand side of (2.18) contains the specific enthalpy of vapor, where

hv= λ0+ cpvTi (2.19)

Here λ0stands for the latent heat at temperature equal to 0 K and cpvdenotes the specific heat

(27)

In order to close the system of equations we require an expression for dmi/dt, the rate of

evaporation and condensation of a droplet. Assuming a vapor film around the droplets of uniform width δ , we may use the following equation for the change of the droplet mass, mi

(Bird et al. 1960): dmi dt = − miSh 3 τd,iSc ln 1 −Yv,δ ,i 1 −Yv,s,i (2.20)

where Yv,δ and Yv,s are the vapor mass fractions at a distance δ from the surface of the

droplet and on the surface of the droplet, respectively. Here, we adopted the Schmidt number Sc = µ/(ρlD) and the Sherwood number Sh = di/δ (Bird et al. 1960). In addition, for the

Sherwood number we use the correlation for a sphere (Bird et al. 1960):

Sh = 2 + 0.6Re1/2_d Sc1/3 (2.21) The driving quantity in (2.20) is the difference Y_{v,δ ,i}− Yv,s,i. The vapor mass fraction at the

surface, Yv,scorresponds to thermodynamic equilibrium at the surface, which determines the

saturation vapor pressure. The saturation pressure pv,satis calculated using Antoine’s relation

(Antoine 1888): pv,sat= 105exp A− B C+ T_i in Pa (2.22)

where A = 11.6834, B = 4089.59 K and C = 500.02 K. Subsequently, the ideal gas law for water vapor is used to calculate the vapor mass fraction at the surface of the droplet:

Yv,s,i=

pv,sat(Ti)

ρ RwaterTi

(2.23)

In order to find Yv,δ ,i, i.e., the vapor fraction in the close vicinity of the droplet, the assumption

of a point-particle approach motivates to calculate Y_{v,δ ,i} as the vapor mass fraction of the carrier gas at the position of the droplet, i.e., Y_{v,δ ,i}= Yv(xi,t).

2.2.3 Coupling terms

The governing equations for the carrier phase contain two-way coupling terms. These terms appear because of the mass, momentum and energy exchange between the carrier and the dispersed phases. The interaction between the two phases conserves mass, momentum and energy. The contributions from the dispersed phase in the equations of mass, momentum, energy and water vapor for the carrier gas are given by:

Qm= −

∑

i dmi dt δ (x − xi) (2.24) Qmom= −

_∑

i d dt(mivi)δ (x − xi) (2.25)

(28)

Qe= −

_∑

i d dt clmiTi+ 1 2mi|vi| 2 ! δ (x − xi) (2.26)

where the sums are taken over all droplets in the domain. The delta-function expresses that the coupling terms act only at the positions of the droplets, consistent with the underlying point-particle assumption.

The coupling term (2.24) in the mass density equation shows that mass transfer between the phases occurs because of evaporation and condensation. In order to calculate this term equation (2.20) will be used. Momentum transfer between the phases as expressed by (2.25), consists of two mechanisms: the drag force between the droplets and the carrier gas and the momentum transfer due to mass transfer arising from phase changes. In order to calculate this coupling term equation (2.16) will be used, rewritten to provide d(mivi)/dt as required.

Finally, in (2.26) the first term on the right-hand side represents the convective heat transfer and the energy transferred to the carrier phase by vapor. The last term in (2.26) appears as the contribution from the kinetic energy of the droplets. This contribution can be evaluated from the right-hand side of (2.18).

This completes the description of the governing equations for the two phases. In the next section we will describe the numerical approaches used for the treatment of the governing equations for the two phases and the coupling terms.

2.3 Numerical method

In this section the numerical approaches applied to the compressible and incompressible for-mulations will be described. In addition, the numerical methods for the dispersed phase and implementation of the coupling terms will be discussed.

The incompressible and compressible models are simulated using different numerical meth-ods. The incompressible formulation adopts a pseudo-spectral method for the carrier phase as described in Russo et al. (2014) and Kuerten (2006). In the two periodic directions, x and z, a Fourier-Galerkin approach is used, while a Chebyshev-collocation approach is applied in the wall-normal direction, y. The droplet equations are integrated in time with the same three-stage Runge-Kutta method as used for the convective terms in the carrier phase equa-tions. The viscous terms and pressure are integrated in time with the implicit Crank-Nicolson method. The pressure is found at each time level from the divergence-free condition for the velocity, following the approach proposed by Kleiser and Schumann (1980).

The numerical model for the compressible carrier phase adopts a second-order accurate ex-plicit four-stage low-storage Runge-Kutta time-stepping with weights equal to 1/4; 1/3; 1/2 and 1. A second order finite volume discretization (Vreman et al. (1992)) is used. The geom-etry is divided into rectangular cells. A uniform grid is used in the two periodic directions. The grid spacing is denoted by ∆ x and ∆ z for the streamwise and the spanwise direction, respectively. In the wall-normal direction a non-uniform grid is applied which is finer near the walls. We use 128 grid points in each direction. This choice for the grid resolution is motivated in detail in Marchioli et al. (2008) where simulation results obtained with various numerical methods were compared. In fact, simulation of a particle-laden turbulent channel flow at Reτ= 150 using second-order central differences, similar to our method, on the same

(29)

gas and particle statistics with other studies based on pseudo-spectral methods. In fact, an agreement on the order of 1% in the mean velocity profile and on the order of 2% in RMS of fluctuating quantities was observed. This demonstrates the reliability of the DNS results reported in this paper and shows that the basic mechanisms are captured well enough to dis-tinguish compressibility effects in multiphase flows.

The variables are stored in the centers of the cells. For each of the cells we integrate the gov-erning equations (2.1) - (2.4) for the carrier phase in conservative form over the volume of the cell. For the integrals of the coupling terms we proceed as follows. As an example, we consider the integration over a grid cell of the coupling term for the mass density equation. Substituting the explicit expression for it, (2.24), we obtain:

Z Vcell QmdV= −

_∑

i Z Vcell dmi dt δ (x − xi)dV = −_i∈cell

∑

dmi dt (xi) (2.27)

where the last sum is taken over all droplets within the cell. At any time, a droplet contributes only to the flux in the grid cell where it is located. The size of the droplet is everywhere smaller than the smallest grid size. The smallest grid size in wall units is equal to 0.6 while the droplet diameter in wall units equals to 0.5.

The system of ordinary differential equations for the droplets is integrated in time applying the same four-stage compact Runge-Kutta method as used for the gas equations. The right-hand sides of the droplet equations contain gas properties at droplet location. In order to determine these, tri-linear interpolation is applied (Marchioli et al. 2008).

2.4 Initial conditions of the simulations

In this section the initial conditions for the different test cases will be described and moti-vated, both for the compressible and the incompressible formulations.

In this study two sets of simulations are chosen. The first set considers applying three dif-ferent values of the Mach number, 0.2, 0.1 and 0.05 to a system under atmospheric pressure and room temperature. We will study convergence of the results to the incompressible results when decreasing the Mach number. The value of the Mach number based on the reference conditions is about 5 × 10−3. A Mach number this small leads to a further reduction of the time step for the chosen explicit numerical method, which would render the simulations in-feasible.

In all the cases the simulations were started from a turbulent velocity field, obtained from a simulation without droplets at the correct values of the Mach number and the Prandtl number. These simulations were done with adiabatic boundary conditions at the walls until the statis-tically steady state was reached. The resulting temperature profiles are shown in Figure 2.2 for the three different values of the Mach number as well as for the incompressible case for which the initial temperature is uniform. The initial mean temperature for these simulations was chosen equal to 293.15 K. The highest Mach number leads to variations in temperature on the order of 1%. It is seen that reduced compressibility implies that the initial conditions for the compressible cases converge to those of incompressible flow.

(30)

Rela-tive humidity φ characterizes the amount of vapor in the channel. It is the ratio of the partial pressure pv of water vapor in the air-vapor mixture to the saturated vapor pressure at the

prescribed temperature, φ = pv/pv,sat. The condition of φ = 100% was selected in order to

have very ‘mild’ coupling conditions for the turbulent flow in the channel as far as mass transfer is concerned, resulting in the smallest differences between the compressible and the incompressible formulations. From the state equation for water vapor we find the vapor mass fraction in the channel:

Yv=

pv,sat(Tre f)

ρ RwaterTre f

(2.28)

where pv,sat(T ) follows from (2.22). Subsequently, the total energy density is found taking

into account the contribution from the added water vapor according to (2.9).

We randomly distribute Ndrop= 2, 000, 000 identical droplets over the volume of the channel

and apply a positive heat flux at the top wall and a negative flux of equal magnitude at the bottom wall. For both the incompressible and the compressible simulations the initial droplet diameter is given by di/H = 3.09 × 10−3. In this way, the initial volume fraction of droplets

is on the order of 10−4. The Kolmogorov length varies along the wall-normal direction from a minimum value of 1.6 near the wall to a maximum value of 3.6 in wall units at the centerline (Marchioli et al. (2008)). Consequently, the maximum ratio of the droplet diameter, which is in wall units equals to 0.5, and the Kolmogorov scale is equal to approximately 0.3. This ratio shows that it is allowed to use the point-particle approach, (Marchioli et al. (2006)), (Elghobashi (1994)). Velocity and temperature of the droplets are initialized using the carrier gas values at the droplet locations.

The second set of simulations reflects physical situations at φ = 100% for which the effect of compressibility is expected to be more important. We will investigate the effects of a higher heat flux through the walls and a higher mean temperature. In the first case we will apply a five times larger heat flux through the walls and in the second the mean temperature will be increased to 323.15 K. Both these changes are quite significant while keeping the flow and process regime comparable, e.g., not inducing effects of boiling. A description of all test cases can be found in Table 2.2.

The slightly higher heat flux in case 5 of about 9% is chosen in order to obtain the same tem-perature gradient at the walls as in cases 1-3, compensating for the temtem-perature dependence of the thermal conductivity of the carrier gas, which, according to Table 2.1 changes by the same amount. In the next section we will discuss the results obtained from these two sets of simulations.

case number Ma Q˙(W/m2) Tmean(K)

1 0.2 31.8987 293.15

2 0.1 31.8987 293.15

3 0.05 31.8987 293.15

4 0.05 159.49 293.15

5 0.05 34.8612 323.15

(31)

−0.02 −0.01 0 0.01 0.02 292.5 293 293.5 294 294.5 295 295.5 y [m] <T> [K]

Fig. 2.2: : Initial fluid temperature averaged over the homogeneous directions as a function of the wall-normal coordinate; dashed: incompressible, circles: Ma = 0.2; dash-dotted: Ma = 0.1; solid: Ma = 0.05.

2.5 Results

In this section the results obtained with the compressible and incompressible formulations will be compared. The main purpose of the comparison is to see whether and which com-pressible results converge to the incomcom-pressible results if we decrease the value of the Mach number. We concentrate on fluid flow aspects, heat transfer properties of the system and on the behavior of the droplets and adopt a spatial resolution of 128 × 128 × 128 with both the finite volume and the spectral method. Four subsections will be presented: in 2.5.1 the agree-ment in the kinematics of the flow between the two formulations is shown, 2.5.2 is dedicated to thermal properties of the carrier phase, in 2.5.3 results for the water vapor are presented and finally, in 2.5.4 we present results of the dispersed phase.

2.5.1 Velocity properties

In this study we investigate compressibility in the form of changes in the mass density due to phase transition taking liquid water from the dispersed water droplets to the vapor and vice versa. This compressibility is expressed by a non-zero divergence of velocity and a non-zero gradient of the mass density. In order to quantify this explicit compressibility we compared the RMS of the term ρ∇ · u with the RMS of the term u · ∇ρ in the continuity equation. In Figure 2.3 we display these as functions of the wall-normal coordinate. We observed that both during the initial stages of the simulations as well as in a well developed stage the RMS of these two terms is of the same order of magnitude, accounting for changes in the local density on the order of a few percent. At initial conditions with stronger heat flux and higher

(32)

mean system temperature these explicit compressibility effects are more pronounced. First, we establish the level of agreement in the kinematics of the turbulent flow in both formulations. In Figure 2.4a the streamwise component of the carrier gas velocity averaged over the periodic directions and over time is presented. The maximum difference between the two formulations at Ma = 0.2 is on the order of 2%. In addition, the bulk Reynolds number history as defined in (2.10) is presented in Figure 2.4b. The difference is on the order of 1%. For both formulations the bulk Reynolds number increases in time. The forcing of the flow is chosen in such a way that the total momentum of the system, the sum of the momentum of the carrier and dispersed phases, is conserved. Droplets tend to migrate to the walls (see section 2.5.4), where they will gradually obtain a lower streamwise velocity. That implies that the carrier gas will gain momentum and this causes the gradual increase of the bulk Reynolds number. At reduced Ma the correspondence between the two formulations becomes even closer and the maximum difference in the averaged streamwise velocity is about 1% at Ma = 0.1 and Ma = 0.05. This also establishes the level of correspondence between the two independent direct numerical simulations at a resolution of 1283, employing very different numerical methods for spatial discretization as well as time integration.

−0.020 −0.01 0 0.01 0.02 0.02 0.04 0.06 0.08 0.1 0.12 0.14 y [m] [kg/s m 3 ] (a) −0.020 −0.01 0 0.01 0.02 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 y [m] [kg/s m 3] (b)

Fig. 2.3: Solid: RMS of ρ∇ · u, dashed: RMS of u · ∇ρ as a function of the wall-normal coordinate. (a): at 0.02s, (b): at 4s .

2.5.2 Heat transfer properties

We apply a heat flux to the walls, which induces a temperature gradient in the wall-normal direction: a region with higher temperatures forms near the warmer top wall and a region with lower temperatures near the colder bottom wall. During the initial transient period the carrier phase adapts to the applied boundary conditions and an asymmetric profile of the gas temperature develops from the initial nonuniform temperature field, cf. Figure 2.2. The evo-lution of the profile of the mean temperature, averaged over the homogeneous directions, is

(33)

−0.020 −0.01 0 0.01 0.02 0.5 1 1.5 2 2.5 y [m] <u x > [m/s] (a) 0 2 4 6 8 2280 2300 2320 2340 2360 2380 2400 2420 2440 2460 t [s] Re b (b)

Fig. 2.4: (a) Streamwise velocity component averaged over the homogeneous directions and over time [2s; 8s] as a function of the wall-normal coordinate. (b): The bulk Reynolds num-ber history. Compressible simulations (solid lines) are obtained at Ma = 0.2, incompressible model is shown with dashed lines, following the settings of Case 1.

presented in Figure 2.5 for two different values of the Mach number. It is seen that reduced compressibility implies a more symmetric profile of the mean temperature at later times. In Figure 2.6 the asymptotic mean temperature profile is presented for different values of the Mach number and for the incompressible formulation. The carrier gas reaches a statistically steady state for which the mean temperature profiles can be determined after approximately 2s. The results are averaged over time interval [2s; 8s] and over the homogeneous directions. From figure 2.6 it can be concluded that a reduced compressibility in the compressible for-mulation leads to a strongly increased agreement with the results obtained with the incom-pressible model. Averaging over a longer time interval does not change the conclusions. In Figure 2.6 one can observe significant differences between the temperature profiles for Ma = 0.2 and for Ma = 0.05. These differences are characteristic for the superposition of two effects. On the one hand, in the statistically steady state the adiabatic simulations yield a non-uniform temperature profile. For the larger Mach number the temperature profile is more non-uniform while remaining symmetric. This reflects the property that not the mean temperature is constant in compressible turbulent channel flow, but rather the stagnation tem-perature. On the other hand, the application of the heat flux to the walls at Ma = 0 leads to an anti-symmetric temperature profile. Combined, the symmetric (no heat flux) and antisym-metric (Ma = 0) profiles lead to the characteristic profile shown. This superposition is present at all Mach numbers, but at Ma = 0.2 both effects are of the same magnitude, explaining the apparent qualitative difference.

Figure 2.7a shows how the root-mean square of the gas temperature evolves at Ma = 0.2. The maximum value of oscillations is equal to approximately 0.4 K which is around 10% of the temperature difference between the two walls.

Figure 2.7b shows the agreement of the root-mean square of the gas temperature for the in-compressible and in-compressible formulations. It can be seen that the agreement between the two formulations increases if the Mach number is reduced. The level of variation in Trmsat

(34)

low Mach numbers, compared to the incompressible results is expected to represent slight differences resulting from differences in the numerical methods at the selected spatial resolu-tion - currently, further grid refinement is not feasible.

To express the efficiency of the heat transfer between the channel walls we consider the Nus-selt number. It has been shown that the NusNus-selt number can be increased by more than a factor of two if solid particles with high heat capacity are added to a turbulent channel flow (Kuerten et al. 2011). The Nusselt number is defined in the following way:

Nu = dTg dy    wall _. ∆ Tg 2H (2.29)

where bars denote averaging over the two homogeneous directions, ∆ Tg is the difference

in gas temperature between the two walls and the derivative with respect to the wall-normal coordinate in the numerator is the average over both walls. Figure 2.8a shows how the Nusselt number depends on time for both formulations. We observe a close agreement with a relative difference of up to about 10%. The compressible formulation predicts a slightly higher value than the incompressible model.

In the initial phase of the simulations the compressible and incompressible formulation result in differences in the value of the Nusselt number up to 3. Kuerten et al. (2011) showed that for incompressible flow with solid particles three contributions to the Nusselt number can be distinguished:

Nu = Nu_lam+ Nuturb+ Nupart (2.30)

which correspond to the Nusselt number for laminar flow, a contribution from turbulence and a contribution from particles. For turbulent flow with droplets there are some additional contributions, but they can be shown to be negligible (Russo et al. 2014b). An important difference between compressible and incompressible flow is that Nuturbin compressible flow

can be written as:

Nuturb= − − 1 α ∆ T Z H −H u2T (2.31)

can be written as:

Nuturb= − 1 α ∆ T Z H −Hu 0 2T0− 1 α ∆ T Z H −Hu2T (2.32)

where · means averaging over the periodic directions and α is the thermal diffusivity. For incompressible flow the second contribution in (2.32) is equal to zero, because the mean of the wall-normal component of the velocity equals zero in incompressible channel flow. Our simulation results for the compressible formulation show that this contribution of the mean wall-normal velocity component is larger than the other contribution by a factor of 3 in the initial phase of case 3. This explains the initial difference in Nusselt number shown in Figure 2.8a. The presence of the additional term in Nuturb in the compressible formulation

reflects the additional net transfer of the gas from the hot wall to the cold wall during the initial stages which increases the value of the Nusselt number.

Figure 2.8b shows the results of the two formulations for cases 3, 4 and 5. The Nusselt number obtained from the compressible simulations is slightly higher than the corresponding

(35)

incom-pressible result for all three cases considered. The case with the higher mean temperature gives a significant increase in the value of the Nusselt number by more than a factor of 2. This is related to the difference in mean temperature between the two walls, figure 2.9. In case 5 the temperature difference between the two walls is smaller by more than a factor of two than in case 3. In the statistically steady state of cases 3 and 5 the droplets near the walls have the same diameter, see Section 2.5.4, but in the case of a higher mean temperature the mass transfer rate from the droplets to the carrier gas is higher (Russo et al. 2014), adding correspondingly to the transfer of heat. This explains the higher Nusselt number in case 5. Figure 2.8b shows that the compressible formulation becomes more important in case 5. The results for case 3 and 4 show very close agreement - the solid lines and square-marked lines practically coincide in Figure 2.8b. The independence of the Nusselt number on the value of the applied heat flux is further illustrated by the temperature difference between the walls: this difference is approximately 5 times larger in case 4 than in case 3, Figure 2.9. Thereby the increase of the heat flux through the walls by a factor 5 is almost completely compen-sated by the same factor of increase in the temperature difference, hence leading to a nearly unchanged Nu in case 4, compared to the reference case 3.

We also performed simulations with the two codes for the case of solid particles and com-pared the agreement in the Nusselt number history for these simulations with the case of droplets, for case 5 of 323.15 K mean temperature. In case of particles dispersed in the flow, changes in mass density of the carrier phase due to evaporation and condensation are absent and the agreement between the compressible and the incompressible formulation is expected to be closer than in case of droplets. Although we investigate in this paper situations with rel-ative humidity of 100% for which the mass transfer between the liquid and gaseous phases is quite small, we do observe that the agreement between the two formulations is indeed better for solid particles than for evaporating droplets. In Figure 2.10 we compare the Nusselt num-ber obtained with the two formulations and confirm that the differences in case particles are dispersed are considerably smaller than in case of droplets. As may be inferred, the Nusselt number in the well-developed stages the droplets cases differ in absolute value by about 3-4, while the differences in the particle cases is reduced to well bellow 1.

The viscous heating term which is present in the compressible formulation and absent in the incompressible transfers kinetic energy into heat. In case when all kinetic energy would be transferred into heat, the gas temperature would change on the order of 0.05%. We also per-formed a simulation without the viscous heating term in the compressible model. The results of this simulation demonstrate that the viscous heating term is negligible when it comes to assessing the agreement between the two formulations.

2.5.3 Water vapor

In this subsection characteristics of the water vapor distribution will be shown. The temper-ature gradient in the wall-normal direction causes evaporation of droplets in the region near the hot wall and condensation of vapor on the droplets in the region near the cold wall. To understand this more precisely we consider the region near the hot wall. The temperature increase in this region causes an increase of the saturation pressure according to Antoine’s relation. Consequently, the vapor mass fraction at the surface of the droplet increases as well and this leads to evaporation. The opposite situation, i.e., condensation of water vapor, is

(36)

−0.02 −0.01 0 0.01 0.02 292 292.5 293 293.5 294 294.5 295 295.5 296 296.5 y [m] < T > [K ] (a) −0.02 −0.01 0 0.01 0.02 291.5 292 292.5 293 293.5 294 294.5 295 y [m] < T > [K ] (b)

Fig. 2.5: Gas temperature averaged over the homogeneous directions as a function of the wall-normal coordinate at different times. Dashed lines: initial profile obtained without droplets, solid lines: evolution of the profile in time, shown at t = n∆t with ∆t = 1s. The arrows show the direction of increasing time. (a): Ma = 0.2, (b): Ma = 0.05.

−0.02 −0.01 0 0.01 0.02 291 291.5 292 292.5 293 293.5 294 294.5 295 295.5 y [m] <T> [K]

Fig. 2.6: Gas temperature averaged over the homogeneous directions and over time [2s; 8s] as a function of the wall-normal coordinate. Dashed: incompressible, circles: Ma = 0.2; dash-dotted: Ma = 0.1; solid: Ma = 0.05.

obtained in the region near the cold wall. This creates a vapor mass density gradient in the wall-normal direction of the channel, as is illustrated in Figure 2.11. Moreover, both water vapor mass density, Figure 2.11a, and its root-mean square, Figure 2.11b, show that decreas-ing compressibility through a decrease of Ma leads to increased correspondence with the incompressible results.

Heat and mass transfer in turbulent multiphase channel flow