# Numerical simulations of bubbly turbulent convection in cubical geometries

J Marichal, Pierre Ruyer, Y Bartosiewicz

► To cite this version:

J Marichal, Pierre Ruyer, Y Bartosiewicz. Numerical simulations of bubbly turbulent convection in cubical geometries. Eurotherm 2024 - 9th European Thermal Sciences Conference, University of Ljubljana, Faculty of Mechanical Engineering, Slovenia, Jun 2024, Bled (Slovénie), Slovenia. pp.012030, 10.1088/1742-6596/2766/1/012030. irsn-04664110

HAL Id: irsn-04664110

https://asnr.hal.science/irsn-04664110v1

Submitted on 9 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical simulations of bubbly turbulent convection in cubical geometries

J Marichal¹, P Ruyer² and Y Bartosiewicz¹

¹ Université catholique de Louvain, Institute of Mechanics, Materials and Civil Engineering, Louvain-la-Neuve 1348, Belgium

² Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-RES/SEMIA/LSMA, Cadarache, 13115 Saint-Paul-Lez-Durance, France

E-mail: joauma.marichal@uclouvain.be

**Abstract.** In this work we present numerical results of pool boiling flow in a turbulent Rayleigh-Bénard convection configuration, using our in-house code in a cubical geometry. The problem in hand is encountered in various natural phenomena as well as in industrial applications. An Eulerian-Lagrangian approach is developed for the mixture of liquid water and vapor bubbles. The liquid mean temperature is close to the saturation temperature and is governed by the quasi-incompressible Navier-Stokes equations that are solved using DNS standards. The motion and growth/shrinkage of each individual vapor bubble is modeled and the backward effect of the bubbles on the fluid is accounted via momentum and energy exchanges between the two phases (two-way coupling), as well as variations in fluid phases volumetric fractions (volumetric coupling). In such pool configuration the non-dimensional parameters governing the flow are both those relative to Rayleigh-Bénard convection, namely Rayleigh number, Prandtl number and aspect ratio, and those to boiling, namely the overall vapor volumetric fraction, the bubble size based Reynolds number, the degree of superheat of the liquid scaled by the overall temperature difference, and the Jakob number (sensible heat to latent heat ratio).

At first, we describe the model used and its corresponding validation, involving natural convection, isolated bubble dynamics and coupling between bubbly and bulk flows. In the second part, we consider the study of the relationship between heat transfer through the pool (Nusselt number) and the flow topology for different settings of the bubbly configurations.

## 1. Introduction

Rayleigh-Bénard convection is a well-known buoyancy driven flow of a fluid enclosed between a lower hot wall, at temperature, $T_{\text{hot}}$ and an upper cold one, at temperature $T_{\text{cold}}$. It has been extensively studied in the literature; see for example [2, 3, 10, 14, 16, 17, 24] and references therein. These studies have shown that the problem is governed by three dimensionless groups: the Rayleigh $Ra$ and Prandtl $Pr$ numbers as well as the aspect ratio $\Gamma$. Research on the subject has also provided important information and physical insight on the flow topology, the properties of thermal and hydrodynamic boundary layers and the turbulence statistics. Correlations between the Nusselt and Rayleigh numbers have also been developed [14, 9].

This study aims at investigating the same problem in the presence of boiling, which is a common phenomenon in various technological applications. We therefore perform Direct Numerical Simulations (DNS) of a two-phase flow composed of water, close to its saturation

<figure><img src="image_2.png" /></figure>
temperature, and vapor bubbles, namely a pool boiling configuration. We rely on an Eulerian-Lagrangian approach which, on the one hand, solves the low-Mach number Navier-Stokes equations for the fluid phase and, on the other, tracks each bubble movement and size. The volume occupied by the bubbles is perceived by the liquid via a space- and time-dependent volumetric fraction (volumetric coupling) while both the force that a bubble applies on, and the energy that a bubble exchanges with the surrounding liquid are accounted for as additional source terms (two-way coupling). The addition of bubbles leads to an additional dimensionless parameter which controls the bubble growth rate, the Jakob number, $Ja = \frac{\rho c_p (T_{hot} - T_{sat})}{\rho_b h_{lv}}$, where $T_{sat}$, $\rho_b$ and $h_{lv}$ are the saturation temperature, density of bubbles and the latent heat of vaporization at a constant pressure.

The next section will describe the governing equations that have been implemented in our in-house code. We will then focus on the validation procedure before analyzing the flow topology and turbulence statistics for various Jakob numbers.

## 2. Governing equations

Boiling Rayleigh-Bénard convection in a closed rectangular geometry is studied by adopting an Eulerian-Lagrangian approach. The fluid is governed by the quasi-incompressible [8, 19] Navier-Stokes equations [12, 18]. The effects of vapor bubbles on the fluid are reflected through the liquid phase volume fraction and additional momentum and energy terms [4, 22]. Accordingly, a general transport equation for a generic quantity $\eta$ reads, in dimensional form,

$$ \frac{\partial(\phi\eta)}{\partial t} + \nabla \cdot (\phi\eta\mathbf{u}) = \nabla \cdot (\phi\Gamma\nabla\eta) + \phi S + f_b, \quad (1) $$

where $\Gamma$ is the diffusion coefficient for $\eta$. To obtain the mass, momentum and energy equations, $\eta$ is set equal to $\rho$, $\rho\mathbf{u}$ and $\rho c_p T$ respectively where $\rho$, $\mathbf{u}$ and $T$ stand for the density, the velocity vector and the temperature of the fluid. The fluid specific heat capacity is denoted by $c_p$ and the additional source term due to the presence of bubbles by $f_b$. The liquid volumetric fraction, $\phi$ is coupled to the bubble void fraction, $\phi_b$ through the following relation: $\phi + \phi_b = 1$.

The bubbles are subjected to the buoyancy force, the drag force, the lift force, the added mass force and a force representing the momentum transfer due to radius variations [6, 11, 13, 15, 22]. The bubble radius is obtained by balancing the heat due to evaporation/condensation and the heat exchanged with the liquid [11, 15] while assuming a constant bubble density (at saturation value) and that the bubble keeps a spherical shape. Equations for the bubble position, velocity and radius can be found in the literature [11, 15, 22].

An interpolation operator is needed to couple the two reference frames (respectively Eulerian and Lagrangian) in a smooth manner, ensuring the conservation of transferred variables [4, 22]. The mapping from the Lagrangian to the Eulerian frame relies on Lagrange polynomials. For the reverse transformation, we follow the work of [4] and use a clipped fourth-order polynomial function, denoted by $\omega(n)$. The additional momentum and energy terms accounting for the backward effects of bubbles on the fluid and the fluid void fraction are therefore expressed as,

$$ f_u(j) = \frac{\sum_n f_n \int_{\Omega_j} \omega(n) d\Omega}{V_{cell}} \quad (2) $$

$$ f_T(j) = \frac{\sum_n Q_n \int_{\Omega_j} \omega(n) d\Omega}{V_{cell}} \quad (3) $$

$$ \phi(j) = 1 - \frac{\sum_n V_n \int_{\Omega_j} \omega(n) d\Omega}{V_{cell}} \quad (4) $$

where $f_n$ is the force applied by the $n^{th}$ bubble on the surrounding fluid (cell $j$) and $Q_n = (4\pi r_n^2 h_b (T_{sat} - T_f|_n))$ is the energy given/taken by the $n^{th}$ bubble on the neighboring
fluid due to phase change. $V_n$ and $V_{cell}$ denote the $n^{th}$ bubble's volume and the volume of cell $j$ respectively.

### 3. Validation

A series of test cases are carried out to validate the code used. A first test case is performed to analyze the motion of a single bubble in a two-dimensional Taylor-Green vortex [23]. The single bubble is initially located at the corner of one of the vortices and is subject to drag, lift and added-mass forces. The bubble is found to follow a spiral trajectory towards the center of the domain [21], driven mainly by the drag and added-mass forces.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Instantaneous contours of the vorticity $\omega/\omega_0$ at time $t_{ff} = 0.2$ for two-way coupling.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Normalized bubble void fraction at the vortex center for one-way coupling (—) and two-way coupling (---). Results for one-way coupling and two-way coupling from [7] (---- and ---) and [22] (--- and ---) are added for comparison.</figcaption></figure>

A second case is performed to study the effects of bubbles on the liquid phase: 230 400 bubbles are evenly distributed in the same 2D Taylor-Green vortex as before, leading to a uniform bubble void fraction of 0.01. The bubbles migrate towards the center of the vortexal structures and modify the flow significantly. Figure 1 illustrates the latter by displaying the instantaneous contours of the vorticity. Figure 2 shows the bubble void fraction at the center of the domain with time for one-way coupling and two-way coupling (volumetric coupling and additional momentum). Our results agree well with the ones obtained by two previous studies [7, 21]. As discussed by [21], the slight differences are explained by difference of numerical schemes and interpolation methods.

The previously described validation tests considered neither buoyancy effects nor bubble growth. To appraise the gravitational forces acting on a bubble, a single bubble of fixed radius is introduced in a quiescent fluid. The bubble is shown to follow a rectilinear vertical trajectory [1] and its vertical velocity tends towards the theoretical terminal speed [21], which is illustrated on figure 3. The change in bubble radius due to pressure is then accounted for, relying on an ideal gas law model for vapour. As time passes, the bubble rises, hydrostatic pressure decreases resulting in a radius increase, thereby leading to a higher velocity. As shown in figure 3, the actual speed stays slightly behind the theoretical terminal velocity.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Velocity of a single bubble in a quiescent liquid: $u_x$ (—), $u_y$ (—) and $u_z$ (---) and the theoretical terminal velocity (—) for a bubble (a) without and (b) with radius change.</figcaption></figure>

Finally, to analyze the thermal aspect related to bubbles, we consider a single bubble at rest in a super-heated liquid, also at rest and in absence of buoyancy. This problem has been studied theoretically and experimentally by [20] and [5] and has been accurately reproduced by our code, as displayed by figure 4.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Radius of a single bubble at rest in a quiescent fluid (—) compared against the results obtained by [20] (—).</figcaption></figure>

**4. Results**

Various simulations are performed at a saturation temperature of $T_{\text{sat}} = 373.15$ K and a latent heat of $h_{lv} = 2257.92 \times 10^3$ J/kg. The Rayleigh number is set to $2 \times 10^5$ and the reference temperature, defined as the mean between the imposed values at the top and bottom, equal to the saturation one. Several Jakob number values are imposed to investigate the effect of this parameter on the flow.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** $Ja = 0.4$ (a) time-averaged streamlines colored with the time-averaged normalized velocity $\langle \frac{u}{u_{ref}} \rangle$, (b) clipped at mid-height, (c) contour plots of the mean amplitude of the in-plane time-averaged normalized velocity, with superimposed in-plane velocity vectors on a diagonal plane and (d) contour plots of time-averaged mean bubble void fraction, $\langle \phi_b \rangle$ on a diagonal plane.</figcaption></figure>

The addition of bubbles to the flow tends to modify the well-known Large-Scale circulation (of classical Rayleigh-Bénard convection in cubical geometries) to a "nearly toroidal roll" structure [16], where the fluid rises along the corners and sinks in the middle, leading to the creation of four convection cells, as shown in figure 5. Indeed, bubbles, with their low density, have a clear tendency of rising therefore resulting in an additional bubble-induced buoyant forcing [15]. On the other hand, bubbles tend to absorb/evacuate heat if the surrounding fluid's temperature is greater/lower than the saturation one, therefore lessening the temperature gradients ergo the liquid buoyancy. Lastly, the bubbles grow/shrink (more or less depending on the Jakob number) in hot/cold liquid regions which again has a direct impact on the fluid's flotability. As the Jakob number is increased, these effects take various proportions relative to one another which is reflected on the flow structure and the heat transfer through the domain.

In all cases, the mean bubble void fraction, $\langle \phi_b \rangle$, is concentrated at the corners of the domain, coinciding with the rising columns of fluid. When the bubbles are not allowed to grow/shrink ($Ja = 0$), the profile does not vary in the vertical direction and the center of the domain encounters very few bubbles. As bubbles are affected by the surrounding temperature, they concentrate on very thin columns in the corners as well as near the bottom wall, coinciding with the hot fluid zones. In the center of the domain, $\langle \phi_b \rangle$ is bigger than zero, which may be explained by the fact that a few bubbles succeed into rising against the flow direction, however slowly, thus affecting the mean value.

An investigation of turbulence statistics has shown that the thickness of the thermal boundary layer decreases as bubbles of fixed radius ($Ja = 0$) are added to the single-phase case and this trend intensifies as bubbles are allowed to grow more intensively ($Ja$ increase) [11]. As explained by Lakkaraju [11], this effect is due to the fact that bubbles are at the saturation temperature and therefore cool the liquid at the bottom/hot side and warm it near the top/cold one.
**5. Conclusion**

For this study, we implemented an Eulerian-Lagrangian method to track vapor bubbles inside a fluid close to saturation in our in-house code. The method as well as the validation procedure put in place have been briefly described and will be further developed during the presentation. An introductory analysis of the influence of bubbles, and more specifically of the Jakob number on the flow fields, highlighted not only the change in flow structure, but also the distribution of bubbles in the domain and the thinning of the thermal boundary layer.

**References**

[1] Bonnefis P, Fabre D and Magnaudet J 2023 *P. Natl. Acad. Sci. USA* **120** e2300897120

[2] Cioni S, Ciliberto S and Sommeria J 1996 *Dynam. Atmos. Oceans* **24** 117-27

[3] Cioni S, Ciliberto S and Sommeria J 1997 *J. Fluid Mech.* **335** 111-40

[4] Darmana D, Deen N G and Kuipers J A M 2006 *J. Comput. Phys.* **220** 216-48

[5] Dergarabedian P 1960 *J. Fluid Mech.* **9** 39-48

[6] Elghobashi S 1994 *Appl. Sci. Res.* **52** 309-29

[7] Ferrante A and Elghobashi S E 2007 *J. Fluid Mech.* **572** 145-77

[8] Georgiou M and Papalexandris M V 2018 *J. Fluid Mech.* **845** 581-614

[9] Hay W A, Martin J, Migot B, Papalexandris M V 2021 *Phys. Fluids* **33** 015104

[10] Krishnamurti R and Howard L N 1981 *P. Natl. Acad. Sci. U.S.A.* **78** 1981-5

[11] Lakkaraju R 2012 *PhD Thesis* University of Twente, The Netherlands

[12] Lessani B and Papalexandris M V 2006 *J. Comput. Phys.* **212** 218-46

[13] Mazzitelli I M, Lohse D and Toschi F 2003 *J. Fluid Mech.* **488** 283-313

[14] Niemela J J, Skrbek L, Sreenivasan K R and Donnelly R J 2000 *Nature* **404** 837-40

[15] Oresta P, Lohse D and Prosperetti A 2009 *Phys. Rev. E* **80** 026304

[16] Pallarès J, Arroyo M, Grau F X and Giralt F 2001 *Exp. Fluids* **31** 208-18

[17] Pallarès J, Cuesta I, Grau F X and Giralt F 1996 *Int. J. Heat Mass Tran.* **39** 3233-47

[18] Paolucci S 1982 *NASA STI/Recon Technical Report N* **83** 26036

[19] Papalexandris M V 2020 *Continuum Mech. Therm.* **32** 1245-9

[20] Plesset M S and Zwick S A 1954 *J. Appl. Phys.* **25** 493-500

[21] Shams E 2009 *PhD Thesis* Oregon State University, USA

[22] Shams E, Finn J and Apte S V 2011 *Int. J. Numer. Meth. Fl.* **67** 1865-98

[23] Taylor G I and Green A E 1937 *P. Roy. Soc. Lond. A Mat.* **158** 499-521

[24] Zhou Q and Xia K 2013 *J. Fluid Mech.* **721** 199-224