1 Introduction
For industrial applications on complex two-phase flows, interfacial forces have to be characterised and modelled. A particular concern of the nuclear industry is the presence of bubbles at the wall where boiling characteristics strongly depend on the distribution and trajectories of bubbles. Hence, the prediction of void fraction, which relies on interfacial force models, is a major issue. In two-phase flows, limitations are related to our lack of understanding of local phenomena such as the bubble wake structures. A powerful tool which may be used to provide further local information is direct numerical simulation (DNS) used as a form of ‘numerical experiment’. An up-scaling process (Bois Reference Bois2017; du Cluzeau, Bois & Toutant Reference du Cluzeau, Bois and Toutant2019) allows the study of two-phase flows towards an improvement of turbulence and interfacial force modelling. This process relies on our capacity to link the exact local formulation of the Navier–Stokes equations in DNS to several notions such as the turbulent dispersion or lift forces which come from the consideration of bubbles as point particles. Splitting interfacial forces as in a particle approach (into drag, lift, added mass etc.) is not straightforward from the local and continuous standpoints, where only the total interfacial force is defined as an integral of the constraint over the interface. It requires an assumption on the topology of the flow (small, point-size bubbles). Besides, even though lift and drag forces can be separated by projecting the integral in the cross-flow and streamwise components, no local definition allows us to separate forces acting in the same direction such as lift and turbulent dispersion forces or drag and added-mass forces.
One of the main concerns of the two-fluid model is to define the content of the momentum transfer between phases (related to the interfacial forces on bubbles). Some properties of this transfer have been demonstrated by Geurst (Reference Geurst1986) and Wallis (Reference Wallis1990). Since their works, we know that the momentum transfer is comprised of a term of added mass and of combined effects of the relative velocity between phases and mean shears (lift force, drag force). In the simplified case of spherical bubbles, Drew & Passman (Reference Drew and Passman1999) propose a fully defined decomposition of these forces. For complex industrial cases, these works thus give the form of the closures for the momentum transfer but not their actual definitions, because these closures must be completed by unknown coefficients (lift, drag, added-mass coefficients etc.). Based solely on these works, it is for instance impossible to measure forces experimentally or numerically in complex configurations, since there is no all-regime definition of the coefficients.
These coefficients have been extensively studied over the past 20 years. Experiments have first brought to light empirical relations for their closures. The lift (Tomiyama et al. Reference Tomiyama, Tamai, Zun and Hosokawa2002) and the drag (Ishii & Zuber Reference Ishii and Zuber1979) forces have been particularly studied for all kind of flow regimes (laminar, turbulent, transitional; for spherical, ellipsoidal or cap bubbles). Despite all these studies, the reversal of the lift force for instance is still poorly modelled. Historically, these coefficients have been highlighted in simplified cases. For instance, the lift force can be studied because it appears alone in the transverse direction in laminar shear flows (Legendre & Magnaudet Reference Legendre and Magnaudet1998; Tomiyama et al. Reference Tomiyama, Tamai, Zun and Hosokawa2002). Then, the compilation of elementary forces gives a trajectory equation of the bubbles (see for instance Legendre, Borée & Magnaudet (Reference Legendre, Borée and Magnaudet1998), Magnaudet & Eames (Reference Magnaudet and Eames2002)) which gives very good predictions on a lot of studies and industrial applications. With the growing complexity of the studied flows, the scientific community tries nowadays to characterise increasingly complex phenomena. In such configurations, intricate interactions arise between turbulence, geometry and nearby bubbles. With this train of thought, the present work proposes a study of the interfacial force balance at the statistical steady state of a rather complicated flow (turbulent bubbly channel flow). In this case, one of the main challenges is that the averaged interfacial force acting on a bubble at statistical steady state is zero. Because of this, as long as interfacial forces (lift, drag etc.) do not rely on local definitions, statistical analysis will not help in understanding the complex equilibrium of forces. Our first objective is then to use the exact local and continuous equations to propose new definitions for the interfacial forces which reflect the meaning usually assigned to the forces in the particle approach (Geurst Reference Geurst1986; Wallis Reference Wallis1990; Drew & Passman Reference Drew and Passman1999).
A second issue raised in this work is to verify the relevance of the particle hypothesis (assimilation of bubbles to point particles) which is generally used in the two-fluid model. Indeed, to the best of our knowledge, strong hypotheses on surface tension and vapour pressure are always made in Reynolds averaged Navier–Stokes (RANS) applications since the work of Prosperetti & Jones (Reference Prosperetti and Jones1984) later taken up by Ishii & Hibiki (Reference Ishii and Hibiki2006). Indeed, they demonstrate that the impact of the vapour pressure is systematically compensated by the action of the surface tension under several hypotheses. In this article, we highlight cases in which these hypotheses are obsolete. Indeed, all theoretical works on the two-fluid model consider, more or less directly, the size of the bubbles as negligible. For instance, the result of Prosperetti & Jones (Reference Prosperetti and Jones1984) requires a generalised Laplace law for which the local pressure jump at the interface is extended to the entire vapour phase. This assumption seems correct for small bubbles, but what about larger ones? For larger bubbles, the averaged momentum transfer term is not directly related to the interfacial force acting on bubbles but to the mean force acting on a fluid element of the vapour phase. Local effects of surface tension and pressure can then arise, as in the work of Moraga et al. (Reference Moraga, Larreteguy, Drew and Lahey2006) in which a desire to go further than the particle hypothesis is present. Indeed, to correct the wall lubrication force, Moraga et al. (Reference Moraga, Larreteguy, Drew and Lahey2006) need to redefine the phase indicator because of the use of the particle hypothesis in a region where it is not valid (in the near-wall region). Their work thus supports the idea that the force acting on each bubble as an entity is not sufficient to describe all the dynamics of the flow. In our work, we are wondering whether this kind of approach (see also Lubchenko et al. (Reference Lubchenko, Magolan, Sugrue and Baglietto2017)) can be generalised and whether the effects of surface tension can be found in other areas of the flow.
A possible application of our proposal concerns the appropriate closure of the turbulent dispersion force for which a lot of different models coexist (Lopez de Bertodano Reference Lopez de Bertodano1998; Chahed, Roig & Masbernat Reference Chahed, Roig and Masbernat2003; Laviéville et al. Reference Laviéville, Mérigoux, Guingo, Baudry and Mimouni2015). The abundance of closure relations for the dispersion force suggests that hidden variables play a role in that phenomenon and that turbulent fluctuations are not alone in driving the bubble spreading. We will see that several additional forces are commonly neglected because of the particle hypothesis. We aim to study those forces to complete the modelling of the dispersion force. In a previous work (du Cluzeau et al. Reference du Cluzeau, Bois and Toutant2019), we have demonstrated that different levels of turbulence can be achieved for the same bubble dispersion. This effect has been linked to surface tension forces which may be responsible for a laminar dispersion force (see du Cluzeau et al. Reference du Cluzeau, Bois and Toutant2019). This study suggests that a non-turbulent dispersion force (comprised of pressure and surface tension terms) has an essential impact on the spatial distribution of bubbles. Indeed, in purely potential flows, the numerical work by Smereka (Reference Smereka1993) and the theoretical development by Biesheuvel & Wijngaarden (Reference Biesheuvel and Wijngaarden1984) have demonstrated that spherical bubbles tend to create horizontal alignment due to the modification of the pressure field by the potential flow around the bubbles. Otherwise, deformable bubbles have a naturally oscillating trajectory which can be seen as a dispersion effect. Those phenomena are not associated with standard turbulence. They may explain the lack of reliability of model prediction and the difficulties encountered with most closures for dispersion forces which only take into account the turbulent contribution.
In § 2, the DNSs are briefly described and some general remarks on bubble distribution are made. Then, in § 3, the complete content of the momentum transfer term is defined from the Navier–Stokes equations without relying on the particle hypothesis. This content is then connected to the trajectory equation of a fluid element of the vapour phase to make the link with the forces acting on a bubble in the particle approach. The resulting equation of forces is profusely discussed and forces are identified. In particular, a new laminar dispersion force is highlighted. In § 4, this new formulation is used to study the different forces from DNS data, in particular, the lateral forces which are predominant in the void-fraction distribution in turbulent channel flows (lift, dispersion forces). The laminar dispersion force is found to be particularly important for the spatial distribution of bubbles. In § 5, we apply our approach to the two-fluid model to identify possible improvements of the model for future Euler–Euler RANS applications.
2 Direct numerical simulations
2.1 Numerical set-up
Physical behaviour of two-phase bubbly flows is strongly related to local phenomena occurring in the surroundings of the bubbles (an increase of the dissipation close to the interface, boundary layer detachment, wake’s structure and interactions). Classical experiments struggle to clearly isolate the individual role of each of these phenomena. DNS used as a numerical experiment is an alternative. In the present work, five calculations are studied for different physical conditions. The computational domain is a
$2h\unicode[STIX]{x03C0}\times 2h\times h\unicode[STIX]{x03C0}$
channel at
$Re_{\unicode[STIX]{x1D70F}}=180$
(where
$Re_{\unicode[STIX]{x1D70F}}=hu_{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D708}_{l}$
) except for the case D127, which is twice smaller in every direction (see table 1) between two vertical walls perpendicular to the direction 2 (
$y$
). The distance between the walls is
$2h=2m$
where
$h$
is the characteristic length of the channel. There are periodic conditions in the directions 1 (
$x$
, streamwise) and 3 (
$z$
, spanwise). The flow is driven upward by a mean pressure gradient calculated to satisfy an imposed mean velocity gradient at the wall, while the acceleration due to gravity acts downwards (along the
$x$
-axis). In all simulated cases, a population of bubbles (
$d_{b}=0.3h$
) is added in the flow with a void fraction of 3 %.
Table 1. Numerical and physical parameters for calculations at
$Re_{\unicode[STIX]{x1D70F}}=127$
and
$Re_{\unicode[STIX]{x1D70F}}=180$
for
$h=1~\text{m}$
(the characteristic length scale of the channel). The bubble diameter is
$d_{b}$
,
$N_{b}$
is the number of bubbles,
$\unicode[STIX]{x0394}x^{+}=(L_{x}/N_{x})Re_{\unicode[STIX]{x1D70F}}/h$
, where
$L_{x}$
and
$N_{x}$
are respectively the length and the number of cells in the
$x$
direction. The Eötvös number is
$Eo=\unicode[STIX]{x1D70C}_{l}gd_{b}^{2}/\unicode[STIX]{x1D70E}$
. The parietal Reynolds number is defined as
$Re_{\unicode[STIX]{x1D70F}}=hu_{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D708}_{l}$
, with
$u_{\unicode[STIX]{x1D70F}}=\sqrt{(\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D70C}_{l})(\unicode[STIX]{x2202}\overline{u}/\unicode[STIX]{x2202}y|_{y=0})}$
;
$Re_{c}=D_{h}\langle u\rangle /\unicode[STIX]{x1D708}_{l}$
is the channel Reynolds number, with the hydraulic diameter
$D_{h}=4h$
;
$Re_{b}=d_{b}\langle u_{r}\rangle /\unicode[STIX]{x1D708}_{l}$
is the bubble Reynolds number, with
$u_{r}=\overline{u_{v}}^{v}-\overline{u_{l}}^{l}$
, the mean upstream relative velocity of the bubbles;
$\unicode[STIX]{x0394}t_{ave}$
,
$\unicode[STIX]{x0394}t_{ave}^{\unicode[STIX]{x1D70F}}$
are the time intervals on which statistical results have been measured. They are expressed respectively in the time unit
$t.u.=h/\langle u\rangle$
,
$t.u.\unicode[STIX]{x1D70F}=h/u_{\unicode[STIX]{x1D70F}}$
.
TrioCFD, through its front-tracking algorithm (Mathieu Reference Mathieu2003), resolves the one-fluid equations of Kataoka (Reference Kataoka1986) as written for channel up-flow by Lu & Tryggvason (Reference Lu and Tryggvason2008),
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}\boldsymbol{u}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{u}\boldsymbol{u}) & = & \displaystyle -\unicode[STIX]{x1D735}P^{\prime }+(\unicode[STIX]{x1D70C}-\langle \unicode[STIX]{x1D70C}\rangle )\boldsymbol{g}-\unicode[STIX]{x1D6FD}\boldsymbol{e}_{x}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}^{\text{T}}\boldsymbol{u})]+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\boldsymbol{n}_{v}\unicode[STIX]{x1D6FF}^{i},\end{eqnarray}$$
where each of the one-fluid variables is defined as a mixture of phase variables:
$\unicode[STIX]{x1D719}=\sum _{k}\unicode[STIX]{x1D712}_{k}\unicode[STIX]{x1D719}_{k}$
where
$\unicode[STIX]{x1D719}_{k}$
can be
$\boldsymbol{u}_{k}$
,
$\unicode[STIX]{x1D70C}_{k}$
,
$\unicode[STIX]{x1D707}_{k}$
or
$P_{k}$
, respectively the velocity, the density, the viscosity and the pressure in phase
$k$
;
$\unicode[STIX]{x1D712}_{k}$
is the phase indicator function, which is equal to 1 in phase
$k$
and 0 otherwise. The transport equation is then
$\unicode[STIX]{x2202}\unicode[STIX]{x1D712}_{v}/\unicode[STIX]{x2202}t+\boldsymbol{u}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}=0$
, where
$\boldsymbol{u}_{i}$
is the interfacial velocity of the front. Here, in the absence of phase change, the velocity is continuous across the interface and
$\boldsymbol{u}_{i}=\boldsymbol{u}$
. Following the proposal of Lu & Tryggvason (Reference Lu and Tryggvason2008), the pressure gradient
$\unicode[STIX]{x1D735}P$
is split into mean
$\langle \unicode[STIX]{x1D735}P\rangle$
and fluctuating
$\unicode[STIX]{x1D735}P^{\prime }$
parts (
$\langle \unicode[STIX]{x1D719}\rangle$
is the space average of
$\unicode[STIX]{x1D719}$
over the whole domain). The parameter
$\unicode[STIX]{x1D6FD}$
is a constant source term containing the spatially averaged weight of the mixture and the driving pressure gradient so that
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D735}\langle P\rangle -\langle \unicode[STIX]{x1D70C}\rangle \boldsymbol{g}$
. It corresponds to an imposed shear stress at the wall. The parameter
$\unicode[STIX]{x1D70E}$
is the surface tension,
$\unicode[STIX]{x1D705}=-\unicode[STIX]{x1D735}_{S}\boldsymbol{\cdot }\boldsymbol{n}_{v}$
is the local curvature, usually negative for bubbles, and
$\boldsymbol{n}_{v}$
is the interface normal defined by
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}=-\boldsymbol{n}_{v}\unicode[STIX]{x1D6FF}^{i}$
, where
$\unicode[STIX]{x1D6FF}^{i}$
is the Dirac impulse at the interface
$i$
.
Figure 1. Instantaneous illustrations of the flows. Interfaces are coloured to indicate the local curvature. The colour scale represents isovalues of the
$\unicode[STIX]{x1D706}_{2}$
criterion (Jeong & Hussain Reference Jeong and Hussain1995). Instantaneous illustration of the case D127 is provided in du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019).
Following the proposal of Tryggvason et al. (Reference Tryggvason, Bunner, Esmaeeli and Al-Rawahi2003), a front-tracking method is used to solve this set of equations in the whole computational domain, including both the vapour and liquid phases. The interface is followed by moving connected marker points. The Lagrangian markers are advected by the velocity field interpolated from the Eulerian grid. In order to preserve the mesh quality and to limit the need for remeshing operations, only the normal component of the velocity field is used in the marker transport. After transport, the front is used to update the phase indicator function, the density and the viscosity at each Eulerian grid point. The Navier–Stokes equations are then solved by a projection method (Puckett et al. Reference Puckett, Almgren, Bell, Marcus and Rider1997) using fourth-order central differentiation for evaluation of the convective and diffusive terms on a fixed, staggered Cartesian grid. Fractional time stepping leads to a third-order Runge–Kutta scheme (Williamson Reference Williamson1980) (see Toutant Reference Toutant2006). In the two-step prediction–correction algorithm, a surface tension source is added to the main flow source term and to the evaluation of the convection and diffusion operators in order to obtain the predicted velocity (see Mathieu (Reference Mathieu2003) for further information). Then, an elliptic pressure equation is solved by an algebraic multigrid method to impose a divergence-free velocity field. TrioCFD has already been widely used for two-phase (Toutant et al. Reference Toutant, Labourasse, Lebaigue and Simonin2008, Reference Toutant, Chandesris, Jamet and Lebaigue2009; Toutant, Mathieu & Lebaigue Reference Toutant, Mathieu and Lebaigue2012; Bois, Fauchet & Toutant Reference Bois, Fauchet and Toutant2016; Bois Reference Bois2017; Bois et al. Reference Bois, du Cluzeau, Toutant and Martinez2017) and single-phase (Chandesris & Jamet Reference Chandesris and Jamet2006, Reference Chandesris and Jamet2009; Chandesris et al. Reference Chandesris, D’Hueppe, Mathieu, Jamet and Goyeau2013; Dupuy, Toutant & Bataille Reference Dupuy, Toutant and Bataille2018) flow studies.
Then, averaged information is extracted from the DNS data. Because of the low void fraction in our simulations, break-up and coalescence are neglected. Since the present focus is on dynamical aspects (interfacial forces and bubble-induced turbulence), the channel is considered isothermal with no phase change or boiling at the wall and constant properties within each phase. Physical and numerical set-ups are summarised in table 1. In order to study statistical profiles, the ensemble averaging has been assimilated to a temporal averaging, particularised to a space and time average by application of the ergodicity hypothesis to the periodic directions of the flow,
$$\begin{eqnarray}\overline{\unicode[STIX]{x1D719}(y)}=\frac{1}{\unicode[STIX]{x0394}t_{ave}L_{x}L_{z}}\int _{t-\unicode[STIX]{x0394}t_{ave}/2}^{t+\unicode[STIX]{x0394}t_{ave}/2}\int _{0}^{L_{x}}\int _{0}^{L_{z}}\unicode[STIX]{x1D719}(x,y,z,\unicode[STIX]{x1D70F})\,\text{d}x\,\text{d}z\,\text{d}\unicode[STIX]{x1D70F},\end{eqnarray}$$
where
$L_{x}$
and
$L_{z}$
are respectively the length and depth of the channel and
$\unicode[STIX]{x0394}t_{ave}$
is the time interval of the average expressed in the time unit
$t.u.=h/\langle u\rangle$
. Validations of the code are provided in du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019).
2.2 General remarks
In figure 1, instantaneous illustrations of the flows are shown for the 4 cases at
$Re_{\unicode[STIX]{x1D70F}}=180$
. Cases S180 and D180 are channel bubbly up-flows in low gravity conditions (
$Re_{b}\in [90,140]$
), respectively with spherical (
$Eo=0.45$
) and deformable (
$Eo=3.6$
) bubbles. Cases S180g8 and D180g8 are similar cases with increased gravity, representative of a more realistic range of bubble Reynolds number (
$Re_{b}\in [450,600]$
) reflecting standard air/water conditions for 2 mm bubbles as observed in reference experiments (Riboux, Legendre & Risso Reference Riboux, Legendre and Risso2013). The fifth case introduced in du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019), is similar to case D180 but at a lower Reynolds number
$Re_{\unicode[STIX]{x1D70F}}=127$
. These cases have been especially designed for the study of the transversal migration of bubbles (i.e. the void fraction profiles). As can be seen from figure 1 and more precisely from figure 2(a), all the typical void fraction profiles are represented. The void-fraction profiles have two important features. The first one is the position of the void-fraction peak related to the action of the lift force on the bubbles. The second is its spreading due to the dispersion forces. By spreading, we mean the thickness of the transitional region between the void-fraction peak and the low void-fraction area.
Figure 2. Statistical profiles of void fraction and Weber number from DNS. (a) Void fraction. (b) Weber number.
2.2.1 Void-fraction peak
The strongest effect of an increasing Eötvös number is the reversal of the lift force starting from a critical value
$Eo_{c}(Re_{\unicode[STIX]{x1D70F}}=127)\approx 2.5$
(Tryggvason et al.
Reference Tryggvason, Dabiri, Aboulhasanzadeh and Lu2013). Because of the spherical shape of the bubbles in case S180 (
$Eo=0.45<Eo_{c}$
), the bubbles are pushed against the wall under the action of the lift force to reach a typical wall-peaked profile. Similarly, case S180g8 shows a wall-peaked profile but it also has the strongest diffusion of bubbles to the bulk. Actually, the diffusion is strong enough to strongly decrease the maximum of void fraction and almost reduce it to the constant value obtained in the core of the flow (
$\unicode[STIX]{x1D6FC}\backsimeq 3\,\%$
for
$y/h>0.4$
). For the deformable cases D127 and D180g8 for which
$Eo=4.5-3.6>Eo_{c}$
, the lift force is reversed. Hence, the bubbles are pushed towards the bulk of the channel, and the flow reaches a classical core-peaked profile. However, the intermediate behaviour in case D180 needs further explanation. The intermediate peak of the void-fraction profile in case D180 suggests that the
$Eo$
number is not sufficient to characterise the orientation of the lift force. Indeed, the Eötvös number is independent of the location whereas bubbles located at
$y/h<0.3$
are pushed to the right (away from the wall) when bubbles located in
$y/h>0.3$
are pushed to the left. The Eötvös number does not take into account the bubble relative velocity which may have an impact on the lift orientation (Adoua, Legendre & Magnaudet Reference Adoua, Legendre and Magnaudet2009). In this context, the Weber number, defined as
$We=\unicode[STIX]{x1D70C}_{l}u_{r}^{2}d_{b}/\unicode[STIX]{x1D70E}$
, seems more appropriate. The Weber profile is plotted in figure 2(b). The value of the critical Weber number
$We_{c}=3$
is found at
$y/h=0.3$
for case D180. This new criterion can explain the intermediate peak of case D180 and agrees with the core-peaked profile of cases D127 and D180g8 for which
$We>We_{c}$
and with the wall-peaked profile of cases S180 and S180g8 for which
$We<We_{c}$
. It clearly separates positive and negative lift forces for all simulated cases.
2.2.2 Void-fraction spreading
The void-fraction spreading arises from a balance between the magnitudes of the lift force and the dispersion forces. For RANS applications, the dispersion forces are often considered as turbulent dispersion forces only. Indeed, the mixing induced by turbulence is an essential source of dispersion. Nevertheless, the present results suggest the presence of hidden variables responsible for other non-turbulent dispersion. In low gravity conditions (i.e. small bubble Reynolds numbers), the averaged bubble wakes are stable for cases S180 and D180 so that turbulent structures (shown in figure 1) are mostly related to the turbulence induced near the wall by the local shear (see the elongated streaks on case D180). The flow rate in case S180 is strongly decreased because of the presence of bubbles against the wall in order to satisfy the criterion
$Re_{\unicode[STIX]{x1D70F}}=180$
. This effect explains the decrease, from 12 100 to
$5400$
, of the channel Reynolds number between cases D180 and S180 and hence the small number of turbulent structures in case S180 (see Dabiri, Lu & Tryggvason (Reference Dabiri, Lu and Tryggvason2017), Lu & Tryggvason (Reference Lu and Tryggvason2008) or du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019) for additional explanation on the flow rate reduction). The small wall-normal Reynolds stresses in case S180 (see figure 3
b) may explain why bubbles are stuck to the wall. The slight increase of the wall-normal Reynolds stresses in case D180 and D127 in comparison with case S180 (because turbulence freely develops in the boundary layer in the absence of bubble) allows for a spreading of the void fraction.
Figure 3. Statistical profiles of velocities and Reynolds stresses from DNS. (a) Cross-correlation, (b) wall-normal Reynolds stresses, (c) liquid velocity, (d) relative velocity.
A small difference of turbulent kinetic energy between cases S180 and D127 induces a significant difference of spreading. This fact can be explained by a modification of the lift force magnitude, which is possibly stronger in case S180 due to a different wake structure. In normal conditions of gravity (i.e. for cases S180g8 and D180g8), turbulent structures related to the interaction and instabilities of the wakes are essential. This effect, combined with an increase of bubble relative velocity, induces stronger Reynolds stresses in cases D180g8 and S180g8 (see figure 3
b). Under the action of the turbulent dispersion force, this increase of the wall-normal Reynolds stresses could explain the substantial increase of spreading between cases S180 and S180g8. However, even if the wall-normal Reynolds stress is six times larger for case D180g8 than for case D127, figure 2(a) shows that both cases have similar void fraction profiles. This observation suggests that the role of turbulence in the spreading process is overvalued and not alone. The Reynolds stresses magnitudes follow the same hierarchy for other components (not shown here). Furthermore, the spreading in case S180g8 may be explained by the strong kinetic energy of bubbles which bounce violently on the walls and are ejected towards the bulk. This effect is visible in a movie of the flow (see supplementary material and statistical data at http://triocfd.cea.fr/Pages/Research_directions/FT_database.aspx). The overrating of the role of turbulence in the dispersion process is in agreement with our previous findings in du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019) and can be explained physically. For instance, in purely laminar flows, Biesheuvel & Wijngaarden (Reference Biesheuvel and Wijngaarden1984) demonstrates that spherical bubbles tend to create horizontal alignment due to the modification of the pressure field by the potential flow around the bubbles. Bunner & Tryggvason (Reference Bunner and Tryggvason2002) show numerically that binary interaction between two bubbles induces an horizontal alignment of bubbles. Furthermore, for bubbly flows with low void fraction (as in our cases
${<}5\,\%$
), Batchelor (Reference Batchelor1972) and Wijngaarden & Kapteyn (Reference Wijngaarden and Kapteyn1990) assume that the flow dynamics is dominated by those binary interactions. Otherwise, deformable bubbles have a naturally oscillating trajectory which can be seen as a dispersion mechanism. To the best of our knowledge, these effects are never taken into account in RANS models. Finally, in our previous work (du Cluzeau et al.
Reference du Cluzeau, Bois and Toutant2019), we showed that an impact of the surface tension on dispersion forces is possible.
3 Force identification method
In this section, a new method of identifying forces is proposed. In § 3.1, the complete content of the momentum transfer term is defined from the Navier–Stokes equations. This definition differs from the classical one because it does not rely on the particle hypothesis. As mentioned in § 2.2, a non-turbulent dispersion force can come from pressure or surface tension effects. Hence, the derivation of the present formulation does not resort to any hypothesis for the surface tension force or the pressure field. Section 3.2 links the momentum transfer definition to the trajectory equation of a fluid element of the vapour phase to make the connection with the forces acting on a bubble in the particle approach. The resulting trajectory equation is then studied in § 3.3 for an ideal laminar shear flow to demonstrate the link between the new local definition of the lift force and its classical closure in a particle approach. These analytical developments show that the new definitions can coincide with classical definitions by integrating the volume of the bubbles. The trajectory equation is then profusely discussed in § 3.4 and some forces are identified. In particular, a new laminar dispersion force is highlighted.
3.1 Evolution equation of momentum transfer
In this section, the averaged formulation of the Navier–Stokes equations is used to find a definition for the momentum transfer in a general case (without any hypotheses). Let us first consider the averaged Navier–Stokes equations in each phase
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}_{k}={\displaystyle \frac{\text{D}(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k})}{\text{D}t}}+\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{k}(\overline{p_{k}}^{k}-\widetilde{p_{l}})]-\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\boldsymbol{g}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FC}_{k}\overline{\unicode[STIX]{x1D70F}_{k}}^{k}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}^{\prime }\boldsymbol{u}_{k}^{\prime }}^{k}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{with}~k\in [l,v]\quad \text{and}\quad \boldsymbol{M}_{l}+\boldsymbol{M}_{v}=\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}},\quad \boldsymbol{M}_{k}=\overline{p_{k}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{k}-\unicode[STIX]{x1D749}_{k}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{k}}, & \displaystyle\end{eqnarray}$$
$\text{D}(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k})/\text{D}t=\unicode[STIX]{x2202}(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k})/\unicode[STIX]{x2202}t+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k}\overline{\boldsymbol{u}_{k}}^{k})=\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\text{d}\overline{\boldsymbol{u}_{k}}^{k}/\text{d}t$
with
$\text{d}/\text{d}t$
the material derivative (see § A.1). Indices
$l$
and
$v$
are respectively related to the liquid and vapour phases. The velocity
$\boldsymbol{u}_{k}$
is decomposed into mean and fluctuating parts
$\boldsymbol{u}_{k}=\overline{\boldsymbol{u}_{k}}^{k}+\boldsymbol{u}_{k}^{\prime }$
and
$\overline{\unicode[STIX]{x1D719}_{k}}^{k}=\overline{\unicode[STIX]{x1D712}_{k}\unicode[STIX]{x1D719}_{k}/}\overline{\unicode[STIX]{x1D712}_{k}}$
is the phase average. The pressure
$p$
is defined relatively to a reference pressure
$\widetilde{p_{l}}$
;
$\overline{\unicode[STIX]{x1D70F}_{k}}^{k}$
is the mean viscous stress. Then, considering (3.1b
) and 
, a more detailed content of the interfacial momentum transfer is found, 
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{v}=\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\boldsymbol{M}_{l} & = & \displaystyle \unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D6FC}_{l}(\unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{v})\boldsymbol{g}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}_{v}{\displaystyle \frac{\text{D}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}}{\text{D}t}}+\unicode[STIX]{x1D6FC}_{l}{\displaystyle \frac{\text{D}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}}{\text{D}t}}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FC}_{l}\overline{\unicode[STIX]{x1D70F}_{l}}^{l}-\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70F}_{v}}^{v}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l})+\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{l}(\overline{p_{l}}^{l}-\widetilde{p_{l}})]+\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\widetilde{p_{l}})].\end{eqnarray}$$
As explained in § A.1,
$\boldsymbol{M}_{v}$
is an interfacial momentum transfer term and not the sum of external forces on the vapour phase. Therefore it has a non-zero mean value (
$\boldsymbol{M}_{v}\neq 0$
) even at statistical steady state. One of the challenges of the two-fluid approach is the modelling of this momentum transfer. As we saw in the Introduction, it is classically composed of added-mass, drag and lift forces (Geurst Reference Geurst1986; Wallis Reference Wallis1990; Drew & Passman Reference Drew and Passman1999). As it stands, classical expressions for these forces are not found in (3.2) because they are in a different form, written with local variables. However, these classical terms are inevitably there, since (3.2) is exact. To determine if the classical terms can be identified and separated, we use in the following the trajectory equation of a fluid element of the vapour phase, in order to have all the forces and to make the link with the trajectory equation of a particle.
3.2 Trajectory equation of a fluid element of the vapour phase
In order to understand the equilibrium of forces which leads to bubble migration, equation (3.1a
) is used again. By replacing
$\boldsymbol{M}_{v}$
in (3.2) by its expression in (3.1a
), an expression for the total sum of forces on the vapour phase
$\boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}$
is reached,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\frac{\text{d}\overline{\boldsymbol{u}_{v}}^{v}}{\text{d}t} & = & \displaystyle \boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70C}}\boldsymbol{g}-\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}}^{l}-\tilde{p_{l}})\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}_{v}{\displaystyle \frac{\text{D}(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l})}{\text{D}t}}+\unicode[STIX]{x1D6FC}_{l}{\displaystyle \frac{\text{D}(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v})}{\text{D}t}}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\overline{\unicode[STIX]{x1D70F}_{l}}^{l})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70F}_{v}}^{v})]\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v})]\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})].\end{eqnarray}$$
For additional clarity, the averaged weight of the mixture included in the pressure gradient (reflecting the hydrostatic equilibrium in the longitudinal direction) can be extracted,
$$\begin{eqnarray}\displaystyle \overline{p_{l}}^{l}=\overline{p_{l}^{\prime }}^{l}-\langle \unicode[STIX]{x1D70C}\rangle \boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{x}. & & \displaystyle\end{eqnarray}$$
To understand the physical meaning of the different terms, divergence terms are split in order to discriminate dispersion forces (
$\propto \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}$
) from other forces:
For the purpose of giving a local definition to forces, it has to be noticed that this equation is not the trajectory equation of a point particle (i.e. Newton’s second law in a particle approach) but the trajectory equation of an elementary volume of fluid in the vapour phase. Then, the right-hand side of (3.6) is the total force acting on the elementary volume and not on the complete bubble. At the statistical steady state, the flow reaches a spatial distribution of bubbles (i.e. void fraction and relative velocity profiles) for which there is an equilibrium between all the forces of the equation. In (3.6), several contributions are clearly identified;
$\boldsymbol{M}^{\unicode[STIX]{x1D72B}}$
is the buoyancy force resulting from a density difference compensated by the drag force. The turbulent dispersion force
$\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}$
(proportional to
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}$
and to the Reynolds stresses) also arises in (3.6). As will be discussed in § 4.2, the term
$\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}$
, coined as turbulent dispersion, is not a purely turbulent contribution. Other forces
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}=\boldsymbol{M}^{\unicode[STIX]{x1D748}}+\boldsymbol{M}^{\boldsymbol{P}}\propto \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}$
are coined as ‘laminar dispersion forces’. The viscous dispersion and viscous force are named
$\boldsymbol{M}^{\unicode[STIX]{x1D749}\boldsymbol{D}}$
and
$\boldsymbol{M}^{\unicode[STIX]{x1D749}}$
. In the following, assessments of
$\boldsymbol{M}^{\unicode[STIX]{x1D749}\boldsymbol{D}}$
and
$\boldsymbol{M}^{\unicode[STIX]{x1D749}}$
from DNS show that they are negligible in comparison with the other forces. Therefore, they will not be discussed in the subsequent analysis. In order to clarify the role of the pressure gradients, of the instationary terms and of the last term of (3.6), the simplified case of an isolated bubble in a laminar shear flow is considered.
3.3 Lift force on an isolated bubble in a laminar shear flow
Figure 4. Scheme of a single quasi-spherical bubble in a shear flow.
In (3.6), terms related to the dispersion process in the lateral direction are proportional to
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}$
. The lift force, commonly expressed in terms of
$\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D6FC}_{v}u_{r}\unicode[STIX]{x2202}\overline{u_{l}}^{l}/\unicode[STIX]{x2202}y$
in the particle approach, does not appear immediately in (3.6). In order to identify it, let us consider an isolated bubble at low bubble Reynolds number in a laminar shear flow. In these conditions, only the lift force
$M^{L}$
leads the lateral trajectory of the bubble (see Tomiyama et al.
Reference Tomiyama, Tamai, Zun and Hosokawa2002),
$$\begin{eqnarray}\boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}\boldsymbol{\cdot }\boldsymbol{e}_{y}=M^{L}.\end{eqnarray}$$
Let us consider an isolated bubble with
$\unicode[STIX]{x1D6FA}(t)$
the domain in which the void fraction is non-zero at a given time
$t$
(see figure 4). For the statistics of a channel flow, the statistical data defined as in (2.3) are relevant because each part of the finite-size bubble (right, left, central parts) passes through the plane
$y/h=Cte$
if a sufficient averaging time is chosen (for
$y>d_{b}/2$
). Then, the statistical average is representative of the whole bubble and statistical quantities can be regarded as particle quantities. For the present derivation, at a given time
$t$
, the statistical average representative of the point-size bubble is
$$\begin{eqnarray}\displaystyle \overline{\overline{\unicode[STIX]{x1D719}(t)}}=\frac{1}{d_{b}L_{x}L_{z}}\int _{\unicode[STIX]{x1D6FA}(t)}\unicode[STIX]{x1D719}\,\text{d}x\,\text{d}y\,\text{d}z & & \displaystyle\end{eqnarray}$$
and the phase averages are thus redefined as
$\overline{\unicode[STIX]{x1D719}_{k}}^{k}=\overline{\overline{\unicode[STIX]{x1D712}_{k}\unicode[STIX]{x1D719}_{k}}}/\overline{\overline{\unicode[STIX]{x1D712}_{k}}}$
(only for the present section). In this configuration, the dispersion forces in (3.6) are zero because a single quasi-spherical bubble cannot be dispersed (there is no turbulent dispersion, oscillating path or swarm effect). The surface tension source term is zero by definition because it is integrated on the whole bubble. The viscous forces are neglected after observations of DNS statistics (see figure 5). Then (3.6) is reduced in the transversal direction to,
$$\begin{eqnarray}\displaystyle M^{L} & {\approx} & \displaystyle -\frac{\unicode[STIX]{x1D6FC}_{v}}{d_{b}}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\overline{p_{l}^{\prime }}^{l}-\tilde{p_{l}})\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle -\,\frac{\unicode[STIX]{x1D6FC}_{v}}{d_{b}}\left[\unicode[STIX]{x1D6FC}_{l}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\unicode[STIX]{x1D70C}_{l}\overline{v_{l}^{\prime }v_{l}^{\prime }}^{l}\,\text{d}y)+\unicode[STIX]{x1D6FC}_{v}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\unicode[STIX]{x1D70C}_{v}\overline{v_{v}^{\prime }v_{v}^{\prime }}^{v}\,\text{d}y)\right],\end{eqnarray}$$
where
$v$
is the transversal component of the velocity
$\boldsymbol{u}$
. All terms involving the Reynolds stresses could be misunderstood as having no laminar origin. For instance, equation (3.9) suggests that the lift force is partly proportional to the Reynolds stresses. Nevertheless, lift forces exist even for laminar or potential flows. In fact, this apparent contradiction is explained by the dual nature of fluctuations in two-phase flows. Indeed, a part of the Reynolds stresses in bubbly flows is induced by local inhomogeneity of the velocity field due to the presence of bubbles. In particular, the averaged wakes of the bubbles and the potential flow around them are related to non-turbulent spatial fluctuations, coined as wake-induced fluctuations. Hence, a part of the Reynolds stresses is non-turbulent (see Riboux et al. (Reference Riboux, Legendre and Risso2013), Amoura, Besnaci & Risso (Reference Amoura, Besnaci and Risso2017), Risso (Reference Risso2018) and du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019) for further details). Following these works, the Reynolds stresses decomposition leads to,
$$\begin{eqnarray}\displaystyle \overline{u_{k}^{\prime }u_{k}^{\prime }}^{k}=\overline{u_{k}^{\prime }u_{k}^{\prime }}^{k}|_{turb}+\overline{u_{k}^{\prime }u_{k}^{\prime }}^{k}|_{lam}. & & \displaystyle\end{eqnarray}$$
Here,
$\overline{u_{k}^{\prime }u_{k}^{\prime }}^{k}|_{turb}$
is comprised of turbulent fluctuations related to shear-induced turbulence or to instabilities and interactions of bubble wakes;
$\overline{u_{k}^{\prime }u_{k}^{\prime }}^{k}|_{lam}$
contains all the non-turbulent spatial fluctuations related to the potential flow around the bubbles and to their averaged wakes. In the present case of laminar shear flow, equation (3.9) becomes,
$$\begin{eqnarray}\displaystyle M^{L} & {\approx} & \displaystyle \underbrace{-\frac{\unicode[STIX]{x1D6FC}_{v}}{d_{b}}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\overline{p_{l}^{\prime }}^{l}-\tilde{p_{l}})\,\text{d}y}_{M_{\unicode[STIX]{x1D735}P}^{L}}\nonumber\\ \displaystyle & & \displaystyle -\,\underbrace{\frac{\unicode[STIX]{x1D6FC}_{v}}{d_{b}}\left[\unicode[STIX]{x1D6FC}_{l}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\unicode[STIX]{x1D70C}_{l}\overline{v_{l}^{\prime }v_{l}^{\prime }}^{l}|_{lam}\,\text{d}y)+\unicode[STIX]{x1D6FC}_{v}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\unicode[STIX]{x1D70C}_{v}\overline{v_{v}^{\prime }v_{v}^{\prime }}^{v}|_{lam}\,\text{d}y)\right]}_{M_{\text{Re}}^{L}}.\end{eqnarray}$$
With this equation, the lift force seems to be comprised of a pressure gradient lift force
$M_{\unicode[STIX]{x1D735}P}^{L}$
and of a Reynolds stress lift force
$M_{\text{Re}}^{L}$
. The following sections demonstrate that both
$M_{\unicode[STIX]{x1D735}P}^{L}$
and
$M_{\text{Re}}^{L}$
can be approximated by a classical lift force model proportional to
$\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{l}u_{r}(\unicode[STIX]{x2202}\overline{u_{l}}^{l}/\unicode[STIX]{x2202}y)$
.
3.3.1 Pressure gradient lift force
Applying Green’s theorem,
$M_{\unicode[STIX]{x1D735}P}^{L}$
becomes proportional to the pressure difference between
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}$
and
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}$
,
$$\begin{eqnarray}\displaystyle M_{\unicode[STIX]{x1D735}P}^{L}=\frac{\unicode[STIX]{x1D6FC}_{v}}{d_{b}}(\overline{p_{l}}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}-\overline{p_{l}}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}). & & \displaystyle\end{eqnarray}$$
Let us consider the averaged pressure in the liquid phase as a sum of a ‘single-phase pressure’
$\overline{p_{l}^{SP}}^{l}$
in the absence of bubbles and of a pressure induced by the bubble through surface tension effects
$\overline{p_{l}^{b}}^{l}$
(this decomposition is close to that proposed by Prosperetti & Jones (Reference Prosperetti and Jones1984)),
$$\begin{eqnarray}\displaystyle \overline{p_{l}}^{l}=\overline{p_{l}^{SP}}^{l}+\overline{p_{l}^{b}}^{l}. & & \displaystyle\end{eqnarray}$$
In a single-phase shear flow, there is no pressure gradient in the transverse direction (
$\unicode[STIX]{x1D735}\overline{p_{l}^{SP}}^{l}=0$
). Hence, the only possible pressure gradient in a shear flow is induced by the bubble itself through the action of surface tension. Indeed, a pressure difference between
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}$
and
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}$
can come from an asymmetrical shape of the bubble induced by the shear. In the absence of shear, the pressure difference between
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}$
and
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}$
is zero by symmetry. Considering a constant pressure inside the bubble, we get,
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{p_{v}}^{v}-\overline{p_{l}^{b}}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}\propto \unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{p_{v}}^{v}-\overline{p_{l}^{b}}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}\propto \unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}. & \displaystyle\end{eqnarray}$$
Then, we consider that
$\unicode[STIX]{x1D705}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}$
and
$\unicode[STIX]{x1D705}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}$
are given by the maximum Gauss curvature of a local osculating three-dimensional ellipsoid which fits the bubble interface at the right and left sides. Considering the properties of an ellipsoid, its maximal curvature can be linked to its aspect ratio
$\unicode[STIX]{x1D6FE}$
, then
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{p_{v}}^{v}-\overline{p_{l}^{b}}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}\propto \unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}=\frac{\unicode[STIX]{x1D70E}}{r_{0}}(\unicode[STIX]{x1D6FE}^{-4}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}+\unicode[STIX]{x1D6FE}^{5}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{p_{v}}^{v}-\overline{p_{l}^{b}}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}\propto \unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}=\frac{\unicode[STIX]{x1D70E}}{r_{0}}(\unicode[STIX]{x1D6FE}^{-4}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}+\unicode[STIX]{x1D6FE}^{5}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}), & \displaystyle\end{eqnarray}$$
where
$r_{0}$
is the radius of the equivalent spherical bubble. Here,
$\unicode[STIX]{x1D6FE}$
is related to the local Weber number (
$We=\unicode[STIX]{x1D70C}_{l}u_{r}^{2}d_{b}/\unicode[STIX]{x1D70E}$
) where the relative velocity is defined as the difference between the local liquid velocity and the mean bubble velocity. A general law for the aspect ratio as a function of the Weber number has been proposed by Moore (Reference Moore1965) but for the purpose of analytical developments, the relation for small Weber numbers (
$We\ll 1$
) is chosen here. A Taylor expansion gives,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}=1+\unicode[STIX]{x1D6FD}We+O(We^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}^{-4}+\unicode[STIX]{x1D6FE}^{5}=2+\unicode[STIX]{x1D6FD}We+O(We^{2}), & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FD}=9/6.4$
. We seek here to characterise the difference in deformation between the two sides of a bubble that is observed for example in the high-shear region of case D180 in figure 1. Since the Weber number characterises the deformation of a bubble, we can naturally link the deformation asymmetry to the variation of relative velocity between the two sides of a bubble in a shear (for a given surface tension and a given diameter). Even if expression (3.18) was obtained for a whole bubble, it characterises the effect of a relative velocity on the deformation. It is therefore assumed that this expression takes into account the impact of shear on the asymmetries of deformation. Equations (3.12), (3.16), (3.17) and (3.19) then give,
$$\begin{eqnarray}\displaystyle M_{\unicode[STIX]{x1D735}P}^{L}\propto \frac{2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70E}}{d_{b}^{2}}(We|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}-We|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}). & & \displaystyle\end{eqnarray}$$
Finally, the definition of the Weber number is used,
$$\begin{eqnarray}\displaystyle M_{\unicode[STIX]{x1D735}P}^{L}\propto 2\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}_{v}(u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}+u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}})\frac{(u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}-u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}})}{d_{b}}. & & \displaystyle\end{eqnarray}$$
The difference of relative velocity between
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}$
and
$\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}$
is equal to the difference of liquid velocity. Then, using
$u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}+u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}=2\overline{u_{r}}$
,
$$\begin{eqnarray}\displaystyle M_{\unicode[STIX]{x1D735}P}^{L}\propto 4\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D6FC}_{v}\overline{u_{r}}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{u_{l}}^{l}}{\unicode[STIX]{x2202}y}}. & & \displaystyle\end{eqnarray}$$
These developments show that the pressure disturbance induced by the bubbles leads to a classical formulation of the lift force closure in the presence of shear. We assume that this reasoning can be extended to other regimes of Weber number (for strongly deformable bubbles for instance) as long as the lift coefficient is adapted. For a more complex flow (channel flow for instance), the pressure gradient induced by bubbles (
$\unicode[STIX]{x1D735}\overline{p_{l}^{b}}^{l}$
) and the classical pressure gradient (
$\unicode[STIX]{x1D735}\overline{p_{l}^{SP}}^{l}$
) cannot be distinguished so that this result is difficult to check in complex configurations. The use of the Weber number instead of the Eötvös number has to be noticed. Indeed, based on the Eötvös number which is independent of the flow, an asymmetrical pressure cannot be obtained as it essentially comes from velocity gradients in conjunction with surface tension effects. The use of the Weber number is mandatory. This is in agreement with the observation made in § 2.2.
3.3.2 Reynolds stress lift force
The second remaining force in the case of a bubble in a shear flow relies on Reynolds stresses,
$$\begin{eqnarray}\displaystyle M_{\text{Re}}^{L}=\frac{\unicode[STIX]{x1D6FC}_{v}}{d_{b}}\left[\unicode[STIX]{x1D6FC}_{l}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\unicode[STIX]{x1D70C}_{l}\overline{v_{l}^{\prime }v_{l}^{\prime }}^{l}|_{lam}\,\text{d}y)+\unicode[STIX]{x1D6FC}_{v}\int _{y\in \unicode[STIX]{x1D6FA}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(\unicode[STIX]{x1D70C}_{v}\overline{v_{v}^{\prime }v_{v}^{\prime }}^{v}|_{lam}\,\text{d}y)\right]. & & \displaystyle\end{eqnarray}$$
Then, because
$\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\ll \unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}$
, the last term is neglected and Green’s theorem is applied,
$$\begin{eqnarray}\displaystyle M_{\text{Re}}^{L}\approx -\frac{\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}}{d_{b}}(\overline{v_{l}^{\prime }v_{l}^{\prime }}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}-\overline{v_{l}^{\prime }v_{l}^{\prime }}^{l}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}}). & & \displaystyle\end{eqnarray}$$
In the cross-flow direction, the laminar Reynolds stresses induced by the bubbles are mainly related to the potential flow around them (see Risso Reference Risso2016; Amoura et al. Reference Amoura, Besnaci and Risso2017). Biesheuvel & Wijngaarden (Reference Biesheuvel and Wijngaarden1984) have calculated the contribution of the potential flow to the Reynolds stresses in the case of spherical bubbles. In the transverse directions, they have found,
$$\begin{eqnarray}\overline{v_{l}^{\prime }v_{l}^{\prime }}^{l}|_{lam}={\textstyle \frac{3}{20}}\unicode[STIX]{x1D6FC}_{v}u_{r}^{2}.\end{eqnarray}$$
Thus, using (3.24) and (3.25) with the definition of the relative velocity given in § 3.3.1, we get,
$$\begin{eqnarray}\displaystyle M_{\text{Re}}^{L} & = & \displaystyle -\frac{3\unicode[STIX]{x1D6FC}_{v}^{2}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}}{20}(u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}+u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}})\frac{(u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{+}}-u_{r}|_{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}^{-}})}{d_{b}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle -\frac{6\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D6FC}_{l}}{20}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{l}\overline{u_{r}}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{u_{l}}^{l}}{\unicode[STIX]{x2202}y}}.\end{eqnarray}$$
These developments show that the disturbance induced by the bubbles through the potential flow leads to a lift force in the presence of shear. Here, the part of the Reynolds stresses induced by the averaged wake of the bubble has been neglected because the development was focused on the transverse direction where the potential flow is larger. However, if the wake is slightly tilted (with respect to the liquid mean flow) due to the trajectory of the bubble, it may also contribute to the transversal Reynolds stresses. Finally, the total lift force is,
$$\begin{eqnarray}\displaystyle M^{L}=M_{\unicode[STIX]{x1D735}P}^{L}+M_{\text{Re}}^{L}\approx -\underbrace{\left(\frac{6\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D6FC}_{l}}{20}-4C_{\unicode[STIX]{x1D735}P}\unicode[STIX]{x1D6FD}\right)}_{C_{L}}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{l}\overline{u_{r}}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{u_{l}}^{l}}{\unicode[STIX]{x2202}y}}, & & \displaystyle\end{eqnarray}$$
where
$C_{\unicode[STIX]{x1D735}P}$
is a parameter. Pressure gradient and Reynolds stress lift forces have opposite behaviours. With the previous derivation, the pressure gradient lift force is connected to the bubble deformation whereas the Reynolds stress lift force is linked to the flow dynamics induced by the shear. This behaviour is in agreement with the physical explanation given by Adoua et al. (Reference Adoua, Legendre and Magnaudet2009). In this study, a dual origin of the lift force has been demonstrated. A part of the lift force comes from the vorticity induced by the liquid shear (flow dynamics), whereas the other part is related to vorticity created at the bubble interface (bubble deformation). In the case of spherical bubbles, experiments show that
$C_{L}$
is always positive (Tomiyama et al.
Reference Tomiyama, Tamai, Zun and Hosokawa2002) because the contribution of the liquid dynamics is stronger than the vorticity production due to small bubble deformation. In (3.28), obtained for quasi-spherical bubbles, the coefficient
$C_{L}$
must be positive. The parameter
$C_{\unicode[STIX]{x1D735}P}$
must therefore satisfy
$C_{\unicode[STIX]{x1D735}P}\leqslant 6\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D6FC}_{l}/80\unicode[STIX]{x1D6FD}$
. It seems then that the coefficient
$C_{\unicode[STIX]{x1D735}P}$
is related to the void fraction. Indeed, the equations (3.16) and (3.17) reflect the local action of surface tension on the liquid pressure (around the interface). Far from the bubble (in most of the liquid phase at moderate void fraction), the pressure is not disturbed by the interface. On average, the pressure in the liquid phase is therefore much lower than the local pressure around the interface. The parameter
$C_{\unicode[STIX]{x1D735}P}$
reflects this effect and that is why it can be expected to be very small and related to the void fraction (it depends on the volume on which the pressure is disturbed by the interface and on the decrease of the pressure in this volume).
Extending this partial result to complex flows, the reversal of the lift force for deformable bubbles may be explained by an increase of the pressure gradient lift force (related to the Weber number). Even if the preceding development is carried out for quasi-spherical bubbles, we assume (although the demonstration is more complex), that the contribution related to the effects of surface tension (characterised by
$C_{\unicode[STIX]{x1D735}P}$
) increases as the bubbles become more deformable until the inversion of the sign of the lift coefficient is reached. This interpretation is also qualitatively in agreement with the observations made in § 2.2. These derivations serve a pedagogical purpose in order to confirm the new definition of the lift force but do not aim at a new model for the coefficient.
We have identified the lift force in the simplified case of almost spherical bubbles in laminar shear flows. We expect that this definition of the lift force can be extended further to more complex flow regimes (turbulent channel flow etc.). This section dedicated to the lift force has given a local definition for the lift force and has confirmed that its common closure is adequate and consistent with the new definition proposed.
3.4 Identification of forces
Finally, using the pressure decomposition as in (3.13) and the Reynolds stresses, splitting as in (3.10), the trajectory equation of a fluid element in the vapour phase (3.6) becomes,
$$\begin{eqnarray}\displaystyle \begin{array}{@{}rcll@{}}\displaystyle \unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}{\displaystyle \frac{\text{D}\overline{\boldsymbol{u}_{v}}^{v}}{\text{D}t}}\ & =\ & \boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}} & \\ \ & =\ & \unicode[STIX]{x1D6FC}_{v}(\unicode[STIX]{x1D70C}_{v}-\langle \unicode[STIX]{x1D70C}\rangle )\boldsymbol{g} & \boldsymbol{M}^{\unicode[STIX]{x1D72B}}\\ \ & \ & -\displaystyle \,\unicode[STIX]{x1D6FC}_{v}{\displaystyle \frac{\text{D}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}}{\text{D}t}}+\unicode[STIX]{x1D6FC}_{l}{\displaystyle \frac{\text{D}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}}{\text{D}t}}-\unicode[STIX]{x1D6FC}_{v}^{2}\unicode[STIX]{x1D735}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l}) & \boldsymbol{M}^{\boldsymbol{A}\boldsymbol{M}}\\ \ & \ & +\,\unicode[STIX]{x1D6FC}_{v}[\overline{\unicode[STIX]{x1D70F}_{v}}^{v}-\overline{\unicode[STIX]{x1D70F}_{l}}^{l}]\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FC}_{v} & \boldsymbol{M}^{\unicode[STIX]{x1D749}\boldsymbol{D}}\\ \ & \ & +\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D70F}_{l}}^{l}+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{\unicode[STIX]{x1D70F}_{v}}^{v}] & \boldsymbol{M}^{\unicode[STIX]{x1D749}}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v}|_{turb}-\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}|_{turb}]\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v} & \boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}|_{turb}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v}|_{lam}-\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}|_{lam}]\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v} & \boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}|_{lam}\\ \ & \ & +\,\unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v} & \boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}}) & -\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\overline{p_{l}^{SP}}^{l}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\tilde{p_{l}})\boldsymbol{\cdot }\boldsymbol{e}_{x}\boldsymbol{e}_{x}+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D6FC}_{l}(\unicode[STIX]{x1D70C}_{l}-\unicode[STIX]{x1D70C}_{v})\boldsymbol{g} & \boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{D}}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}|_{turb})+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v}|_{turb})]\boldsymbol{\cdot }\boldsymbol{e}_{x}\boldsymbol{e}_{x} & \boldsymbol{M}_{\text{Re}}^{\boldsymbol{D}}|_{turb}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}|_{lam})+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v}|_{lam})]\boldsymbol{\cdot }\boldsymbol{e}_{x}\boldsymbol{e}_{x} & \boldsymbol{M}_{\text{Re}}^{\boldsymbol{D}}|_{lam}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\tilde{p_{l}})\boldsymbol{\cdot }\boldsymbol{e}_{y}\boldsymbol{e}_{y} & \boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}|_{turb})+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v}|_{turb})]\boldsymbol{\cdot }\boldsymbol{e}_{y}\boldsymbol{e}_{y} & \boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}|_{turb}\\ \ & \ & -\,\unicode[STIX]{x1D6FC}_{v}[\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}|_{lam})+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v}|_{lam})]\boldsymbol{\cdot }\boldsymbol{e}_{y}\boldsymbol{e}_{y} & \boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}|_{lam},\end{array} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{M}^{\unicode[STIX]{x1D72B}}$
and
$\boldsymbol{M}^{\boldsymbol{A}\boldsymbol{M}}$
are respectively the buoyancy and the added-mass forces;
$\boldsymbol{M}^{\unicode[STIX]{x1D749}\boldsymbol{D}}$
and
$\boldsymbol{M}^{\unicode[STIX]{x1D749}}$
are the viscous dispersion and the viscous forces (both negligible);
$\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}|_{turb}$
and
$\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}|_{lam}$
are the turbulent and laminar parts of the inappropriately called turbulent dispersion;
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}$
is the new laminar dispersion force. The single-phase pressure gradient is simplified because the reference pressure
$\widetilde{p_{l}}$
is constant. The drag force
$\boldsymbol{M}^{\boldsymbol{D}}$
is comprised of a pressure contribution
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{D}}$
and Reynolds stress parts (
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{D}}|_{lam}$
for the laminar fluctuations and
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{D}}|_{turb}$
for the turbulent ones). Similarly, the lift force
$\boldsymbol{M}^{\boldsymbol{L}}$
is comprised of a pressure contribution
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
and Reynolds stress parts (
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}|_{lam}$
for the laminar fluctuations and
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}|_{turb}$
for the turbulent ones):
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}^{\boldsymbol{D}}=\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{D}}+\boldsymbol{M}_{\text{Re}}^{\boldsymbol{D}}|_{lam}+\boldsymbol{M}_{\text{Re}}^{\boldsymbol{D}}|_{turb}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}^{\boldsymbol{L}}=\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}+\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}|_{lam}+\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}|_{turb}. & \displaystyle\end{eqnarray}$$
In addition to the lift force definition, equation (3.29) also presents a definition for the drag force comprised of pressure and turbulent parts. Validation elements for this definition are provided in § 4.3. In (3.29), the total liquid pressure
$\overline{p_{l}^{\prime }}^{l}$
is split into an equivalent single-phase pressure in the absence of bubbles
$\overline{p_{l}^{SP}}^{l}$
and a pressure induced by the bubbles through surface tension effects
$\overline{p_{l}^{b}}^{l}$
. In practice, these two pressures are indistinguishable if
$\overline{p_{l}^{SP}}^{l}\neq 0$
; the decomposition of the Reynolds stresses into turbulent and laminar parts is not possible either because spatial and temporal fluctuations cannot be distinguished in the present configuration (see du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019), Amoura et al. (Reference Amoura, Besnaci and Risso2017) for adapted configurations). These limits complicate the interpretation of the interfacial forces.
In (3.29), a new definition of the added-mass force is given. Physically, this force is related to unsteady effects when both phases have different accelerations (Geurst Reference Geurst1985). Further investigations are necessary to analyse the link to the classical formulation of the added-mass force. However, we believe that a procedure similar to the analysis of the lift force can be performed. Firstly, part of the added-mass force may be related to the term
$-\unicode[STIX]{x1D6FC}_{v}^{2}\unicode[STIX]{x1D735}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})$
. Indeed, the acceleration of the liquid or vapour is necessarily related to a pressure gradient. Physically, the difference in the pressure gradients between liquid and vapour seems to be closely related to an added-mass effect for an unsteady flow. The connection between a pressure difference and the added-mass force has already been made by Prosperetti & Jones (Reference Prosperetti and Jones1984) for instance. Note that connecting the liquid pressure gradient
$-\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}\overline{p_{l}^{SP}}^{l}$
to an acceleration can similarly be interpreted as a contribution to the Tchen force. Thus, by combining the unsteady terms involving material derivatives and the terms comprised of a pressure gradient, we believe that it is possible to recover the conventional added-mass force. Ultimately, the form of material derivatives will certainly remain a matter of debate because the classical formulation of the added-mass force involves a convective derivative related to the mean bubble velocity (and not to the mean velocity of the fluid element as in (3.29)) (Magnaudet & Eames Reference Magnaudet and Eames2002). At this point, it should be noted that (3.29) is not the trajectory equation of the bubble itself but the trajectory equation of a fluid element in the vapour phase (there is no integration on the volume of a bubble). In doing so, the definition of the added-mass force in the particle approach and the definition proposed here could be slightly different while characterising the same effect. In (3.29), some other forces, such as the history force, are not directly observable. For instance, we can expect to find the history force in a term involving viscosity, such as
$\boldsymbol{M}^{\unicode[STIX]{x1D749}}$
or
$\boldsymbol{M}^{\unicode[STIX]{x1D749}\boldsymbol{D}}$
, which are negligible here (as is the history force in the case of bubbles at high Reynolds numbers (Magnaudet, Fabre & Rivero Reference Magnaudet, Fabre and Rivero1995)).
It is worth noting that the interpretation of (3.29) is not complete. In this paper, the work is focused on the lateral forces at the statistical steady state which drive the void-fraction profile (lift and dispersion forces). However, all the classical terms are inevitably retrievable since (3.29) is exact. We further believe that these forces can be identified and separated. Further work is necessary to formally highlight unsteady forces (added-mass, Tchen and history forces) as has been done here for the lift force. Successfully identifying all the known forces in simplified cases in this equation could finally give them a more general definition. This definition would not challenge the previous work on the added mass or on the drag force and so on (Geurst Reference Geurst1986; Wallis Reference Wallis1990; Drew & Passman Reference Drew and Passman1999), on the contrary, the two approaches are complementary. An important contribution of this new formulation is, for example, to measure the different forces in complex situations and to prevent the effects of error compensation or the shortcomings of models that assemble closures obtained from academic cases. Finally, a first definition of the forces from local quantities is proposed here (see table 2 for a summary). Future adjustments may be necessary for local definitions to coincide with formulations proposed in the literature. This equation can become an effective tool to understand and model forces. In the following sections, this equation is used to study the lift force, the dispersion forces and the perspectives for the modelling in RANS Euler–Euler calculation codes.
Table 2. Local definitions of forces without the particle hypothesis versus models of interfacial forces from the particle approach. Some identifications between the two viewpoints have been demonstrated (for the lift force in particular). The closure identification for the turbulent dispersion force makes sense (except for the coefficient). For the others (drag and added-mass forces), nothing has been formally demonstrated but the identification seems physically consistent, as discussed previously.
Figure 5. Statistical results from DNS: contributions to (3.6) projected on the
$y$
-axis for all simulated cases. Comparison between the definition of the lift force from local quantities and the classical form provided in table 2 with adjusted parameters
$C_{L}^{\star }$
. (a) D127,
$C_{L}^{\ast }=0.017$
, (b) D180,
$C_{L}^{\ast }=0.011$
, (c) D180g8,
$C_{L}^{\ast }=0.054$
, (d) S180,
$C_{L}^{\ast }=0.01$
, (e) S180g8,
$C_{L}^{\ast }=0.12$
.
4 Results
In the previous section, we proposed a new definition of the forces that appear in the momentum transfer term in two-fluid models. This section applies the previously described principles to DNS results to understand the physics of these forces and their importance. The § 4.1 presents the study of the lift force in the complex configuration of a channel flow. Based on the DNS data described in § 2, § 4.2 shows that the impact of the new non-turbulent dispersion force is tremendous in the migration process. Section 4.3 strengthens the definition of the drag force by comparing it with its most common closure.
4.1 Lift force in complex flows
Section 3 highlighted the content of the lift force
$\boldsymbol{M}^{\boldsymbol{L}}$
comprised of a pressure lift force
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
and a Reynolds stress lift force
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
.
The inseparability of
$\overline{u_{l}^{\prime }u_{l}^{\prime }}^{l}|_{lam}$
and
$\overline{u_{l}^{\prime }u_{l}^{\prime }}^{l}|_{turb}$
complicates the interpretation of
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
in a turbulent channel flow. To date, the only configurations found to study separately those fluctuations are fixed arrays of solid spheres (or bubbles), in which the gradient of the Reynolds stress
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
and the gradient of pressure
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
are zero if the distribution is spatially uniform. Figure 5 compares the lift definition introduced in (3.29) and the classical closure
$C_{L}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{l}u_{r}(\unicode[STIX]{x2202}\overline{u_{l}}^{l}/\unicode[STIX]{x2202}y)$
(normalised with an appropriate
$C_{L}=\pm cte$
for each case). Case S180 is almost laminar (see figure 1). Figure 5(d) shows that
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
is non-zero due to the non-turbulent fluctuations. In that case, as expected from the theoretical development in § 3.3, the classical closure gives satisfactory results. For case S180g8, in figure 5(e),
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
is bigger than in case S180 because of an increase of the relative velocity and of the turbulent part of the lift (turbulent fluctuations induced by interactions and instabilities of wakes are observable in figure 1). Even in the presence of the turbulent lift force, the classical model works very well for case S180g8 with an appropriate
$C_{L}$
. Hence, the modelling of
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
can be generalised
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}}=(C_{L}^{\text{Re}}|_{lam}+C_{L}^{\text{Re}}|_{turb})\unicode[STIX]{x1D70C}_{l}\unicode[STIX]{x1D6FC}_{v}u_{r}{\displaystyle \frac{\unicode[STIX]{x2202}\overline{u_{l}}^{l}}{\unicode[STIX]{x2202}y}}, & & \displaystyle\end{eqnarray}$$
where
$C_{L}^{\text{Re}}|_{lam}$
and
$C_{L}^{\text{Re}}|_{turb}$
are respectively the laminar and turbulent lift force coefficients. To complete the analysis, the role of the pressure lift force needs investigations. The inseparability of
$p_{l}^{b}$
and
$p_{l}^{SP}$
makes it very difficult to compare the pressure lift force to its adapted closure in the turbulent channel. The disturbance of the pressure gradient induced by the single-phase flow prevents us from making a comparison with the classical lift model. This can explain why the pressure gradient is not proportional to
$\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{l}u_{r}\unicode[STIX]{x2202}\overline{u_{l}}^{l}/\unicode[STIX]{x2202}y$
in the spherical cases S180 and S180g8 in figures 5(d) and 5(e).
For deformable bubbles (cases D127, D180 and D180g8), the lift force (
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}+\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
) changes its orientation (see § 2.2). The reversal of the total lift force is concurrent with the reversal of
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
and
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
(see figure 5). For
$y/h<0.2$
,
$C_{L}^{\text{Re}}=C_{L}^{\text{Re}}|_{lam}+C_{L}^{\text{Re}}|_{turb}>0$
whereas for
$y/h>0.7$
, the result suggests that
$C_{L}^{\text{Re}}<0$
. In between these two regions, for
$0.2<y/h<0.7$
, the sign of
$C_{L}^{\text{Re}}$
changes;
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
works in the opposite direction to
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
(as shown by the theoretical development in the previous section). Nevertheless, the sum of
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
and
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
is not zero and constitutes the effective lift force. Then, the lift orientation is governed by the relative importance of the lift coefficients
$C_{L}^{\text{Re}}|_{lam}$
,
$C_{L}^{\text{Re}}|_{turb}$
and
$C_{L}^{\unicode[STIX]{x1D735}P}$
and their behaviours regarding the bubble deformability. For now, those complex behaviours of
$\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}$
and
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
for deformable bubbles are not explained. The present analysis is not sufficient to fully understand the effective mechanism. Pending new ways to separate physical phenomena in simplified cases, global studies of the lift force, as done in § 2.2 or in Tomiyama et al. (Reference Tomiyama, Tamai, Zun and Hosokawa2002), remain the only way to model it (with a particle approach). However, the new balance equation of forces is an alternative to this classical particle approach. It shows that only the Reynolds stress tensor and surface tension effects need models for the purposes of lift force prediction. Further research (for instance on a non-uniform array of fixed bubbles) based on the local definitions given in (3.29) should be considered. The new balance equation can become an effective tool to model the reversal of the lift force in complex flows.
4.2 Dispersion forces
Wall-normal components of the interfacial forces given in 3.6 are plotted on figure 5 for all simulated cases at the statistical steady state. First of all, the residue of the equation is negligible for all simulated cases except for case S180 for which slight variations remain. This means that the flow has reached a bubble distribution (i.e. a void-fraction profile) for which there is a balance between all the forces of the equation. The void fraction profiles have two essential features. The first one is the position of the void-fraction peak resulting from the action of the lift force on the bubbles. The second one is the spreading due to the dispersion forces that are the subject of this section.
Dispersion forces, comprised of laminar and turbulent parts, have been identified in § 3. In the interest of readability, the turbulent dispersion force refers in the following to its definition in table 2 even if it contains a part of laminar dispersion force due to the separation of the Reynolds stresses. Then, the laminar dispersion force is defined only by the difference between the interfacial pressure force
$\boldsymbol{M}^{\boldsymbol{P}}$
and the surface tension force
$\boldsymbol{M}^{\unicode[STIX]{x1D748}}$
(see 3.6). Both laminar and turbulent dispersion forces have an impact on the spreading of the void-fraction profile. In figure 5, the turbulent dispersion force (
$\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}$
) is always much smaller than the laminar dispersion force (
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}$
). This observation holds for all simulated cases whether bubbles are deformable or not. The relative importance of the turbulent dispersion, compared to the laminar part, does not increase significantly with the turbulence magnitude. Indeed, even for the most turbulent case (D180g8), its contribution is negligible (as in cases S180 and S180g8). At its maximum, for case D180, we have
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}\approx 3\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}$
.
In the following, ‘non-Eulerian effect’ refers to a phenomenon which is difficult to describe in a point-particle formalism because it occurs at scales below the bubble diameter and it relies on the finite size of the bubbles. In the two cases which present non-Eulerian effects due to the presence of bubbles in the near-wall region (S180 and S180g8), the laminar dispersion force becomes stronger than the pressure and lift forces. Nevertheless, the interpretation of the forces in this region changes. For instance, in figure 2(a), the void-fraction peak at
$y/h\approx 0.15$
for case S180 can easily be interpreted as the expression of a repulsive force from the wall (see, e.g., Antal, Lahey & Flaherty Reference Antal, Lahey and Flaherty1991). Actually, the position of the void-fraction peak is determined by the finite size of the bubbles and corresponds here to the location of the centre of gravity of wall-peaking bubbles (dotted line in figure 2(a) corresponds to the bubble diameter). Thus, the apparent repulsive effect is not a hydrodynamic force but a manifestation of the surface tension force (Moraga et al.
Reference Moraga, Larreteguy, Drew and Lahey2006). It comes from the averaging operator and the use of a particle approach (Lubchenko et al.
Reference Lubchenko, Magolan, Sugrue and Baglietto2017). The bubble as a whole is not necessarily subjected to a repulsive force from the wall. Surface tension acts to restore equilibrium and prevents high void-fraction values below a bubble radius. The same kind of interpretation could happen for the laminar dispersion force in the near-wall region. For
$y<d_{b}$
,
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}$
is a surface tension force related to a non-Eulerian effect due to the finite size of the bubbles.
As a conclusion, the non-turbulent dispersion force is one of the most important terms for the void-fraction profile. Hence, it needs to be modelled, whereas it is totally neglected in the most common models. The laminar dispersion force is comprised of two terms: the interfacial pressure force
$\boldsymbol{M}^{\boldsymbol{P}}$
and the surface tension force
$\boldsymbol{M}^{\unicode[STIX]{x1D748}}$
. Concerning the interfacial pressure force, its definition is,
$$\begin{eqnarray}\boldsymbol{M}^{\boldsymbol{P}}=-\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}.\end{eqnarray}$$
This force is related to an anti-dispersion process because it is proportional to
$-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}$
. This behaviour is quite uncommon. Indeed,
$\overline{p_{v}}^{v}-\overline{p_{l}}^{l}>0$
due to the increase of pressure inside the bubble (i.e. the Laplace law). A similar force has already been studied by Stuhmiller (Reference Stuhmiller1977), Pauchon & Banerjee (Reference Pauchon and Banerjee1986) and more recently by Vaidheeswaran & de Bertodano (Reference Vaidheeswaran and de Bertodano2016) for ellipsoids within a potential flow. These authors have understood that pressure forces may be responsible for the dispersion process. An anti-dispersion process is a surprising result and no experiment has isolated this behaviour. However, the laminar dispersion force is also comprised of the surface tension term which is a dispersive force always larger than the interfacial pressure force. Hence, the combined effect of both of them is always dispersive. The surface tension force is defined as,
$$\begin{eqnarray}\boldsymbol{M}^{\unicode[STIX]{x1D748}}=\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70E}\overline{\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}.\end{eqnarray}$$
This force is related to a dispersion process and is partly reduced by the interfacial pressure force. Indeed, for monodisperse spherical bubbles at mechanical equilibrium, the curvature is constant. The pressure jump between the phases is locally dictated by the Laplace law:
$p^{+}-p^{-}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}$
where
$p^{+}$
is the pressure next to the interface in the vapour phase and
$p^{-}$
in the liquid phase. If the bubbles are sufficiently small and the liquid pressure field close to the interface can be assimilated to the pressure of the liquid phase (
$\overline{p_{l}}^{l}=p^{-}$
), the average pressure jump between the phases is dictated by a Laplace law generalised to the entire phase (
$\overline{p_{v}}^{v}-\overline{p_{l}}^{l}=\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}=Cte$
). In these conditions we get a zero laminar dispersion force,
$$\begin{eqnarray}\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}=\boldsymbol{M}^{\unicode[STIX]{x1D748}}+\boldsymbol{M}^{\boldsymbol{P}}=\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70E}\overline{\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}=0.\end{eqnarray}$$
In our cases, we show that the laminar dispersion force is not negligible. It can be explained by the fact that, in our cases, the Laplace law cannot be generalised to the whole vapour phase when the bubbles are of significant size. We therefore clearly observe a limitation of the particle hypothesis (or of the hypothesis of small bubbles). Our results show that the difference between
$\boldsymbol{M}^{\boldsymbol{P}}$
(comprised of liquid and vapour pressures) and
$\boldsymbol{M}^{\unicode[STIX]{x1D748}}$
is at the origin of a dispersion force. This conclusion contrasts with Prosperetti & Jones (Reference Prosperetti and Jones1984) who show that surface tension effects are compensated by the pressure effects of the vapour phase, using an extended Laplace law, and that the bubble dynamics is entirely determined by the liquid phase. The work of Prosperetti & Jones (Reference Prosperetti and Jones1984) is in agreement with the classical formulation of the interfacial pressure force described in Pauchon & Banerjee (Reference Pauchon and Banerjee1986) and Vaidheeswaran & de Bertodano (Reference Vaidheeswaran and de Bertodano2016), which thus takes into account only the liquid pressure gradient. These studies therefore assume that vapour pressure has no influence on the bubbles dynamics. However, for large-sized bubbles in which the pressure is not constant, pressure effects generated by a stress on one side of the bubble can be transmitted to the other side via the pressure in the vapour phase. Thus, it appears that, even for small density and viscosity in the vapour, the flow dynamics is not entirely determined by the liquid phase. In our particular cases, the interfacial transfer is affected by the pressure distribution within the bubble.
4.3 Drag force
In figure 6, all the streamwise forces of (3.6) projected in the axial direction have been plotted for all the DNS cases at the statistical steady state. The residue of the equation is insignificant for all simulated cases. In the absence of non-Eulerian effects (all cases except S180), the drag force is widely predominant and compensates for the buoyancy of the vapour phase. Classically, the drag force is related to the square of the relative velocity,
$$\begin{eqnarray}\displaystyle \boldsymbol{M}^{\boldsymbol{D}}=-\frac{3}{4d_{b}}\unicode[STIX]{x1D6FC}_{v}C_{D}\unicode[STIX]{x1D70C}_{l}|\boldsymbol{u}_{r}|\boldsymbol{u}_{r}, & & \displaystyle\end{eqnarray}$$
where
$C_{D}$
is the drag coefficient. In figure 6, this kind of closure fits very well with the drag force for all simulated cases (with an adjusted coefficient
$C_{D}^{\ast }$
for each case). It confirms the definition given in (3.29).
Figure 6. Statistical results from DNS: contributions to (3.6) projected onto the
$x$
-axis for all simulated cases. Comparison between the definition of the drag force from local quantities and its current model. (a) D127,
$C_{D}^{\ast }=0.1$
, (b) D180,
$C_{D}^{\ast }=0.05$
, (c) D180g8,
$C_{D}^{\ast }=0.44$
, (d) S180,
$C_{D}^{\ast }=0.22$
, (e) S180g8,
$C_{D}^{\ast }=1.73$
.
5 RANS Euler–Euler application
The lessons learned from previous sections are finally applied here to the RANS Euler–Euler formalism to provide some insight into the improvement of RANS models and computations.
The averaged Navier–Stokes equations can be written in the Euler–Euler RANS formalism (see (A 21a ) and (A 21b ) from §§ A.1 and A.2),
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}\overline{\boldsymbol{u}_{v}}^{v}) & = & \displaystyle -\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\boldsymbol{g}+\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}\overline{\boldsymbol{u}_{l}}^{l}) & = & \displaystyle -\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}-\unicode[STIX]{x1D6FC}_{l}\overline{\unicode[STIX]{x1D70F}_{l}}^{l})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\boldsymbol{g}+\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=-\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}+\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})]-\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}}). & & \displaystyle\end{eqnarray}$$
Equations (5.1a
) and (5.1b
) are written in the Euler–Euler RANS two-fluid one pressure formalism (with a liquid pressure gradient in the vapour momentum equation) respectively in the vapour and liquid phases. In the RANS Euler–Euler approach, the liquid pressure in the absence of bubbles
$\overline{p_{l}^{SP}}^{l}$
is included in the resolution of the system whereas the pressure inside the bubbles
$\overline{p_{v}}^{v}$
and the liquid pressure induced by bubbles through surface tension effects
$\overline{p_{l}^{b}}^{l}$
have to be closed. The interpretation of the resolved pressure field in a classical RANS Euler–Euler calculation is thus tricky. The part of pressure due to surface tension has an impact on the balance equation of forces (mainly on the lift force, see § 3.3). Thus, this part, which is not directly solved, is considered through the interfacial forces.
In this two-fluid model, the interface is not explicitly represented. This means that the presence of bubbles is accounted for by an averaged void fraction so that each bubble position is lost during the averaging procedure. Despite this formalism, equations (5.1a
) and (5.1b
) are equivalent to the exact averaged Navier–Stokes equations under conditions depicted in §§ A.1 and A.2;
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
and
$\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
are the interfacial momentum transfer terms in the vapour and liquid phases respectively. Because of surface tension, the action–reaction principle cannot be applied in a straightforward fashion and
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}\neq -\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
(5.1c
). This disequilibrium is taken into account here through pressure and surface tension terms. Nevertheless, this formulation can tend to
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=-\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}-\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}})$
by assuming an extended Laplace law:
$\unicode[STIX]{x1D70E}\overline{\unicode[STIX]{x1D705}}=\overline{p_{v}}^{v}-\overline{p_{l}}^{l}=Cte$
(but this classical approximation is a source of serious bias so it will not be used in the following). In the limit of this local equilibrium, the interfacial pressure jump and the surface tension forces are balanced. If all surface tension effects are neglected, this formulation tends to the classical
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=-\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
. Indeed, without a surface tension effect, a local equilibrium between liquid and vapour pressures is acceptable
$\overline{p_{v}}^{v}-\overline{p_{l}}^{l}=0$
; the pressure induced by surface tension
$\overline{p_{l}^{b}}^{l}$
is proportional to
$\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}$
and it vanishes with the surface tension along with the source term
$\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}$
. Hence, the difference between phases is a pure surface tension effect which characterises the impact of interfaces on the momentum budget. For instance,
$\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}$
transfers the momentum from the region where
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}>0$
(
$\unicode[STIX]{x1D705}<0$
by definition) to the region where
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}<0$
(but its integrated value over the whole domain is zero). Similar behaviours of
$\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}})$
and
$\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})]$
are expected. Thus the difference between phases can locally act as a momentum source or as a momentum sink. Removing this term from the momentum budget would result locally in a poor assessment of the velocity field.
In the previous section, we have shown that the total interfacial force is given by (3.6) which has a short form,
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}=-\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}})+\boldsymbol{M}^{\unicode[STIX]{x1D72B}}+\boldsymbol{M}^{\boldsymbol{D}}+\boldsymbol{M}^{\boldsymbol{A}\boldsymbol{M}}+\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}+\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}+\boldsymbol{M}^{\boldsymbol{L}}. & & \displaystyle\end{eqnarray}$$
The commonly used models to close equation (5.2) are given in table 2. The lift closure of Tomiyama et al. (Reference Tomiyama, Tamai, Zun and Hosokawa2002) catches a part of the phenomenon due to the pressure
$\overline{p_{l}^{b}}^{l}$
induced by surface tension (as discussed in § 3.3). Nevertheless, it does not take into account the turbulent part of the lift force. For the turbulent dispersion force, the closure proposed by Laviéville et al. (Reference Laviéville, Mérigoux, Guingo, Baudry and Mimouni2015) has been chosen because it takes into account most of the turbulent contribution of the classical forces. Nevertheless, contrary to what the definition in table 2 suggests, the turbulent fluctuations in the vapour phase are neglected in the model of Laviéville et al. (Reference Laviéville, Mérigoux, Guingo, Baudry and Mimouni2015). A lot of closures exist for the drag force in different configurations (spherical bubbles, deformable bubbles, highly turbulent flows, swarm effects...). For instance, Ishii & Zuber (Reference Ishii and Zuber1979) propose a generic model which predicts the drag force accounting for all these phenomena. Finally, the main lessons of this table are the lack of a model for the laminar dispersion force and the turbulent part of the lift force which both have a significant impact on the bubble lateral migration (see §§ 4.1 and 4.2).
Furthermore, the interfacial momentum transfer in the RANS Euler–Euler formalism (
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
) is classically comprised of drag, lift, added-mass and turbulent dispersion forces. To check whether this common practice is adequate, the content of
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
is studied based on our DNS configurations. In order to rebuild
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
from
$\boldsymbol{M}_{v}$
, and
$\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
from
$\boldsymbol{M}_{l}$
, relations are demonstrated in § A.2,
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}_{v}=\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}})+(\overline{p_{l}}^{l}-\widetilde{p_{l}})\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}+\unicode[STIX]{x1D735}\left[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}_{l}=\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}+\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}})-(\overline{p_{l}}^{l}-\widetilde{p_{l}})\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}. & \displaystyle\end{eqnarray}$$
Then, by using (3.2), (3.4) and (5.3), we get,
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}} & = & \displaystyle \boldsymbol{M}^{\boldsymbol{D}}+\boldsymbol{M}^{\boldsymbol{A}\boldsymbol{M}}+\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}+\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}+\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}+\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}+\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}})+\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}-\boldsymbol{M}^{\unicode[STIX]{x1D72B}},\end{eqnarray}$$
where
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{u_{v}^{\prime }u_{v}^{\prime }}^{v})$
is related to the Reynolds stresses of the dispersed phase. For the liquid phase, by using (5.1c
) and the definitions of the interfacial forces, we get
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=-\boldsymbol{M}^{\boldsymbol{D}}-\boldsymbol{M}^{\boldsymbol{A}\boldsymbol{M}}-\boldsymbol{M}^{\boldsymbol{T}\boldsymbol{D}}-\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}-\boldsymbol{M}_{\text{Re}}^{\boldsymbol{L}}+\frac{\unicode[STIX]{x1D6FC}_{l}}{\unicode[STIX]{x1D6FC}_{v}}(\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}+\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}), & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
and
$\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
are not the total interfacial forces applied on their phases. Their proper definitions are given by (5.5) and (5.7). In practice in the Euler–Euler framework, the hypothesis
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=-\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=\boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}-\boldsymbol{M}^{\unicode[STIX]{x1D72B}}+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}})$
is always assumed. By doing so, the impact of the additional term
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}$
related to velocity fluctuations in the vapour phase is neglected (see (5.6)) in the momentum budgets (necessary for the prediction of the velocity field). Even for
$\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
,
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}$
is related to the fluctuations in the vapour phase while fluctuations in the liquid phase might have been expected. This happens because the forces have been defined with the balance equation of forces for the dispersed phase (and not for the liquid phase). This choice is in agreement with the literature on interfacial forces (the closure of interfacial forces in table 2 are proportional to the void fraction
$\unicode[STIX]{x1D6FC}_{v}$
because they are written for the vapour phase). However, starting from the equilibrium equation of forces in the liquid phase, a formulation equivalent to (5.7) can be obtained. It involves fluctuations in the liquid phase instead of fluctuations in the vapour phase for
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}$
. In this formulation, the forces must be corrected with a multiplicative factor
$\unicode[STIX]{x1D6FC}_{l}/\unicode[STIX]{x1D6FC}_{v}$
. The choice between these two equivalent formulations is based on our ability to model the fluctuations in the liquid or in the vapour phase.
Figure 7. Contributions to (5.5) and (5.7) in the wall-normal direction for case D127. Comparison with the pressure gradient in the liquid phase (
$-\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\overline{p_{l}}^{l}=-\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}\overline{p_{l}^{SP}}^{l}-(\unicode[STIX]{x1D6FC}_{l}/\unicode[STIX]{x1D6FC}_{v})\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
). (a) Liquid: equation (5.7), (b) vapour: equation (5.5).
Equations (5.5) and (5.7) are plotted in figure 7 projected in the wall-normal direction for case D127. The balance between the forces depends on the phase because
$\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}\neq -\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
due to the impact of surface tension. We have already said that we cannot distinguish
$\overline{p_{l}^{SP}}^{l}$
and
$\overline{p_{l}^{b}}^{l}$
in the present configuration. For this reason,
$\boldsymbol{M}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
and
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
are plotted together in figure 7.
In the liquid phase, the budget is dominated by the difference between the phases attributed to the surface tension effect responsible for the laminar dispersion force
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}$
and for the gradient of bubble pressure
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
(the other cases are not shown but have similar behaviour). Unfortunately, because we cannot distinguish
$\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
and
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
, we cannot know whether
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
compensates for the action of
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}$
or not. For comparison, the total lateral pressure gradient is plotted in figure 7 (dotted line). The order of magnitude of this term suggests that
$\boldsymbol{M}_{\unicode[STIX]{x1D735}\boldsymbol{P}}^{\boldsymbol{L}}$
may be in the same range as
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}$
. If so, the difference between the phases due to the laminar dispersion and to the gradient of bubble pressure could possibly balance each other out. In such a case, we would recover
$\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}\approx -\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
, as classically stated in the RANS Euler–Euler framework, and the surface tension effects could be neglected in the difference between phases in (5.1c
) (but not in the definitions of the interfacial forces, as already described in § 4.2). Pending new tools to separate the undisturbed pressure and the pressure induced by surface tension effects, this discussion cannot be concluded and we are far from being able to neglect this disequilibrium.
On the other hand, for the vapour phase,
$\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}$
is a balance between lift, laminar and turbulent dispersion forces as well as
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}$
, which is not negligible. For the other cases, the significance of
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}$
in the momentum budget is also confirmed (not shown here). This force has to be considered in the vapour momentum equation. If neglected, the estimation of the momentum magnitude (and then of the mean velocity) could be inaccurate. To the best of our knowledge,
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}$
is classically neglected in all practical applications.
As a conclusion, the new terms
$\boldsymbol{M}^{\boldsymbol{e}\boldsymbol{x}\boldsymbol{t}\boldsymbol{r}\boldsymbol{a}}$
and
$\boldsymbol{M}^{\boldsymbol{L}\boldsymbol{D}}$
, which are routinely passed over in RANS applications, should be considered for further debate and modelling.
6 Conclusion
Transversal forces have been investigated from DNS data. Five DNSs of turbulent bubbly flows have been designed for the study of bubble migration. Very different situations are discussed. Core, wall and intermediate peaking are depicted. A new criterion has been proposed for the reversal of the lift force orientation which is classically based on the Eötvös number. Based on the DNS data, the Weber number seems more appropriate to describe this inversion. Based on the local and averaged formulations of the Navier–Stokes equations, a new balance equation for the forces is established. We propose a complementary approach to the paradigm of the particle hypothesis (assimilation of a bubble to a point-size particle or to a bubble of negligible size) which neglects in particular surface tension effects. For large-sized bubbles, the averaged momentum transfer term that requires closure is no longer directly derived from the average interfacial force acting on the bubbles. Instead, it is related to the mean force acting on a fluid element of the vapour phase. Then, we do not study here the trajectory equation of the entire bubble but of a fluid element of the vapour phase. This equation gives a local definition for all the forces acting on this element. Some of these definitions have been shown to reflect the usual meaning of interfacial forces (lift, drag and turbulent dispersion) in the classical approach. Other forces, especially unsteady forces, still need work to be identified (added-mass, Basset or Tchen forces). New forces are introduced and their physical meanings are explained. They are particularly relevant for the prediction of the void-fraction distribution. Furthermore, the new balance equation of forces is complementary to the classical particle approach. In particular, based on the local definitions, it becomes possible to measure experimentally or numerically some of the forces in all cases (assuming we have access to the Reynolds stresses or the pressure decomposition for some terms). This balance equation also shows that only the Reynolds stress tensor for each phase and surface tension effects need models. Hence, several strategies are possible:
(i) To bridge the gaps between the continuous and the particle approach and to benefit from the numerous literature on force modelling (classically based on particle hypothesis).
(ii) To use solely the new local formulation as a whole. The closure problem would then be transferred from force modelling to the developments of models for the Reynolds stresses of both phases, for pressures
$\overline{p_{l}^{b}}^{l}$
,
$\overline{p_{v}}^{v}$
and for
$\overline{\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}$
.(iii) To develop an intermediate approach by completing the existing models for interfacial forces with the understanding allowed by the new definitions. We could for instance benefit from the useful literature on unsteady force closures and propose an improvement of the lift force closure based on the force definitions established in this paper. More generally, drag force, lift force and Reynolds stresses are linked by the new balance equation established in this paper. Equation (3.29) is then a compatibility constraint for the models which should be adjusted to satisfy this criterion.
We also have shown that the new definition of the lift force is in agreement with its classical closure by analytical derivations under the simplified condition of an isolated bubble in a laminar shear flow, and that it is naturally capable of transitioning from positive to negative values depending on the bubble deformability and background turbulence. The lift force can be split into a turbulent and a laminar contribution. In the literature, there is still a lack of reliability concerning the closure of what we have called the turbulent lift coefficient. The study of the turbulent lift force is complicated due to the inseparability of turbulent and non-turbulent Reynolds stresses. Future calculations with fixed bubbles (Amoura et al. Reference Amoura, Besnaci and Risso2017) could be considered to study the reversal of the lift force for instance.
The definition of the non-turbulent dispersion force has also been provided. To the best of our knowledge, the laminar dispersion force has never been studied whereas it is a significant contribution to the migration process. In comparison with this laminar dispersion force, the turbulent dispersion is always smaller, if not negligible, in the configurations considered. Different known physical behaviours could be modelled via this laminar dispersion force: the horizontal clustering of spherical bubbles in laminar flows and the oscillating trajectories of deformable bubbles. This force is comprised of surface tension and interfacial pressure terms. The interfacial pressure term acts as an anti-dispersion process whereas the surface tension force acts as a stronger dispersion. The combined effect of both is hence a dispersion force. This force will have to be modelled to be integrated into averaged calculation codes. This goal will be the focus of future investigations.
Acknowledgements
The authors would like to convey their sincere thanks to GENCI and the TGCC for providing the necessary computational resources for performing the DNS calculations. This work was granted access to the HPC resources of CINES under the allocations A0022b7712 and A0042b7712 made by GENCI. Also, we would like to acknowledge the CEA/DEN for the HPC resources allocated on COBALT super-computer at TGCC.
Appendix A. From the Euler–Euler RANS two-fluid formulation to continuous equations
A.1 Euler–Euler RANS two-fluid formulation
This appendix presents the derivation of the two-fluid model with specific care to preserve all interfacial forces responsible for momentum transfer. Our goal is to clarify the definitions of the total force applied on the vapour phase and the interfacial momentum transfer term in the momentum equations. The classical two-fluid model can be derived from a particle approach. It considers the bubble as a point-size particle of mass
$m_{p}$
on which different forces are applied. The trajectory of a bubble is then driven from a Lagrangian viewpoint by the trajectory equation (A 1)
$$\begin{eqnarray}\displaystyle m_{p}\frac{\text{d}\boldsymbol{u}_{p}}{\text{d}t}=m_{p}\boldsymbol{g}+\boldsymbol{F}^{\boldsymbol{S}_{\boldsymbol{b}}}, & & \displaystyle\end{eqnarray}$$
where
$\text{d}/\text{d}t$
is the material derivative, with
$$\begin{eqnarray}\displaystyle \boldsymbol{F}^{\boldsymbol{S}_{\boldsymbol{b}}}=\oint _{S_{b}}(-p_{l}\boldsymbol{n}+\unicode[STIX]{x1D749}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}S, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{F}^{\boldsymbol{S}_{\boldsymbol{b}}}$
are the surface forces;
$\boldsymbol{F}^{\boldsymbol{S}_{\boldsymbol{b}}}$
is comprised of particle forces such as lift, drag, added-mass, pressure and viscous forces. Surface tension forces are classically neglected by arguing that their integrated value over a closed surface is equal to zero. In order not to neglect several crucial quantities from the outset,
$\boldsymbol{F}^{\boldsymbol{S}_{\boldsymbol{b}}}$
is considered as an unknown in the following.
Simonin (Reference Simonin2000) shows the following transport equation for any variable
$\unicode[STIX]{x1D731}$
,
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\unicode[STIX]{x1D731}}^{\dot{v}}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{p}\unicode[STIX]{x1D731}}^{\dot{v}})=\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\frac{\text{d}\boldsymbol{u}_{p}}{\text{d}t}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D731}}{\unicode[STIX]{x2202}\boldsymbol{u}_{p}}}}^{\dot{v}},\end{eqnarray}$$
where
$\dot{\unicode[STIX]{x1D6FC}}_{v}$
is the averaged void fraction,
$\overline{\unicode[STIX]{x1D731}}^{\dot{v}}$
is the particle phase average of
$\unicode[STIX]{x1D731}$
(see table 3 for definitions, all the variables with a dot are linked to the point-size particle). For
$\unicode[STIX]{x1D6F7}=\boldsymbol{u}_{p}$
(the particle velocity), the right-hand side of (A 3) is given by 
. Then, with the decomposition
$\boldsymbol{u}_{p}=\overline{\boldsymbol{u}_{p}}^{\dot{v}}+\boldsymbol{u}_{p}^{\prime }$
, equation (A 3) gives the averaged momentum equation of the ‘particle phase’ (A 4) from a Eulerian viewpoint
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{p}}^{\dot{v}}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{p}}^{\dot{v}}\overline{\boldsymbol{u}_{p}}^{\dot{v}})=\underbrace{\underbrace{-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{p}^{\prime }\boldsymbol{u}_{p}^{\prime }}^{\dot{v}})}_{\boldsymbol{M}^{\text{Re}}}+\underbrace{\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\boldsymbol{g}}_{\boldsymbol{M}^{\boldsymbol{g}}}+\underbrace{\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\left(\frac{\boldsymbol{F}_{\boldsymbol{v}}^{\boldsymbol{S}_{\boldsymbol{b}}}}{m_{p}}\right)}^{\dot{v}}}_{\boldsymbol{M}^{\boldsymbol{S}_{\boldsymbol{b}}}}}_{\boldsymbol{M}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}}, & & \displaystyle\end{eqnarray}$$
where,
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\text{D}\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k}}{\text{D}t}} & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k}\overline{\boldsymbol{u}_{k}}^{k}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\left[{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\boldsymbol{u}_{k}}^{k}}{\unicode[STIX]{x2202}t}}+\overline{\boldsymbol{u}_{k}}^{k}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{\boldsymbol{u}_{k}}^{k}\right]+\overline{\boldsymbol{u}_{k}}^{k}\left[\underbrace{{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\overline{\boldsymbol{u}_{k}}^{k})}_{=0}\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & = & \displaystyle \unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D70C}_{k}\left[{\displaystyle \frac{\unicode[STIX]{x2202}\overline{\boldsymbol{u}_{k}}^{k}}{\unicode[STIX]{x2202}t}}+\overline{\boldsymbol{u}_{k}}^{k}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\overline{\boldsymbol{u}_{k}}^{k}\right]=\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{k}\frac{\text{d}\overline{\boldsymbol{u}_{k}}^{k}}{\text{d}t},\end{eqnarray}$$
with
$\text{d}/\text{d}t$
the material derivative. Hence, equation (A 4) is a trajectory equation of a set of particles in the Eulerian formalism, comparable to Newton’s second law;
$\boldsymbol{M}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}$
is the sum of external forces applied on particles. At the statistical steady state,
$\boldsymbol{M}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}=0$
. Equation (A 4) also looks like the averaged Navier–Stokes equations for the vapour phase without pressure and viscous terms. Because of the particle approach, the viscous term inside the ‘particle’ does not exist. To approach the Navier–Stokes equations for the vapour phase (A 18), a pressure gradient can be extracted from the surface force
$\boldsymbol{F}^{\boldsymbol{S}_{\boldsymbol{b}}}$
(see Gatignol Reference Gatignol1983).
Then, equation (A 4) can be reformulated to introduce a pressure gradient,
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{p}}^{\dot{v}}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{p}}^{\dot{v}}\overline{\boldsymbol{u}_{p}}^{\dot{v}}) & = & \displaystyle -\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}}^{\dot{l}}-\widetilde{p_{l}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{p}^{\prime }\boldsymbol{u}_{p}^{\prime }}^{\dot{v}})+\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D70C}_{v}\boldsymbol{g}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{\boldsymbol{M}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}-\boldsymbol{M}^{\boldsymbol{g}}-\boldsymbol{M}^{\text{Re}}+\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}}^{\dot{l}}-\widetilde{p_{l}})}_{\dot{\boldsymbol{M}}_{v}}.\end{eqnarray}$$
Equation (A 8) now looks like the continuous Navier–Stokes equations of the vapour phase with particle quantities (see (A 18) in § A.2 for a comparison between the two formulations). In the two-fluid formulation,
$\dot{\boldsymbol{M}}_{v}$
is often named the ‘total interfacial force’. Previous equations show that this current appellation is misleading. The real ‘total interfacial force’ is given by the right-hand side of the trajectory equation (A 4). Thus,
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}=\dot{\boldsymbol{M}}_{v}+\boldsymbol{M}^{\boldsymbol{g}}+\boldsymbol{M}^{\text{Re}}-\dot{\unicode[STIX]{x1D6FC}}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}}^{\dot{l}}-\widetilde{p_{l}}), & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}$
is the total interfacial force and
$\dot{\boldsymbol{M}}_{v}$
is the interfacial momentum transfer term between phases. The paper proposes a clarification of the content of
$\boldsymbol{M}_{v}^{\boldsymbol{t}\boldsymbol{o}\boldsymbol{t}}$
and
$\dot{\boldsymbol{M}}_{v}$
.
Table 3. Definitions of particle and continuous functions and phase averages;
$\unicode[STIX]{x1D6FF}$
is a Dirac impulsion and
$H$
is the Heaviside function.
A.2 Particle versus continuous equations
From the local balance equations, the two-phase exact formulation of the vapour and liquid phases momentum equations without phase change is (Delhaye Reference Delhaye2008),
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}\overline{\boldsymbol{u}_{v}}^{v}) & = & \displaystyle -\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\widetilde{p_{l}})]-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\boldsymbol{g}+\boldsymbol{M}_{v}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70F}_{v}}^{v}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}\overline{\boldsymbol{u}_{l}}^{l}) & = & \displaystyle -\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{l}(\overline{p_{l}}^{l}-\widetilde{p_{l}})]-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\boldsymbol{g}+\boldsymbol{M}_{l}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\overline{\unicode[STIX]{x1D70F}_{l}}^{l}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}_{k}=\overline{(-p_{k}\boldsymbol{I}+\unicode[STIX]{x1D70F}_{k})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{k}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}_{l}+\boldsymbol{M}_{v}=\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}, & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70F}_{k}=\unicode[STIX]{x1D707}_{k}(\unicode[STIX]{x1D735}\boldsymbol{u}_{k}+\unicode[STIX]{x1D735}^{\text{T}}\boldsymbol{u}_{k})$
and
$\widetilde{p_{l}}$
is the reference pressure. In order to compare the formulation (A 8) (written for particles) with the continuous and exact two-phase formulation (A 10), the formal link between the particle and continuous viewpoints has to be established. Even if the relation between the two columns of table 3 is not obvious, a Taylor expansion is proposed by Lhuillier, Morel & Delhaye (Reference Lhuillier, Morel and Delhaye2000). It can be written for any volume variable
$\unicode[STIX]{x1D719}$
if this variable has a typical length scale bigger than the diameter of a bubble,
$$\begin{eqnarray}\displaystyle \overline{\unicode[STIX]{x1D712}_{k}\unicode[STIX]{x1D719}}=\overline{\dot{\unicode[STIX]{x1D712}}_{k}\int \unicode[STIX]{x1D719}dv}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\overline{\dot{\unicode[STIX]{x1D712}}_{k}\int \boldsymbol{r}\unicode[STIX]{x1D719}dv}\right)+\frac{1}{2}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left[\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\overline{\dot{\unicode[STIX]{x1D712}}_{k}\int \boldsymbol{r}\boldsymbol{r}\unicode[STIX]{x1D719}dv})\right]-\cdots \,. & & \displaystyle\end{eqnarray}$$
With this relation we get, at first order,
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D6FC}}_{k}=\unicode[STIX]{x1D6FC}_{k}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D719}}^{\dot{v}}=\overline{\unicode[STIX]{x1D719}}^{v}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{\unicode[STIX]{x1D719}}^{\dot{l}}=\overline{\unicode[STIX]{x1D719}}^{l}. & \displaystyle\end{eqnarray}$$
Then, expressions (A 15), (A 16) and (A 17) applied to (A 8) give, at first order, the two-fluid equation written for the vapour phase (A 18),
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}\overline{\boldsymbol{u}_{v}}^{v}) & = & \displaystyle -\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}}^{l}-\widetilde{p_{l}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\boldsymbol{g}+\dot{\boldsymbol{M}}_{v}.\end{eqnarray}$$
The comparison between (A 18) and (A 10) confirms the physical meaning of the particle hypothesis and allows us to express
$\boldsymbol{M}_{v}$
as a function of
$\dot{\boldsymbol{M}}_{v}$
(the viscous effect inside the bubbles
$\unicode[STIX]{x1D6FC}_{v}\overline{\unicode[STIX]{x1D70F}_{v}}^{v}$
is neglected). The relation between the two viewpoints is then, at first order, by assimilation between (A 18) and (A 10),
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{v}=\dot{\boldsymbol{M}}_{v}+(\overline{p_{l}}^{l}-\widetilde{p_{l}})\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}+\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})]. & & \displaystyle\end{eqnarray}$$
Taking (A 19) in the averaged Navier–Stokes equations (A 11) on the liquid phase gives with (A 13),
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}\overline{\boldsymbol{u}_{l}}^{l}) & = & \displaystyle -\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}(\overline{p_{l}}^{l}-\widetilde{p_{l}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}+\unicode[STIX]{x1D6FC}_{l}\overline{\unicode[STIX]{x1D70F}_{l}}^{l})+\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\boldsymbol{g}\nonumber\\ \displaystyle & & \displaystyle \underbrace{-\dot{\boldsymbol{M}}_{v}+\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})]}_{\dot{\boldsymbol{M}}_{l}}.\end{eqnarray}$$
In addition to classical terms of the two-fluid formulation, equation (A 20) presents a surface tension part
$\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}$
and a pressure term
$-\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})]$
. They are essential to our analysis; they are at the origin of the laminar dispersion force introduced in this work. This formulation tends to the classical relation
$\dot{\boldsymbol{M}}_{l}=-\dot{\boldsymbol{M}}_{v}$
assuming that the interfacial jump condition in the momentum equation is led, at first order and for averaged quantities, by the pressure jump only in a form of the averaged Laplace law
$\overline{p_{v}}^{v}-\overline{p_{l}}^{l}=\unicode[STIX]{x1D70E}\overline{\unicode[STIX]{x1D705}}$
. In this case, the pressure jump at the interface compensates exactly the surface tension term and
$\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})]=0$
.
In practice, for RANS Euler–Euler computations, only the undisturbed liquid pressure in the absence of bubbles
$\overline{p_{l}^{SP}}^{l}$
is solved. This means that the part of the liquid pressure induced by the surface tension
$\overline{p_{l}^{b}}^{l}$
is considered as a momentum source term responsible for interfacial forces. Thus, equations (A 18) and (A 20) become,
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}}^{v}\overline{\boldsymbol{u}_{v}}^{v}) & = & \displaystyle -\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\overline{\boldsymbol{u}_{v}^{\prime }\boldsymbol{u}_{v}^{\prime }}^{v})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{v}\boldsymbol{g}+\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}}^{l}\overline{\boldsymbol{u}_{l}}^{l}) & = & \displaystyle -\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}(\overline{p_{l}^{SP}}^{l}-\widetilde{p_{l}})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\overline{\boldsymbol{u}_{l}^{\prime }\boldsymbol{u}_{l}^{\prime }}^{l}-\unicode[STIX]{x1D6FC}_{l}\overline{\unicode[STIX]{x1D70F}_{l}}^{l})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D70C}_{l}\boldsymbol{g}+\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=-\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}+\overline{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D712}_{v}}-\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})]-\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}}), & & \displaystyle\end{eqnarray}$$
$\boldsymbol{M}_{k}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}=\dot{\boldsymbol{M}_{\boldsymbol{k}}}-\unicode[STIX]{x1D6FC}_{k}\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}})$
. Thus, the relations between the DNS and the RANS viewpoints are given by (A 19), 
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}_{v}=\boldsymbol{M}_{v}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}+\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}})+(\overline{p_{l}}^{l}-\widetilde{p_{l}})\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}+\unicode[STIX]{x1D735}[\unicode[STIX]{x1D6FC}_{v}(\overline{p_{v}}^{v}-\overline{p_{l}}^{l})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{M}_{l}=\boldsymbol{M}_{l}^{\boldsymbol{R}\boldsymbol{A}\boldsymbol{N}\boldsymbol{S}}+\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D735}(\overline{p_{l}^{b}}^{l}-\widetilde{p_{l}})-(\overline{p_{l}}^{l}-\widetilde{p_{l}})\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}_{v}. & \displaystyle\end{eqnarray}$$
