Investigation of the force closure for Eulerian-Eulerian simulations: a validation study of nine gas-liquid flow cases

TL;DR

This study evaluates the impact of different force models in Eulerian-Eulerian gas-liquid flow simulations, identifying key forces necessary for accurate predictions across various industrial applications and validating results against experimental data.

Contribution

It provides a comprehensive validation of force closure models in Eulerian-Eulerian simulations for diverse industrial bubbly flows, establishing guidelines for numerical settings.

Findings

Drag and turbulent dispersion forces are essential for accurate predictions.

Lift and wall lubrication forces are important near walls, especially in pipe flows.

Neglecting lateral forces is acceptable in some flow configurations.

Abstract

Gas-liquid flows can be simulated by the Eulerian-Eulerian (E-E) method. Whether to include a specific momentum interfacial exchange force model remains as a question with no answer. In this work, our aim is to seek a general numerical settings for the E-E method, which can provide competent results for industrial bubbly flows with different geometries under different operations. They were selected from different industries including chemical, nuclear, bio-processing and metallurgical engineering. Simulations were launched by the OpenFOAM solver reactingTwoPhaseEulerFoam, in which the E-E method was implemented with sophisticated numerical techniques to ensure stabilities. Predictions were compared against experimental data. It was found that the drag force and turbulent dispersion force play the most important role on the predictions and should be included for all simulations. The first one accounts for the two-way coupling while the second one accounts for the turbulence effect and ensures the E-E equations to be well-posed. The lift force and wall lubrication force should be included to address the phase fraction accumulation in the vicinity of the wall, especially for pipe flows with large aspect ratio. In other cases the lateral forces can be safely neglected. All the test case are open-sourced and are available as supplementary data for anyone to download as baseline test cases.

Tables

Table 1. Table 1: Numerical configurations used in the test cases: ψ 𝜓 \psi denotes a generic variable. ( … ) f subscript … 𝑓 (...)_{f} is the face interpolation operator. ∇ ⟂ superscript ∇ perpendicular-to \nabla^{\perp} is the surface-normal gradient. The number “1” indicates the compliance of the scheme with the definition of TVD scheme. A value of 1 indicates full TVD compliance.

Term	Configuration
$\partial /\partial t$	Euler implicit
$\nabla \psi$	Gauss linear
$\nabla p$	Gauss linear
$\nabla \cdot \left(\psi \mathrm{𝐔𝐔}\right)$	Gauss limitedLinearV 1;
$\nabla \cdot \left(𝐔\psi \right)$	Gauss limitedLinear 1
$\nabla \cdot \tau$	Gauss linear
${\nabla }^2\psi$	Gauss linear uncorrected
${\nabla }^{⟂}\psi$	uncorrected
${\psi }_f$	linear

Table 2. Table 2: Solvers and related settings used in the test cases.

	Solver	Preconditioner	Rel. tol.	Final tol.
$p$	PCG	DIC	0.01	1e-7
$k$	PBiCGStab	DILU	-	1e-7
$\epsilon$	PBiCGStab	DILU	-	1e-7

Table 3. Table 3: Main boundary conditions adopted in the simulations.

Test case	Operating conditions
A.1	Superficial gas velocity: 0.0024 m/s Inlet liquid velocity: 0 m/s Bubble diameter: 0.00505 m
A.2	Superficial gas velocity: 0.04 m/s Inlet liquid velocity: 0 m/s Bubble diameter: 0.004 m
A.3	Inlet gas velocity: 1.87 m/s Inlet liquid velocity: 1.57 m/s Bubble diameter: 0.002 m
B.1	Superficial gas velocity: 0.0115 m/s Superficial liquid velocity: 1.0167 m/s Bubble diameter: 0.0048 m
B.2	Superficial gas velocity: 0.0151 m/s Superficial liquid velocity: 1.017 m/s Bubble diameter: 0.0046 m
B.3	Superficial gas velocity: 0.005 m/s Superficial liquid velocity: 0.43 m/s Bubble diameter: 0.006 m
B.4	Inlet gas velocity: 0.0087 m/s Inlet liquid velocity: 0 m/s Bubble diameter: 0.0042 m
C.1	Inlet gas velocity: 4 cm3/s Inlet liquid velocity: 5 l/s Bubble diameter: 0.005 m
C.2	Inlet gas velocity: 50 ml/s Inlet liquid velocity: 0 m/s Bubble diameter: 0.006 m
Gas velocity at walls: slip. Liquid velocity at walls: no-slip. $k$ and $\epsilon$ at walls: wall function. Outlet: zero-gradient.

Table 4. Table 4: Available experimental data and features of test case A.1 to A.3.

Test case	Exp. data	Features	Pseudo-steady-state
A.1	Plume oscillating period Gas holdup	Periodic flow field	No
A.2	Phase frac. distri.	High phase fraction	No
A.3	Upward liquid vel. Turb. kinetic energy	Stagnant vortex	Yes

Table 5. Table 5: Comparison of the gas holdup and POP predicted by different force closure combination with experimental data for test case A.1.

	Drag force	Drag & dispersion forces	Drag & dispersion & lift and wall forces	Exp.
Gas holdup	0.00613	0.0065	0.0075	0.0069
Gas holdup rel. error	-11%	-6%	+8%	-
POP	9.37	8.87	8.98	11.38
POP rel. error	-17%	-22%	-21%	-

Table 6. Table 6: Available experimental data and features of test case B.1 to B.4.

Test case	Exp. data	Features	Pseudo-steady-state
B.1	Phase frac. distri. Upward gas vel.	Wall peak	Yes
B.2	Phase frac. distri.	Double peak	Yes
B.3	Phase frac. distri.	Central peak	Yes
B.4	Global gas holdup	Non-regular components	No

Table 7. Table 7: Predicted gas holdup compared with the experimental data for case B.4.

	Drag	Drag & Dis.	Drag & Dis. & Lift & Wall	Exp.
Gas holdup	0.0226	0.02318	0.02358	0.0287

Table 8. Table 8: Available experimental data and features of test case C.1 to C.2.

Test case	Exp. data	Features	Pseudo-steady-state
C.1	Axial liquid vel. Radial liquid vel.	Non-regular flow	Yes
C.2	Axial liquid vel. Turbulent kinetic energy	Central bubble plume	Yes

Equations

$$\frac{\partial \left({\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}\right)}{\partial t}+\nabla \cdot \left({\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}}\right)=0,$$

$$\frac{\partial \left({\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}\right)}{\partial t}+\nabla \cdot \left({\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}\right)=0,$$

$$\begin{gathered} \frac{\partial \left({\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}}\right)}{\partial t}+\nabla \cdot \left({\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}}\right)- \\ \nabla \cdot \left({\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}{\tau }_{\mathrm{a}}\right)=-{\alpha }_{\mathrm{a}}\nabla p_{\mathrm{a}}+{\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}\mathbf{g}-\mathbf{M}, \end{gathered}$$

$$\begin{gathered} \frac{\partial \left({\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}\right)}{\partial t}+\nabla \cdot \left({\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}\right)- \\ \nabla \cdot \left({\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\tau }_{\mathrm{b}}\right)=-{\alpha }_{\mathrm{b}}\nabla p_{\mathrm{b}}+{\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}\mathbf{g}+\mathbf{M}, \end{gathered}$$

$${\mathbf{M}}_{\mathrm{drag}}=\frac{3}{4}{\alpha }_{\mathrm{a}}{\rho }_{\mathrm{b}}{C}_{\mathrm{D}}\frac{1}{d_{\mathrm{a}}}\left|{\mathbf{U}}_{\mathrm{b}}-{\mathbf{U}}_{\mathrm{a}}\right|\left({\mathbf{U}}_{\mathrm{b}}-{\mathbf{U}}_{\mathrm{a}}\right),$$

$${\mathbf{M}}_{\mathrm{lift}}={\alpha }_{\mathrm{b}}{C}_{\mathrm{L}}{\rho }_{\mathrm{a}}{\mathbf{U}}_{\mathrm{r}}\times \left(\nabla \times {\mathbf{U}}_{\mathrm{b}}\right),$$

$${\mathbf{M}}_{\mathrm{wall}}=C_{\mathrm{wall}}{\rho }_{\mathrm{b}}{\alpha }_{\mathrm{a}}|{\mathbf{U}}_{\mathrm{r}}{|}^2\cdot \mathbf{n},$$

$${\mathbf{M}}_{\mathrm{turb}}=C_{\mathrm{T}}{\rho }_{\mathrm{b}}{k}_{\mathrm{b}}\nabla {\alpha }_{\mathrm{a}},$$

$$\begin{gathered} \frac{\partial {\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}}}{\partial t}+\nabla \cdot \left({\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}}\right)- \\ \nabla \cdot \left({\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}{\tau }_{\mathrm{a}}\right)=-\mathbf{Kd}{\mathbf{U}}_{\mathrm{a}}-{\mathbf{M}}_{\mathrm{lift}}-{\mathbf{M}}_{\mathrm{turb}}-{\mathbf{M}}_{\mathrm{wall}}, \end{gathered}$$

$$\begin{gathered} \frac{\partial {\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}}{\partial t}+\nabla \cdot \left({\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}\right)- \\ \nabla \cdot \left({\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\tau }_{\mathrm{b}}\right)=-\mathbf{Kd}{\mathbf{U}}_{\mathrm{b}}+{\mathbf{M}}_{\mathrm{lift}}+{\mathbf{M}}_{\mathrm{turb}}+{\mathbf{M}}_{\mathrm{wall}}, \end{gathered}$$

$$\mathbf{Kd}=\frac{3}{4}{\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}{\mathrm{C}}_{\mathrm{D}}\frac{1}{d_{\mathrm{b}}}|{\mathbf{U}}_{\mathrm{a}}-{\mathbf{U}}_{\mathrm{b}}|.$$

$$A_{\mathrm{a},\mathrm{P}}{\mathbf{U}}_{\mathrm{a},\mathrm{P}}+\sum A_{\mathrm{a},\mathrm{N}}{\mathbf{U}}_{\mathrm{a},\mathrm{N}}=S_{\mathrm{a},\mathrm{P}},$$

$$A_{\mathrm{b},\mathrm{P}}{\mathbf{U}}_{\mathrm{b},\mathrm{P}}+\sum A_{\mathrm{b},\mathrm{N}}{\mathbf{U}}_{\mathrm{b},\mathrm{N}}=S_{\mathrm{b},\mathrm{P}},$$

$${\mathbf{HbyA}}_{\mathrm{a},\mathrm{P}}=\frac{1}{A_{\mathrm{a},\mathrm{P}}}(-\sum A_{\mathrm{a},\mathrm{N}}{\mathbf{U}}_{\mathrm{a},\mathrm{N}}+S_{\mathrm{a},\mathrm{P}}),$$

$${\mathbf{HbyA}}_{\mathrm{b},\mathrm{P}}=\frac{1}{A_{\mathrm{b},\mathrm{P}}}(-\sum A_{\mathrm{b},\mathrm{N}}{\mathbf{U}}_{\mathrm{b},\mathrm{N}}+S_{\mathrm{b},\mathrm{P}}),$$

$$\nabla p_{\mathrm{rgh}}=\nabla p-\rho \mathbf{g}-\mathbf{g}\cdot \mathbf{h}\nabla \rho ,$$

$$\begin{gathered} -{\alpha }_{\mathrm{a}}\nabla p+{\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}\mathbf{g}=-{\alpha }_{\mathrm{a}}\nabla p_{\mathrm{rgh}}-{\alpha }_{\mathrm{a}}\rho \mathbf{g}-{\alpha }_{\mathrm{a}}\mathbf{g}\cdot \mathbf{h}\nabla \rho +{\alpha }_{\mathrm{a}}{\rho }_{\mathrm{a}}\mathbf{g} \\ ={\alpha }_{\mathrm{a}}\nabla p_{\mathrm{rgh}}-{\alpha }_{\mathrm{b}}{\alpha }_{\mathrm{a}}({\rho }_{\mathrm{b}}-{\rho }_{\mathrm{a}})\mathbf{g}-{\alpha }_{\mathrm{a}}\mathbf{g}\cdot \mathbf{h}\nabla \rho , \end{gathered}$$

$$\begin{gathered} -{\alpha }_{\mathrm{b}}\nabla p+{\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}\mathbf{g}=-{\alpha }_{\mathrm{b}}\nabla p_{\mathrm{rgh}}-{\alpha }_{\mathrm{b}}\rho \mathbf{g}-{\alpha }_{\mathrm{b}}\mathbf{g}\cdot \mathbf{h}\nabla \rho +{\alpha }_{\mathrm{b}}{\rho }_{\mathrm{b}}\mathbf{g} \\ ={\alpha }_{\mathrm{b}}\nabla p_{\mathrm{rgh}}-{\alpha }_{\mathrm{a}}{\alpha }_{\mathrm{b}}({\rho }_{\mathrm{a}}-{\rho }_{\mathrm{b}})\mathbf{g}-{\alpha }_{\mathrm{b}}\mathbf{g}\cdot \mathbf{h}\nabla \rho , \end{gathered}$$

$$\nabla \cdot ({\alpha }_{\mathrm{b}}{\mathbf{U}}_{\mathrm{b}}+{\alpha }_{\mathrm{a}}{\mathbf{U}}_{\mathrm{a}})=0.$$

$$\sum ({\alpha }_{\mathrm{a},\mathrm{f}}{\mathbf{U}}_{\mathrm{a},\mathrm{f}}+{\alpha }_{\mathrm{b},\mathrm{f}}{\mathbf{U}}_{\mathrm{b},\mathrm{f}})\cdot {\mathbf{S}}_{\mathrm{f}}=0.$$

$${\mathbf{U}}_{\mathrm{a},\mathrm{P}}=\mathbf{Hby}{\mathbf{A}}_{\mathrm{a},\mathrm{P}}+\frac{{\alpha }_{\mathrm{a},\mathrm{P}}}{{\mathbf{A}}_{\mathrm{a},\mathrm{P}}}(\nabla p_{\mathrm{rgh},\mathrm{P}}-{\alpha }_{\mathrm{b},\mathrm{P}}({\rho }_{\mathrm{b}}-{\rho }_{\mathrm{a}})\mathbf{g}-\mathbf{g}\cdot {\mathbf{h}}_{\mathrm{P}}\nabla {\rho }_{\mathrm{P}})+\frac{\mathrm{K}{\mathrm{d}}_{\mathrm{P}}}{{\mathbf{A}}_{\mathrm{a},\mathrm{P}}}{\mathbf{U}}_{\mathrm{b},\mathrm{P}},$$

$${\mathbf{U}}_{\mathrm{b},\mathrm{P}}=\mathbf{Hby}{\mathbf{A}}_{\mathrm{b},\mathrm{P}}+\frac{{\alpha }_{\mathrm{b},\mathrm{P}}}{{\mathbf{A}}_{\mathrm{b},\mathrm{P}}}(\nabla p_{\mathrm{rgh},\mathrm{P}}-{\alpha }_{\mathrm{a},\mathrm{P}}({\rho }_{\mathrm{a}}-{\rho }_{\mathrm{b}})\mathbf{g}-\mathbf{g}\cdot {\mathbf{h}}_{\mathrm{P}}\nabla {\rho }_{\mathrm{P}})+\frac{\mathrm{K}{\mathrm{d}}_{\mathrm{P}}}{{\mathbf{A}}_{\mathrm{b},\mathrm{P}}}{\mathbf{U}}_{\mathrm{a},\mathrm{P}}.$$

$$\alpha_{\rm b,f}\phi_{\rm b,f}+\alpha_{\rm a,f}\phi_{\rm a,f}=\nabla\cdot\Bigg{(}\Big{(}\alpha_{\rm b,P}\frac{\alpha_{\rm b,P}}{\mathbf{A}_{\rm b,P}}+\alpha_{\rm a,P}\frac{\alpha_{\rm a,P}}{\mathbf{A}_{\rm a,P}}\Big{)}\nabla p_{\rm rgh,P}\Bigg{)},$$

$${\varphi }_{\mathrm{a},\mathrm{f}}=(\mathbf{Hby}{\mathbf{A}}_{\mathrm{a},\mathrm{f}}+\frac{{\alpha }_{\mathrm{a},\mathrm{f}}}{{\mathbf{A}}_{\mathrm{a},\mathrm{f}}}(-{\alpha }_{\mathrm{b},\mathrm{f}}({\rho }_{\mathrm{b}}-{\rho }_{\mathrm{a}})\mathbf{g}-\mathbf{g}\cdot {\mathbf{h}}_{\mathrm{f}}\nabla {\rho }_{\mathrm{f}})+\frac{\mathrm{K}{\mathrm{d}}_{\mathrm{f}}}{{\mathbf{A}}_{\mathrm{a},\mathrm{f}}}{\mathbf{U}}_{\mathrm{b},\mathrm{f}})\cdot {\mathbf{S}}_{\mathrm{f}},$$

$${\varphi }_{\mathrm{b},\mathrm{f}}=(\mathbf{Hby}{\mathbf{A}}_{\mathrm{b},\mathrm{f}}+\frac{{\alpha }_{\mathrm{b},\mathrm{f}}}{{\mathbf{A}}_{\mathrm{b},\mathrm{f}}}(-{\alpha }_{\mathrm{a},\mathrm{f}}({\rho }_{\mathrm{a}}-{\rho }_{\mathrm{b}})\mathbf{g}-\mathbf{g}\cdot {\mathbf{h}}_{\mathrm{f}}\nabla {\rho }_{\mathrm{f}})+\frac{\mathrm{K}{\mathrm{d}}_{\mathrm{f}}}{{\mathbf{A}}_{\mathrm{b},\mathrm{f}}}{\mathbf{U}}_{\mathrm{a},\mathrm{f}})\cdot {\mathbf{S}}_{\mathrm{f}}.$$

$${\mathbf{U}}_{\mathrm{a},\mathrm{P}}^{\mathrm{s}}=\mathbf{Hby}{\mathbf{A}}_{\mathrm{a},\mathrm{P}}+\frac{{\alpha }_{\mathrm{a},\mathrm{P}}}{{\mathbf{A}}_{\mathrm{a},\mathrm{P}}}(\nabla p_{\mathrm{rgh},\mathrm{P}}-{\alpha }_{\mathrm{b},\mathrm{P}}({\rho }_{\mathrm{b}}-{\rho }_{\mathrm{a}})\mathbf{g}-\mathbf{g}\cdot {\mathbf{h}}_{\mathrm{P}}\nabla {\rho }_{\mathrm{P}}),$$

$${\mathbf{U}}_{\mathrm{b},\mathrm{P}}^{\mathrm{s}}=\mathbf{Hby}{\mathbf{A}}_{\mathrm{b},\mathrm{P}}+\frac{{\alpha }_{\mathrm{b},\mathrm{P}}}{{\mathbf{A}}_{\mathrm{b},\mathrm{P}}}(\nabla p_{\mathrm{rgh},\mathrm{P}}-{\alpha }_{\mathrm{a},\mathrm{P}}({\rho }_{\mathrm{a}}-{\rho }_{\mathrm{b}})\mathbf{g}-\mathbf{g}\cdot {\mathbf{h}}_{\mathrm{P}}\nabla {\rho }_{\mathrm{P}}).$$

$${\mathbf{U}}_{\mathrm{a},\mathrm{P}}={\mathbf{U}}_{\mathrm{a},\mathrm{P}}^{\mathrm{s}}+{\mathrm{D}}_{\mathrm{a},\mathrm{P}}{\mathbf{U}}_{\mathrm{b},\mathrm{P}},$$

$${\mathbf{U}}_{\mathrm{b},\mathrm{P}}={\mathbf{U}}_{\mathrm{b},\mathrm{P}}^{\mathrm{s}}+{\mathrm{D}}_{\mathrm{b},\mathrm{P}}{\mathbf{U}}_{\mathrm{a},\mathrm{P}},$$

$${\mathbf{U}}_{\mathrm{P}}={\alpha }_{\mathrm{b},\mathrm{P}}({\mathbf{U}}_{\mathrm{b},\mathrm{P}}^s+{\mathrm{D}}_{\mathrm{b},\mathrm{P}}{\mathbf{U}}_{\mathrm{b},\mathrm{P}})+{\alpha }_{\mathrm{a},\mathrm{P}}({\mathbf{U}}_{\mathrm{a},\mathrm{P}}^{\mathrm{s}}+{\mathrm{D}}_{\mathrm{a},\mathrm{P}}{\mathbf{U}}_{\mathrm{b},\mathrm{P}}),$$

Full text

††journal: Industrial & Engineering Chemistry Research

1 Introduction

Gas-liquid flows are encountered in a variety of applications and can be simulated by different models, such as the Eulerian-Eulerian (E-E) method, the Eulerian-Lagrangian (E-L) method and the direct method which is also called multi-phase direct numerical simulation (DNS). In the E-E method [1, 2, 3, 4, 5, 6, 7, 8, 9], each fluid phase is considered as a continuum in the computational domain under consideration which can interpenetrate with the other fluid phase. Several averaging methods (such as volume or ensemble averaging) are used to formulate basic governing equations. In the E-L approach [2, 10, 11, 12], the continuous phase is treated in an Eulerian framework whereas the motion of individual bubbles is simulated by solving the force balance equation. The trajectories of these bubbles are computed in the control volume and averaged at the computational level. In the multi-phase DNS [13, 14, 15, 16], one employs the Navier-Stokes equations directly, without further manipulation and the topology of the interface between the two-phases is determined as part of the solution. No additional modelling assumptions are introduced. DNS requires very high resolution in order to resolve a broad range of temporal and spatial scales. These scales are associated with the topology of the interface, e.g. the size of the dispersed particles, or with the fluid motion, e.g. the eddies encountered in the turbulent motion. Resolving these scales is computationally expensive both in terms of computer memory size execution time. Therefore, DNS is restricted to only low Reynolds numbers and a few bubbles/particles due to its high computational cost.

Thanks to the computational economics, many works used the E-E method to simulate gas-liquid flows. However, as a macro-scopic method, the E-E method 1) requires constitutive models (e.g., the solid-phase stress model); 2) requires models to address the momentum and energy exchange terms (e.g., the drag model); 3) loses the characteristic feature for those scales which is much smaller than the mesh resolution. To obtain reasonable results, suitable sub-models and parameters should be adjusted. Moreover, the observed flow field is also different for the gas-liquid flows in different industries. For example, bubble columns are quite common in chemical engineering industry. The liquid flow field in the bubble column tends to form a circulation model due to the movement of the bubbles, which is commonly referred as “Gulf-stream” or “cooling tower” flow regime. In the nuclear engineering, the phase fraction of the bubbles in vertical upward pipe shows a wall peak due to the lateral movement of the small bubbles [17, 18, 19]. Without loss of the generality, it was proven that the E-E method is able to capture the typical flow field information. However, amount of work needs to be studied for the sub-models, especially for the momentum interfacial exchange models.

The drag force and the buoyancy force, as the most important momentum interfacial exchange forces, determine the bubble terminal velocity and address the coupling between the disperse phase with the continuous phase. The lift force and wall lubrication force address the lateral movements of the bubbles. A correct description of the lift coefficient and wall force coefficient is crucial in order to model this transversal force correctly. The turbulence dispersion force acts as a driving force for bubbles to move from areas with higher phase fractions to areas of lower phase fraction. It arises due to the pressure variations in the continuous phase that are not resolved at meso-scale. It is also shown by Panicker et al. [20] and Vaidheeswaran and Lopez de Bertodano [21] that the addition of a dispersion term ensures the hyperbolicity of the PDE, and prevents the non-physical instabilities in the predicted multiphase flows upon grid refinement.

Although there is universal consensus in the literature on the need to incorporate drag into any model of a bubble column, many works admit that there is still no agreement in the community on the other forces to be used at best [22]. The purpose of the present contribution is to simulate large amount of gas-liquid flows with different force closure combinations to investigate their importance. Our aim is to seek a general numerical settings for the E-E method, which can provide competent results for industrial bubbly flows.

2 Description of models

Within the Eulerian framework, two sets of Navier-Stokes equations are ensemble-averaged, and the effects of turbulence and inter-phase phenomena are taken into account using closure models. The conservation of mass for phase $a$ and phase $b$ can be expressed by [23] :

[TABLE]

[TABLE]

where $\alpha_{\rm a}$ and $\alpha_{\rm b}$ are the phase fraction of phase $a$ and phase $b$, $\rho_{\rm a}$ and $\rho_{\rm b}$ are the density for phase $a$ and phase $b$, and $\mathbf{U}_{\rm a}$ and $\mathbf{U}_{\rm b}$ are the average velocity for phase $a$ and phase $b$, respectively. The average velocities $\mathbf{U}_{\rm a}$ and $\mathbf{U}_{\rm b}$ can be calculated by solving the corresponding phase momentum equations [23]:

[TABLE]

[TABLE]

where $p_{\rm a}$ and $p_{\rm b}$ are the pressure for each phase, $\mathbf{\tau}_{\rm a}$ and $\mathbf{\tau}_{\rm b}$ are the effective Reynolds stress tensors, $\mathbf{g}$ is the gravitational acceleration vector, and $\mathbf{M}$ is the interfacial force exchange term. It is common to break it down into the axial drag force $\mathbf{M}_{\mathrm{drag}}$; the lateral lift force $\mathbf{M}_{\mathrm{lift}}$, which acts perpendicular to the direction of the relative motion of the two phases [24]; the wall lubrication force $\mathbf{M}_{\mathrm{wall}}$ which acts to drive the bubble away from the wall [25]; and the turbulent dispersion force $\mathbf{M}_{\mathrm{turb}}$, which is the result of the turbulent fluctuations of the liquid velocity [26]. Specifically, the drag force can be calculated as follows:

[TABLE]

where $C_{\mathrm{D}}$ is the drag force coefficient and $d_{\rm a}$ is the diameter of phase $a$. The lift force can be calculated as follows [27]:

[TABLE]

where $C_{\mathrm{L}}$ is the lift force coefficient, $\mathbf{U}_{\rm r}$ is the relative velocity which equals $\mathbf{U}_{\rm b}-\mathbf{U}_{\rm a}$. The wall lubrication force can be calculated as follows [25]:

[TABLE]

where $C_{\rm wall}$ is the wall lubrication force coefficient. The turbulence dispersion force can be calculated as follows [26]:

[TABLE]

where $C_{\rm T}$ is the turbulence dispersion force coefficient. A comprehensive discussion of the E-E method can be found in other latest work [28, 29]. Readers are referred to the our previous work [30] and the Appendix for the finite volume discretization of the E-E method.

In the E-E method the interaction between the bubbles and the liquid phase is modelled through the momentum interfacial exchange forces in the momentum equation of the liquid and the gas phase. There is still no agreement in the community on the closures to be used at best. Physically, the effect of these forces are reported in Fig. 1. In a quiescent environmental, the buoyancy force pushes the bubbles upwardly. On the other hand, the drag force tends to exert on the opposite direction of the bubble movement. The direction of the drag force is equivalent to the direction of relative velocity. In the simplest cases, these two forces can reach to a balance and the terminal bubble velocity can be determined.

The lift force arises due to the presence of shear in the liquid phase and acts perpendicular to the bubble rise direction. Its direction depends on the value of sign of $C_{\rm L}$. Different models can be used, such as the constant lift coefficient model and the model developed by Tomiyama [24]. Experiments show that the movement of large bubbles and small bubbles is opposite, especially in pesudo-steady-state pipe flows. It can only be predicted by employing a non-constant lift coefficient. Thus, the lift force coefficient developed by Tomiyama is more reasonable. Moreover, the largest shear exists not only in the vicinity of no-slip walls, it may be also large when one phase was injected into another phase with high velocity. The wall lubrication force exists only in the vicinity of solid objects. It pushes the bubbles away from the wall avoiding bubble accumulation which is not observed from the experiments. It should be noted that this force exists only near the walls. Thus, the distance between the bubble and the nearest wall should be calculated in the wall force model. The turbulent dispersion force acts on the gradient of the phase fraction. It can be seen as a diffusion and it accounts for the random bubble movement in turbulent flows. It was found in many works that the results predicted by the E-E method is highly mesh depended and it was due to the mathmatical characteristic of the E-E equations [21, 20]. Including the turbulent dispersion force avoids the ill-posed problems and improves the stability of the E-E model.

A large amount of force models were developed in the literature and it is not possible to compare all of them for all the test cases investigated in this work. Therefore, unless stated on intention, we employ the classic drag model by Ishii and Zuber [31], the lift force model by Tomiyama et al. [24], the wall force by Antal et al. [25] and the turbulent dispersion force by Loptez De Bertodano [26] in our simulations. Since we focus on the validation between numerical results and experiments, readers are suggested to Appendix for the discretization of the governing equations.

3 Test cases

As previously mentioned, the lack of an ideal experiment, with a good mix of local and global measurements performed under a wide spectrum of operating conditions, necessary for a complete model validation, brought us to consider many different test cases with different geometries. Only by using a large amount of test cases can we find a general numerical settings. These test cases were selected from different works. The schematic representation of these test cases is reported in Fig. 2. All of them were investigated by experiments which implies they are suitable benchmark test cases for numerical investigations. Moreover, they were observed with different features as listed as follows:

Test case A.1: investigation of bubble plumes in a rectangular bubble column of Díaz et al. [32]. The experimental equipment consists of a 0.2 m wide, 1.8 m high and 0.04 m deep partially aerated bubble column, filled with tap water up to 0.45 m from the bottom at room temperature and atmospheric pressure, while the air is fed through a sparger composed of eight centered holes of 1 mm of diameter and 6 mm pitch. This test case was proven to be a very interesting test case because the liquid vortices generated by the bubble plumes are a favorable factor for mixing and, therefore, for speeding up all transport processes [33]. Additionally, the existence of flow structures showing unsteady liquid recirculation is a typical phenomenon in industrial-scale bubble columns.

• 2.

Test case A.2: investigation of gas-liquid flows in an industrial bubble column of McClure et al. [34]. It is a partially aerated cylindrical bubble column of 0.19 m diameter with a multi-point sparger filled with water up to 1 m from the bottom. The aspect ratio ($L/D$) is about 5. In the bio-processing industry, the bubble column height to diameter typically ranges from 2 to 5 and is usually operated in the heterogeneous flow regime to speed up mass and heat transfer.

• 3.

Test case A.3: investigation of sudden enlargement pipe flow of Bel F’dhila [35]. The test case studied here is that of a bubbly air/water upward flow through a pipe with a sudden enlargement. The diameters of the two pipe sections are 50 mm and 100 mm, respectively. The inlet phase fraction is characterised by a wall peak in experiments. The feature of this test case is that there is a large separation zone at the bottom of the pipe. This case has been employed extensively to verify the E-E method and implementation [36, 37, 38, 39].

• 4.

Test case B.1: investigation of bubbly flow in a cyliderial pipe of Lucas et al. [40]. Mixture of the gas and liquid is injected from the bottom pipe of 0.0256 m diameter and 3.53 m height. For the vertical upward flow smaller bubbles tend to move towards the wall. A wall peak of the gas phase fraction occurs at the position of high $L/D$. This was also observed for single bubbles by Tomiyama et al. [41]. In the case of vertical co-current pipe flow the radial flow field is symmetrically stable over a long distance. Therefore, this type of flow is well suited for the investigation of the non-drag forces.

• 5.

Test case B.2: investigation of bubbly flow in a circular pipe of Banowski et al. [42]. The experiment comprises a vertical pipe with 54.8 mm inner diameter and 6 m length. Gas and liquid mixture is injected from the bottom. It is similar with the test case B.2. However, besides the typical wall peak formed due to the smaller bubbles movement, a double peak of the phase fraction can also be observed due to the existence of large bubbles because the large bubbles tend to move to the center.

• 6.

Test case B.3: investigation of bubbly flow in a rectangle pipe of Žun [43]. This test case is similar with test case B.3 with slightly difference of the geometry. Mixture of the gas and liquid is injected from the bottom of a rectangle channel of 0.0254 m length and 2 m height. Only the central peak of the gas phase fraction was observed.

• 7.

Test case B.4: investigation of bubbly flow of Besagni et al. [44]. The gas-liquid flow in an annular gap bubble column with two non-regular internal pipes was investigated. It consists of a non-pressized vertical column with an inner diameter 0.24 m and a height 5.3 m. In simulations the height of the domain is limited to 5 m. Two internal pipes are placed inside the column: one centrally positioned (with an external diameter of 0.06 m) and one asymmetrically positioned (with an external diameter of 0.075 m). The sparger is modeled as a uniform cylindrical surface with a height of 0.01 m placed on the lateral inner pipe at the vertical position of 0.3 m from the bottom of the domain. The aspect ratio of the geometry is small. Due to the existence of the non-irregular components, the gas phase fraction distribution developed to be quite flatten and no wall peak was observed in the experiments, even it can be seen as a pipe flow.

• 8.

Test case C.1: investigation of gas-liquid flows in a continuous casting molds of Iguchi and Kasai [45]. The geometry employed in this test case is quite different with previous ones. The gas and liquid is injected into a rectangular vessel of 0.3 m length and 0.15 m width. In the experiments, it was observed that larger bubbles are lifted towards the liquid surface due to the buoyancy force acting on them, while smaller bubbles are carried deep. Such phenomenon is also known as phase segregation or poly-dispercity in other works, which was proven as a tough work for the E-E method.

• 9.

Test case C.2: investigation of gas-liquid flows in a continuous casting molds of Sheng and Irons [46]. The gas phase is injected from the bottom of the vessel of 0.76 m height and 0.5 m diameter. In this test case, the measured turbulence fields, gas phase fraction distribution, gas/liquid velocities in the plume zone were used for validation of various turbulence models. It can be seen as a suitable test case to validate the multi-phase turbulence model against experimental data.

The reader may find these test cases are sorted by different groups. The bubble columns investigated in test case A.1 - A.3 are mainly encountered in chemical and bio-processing engineering. The aspect ratio of the geometry is small. They are usually used to speed up mixing and heat/mass transfer and are typical operated in medium or high superficial velocity. A strong coupling between the disperse phase and continuous phase is observed. The local phase fraction may be high. Meanwhile, the liquid flow in these bubble columns may be highly transient and full of chaos with large-scale vortices. The pipe flow investigated in test case B.1 - B.4 are mainly encountered in nuclear engineering process. The aspect ratio of these pipes are relatively large. Different flow types exist depending on the gas flow rate. Dilute bubbly flows are found when operated at low gas superficial velocity. In bubbly flows, the liquid flow is rather stable and a steady-state can be achieved. The shape of the phase fraction distribution develops gradually to a stable distribution along the pipe axial direction, since the aspect ratio is quite large which implies the gas phase has enough time to develop. The bubbly flows investigated in test case C.1 - C.2 are mainly encountered in metallurgical engineering. The scale of the gas-liquid reactor is the same with that employed in the chemical engineering. The aspect ratio is also small. However, non-irregular design is quite common and the flow field is complex due to the existence of the complex geometry. Meanwhile, the small bubbles in the liquid phase, due to the phase segregation movement, are usually seen as impurity which need to be removed. At last, it should be stressed that although these test cases coming from different industries are full of different features, the core problem lies in the investigation of the gas-liquid flows by the numerical method.

The computational grids employed in these test cases are shown in Fig. 3. 3D hexahedral cells were generated for test case A.1, A.2 and C.1. 2.5D axi-symmetric wedge cells were employed for test case A.3, B.1, B.2 and C.2. 2D hexahedral cells were employed for test case B.3. Non-regular gambit-paving meshes was employed for test case B.4 due to the existence of the non-regular geometries. Our grid independence investigation, though not shown herein for brevity, suggested that the predicted results are not sensitive to the grid resolution. The base setup for these test cases are listed in Table 1 and 2. In our preliminary investigations such numerical setup compromises between stability and accuracy. Therefore, the remainder of the simulations in this study will utilize this base setup. The main boundary conditions are listed in Table 3. Since large amount test cases were employed in this work, it is not possible to document all the settings for brevity. Readers are referred to the source code of the test cases for further details.

All the test cases were simulated by three basic settings. The first one only consists of the drag force. The second one consists of the drag force and turbulence dispersion force. The third one consists of the drag force, turbulence dispersion force, lift force and wall force. Our aim is to verify which force closure combination is able to provide the most universal settings.

4 Results and discussions

4.1 Test case A.1 - A.3

In this section, the numerical predicted results for test case A.1 to A.3 were investigated. The available experimental data and features are listed in Table 4. These test cases are common in chemical and bio-processing engineering. In test case A.1, the gas phase was injected into the column at a relatively small superficial velocity. A “cooling tower” flow regime was formed due to the existence of periodic large vortex. In test case A.2, the gas phase was injected into the cylinder column at high superficial velocities. The flow fields were highly transient and full of vortexes. In test case A.3, the gas phase was injected into a sudden enlargement pipe and a steady-state stagnant vortex was formed at the pipe bottom. In experiments, the turbulence kinetic energy was monitored, which can be used to validate the turbulence model employed in the E-E method.

Fig. 4 show the horizontal liquid velocity predicted by different force closure combination for test case A.1. It can be seen that all the force closure combination can predict the periodically “Gulf-stream” phenomenon, which proves that the periodically vortex in the bubble column is not sensitive to the force closure. Table 5 shows the predicted gas holdup and the plume oscillation period (POP). Compared with the experimental data, the POP predicted by all these force closure combination was under-estimated. Results can be improved by adjusting the turbulence model and inclusion of the bubble induced turbulence as was shown in our previous work [30]. On the other hand, all these three force closure combination is able to predict good results of the gas holdup with errors smaller than 11%. When the turbulence dispersion force and wall forces were included, the prediction of the gas holdup was slightly improved. However, the addition of wall forces and turbulence dispersion forces cannot improve the results of POP. It can be explained by the fact that the phase shear rate in this test case is too small to present differences, and the occurence of POP (or the pediodically vortex) is mainly because of the drag.

Fig. 5 shows the averaged phase fraction predicted by different force closure combinations for test case A.2. It can be seen that the phase fraction is highly dependent on the employed force models and the difference is quite large. If only the drag force was included, the predicted phase fraction was accumulated at the center of the column since no lateral forces and diffusion force were used. If the turbulence dispersion force was included, the bubbles at the center of the column were diffused along the gradient of the phase fraction and the bubble accumulating problem was prevented. Although the predicted phase fraction was flattened, it was much better than the predicted results when only the drag force was included. Moreover, the flattened phase fraction can be handled by using a smaller the turbulence dispersion force coefficients. The addition of the wall forces cannot improve the results. Instead, the wall forces predict a peak in the vicinity of the wall, which was not observed in the experiments.

Fig. 6 shows the predicted mean axial liquid velocities and the turbulent kinetic energy for test case A.3. It can be seen that the prediction of the axial liquid velocity agree well with experiments when the drag force and turbulence dispersion force were included. However, including the wall forces cannot predict reasonable results. Meanwhile, accurate prediction of the turbulent kinetic energy was proven to be quite difficult. Although the combination of the drag force and turbulence dispersion force improves the results compared with that predicted by addition of the wall forces, the predicted turbulent kinetic energy was under-estimated. Further research is needed to quantify these errors and possibly employ a bubble induced turbulence model to correct the deficiencies [47].

4.2 Test case B.1 - B.4

In this section, the numerical predicted results for test case B.1 to B.4 were investigated. The available experimental data and features are listed in Table 6. These test cases are common in nuclear engineering. In test case B.1 to B.3, the aspect ratios of the geometry are quite large, which implies that there is enough time for the bubbles and the phase fraction to develop along the vertical pipe direction. Pseudo-steady-state can also be reached due to small gas phase fraction and low gas holdup. In test case B.4, non-regular internal pipes are placed in the computational domain and the aspect ratio is relatively small. Meanwhile, it was found in the experiments that the flow field in the pipe is transient due to the existence of the non-regular internal pipes. These test cases can be distinguished by different features. In our preliminary investigations as mentioned previously, we found that the non-drag forces are crucial to obtain reasonable radial phase fraction distribution, especially for the pseudo-steady-state test cases (e.g., case B.1 to B.3). Therefore, in the following, we will turn our attention to the non-drag forces.

Fig. 7 shows the comparison of the upward gas velocity and phase fraction compared with experimental data. It can be seen that the bubble upward velocity is not sensitive to the force closure. These force combination predict similar results. Accurate predictions of the upward gas velocity can be also obtained if only the drag force was included. However, reasonable phase fraction can only be predicted when all the forces was included. Meanwhile, we found that the predicted phase fraction was highly dependent on the wall lubrication force model. For the wall lubrication model developed by Antal et al. [25], using a smaller $C_{2}$ improved the results. Otherwise many bubbles are pushed away when a large $C_{2}$ was used.

Fig. 8 shows the radial profiles of the predicted phase fraction for test case B.2. The difference between test case B.1 and B.2 is that a double peak was observed in the latter one due to the existence of bubbles with different diameters. In the mono-disperse plot, the phase fraction of the small bubbles ($d<5.6$ mm) was reported. All these small bubbles tend to move towards the wall and the wall peak was observed. As mentioned previously, reasonable phase fraction can only be predicted when all the forces was included. By using a smaller $C_{2}$, the wall peak of the phase fraction can be improved. However, none of these models predicted the central peak. It comes from the drawback of the E-E method. As a mono-disperse mathematical model, it can only be used for mono-disperse test cases. Improvements can be obtained by using a higher-order moments methods as it was done in our previous works [48, 49]. Fig. 9 shows the radial profiles of the predicted phase fraction for test case B.3. In this test case, a large bubble diameter was used. Therefore, they move towards the pipe center due to negative lift force coefficients. It can be succeffully predicted by the lift force model developed by Tomiyama et al. [24]. Again, reasonable phase fraction can only be predicted when all the forces was included.

Table 7 reported the predicted gas holdup compared with the experimental data for test case B.4. It can be seen that the predictions were slightly under-estimated. Results may be improved by taking polydispersity into account as was done in the work of Besagni et al. [44]. However, even the gas holdup was under-estimated, it can be seen that the prediction was improved when all the forces were included. If only the drag force was included, it results in the largest error.

4.3 Test case C.1 - C.2

Test case C.1 and C.2 are common in metallurgical industry. The available experimental data and features are listed in Table 8. In test case C.1, the gas phase is injected into the mold by a submerged entrance nozzle (SEN). Large bubbles are lifted toward the meniscus due to the buoyancy force acting on them and removed from the mold, which smaller bubbles are carried deep into the mold. These small bubbles are trapped in the steel and cause pin hold defects. In the field of numerical simulation, this is also called phase segregation which is very difficult to address. In test case C.2, the gas phase was injected to remove the non-metallic inclusions in metallurgical reactor. This test case is important because it is the only one for which the turbulent kinetic energy was reported in the experiments.

Fig. 10 shows the predicted axial and radial liquid velocities. It can be seen that all the force closure combination can predict accurate axial liquid velocity compared with experimental data. The predicted axial liquid velocity agree well with the experimental data. On the other hand, none of the models can predict satisfactory radial liquid velocity. Such inconsistency may come from the mono-disperse assumption of the E-E method and can be handled by using a poly-disperse model. The predicted racial liquid velocities at the right tail of the plot were smaller than 0, which implies the liquid is moving toward the bottom of the mold and a large vortex at the right bottom of the mold is formed. Using a larger bubble diameter improves the predictions due to the effect of large buoyancy. Similar unrealistically predictions for low gas flow rate were also found by Liu and Li [50], in which the flow fields were investigated by large eddy simulation.

Fig. 11 shows the predicted axial liquid velocities, phase fraction and the turbulent kinetic energy. It can be seen that the predicted axial liquid velocity agree well with experiments if the drag force and turbulence dispersion force were included. Predictions became worse if the turbulence dispersion force was neglected. Including the lift and wall forces cannot improve the results, which implies their effects can be omitted. On the contrary with the test cases investigated previously, when only the drag force was included, the predicted turbulence kinetic energy was better than that predicted by all forces combination. However, the predicted phase fraction was seriously over-estimated, which implies the bubbles were not diffused and accumulated at the center line of the reactor. Fig. 12 shows the predicted phase fraction compared with experiments investigated by Castillejos et al. usting a similar equipment [51]. The bubble flow rate is 257 N cm3/s. It can be seen that the predicted results agree well with the experiments if the drag force and turbulence dispersion force were included.

5 Conclusions

In this work, we employed different force combinations to simulate industrial bubbly flows. Instead of investigating the effects of momentum interfacial exchange force by a specific test case, nine different test cases were employed. Our aim is to seek a general numerical settings for the Eulerian-Eulerian model, which can provide competent results for industrial bubbly flow simulations. These test cases were selected from different industries including chemical, nuclear, bio-processing and metallurgical engineering. Simulations were launched by the OpenFOAM solver reactingTwoPhaseEulerFoam. The drag force developed by Ishii and Zuber [31], turbulent dispersion force developed by Lopez De Bertodano [26], lift force developed by Tomiyama et al. [24] and wall force developed by Antal et al. [25] were employed. Predictions were compared against experimental data. It was found that the drag force and turbulent dispersion force play the most important role on the predictions and should be included for all simulations. Otherwise the bubbles tend to accumulate since the spreading effect of the lift force is weak. For the pipe flows with large aspect ratio, the lift force and wall lubrication force should be included to address the phase fraction accumulation in the vicinity of the wall. In other cases these lateral forces can be safely neglected.

At last, the test case library is open-sourced and are available as supplementary data for anyone to download as baseline test cases for further investigations.

Contributions

Yulong Li: wrote the paper; post-processed the data; performed the analysis. Dongyue Li: conceived the analysis; corrected the paper; prepared the OpenFOAM test cases.

Acknowledgement

One of the author (Dongyue Li) wants to acknowledge the CFD-China community for the fruitful discussions of gas-liquid flows.

Appendix

The finite volume discretization for the E-E method was implemented in OpenFOAM-6 by the OpenFOAM foundation with novel techniques to ensure stabilities. To the authors’ knowledge, these techniques were not published. Therefore, we summarize the procedures in this section. After omitting the buoyancy and pressure terms, Eq. (1) and (2) can be written as follows:

[TABLE]

[TABLE]

where

[TABLE]

The discretized form of Eq. (9) and (10) can be written as follows:

[TABLE]

[TABLE]

where the subscript P denotes the owner cell, the subscript N denotes the neighbour cells, $A_{\rm P}$ and $A_{\rm N}$ are the matrix diagonal coefficients contributed by the owner cells and neighbour cells, $S$ are source terms. The predicted velocities can be obtained by solving the discretized form of Eq. (12) and (13):

[TABLE]

[TABLE]

where ${\mathbf{HbyA}_{{\rm a},{\rm P}}}$ and ${\mathbf{HbyA}_{{\rm b},{\rm P}}}$ are the predicted velocities for phase a and phase b, respectively. At this stage, the predicted velocity fields are not divergence free and a pressure equation needs to be constructed. In order to simplify the designation of the pressure boundary conditions and to reduce the spurious currents caused by hydrostatic pressure in the non-orthogonal grids, the pressure term and buoyancy term are usually combined together, and the pressure without the hydrostatic part, $p_{\rm rgh}=p-\rm\rho\mathbf{g}\cdot\mathbf{h}$, is used. The gradient of $p_{\rm rgh}$ can be calculated by

[TABLE]

Substituting the definaiton of $p_{\rm rgh}$ and Eq. (16) into the buoyance and pressure terms leads to:

[TABLE]

[TABLE]

In this manner, a $p_{\rm rgh}$ equation can be constructed.

For incompressible flows, summarizing Eq. (1) and Eq. (2) produces the divergence constraint equation:

[TABLE]

The discretized form of Eq. (19) can be written as follows:

[TABLE]

where the subscript f implies variables defined at the cell faces. Eq. (20) is usually employed as a restrictive condition to construct pressure Poisson equation. Since the pressure gradient and buoyancy term are not considered in Eq. (14) and (15), adding these terms produces the predicted velocities:

[TABLE]

[TABLE]

Substituting the interpolated cell face velocities into Eq. (20), the phase fluxes can be written as follows:

[TABLE]

where

[TABLE]

[TABLE]

Eq. (23) can be used to calculate the $p_{\mathrm{rgh}}$ in the PISO loops [52]. After the pressure was updated, a divergence-free velocities can be updated.

When the drag coefficient $\mathbf{K}\mathbf{d}$ is much larger than the diagonal coefficient $\mathbf{A}$ at certain cells, the classical semi-implicit algorithm discussed previously leads to very large relative velocity Courant number and the time step will be very small. In OpenFOAM-6, a robust approach was implemented by Weller [53] to handle the pressure-velocity coupling and it was not published. In this method, the velocities without contributions of the drag forces are defined as follows:

[TABLE]

[TABLE]

The phase velocities can be written as

[TABLE]

[TABLE]

[TABLE]

where $\rm D_{\rm b,P}=\frac{\rm Kd_{p}}{\mathbf{A}_{\rm b,p}}$, $\rm D_{\rm a,P}=\frac{\rm Kd_{p}}{\mathbf{A}_{\rm a,p}}$. The relative velocity is defined by

[TABLE]

Subtracting $\rm D_{\rm b,P}\rm D_{\rm a,P}(\mathbf{U}_{\rm b,P}-\mathbf{U}_{\rm a,P})$ from the L.H.S. and the R.H.S. of Eq. (31) leads to:

[TABLE]

After arrangement, Eq. (32) can be written as follows:

[TABLE]

Once the relative velocity $\mathbf{U}_{\rm P}^{\rm r}$ is updated, the velocity of discrete phase and continuous phase can be calculated by the following equations, respectively:

[TABLE]

[TABLE]

For the flows with strong swirl or strong body forces, the algorithm discussed above tends to produce phase fraction oscillations [54]. This issue can be mitigated by including these body forces on the cell faces to construct the face-based pressure equation. In this manner, Eq. (9) and Eq. (10) can be re-written as follows:

[TABLE]

[TABLE]

It can be seen that in Eq. (36) and Eq. (37) the contribution of the non-drag forces is neglected. The contribution needs to be considered to update the predicted velocities. Therefore, Eq. (21) and Eq. (22) can be re-written as follows:

[TABLE]

[TABLE]

The corresponding fluxes can be written as follows:

[TABLE]

[TABLE]

Substituting $\phi_{\rm a,f}$ and $\phi_{\rm b,f}$ into Eq. (23) leads to the pressure Poisson equation which can avoid the oscillating results due to strong body forces.