Sub-grid Scale Modelling and a-Posteriori Tests with a Morphology Adaptive Multifield Two-Fluid Model Considering Rising Gas Bubbles

Abstract

The predictive simulation of gas–liquid multiphase flows at industrial scales reveals the challenging task to consider turbulence and interfacial structures, which span a large range of length scales. For simulation of relevant applications, a hybrid model can be utilised, which combines the Euler–Euler model for the description of small interfacial structures with a volume-of-fluid model as a scale-resolving multiphase approach. Such a hybrid model needs to be able to simulate interfaces, which are hardly resolved on a coarse numerical grid. The goal of this work is to improve the prediction of interfacial gas–liquid flows on a numerical grid with comparably large grid spacing. From the low-pass filtering of the two-fluid model five unclosed sub-grid scale terms arise. The convective and the surface tension part of the aforementioned contributions are individually modelled with multiple closure formulations. Those models are a-posteriori assessed in cases of two- and three-dimensional gas bubbles rising in stagnant liquid. It is shown, that the chosen closure modelling approach is suitable to improve the predictive power of the numerical model utilised in this work. Hence, simulation results on comparably coarse grids are changed towards results obtained with higher spatial resolution.

1 Introduction

In a large number of applications and processes, e.g., in renewable energy, mining, processing or nuclear power industry, gas–liquid multiphase flows have a major impact on system performance and reliability. In particular, interactions between interfacial and turbulent dynamics in the contact region between individual immiscible phases are of great importance as they significantly influence the behaviour of the whole system. In order to develop and improve these processes in terms of safety and efficiency, a way for reliably predicting actual physical occurrences at industrial scales with numerical tools is urgently needed. The choice of a method has always been a trade-off between the claim for precision and computational expenses. In recent years a strong trend towards hybrid multiphase models is observed, where interface resolving approaches are combined with ones, that describe interfacial structures in a statistical manner. The long-term objective is to simulate both large and small features of interfaces appropriately and to switch between different types of flow description depending on the actual type of flow. For this purpose, Hänsch et al. (2012) proposed a hybrid model, which involves an algebraic volume-of-fluid-like approach and the Euler–Euler model for the description of large-scale and small-scale interfacial structures, respectively. This concept is entirely formulated in the spirit of the multifield two-fluid model, which forms the mathematical basis for the Euler–Euler model (Drew and Passman 1999). Under the condition of strong interfacial coupling, the two-fluid model behaves like the one-fluid model by equalising all phase-specific velocity fields (Yan and Che 2010). Meller et al. (2021) apply compact momentum interpolation (Cubero et al. 2014) and interfacial drag coupling (Štrubelj and Tiselj 2011) to the hybrid model of Hänsch et al. (2012). This numerical framework forms the basis for the present work. It is noted that the term multifield two-fluid model expresses the applicability to multiphase flows rather than being restricted to two-phase flows (Meller et al. 2021). Even in the case of the presence of two distinct physical phases, different multiphase morphologies are handled by means of individual sets of equations, hence, the term multifield. Such situations are beyond the scope of the present work.

Both ways to describe multiphase flows—Eulerian–Eulerian (E–E) and algebraic volume-of-fluid modelling (VOF)—are based on different assumptions regarding the ratio of length scales of interfacial structures, e.g., representative bubble diameter $D_{\mathrm{b}}$, and spacing of the computational grid $\Delta _{x}$. In the first case, dispersed structures, such as bubbles, droplets or particles, are modelled to be much smaller compared to the size of a grid cell and vice-versa for the second case. Besides a schematic of both model descriptions, combined in the hybrid multiphase model, Fig. 1 shows their respective positions on the axis of length scale ratio $D_{\mathrm{b}}/\Delta _{x}$. It turns out, that there is a gap in the applicability of underlying assumptions of the basic multiphase methods in the range of $10^{0} \lesssim D_{\mathrm{b}}/\Delta _{x}\lesssim 10^{1}$. Aiming at a fully applicable hybrid multiphase model, a description of interfacial structures of every size, based on legitimate requirements, is crucial. The present work focuses on the representation of coarsely resolved interfacial and turbulence structures via algebraic VOF.

Fig. 1

Range of scales of interfacial structures in relation to spatial resolution for a hybrid multiphase model

For this particular purpose, the concept of large-eddy simulation (LES) has been expanded to multiphase flows in the context of the one-fluid model several times in the past: Most often spatial low-pass filtering of governing equations (Sagaut 2006) forms the theoretical basis for this type of simulations. This procedure gives rise to numerous different unclosed terms in the so-called filtered Navier–Stokes equations. The influence of those terms has been assessed with a-priori investigations (Labourasse et al. 2004, 2007; Toutant et al. 2006, 2008, 2009a, b; Vincent et al. 2008; Larocque et al. 2010; Chesnel et al. 2011; Liovic and Lakehal 2007a; Ketterl and Klein 2016, 2018; Klein et al. 2019; Mimouni et al. 2017; Hasslberger et al. 2020; Saeedipour et al. 2019a; Saeedipour and Schneiderbauer 2019). In only a few investigations model formulations for those unclosed terms have been applied in actual multiphase simulations for a-posteriori model assessment (Liovic and Lakehal 2007a, b, 2012; Saeedipour et al. 2019a, b; De Villiers et al. 2004; Aniszewski et al. 2012; Herrmann 2013, 2015; Ketterl et al. 2019). All of the previously listed investigations have been carried out in the context of the one-fluid formulation. In contrast, Fleau (2017) and Mimouni et al. (2017) adopted the spatial low-pass filter formalism to multifield two-fluid models. Schneiderbauer (2017), Cloete et al. (2018) and Sarkar et al. (2016) apply a similar approach of filtering the two-fluid model to gas-sold flows. In those investigations, the focus is on the sub-grid interfacial drag force in the frame of dispersed flows.

In the present work different closure models for sub-grid scale terms are applied to gas–liquid multiphase flows considering single rising bubbles with the hybrid multiphase model mentioned above (Meller et al. 2021). The focus of the present work is to create a basis for sub-grid scale modelling in the context of the filtered multiphase two-fluid model equations. It is the goal to improve the prediction of gas–liquid interfacial flows with coarse computational grids. For this purpose, the basic equations and the modelling strategy are given in Sect. 2 including closure formulations for unclosed convective and surface tension sub-grid scale (SGS) terms being reviewed and adapted to the present framework. The numerical method is briefly discussed in Sect. 3. Simulation results of applied SGS closure models together with detailed discussions are presented in Sect. 4 considering two- and three-dimensional rising gas bubbles in stagnant liquid. The work is finally concluded in Sect. 5.

2 Modelling

2.1 Basic Equations

The hybrid multiphase model is based on the multifield two-fluid model equations (Drew and Passman 1999). The methodology of spatial low-pass filtering is applied to this equation system, which is realised with a convolution integral (Labourasse et al. 2007; Sagaut and Germano 2005):

$$\begin{aligned} \overline{\varphi }({\mathbf{x}}, t) = G * \varphi = \int \varphi ({\mathbf{y}}, t)\,G({\mathbf{x}} - {\mathbf{y})}\,{\mathrm{d}}{\mathbf{y}}. \end{aligned}$$

(1)

The variable $\varphi$ represents a general scalar, vector or tensor field quantity. Mimouni et al. (2017) and Fleau (2017) present the application of this formalism to the multifield two-fluid model, which results in the filtered multifield two-fluid model equations. Under the assumption, that density and molecular viscosity are constant for each individual phase $\alpha$, this includes conservation equations of phase volume fraction $r_\alpha$

$$\begin{aligned} \frac{\partial \overline{r}_\alpha }{\partial t} + \frac{\partial }{\partial x_{j}}\left( \overline{r}_\alpha \,\overline{u}_{\alpha ,j}\right) = - \frac{\partial \tau _{ru,\alpha ,j}}{\partial x_{j}} \end{aligned}$$

(2)

and momentum

$$\begin{aligned} \begin{aligned}&\rho _\alpha \left( \frac{\partial }{\partial t}\left( \overline{r}_\alpha \overline{u}_{\alpha ,i}\right) + \frac{\partial }{\partial x_{j}}\left( \overline{r}_\alpha \overline{u}_{\alpha ,i} \overline{u}_{\alpha ,j}\right) \right) = - \overline{r}_\alpha \frac{\partial \overline{p}}{\partial x_{i}} + \frac{\partial }{\partial x_{j}}\left( \overline{r}_\alpha 2 \mu _\alpha \overline{S}_{\alpha ,ij} \right) \\&\quad + \overline{r}_\alpha \rho _\alpha g_{i} + \overline{r}_\alpha \sum _{\beta \ne \alpha }\overline{\sigma }_{\alpha \beta }\overline{\delta }|^{\mathrm{S}}_{\alpha \beta }\overline{\kappa }|^{\mathrm{S}}_{\alpha \beta }\overline{n}|^{\mathrm{S}}_{\alpha \beta ,i} + \overline{f}_{\alpha ,i}^{\mathrm{D}}\\&\quad - \rho _\alpha \frac{\partial }{\partial t}\left( \tau _{ru,\alpha ,i}\right) + \frac{\partial }{\partial x_{j}}\left( - \rho _\alpha \,\tau _{ruu,\alpha ,ij} + 2 \mu _\alpha \tau _{rS,\alpha ,ij} \right) \\&\quad - \tau _{rp,\alpha ,i} + \sum _{\beta \ne \alpha }\tau _{rnn,\alpha \beta ,i}. \end{aligned} \end{aligned}$$

(3)

Linearity of the filter operation as well as a negligible commutation error of derivation and filtering are assumed. Einstein summation convention applies wherever two identical indices occur in a single term. Phase-specific density, molecular viscosity and velocity are denoted with $\rho _\alpha$, $\mu _\alpha$ and ${\mathbf{u}}_\alpha$, respectively. Pressure p is shared between all phases. This assumption is fundamental to the formulation of a single pressure equation. This differential equation is derived from mass conservation of the fluid mixture and phase-specific momentum equations, which are weighted with the phase volume fraction $r_\alpha$ (Ferziger et al. 2002). The gravitational force is ${\mathbf{g}}$. Surface tension force is modelled as a continuum surface force according to Brackbill et al. (1992). For a pair of phases $\alpha$ and $\beta$ the symbols $\sigma _{\alpha \beta }$, $\delta |^{\mathrm{S}}_{\alpha \beta }$, $\kappa |^{\mathrm{S}}_{\alpha \beta }$ and ${\mathbf{n}}|^{\mathrm{S}}_{\alpha \beta }$ denote the surface tension coefficient, interface-Dirac function, interface curvature and interface normal vector, respectively. The interfacial drag force acting on phase $\alpha$ is denoted by ${\mathbf{f}}_\alpha ^{\mathrm{D}}$. The filtered strain rate tensor is defined as

$$\begin{aligned} \overline{S}_{\alpha ,ij} = \frac{1}{2}\left( \frac{\partial \overline{u}_{\alpha ,i}}{\partial x_{j}} + \frac{\partial \overline{u}_{\alpha ,j}}{\partial x_{i}} \right) - \delta _{ij} \frac{1}{3}\frac{\partial \overline{u}_{\alpha ,k}}{\partial x_{k}}. \end{aligned}$$

(4)

Due to the filtering procedure sub-grid scale (SGS) contributions appear. Those terms are similarly named as in the work of Ketterl and Klein (2018), who applied the spatial low-pass filtering to the one-fluid model:

$$\begin{aligned} \begin{aligned}&{\text { SGS convection:}}&\tau _{ruu,\alpha ,ij} =&\ \overline{r_\alpha u_{\alpha ,i} u_{\alpha ,j}} - \overline{r}_\alpha \overline{u}_{\alpha ,i} \overline{u}_{\alpha ,j}, \\&{\text { SGS diffusion:}}&\tau _{rS,\alpha ,ij} =&\ \overline{r_\alpha S_{\alpha ,ij}} - \overline{r}_\alpha \overline{S}_{\alpha ,ij}, \\& &\ &\ {\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!- \overline{r}_\alpha \overline{\sigma }_{\alpha \beta }\overline{\delta }_{\alpha \beta }\overline{\kappa }|^{\mathrm{S}}_{\alpha \beta }\overline{n}|^{\mathrm{S}}_{\alpha \beta ,i}}, \\&{\text { SGS surface tension:}}&\tau _{rnn,\alpha \beta ,i} =&\ \overline{r_\alpha \sigma _{\alpha \beta }\delta _{\alpha \beta }\kappa |^{\mathrm{S}}_{\alpha \beta }n|^{\mathrm{S}}_{\alpha \beta ,i}} \\& {\text { SGS interfacial term:}}&\tau _{ru,\alpha ,i} =&\ \overline{r_\alpha u_{\alpha ,i}} - \overline{r}_\alpha \overline{u}_{\alpha ,i}, \\&{\text { SGS pressure:}}&\tau _{rp,\alpha ,i} =&\ \overline{r_\alpha \frac{\partial p}{\partial x_{i}}} - \overline{r}_\alpha \frac{\partial \overline{p}}{\partial x_{i}}. \end{aligned} \end{aligned}$$

(5)

It is worth noting, that each individual one of the five SGS terms is affected by the filtered gradient of phase fraction $r_\alpha$, which appears in the interface region. Therefore, all SGS contributions are to be expected in two-phase flows. While the classical SGS stress, as known from single-phase flows, is purely sensitive to gradients in the velocity field, the corresponding SGS convection term additionally considers filtered phase fraction gradients. Hence, an SGS convection contribution is not limited to turbulent flows but appears in (under-resolved) laminar two-phase flows as well.

In regimes of large-scale interfacial structures, different phases must not interpenetrate each other. In order to reproduce such dynamics with the underlying multifield two-fluid model, an interfacial no-slip condition is enforced via the drag model formulation of Štrubelj and Tiselj (2011):

$$\begin{aligned} \overline{{\mathbf{f}}}_{\alpha \beta }^{\mathrm{D}}= K_{\alpha \beta }^{\mathrm{D}}(\overline{\mathbf{u}}_\beta - \overline{\mathbf{u}}_\alpha ) \quad {\text {with}}\quad K_{\alpha \beta }^{\mathrm{D}}= \overline{r}_\alpha \overline{r}_\beta \rho _{\alpha \beta }\frac{1}{\tau _{r}} \end{aligned}$$

(6)

with mixture density

$$\begin{aligned} \rho _{\alpha \beta }= \frac{\rho _\alpha \,\overline{r}_\alpha + \rho _\beta \,\overline{r}_\beta }{\overline{r}_\alpha + \overline{r}_\beta }. \end{aligned}$$

(7)

Relaxation time $\tau _{r}$ is chosen to be several orders of magnitude smaller compared to the numerical time step size $\Delta _{t}$. If the interfacial drag coupling is very strong, the behaviour of the one-fluid model will be recovered (Yan and Che 2010; Meller et al. 2021) equalising all fields of phase-specific velocities. Therefore, index $\alpha$ specifying the phase of a phase-specific velocity will be skipped, when investigating velocity fields in simulation results. The fluid dynamics are mathematically described within the two-fluid model at any time. Due to the drag coupling described above, both individual momentum equations effectively collapse onto each other delivering an identical behaviour, which is identical to the homogeneous model. Hence, the numerical method can be described as an algebraic volume-of-fluid-like method.

2.2 Sub-grid Scale Modelling

The aforementioned five unclosed SGS terms are present for each individual phase $\alpha$. Via an a-priori model analysis considering primary atomisation of a liquid jet, Ketterl and Klein (2016) estimated the contribution of SGS convection to have the strongest influence among all present SGS terms. However, the absolute value for the resulting force of the SGS surface tension contribution is comparably small. According to Herrmann and Gorokhovski (2008), the contribution of this unclosed term might nevertheless be of great importance due to its dependence on interface curvature. For these reasons, the present work focuses on the investigation of closure models for both SGS convection and the SGS surface tension.

2.2.1 Convective Sub-grid Scale Term

This unclosed term is object of numerous investigations since many decades and different modelling strategies have proven to be reasonable, such as functional and structural modelling. All model formulations addressed in this section are originally formulated in the context of the one-fluid model. They are adapted to the multifield two-fluid formulation by accounting for the filtered phase volume fraction $\overline{r}_\alpha$. Some convective SGS closure models rely on the filter length $\overline{\Delta }$. In the course of this work, this length scale is estimated to be $\overline{\Delta }\approx \sqrt{A_{\mathrm{CV}}}$ or $\overline{\Delta }\approx \root 3 \of {V_{\mathrm{CV}}}$ in two- and three-dimensional space, respectively. The surface area of a control volume is referred to as $A_{\mathrm{CV}}$ and $V_{\mathrm{CV}}$ denotes its volume.

Functional models for the closure of the convective SGS term usually attempt to mimic contributions of non-resolved exchange of momentum by increasing viscosity. The most prominent SGS model of such type is the one proposed by Smagorinsky (1963), hereafter denoted as SMAG. Another SGS closure model of eddy-viscosity type is the sigma model (SIG) proposed by Nicoud et al. (2011). Due to its ability to generate zero eddy-viscosity for two-dimensional or two-component flows, as well as for axisymmetric and isotropic compressions/expansions it delivers the proper cubic behaviour in near-wall regions (Nicoud et al. 2011). Formulations for all individual closure models, which are subject to this section, are listed in Table 1.

Table 1 Overview over models for the convective sub-grid scale term $\varvec{\tau }_{ruu,\alpha }$

In contrast, structural models aim at the estimation of all individual elements of the convective SGS stress tensor. Clark et al. (1979) proposed the gradient model (GM), which is based on the idea, that SGS structures of the velocity field may be approximated by unfolding the filtered field via the leading term of a Taylor-expansion. The model has been extended by Ketterl and Klein (2018) to take varying fluid properties due to spatial distribution of multiple phases into account and therefore is called extended gradient model (EGM). Approximation of SGS structures via a test filter operation $\check{(\cdot )}$ has been proposed by Liu et al. (1994) and, as it is based on the assumption of scale-similarity (Bardina et al. 1980), is denoted as scale-similarity model (SSM). The test filter is realised as the so-called simple filter (The OpenFOAM Foundation Ltd. 2020), which is constructed as a linear interpolation of values from cell centres to cell faces and back via a surface integral, weighted with the surface area of individual cell faces.

Contrary to most functional models, which increase effective viscosity and therefore act dissipatively, the proposed structural models may allow for anti-dissipative behaviour (Bardina et al. 1980). Although this might correctly describe physical processes, it may be harmful to the stability of the numerical procedure (Anderson and Domaradzki 2012). In order to counteract this drawback, Bardina et al. (1980) linearly combined a structural model with a functional model, SMAG in this case. Hereafter, resulting model combinations are called mixed model (M-mod-SMAG). Generally, mixed models can be a combination of structural models with different individual eddy-viscosity type models, which would also allow for a combination with SIG. It turns out that the models SMAG and SIG reveal nearly identical results in the investigations carried out in the course of the present work. In order to limit the number of model combinations, mixed models are purely formulated with SMAG as eddy-viscosity model formulation.

Apart from linear combination with functional models, Kobayashi (2018) proposed a regularisation procedure, which is based on the idea to use arbitrary structural models $\varvec{\tau }_{ruu,\alpha }^{\mathrm{mod}}$, while preventing anti-diffusion. This modelling approach is modified by Klein et al. (2020) to be applicable without solving a transport equation for turbulent kinetic energy and is referred to as KOB. Additionally, Klein et al. (2020) proposed a simplification of the regularisation originating from the concept of adding just the right amount of dissipation in case the structural model itself delivers anti-dissipation. In this way the net energy transfer is set to zero, whenever anti-dissipation is detected. The latter formulation is denoted with KL.

2.2.2 Surface Tension Sub-grid Scale Term

Analogous to the convective SGS stress, functional as well as structural modelling approaches have been proposed for approximation of the force resulting from SGS surface tension. The application of the scale similarity concept to this unclosed term is formulated in terms of explicit volume filtering. In contrast to that, the test filter could be applied based on the concept of surface filtering as well (Hasslberger et al. 2020). However, this work focuses on volume-based filtering for the test filter operation. The resulting model formulation is proposed by Ketterl and Klein (2018) and is referred to as nn-SSM. All SGS surface tension closure formulations are listed in Table 2.

Table 2 Overview over models for surface tension sub-grid scale term $\varvec{\tau }_{rnn,\alpha \beta }$

Most functional modelling approaches for closure of the convective SGS stress are based on the argumentation, that on a computational grid with given spatial resolution only a limited effective Reynolds number can be reliably reproduced. Therefore, the viscosity is artificially increased by means of a modelled turbulent viscosity. Ketterl and Klein (2018) adapted this formalism to SGS surface tension via the concept of a maximum effective Weber-number $We _{\mathrm{crit}}$ below which an interfacial structure is not expected to breakup. The condition of a maximum effective Weber-number is ensured via an additional surface tension coefficient $\sigma _{\mathrm{SGS},\alpha \beta }$, which relies on the volume fraction weighted velocity $\overline{\mathbf{u}}_{\alpha \beta }$ of phase pair $\alpha \beta$. This quantity is determined analogously to Eq. (7). The quantity $\sigma _{\mathrm{SGS},\alpha \beta }$ is based on the definition of a fixed critical Weber-number, which was originally assumed to be $We _{\mathrm{crit}}=10$ for liquid droplets in gas (Ketterl and Klein 2018). For gas bubbles in liquid Miksis (1981) proposed a value of $We _{\mathrm{crit}}=2.76$. As it will turn out later, this requirement is too weak to deliver any significant model influence in the cases under investigation. Hence, smaller values will be assumed for $We _{\mathrm{crit}}$ as well. In favour of numerical stability, the Weber-number model is reformulated in such a way, that

• 1. The formulation is based on the mixture density $\rho _{\alpha \beta }$ of a given phase pair $\alpha \beta$ (see Eq. (7)), and

• 2. The effective surface tension coefficient is defined to be the maximum of its molecular value and the SGS model value: $\sigma _{\mathrm{SGS},\alpha \beta }^{\mathrm{eff}}=\max (\sigma _{\alpha \beta }, \sigma _{\mathrm{SGS},\alpha \beta }^{ We -{\mathrm{mod}}})$.

Hereafter, the Weber-number model in connection with a fixed critical Weber-number is referred to as We-KE.

The critical Weber-number is a model parameter and appropriate values may vary among different cases. Therefore, a formulation, which is distinct from a fixed value is highly desirable. For this purpose, the functional relationship between Reynolds-number and critical Weber-number by Ryskin and Leal (1984) is proposed here to be used in conjunction with the general approach of the Weber-number model. The local Reynolds-number $Re _{\mathrm{loc}}$ relies on the volume fraction weighted kinematic viscosity $\overline{\nu }_{\alpha \beta }$ of phase pair $\alpha \beta$, which is also calculated from phase-specific quantities analogously to Eq. (7). The resulting model formulation is denoted with We-Re.

3 Numerical Method

The system of differential equations is solved via a finite-volume method on an unstructured computational grid, where solution variables are stored in the centre position of each grid cell. A segregated solution procedure is used including a projection method for pressure-velocity coupling. All terms are spatially discretised with second order accuracy including flux limiting schemes (Hirsch 1990) for the convective terms of both phase fraction and phase-specific momentum transport equations. The transient terms are discretised first order explicitly and implicitly in phase fraction equations and in phase-specific momentum transport equations, respectively. Thus, the overall solution procedure is characterised by semi-implicit discretisation in time.

With an algebraic volume-of-fluid-like method as it is utilised in this work, gradients of fields, e.g., of phase fraction $r_\alpha$, are by definition finite in regions of resolved interfaces. Nevertheless, magnitudes of gradients of field quantities may become very large. In order to accurately capture this in a simulation, while maintaining numerical stability the compact momentum interpolation (CMI), which is proposed by Cubero et al. (2014) is applied. From interfacial drag coupling via the drag model, described in Sect. 2.1 results a strong coupling between momentum equations of individual phases and, hence, a high stiffness of the resulting coupled system of equations. For segregated solution algorithms this implies a bad convergence rate over iterations. In order to overcome this drawback, the partial elimination algorithm (Spalding 1981) is employed.

The implementation is realised in the framework of the multiphaseEulerFoam solver, which is part of the C++ software library OpenFOAM (The OpenFOAM Foundation Ltd. 2020). Details on the numerical model as well as the public source code are provided by Meller et al. (2021) and Schlegel et al. (2021), respectively.

4 Simulation Results

The goal is to identify SGS closure models, which allow for improved predictions of interfacial flows on coarse numerical grids. Hence, the obtained simulation results shall ideally be similar to ones, which are obtained on grids with higher spatial resolution.

Both test cases considered in this work are characterised by the liquid phase being stagnant, which results in a zero turbulence intensity of the liquid flow approaching the gas bubbles. Hence, the SGS closure models react to the implicitly filtered interfacial region and the attached velocity boundary layers rather than turbulent eddies of chaotic nature. The numerical setups for both test cases are provided by Hänsch et al. (2021) and can be executed with the public source code (Schlegel et al. 2021).

4.1 Convective Sub-grid Scale Term

4.1.1 Three-Dimensional Rising Gas Bubble

The test case under investigation considers a single three-dimensional gas bubble, which rises in stagnant liquid. Material properties and constraints are selected, such that a wobbling bubble regime is achieved. The system is characterised by several dimensionless numbers: Reynolds number $Re _{g} = \rho _{\mathrm{L}}U_{g} D_{\mathrm{b}}/\mu _{\mathrm{L}}$, Eötvös number $Eo = \rho _{\mathrm{L}}\left| \mathbf{g}\right| D_{\mathrm{b}}^{2}/\sigma _{\mathrm{GL}}$ and Morton number $Mo =\left| \mathbf{g}\right| \mu _{\mathrm{L}}^{4} (\rho _{\mathrm{L}}- \rho _{\mathrm{G}})/(\rho _{\mathrm{L}}^{2} \sigma _{\mathrm{GL}}^{3})$ as well as by ratios of phase-specific density and kinematic viscosity values. Material properties specific for liquid and gas phases are denoted with L and G, respectively, whereas index GL indicates the pair of both phases. The gas bubble is initialised as a sphere of diameter $D_{\mathrm{b}}$. Furthermore, gravitational time $t_{g} = \sqrt{D_{\mathrm{b}}/\left| \mathbf{g}\right| }$ and gravitational velocity $U_{g} = \sqrt{\left| \mathbf{g}\right| D_{\mathrm{b}}}$ are defined. Quantities of length, time and velocity are made dimensionless via scaling with $D_{\mathrm{b}}$, $t_{g}$ and $U_{g}$, respectively. This operation is denoted by $\widehat{(\cdot )}$. The corresponding values of all characteristic dimensionless parameters for the test case under investigation are listed in Table 3. The cuboid computational domain of size $4D_{\mathrm{b}}\times 4D_{\mathrm{b}}\times 6D_{\mathrm{b}}$ is specified to be bounded by periodic conditions in both horizontal directions x and y. At the top of the domain, pure liquid enters with a uniform downward velocity $u_{\mathrm{inlet}}$. The bottom boundary is set up to serve as an outlet for pure liquid. The gas bubble is initialised with the centre of gravity at the target height, which is located in a distance of $2D_{\mathrm{b}}$ to the inlet boundary at the top. The deviation of the centre of gravity of the gas bubble from the target height is denoted with $\Delta _{Z_{\mathrm{b}}}$. Via a proportional-integral-derivative (PID) type controller the inlet velocity $u_{\mathrm{inlet}} = f_{\mathrm{PID}}(\Delta _{Z_{\mathrm{b}}}(t))$ is manipulated with the control value at time step n:

$$\begin{aligned} f_{\mathrm{PID}}(x_{n}) = K_{P}\,x_{n} + \sum \limits _{m=1}^{n}\,K_{I} \Delta _{t,m}\,x_{m} + K_{D} \frac{x_{n} - x_{n-1}}{\Delta _{t,n}}. \end{aligned}$$

(8)

The controller coefficients for proportional, integral and differential contributions are $K_{P}=-{0.8}\hbox { s}^{-1}$, $K_{I}=-{0.05}{\hbox { s}^{-2}}$ and $K_{D}=-{0.2}$, respectively. In this way, the bubble maintains its original vertical position, such that the aforementioned deviation is negligible. In that sense, the gas bubble is simulated in a frame of reference, which is attached to its vertical position, similarly to the approach of Fang et al. (2013). Therefore, inlet velocity $u_{\mathrm{inlet}}(t)$ is identical to the bubble rising velocity $U_{\mathrm{b}}(t)$. The centre of gravity itself is calculated via (Chen et al. 2004)

$$\begin{aligned} {\mathbf{x}}_{\mathrm{b}}= \begin{pmatrix}X_{\mathrm{b}}\\ Y_{\mathrm{b}}\\ Z_{\mathrm{b}}\end{pmatrix} = \frac{\int \limits _\Omega \overline{r}_{\mathrm{G}}{\mathbf{x}}\,{\mathrm{d}}V}{\int \limits _\Omega \overline{r}_{\mathrm{G}}\,{\mathrm{d}}V}. \end{aligned}$$

(9)

A schematic of the setup is shown in Fig. 2. In order to compare average field quantities, time averaging $\left\langle \cdot \right\rangle _{t}$ for the time period $14.38<\widehat{t}<\widehat{T}$ is applied. Additionally, field variables are averaged over circumferential angle $\Phi$ for each radial position r relative to the bubble’s centre of gravity (see Fig. 2):

$$\begin{aligned} \left\langle \varphi \right\rangle _\Phi = \frac{1}{2\pi }\int \limits _{0}^{2\pi }\varphi (\Phi )\,{\mathrm{d}}\Phi . \end{aligned}$$

(10)

Via an averaging procedure with respect to time t and angle $\Phi$, average radial profiles of vertical velocity component $\left\langle \overline{u_{z}}\right\rangle _{t,\Phi }(r)$ are obtained. This is equivalent to mapping the data to the lateral position of the bubble’s centre of gravity.

Table 3 Dimensionless parameters for test case of three-dimensional rising gas bubble

Fig. 2

Overview over test case of three-dimensional rising gas bubble with computational domain and boundary conditions; the gas bubble is shown in a position, which vertically deviates from this initial height (dashed line)

4.1.2 Sensitivity to Spatial Resolution

Prior to the assessment of model influences, the effect of spatial discretisation shall be investigated. For this purpose, four different levels of grid refinement are defined, namely G1 to G4. The corresponding quantities concerning the number of grid cells, grid spacing $\widehat{\Delta }_{x}$, time step size $\widehat{\Delta }_{t}$ and simulated time $\widehat{T}$ are listed in Table 4. In order to give an impression of bubble deformation and spatial resolution, Fig. 3 shows the bubble surface in front of the computational grid for refinement levels G1 to G4. First of all the focus will be on the trajectory of the centre of gravity. The corresponding curves in 3D space as well as their projection into the xy-plane are depicted in Fig. 4. Every time the gas bubble passes a pair of periodic boundaries, the trajectory is extended to show a continuous behaviour. In that way, discontinuities are prevented indicating the jump of the bubble centre of gravity from one side of the domain to the opposite one. Hence, the horizontal coordinates $\widehat{x}$ and $\widehat{y}$ used for the representation of the bubble trajectory may take values, which are formally beyond the bounds of the computational domain. When the gas bubble is cut by a pair of periodic boundaries, the centre of gravity is individually calculated for each gas structure and, subsequently, from those the common value is evaluated. Furthermore, the lengths of the different trajectories differ according to the individual maximum simulation times $\widehat{T}$ as listed in Table 4. It turns out, that the coarse spatial resolution of G1 leads to a zig-zag type rising path with a fixed orientation of the lateral movement. With grid level G2 the zig-zag path is still apparent but the orientation of lateral movements stays fixed only for short sections before it is slightly tiled in the z direction. With G3 the gas bubble shows a rather chaotic movement with short sections of the shape of a flattened helix, whereas a stable flattened helical rise is observed with G4. According to Cano-Lozano et al. (2016) the latter type of path is physical for the investigated bubble regime. Hence, the data obtained with G4 serve as a reference in the course of this work.

Table 4 Levels of grid refinement and corresponding characteristic numerical parameters

Fig. 3

Side view of the instantaneous bubble surface in front of computational grid for four different grid refinement levels

Fig. 4

Trajectories of rising gas bubble with different spatial resolutions as three-dimensional path and as projection into the plane $\widehat{z}=0$

Additionally to the trajectory, the bubble rising velocity $\widehat{U}_{\mathrm{b}}$ is an important characteristic. The corresponding time averaged values are listed in Table 5. As it turns out, a zig-zag type trajectory is connected to a lower rising velocity $\widehat{U}_{\mathrm{b}}$ compared to a chaotic path. This trend is visible for refinement levels G1 to G3. For the last refinement step from G3 to G4, the resulting averaged rising velocity decreases gently. Hence, the chaotic bubble movement reveals the highest bubble rising velocity $\widehat{U}_{\mathrm{b}}$ among the different types of trajectories.

Table 5 Temporally averaged dimensionless bubble rising velocity for different spatial resolution

The radial profile of vertical velocity component $\left\langle \widehat{\overline{u}}_{z}\right\rangle _{t,\Phi }$ allows for a more detailed look into the flow inside and outside the gas bubble as well as across the gas–liquid interface. Those distributions are shown in Fig. 5. A nearly uniform downward velocity is observed in the liquid flow in locations of large radial distance from the bubble centre of gravity. Those values differ slightly for the individual grid refinement levels, which corresponds to the different bubble rising velocities discussed above. Nevertheless, the absolute values are slightly larger compared to the average inflow velocity. This behaviour results from the fact, that a significant portion of the cross section of the computational domain is blocked by the gas bubble. In the interface region at $\widehat{r}\approx {0.5}$ the velocity increases towards the bubble centre of gravity and reaches a global maximum of positive value for all grid level corresponding to an internal gas flow, which is directed upwards. For G2 to G4 the value of the local velocity maximum reduces with increasing spatial resolution, while the radial velocity profile becomes flatter in this region. That implies, that the maximum upward velocity rises with larger grid spacing. This trend does not apply to G1 and it is assumed, that the grid spacing is just too coarse for an even sharper velocity peak than obtained with G2 to be resolved in the bubble centre. Furthermore, the radial velocity gradient rises in the interface region ($\widehat{r}\approx {0.5}$) with increasing refinement, at least for G1 to G3. The chaotic trajectory and the high rising velocity distinguish the result with refinement level G3 from the remaining ones. This behaviour is explained with the connection between bubble movement and an unique flow pattern inside and outside the gas bubble at this particular level of spatial resolution.

Fig. 5

Radial profiles of vertical velocity component for different mesh resolutions; the interface position for a volume-equivalent sphere is marked with a black solid vertical line

4.1.3 Model Assessment

In the following the individual SGS closure models presented in Sect. 2.2.1 are applied to the aforementioned test case. All a-posteriori investigations are carried out on spatial grid refinement level G2 and the time period for gathering transient and statistical data corresponds to simulation G2-O (see. Table 4). A comparison is drawn to the numerical results obtained with G2 and with G4 without application of any SGS closure model. The corresponding results are referred to as G2-O and G4-O, respectively.

The time averaged values of rising velocity $\widehat{U}_{\mathrm{b}}$ obtained with different convective SGS closure models obtained with G2 as well as the results without any model formulation (G2-O) are listed in Table 6. Those values are further expressed as relative deviation from the reference result G4-O in percent. It is observed, that the general trend of a lower rising velocity on G2 compared to G4 is maintained, regardless of the choice of closure model. The eddy-viscosity type models SMAG and SIG lead to slightly increased rising velocities with respect to G2-O. With the structural closure models the value of average rising velocity varies stronger compared to the functional model approaches. Among mixed and Klein-regularised models, M and KL, the gradient model (GM) and the extended gradient model (EGM) show small and moderate positive influence on the average rising velocity, respectively. However, the scale similarity model (SSM) results in a weaker bubble acceleration (KL-SSM) or even a deceleration (M-SSM-SMAG) compared to reference G2-O. Generally, structural models with Klein’s regularisation (KL) lead to more precise predictions of $\left\langle \widehat{U}_{\mathrm{b}}\right\rangle _{t}$ compared to the results without any model, with eddy-viscosity models and with the mixed approach. Among the three Kobayashi model configurations (KOB) the gas bubble rises fastest with the gradient model KOB-GM, which represents the most precise prediction. KOB-SSM leads to a slightly slower bubble rise, which is still an improvement compared to the result of G2-O. KOB-EGM decreases the value of the predicted velocity and, therefore, reveals a worse prediction in this context. The best velocity prediction among all closure models is achieved with the Klein-regularised extended gradient model (KL-EGM).

Table 6 Dimensionless rising velocity obtained with different closure models for convective SGS term $\varvec{\tau }_{ruu,\alpha }$ averaged over time in absolute values and as relative deviation from reference result of G4-O

In order to get an insight into the dynamical behaviour of the gas bubble, the trajectories predicted by the individual SGS closure models are projected into the horizontal plane. They are shown in Fig. 6 side-by-side together with the results without closure models. Again, discontinuities in the trajectories due to the gas bubble passing a pair of periodic boundaries are avoided as described in Sect. 4.1.1. Furthermore, it is noted, that the different periods of averaging for individual levels of grid refinement according to Table 4 result in trajectories of different length. The SGS models SMAG and SIG result in stable zig-zag bubble trajectories, which maintain their spatial orientation during the whole simulation time (Fig. 6e, f). Those results are similar to G1-O (Fig. 6a). Apparently, both eddy-viscosity models induce a dissipation by increasing the effective viscosity. This has a similar effect on the discretisation error in connection with a low spatial resolution. This behaviour is not desired. Instead an SGS closure model shall allow for simulation results comparable to a numerical grid, which is sufficiently fine to produce reliable results without any sub-grid scale model. Such a tendency is indeed observed for the investigated structural models in a sense, that the bubble rising path tends to be more chaotic compared to the results of G2-O (Fig. 6b), just as observed in G3-O (Fig. 6c). With all three mixed models as well as with KOB-EGM and KL-GM sections of stable zig-zag rising behaviour alternate with events of changing horizontal bubble position and path orientation. Looking at the model combinations, KOB-GM, KOB-SSM and KL-SSM the gas bubble tends to temporarily follow a straight diagonal rising path and from time to time changes the orientation of that trajectory. KL-EGM results in a bubble movement which comes closest to a chaotic change of rising direction among the investigated closure models. Only a single section of nearly straight motion is observed. From this point of view the latter model combination is evaluated to perform best regarding the shape of the trajectory. A flattened helical path as observed in G4-O (Fig. 6d) cannot be recovered with any of the investigated convective SGS closure models.

Fig. 6

Trajectories of rising gas bubble for different closure models for convective SGS term $\varvec{\tau }_{ruu,\alpha }$ projected into the horizontal plane

In order to gain more insight into flow inside and outside the gas bubble, radial profiles of averaged vertical velocity component $\left\langle \widehat{\overline{u}}_{z}\right\rangle _{t,\Phi }$ in the horizontal plane at the height of the bubble centre are shown in Figs. 7 and 8. Data resulting from application of closure models are compared to reference solutions G2-O and G4-O. It is evident from the radial profiles, that the eddy-viscosity models SMAG and SIG as well as all model combinations incorporating the gradient model GM show only a minor influence on the velocity fields. This result most likely stems from the low overall turbulent intensity of the flow. All the SGS models have in common, that they do not take into account the structure of the large-scale interface, which is expressed via $\overline{r}_\alpha$. Contrary, all model combinations concerning the extended gradient model (EGM) result in a lower vertical velocity inside the gas bubble. Furthermore they lead to steeper velocity gradients in the interface region at $\widehat{r}\approx {0.5}$. This results in smaller deviations from reference solution G4-O. For models combinations with scale similarity model SSM both effects are even more pronounced. As a consequence, the radial velocity profile obtained with M-SSM-SMAG is close to G4-O, while the ones resulting from KOB-SSM and KL-SSM are nearly identical to the reference. Note the small velocity offset in Figs. 7 and 8 for $\widehat{r}>1$, which reflects the different rising velocities.

Fig. 7

Radial profiles of vertical velocity component for eddy-viscosity and mixed models for the SGS convective term; the interface position for a volume-equivalent sphere is marked with a black solid vertical line

Fig. 8

Radial profiles of vertical velocity component for Kobayashi and Klein regularised models for the SGS convective term; the interface position for a volume-equivalent sphere is marked with a black solid vertical line

In the frame of modelling of the convective SGS term, structural models, which take the interfacial structure directly into account, turn out to improve overall simulation results. Regarding rising velocity and trajectory of the gas bubble, the model combination KL-EGM delivers the best predictions, while KOB-SSM and KL-SSM lead to improved radial velocity profiles. The results obtained with those models are qualitatively closer to high resolution simulation results compared to simulations without closure model.

4.2 Surface Tension Sub-grid Scale Term

In order to assess the SGS surface tension closure models, presented in Sect. 2.2.2, a test case with strong curvature of the interface is selected. For this purpose, the investigation of model influence is carried out in the two-dimensional benchmark case of Hysing et al. (2009), which is referred to as case 2 in the reference. Convergence and consistency of the underlying numerical framework in general as well as for this particular case is demonstrated by Meller et al. (2021). Deviating from the original configuration, the case is adapted in such a way, that both ratios of phase-specific density values and molecular viscosity values are reduced to 50. The reason for this choice is, that the method of Brackbill et al. (1992) for modelling surface tension is known to induce spurious currents. In the current case, those numerical artefacts are effectively dampened by increased gas density and viscosity. As the Weber-number based SGS models effectively increase the surface tension coefficient, in that way model influences can be investigated without interfering with the numerical deficiencies named before. This test case features a rectangular computational domain with slip conditions on both top and bottom as well as no-slip conditions on left and right boundaries. The gas bubble is initialised as a circle of diameter $D_{\mathrm{b}}$ with its centre $1D_{\mathrm{b}}$ above the lower boundary. The domain has a width and a length of $2D_{\mathrm{b}}$ and $4D_{\mathrm{b}}$, respectively. The characteristic dimensionless parameters are listed in Table 7. Two different spatial refinement levels with homogeneous, orthogonal numerical grids are investigated: with G1 the sphere-equivalent bubble diameter $D_{\mathrm{b}}$ is resolved with 10 grid cells, while with G2 $D_{\mathrm{b}}$ is equivalent to 160 grid cells. As reported by Meller et al. (2021) spatial resolution G2 is sufficiently fine to serve as a reference solution. Dimensionless time and velocity, $t_{g}$ and $U_{g}$, are defined as described in Sect. 4.1.1.

Table 7 Dimensionless parameters of two-dimensional test case considering rising gas bubble (G) in liquid (L)

The performance of SGS surface tension closure models is compared to simulation data obtained without model formulation on resolution G1 and G2, from now on referred to as G1-O and G2-O. No model influence is observed for model We-KE with $We _{\mathrm{crit}}={2.75}$ (We-KE-2.76). Therefore, the model is additionally assessed with critical Weber number $We _{\mathrm{crit}}={0.25}$ (We-KE-0.25) and the corresponding results are shown in Fig. 9. In reference solution G2-O thin ligaments form on both left and right sides of the bottom of the gas bubble due to the shear in the surrounding liquid flow, which results from the narrow computational domain. On grid G1 these structures are predicted to be comparably thick, because in particular small details cannot be resolved. Nevertheless, the necking, which is observed in result G2-O, is not recovered without SGS surface tension closure model on the coarse grid (G1-O). The model nn-SSM results in a qualitatively identical behaviour just as We-KE-2.76 does. In contrast, Weber-number models with a stricter Weber-number criterion, namely We-KE-0.25 and We-Re in this case, do allow for this particular process. Apart from that, results on G1 don’t show significant deviation among different SGS surface tension closure formulations.

Fig. 9

Interface position $\overline{r}_{\mathrm{G}}={0.3}$ at $\widehat{t}={4.2}$ for different SGS surface tension models with G1 as well as results with G1 and G2 without any SGS closure model; a whole gas bubble and b a zoom to ${1.5}<\widehat{x}<{1.8}$ and ${1.2}<\widehat{y}<{1.8}$

In order to understand, why the necking is recovered with Weber-number models in general, the distribution of model forces is shown in Fig. 10. At both $\widehat{t}={2.8}$ and $\widehat{t}={4.2}$ the scale similarity model nn-SSM delivers strong contributions in regions of higher gas volume fractions $\overline{r}_{\mathrm{G}}$ with the resulting force pointing in the direction of the liquid phase. This explains the negligible influence of this closure formulation in the present case. It is worth noting, that in this particular test case the interface defined as $0<\overline{r}_\alpha <1$ is predicted to be comparably thick due to the coarse spatial resolution G1. Hence, it is possible to observe contributions of the SGS model force in such a large area. The distribution of SGS surface tension over the interfacial region is considered in more detail by Hasslberger et al. (2020) in a-priori analysis. The Weber-number models We-KE-0.25 and We-Re reveal resulting model forces, which mainly appear at the lower side of the gas bubble acting in a direction, such that they counteract the curvature of the interface. In this way, the behaviour is reproduced, which was the theoretical starting point for formulating the Weber-number model (Ketterl et al. 2019). At $\widehat{t}={2.8}$ both Weber-number models reveal similar model forces, whereas at $\widehat{t}={4.2}$ We-KE-0.25 shows a smaller effect compared to We-Re, especially in the region and the ends of the ligaments. At the same time, We-Re shows a tendency to compress those slender gas structures. The effect of enhanced necking of the ligaments is common to both Weber-number models, We-KE-0.25 and We-Re. The model force contribution resulting from We-KE-2.76 is negligibly small and, therefore, is not shown here.

Fig. 10

Distribution of forces at $\widehat{t}={2.8}$ a–c and at $\widehat{t}={4.2}$ d–f resulting from G1 with different closure formulations for SGS surface tension $\varvec{\tau }_{rnn,\alpha \beta }$; the interface position $\overline{r}_{\mathrm{G}}={0.3}$ is marked with the black line

In this two-dimensional test case, the scale similarity model did not show any influence on simulation results, which applies to We-KE-2.76 as well. Contrary, We-KE-0.25 and We-Re formulations allow for the prediction of features of gas–liquid flows on a coarse numerical grid, which are typically reproduced on computational grids with high spatial resolutions. Those closure models help to predict a bubble shape, which is qualitatively similar to the results of a finer numerical grid without SGS modelling.

4.3 Combined Application of Convective and Surface Tension Sub-grid Scale Terms

After closure models for convective and surface tension SGS terms have been individually assessed previously, the combined effect of models for both terms is investigated. For this purpose, a single closure formulation for each individual SGS term is selected, namely We-Re for modelling the surface tension SGS term and KL-SSM for the convective SGS term. Both models perform very good, when applied and tested in the two- or three-dimensional test cases described earlier. In the following both closure formulations are simultaneously applied to both test cases.

4.3.1 Two-Dimensional Case

At first the test case from Sect. 4.2 is used to assess the closure models. For that purpose, interface position $\overline{r}_{\mathrm{G}}={0.3}$ at $\widehat{t}={4.2}$ is shown in Fig. 11 for the whole gas bubble and as zoom showing the right ligament of the gas bubble. Besides the results G2-O, G1-O and We-Re-G1, which are already presented in Fig. 9, KL-SSM is applied exclusively as well as in combination with We-Re. It turns out that KL-SSM shows a very minor overall influence on the shape of the gas structure. The result obtained purely with this closure formulation hardly differs from G1-O without model for any SGS term. Solely in the region of the ligament, the gas structure tends to be slightly more elongated downwards at the tip and to have a slight tendency towards enhanced necking at $\widehat{y}\approx {1.5}$, compared to G1-O. The same influence is observed, when investigating the application of KL-SSM together with We-Re. With both models combined the ligament tends to be slightly longer than with We-Re applied only. However, the overall influence of KL-SSM is small compared to the positive performance of We-Re in this two-dimensional test case.

Fig. 11

Interface position $\overline{r}_{\mathrm{G}}={0.3}$ at $\widehat{t}={4.2}$ for individual application of KL-SSM, We-Re and the combination of both models with G1 as well as results with G1 and G2 without any SGS closure model; a whole gas bubble and b a zoom to ${1.5}<\widehat{x}<{1.8}$ and ${1.2}<\widehat{y}<{1.8}$

4.3.2 Three-Dimensional Case

For the test case from Sect. 4.1 radial profiles of vertical velocity component are shown in Fig. 12. Besides the results G4-O, G2-O and Kl-SSM-G2, which are already presented in Fig. 8, the results from isolated application of We-Re as well as in combination with KL-SSM are presented here. It is evident, that the Weber-number model does not show any influence on the radial velocity profile. Hence, the results of We-Re-G2 and G2-O are nearly identical to each other. The same holds for results We-Re-G2 and KL-SSM+We-Re-G2. It turns out, that in this test case the gas bubble is not deformed strong enough, such that the interface shows a curvature, which is large enough for the Weber-number criterion of We-Re to activate the model.

Fig. 12

Radial profiles of vertical velocity component for individual application of KL-SSM, We-Re and the combination of both models with G2 as well as results with G2 and G4 without any SGS closure model; the interface position for a volume-equivalent sphere is marked with a black solid vertical line

In order to check the consistency of those results, the bubble rising paths as projection into the horizontal plane are investigated. Those are shown in Fig. 13 for refinement level G2 and for the configurations of SGS closure models named above. Indeed, the observation of minor influence of We-Re in this case is confirmed. Neither with or without the usage of the convective SGS model KL-SSM, the Weber-number models reveals an influence on the quality of the results. In case no convective SGS closure model is applied, the gas bubble describes a zig-zag path with moderate variance of spatial orientation of the lateral movement. With KL-SSM, the bubble shows longer distances of straight trajectories between random changes of orientation. Both statements hold, regardless of the application of We-Re.

Fig. 13

Trajectories of rising gas bubble projected into the horizontal plane; results are shown for different closure models for the convective and the surface tension SGS terms

Following the basic idea of the Weber-number model of Ketterl and Klein (2018) together with the specific formulation of the effective surface tension coefficient proposed in this work, this closure model must only be active in regions with extreme interface curvatures with respect to the numerical grid. As this situation is not met in the present test case, the model works as intended.

For both KL-SSM and We-Re it is concluded, that both individual SGS models improve the prediction in cases, where the respective unclosed SGS terms are of particular importance. At the same time both models do not change the results qualitatively in the opposite case.

5 Conclusion

In the present work the multifield two-fluid model equations are spatially low-pass filtered. In this way, a hybrid multiphase model combining Euler–Euler and algebraic volume-of-fluid methods in a single numerical framework is extended for scale resolving simulations of multiphase flows with large scale interface structures, particularly for gas–liquid flows. The prediction of such flows utilising coarse computational grids is successfully improved. From the five different unclosed sub-grid scale terms, arising from the filter formalism, convective and surface tension sub-grid scale contributions are selected to be modelled with several closure formulations each. The selected models are either of functional or of structural type and are adapted for the formulation in the hybrid multiphase model. Functional closure formulations for the convective sub-grid stress are additionally regularised in order to prevent anti-diffusivity and, therefore, numerically destabilising behaviour. For sub-grid surface tension, the functional Weber-number model is extended to rely on a Reynolds-number dependent critical Weber-number model, which avoids case specific definition of the latter quantity.

All models are individually tested in an a-posteriori fashion in cases of two- and three-dimensional rising gas bubbles in liquid. Among all investigated closure formulations for the convective sub-grid stress, structural models perform best. The reason is that they account for the structure of the flow field as well as the interfacial structure itself, which is expressed in terms of phase-specific volume fraction. Especially, the extended gradient model (Ketterl and Klein 2018) (EGM) together with regularisation according to Klein et al. (2020) (KL) delivers promising results overall. Regarding sub-grid surface tension, the functional approach of the Weber-number model successfully limits interface curvature according to the level of actual spatial resolution. Together with the criterion of a Reynolds-number dependent critical Weber-number, this approach delivers promising results in the two-dimensional test case without demanding for a case specific tuned model parameter. The combined application of one convective and one surface tension SGS model shows, that the prediction of the simulation is improved in situations, where the respective SGS term is of great relevance, while not changing the result qualitatively in other cases.

While the present work focuses on establishing closure models for unclosed SGS contributions arising from filtering across the phase interface in a-posteriori analysis, the investigated setups are comparatively simple with low and moderate degrees of turbulence. Future work shall focus on more complex bubble flow configurations with bubble swarms under more turbulent flow conditions considering interactions of interfaces with stronger velocity fluctuations. Furthermore, finding suitable combinations of convective and surface tension SGS closures shall be also part of future endeavours.