DNS and Highly-Resolved LES of Heat and Mass Transfer in Two-Phase Counter-Current Condensing Flow

Abstract

A comprehensive study of direct-contact condensation heat transfer for turbulent, counter-current, liquid/vapour flow in a nearly horizontal channel at high pressure (i.e. 5 MPa) has been carried out based on Direct Numerical Simulation (DNS) and highly-resolved Large Eddy Simulation (LES) approaches. To simulate the two-phase flow situation, driven in this case by a constant pressure gradient, a single set of Navier–Stokes equations, coupled with an enthalpy conservation equation, have been employed. The interfacial mass transfer, seen in this case to be dominated by condensation, has been calculated directly from the heat flux at the liquid/vapour interface. To investigate the effect of condensation on the turbulence phenomena, and vice versa, cases have been considered involving two friction Reynolds numbers: namely $Re_{*}=u_{*}h/\nu =178$ and $Re_{*}=u_{*}h/\nu =590$ ($u_{*}=(h\varDelta P/\rho )^{1/2}$). At the lower Reynolds number, three levels of water subcooling—0 K, 10 K and 40 K—have been investigated. The use of water subcooling of 0 K has enabled the validation and verification procedures associated with the numerical approach to be compared against experimental and numerical data reported in the literature. The choice of the maximum degree of water subcooling is dictated by the need to justify the periodic boundary conditions applied in this numerical study. In the simulation for the higher Reynolds number, only the case of 10 K subcooling has been included, as a consequence of the very high computation effort involved. A detailed statistical analysis of the DNS and LES data obtained from the application of the well-known wall laws has also been assessed. In the vicinity of the liquid/vapour interface, the characteristics of the turbulent motions appear somewhat diverse, depending on whether the interface is basically flat or wavy in character. For a flat interface, some damping effect of the presence of the interface on the turbulence intensity has been observed, a feature which becomes enhanced as the level of liquid subcooling is increased. In the case of a wavy interface, the damping effect is predicted as considerably less pronounced.

1 Introduction

Direct Contact Condensation (DCC) denotes the heat transfer process, including phase change, that occurs between two or more fluids at any common free interface between them. DCC is an important process in a multitude of diverse industrial applications, such as in cooling towers, direct contact condensers, contact feed-water heaters, de-aerators and water desalination units. Condensation of vapour on a cold liquid surface is also of major importance for safety analyses in the context of nuclear reactor systems, in particular during a (postulated) two-phase Pressurised Thermal Shock (PTS) event (Lucas et al. 2009a, b), in which cold Emergency Core Cooling (ECC) water is injected into a feed-water pipe of the primary circuit of a Pressurized Water Reactor (PWR) and partially mixes with the hot coolant which is still present. Another safety relevant phenomenon involving DCC is condensation-induced water hammer (CIWH) in piping systems, which may lead to destructive pressure overloads (Kirsner 1999; Barna et al. 2010; Strubelj et al. 2010). For example, several such events have been recorded during bitumen production at the Athabasca Oil Sands in Canada, where the Steam-Assisted Gravity Drainge (SAGD) technique is employed (ERCB 2008; Urban and Schlüter 2014).

Computational Fluid Dynamics (CFD) approaches, which are, in principle, able to take into account the 3-D nature of the phenomena, are currently being used to estimate the damage potential (Strubelj et al. 2010; Ceuca and Macian-Juan 2012) in such events. The essential closure model for the simulation of such DCC phenomena under turbulent flow conditions is represented by the overall liquid heat transfer coefficient, $HTC_L$. This parameter is essential in estimating the interfacial heat flux, from which the mass transfer rate is calculated, and hence upon the pressure propagation characteristics. Current practice is based on the use of empirical correlations, obtained for a specified geometry, and within a limited range of operating conditions (i.e. as determined by the pressure, temperature, etc.). The application of such correlations to other situations for which they have been formulated is not always straightforward, could be questionable, or even impossible to substantiate (Boehm and Kreith 1988). Moreover, the correlations are usually expressed in terms of integral quantities, such as the bulk liquid velocity and average temperature, parameters which are not obviously appropriate for the 3-D CFD modelling of $HTC_L$.

Over the last few decades, multiple experimental studies on direct-contact condensation heat and mass transfer under two-phase flow conditions have been performed, and a considerable number of empirical correlations for the all-important condensation heat transfer coefficient have been proposed (Linehan et al. 1970; Segev et al. 1981; Kim and Bankoff 1983; Bankoff and Lee 1987; Lee et al. 2006). All correlations are expressed in the form of a Nusselt number $Nu= HTC_L L/ \lambda _L$ (where L is a characteristic length scale and $\lambda _L$ is the thermal conductivity); the Nusselt number expresses the ratio of convective to conductive heat transfer. Typical correlations are expressed as follows:

$$\begin{aligned} Nu=C_1Re_t^{C_2}Pr^{C_3} \quad \text{ or } \quad Nu=C_4Re_L^{C_5}Re_V^{C_6}Pr^{C_7}, \end{aligned}$$

(1)

in which Pr is the molecular Prandtl number, and $C_1$—$C_7$ are model constants, all to be determined empirically.

The constants represented in the expression (1) vary considerably from one author to another, as one might expect (Fulgosi 2004). It may be supposed that the disagreement reflects a general lack of understanding of the fundamental liquid-side interfacial mass, momentum and energy transport mechanisms. Perhaps part of the scatter in the data may be attributed to often undisclosed experimental uncertainties in measuring temperatures and turbulence parameters near to or at the interface (Celata et al. 1989), or just a fundamental lack of available measurements or uncertainty qualification in this critical region of the flow field (Peturaud et al. 2011). Certainly, the presence of high interfacial and wall shear stresses, together with turbulent motions at small scales, make measurements close to walls and interfaces intrinsically very difficult in many circumstances (Kim et al. 1987). To emphasise the point, Apanasevich et al. (2014) has presented pre-test blind simulations of a steady-state TOPFLOW-PTS, steam/water test with condensation. The TOPFLOW-PTS test facility of the Helmholtz-Zentrum Dresden-Rossendorf (HZDR) centre was specifically constructed to investigate such two-phase PTS situations (Peturaud et al. 2011). Three CFD codes ANSYS CFX (ANSYS 2018a), ANSYS FLUENT (ANSYS 2018) and NEPTUNE_CFD (Méchitoua et al. 2003) were featured in this benchmark analysis. Notwithstanding the different approaches, and software options, considerable discrepancies in the numerical estimates of the overall heat transfer were predicted. Granted that subsequently to these analyses progress has been achieved in the modelling of some of the key flow phenomena (Lakehal and Labois 2011; Coste 2013), it appears that CFD modelling of two-phase flows with phase change still remains far from satisfactory, due both to the limited theoretical knowledge of the actual DCC process, the lack of CFD-grade experimental data, and a common approach to predicting such data.

Direct Numerical Simulation (DNS) of turbulent, single- and two-phase flows at the micro-scale has become a newly discovered approach to extend the knowledge base of our understanding of the governing physical phenomena, and, based on this, to derive improved correlations to be applied at the meso-scale (i.e. between macro and micro scales). DNS has already been successfully applied to study the fundamentals of single-phase flow situations, including near-wall turbulence and turbulent heat transfer (Kawamura et al. 1998; Tiselj et al. 2001; Hoyas and Jiménez 2006; Bergant and Tiselj 2007). Komori et al. (1993) performed a DNS study of three-dimensional, open-channel flow, involving a zero-shear, gas/liquid interface using a Boundary Fitting Method (BFM) to identify the location of the interface, and were able to reproduce most of the experimentally observed turbulence quantities. Specifically, it was recognized that the relationship between scalar transport of physical quantities across a phasic interface, together with the turbulent motions occurring near the interface, should together be studied in detail. In their DNS computations of open-channel flows, Lam and Banerjee (1992), Pan and Banerjee (1995) and De Angelis (1998) treated the liquid/vapour interface as a rigid, free-slip wall. However, this simplistic approach restricts the application areas to situations in which interface deformations, and their contributions to the transfer processes, can be neglected, and this is not always the case, such as in the case of stratified wavy flows (Cohen and Hanratty 1968; Nordsveen and Bertelsen 1997). It should be noted that in these studies only one phase (liquid) was taken into consideration: this is a restriction in itself. In particular, any possible interaction between the turbulent structures generated separately in the gas and liquid phases on the two sides of the interface would not be represented, nor the interaction between them. To overcome this limitation, Lombardi et al. (1996) suggested another simulation approach (based on BFM) that would specifically account for the coupling of the turbulent motions between the phases, and successfully applied their model to a counter-current air/water flow situation, in a stably stratified configuration, with a flat interface between the phases. Fulgosi et al. (2003), following up on this earlier initiative, used this same coupled simulation technique to investigate turbulence near an air/water interface under conditions in which capillary waves would be present. Results create the impression that the gas phase “sees” the interface as a rigid wall, while, for the liquid phase, the interface appears to act like a free-slip boundary. These are important observations in regards to the directions future numerical developments should take.

Three-dimensional DNS studies of two-phase flows involving phase change (due to condensation or vaporisation), in which interfacial heat and mass transfer are to be directly simulated rather than be modelled by empirical correlations, are somewhat rare. Such simulations make high demands on available computer hardware, since important physical phenomena need to be represented (or modelled) at very different scales (Hoyas and Jiménez 2006). Fulgosi (2004) proposed that the energy equation be added to the method of Lombardi et al. (1996), and performed DNS computations of condensing, counter-current, stratified water/steam flow in a rectangular channel (at constant Reynolds number) to support the argument. Unfortunately, details of the liquid temperature field were not reported, accurate representation of which is very important for model development. Further, it should be emphasized that the use of BFM is restricted to interfaces for which it is assumed that strong topology changes do not take place, or at least have negligible effect on the phasic exchange processes.

This particular shortcoming can be overcome by the use of direct interface-capturing techniques. Recently, Sato and Ničeno (2013) have developed a three-dimensional, mass-conservative, phase-change model for such direct interface capturing, and implemented the model into the in-house 3-D CFD code PSI-Boil: an acronym of Parallel SImulator of BOILing phenomena (Sato and Ničeno 2012). In this approach, the mass transfer rate is calculated directly from the separate heat fluxes arriving at the interface from the liquid and vapour phases. Another key component of the model is the capturing of the sharpness of the mass transfer rate at the liquid/vapour interface: as modelled, the phase change is deemed to occur only in those computational cells directly intersected by the interface. This strategy enables the velocity jump across the interface to be localised, while maintaining the accuracy in the transfer processes. For further discussion of two-dimensional simulation approaches to phase-change phenomena, the interested reader is referred, for example, to the papers of Jamet et al. (2001), Tomar et al. (2005) and Badillo (2012).

The objective of this work is to continue the investigation of turbulent, counter-current, condensing liquid/vapour flow in a horizontal channel, but at high pressure, namely 5 MPa. To this purpose, 3-D DNS and highly-resolved Large Eddy Simulation (LES) calculations have been undertaken using the CFD code PSI-Boil. To overcome the limitations of the numerical study of Fulgosi (2004) the two-phase flow configuration considered here includes both the interfacial and wall regions on both sides of the phasic interface. The reported DNS studies of turbulent, adiabatic (De Angelis 1998; Lombardi et al. 1996) and condensing (Fulgosi 2004) counter-current liquid/vapour flows have been performed only for a single friction Reynolds number $Re_*< 180$, for which significant low Reynolds number effects are apparent. In the present study, two friction Reynolds numbers, $Re_{*}=178$ and $Re_{*}=590$, are considered, specifically to investigate the effect of condensation on the turbulence levels, and vice versa.

In the following sections, we will present the mean statistical properties of turbulence in the interfacial and wall regions for both phases, and the associated liquid temperature field. By this means, we intend to contribute to the ongoing discussions regarding the interplay between the turbulence levels and the phase-change phenomena as direct functions of the Reynolds number and the degree of liquid subcooling. Due to the lack of CFD-grade experimental data on DCC, specifically regarding velocity and temperature measurements close to, and at, the liquid/vapour interfaces, the newly established database may be of considerable value for the development and testing of turbulence and heat transfer closure models in the general context of two-phase interacting flows.

2 Numerical Procedure

The numerical method used for the simulation of two-phase flow situations is based on the finite-volume approach, within a Cartesian grid system, in a staggered-variable arrangement, which dates back to the pioneer work of Harlow and Welch (1965). For the coupling of the velocity and pressure fields, the projection method of Chorin (1968) is incorporated. The interface between the two phases is simulated using an interface capturing procedure proposed by Sato and Ničeno (2012). Within the two-phases methodology, a single set of Navier–Stokes equations describing the fluid flow conditions is coupled with an enthalpy conservation equation governing the heat transfer between them, and a sharp-interface model to represent the phase-change processes taking place at the liquid/vapour interfaces themselves.

2.1 DNS Method

The conservation equations for mass, momentum and enthalpy are generally expressed in the following forms (Sato and Ničeno 2013):

$$\begin{aligned}&\nabla \cdot \mathbf {U} = \varGamma \left( \frac{1}{\rho _V}-\frac{1}{\rho _L} \right) , \end{aligned}$$

(2)

$$\begin{aligned}&\frac{\partial \left( \rho \mathbf {U} \right) }{\partial t} + \nabla \cdot \left( \rho \mathbf {U} \mathbf {U} \right) =-\nabla P + \nabla \cdot \mathbf {\tau } + \rho \mathbf {g} + \mathbf {f}, \end{aligned}$$

(3)

$$\begin{aligned}&\rho c_p \left( \frac{\partial T}{\partial t} + \mathbf {U} \cdot \nabla T \right) = \nabla \cdot \left( \lambda \nabla T\right) , \end{aligned}$$

(4)

in which $\mathbf {U}$ is the flow velocity vector (m/s), $\rho$ the density (kg/m$^3$), t the physical time (s), P the pressure (Pa), $\mu$ the molecular dynamic viscosity (Pa s), $\tau$ the molecular viscous stress tensor (Pa), $\mathbf {g}$ the gravitational acceleration (m/s), $\mathbf {f}$ any appropriate body force (N/m$^3$), $\varGamma$ the overall mass transfer rate (kg/m$^3$ s), $\lambda$ the thermal conductivity (W/mK), and $c_p$ the specific heat capacity at constant pressure (J/kg K). The body force $\mathbf {f}$ incorporates the surface tension force at a liquid/vapour interface and is here defined in terms of the Continuum Surface Force (CSF) model of Brackbill et al. (1992) as:

$$\begin{aligned} \mathbf {f}=\sigma \kappa _\sigma \nabla \phi , \end{aligned}$$

(5)

where $\sigma$ is the molecular surface tension coefficient (N/m), $\kappa _\sigma$ represents the local curvature (1/m), and $\phi$ is the colour function, defined in this formulation as the liquid volume fraction inside a particular cell containing the interface. The curvature $\kappa _\sigma$ is calculated from the signed distance function d instead of the color function, because the curvature calculation from the color function tends to include numerical error Yabe et al. (2007):

$$\begin{aligned} \kappa _\sigma =-\nabla \cdot \left( \frac{\nabla d}{|\nabla d|} \right) . \end{aligned}$$

(6)

The signed distance function is calculated from the color function using the method proposed by Russo and Smereka (2000). The liquid volume fraction is commonly used to define the mean density and viscosity of the fluid occupying the cell as follows:

$$\begin{aligned} \rho= & {} \phi \rho _L + \left( 1- \phi \right) \rho _V, \end{aligned}$$

(7)

$$\begin{aligned} \mu= & {} \phi \mu _L + \left( 1- \phi \right) \mu _V, \end{aligned}$$

(8)

in which the subscripts L and V here denote the liquid and vapour phases, respectively. In contrast to the definition of the density and viscosity, the thermal conductivity and the specific heat capacity are not considered to linearly depend on the colour function but are defined as follows:

$$\begin{aligned} \lambda = \left\{ \begin{array}{ll} \lambda _L, &{} \phi \ge 0.5 \\ \lambda _V, &{} \phi< 0.5 \end{array} \right. , \quad c_p = \left\{ \begin{array}{ll} c_{p,L}, &{} \phi \ge 0.5 \\ c_{p,V}, &{} \phi < 0.5 \end{array} \right. . \end{aligned}$$

(9)

This definition of the molecular coefficients in the energy equation allows the jump condition in the fluid properties to be captured during the numerical solution procedure.

The mass and momentum equations (2, 3) are solved within the code PSI-Boil using the semi-implicit projection method of Chorin (1968). That is, the convective terms are discretised using the two-level, explicit Adams–Bashforth scheme (Ferziger and Peric 2002) in respect to the time discretisation and the implicit Crank–Nicolson method (Ferziger and Peric 2002) is employed for the time discretisation of the diffusion terms. The solutions in space are derived via a finite-volume procedure, incorporating a staggered variable arrangement for the vector and scalar fields. A second-order, central-differencing scheme (Ferziger and Peric 2002) is applied to the advection and diffusion terms.

The location of the interface between the liquid and vapour phases is estimated from the solution of the following transport equation of the colour function:

$$\begin{aligned} \frac{\partial \phi }{\partial t} + \nabla \cdot \left( \phi \mathbf {U}\right) =-\varGamma \frac{1}{\rho _L}. \end{aligned}$$

(10)

For the derivation and discretisation of Eq. 10 the reader is referred to Sato and Ničeno (2013). In the second step of the interface-capturing method, an interface sharpening procedure is invoked in order to prevent the smearing of the colour function, and keep the interface thickness sharp (Olsson and Kreiss 2005), as is common in other works (Sato and Ničeno 2012).

The interfacial mass transfer rate is computed directly from the heat flux at the liquid/vapour interface as follows (Sato and Ničeno 2013):

$$\begin{aligned} \dot{m}=\frac{\left( q_L + q_V \right) }{H_{PC}}, \end{aligned}$$

(11)

in which

$$\begin{aligned} q_L= & {} \lambda _L\left( \nabla T_L\right) \mathbf {n}, \end{aligned}$$

(12)

$$\begin{aligned} q_V= & {} -\lambda _V\left( \nabla T_V\right) \mathbf {n}. \end{aligned}$$

(13)

In this set of Eqs. 11–13, $\dot{m}$ presents the interfacial mass transfer rate per unit area (kg/m$^2$ s); $\mathbf {n}$ is the normal vector to the liquid/vapour interface pointing from the vapour phase to the liquid phase; $H_{PC}$ is the latent heat of phase change from water to vapour (J/kg); $q_V$ (W/m$^2$) and $q_L$ (W/m$^2$) denote the respective heat fluxes to the interface from the vapour and liquid sides, while $\nabla T_L$ and $\nabla T_V$ are the temperature gradients on the liquid and vapour phases of the interface, respectively.

The mass transfer rate per unit volume $\varGamma$ (kg/m$^3$ s) may be expressed as follows:

$$\begin{aligned} \varGamma = \dot{m}A=\dot{m}\frac{S_{i}}{V_{c}}, \end{aligned}$$

(14)

where the interfacial area density A is defined as the ratio of the area of the liquid/vapour interface in the given cell, $S_{i}$ (m$^2$), to the total cell volume $V_{c}$ (m$^3$). The liquid/vapour interface area, A, is here calculated using a high-resolution 3D surface construction algorithm (Lorensen and Cline 1987). From Eq. 14, it is evident that $\varGamma$ is exactly zero for all non-interface cells.

For the discretisation of the enthalpy conservation equation in time, forward and backward Euler techniques are employed for the advection and diffusion terms, respectively (Ferziger and Peric 2002). Details of the space discretisation, as well as of the verification and validation of the method employed can be found in Sato and Ničeno (2013).

2.2 LES Method

Since DNS of two-phase flows at even moderate Reynolds numbers remains very CPU intensive, a highly-resolved LES method (cf. profiles of the modelled viscosity and thermal conductivity in App. (A) is used instead in the present work to investigate condensing, two-phase flow at a moderate friction Reynolds number.

In the LES approach, any flow variable $\psi$ is decomposed into a resolved component $\overline{\psi }$ and an unresolved component $\psi '$. This decomposition is the result of the application of a filtering operation, as originally proposed by Leonard (1975). One important quantity in the filtering algorithm is the filter width, which is usually denoted by $\varDelta$. Scales larger than $\varDelta$ are retained in the flow field, and are referred to as resolved or super-grid scales (GS), while the effects of scales smaller than $\varDelta$ need to be modelled. The unresolved part is called the sub-grid-scale (SGS) model. In the current study, $\varDelta$ is specified by $\varDelta =(\varDelta x \varDelta y \varDelta z)^{1/3}$.

Application of the filtering operation to the Navier–Stokes and energy equations, Eqs. 2–4, enables to derive the equations for the resolved part of the flow (Lesieur et al. 2005). In the current work, the Wall-Adapting Local Eddy-Viscosity (WALE) model, proposed by Nicoud and Ducros (1999), is used to calculate the SGS dynamic viscosity $\mu _t$. This model produces the correct asymptotic behaviour of the eddy viscosity in the near-wall region without the use of a damping function. The WALE model constant $C_W$ is set to 0.45 Nicoud and Ducros (1999). The modelling of SGS thermal conductivity $\lambda _t$ is simply via the concept of a turbulent Prandtl number, $Pr_t$:

$$\begin{aligned} Pr_t=\frac{\mu _t c_p}{\lambda _t}. \end{aligned}$$

(15)

Several DNS studies indicate that the turbulent Prandtl number ranges between 0.9 and 1.0 (Kim et al. 1987; Lyons et al. 1991; Kasagi et al. 1992; Nicoud 1998). Based on these findings, $Pr_t$ is set to the constant value of 0.9 in the current work. The choice of the simple SGS model for the computation of the unresolved heat flux is motivated by an expectation that, for a smaller modelled contribution, there is also a smaller potential for error. Hence, satisfactory results can be obtained using even relatively simple SGS models.

3 Computational Set-Up

3.1 Physical Problem

A schematic of the considered flow situation is given in Fig. 1. Two incompressible fluids, flowing in opposite directions, are separated by a (deformable) phase interface, and confined laterally by two non-slip walls: subcooled liquid is situated beneath the interface and saturated vapour lies above, providing a gravitationally stable configuration, cf. Fig. 3. The respective phases are driven in counter-current flow by imposed (constant) mean horizontal pressure gradient and the gravitational force. The vapour and liquid phases are here saturated steam and subcooled water, respectively. The physical properties of the working fluids are evaluated at a system pressure of 5 MPa, and an associated saturation temperature $T_{Sat}$ of 537.09 K, as listed in Table 1.

Fig. 1

Schematic of the physical situation represented

Table 1 Properties of working fluids at 5 MPa

To ensure that steady-state conditions prevail, the mean interfacial shear stress on the liquid side needs to converge to that on the gas side: i.e.

$$\begin{aligned} \tau _{i,L}=\tau _{i,V}, \end{aligned}$$

(16)

which necessitates that the pressure at the interface is equal for both phases:

$$\begin{aligned} P_{i,V}= & {} P_{i,L}, \end{aligned}$$

(17)

$$\begin{aligned} P_{i,V}= & {} P_{x} + \frac{\partial P_V}{\partial x}x, \end{aligned}$$

(18)

$$\begin{aligned} P_{i,L}= & {} P_{x} + \frac{\partial P_L}{\partial x}x. \end{aligned}$$

(19)

From Eqs. 18 and 19 it is evident that the condition (17) is valid only if the horizontal pressure gradients $\partial P_V/\partial x$ and $\partial P_L/\partial x$ are equal. This means that one phase flowing in the positive x direction (in the present case it is the liquid phase) needs to overcome this positive pressure gradient. This is achieved by tilting the channel at an angle $\alpha$. In order to determine the angle $\alpha$, the shear stress continuity condition, Eq. 16, is used:

$$\begin{aligned} \mu _V \left( \frac{1}{\mu _V} \frac{\partial P}{\partial x}-\frac{1}{\nu _V}g \sin \alpha \right) h=-\mu _L \left( \frac{1}{\mu _L} \frac{\partial P}{\partial x}-\frac{1}{\nu _L}g\sin \alpha \right) h, \end{aligned}$$

(20)

where $\nu _V$ and $\nu _L$ are the vapour and liquid kinematic viscosities, respectively. From Eq. 20, it can be derived that

$$\begin{aligned} \sin \alpha = 2 \left| \frac{\partial P}{\partial x} \right| \frac{1}{g(\rho _L + \rho _V)}. \end{aligned}$$

(21)

Since the vapour density is considerably smaller than the liquid density, the angle $\alpha$ is very small ($\sim O\left( 10^{-4\circ }\right)$) for the given friction Reynolds numbers.

To study the effects of condensation on turbulence, and vice versa, three degrees of water subcooling are investigated in this work: $\varDelta T_{sub1}=0$ K, $\varDelta T_{sub2}=10$ K and $\varDelta T_{sub3}=40$ K. In the first case, phase change is suppressed by imposing saturation conditions for both phases. One is then able to examine the dynamics of the flow in the absence of phase-change complications; the base case. To investigate the effects of surface shear, two Reynolds numbers, as defined in Eq. 25, are selected: $Re_*=178$ and $Re_*=590$. The condensing, two-phase flow situation at the lower Reynolds number is studied by means of DNS, since this is computationally achievable, whereas the higher-Re case is investigated using highly-resolved LES. Since the variations in the water properties (density, thermal conductivity, specific heat capacity and surface tension) at the given levels of subcooling are small, all molecular properties are taken at saturation temperature. The density ratio $R=\sqrt{\rho _L/\rho _V}$ is set equal to 5.54 in the simulations. This quantity, not only representative of reality in many steam/water situations of physical relevance, can also be considered as an indicator for strong dynamic coupling between the phases, and thereby the magnitude of interfacial momentum transfer (Lombardi et al. 1996).

The simulated two-phase flow situation is often characterised in terms of the following dimensionless parameters (Lakehal et al. 2008):

$$\begin{aligned}&We = \frac{\rho _L h u_{*L}^2}{\sigma }=\frac{\rho _V h u_{*L}^2}{\sigma }, \end{aligned}$$

(22)

$$\begin{aligned}&Ja_L=\frac{c_{p,L} \varDelta T_{sub}}{H_{PC}}, \end{aligned}$$

(23)

$$\begin{aligned}&Fr = \frac{ u_{*L} \sqrt{\rho _L}}{\sqrt{gh(\rho _L - \rho _V)}}=\frac{u_{*V} \sqrt{\rho _V}}{\sqrt{gh(\rho _L - \rho _V)}}, \end{aligned}$$

(24)

$$\begin{aligned}&Re_* = \frac{u_{*L} h}{\nu _L} = \frac{u_{*V} h}{\nu _V}, \end{aligned}$$

(25)

$$\begin{aligned}&Pr_L=\frac{\nu _L}{a_L}=\frac{\mu _L c_{p,L}}{\lambda _L}. \end{aligned}$$

(26)

The Weber number We describes the ratio of inertia forces to that of surface tension $\sigma$. The Jakob number Ja represents the ratio of the sensible heat to latent heat in heat transfer processes. The Froude number Fr is the ratio of inertial to gravitational forces, and the shear Reynolds number $Re_*$ is defined, as always, in terms of the imposed shear velocity $u_*=\sqrt{h/\rho \ \partial P/\partial x}$, being equal for the two phases. In Eq. 26, $a_L$ denotes the thermal diffusivity. Table 2 lists explicitly the values adopted in the present work. For the sake of clarity, the simulation name is identified in terms of the shear Reynolds number and degree of water subcooling: for example, the simulation for $Re_*=178$ and $\varDelta T_{sub1}=10$ K is labelled as 178_10K.

Table 2 Geometry and flow parameters

To present and discuss the simulation results (cf. Sect. 4) a suitable coordinate system, and appropriate flow regions, both need to be defined. As shown in Fig. 2, the z-axis represents the interface/wall normal (or vertical) direction. Each subdomain, of depth 2h, occupied by either the vapour or liquid in a stable stratification, is divided into two parts, an interface region and a wall region: for their definition, the dimensionless coordinate z/h is used. The wall region is that between the bottom/top wall (z/h = 2) and the centre lines (z/h = 1) within the respective vapour/liquid regions, while the interface region is that between the interface itself (i.e. the red line, z/h = 0) and the two centre lines. Note that sometimes z/h = − 2 and z/h = 2 are used in the literature to indicate the bottom and top walls, respectively. This notation has been avoided here: the vapour velocity is more than five times greater than that of the liquid (as a consequence of the density difference, given that the same pressure gradient is being applied to drive the flows). Hence, it is advantageous to consider the mean profiles of the variables, principally the mean velocity and the velocity fluctuations, on each side of the interface separately.

Fig. 2

Definition of the coordinate system, wall and interface regions used in the present work

Results from a statistical analysis are presented here in normalised form: the normalisation defined in terms of the shear velocities according to $u_{\tau ,i}=\sqrt{\tau _i/\rho }$ or $u_{\tau ,w}=\sqrt{\tau _w/\rho }$, as appropriate. Note that, for both the interface and wall regions, the dissipation temperatures are defined as follows: $T_{\tau ,i}=q_i/(\rho _L c_{p,L} u_{\tau ,i})$ or $T_{\tau ,w}=q_w/(\rho _L c_{p,L} u_{\tau ,w})$; the vertical coordinate has been normalised according to $\nu /u_{\tau }$. In the wall region, the mean profiles are plotted as functions of the dimensionless distance from the wall $z_{\tau ,w}^+$, for which $z_{\tau ,w}^+=0$ indicates either the bottom or top wall, as appropriate. For the presentation of the results in the interface region, a similar procedure is followed: the mean profiles are displayed as functions of the dimensionless distance from the interface $z_{\tau ,i}^+$, where $z_{\tau ,w}^+=0$ designates the interface itself. Note that $z_{\tau ,i}^+$ is not exactly zero at the interface, since $\varDelta z_{i}^+ \ne 0$ there, but it is nonetheless sufficiently small ($z_{\tau ,i}^+<1$) to enable the viscous sublayers surrounding the interface to be properly resolved.

3.2 Computational Domain and Grid

In the DNS study, the domain lengths for the computation are given by $L_x=3.82\pi h$, $L_y=1.33\pi h$, $L_z=4h$. In shear-based units ($L^+=Lu_{*L}/\nu _L$), the product $L_x^+\times L_y^+\times L_z^+$ is $2136\times 748 \times 712$. The calculations are performed on a grid with $192\times 128 \times 512$ cells, resulting in a resolution of $\varDelta x^+=11.12$, $\varDelta y^+=5.84$ and $\varDelta z^+=0.56$ to 1.81. The choice of this grid resolution is based on an analysis of our own single-phase DNS of turbulent channel flow at $Re_*=178$ and those from the literature (Kim et al. 1987; Komori et al. 1993).

The time step $\varDelta t$, which is limited both by the Courant–Friedrich–Levy (CFL) condition and the stability condition for the surface tension Brackbill et al. (1992), is given by:

$$\begin{aligned} \varDelta t=min \left( C_{CFL} \frac{\varDelta x }{|\mathbf {U}|_{max}},C_{\sigma } \left( \frac{\rho _V x^3}{2 \pi \sigma }\right) ^{1/2} \right) , \end{aligned}$$

(27)

where $C_{CFL}$ and $C_{\sigma }$ are dimensionless constants, representing safety factors. The first term describes the limitation imposed by the CFL condition, and $C_{CFL}=0.25$ is used for all the simulations reported in the paper. The second term represents the one due to the surface tension effect. In the original CSF model Brackbill et al. (1992), the term is defined as: $\left( \frac{0.5(\rho _V+\rho _L) x^3}{2 \pi \sigma }\right) ^{1/2}$. However, according to our experience, numerical instability ensues when this definition is used. Thus, in the present work, we use $\rho _V$ instead of in Eq. 27. Since $\rho _V \ll \rho _L$, this resolves the issue of numerical instability. Nonetheless $\varDelta t$ limited by the factor $\left( \frac{\rho _V x^3}{2 \pi \sigma }\right) ^{1/2}$ imposes an unaffordable restriction on computer resources. Consequently, we increase $\varDelta t$ by multiplying a factor $C_{\sigma }=10$, which is subsequently used for all the computations reported in this paper. In the simulation without condensation, a time step of $\varDelta t=0.0052\nu _L/u_{*L}^2$ in dimensionless time wall units is employed, while, in the case with condensation, $\varDelta t=0.0024\nu _L/u_{*L}^2$. The corresponding maximum CFL limits are 0.12 and 0.06, respectively. For comparison, a typical time step in a DNS simulation of single-phase, turbulent channel flow at a similar friction Reynolds number would be about $0.1\nu /u_{*}^2$, resulting in a CFL limit less than or equal to 0.3 (Bergant and Tiselj 2007; Komori et al. 1993). This shows that two-phase simulations involving phase change at moveable interfaces are intrinsically less stable, and significantly more time consuming than their single-phase counter parts. Here, transient data are collected over a dimensionless time interval of approximately $6000\nu _L/u_{*L}^2$. For the liquid phase, it corresponds to a number of flowthroughs of 45. A sampling frequency of $2.44\nu _L/u_{*L}^2$ is used for all the simulations. Each calculation has been performed on a cluster computer system using 192 cores (node configuration: 2x Intel(R) Xeon(R) CPU E5-2690 (8 cores), 2.90 GHz, 2 GB RAM per core). Each simulation took about 5 months to reach convergence, which illustrates the computational overhead in applying DNS to problems of this type.

In the LES investigation of stratified condensing two-phase flow at the moderate Reynolds number of $Re_*=590$, a domain size is $3895\times 2478 \times 2360$ (in wall units) was employed. This particular choice of grid resolution is based on the aim to resolve as much of the two-phase flow field as possible. Always, the grid resolution should be as fine as possible, given the constraints of the available computational resources. Based on these considerations, the computational domain is discretised using 256 cells in the streamwise and spanwise directions, and 768 cells in the vertical direction. This choice results in a grid spacing of $\varDelta x^+=15.2$, $\varDelta y^+=12.9$ and $\varDelta z^+=0.6-7.9$. The grid is stretched in the wall/interface-normal direction by a factor less than 1.05, as recommended by Fröhlich (2006). The total number of cells is then about 50.3 million. Such a grid resolution avoids the use of wall functions. The time step is chosen as $\varDelta t=0.0043\nu _L/u_{*L}^2$, to maintain the numerical stability of the simulation. The transient data, collected with a sample frequency of $7\nu _L/u_{*L}^2$, are averaged over a time interval of $5800\nu _L/u_{*L}^2$ in wall-time units. For the liquid phase, it is equivalent to a number of flowthroughs of 24. The calculation has been performed using 512 processors (node configuration: 2x Intel Haswell Xeon CPU E5-2680 (12 cores), 2.5 GHz, 2.5 GB RAM per core). The simulation took approximately 4 months to complete.

3.3 Boundary and Initial Conditions

In the DNS study of Fulgosi (2004), who first simulated condensing, counter-current, stratified water/steam flow in a horizontal channel, periodic boundary conditions were employed in the streamwise and spanwise directions. At the outer boundaries, a free-slip condition was imposed for each velocity component (i.e. $\partial U / \partial z = \partial V / \partial z = \partial W / \partial z = 0$) in the case without phase change. In the presence of phase change phenomena, the problem becomes more complicated. The total volume of each phase changes as a consequence of condensation. In order to maintain steady-state conditions, and to satisfy mass balance, the condensing mass flux was injected at the top of domain and withdrawn at the bottom of the domain (Lakehal et al. 2008). To this purpose, the boundary condition for the vertical velocity W was changed to one of the Dirichlet type, according to:

$$\begin{aligned} W_V= & {} -\frac{\sum \nolimits _{i,j=1}^{N_x N_y}\dot{m}_{ij}V_{c,ij}}{\rho _V}, \end{aligned}$$

(28)

$$\begin{aligned} W_L= & {} -\frac{\sum \nolimits _{i,j=1}^{N_x N_y}\dot{m}_{ij}V_{c,ij}}{\rho _L}, \end{aligned}$$

(29)

in which $N_x$ and $N_y$ are the number of cells in the streamwise and spanwise directions, respectively.

In the present study, periodic boundary conditions are used in the streamwise and spanwise directions. In the simulation without condensation, contrary to the DNS study of Fulgosi (2004), a non-slip (wall) condition is imposed for the velocity components at the outer boundaries, i.e.

$$\begin{aligned} U_{L,V}=V_{L,V}=W_{L,V}=0 \end{aligned}$$

(30)

In the presence of phase change, the boundary conditions for the vertical velocity at the top and bottom boundaries are set according to Eqs. 28 and 29, respectively. The velocities at these boundaries ($W_V \sim O(10^{-5} m/s), W_L \sim O(10^{-6} m/s)$) are updated at each time step. Since the values are very close to zero, the boundary conditions given in Eq. 28 and Eq. 29 have the same effect on the flow as the imposition of a rigid wall boundary condition. Therefore, the periodic boundary conditions imposed in the streamwise and spanwise directions are not influenced by the re-injection and extraction of the condensing mass flux.

In order to simulate the phase change effect due to condensation, a thermal equilibrium condition is imposed at the interface: i.e. the vapour and liquid phases are assumed to be at saturation temperature (cf. Fig. 3):

$$\begin{aligned} T_i = T_{Sat,V}=T_{Sat,L}=T_{Sat}. \end{aligned}$$

(31)

Water subcooling is maintained at the constant temperature of $T_{Sat}-\varDelta T_{sub}$ at the lower boundary of the domain, while the saturation temperature is prescribed at the upper boundary.

Fig. 3

Sketch of boundary conditions

Appropriate choice of the initial conditions can significantly reduce the computational effort. For a simulation involving no phase change, a method for generating a Gaussian random field is used to initialise the spanwise and vertical velocity components (Kraichnan 1970; Li et al. 1994). A velocity field calculated in such a way resembles pseudo-isotropic turbulence, which aids convergence. The streamwise velocity component is stipulated using Reichardt’s wall law (Reichardt 1951):

$$\begin{aligned} \frac{U}{u_*}=2.5 \ln \left( 1+0.4z^+\right) +7.8\left( 1-e^{-z^+/11}-\frac{z^+}{11}e^{-0.33z^+}\right) , \end{aligned}$$

(32)

where $z^+=zu_*/\nu$ is the dimensionless distance from the wall or interface, as calculated from the mean value of the imposed shear velocity. The results for the isothermal case (zero subcooling) are used as initial conditions for the calculations involving phase change. Additionally, the temperature field is initialised assuming a linear temperature gradient prevails between the subcooled wall at the base of the channel and the saturated interface above.

For continuity purposes, in the LES study, the boundary conditions are the same as those used in the DNS simulations. The velocity and temperature fields are initialised in the same manner, as described above. A critical evaluation of the two approaches will then be made possible.

4 Results and Discussion

For a detail description of the verification and validation of the DNS and LES methodologies used here, the reader is referred to Appendix B and Appendix C, respectively.

4.1 Determination of the Mean Quantities

As set up, the two-phase, turbulent flow situation that has artificially been constructed here is, statistically, a stationary process, meaning that all variables do not change as functions of time. Fully developed, turbulent conditions are achieved following an initial transient period of approximately $(2000-2200)\nu _L/u_{*,L}^2$ time units. The statistically stationary-state is confirmed from examination of the viscous and Reynolds stress profiles.

In the presence of interfacial waves (wave slope about 0.02), case $Re_*=590$, a control volume (or cell) in the vicinity of the interface might be occupied by liquid or vapour, or a combination of both. Evidently then, any averaging procedure involved in the numerical algorithm has to take into account the presence of the waves associated with the liquid and vapour phase motions. To define an appropriate mean value, a phase-weighted average can be used. For any given variable F (e.g. velocity or temperature) of the kth phase, the phase-weighted average can be expressed as:

$$\begin{aligned} \langle F_k \rangle = \langle \langle \varPhi _k F_k \rangle \rangle / \langle \langle \varPhi _k \rangle \rangle , \end{aligned}$$

(33)

where the double brackets $\langle \langle \cdots \rangle \rangle$ represent time averaging (Ishii and Hibiki 2011). In the present study, the colour function $c_k$ has been employed as a weighting function, represented by $\varPhi _k$ in Eq. 33.

The velocity profiles in the interface/wall-normal direction were averaged over time, once a statistically steady-state condition had been achieved:

$$\begin{aligned} \langle U \rangle (x,y,z) = \langle \langle c_L U \rangle \rangle / \langle \langle c_L \rangle \rangle . \end{aligned}$$

(34)

In order to reduce statistical uncertainties, the velocities were additionally averaged across the horizontal planes of the configuration (i.e. in the homogeneous directions). Root-mean-square fluctuations of the streamwise velocity component were computed for both phases according to the following prescription:

$$\begin{aligned}&u'^2 (x,y,z) = \frac{1}{T} \sum \limits _{t=0}^{T} \left( U(x,y,z,t) - \langle U \rangle (x,y,z)\right) ^2, \end{aligned}$$

(35)

$$\begin{aligned}&u_{rms} (z) = \sqrt{u'^2 (z)}. \end{aligned}$$

(36)

The fluctuations of the spanwise and vertical velocity components were calculated in the same manner.

According to the literature, related studies of heat transport in single-phase flow using DNS, e.g. Bergant and Tiselj (2007), the dimensional temperature T is usually transformed to dimensionless form using

$$\begin{aligned} \theta ^+ = (T_{ref} - T) / T_{\tau }, \end{aligned}$$

(37)

where $T_{ref}$ is a reference temperature (usually the wall temperature) and $T_{\tau }$ is the dissipation temperature. In the present study, the interface temperature $T_i$ is used to represent $T_{ref}$. The mean temperature profiles are computed in the same way as those for the velocities. In the case of $Re_*=590$, all the profiles presented in the figures below relate to the resolved or filtered quantities, as appropriate. For the sake of clarity, all the filtered variables are written without an overbar.

4.2 Velocity Field

Figure 4a, c display the liquid mean velocity profiles in the interface region. If the interfacial velocity is subtracted out, all the profiles for $Re_*= 178$ collapse onto a single curve in the viscous, buffer and logarithmic sublayers, as revealed by Fig. 4b. The velocity profile for $Re_*= 590$ differs somewhat from those for $Re_*=178$, with the maximum offset being observed in the logarithmic sublayer $z_{\tau ,i}^+>30$, cf. Fig. 4d.

Fig. 4

Liquid interface region: dimensionless liquid mean velocity profiles for $Re_*=178$ (a, b), dimensionless liquid mean velocity profiles for $Re_*=590$ (c, d), the diagnostic quantity $\gamma$ for a velocity log law description (e), and dimensionless mean velocity profile for $Re_*=590$ showing the law of the wake (f)

In addition, the variation of the mean velocity profiles for increased water subcooling and Reynolds number can be distinguished in terms of the more sensitive parameter $\gamma =z_\tau ^+\left( \partial U / \partial z \right) ^+$, with a value of $1/\kappa$ ($\kappa$ is the von Karman constant) in the log-law range. From $z_{\tau ,i}^+=50$ to $z_{\tau ,i}^+=100$, $\gamma$ increases slightly for the $Re_*=178$ case, the corresponding values of $1/\kappa$ varying from 2.43 to 2.93. The $\kappa$ values lie between 0.41 and 0.34 in the interface region, cf. Fig. 4e. For the case $Re_*=590$, lower values of $\gamma$ are predicted, varying from 1.75 to 2.95 in the range $30<z_{\tau ,i}^+<250$. Closer analysis of the diagnostic quantity $\gamma$ clearly shows that, in the simulations for $Re_*=178$, the log sublayer is very narrow. It is also important to note that in the case $Re_*=590$, the liquid mean velocity deviates significantly from the log law for $z_{\tau ,i}^+>250$.

According to Coles (1956), this region is known as the “defect layer”, characterised by a universal function $f \left( z/h \right)$ and a flow-dependent wake parameter $\varPi$. The former usually takes the form of $f\left( z/h\right) =2sin^2\left( \frac{\pi }{2}\frac{z}{h} \right)$. In order to determinate $\varPi$, the mean velocity profile is formulated in the form of the velocity defect law. In the present circumstances, a value of $\varPi =0.45$ may be estimated within the interface region of the liquid phase. This is an interesting result since the wake strength should be very small (i.e. $\varPi \approx 0$) for channel/pipe flows (Pope 2000). As a result, the mean velocity can be effectively represented as the sum of the log and wake laws, as illustrated in Fig. 4f, which depicts collectively the actual log law (solid black line), the wake contribution $\varPi f\left( z/h\right) /\kappa$ with $\varPi =0.45$ and $1/\kappa =1.85$ (black dashed line), as well as the sum of the log-law and the wake contributions. Obviously, the sum agrees very well with the liquid mean velocity profile calculated here for the 590_10K case.

The liquid mean velocity profiles in the near-wall region are presented in Fig. 5a, b. The DNS data for turbulent, single-phase channel flow at $Re_*=180$ and $Re_*=590$ obtained by Kim et al. (1987) and Moser et al. (1999) respectively, are also included for comparison. Obviously, the calculated profiles coincide with the reference data in the viscous and buffer sublayers, but for $z_{\tau , w}^+>20$ they begin to differ. This may be seen more clearly in Fig. 5c, which shows the diagnostic quantity $\gamma$ for a log law formulation. For $Re_*=178$, the mean profiles of $\gamma$ reach a plateau between $z_{\tau ,i}^+=50$ and $z_{\tau , i}^+=100$, but the duration of the plateau decreases as the water subcooling is increased. For $50<z_{\tau ,i}^+<100$, $\gamma$ is approximately equal to 2.6 (for $\kappa = 0.38$). In contrast to the interface region, however, the defect layer is less pronounced in the wall region, for which a value of $\varPi =0.15$ is obtained. Since the wake strength in the wall region is significantly smaller than that derived for the interface region, the wake contribution to the magnitude of the liquid mean velocity is considerably smaller there. This is clearly illustrated in Fig. 5d.

Fig. 5

Liquid wall region: dimensionless liquid mean velocity profiles for $Re_*=178$ (a); dimensionless liquid mean velocity profiles for $Re_*=590$ (b); the diagnostic quantity $\gamma$ based on a velocity log law description (c); and dimensionless mean velocity profile for $Re_*=590$ showing the law of the wake (d)

Turning now to the vapour side, Fig. 6a, c displays the dimensionless vapour mean velocity profiles in the interface region. It can be clearly seen that the mean vapour interfacial velocity increases with increase in the degree of water subcooling (cf. “178_0K”, “178_10K” and “178_40K”) and for increased Reynolds number $Re_*$ (i.e “178_10K” and “590_10K”). If the averaged interfacial velocity is subtracted out, as depicted in Fig. 6b, all three profiles for $Re_*=178$ coincide for $z_{\tau ,i}^+<5$, but start to deviate noticeably from each other for $z_{\tau ,i}^+>50$. In Fig. 6e, the profiles of $\gamma =z_\tau ^+\left( \partial U / \partial z \right) ^+$ are shown for the vapour interface region. For the “178_40K” case there is no distinctive plateau for $z_{\tau ,i}^+>50$, which indicates that the log layer is considerably narrower for higher degrees of water subcooling. The wake strength, computed for the vapour phase in the interface region, is given by $\varPi =0.15$. Adopting this value, very good agreement between the “590_10K” profile and the sum of the log and wake laws is achieved, as revealed in Fig. 6f.

Fig. 6

Vapour interface region: dimensionless vapour mean velocity profiles for $Re_*=178$ (a, b), dimensionless vapour mean velocity profiles for $Re_*=590$ (c, d), the diagnostic quantity $\gamma$ for a velocity log-law description (e), and dimensionless mean velocity profile for $Re_*=590$ showing the law of the wake (f)

In the region very close to the upper wall, the vapour mean velocity profile agrees well with the linear form of the wall law (not shown here), as well as with the data of Kim et al. (1987) and Moser et al. (1999); cf. Fig. 7a and b, respectively. The diagnostic quantity $\gamma$ is drawn in Fig. 7c. For the lower Reynolds number case, $Re_*=178$, $\gamma \approx 2.65$ for $50<z_{\tau ,w}^+<100$. With an increase of water subcooling from 10K to 40K, however, $\gamma$ increases to 2.89. For the higher $Re_*=590$ case, the mean profile of $\gamma$ agrees well with the “Moser_590” profile (Moser et al. 1999): $\gamma \approx 2.3$. It is interesting to note that a larger wake region can be discerned in the vapour wall region, characterised by a wake strength of $\varPi =0.23$. The sum of the log and wake contributions is given in Fig. 7d, which is in good agreement with the calculated vapour mean velocity.

Fig. 7

Vapour wall region: dimensionless vapour mean velocity profiles for $Re_*=178$ (a), dimensionless vapour mean velocity profiles for $Re_*=590$ (b), the diagnostic quantity $\gamma$ for a velocity log-law description (b), and dimensionless mean velocity profile for $Re_{*}=590$ showing the law of the wake (c)

Tables 3 and 4 summarise the log laws obtained from the DNS and LES data for the liquid and vapour phases, respectively. Analysis of the velocity field of condensing, two-phase flow at $Re_*=178$ and $Re_*=590$ clearly shows that the universal wall law is valid for the upper and lower wall regions. At the interface region itself, the situation is more complex. In the case $Re_*=590$, for which the interface is wavy, the mean velocity profile does not follow the linear form of the wall law in either the liquid or the vapour phases. In the absence of interfacial waves (i.e. the cases for which $Re_*=178$), this tendency is observed only on the vapour side; on the liquid side, the law remains valid.

Table 3 Log laws obtained for the liquid phase

Table 4 Log laws obtained for the vapour phase

The RMS velocity fluctuations in the interface region on the liquid side are plotted in Fig. 8. At the interface ($z_{\tau ,i}^+=0$), the $u_{rms,i}^+$ values decrease to 2.38, 2.27 and 1.96 for the “178_0K”, “178_10K” and “178_40K” cases, respectively; cf. the magnified view of the $u_{rms,i}^+$ profiles given in Fig. 8a. The same behaviour is apparent for the spanwise velocity fluctuations. These findings reveal that the damping effect of the interface on the turbulence levels is enhanced in the presence of phase change phenomena. For $Re_*=590$, significant differences can be discerned in comparison with the $Re_*=178$ case. Namely, the fluctuations of the streamwise velocity increase significantly in the vicinity of the interface with increasing $Re_*$. It is also evident that, in contrast to the lower $Re_*=178$ case, the fluctuations of the vertical velocity component $w_{rms}$ are enhanced as the interface is approached.

Fig. 8

Dimensionless RMS velocity fluctuations in the interface region of the liquid phase, both in shear-based coordinates (a, c, e) and in global coordinates (b, d, f)

Figure 9 shows the vapour RMS velocity fluctuations in the interface region. For the lower $Re_*$, it can be seen that both $u_{rms}$ and $v_{rms}$ have considerably smaller values at the interface compared to those on the liquid side. Moreover, $u_{rms,i}^+$ and $v_{rms,i}^+$ decrease with increase in the degree of water subcooling. The variations in the maxima of $v_{rms}$ and $w_{rms}$ are themselves sensitive to the degree of water subcooling. Comparison of the “178_10K” and “590_10K” cases reveals that with increase of Reynolds number the damping effect on the turbulence due to the presence of the interface is significantly weakened. Thus, it appears that the statement that the vapour phase “sees” a moving wall is only of limited validity for the case $Re_*=590$.

Fig. 9

Dimensionless RMS velocity fluctuations in the interface region of the vapour phase, both in shear-based coordinates (a, c, e) and in global coordinates (b, d, f)

In the wall regions for both phases, Figs. 10 and 11, the profiles are compared with the single-phase DNS data of Moser et al. (1999): good agreement is observed for the “590_10K” case. Inspection of the RMS profiles reveals that the ratio of the maximum value of $u_{rms}$ in the interface region to that in the wall region remains almost constant as the Reynolds number is increased: namely, 1.10 vs. 1.14 in the cases $Re_*=178$ and $Re_*=590$, respectively.

Fig. 10

Dimensionless root-mean-square velocity fluctuation in the wall region of the liquid phase: in shear-based coordinates (a, c, e), and in global coordinates (b, d, f)

Fig. 11

Dimensionless RMS velocity fluctuations in the wall region of the vapour phase, both in shear-based coordinates (a, c, e) and global coordinates (b, d, f)

The above analysis of the velocity fluctuations reveals that with increase of water subcooling, the $v_{rms}$ and $w_{rms}$ values of the liquid phase are enhanced in the interface region, but reduced in the wall region. On the vapour side, the opposite is true. Note that all the profiles of velocity fluctuations are asymmetric with the respect to the centre line.

The dimensionless liquid Reynolds stress profiles in the interface and wall regions are presented in Fig. 12; the data are normalised in terms of either $\rho u_{\tau ,i}^2$ or $\rho u_{\tau ,w}^2$, as appropriate. In the interface region, Fig. 12a, b, the Reynolds stress increases with increase in water subcooling, while the opposite trend is observed in the wall region: Fig. 12c, d. It can be also clearly discerned that the Reynolds stress significantly increases with increase in $Re_*$. It should be pointed out that the viscous (not shown here) and turbulent stresses equalize at $z_{\tau ,i}^+=8$ in the interface region for both Reynolds numbers. In the wall region, the data of Moser et al. (1999) are included for comparison. It is worth noting that the location of the stress equilibrium point is at $z_{\tau ,w}^+=12$ for the “178_0K”, “178_10K” and “590_10K” cases, and at $z_{\tau ,w}^+=13$ for the “178_40K” case. In each simulation, the location at which the stresses are balanced agrees with the location of the maximum in the production rate. These findings are in line with the theory of single-phase, wall-bounded turbulent flow (Pope 2000).

Fig. 12

Profiles of Reynolds stresses in dimensionless form on the liquid side, both for interface (a, b) and wall (c, d) regions

The vapour Reynolds stress profiles in the interface and wall regions are shown in Fig. 13. The dependency of $\langle u'w' \rangle ^+$ on the degree of subcooling in the interface and wall regions is opposite to those seen on the liquid side. In the interface region, the molecular viscous shear stress (not shown here) equalises with that of the turbulent one at $z_{\tau ,i}^+=13.5$, $z_{\tau ,i}^+=14.6$, $z_{\tau ,i}^+=16.6$ and $z_{\tau ,i}^+=15.3$ for the “178_0K”, “178_10K”, “178_40K” and “590_10K” cases, respectively. In the near-wall region, the cross-over point is located at $z_{\tau ,w}^+=12$ for the “178_0K”, “178_10K” and “590_10K” cases, and at $z_{\tau ,w}^+=11$ for the “178_40K” case. At the higher Reynolds number, it is to be noted that $\langle u'w' \rangle _L^+$ and $\langle u'w' \rangle _V^+$ do not vanish at the interface.

Fig. 13

Profiles of Reynolds stresses in dimensionless form on the vapour side, for both the interface (a, b) and wall (c, d) regions

It is important to note that for $Re_*=590$ the Reynolds stress profiles remain linear (for $z_{\tau ,i}^+>100$, $z_{\tau ,w}^+>100$ and $z/h>0.2$) on both sides of the interface. This observation confirms that the interfacial waves influence the mean velocity field characteristics only in the region very close to the interface. The profiles of the Reynolds stress shear $\langle u'w' \rangle ^+$ emphasise again the significant difference between a rigid wall and a moving interface, which results in the latter case to an asymmetry of the profiles around the centre line. Even in the vapour phase, for which we have noted the interface acts like a rigid wall (at least for the $Re_*=178$ cases), the asymmetry of the stress profiles is still apparent. Far from the wall or interface, the profiles of the velocity fluctuations and Reynolds stresses, presented above in terms of global coordinates (z/h) for the two Reynolds numbers, indicate that in the outer region, $z/h>0.2$, mean quantities scale with outer variables rather than with wall variables.

4.3 Temperature Field

The dimensionless temperatures for the two Reynolds numbers considered here are displayed in Fig. 14a as functions of the dimensionless distance from the interface. It can be clearly seen that for $Re_*=178$ the normalized temperature profiles coincide perfectly. Further analysis reveals that close to the interface a diffusive sublayer exists in which the calculated temperature profiles, for both the “178_10K” and “178_40K” cases, agree with the linear relation $\langle \theta \rangle ^+=z_{\tau ,i}^+Pr$. The temperature profile obtained for the $Re_*=590$ case is similar to that of the velocity in the vicinity of the interface: i.e. the mean temperature is essentially constant between $z_{\tau ,i}^+=0$ and $z_{\tau ,i}^+=10$, and does not follow the linear wall law. In the logarithmic region, $z_{\tau ,i}^+>30$, the present results for $Re_*=178$ can be well approximated using the relation $\langle \theta \rangle ^+=2.9 \ln z_{\tau ,i}^+ +1.7$, where $1/\kappa _\theta =2.9$ ($\kappa _\theta$ is the von Karman constant for the thermal field) is the slope and $B_\theta =1.7$ the shift, respectively. The diagnostic quantity $\gamma _{\theta }=1/\kappa _\theta$ for the temperature field, presented in Fig. 14b, is calculated in a similar way to that for $\gamma$ for the velocity field; i.e. $\gamma _{\theta }=z_{\tau }^+ \langle \partial T/ \partial z \rangle ^+$. It is interesting to note that $\gamma _{\theta }$ does not tend to zero far from the interface, though $\gamma$ does (cf. Fig. 4c). This is due to the fact that the mean temperature gradient does not vanish in the liquid bulk. The given values of the slope $1/\kappa _\theta$ and the shift $B_\theta$ agree quite well with the experimental data of Kader (1981), as well as with the DNS results of Kasagi et al. (1992). It can be clearly seen that, for $Re_*=590$, $\gamma _{\theta }$ is considerably smaller than for the $Re_*=178$ case. This implies that $\kappa _\theta$ increases in the logarithmic region with increase of $Re_*$. In the region $30<z_{\tau ,i}^+<200$, the temperature distribution obeys a logarithmic law of the form $\langle \theta \rangle ^+=2.12 \ln z_{\tau ,i}^+ +3.8$. From comparisons of the temperature profiles for both cases, it can be concluded that for the lower $Re_*=178$ case the logarithmic sublayer ($30<z_{\tau ,i}^+<100$) is considerably narrower, and therefore cannot be distinguished clearly from the wake region. For $Re_*=590$, the wake region, which is now more distinctive, begins approximately at $z_{\tau ,i}^+=200$.

Fig. 14

Dimensionless liquid mean temperature profiles in the interface region (a), the diagnostic quantity $\gamma _{\theta }$ for a temperature log-law description in the interface region (b), dimensionless liquid mean temperature profiles in the wall region (c), and the diagnostic quantity $\gamma _{\theta }$ for a temperature log-law description in the wall region (d)

Figure 14c shows the averaged temperature profiles in the wall region for both Reynolds numbers; here temperatures are plotted against the dimensionless distance from the wall. All the profiles coincide perfectly with each other in the conduction and buffer sublayer regions. Evidently then, a diffusive sublayer also exists very close to the wall. Further, it can be seen that the width of the log sublayer is considerably larger in the case of $Re_*=590$. In the wall region, the temperature profiles have a more significant slope and shift (cf. Table 5). The plots of $\gamma _{\theta }$, given in Fig. 14b, d, confirm that the von Karman constant for the thermal field has a smaller value in the wall region than in the interface region. The well-evidenced logarithmic profile, observed in both the wall and interface regions, provide evidence that the two-phase flow situation considered here at $Re_*=590$ is free of low-Reynolds number effects (this is also confirmed by inspection of the velocity field).

Table 5 Comparison of temperature logarithmic laws obtained for the $Re_*=178$ and $Re_*=590$ cases

The behaviour of the fluctuating temperature field in the liquid phase is shown in Fig. 15. In the interface region, Fig. 15a, the maximum of the temperature fluctuations $\langle \theta ' \rangle ^+$ occurs at $z_{\tau ,i}^+=16$, while in the wall region the maximum is observed at $z_{\tau ,w}^+=19$. The DNS results, obtained by Kawamura et al. (2000) for Poiseuille flow at $Re_*=180$ with constant wall temperature difference for $Pr=0.71$ are also included for comparison. From a qualitative point of view, both the two-phase profiles (“178_10K” and “178_40K”) and the single-phase profile, labelled Kawamura_180, exhibit similar behaviour. In the case of single-phase, wall-bounded flow (Kawamura_180), the near-wall maximum of the temperature fluctuations is $\langle \theta ' \rangle ^+=2.63$, which is lower than for the two-phase condensing cases. It is also interesting to note that the maximum of the temperature fluctuations, located in the bulk of the liquid phase, is about $\langle \theta ' \rangle ^+=3$, both for the single-phase and two-phase situations. For the higher $Re_*=590$ case, the temperature fluctuations do not vanish at the interface, due to the presence of interfacial waves. Near the interface, the maximum $\langle \theta ' \rangle ^+=3.24$ occurs at $z_{\tau ,i}^+=10$, but in the wall region, the temperature fluctuations are much smaller.

Fig. 15

Dimensionless liquid temperature fluctuations in the interface (a) and wall (b) regions

Figure 16 shows the dimensionless, time-averaged turbulent heat flux in the streamwise, $\langle u'\theta ' \rangle ^+$, and interface-normal, $\langle w'\theta ' \rangle ^+$, directions versus the non-dimensional distance from the interface/wall. The data are normalised in terms of $u_{\tau ,i}T_{\tau ,i}$ and $u_{\tau ,w}T_{\tau ,w}$, respectively; the DNS data of Kawamura et al. (2000) are also included for comparison. In single-phase channel flow, the peak value of $\left| \langle u'\theta ' \rangle ^+ \right| =5.93$ is located at $z_{\tau ,i}^+=16$; cf. Fig. 16a. For the lower $Re_*=178$ two-phase flow cases, a maximum value $\left| \langle u'\theta ' \rangle ^+ \right| =6.72$ is predicted in both the interface and wall regions. The near-interface maximum is located at $z_{\tau ,i}^+=13$, while in the wall region the peak value occurs at $z_{\tau ,w}^+=15$.

Fig. 16

Liquid streamwise (a, b) and interface/wall-normal (c, d) turbulent heat flux profiles; and comparison with DNS data of Kawamura et al. (2000)

Turning now to the distribution of the normalised turbulent heat flux in the interface-normal direction, $\langle w'\theta ' \rangle ^+$: very close to the interface/wall ($z_{\tau ,i}^+<30$, $z_{\tau ,w}^+<30$), our DNS data agree quite well with those of Kawamura et al. (2000), also derived using DNS, but for single-phase flow. However, away from the interface/wall, the single-phase profile of the normal turbulent heat flux (labelled Kawamura_180) differs from those calculated for the two-phase condensing cases. It can be seen that in the interface region the normal turbulent heat flux $\langle w'\theta ' \rangle ^+$ increases according to the degree of water subcooling. This is in line with the behaviour of the vertical velocity fluctuations and Reynolds stresses, which also increase with water subcooling (cf. Fig. 12a). In the wall region, the opposite trend is seen: the normal turbulent heat flux decreases with increased water subcooling, since $\langle w'w' \rangle$ decays in this region. For the higher $Re_*$ case, a similar trend than the one observed for the temperature fluctuations can be noticed; and the streamwise and interface-normal turbulent heat fluxes are not zero at the interface. Also, the turbulent heat fluxes are enhanced with increased Reynolds number.

Finally, the mass transfer rates due to condensation are displayed in Fig. 17 as functions of the friction Reynolds number. The condensation rate has been normalized according to $m^+=\dot{m}H_{PC}h/\left( \lambda _L \varDelta T_{sub} \right)$, with $h=0.025m$ (cf. Fig. 17b). The dashed line reflects an assumed linear dependency of the condensation rate on Reynolds number for 10K liquid subcooling, which is seen not to be strictly valid. For the LES validation case, presented in App. C, the mass transfer rate was predicted to an accuracy of 6%. Comparison of the two DNS computations, for which the liquid subcooling has been varied from 10K to 40K (this for the $Re_*=178$ case only), clearly indicates that the condensation rate increases monotonically, though not linearly, with the degree of subcooling (cf. Fig. 17a).

Fig. 17

Absolute (a) and dimensionless (b) condensation rates versus friction Reynolds number

4.4 Consequences for a Reynolds-Averaged Navier–Stokes Modelling of Stratified Two-Phase Flow

The interfacial viscous shear stress, normalised in terms of a suitable reference velocity, is commonly referred to as a drag coefficient. For example, using the bulk velocity in our case, it is defined as:

$$\begin{aligned} C_D = \frac{\tau _\nu }{0.5 \rho U_b^2}. \end{aligned}$$

(38)

The magnitude of the interfacial shear stress is not known a priori in a Reynolds-Averaged Navier–Stokes (RANS) simulation. Hence, for such a computation, the drag coefficient needs to be defined, to take into account the effect of the interfacial shear stress. Using the value obtained in the present DNS and LES simulations, a drag coefficient has been calculated, and is plotted in Fig. 18 as a function of both the liquid subcooling and friction Reynolds number. Evidently, $C_D$ increases non-linearly as $\varDelta T_{sub}$ increases. This trend agrees well with the experimental findings of Kim and Bankoff (1983). Figure 18b shows that for the “590_10K” case $C_D$ is slightly less than that for the “178_10K” simulation. This reduction indicates that a part of the flow energy provided by the constant pressure gradient driving the flow is being used to create and sustain interfacial waves. Thus, that part of the flow energy responsible for the interfacial shear stress is also reduced, which results in a reduced value of $C_D$. Proper modelling of the drag force in the present context should incorporate the fact that $C_D$ is a function of both the Reynolds number and the degree of water subcooling.

Fig. 18

Drag coefficient on the liquid side as functions of liquid subcooling (a) and Reynolds number (b)

The turbulent Prandtl number $Pr_t$ plays an important role in the modelling of turbulent heat transfer within a RANS framework (Hanjalic and Launder 2011). This quantity, defined as

$$\begin{aligned} Pr_t = \frac{\nu _t}{a_t}=\frac{\langle u'w' \rangle }{\langle w'T' \rangle } \frac{\langle \partial T / \partial z \rangle }{\langle \partial U / \partial z \rangle }, \end{aligned}$$

(39)

is plotted in Fig. 19, based on our DNS and LES calculations. Note that, for the higher Reynolds number case, $Pr_t$ is computed in terms of the resolved variables according to the underlying LES strategy. It can be readily seen that both the DNS and LES computations have produced very similar distributions for $Pr_t$ for both the $z_{\tau ,i}^+<50$ and $z_{\tau ,w}^+<50$ regions of the flow, with $Pr_t$ close to unity in both regions. This value is in excellent agreement with the DNS analysis of Kasagi et al. (1992), who obtained a value $Pr_t=1.02$ for single-phase, turbulent channel flow at $Re_*=150$. In our case, sufficiently far from the interface and wall, i.e. $z_{\tau ,i}^+>50$ and $z_{\tau ,w}^+>50$, the same trend is observed: that is, $Pr_t$ decreases progressively, but then increases again for $z_{\tau ,i}^+>100$ for $Re_*=178$. In particular, again for the “178_10K” case, $Pr_t$ falls to 0.75 and 0.55 in the interface and wall regions, respectively. For the “590_10K” case, $Pr_t$ appears to approach a value of 0.75 in the bulk of the liquid phase.

Fig. 19

Distribution of the turbulent Prandtl number in the liquid phase for the interface (a) and wall (b) regions

Generally, to improve RANS modelling of interfacial mass, momentum and energy transfer for large-industrial scale applications, it is recommended to model the principal mean two-phase flow quantities, such as velocity, temperature, turbulent kinetic energy and dissipation, at the interface using a so-called interface function approach, similar to the wall function method for single-phase, turbulent flow in wall-adjacent cells. The present data for the velocity and temperature fields at the higher Reynolds number ($Re_*=590$) clearly confirm the existence of a well-pronounced logarithmic sublayer (thickness about 220 in shear-based units) in both the interface and wall regions, and for both phases.

The log laws derived in the current study can be directly applied to specify the values of the velocity and temperature in the first cell next to the interface, assuming that this cell lies within the log-law region (i.e. $30<z_{\tau ,i}^+<100$). For the turbulence quantities, appropriate interface functions still need to be developed. To this purpose, the computed turbulent kinetic energy (TKE) budgets could provide invaluable assistance (Apanasevich 2019). Additionally, the wall-function expressions for TKE and turbulent eddy dissipation rate, which have been successfully applied in CFD simulations of single-phase flow in large-industrial applications (Pope 2000), could also provide a good starting point here.

5 Conclusions

Direct Numerical Simulation (DNS) and highly-resolved Large Eddy Simulation (LES) computations have been performed to investigate condensation heat transfer in turbulent, counter-current, liquid/vapour flow in a horizontal, rectangular channel. The two phases (both assumed incompressible), driven by a constant mean pressure gradient and the gravitational force, are separated from each other by a deformable liquid/vapour interface, enclosed between confining walls, with subcooled liquid below and saturated vapour above. In the present study, two friction Reynolds numbers have been considered: $Re_*=178$ and $Re_*=590$. The choice of the maximum value of the Reynolds number was dictated by the available computational resources. For the lower Reynolds number, water subcoolings of 0K, 10K and 40K have been investigated. In the one simulation at the higher Reynolds number, only the case of 10K subcooling is included, this again as a consequence of very high computation effort involved. The numerical simulations were performed using the in-house 3D CFD code PSI-Boil. A single set of Navier–Stokes equations, coupled with an enthalpy conservation equation, together with a sharp-interface phase change model, have been solved simultaneously in the transient simulations. Within PSI-Boil, the interfacial mass transfer rate across any liquid/vapour boundary is computed directly from the local heat flux using a thermal equilibrium condition. The principal results from the current analysis are summarised below.

Velocity log-laws for both the liquid and vapour phases adjacent to solid boundaries have been confirmed from the DNS and LES data obtained. However, with increase of Reynolds number, the coefficients in the log-laws change: the slope ($\gamma = 1/\kappa$) decreases, while the shift increases. In addition, in the case of the higher Reynolds number ($Re_*=590$), a “defect” sublayer has been identified both in the interface and wall regions, and for both phases. Further analysis has indicated that, in the liquid phase, the wake strength parameter $\varPi$ is 0.45 and 0.15 in the interface and wall regions, respectively. On the vapour side, the corresponding values are 0.15 and 0.23. This is in sharp contrast to single-phase channel/pipe flows, for which the wake strengths are considerably smaller, i.e. $\varPi \approx 0$. This feature is important for the characterization of the mean velocity profiles.

In the interface region on the liquid side, the turbulence characteristics differ, depending on whether the interface is smooth or wavy. For a smooth interface, the fluctuations in the (liquid) streamwise velocity component decrease noticeably as the interface is approached. This damping effect is enhanced as the degree of liquid subcooling is increased. By contrast, for a wavy interface, the damping effect is considerably less pronounced, if it occurs at all. It is also worth noting that, at a wavy interface, Reynolds stresses do not vanish.

In contrast, on the vapour side, the velocity fluctuations are almost zero near the interface for the lower Reynolds number case, indicating that a moving, smooth interface acts on the vapour side much like a non-slip solid boundary. With increased Reynolds number, velocity fluctuations increase significantly at the interface, overcoming somewhat the interface damping effect.

Analysis of the turbulence characteristics suggests that, at the interface, the turbulent kinetic energy (TKE) is a strong function of the Reynolds number. Hence, proper modelling of this parameter within a RANS formulation, should take into account this intrinsic dependency.

Comparison with other DNS data obtained from the literature for this configuration has shown that, near the walls, for both the liquid and vapour phases, the behaviour of the velocity fluctuations is very similar to that observed for single-phase, turbulent, wall-bounded flow.

The temperature log-laws derived on the liquid side reveal that the von Karman constant for the thermal field, $\kappa _\theta$, and the shift or offset both increase with increase of $Re_*$ for both the interface and wall regions. Temperature fluctuations $\langle \theta ' \rangle ^+$, streamwise $\langle u' \theta ' \rangle ^+$ und interface-normal $\langle w' \theta ' \rangle ^+$ heat fluxes also increase significantly close to the interface with increasing Reynolds number.

Also, the present simulations have revealed that, with increase of water subcooling, the drag coefficient, $C_D$, increases non-linearly. In the presence of interfacial waves, $C_D$ is slightly lower than in the corresponding case with a smooth interface. It is suggested that part of the energy of the flow is being used to create and sustain the waves themselves. It is also proposed here that accurate modelling of the total drag force (at least in the present context) should incorporate the fact that $C_D$ is a function of both the Reynolds number and the degree of water subcooling.

Overall, analysis of the present DNS and LES computational data has indicated that the turbulent Prandtl number $Pr_t$ is very close to unity within the viscous interface ($z_{\tau ,i}^+<50$) and viscous wall ($z_{\tau ,w}^+<50$) regions for both the Reynolds numbers considered. Results confirm that the simple assumption of a constant turbulent Prandtl number between 0.7 and 1, which is widely accepted within the RANS turbulence modelling community, appears to be quite reasonable for the present counter-current, two-phase flow case.

In the general context, the DNS and LES databases derived from this study can be used to improve the modelling of stratified two-phase flows in the presence of phase change within a RANS framework.

Appendix

1.1 A Subgrid-Scale Viscosity and Thermal Conductivity in “590_10K” Case

The eddy viscosity profiles for both phases are presented in Fig. 20 in order to discuss the relationship between the modelled and resolved scales. Figure 20a shows the profiles of the resolved $\mu _t^{res}$ and modelled $\mu _t$ eddy viscosities in the interface region of the liquid phase. The black dot-dashed line indicates the molecular viscosity. It is evident that $\mu _t^{res}$ is almost two to three orders of magnitude greater than both the molecular and modelled viscosities. Thus, the contribution of $\mu _t$ deriving from the modelled scales is very small in the interface region. In the wall region of the liquid phase, Fig. 20b, the resolved eddy viscosity is also considerably greater than the molecular and modelled (eddy) viscosities. Inspection of the eddy viscosity on the vapour side shows that in the interface and wall regions the modelled viscosity is considerably smaller than the molecular viscosity.

Fig. 20

“590_10K” case: eddy viscosity profiles in the liquid interface region (a), liquid wall region (b), vapour interface region (c) and vapour wall region (d)

The ratio of the modelled turbulent eddy viscosity to the total viscosity ($\mu + \mu _t$) in the liquid and vapour phases is shown in Fig. 21a and b, respectively. On the liquid side, $\varPsi _{\mu _t}$ reaches a maximum value of 0.13, while in the vapour phase the maximum is 0.12. From this analysis, it can be clearly seen that the contribution of the modelled part is very small in comparison with the resolved part. Figure 21c presents the ratio of the turbulent thermal conductivity to the total thermal conductivity (recall that the latter is the sum of the molecular thermal conductivity and the turbulent thermal conductivity). It can be observed that the ratio, $\varPsi _{\lambda _t}$, is negligibly small at and close to the interface. Away from the interface, it attains a maximum value of 0.125. In the wall region, it is even smaller, with a maximum value of 0.1. These values indicate that the turbulent thermal conduction is significantly smaller than the molecular thermal conduction. Thus, the contribution from the unresolved scales to the total heat transport is also small, and so the present computation can be regarded as a highly-resolved LES of the temperature field.

Fig. 21

“590_10K” case: the ratio of the turbulent viscosity to the total viscosity in the liquid (a) and vapour (b) phases; the ratio of the turbulent thermal conductivity to the total thermal conductivity in the liquid phase (c)

1.2 B Verification and Validation of the Isothermal Case “178_0K”

For the verification of the DNS approach used, and the results obtained, numerical predictions for the isothermal case, “178_0K”, have been compared against data from the literature. To this purpose, the DNS studies of Lombardi et al. (1996) and Kim et al. (1987) were selected. Lombardi et al. (1996) investigated isothermal, counter-current gas/liquid flow in a channel at $Re_*=171$, cf. Fig. 22a. In this work, free-slip conditions were imposed at the bottom and top boundaries of the domain; i.e. the shear wall effect on the velocity field was not simulated. Consequently, these data can only be used for comparisons in the interface region. For the verification of the present results in the wall region, DNS results for turbulent, single-phase, closed-channel flow at $Re_*=180$, obtained by Kim et al. (1987) were utilised, cf. Fig. 22b. In addition, the experimental data of Komori and Ueda (1982) were used for the comparisons in the wall region, cf. Fig. 22c. In this experiment, the velocity fluctuations were measured in an open-channel flow near and at the free surface. In contrast to the case “178_0K”, no shear stress was imposed at the interface; i.e. $\tau _i \approx 0$, in the experiment of Komori and Ueda (1982). The zero-shear interface corresponds to a free-slip condition.

Fig. 22

Two-phase flow configuration in the DNS study of Lombardi et al. (1996) (a); single-phase, closed-channel DNS setup of Kim et al. (1987) (b); and schematic of the open-channel flow configuration studied experimentally by Komori and Ueda (1982) (c)

The auto-correlation functions for the velocity components are shown for the liquid and vapour phases in Figs. 23 and 24, respectively. Three $x-y$ planes are arbitrarily selected for the determination of the auto-correlation functions: the first is located close to the interface, at $z_{\tau ,i}^+=12$; the second near the centre line at $z_{\tau ,i}^+=146$; and the last close to the wall at $z_{\tau ,i}^+=350$ (or $z_{\tau ,w}^+=12$ in terms of dimensionless distance from the wall). It can be clearly seen for both phases that the auto-correlation functions fall off to zero values in the first half of the channel indicating that the computation domain is sufficiently large in the streamwise and spanwise directions. Figure 25 shows the ratio of the Kolmogorov scale to the $\varDelta z$, cell size in the wall normal/interface direction. Note that $z/h=0$ corresponds to the interface, while $z/h=2$ indicates the wall. The minimum value of $\eta _k/\varDelta z$ is about 1.75 in both phases.

Fig. 23

“178_0K” case: the auto-correlation functions for the fluctuations of the liquid velocity components: (a, b, c) streamwise direction and (d, e, f) spanwise direction

Fig. 24

“178_0K” case: the auto-correlation functions for the fluctuations of the vapour velocity components: (a, b, c) streamwise direction and (d, e, f) spanwise direction

Fig. 25

“178_0K” case: the ratio of the Kolmogorov scale $\eta _k$ to $\varDelta z$ (cell size in z direction): liquid phase (a) and vapour phase (b)

The profiles of the mean liquid and vapour streamwise velocities for the case “178_0K” (black line) at the interface and wall regions are presented in Fig. 26. For the sake of clarity, the values of the Reynolds number are not included in the curve labelling. Firstly, the liquid velocity profiles in the interface region, normalised in terms of the liquid interfacial shear velocity $u_{\tau ,i}$, are discussed, see Fig. 26a. We recall that in the present study the density ratio is $R=5.54$. The velocity profiles calculated by Lombardi et al. (1996) for the two density ratios $R=1$ and $R=10$ are given by the red and blue curves, respectively. As can be seen, the profiles of the three cases are almost the same in the buffer and log-law regions of the flow, i.e. from $z_{\tau ,i}^+=10$. Close to the interface, $z_{\tau ,i}^+<10$, an offset occurs, as expected. Lombardi et al. (1996) have reported an interfacial velocity $\langle U \rangle _i^+=2.298$ and $\langle U \rangle _i^+=0.002$ for the $R=10$ and $R=1$ cases, respectively, while in the case “178_0K” a value of $\langle U \rangle _i^+=1.18$ was calculated. Evidently, with the decrease of density ratio (implying an increase of the interactions between the phases), the interfacial velocity decreases on the liquid side. It is important to note that here all three profiles coincide for $z_{\tau ,i}^+<10$, but a velocity offset appears in the logarithmic region $z_{\tau ,i}^+>30$ if the interfacial velocity is subtracted out (cf. Fig. 26b). As explained by Lombardi et al. (1996), this behaviour indicates that the velocity profiles in the logarithmic region have the same slope, but differ in magnitude by the value of the interfacial velocity. Comparisons for the wall region are presented in Fig. 26c. The normalisation is in terms of the wall shear velocity $u_{\tau ,w}$. It is clearly seen that there is no discernible difference between the “0K” velocity profile obtained here and that of Kim et al. (1987). The mean velocity profiles in the vapour phase are displayed in Fig. 26d, e. The profiles are very similar for all simulations in the interface region. In the wall region (cf. Fig. 26e), the “0K” vapour velocity profile collapses onto the single-phase profile produced by Kim et al. (1987).

Fig. 26

“178_0K” case: mean streamwise velocity profiles in the vapour and liquid phases

RMS results from the case “178_0K” are compared again with the DNS data of Lombardi et al. (1996), only for the case “R=10”, and the data of Kim et al. (1987) in the same way, as already shown for the mean velocity profiles. Figure 27 presents the velocity fluctuations in the liquid phase. Normalisation is performed with respect to the liquid interfacial and wall shear velocity in the interface and wall regions, respectively.

Fig. 27

“178_0K” case: liquid root-mean-square velocity fluctuations in the interface and the wall regions

In the interface region (cf. Fig. 27a–c), the calculated quantities show good agreement with the data of Lombardi et al. (1996). However, there are some quantitative and qualitative deviations observed close to the interface and near the centre line that need to be explained. In the reference case “Lombardi R=10”, in the region close to the interface, a higher maximum value of $u_{rms}$ is predicted, Fig. 27a. Additionally, for $z_{\tau ,i}^+<10$, the “0K” profile decreases monotonically, while the “R=10” curve reaches a plateau. This indicates that the interface reduces the rms fluctuations of the streamwise velocity, and this damping effect increases as the density ratio decreases. The profiles of the spanwise velocity fluctuations, shown in Fig. 27b, show a similar tendency: the smaller the density ratio, the stronger damping effect at the interface. Concerning the third quantity, $w_{rms}$, there is a very good agreement between both the “0K” and “Lombardi R=10” profiles in the interface region for $z_{\tau ,i}^+<130$. Due to gravity and surface tension, the vertical velocity fluctuations vanish at the interface. From $z_{\tau ,i}^+=130$, in the “Lombardi R=10” case, the $w_{rms}$ values decrease rapidly to zero at $z_{\tau ,i}^+<170$, whereas in the “178_0K” case the vertical velocity fluctuations reduce considerably slower to a minimum value. The different behaviour of $w_{rms}$ in both cases is ascribed to the fact that in the “Lombardi R=10” simulation a free-slip condition was imposed at the outer domain boundaries (i.e. at $z_{\tau ,i}^+=170$), which forces $w_{rms}=0$ there. The decrease of $w_{rms}$ to zero, and increase of $u_{rms}$ and $v_{rms}$ as $z_{\tau ,i}^+\rightarrow 170$, indicate a redistribution of turbulent energy from the vertical velocity component to the streamwise and spanwise components in the case “Lombardi R=10”. This observation is in line the works of Hunt and Graham (1978) and Brumley and Jirka (1988). In the computation “178_0K”, the wall region is also considered. For this reason, and contrary to the simulation “Lombardi R=10”, the vertical velocity fluctuations $w_{rms}$ do not vanish in the bulk of liquid, and consequently $u_{rms}$ and $v_{rms}$ do not increase near the centre line.

In the wall region, Fig. 27d–f, there is good agreement between the results of the “178_0K” simulation and the reference data for all three fluctuating velocity components. Evidently, near the wall, the maximum values of the velocity fluctuations are lower in the two-phase flow case than those for single-phase channel flow. Note that in the case “178_0K” the interface is a deformable and moving boundary within which the fluid properties change rapidly. This is very different to the single-phase flow conditions studied by Kim et al. (1987), who considered only rigid boundaries. Comparison with the experimental data of Komori and Ueda (1982) (green symbols) offers the possibility for validation of the “178_0K” case in the outer sublayer of the wall region; i.e. for $z_{\tau ,w}^+>85$. Obviously, good agreement can be observed for $u_{rms}$ and $v_{rms}$, whereas $w_{rms}$ agrees well for $z_{\tau ,w}^+<150$. For $z_{\tau ,w}^+>150$, in the experiment of Komori and Ueda (1982), the vertical velocity fluctuations decrease rapidly as the free surface is approached; i.e. the same behaviour as $w_{rms}$ in the “Lombardi R=10” case. The shapes of the all three experimental curves also confirm that the damping of $w_{rms}$ close to the free surface leads to an enhancement of the two other fluctuating velocity components (i.e. a redistribution of turbulent energy between the velocity components).

Changing now to the vapour side, Fig. 28a–f show that the profiles of the turbulence intensities in the interface and near-wall regions behave very similarly. Comparing the “10K” and “Lombardi R=10” results are in very good agreement for $v_{rms}$ and $w_{rms}$, except for the flow region near $z_{\tau ,i}^+=150$ (same reasons for this discrepancy as on the liquid side). For $u_{rms}$, different maxima are recorded due to the different density ratios in the “178_0K” and “Lombardi R=10” runs. In the near-wall region (cf. Fig. 28d–f), in comparison with the liquid side, a similar level of agreement is observed between the “10K” curves and the results of Kim et al. (1987) and data of Komori and Ueda (1982). It is worth noting that in the case “178_0K” the difference between the maxima of the rms values near the interface and the wall is negligibly small compared to the liquid side. This observation confirms that the interface acts like a rigid wall to the vapour phase.

Fig. 28

“178_0K” case: vapour root-mean-square velocity fluctuations in the interface and the wall regions

Comparing the mean velocities and turbulence intensities in the two phases, it may be concluded that the results of the “178_0K” case agree well with the reference data from the literature. The main differences are mostly caused by the use of different boundary conditions and density ratio.

1.3 C Qualification of LES Approach Against the DNS Data

Before simulating condensing two-phase flow at the higher Reynolds number, the potential of the LES method for the simulation of such flows was investigated. To this purpose, a single-phase channel flow at a friction Reynolds number of 590 was simulated. The size of the computational box is $L_x=2.15\pi h$ in streamwise direction and $L_y=1.34\pi h$ in spanwise direction, which is slightly larger than the one used in the reference DNS study (Moser et al. 1999), for which $L_z=2$, where h is the half channel width. This yields a grid resolution of $\varDelta x^+=15.2$,$\varDelta y^+=12.9$ and $\varDelta z^+=0.65-16.6$ in wall units in the three coordinate directions. The total number of cells is then approximately 12 583 000. Thus, the resolution of this grid satisfies the conditions for well-resolved LES, as specified by Piomelli et al. (1990). Note that the grid is stretched in the wall-normal direction by different factors. The stretching factor is less than 1.05 for the entire domain, which follows the recommendations of Fröhlich (2006).

The averaging time of all the calculations presented below is about $3000\nu /u_*^2$ in dimensionless time units. The time step is about $0.2\nu /u_*^2$ in the simulations with the coarse grid, and $0.04\nu /u_*^2$ using the refined grid. This corresponds to a maximum CFL number of approximately 0.3 in all the simulations. The simulation is carried out using 192 cores (Node configuration: 2x Intel(R) Xeon (R) CPU E5-2690 (8 cores), 2.90 Hz, 2 GB RAM per core).

Firstly, the influence of the WALE model constant on the flow structures was investigated. Four values are tested here: $C_W=0.325$, 0.45, 0.5 and 0.6. From analysis of the effect of variation of $C_W$, it can be concluded that only the streamwise fluctuating velocity component $u_{rms}^+$ exhibits a strong dependency on the $C_W$ value. The smallest differences between LES and DNS were observed when using $C_W=0.45$, which is the lower limit of $C_W$ derived by Nicoud (1998). Only these results are presented in the following section.

An important property of the WALE model is its ability to predict the correct asymptotic behaviour of the subgrid viscosity in the wall region. The dynamic eddy viscosity is plotted as function of the dimensionless distance from the wall $z^+$ in Fig. 29a. The solid black line shows the theoretical near-wall scaling of the eddy viscosity ($\mu _t \sim (z^+)^3$), while the dot-dashed line denotes the molecular dynamic viscosity. Evidently, the correct asymptotic behaviour of the subgrid viscosity is reproduced using the WALE model.This quantity is less than the molecular viscosity in the entire domain. The eddy viscosity is also seen to be negligible within a thin, wall-adjusted layer of thickness 10 wall units.

Fig. 29

Computed turbulent viscosity (a) and the ratio of the turbulent viscosity to the total viscosity (b)

To give an impression of the importance of the SGS model, the ratio of the turbulent viscosity to the total viscosity, including both the eddy and molecular viscosity contributions, can be used (Hertel et al. 2013):

$$\begin{aligned} \varPsi _{\mu _t}=\frac{\langle \mu _t \rangle }{\langle \mu _t \rangle + \mu } \end{aligned}$$

(40)

Figure 29b plots the quantity $\varPsi _{\mu _t}$ in global coordinates; $z/h=0$ corresponds to the wall. The diagram indicates that the eddy viscosity contributes only up to 20% to the total viscosity.

The dimensionless mean velocity is displayed in Fig. 30a. The plot shows very good agreement between the LES and DNS: the error is less than about 1.8%. Figure 30b shows the normalised RMS values of the velocity fluctuations. In the near-wall region, the streamwise velocity fluctuations are calculated to an accuracy of better than 1%, while in the outer region the maximum error is about 3%. For the spanwise velocity fluctuations, the error is 9.9% in the near-wall region, and less than 1% in the channel core. As for the spanwise velocity fluctuations, the values of $w_{rms}^+$ differ from the reference values only in the near-wall region; the discrepancies are about 6%. The Reynolds stress variation is given in Fig. 30c. There is excellent agreement between the LES and DNS predictions.

Fig. 30

Mean streamwise velocity profiles (a), root-mean-square velocity fluctuations (b) and Reynolds stress (c)

In a second validation step, LES of condensing counter-current liquid/vapour flow at a friction Reynolds number of 178 is conducted and compared with our own two-phase DNS data. The domain size, flow parameters, and boundary conditions are the same as those for the reference DNS case with a water subcooling of 10K.

For the discretisation of the computational domain, two grid resolutions are examined. The base grid (M1) has 96 cells in the streamwise direction and 64 cells in the spanwise direction. To resolve the near-wall and interface sublayers, the grid is stretched in the wall/interface-normal directions, and consists of 192 cells in total. The grid spacing is $\varDelta x^+=22.2$, $\varDelta y^+=11.6$ and $\varDelta z^+=1.0-5.1$ in the streamwise, spanwise and vertical directions, respectively. Note that a cell size $\varDelta z^+=1.0$ is imposed at the interface and the walls, while $\varDelta z^+=5.1$ refers to the grid resolution at the centre line for each phase. The refined grid (M2) comprises 192 cells in the streamwise direction, 96 cells spanwise, and 256 cells in the vertical direction. In the z direction, the grid is stretched using a factor smaller than 1.05. This results in a grid resolution of $\varDelta x^+=11.1$, $\varDelta y^+=7.8$ and $\varDelta z^+=0.6-7.9$ in three coordinate directions. The distribution of the cells from the interface to the centre line is summarised for both LES grids in Table 6.

Table 6 The cell distributions in the vertical direction

The WALE model constant is set here to $C_W=0.45$ and a turbulent Prandtl number of 0.9 is adopted. In the simulation with the base grid M1, the non-dimensional time step is set to $0.004\nu /u_*^2$, whereas for the refined grid M2 it is around $0.004\nu /u_*^2$ wall units. In order to calculate the turbulence statistics, transient data are collected over $4000\nu /u_*^2$ wall-time units, with a time interval of $2.44\nu /u_*^2$ in both simulations. For the simulations with the base and refined grids, 32 and 64 processors are used, respectively.

Energy spectra for the U and V velocity components and temperature, obtained with the refined grid M2, are shown for the liquid and vapour phases in Figs. 31 and 32, respectively. Three $x-y$ planes are arbitrarily selected for the determination of the spectra: the first is located close to the interface, at $z_{\tau ,i}^+=12$; the second near the centre line at $z_{\tau ,i}^+=146$; and the last close to the wall at $z_{\tau ,i}^+=350$ (or $z_{\tau ,w}^+=12$ in terms of dimensionless distance from the wall). The wave number used is $k=2\pi /b$, where b is the wavelength. The subscripts x and y denote wave numbers in the streamwise and spanwise directions, respectively. For the sake of clarity, when presenting the DNS spectra, only every fourth point is plotted. The energy spectra of the vertical velocity component W are not shown here, since they are very similar to those of the spanwise velocity V. It can be seen that the spectra obtained using LES agree very well with those obtained with DNS. The $E_i^y$ spectra exhibit somewhat faster decay for high wave numbers compared to the spectra obtained using DNS. Furthermore, all the LES spectra clearly show that there is no energy accumulation at the high wave numbers. From minimal observed differences between the DNS and LES spectra, it might be concluded that the eddy viscosity, computed by the SGS model of the LES computation, is of the same order of magnitude as the molecular (laminar) viscosity.

Fig. 31

Energy spectra on the liquid side: the solid lines refer to the LES simulation (refined mesh M2), and the symbols to DNS

Fig. 32

Energy spectra on the vapour side: the solid lines refer to the LES simulation (refined mesh M2), and the symbols to DNS

Figure 33a, b show the SGS viscosity in the liquid phase, as obtained with the base grid M1 (blue solid curve) and the refined mesh M2 (red solid line). The black, dash-dotted line depicts the molecular viscosity for the liquid phase, $\mu _L$. The black, solid curve is the eddy viscosity obtained from the DNS calculation $\mu _t^{DNS}$, defined as:

$$\begin{aligned} \mu _t^{DNS}=\frac{-\rho \langle u'w' \rangle }{\langle \partial U / \partial z \rangle } \end{aligned}$$

(41)

Note that $\mu _t^{DNS}$ is a post-processing quantity in the context of a DNS analysis. The resolved turbulent viscosity $\mu _t^{res}$, shown by the red symbols only for the LES case with the refined grid M2, is calculated in the same manner as for $\mu _t^{DNS}$, using Eq. (41) and the resolved (filtered) variables.

Fig. 33

Eddy viscosity in the liquid and vapour phases

It can be seen from Fig. 33a that the SGS viscosity is negligibly small as the interface is approached ($z_{\tau ,i}^+\rightarrow 0$), but increases progressively with distance from the interface. Evidently, the molecular viscosity is greater than that modelled for both the LES cases. As expected, the magnitude of the eddy viscosity decreases as the degree of grid resolution is increased. Considering the shape of the eddy viscosity profiles, it can be noticed that the grid resolution has some effect on the slope of the curve in the interface region. In the “M1” case, the eddy viscosity increases more steeply with $z_{\tau ,i}^+$, while in the case with the refined grid (“M2”) the SGS viscosity profile has the same slope as the profile of $\mu _t^{DNS}$, cf. Fig.  33a. It is worth noting that the differences between the LES resolved eddy viscosity and the DNS directly calculated eddy viscosity are very small. Obviously, the magnitude of the resolved eddy viscosity is two to three orders of magnitude higher than the modelled one. In the wall region, Fig.  33b, no grid effect on the profile shape can be discerned. Both the LES profiles have the same slope as the DNS one. The resolved eddy viscosity profile obtained from LES collapses onto the DNS profile, and is also two to three orders of magnitude higher than the modelled viscosity.

Turning now to the vapour side, Fig. 33c, d, similar trends can be observed: i.e. the modelled viscosity is considerably smaller than the molecular one, and the resolved eddy viscosity is significantly greater than the modelled (SGS) and molecular viscosities. It is worth noting that LES is able to reproduce the slopes of the eddy viscosity profile on the liquid and vapour sides of the interface.

The ratio of the turbulent viscosity to the total viscosity $\varPsi _{\mu _t}$, Eq. (40), in both phases, across the channel is depicted in Fig. 34. The profiles confirm that the eddy viscosity is zero at both walls, and at the interface. For the coarse grid, the maximum value of the viscosity ratio is about 0.22 in both the liquid and vapour phases. With grid refinement, $\varPsi _{\mu _t}$ reduces to 0.13 in both the phases. From the low values of $\varPsi _{\mu _t}$ predicted over the entire domain, and from consideration of the spectra, it can be concluded that the influence of the SGS model is minimal for the velocity field simulation.

Fig. 34

Ratio of the turbulent viscosity to the total viscosity in the liquid (a) and vapour (b) phases

The mean streamwise velocity profiles on the liquid side are compared in Fig. 35. Considering the interface region, Fig. 35a, it is important to notice that there is a grid effect on the mean velocity profile only for $z_{\tau ,i}^{+}<10$. As expected, the profile obtained using the refined grid is significantly closer to the (assumed exact) DNS profile; note that the refined grid, M2, contains more grid points in the region $z_{\tau ,i}^+<30$, Table 6. For $z_{\tau ,i}^{+}>10$, differences between the LES profiles and that obtained using DNS are negligible. In the wall region, Fig. 35b, on the other hand, the LES profiles collapse onto the DNS profile.

Fig. 35

Mean streamwise velocity profiles on the liquid and vapour sides; LES vs. DNS

In contrast to the liquid side, only a very weak grid effect on the vapour mean velocity profile is observed for $z_{\tau ,i}^{+}<10$ in the interface region, Fig. 35c. In the logarithmic sublayer region, $z_{\tau ,i}^{+}>30$, the velocity differences peak at 5% and 2.5% for the base M1 and refined M2 grids, respectively. In the wall region, Fig. 35d, all the profiles collapse onto a single curve for $z_{\tau ,w}^{+}<30$. In the logarithmic sublayer region, the differences are less than 2% for the M1 and M2 cases, respectively.

The next important consideration is the root-mean-square velocity fluctuations. Figure 36 shows comparisons for the liquid phase. From a qualitative point of view, there are no differences between the LES and DNS results. The LES prediction comes closer to the DNS result as the degree of grid resolution is increased. A quantitative analysis reveals that the quality of the predictions of peaks in $v_{rms}^+$ and $w_{rms}^+$ are more sensitive to the grid resolution. For the base grid, M1, in the interface region, the peaks of $v_{rms}^+$ and $w_{rms}^+$ are underpredicted by 9.6% and 8.1%, respectively. For $u_{rms}^+$, there is an overprediction of about 5.8%. In the wall region, the predictions of $u_{rms}^+$, $v_{rms}^+$ and $w_{rms}^+$ are generally better by 1–2% than those in the interface region. For the refined grid, M2, $v_{rms}^+$ and $w_{rms}^+$ are predicted with an accuracy of 3–4% in both regions, while differences in $u_{rms}^+$ are less than 3%. Evidently, the percentage error (LES vs. DNS) using the refined grid (M2) is reduced by a factor of two at least from that of the coarse grid (M1).

Fig. 36

Root-mean-square velocity fluctuations on the liquid side in the interface (a, c) and wall (b, d) regions

Figure 37 compares the rms velocity fluctuations on the vapour side. For $u_{rms}^+$, it can be seen that in the interface region the accuracy is about 2.5% for both grids, while the near-wall peak is overpredicted by about 8%. Concerning the two other velocity fluctuation components, a similar accuracy is observed as on the liquid side. For the refined grid, the accuracies of the LES predictions of $v_{rms}^+$ and $w_{rms}^+$ are better than 5% in both the interface and wall regions; Fig. 37c, d.

Fig. 37

Root-mean-square velocity fluctuations on the vapour side in the interface (a, c) and wall (b, d) regions

In Fig. 38a, b, comparisons of the liquid Reynolds stress are displayed. It can be clearly seen that the LES results agree very well with those from DNS: indeed, the discrepancies are less than 2% for the base and refined grid cases. Turning now to the vapour phase, Fig. 38c, d give the profiles of the Reynolds stresses in the interface and wall regions, respectively. As expected, the greatest discrepancies (up to 12%) occur in the interface region for the base grid case M1. For the refined grid, this error is reduced to 4%, while in the wall region these discrepancies are reduced: 1.6–0.8%.

Fig. 38

Reynolds stress profiles on the liquid and vapour sides of the interface

Turning now to the analysis of the temperature field, Fig. 39a shows the ratio of the turbulent thermal conductivity to the total thermal conductivity (i.e. the sum of the molecular and SGS thermal conductivities), $\varPsi _{\lambda _t}$. By analogy to $\varPsi _{\mu _t}$, $\varPsi _{\lambda _t}$ provides an estimation of the role of the SGS model used on the computation of the unresolved heat flux. For the coarser resolution case (M1), $\varPsi _{\lambda _t}$ reaches a maximum of 0.16 (16%) and 0.20 (20%) in the wall and interface regions, respectively. In the refined mesh case (M2), the maximum of $\varPsi _{\lambda _t}$ is significantly lower: 0.10 (10%) and 0.12 (12%) in the wall and interface regions, respectively.

Fig. 39

Ratio of the turbulent thermal conductivity to the total thermal conductivity in the liquid phase (a); mean temperature profiles in the liquid phase in the interface (b) and wall (c) regions

The mean temperature profiles are compared in Fig. 39b, c. In the interface region, the DNS and LES profiles collapse on each other in the vicinity of the interface: i.e. for $z_{\tau ,i}^{+}<10$. This illustrates that both LES calculations are able to predict the thickness of the temperature diffusion sublayer (a zone where the temperature changes linearly, and heat is transferred only by molecular thermal conduction), which is about 4-4.5 wall units. For $z_{\tau ,i}^{+}>30$, the temperature predictions are within 12.5% and 4.6% of the DNS values for the coarse and refined grids, respectively. In the wall region, the agreement between LES und DNS curves is better still: differences are less than 2%. It is worth noting that the temperature gradient at the interface is underpredicted by 3.25% and 6.25% for the base M1 and refined M2 grid cases, respectively. The higher temperature gradient at the interface for the coarser grid can be explained by the greater overprediction of the streamwise velocity fluctuations, which results in higher turbulent heat diffusion in the interface region. Due to this, steeper temperature gradients are observed in the interface region for the coarser grid. The predicted condensation rates for the coarser and refined grids are underpredicted by 3.25% and 6.25%, respectively.

Finally, the resolved rms temperature fluctuations and the resolved streamwise and interface/wall-normal turbulent heat fluxes are shown in Fig. 40. Evidently, the intensity of the fluctuating temperature is overestimated in the interface and wall regions for the coarser grid (M1), cf. Fig. 40a, b; the discrepancies are up to 12.4%. For the finer mesh resolution (M2), the temperature fluctuations agree very well with those derived using DNS: differences are just less than 3.3%.

Fig. 40

Root-mean-square temperature fluctuations in the interface (a) and wall (b) regions; distributions of the streamwise turbulent heat flux in the interface (c) and wall (d) regions; and distributions of the interface/wall-normal turbulent heat flux in the interface (e) and wall (f) regions

Examination of the heat flux predictions shows that the streamwise heat fluxes, Fig. 40c, d, are more sensitive to the degree of grid resolution. The maxima occurring near the wall and interface for the base grid (M1) are overestimated by 14% and 17.3%, respectively. For the refined grid, the discrepancies are up to 6.3% and 7.2% in the wall and interface regions, respectively. Considering the turbulent heat flux in the direction normal to the interface/wall, Fig. 40e, f, it can be seen that for the coarse grid (M1) the fluxes are overestimated by approximately 7% in the far-field regions, $z_{\tau ,i}^{+}>30$ and $z_{\tau ,w}^{+}>30$. However, the LES predictions improve with increase in grid resolution (M2): discrepancies are then less than 3%.

The above promising results allow finally to conclude that the LES approach, in combination with the WALE SGS model, can be applied to simulate stratified two-phase flows, including phase change, with acceptable accuracy.