**About the Editors**

**Li Xi** obtained his Ph.D. from the University of Wisconsin-Madison. After completing his postdoctoral training at the Massachusetts Institute of Technology, he joined the Department of Chemical Engineering at McMaster University as a faculty member. He has broad research interests in chemical and materials engineering, with recent focuses on the modeling and analysis of chemical and pharmaceutical reactors and the development of sustainable polymer production processes. In addition, he also actively contributes to fundamental research in turbulent flows, polymer physics, and rheology.

**De-Wei Yin** is a research engineer at the Dow Chemical Company in Midland, Michigan, U.S.A. As a member of the Fluid Mechanics & Mixing Discipline in the Core R&D Engineering & Process Sciences Laboratory at Dow, he has conducted research on porous media flows, multiphase stirred-tank reactor design, flow measurement techniques, characterization and optimization of mixing processes, numerical algorithm development, and the application of dimensional analysis in design of experiments. He is an active member of the American Institute of Chemical Engineers and of the North American Mixing Forum.

**Jae Sung Park** is currently an assistant professor in the Department of Mechanical and Materials Engineering at the University of Nebraska-Lincoln. He was a postdoctoral associate in the Department of Chemical and Biological Engineering at the University of Wisconsin-Madison. He received his Ph.D. in Mechanical Engineering from the University of Illinois at Urbana-Champaign. His research interests encompass a wide range of transport phenomena, including nano/micro/biofluids, electrokinetics, turbulent flows, and flow control. He has developed advanced computational algorithms to investigate various applications in fluid mechanics.

### *Editorial* **Special Issue "CFD Modeling of Complex Chemical Processes: Multiscale and Multiphysics Challenges"**

**Li Xi 1,2,\*, De-Wei Yin <sup>3</sup> and Jae Sung Park <sup>4</sup>**

	- <sup>3</sup> The Dow Chemical Company, Midland, MI 48667-1776, USA; dwyin@dow.com
	- <sup>4</sup> Department of Mechanical and Materials Engineering, University of Nebraska-Lincoln, Lincoln, NE 68588-0526, USA; jaesung.park@unl.edu

After decades of development, computational fluid dynamics (CFD), which solves fluid mechanics and, more generally, transport phenomena problems using numerical analysis, has become a main-stream tool in many areas of engineering practice. In chemical processes, CFD is widely used not only in predicting fluid flow and mixing patterns in chemical reactors, but also in the simulation and modeling of various unit operations where complex flow and transport processes make process outcomes hard to predict and control. Although CFD models typically still require experimental validation, their application brings many advantages in process research and development. First, numerical solutions from CFD contain fully detailed information about the flow field and distribution profiles of other quantities of interest (temperature, pressure, chemical species concentration, etc.), which are often hard to measure directly from experiments, especially in situ. The availability of such detailed data on the fluid dynamics is often critical to process analysis and improvement. Second, a carefully validated CFD model can be used to predict process outcomes in new operating conditions and process design, with which untested process parameters can be prescreened. As such, the number of experiments and trials and errors can be greatly reduced. This not only reduces the cost in, e.g., process design, adaption, and scale-up, but also reduces the risk of potential safety, environmental, and health hazards associated with many chemical processes.

With the increasing power of high-performance computing (HPC) facilities and continued development of numerical techniques, the accuracy and sophistication of CFD models have improved. Obtaining reliable and practically useful predictions, however, is far from a simple click of a button, owing to the intrinsic complexities in many chemical processes. Most chemical processes involve complex coupling of fluid flow with other physicochemical processes, such as heat and mass transfer, chemical reactions, and particle/bubble/droplet dynamics. In addition, dynamics over vastly different length and time scales also interact in nontrivial ways. For one thing, many chemical processes involve flow turbulence, which is intrinsically multiscale. For another, multiphase flow processes where phase separation and fluid flow occur at different scales are also common. Overcoming multiscale and multiphysics challenges is thus often a central theme and requires innovative modeling approaches that best suit the particular systems.

This special issue "CFD Modeling of Complex Chemical Processes: Multiscale and Multiphysics Challenges" highlights some recent advances in the application of CFD in chemical processes. We take a very liberal approach and define chemical processes as those where principles of chemical engineering, such as transport phenomena, are important, even though the direct application may not fall within the traditional scope of the chemical industry. A total of 16 papers are collected, including one review, 14 original research articles, and one correction. The topics range from pure methodology development to the application of existing CFD tools in practical systems, and everything in between.

**Citation:** Xi, L.; Yin, D.-W.; Park, J.S. Special Issue "CFD Modeling of Complex Chemical Processes: Multiscale and Multiphysics Challenges". *Processes* **2021**, *9*, 775. https://doi.org/10.3390/pr9050775

Received: 15 April 2021 Accepted: 16 April 2021 Published: 28 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The review paper by Nakhaei et al. [1] offered a comprehensive and in-depth account of CFD application in gas–solid cyclone separators. It covered the fundamentals of gas– solid cyclone separator systems including their key performance metrics and how they are influenced by operating conditions. This was followed by an extensive review of available approaches for the numerical modeling of such systems and a summary of existing CFD studies at both ambient and elevated temperatures. Capabilities and limitations of existing approaches were also discussed.

Among the original research articles, Vachaparambil and Einarsrud [2] had a focus on modeling methodology. Based on simple benchmark problems, the study compared different surface tension models used in the volume of fluids method for gas–liquid twophase flows. Special attention was given to the occurrence of spurious currents due to model inaccuracies.

Other studies all targeted more complex and practical engineering systems. Several studies applied CFD to different types of chemical and biochemical reactors. Ebrahimi et al. [3] studied the flow and mixing performance of double impellers in a stirred-tank reactor. Effects of impeller configuration and speed were reported. Han et al. [4] studied the combustion performance in a pulverized coal furnace with burners installed on the front wall by coupling the momentum balance (with turbulence model), species transport, and energy balance equations. The goal was to optimize the process for reduced nitrogen emission and air pollution. Reactors with multiphase flow were also represented. Khezri et al. [5] studied a bubbling fluidized-bed reactor for biomass gasification. The gas–solid flow was modeled with an Eulerian–Eulerian approach and various turbulence and drag models were compared. Effects of air distributors with different pore sizes were also examined. Tao et al. [6] again used an Eulerian–Eulerian two-fluid model to study a large-scale high-pressure gas–liquid bubble column reactor, with a focus on the effects of operating parameters on the gas holdup of bubbles of different sizes. An Eulerian–Eulerian two-fluid model was also used in Sarkizi Shams Hajian et al. [7] for the gas–liquid two-phase system in the production process of baker's yeast at the industrial scale. The study combined turbulence modeling with species transport and bio-kinetic models to analyze species distribution in the reactor.

CFD studies of other unit operations are also included. Wang, Z. et al. [8] studied a rotating packed bed process for gas desulfurization, where a liquid solution was dispersed in the gas phase for the selective absorption of H2S. An Eulerian–Lagrangian approach was used where the dynamics of liquid droplets were described by coupling force balance with kinetic models for their coalescence and breakup, with which characteristics of droplets were studied. Qadir et al. [9] numerically simulated the gas separation process in membrane modules and studied the effects of membrane configuration and operating parameters. Landauer and Foerst [10] experimentally studied the triboelectric separation of starch and protein particles. After passing through a charging tube, particles carrying different triboelectric charges settled to different distances in a separation chamber under an electrical field. Effects of particle size and material on the settling distance were studied. The gas flow field in the chamber was modeled with CFD, and simulated particle trajectories in the flow and electrical fields were used to estimate the minimum charge for particle separation. Along the lines of particle processes, Yang et al. [11] used an Eulerian– Lagrangian discrete particle model to simulate the pneumatic conveying of coal particles, studied the effects of operating conditions on the pressure drop, and identified the range of model applicability by comparison with experiments.

Applications in transport processes not in the traditionally-defined scope of chemical engineering, which nonetheless were built on the same principles, are also represented in this special issue. Gao et al. [12] used CFD to model the hydrodynamics in the lateral inlet/outlet part of a pumped hydroelectric storage system. Surrogate models were then built based on CFD results and used for the geometric optimization of the device. Wang, W. et al. [13] measured the air flow, temperature, and relative humidity in the blind heading of an underground mine in situ to study the heat emission by equipment. CFD was conducted to obtain the detailed flow and temperature fields and to optimize the design

of the ventilation system. Also in mining, Zhou et al. [14] numerically solved the air flow and species transport equations in a typical stone-coal mine laneway to study the effects of forced ventilation on the radon concentration. Numerical analysis of coupled fluid flow and transport processes was also the theme of Mousavi et al. [15], where the air flow and heat transfer in a novel sustainable farming compartment (SFC) was modeled. This SFC was designed with an evaporative cooling system where the vaporization of treated wastewater was used to absorb heat and the cooling effect was distributed across the domain by forced convection. An experimental prototype was also set up to demonstrate the concept, and CFD was used to optimize the system design and operating conditions.

Finally, we would like to express our sincere gratitude towards all authors for contributing to this special issue. We look forward to the continued development of CFD methodology and its application in chemical process engineering in years to come.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Review* **CFD Modeling of Gas–Solid Cyclone Separators at Ambient and Elevated Temperatures**

#### **Mohammadhadi Nakhaei 1,2,\*, Bona Lu 3,4, Yujie Tian 3,4, Wei Wang 3,4, Kim Dam-Johansen <sup>1</sup> and Hao Wu <sup>1</sup>**


Received: 15 January 2020; Accepted: 7 February 2020; Published: 15 February 2020

**Abstract:** Gas–solid cyclone separators are widely utilized in many industrial applications and usually involve complex multi-physics of gas–solid flow and heat transfer. In recent years, there has been a progressive interest in the application of computational fluid dynamics (CFD) to understand the gas–solid flow behavior of cyclones and predict their performance. In this paper, a review of the existing CFD studies of cyclone separators, operating in a wide range of solids loadings and at ambient and elevated temperatures, is presented. In the first part, a brief background on the important performance parameters of cyclones, namely pressure drop and separation efficiency, as well as how they are affected by the solids loading and operating temperature, is described. This is followed by a summary of the existing CFD simulation studies of cyclones at ambient temperature, with an emphasis on the high mass loading of particles, and at elevated temperatures. The capabilities as well as the challenges and limitations of the existing CFD approaches in predicting the performance of cyclones operating in such conditions are evaluated. Finally, an outlook on the prospects of CFD simulation of cyclone separators is provided.

**Keywords:** gas–solid; cyclone separator; computational fluid dynamics; elevated temperature process

#### **1. Introduction**

Gas–solid cyclones are frequently used in industrial processes with the primary purpose of two–phase flow separation, i.e., separation of a high-density phase from a lower-density carrier phase, using a turbulent swirling flow. The state-of-the-art industrial cyclone designs are able to operate at elevated temperatures and moderate-to-high loading of solids, while at the same time meeting the required separation efficiency and having low investment and maintenance costs. This has led to the frequent use of cyclones as the only/initial stage of separation and cleaning processes rather than other industrial separators, e.g., bag filters, electrostatic separators, etc. Examples of cyclone separator applications in demanding industrial process conditions are high temperature gas–solid heat exchangers [1], e.g., in the cement industry; gasification and combustion of solid fuels [2–4]; coal pyrolysis [5]; and gas–solid separation in circulating fluidized beds (CFBs) [6].

With the main purpose of gaining improved insight into the flow physics and factors affecting the two important performance parameters of cyclones, i.e., pressure drop and separation efficiency, single–phase and gas–solid flows in pilot-scale cyclones have been extensively studied with the use of experimental methods (for example, see [7–10]). Accordingly, a number of simple semi-empirical/-theoretical models have been proposed to address the flow field and performance

of cyclones (for example see [11–16]), and some of them are still being used in cyclone design and optimization.

In the past two decades, computational fluid dynamics (CFD) has been broadly applied to predict the separation efficiency and pressure drop of cyclone separators and to optimize cyclone designs. Compared to the models based on classical cyclone theory, three-dimensional CFD simulations have the advantage of taking into account the unsteadiness and asymmetry of cyclone flow. On the other hand, the presence of strong swirl and anisotropic turbulent flow as well as adverse pressure gradients in the cyclones has driven the CFD simulation studies to use more advanced turbulence models, e.g., the Reynolds stress transport model (RSTM) and large eddy simulation (LES), as well as higher-order discretization techniques, that are capable of capturing these specific flow physics. These models/techniques, however, are computationally demanding compared to the more commonly used models, e.g., k–*-*. Furthermore, with the addition of extra physics to the CFD simulation of cyclones, e.g., presence of particles at high loadings and gas–solid heat transfer, additional complications will arise with respect to the solution's accuracy and stability.

In this paper, a brief introduction to the basic operation principles of cyclones at ambient and elevated temperatures is initially presented. Subsequently, the existing CFD simulation studies of cyclones are summarized and discussed based on the operating temperature, i.e., ambient and elevated temperature, and with a special focus on the studies with moderate-to-high loading of particles. The paper highlights specific process parameters that are able to be captured by the CFD simulations at these temperatures. Furthermore, the important sub-models utilized as well as specific challenges that may be encountered in such CFD simulations are addressed.

#### **2. Fundamentals of Gas–Solid Cyclone Separators**

The reverse-flow-type gas–solid cyclones are usually equipped with either tangential inlets, which are the main focus in this study, or axial inlets (for example, see [17,18]), usually referred to as swirl tubes [19]. A common configuration of a tangential inlet reverse–flow cyclone is schematically shown in Figure 1a. In this type of cyclone, the swirling gas flow pattern usually consists of a "double vortex": an outer vortex with a downward axial velocity and an upward-moving inner vortex. The tangential flow pattern, caused by the tangential inlet of the gas to the cyclone system, is referred to as a "Rankine vortex", comprised of a near-solid-body rotating flow at the core (i.e., tangential velocity is directly proportional to the radius) and a loss-free vortex flow (i.e., tangential velocity is inversely proportional to the radius) at the walls [19] (see Figure 1b). The outer vortex gradually loses its downward momentum in the bottom regions of the cyclone and changes its direction, i.e., "vortex end", to form the inner vortex. The inner vortex is led to the exit pipe, which is usually extended to the cyclone body to separate the inner vortex from the inlet velocity field. The axial distance from the "vortex end" location to the vortex finder is termed the "natural turning length" [19]. The double vortex structure is inherently unstable due to the presence of a radial pressure gradient imposed by the vortex itself [20]. This phenomenon is referred to as a precessing vortex core (PVC), which causes the location of the vortex axis to oscillate with a specific frequency. In some cyclones, in order to stabilize the vortex end position, a vortex stabilizer is installed either under or above the dust exit location (see Figure 1b), which leads to improved separation efficiency and desirable static pressure below the stabilizer [19].

The static pressure in the cyclone increases monotonically toward the walls to maintain the equilibrium of the rotating flow pattern. This pressure gradient is also present in the boundary layer of the cyclone walls, where the tangential velocity is small [19]. This leads to the presence of a secondary flow, namely inwardly-directed gas flow along the walls, at the cyclone roof, and the conical section, as depicted in Figure 1c.

**Figure 1.** (**a**) Schematic drawing of a conical reverse-flow cyclone separator illustrating the basic operating principle and the presence of a double vortex inside the cyclone [21], reproduced with permission from G. Towler and R. Sinnott, Specification and design of solids-handling equipment, published by Elsevier, 2013. (**b**) Qualitative patterns of axial, tangential, and radial velocity components of the gas-flow field in cyclones (right) [22], reproduced with permission from M. Trefz and E. Muschelknautz, Extended cyclone theory for gas flows with high solids concentrations, published by John Wiley and Sons, 1993. (**c**) The secondary flow pattern caused by the swirl and pressure gradients in the cyclone [19], reproduced with permission from A. Hoffmann and L. Stein, Gas Cyclones and Swirl Tubes: Principles, Design and Operation; published by Springer Nature, 2007.

#### *2.1. Performance Parameters of Cyclones*

For the main purpose of gas–solid separation, two important design parameters of reversed-flow conical cyclones are separation efficiency and pressure drop. These parameters are affected by the gas–solid flow field inside the cyclone, which in turn is influenced by the cyclone's geometrical features as well as the operating conditions. In this section, these performance parameters are discussed and the influences of operating conditions, e.g., solids loading and gas temperature, on cyclone performance are explained.

#### 2.1.1. Separation Efficiency

Cyclone separators are suitable for the separation of solid particles with a size range of 2–2000 microns, which are found in many industries, e.g., heavy industrial smoke, coal dust, cement dust, etc. [19] When solid particles are fed to a cyclone, they are affected by two main forces: the radially outward-directed centrifugal acceleration force, proportional to the cube of particle diameter, and the fluid drag force applied to the particles in the opposite direction; whereas the gravity force is reported to be of minor importance [23]. When the Stokes law applies for the drag force, which is often valid for solid particles in cyclones, the drag force is proportional to the particle diameter. This indicates the dominance of the centrifugal acceleration force for larger particles, leading to their improved separation. The larger particles experience the centrifugal effect as soon as the gas and solids experience the rotational flow at the inlet, and subsequently, they are pushed toward the walls and lose their momentum. Once the particles approach the wall, they start to move downward due to gravity and the drag force applied to them by the downward-directed gas flow. Finally, these particles are separated at the bottom of the cyclone.

Conversely, fine particles, e.g., 1–10 microns, have a greater tendency to follow the gas flow in a cyclone. These particles are poorly separated from the gas phase due to either following the bypass gas to the vortex finder (see Figure 1c) or entrainment to the inner vortex along the inner and outer vortex boundary. The turbulence dispersion of particles intensifies the re-entrainment of fine particles to the inner vortex followed by particle carry over to the vortex finder. On the other hand, the agglomeration of small particles as well as particle attrition can affect the separation efficiency. For example, in the

reported grade efficiencies of a dilutely loaded pilot-scale cyclone, the separation efficiency of particles with a diameter of 0.8–1.0 μm is smaller than the separation efficiency of very fine particles, e.g., 0.3 μm, as well as of large particles [19]. This is known as "fish-hook" behavior [24] and is mainly due to agglomeration of small-sized particles to larger particles, leading to the improved grade efficiency of these particles. The mentioned trend is also observed in other experimental studies [25].

#### 2.1.2. Pressure Loss

One of the main parameters considered in the design of industrial cyclone systems is the energy loss, usually termed the pressure loss. The smaller the pressure loss of the cyclone in a process, the cheaper the process cost. The overall pressure loss in a reversed-flow conical cyclone can be separated into three parts [19,20,26]:


The first contributor to the pressure loss is usually of minor importance. The pressure loss in the separation zone is mainly due to the frictional losses at the wall of the cyclone body. The vortex finder pressure loss is due to the dissipation of the swirl dynamic pressure in the vortex finder, which usually takes place without recovering the dynamic pressure into the static pressure. This pressure drop is the main contributor to the overall pressure drop for single-phase flow or dilutely loaded cyclones, and it is proportional to the square of the tangential velocity magnitude. The vortex finder pressure loss is usually affected by the frictional pressure loss, e.g., if the frictional pressure increases due to higher wall roughness or other operating conditions, the vortex will be weakened, leading to the reduction of the vortex finder pressure loss. This indicates that the contribution of different parts of pressure loss changes once the operating conditions of the cyclone change, e.g., an increase in the temperature or the solids loading, which will be further discussed in subsequent sections.

#### *2.2. Effect of Solids Loading on Cyclone Performance*

At low solid-to-gas mass loading ratios, e.g., 0.01 kg*s*/kg*g*, solid particles are clustered as thin strands on the cyclone walls and are transported in the downward direction in a spiral path, while at high mass loadings, e.g., 10 kg*s*/kg*g*, a major portion of (or the whole) wall surface area is covered with a layer of solid particles (also know as a dense strand in [27]), which is slipped directly into the solids discharge [22]. According to the existing experimental studies of cyclones, both overall and grade separation efficiencies are improved by increasing the inlet solids mass loading ratio [7,8,19]. Hoffmann and Stein [19] reported an overall improvement of the grade efficiency of almost all particle sizes for an increase in the inlet particle mass loading from 3.7×10−<sup>3</sup> to 2.6×10−<sup>2</sup> kg*s*/kg*g*. However, they have stated that the separation efficiency at high mass loading of particles is somewhat dependent on the inlet gas velocity of the cyclone, i.e., the swirl strength [28]. The overall separation efficiency versus inlet solids loading of selected experimental studies is presented in Figure 2. In this figure, the improvement of separation efficiency with the cyclone Reynolds number can be observed from the experimental data of [19] for a fixed value of kinematic response time of particles of *<sup>τ</sup><sup>p</sup>* <sup>=</sup> *<sup>ρ</sup>pd*<sup>2</sup> *p* <sup>18</sup> <sup>μ</sup>*<sup>g</sup>* = 0.1 milliseconds. The improvement is more noticeable for lower mass loadings, i.e., below 0.01 kg*s*/kg*g*. The kinematic response time of particles used in the study of Fassani and Goldstein [29] is somewhat higher compared to the particles of other studies presented in Figure 2, leading to a very high efficiency of separation. For the rest of the studies, no conclusive argument can be stated regarding the effect of this parameter.

Despite available experimental studies on the effect of solids loading on the separation efficiency, there is no consensus regarding the exact mechanism of this improvement [19]. According to the concept of the critical mass loading ratio, *φG*, initially introduced by Muschelknautz and Brunner [8,13], when the inlet particle mass loading is higher than the critical mass loading, the share of particles corresponding to the extra amount of mass loading is ideally separated at the cyclone inlet, while the rest of the particles are separated based on a balance between the centrifugal forces and the turbulent dispersion, commonly referred to as "inner separation" [22]. The concept of Muschelknautz, however, suggests that the separation efficiency is independent of the mass loading for loadings smaller than *φG*, which is not consistent with the experimental data in the literature [28]. On the other hand, based on another concept introduced by Mothes and Löffler [30], the improvement in the separation efficiency at high solid mass loadings is attributed to the agglomeration of smaller particles, i.e., 2 μm diameter and smaller, to larger particles, i.e., 15 μm diameter. Hoffmann and Stein [19] interpreted the separation efficiency of cyclones at high mass loadings to be a combined effect of different contributions. For very fine particles, the agglomeration effect is more dominant, leading to a high efficiency of these particles, especially for a humid carrier gas [19]. Furthermore, they argued that the increased concentration of particles at the inlet can lead to reduced drag force, which improves the separation process.

**Figure 2.** Overall cyclone separation efficiency versus mass loading ratio at the cyclone inlet for selected experimental data in the literature [7,8,19,29,31]. The Reynolds number is calculated based on the inlet velocity and the cyclone body diameter.

Baskakov et al. [32] and Muschelknautz and Brunner [8] have reported a decrease in the pressure drop up to a specific amount of solids loading, followed by an increase. A similar trend is also observed in other experimental [33] and CFD simulation studies [34,35]. This behavior is attributed to the increase of the frictional pressure loss due to the presence of particles near the walls. Although, at the same time as the friction velocity at the cyclone walls increases, the swirl intensity is weakened; consequently, the pressure loss associated with the vortex finder is decreased. Muschelknautz and Brunner [8] reported an 85% reduction in the tangential velocity peak by increasing the solids loading from 0 to around 20 kg*s*/kg*g*. At low mass loadings, the decrease in the vortex finder pressure loss is more pronounced, leading to the reduction of the overall pressure loss, whereas when the particle mass loading is sufficiently high, the vortex finder pressure loss becomes negligible due to a very weak tangential swirl and the frictional pressure drop still increases with the loading of particles [35]. As a consequence, the overall pressure drop starts to increase with the solids loading after reaching a minimum value.

A summary of selected experimental data on pressure drop versus solids loading in pilot-scale cyclones from the literature is presented in Figure 3. Despite the mentioned explanation, there are other studies in the literature that dispute the presence of a minimum pressure drop when the loading of particles changes. For example, in some studies [7,9,10,31], it is reported that by increasing the mass loading of particles in the cyclone, the pressure drop is reduced, although these studies are limited to small mass loadings, e.g., up to 0.5 kg*s*/kg*g* in [31], 0.1 kg*s*/kg*g* in [7], and 0.2 kg*s*/kg*g* in [9]. The mass loading ratio at which the minimum pressure drop occurs is reported to be around 0.25 (at an operating gas temperature of 250 ◦C) and 2.0 (at ambient temperature) in the studies of Baskakov et al. [32] and Muschelknautz and Brunner [8], respectively. Furthermore, in some of the studies, it is reported that by increasing the inlet solids loading, the pressure drop increases [36] or remains at a nearly constant value lower than the pressure drop of the single-phase flow cyclone [29].

**Figure 3.** Experimental data of pressure drop (normalized with the pressure drop of a particle-free cyclone, Δ*P*0) versus mass loading of pilot-scale cyclones from selected studies [7–9,29,31,32,37]. The lines are numerical fits to each set of experimental data. The Reynolds number is calculated based on the inlet velocity and the cyclone body diameter.

#### *2.3. Effect of Operating Temperature on Cyclone Performance*

The use of gas–solid cyclones is popular in high-temperature industrial processes, e.g., for drying, solidification, heating, flue gas cleaning purposes, etc. The high-temperature cyclone separators can be loaded with large amounts of particles, e.g., in gas–solid CFB cyclones, the solids loading can be on the order of 10–100 kg*s*/kg*g* [38]. As briefly mentioned earlier, the cyclone's performance is affected by its operating temperature. In this section, a brief introduction to the effect of cyclone gas temperature on separation efficiency, pressure drop, and gas–solid heat transfer is provided.

In general, as gas temperature rises, gas density and viscosity are, respectively, decreased and increased. At a fixed volumetric gas-flow rate to the cyclone, a direct consequence of this change in the fluid properties is the reduction of the cyclone's Reynolds number and the swirl intensity in the cyclone at elevated operating temperatures. As a result, the cyclone's pressure drop is reduced. This trend is reported in many experimental [39,40] and CFD studies (for example, see [41]) of cyclones operating at elevated temperatures.

For a fixed inlet gas velocity to the cyclone, as the temperature increases, the particle cut-size, i.e., the diameter at which the cyclone separation efficiency is equivalent to 50%, increases, indicating a negative influence of the operating temperature on the separation efficiency [39,40]. This behavior can be partly attributed to an intensified drag force applied to the particles by the carrier fluid as a consequence of increased gas viscosity at elevated temperatures. Furthermore, as stated before, the reduction in the swirl intensity leads to a weaker centrifugal effect on the particles compared to the intensified drag force, and as a result, the separation efficiency is reduced [40]. On the other hand, the gas density is inversely proportional to the temperature; but its effect on particle separation is less significant as only the difference between the gas and particle densities is the determining factor on the buoyancy force applied to the particles.

In the existing experimental studies in the literature, the heat transfer rate to the particles in gas–solid cyclone heat exchangers is positively affected by the loading of particles and the inlet gas velocity, while being inversely influenced by the particle size [1]. Improvement in the heat-transfer rate with solids loading is attributed to the reduced radial gas velocity in the cyclone as a result of stronger gas–solid momentum coupling. Consequently, the amount of gas bypassed to the vortex finder at the entrance of the cyclone is reduced, leading to the presence of a higher-temperature gas at the bottom regions of the cyclone and further improvement in the heat-transfer rate to the solids [1]. As the

particle loading further increases, the rate of increase in the gas–solid heat transfer rate reduces [1,42] and approaches a certain value. Increased gas–solid heat transfer as a result of higher inlet gas velocity as well as the reduction of particle size can be simply explained by an improved driving factor for gas–solid heat transfer, i.e., more availability of the high temperature gas to transfer heat to the solids and improved gas–solid surface area for heat transfer, respectively.

One of the challenges that may exist at the cyclone's internal walls, especially for cyclones operating at elevated temperatures, is erosion due to wall–particle collision and/or corrosion due to the presence of aggressive species and agents in the gas, even in the presence of a strong refractory lining. It is more likely to have erosion on the cyclone wall in front of the gas inlet, caused by the very first impact that solids have with the cyclone wall, followed by a change in the direction of particle movement [43]. An immediate outcome of cyclone wall erosion is more rough surfaces on a cyclone's internal walls, affecting the pressure drop and separation efficiency of particles. A common solution for this challenge is "hardware" changes in the cyclone system during the maintenance period, e.g., installation of wear plates and replacement of eroded materials [19]. Furthermore, cyclone design modification is an option to improve the erosion behavior, e.g., adding a flat-disk vortex stabilizer close to the particle outlet reduces erosion [44]. Apart from erosion and corrosion, the deposition of particles on the cyclone walls can induce operational challenges in the long term. Deposit formation can take place on the outer surface of a vortex finder [45], dipleg [46], or on the internal walls of the cyclone body [47].

#### **3. Approaches for the Numerical Modeling of Gas–Solid Systems**

The use of numerical simulations is undoubtedly a valuable tool for understanding and improving industrial gas–solid systems. However, the presence of a very large span for length and time scales of both phases, which are influenced by each another, makes resolving all of the details of a gas–solid system through computational simulation a challenging task [48]. In the existing numerical methods, the gas–solid flow is resolved at different but limited ranges of length and time scales, and the smaller-scale details, if any, are modeled. A summary of the available approaches for a four-way coupled (to be explained later) gas–solid system is provided in Table 1. In the most fundamental approach, usually referred to as direct numerical simulation (DNS), the boundary layer over the particles surface is resolved. DNS can itself be categorized into two approaches: the resolved Eulerian–Lagrangian model, e.g., immersed boundary model, and the Lagrangian–Lagrangian approach using molecular dynamics, e.g., Lattice Boltzmann, applicable for simulations of gas–solid flows at extremely small scales of 0.01 and 0.001 m, respectively [48]. DNS is usually utilized for fundamental studies of gas–solid flows at small scales in order to develop sub-models, e.g., drag and Brownian motion models, for more simplified approaches.

In the Eulerian–Eulerian (E–E) approach, both phases are described as a continuum medium and solved in an Eulerian domain [49], while the gas–solid interactions are incorporated through a drag model and the particle–particle interactions through models for solid pressure and viscosity, e.g., the kinetic theory of granular flows (KTGF) [50]. A challenge in using the basic formulation of the E–E approach is the difficulty in taking into account poly-disperse particles in the model. However, there are several sub-models, such as the population balance model (PBM), that can be added to the basic E–E equations to account for different particle sizes, with the cost of higher computational overhead since separate momentum and continuity equations should be solved for each size bin [51–53]. In addition, some complex physical phenomena such as particle agglomeration and break-up can be modeled in combination with the PBM [54]. A limitation of the E–E approach is the validity of the continuum hypothesis for both phases, i.e., both phases are adequately present over the whole computational domain. This makes the model more appropriate for gas–solid systems with more homogeneous distributions of particles, mainly dense gas–solid systems such as fluidized beds.

In the (unresolved) Eulerian–Lagrangian (E–L) approach, also referred to as the discrete particle model (DPM), the carrier gas flow field is considered as a continuum medium and solved using an Eulerian formulation, while solid particles are tracked in a Lagrangian platform using Newton's second law. Compared to the resolved E–L model, categorized as DNS in this study, the boundary layer over the particle surface is not resolved in the unresolved E–L and the size of the computational grid is significantly larger than the particle size. One distinct advantage of this model compared to the E–E approach is the direct resolving of particle movement, i.e., fewer closure models are employed for the particle phase, which makes the model more accurate. Furthermore, poly-disperse systems can be readily modeled using the E–L method. The particle–particle interactions in the E–L model are usually incorporated as either hard-sphere or soft-sphere collision models. In the hard-sphere model, the collisions take place instantly and can be employed in both stochastic and deterministic collision approaches; thus it is suitable for dilute flows where binary collisions happen. Conversely, this model is not suitable for cases where particles go through multiple collisions with other particles simultaneously. The soft-sphere model, also referred to as the discrete element model (DEM), is a deterministic particle collision approach that takes into account overlapping particles during the collision and is able to incorporate multiple simultaneous collisions for a particle [55]. As the loading of particles in a gas–solid system increases, the E–L approach becomes computationally expensive, especially if the DEM model is utilized [55].

As an alternative, the particle–particle interactions can be modeled by using hybrid E–E and E–L approaches, such as the multi-phase particle in-cell (MP-PIC) method [56–58] and the dense discrete phase model (DDPM) [59]. The basic idea in the hybrid approach is to solve the carrier phase in an Eulerian frame with the inclusion of the particle phase volume fraction, taken from the Lagrangian tracking, in the conservation equations. Similar to the E–L method, the particles are tracked in a Lagrangian frame, while the particle–particle interactions are modeled using solid stress calculated from the Eulerian grid [58,59]. Since the interparticle interactions are modeled instead of being directly resolved, the computational overhead of the solution is reduced compared to the unresolved E–L approach, while the model is superior over the E–E model (e.g., straightforward handling of poly-dispersed particles).

**Table 1.** Summary of the available approaches for computational fluid dynamics (CFD) simulation of gas–solid flows with the inclusion of closure models to be considered in each approach. The table is inspired from [60].

#### **4. CFD Simulation Studies of Gas–Solid Cyclones at Ambient Temperature**

When solid particles are added to the gas flow in a cyclone separator, the mean and turbulent properties of the flow field may change. According to the classification of gas–solid interactions made by Elghobashi [61] (see Figure 4), for particle volume fractions lower than 10−<sup>6</sup> (equivalent to the interparticle spacing of around 80*dp* [20]), particles have a negligible influence on the turbulent flow field, and the gas–solid coupling is termed as one–way coupling. When the particle volume fraction reaches values higher than 10−<sup>6</sup> but below 10−<sup>3</sup> (equivalent to the interparticle spacing of around 8*dp* [20]), the particles either attenuate or intensify the turbulence in the gas–solid flow, i.e., two-way coupled flow field, depending on the relative kinetic response time of particles to the characteristic flow time-scale [61]. In the context of gas–solid cyclone flows, attenuation is usually the case [20]. If the solids volume fraction is increased further, namely to the four-way coupling region, the interaction between particles becomes important too.

**Figure 4.** Map of gas–solid interaction regimes of particle-laden turbulent flows. *τp* and *τe* are the particle kinetic response time and time-scale of large eddies in a turbulent flow, respectively. Reproduced with permission from S. Elghobashi, On predicting particle-laden turbulent flows; published by Springer Nature, 1994 [61].

There are a vast number of CFD simulation studies of cyclone separators at ambient temperature, which are usually carried out for pilot-scale cyclone designs. In the current study and according to the above classification, these studies are separated into three groups of (I) single-phase gas flow; (II) one-way coupled gas–solid flow; and (III) two- and four-way coupled gas–solid flow in cyclones.

#### *4.1. Group I: Single-Phase Flow CFD Simulations*

As briefly mentioned earlier, the turbulent flow inside a particle-free conical cyclone is anisotropic and inherently unsteady due to the significant swirl and radial shear intensity as well as the adverse pressure gradient, which may induce recirculation regions. Accordingly, the turbulence model utilized in the CFD simulation of cyclones should be capable of correctly resolving or modeling the high curvature and intensive swirl flow and, at the same time, would be appropriate for modeling the adverse pressure gradient and recirculating flows [20].

As a result of being based on an isotropic turbulent flow assumption, the commonly used classical first-order turbulence models of the steady-state and unsteady Reynolds-averaged Navier–Stokes equations (URANS), i.e., zero-, one-, and two-equation models such as the k- model, are usually unable to accurately reproduce the flow field inside cyclones [62–65]. However, there are some attempts that have been carried out to improve the accuracy of the first-order turbulence models with the inclusion of anisotropy effects in their modified versions, e.g., the k- model based on renormalization group (RNG) theory [63,65], the k-curvature correction (cc) model [66], and the model proposed by Meier and

Mori [67]. These models provide improved accuracy in the prediction of single-phase flow in cyclones compared to the standard two-equation models. Alternatively, the second-order URANS closure models, e.g., the Reynolds stress transport model (RSTM) [66,68–70] and scale-resolving turbulence models, e.g., large eddy simulation (LES) [71,72], are able to capture the anisotropy effects in highly swirling flows and have been proven to be accurate in predicting the cyclone gas flow field. Among the two widely used variations of the differential RSTM in modeling the pressure-strain term, i.e., the Launder, Reece, and Rodi (LRR) model [73] and the Speziale, Sarkar, and Gatski (SSG) model [74], the latter provides superior accuracy [66]. With the recent advancements in computer capabilities, the use of LES for the prediction of the flow inside cyclones has become possible, but with a higher computational overhead cost. Using LES, it is possible to study cyclone flow physics in extensive detail and the accuracy of results is superior compared to the RSTM [75–77]. Nonetheless, the RSTM can provide reasonable predictions of the gas flow in cyclones using a relatively coarser computational grid compared to the ones used in LES studies, and at a much lower computational overhead than LES [76].

Apart from the turbulence model, there are important considerations related to the numerical aspects of CFD simulations of gas flow in cyclones. Due to the presence of strong velocity and pressure gradients, the use of high-order (i.e., second order and above) methods in the discretization of advection terms as well as the alignment of the computational cells with the flow direction (i.e., structured grids) are recommended in order to avoid unwanted numerical diffusion (for example, see [78]).

There are a large number of existing CFD studies focusing on single-phase flow characteristics in cyclones that are not fully listed here as it is not the main subject of this review. A summary of some of the previous single-phase flow studies, up to 2007, can be found in [20].

#### *4.2. Group II: One-Way Coupled Gas–Solid Flow Simulations*

In the second category of existing CFD studies of cyclones operating at ambient temperature, category (II), the one-way coupling method is employed for very dilute loading of particles (see Figure 4), and the E–L approach is commonly utilized. The particles are tracked either after the simulation has reached a steady-state solution, as a post-processing step [69,79–87], or in a transient way along the unsteady simulation [75,79,88–93]. After the completion of the particle-tracking procedure, the particle separation efficiencies are predicted and usually compared with the available experimental data. In general, it can be stated that there is a tendency for the particle cut-size diameter to be underpredicted in LES calculations [88,89,93] and overpredicted in RSTM calculations [20,93]. Furthermore, the separation efficiency is predicted more accurately when LES is used [93], due to a better prediction of the mean and fluctuating velocity fields. One must keep in mind that the accuracy of the predicted separation efficiency, however, depends on the validity of the basic assumptions for the simulations, e.g., the accumulation of particles near the cyclone walls may cast doubt on the validity of the one-way coupling assumption.

An important consideration in one-way coupled E–L simulations is the influence of turbulent dispersion, i.e., the movement of particles due to the presence of velocity fluctuations, especially for smaller-sized particles. If the resolved scales have comparable time-scalse as compared to the kinematic response time of particles, e.g., LES and DNS, the turbulent dispersion is resolved directly and no additional turbulent dispersion model is required [79,94], while the effect of sub-grid scale (SGS) structures in LES on particle movement can usually be neglected. For URANS simulations, however, it is common to apply a stochastic turbulent dispersion model [95]. Accurate prediction of the gas velocity fluctuations, i.e., rms velocities, is important for properly resolving particle dispersion due to turbulence. In some of the existing studies, it is reported that even though the mean velocities are properly captured by the RSTM, the rms fluctuations of the velocity field are underpredicted compared to the measurements [76,96], due to a failure in accurately predicting PVC behavior that gives rise to the rms velocities, especially in central regions of a cyclone [76]. Shukla et al. [96] reported

an overprediction of separation efficiency for smaller-sized particles due to an underprediction of velocity fluctuations, as shown in Figure 5. On the other hand, LES is able to capture both mean and rms velocities properly, leading to accurate predictions of grade efficiency compared to the measurements [96].

**Figure 5.** Comparison of predicted mean velocity profiles (**a**,**b**), rms velocities (**c**,**d**), and grade efficiencies (**e**) using the Reynolds stress transport model (RSTM) and large eddy simulation (LES) for a cyclone with a body diameter of 0.29 m and operating at ambient temperature and pressure [96]. The experimental data on separation efficiency are from [97]. Reproduced with permission from S. Shukla et al., The effect of modeling of velocity fluctuations on prediction of collection efficiency of cyclone separators; published by Elsevier, 2013.

#### *4.3. Group III: Two- and Four-Way Coupled Gas–Solid Flow Simulations*

Both the E–E and the E–L methods are employed for CFD simulation of two- and four-way coupled gas–solid flows. Generally speaking, and as discussed earlier, the volume fraction of particles in a cyclone separator is locally inhomogeneous, as the design of such equipment demands. Consequently, in the existing studies of cyclone separators, the E–L approach is preferred over the E–E approach.

#### 4.3.1. E–L and Hybrid Model Simulations

As for the third category of CFD studies, category (III), compared to the previous categories, there are fewer E–L studies with two-way coupling [34,77,98–105] and only a limited number of studies with four-way coupling [35,106–109]. In almost all of the studies (except [99]), due to the lack of available experimental data for the cyclone flow loaded with particles, no comparison of the predicted modification of velocity profiles (compared to the single-phase flow) with the measurements is provided. Instead, the predicted overall/grade efficiencies [34,77,101,102] and/or pressure drops [34,35,98,101,102,106,107,109] are compared with the measurements for certain ranges of particle mass loading. A summary of the existing four-way coupled CFD studies of gas–solid cyclones operating at ambient temperature and with inlet mass loading ratios higher than 0.1 kg*s*/kg*g* is presented in Table 2.



**Table**

16

#### *Processes* **2019**, *8*, 228

 order.

#### Effect of Solids Loading

As the solids loading ratio increases, the particle–particle interactions become important and neglecting these interactions in the cyclone CFD simulation may lead to poor predictions of cyclone performance. For example, Wasilewski and Duda [115] have reported an underprediction of both pressure drop and separation efficiency for two different pilot-scale cyclone geometries, though the overall trend of changes was well-predicted. The authors have claimed that this discrepancy is due to the simplification of the CFD model, neglecting particle–particle collisions and agglomeration. Both CFD–DEM [107–109] and hybrid E–E and E–L approaches [35,106] have been used to consider the effect of particle–particle interactions in existing cyclone separator simulation studies. In the CFD–DEM study of Zhou et al. [107], the effect of large (2000–2800 μm) and "ultra light" particles with a solid to gas density ratio of 33 is investigated for different loadings of particles up to 0.46 kg*s*/kg*g*. They have analyzed the significance of gas–solid and solid–solid forces in different regions of cyclones and concluded that the modification of the velocity field due to the presence of particles is hardly noticeable. On the other hand, based on CFD–DEM calculations, Chu et al. [109] reported a significant modification of the velocity field due to the presence of particles of nearly the same size but much heavier compared to the ones studied in [107]. According to this study, the particle–particle forces, which are much more dominant than the particle–gas forces, increase in magnitude with the loading of particles, with the most active region being inside the particle strands along the cyclone wall.

In most of the existing four-way coupled CFD simulations of cyclones using the hybrid methods [35,106], the focus is on the effect of solids loading on the performance of cyclones. Hwang et al. [35] carried out DDPM–KTGF simulations of a pilot-scale cyclone laden with 2000 μm diameter particles and reported a minimum in the predicted pressure drop versus solids loading. The comparison with available measurements, though, is carried out only for smaller solids loadings, i.e., loadings smaller than 2.5 kg*s*/kg*g*, while the minimum pressure drop takes place at solids loadings of around 6 kg*s*/kg*g*. Conversely, in the hybrid E–E and E–L study of Kozolub et al. [106], an increasing trend in the predicted pressure drop of a cyclone laden with 40–80 μm diameter particles is observed for all of the studied mass loadings up to 2.6 kg*s*/kg*g*. This trend is in line with the available measurements, but the values are somewhat underpredicted.

One of the important features that is captured by two- and four-way coupled CFD simulations of cyclones is the modulation of velocity field and turbulence statistics as compared to the particle-free flow. The swirl velocity is predicted to reduce substantially with the addition of particles to the flow and decrease further as the loading of particles increases [35,102,106,109,116]. The predicted reduction of tangential velocity magnitude as well as turbulent fluctuations by the addition of 0.05 and 0.1 kg*s*/kg*<sup>g</sup>* of particles to a pilot-scale cyclone is shown in Figure 6 based on the LES study of Derksen et al. [116]. The reduction of swirl is reported to be stronger in the loss-free vortex part of the swirl due to the higher concentration of particles in this region. On the other hand, the tangential velocity at regions close to the walls is still high (see Figure 6) due to the presence of particles and partial transfer of tangential momentum to the gas phase. The reduction of swirl can, in turn, reduce the effectiveness of centrifugal force on particle separation in the cyclone, making the particles be less concentrated near the walls. However, at the same time, the turbulence is attenuated by the presence of particles all over the cyclone body [116], leading to the reduction of turbulent dispersion; also, the amount of gas bypassed to the vortex finder at the entrance of the cyclone is reduced. The overall separation efficiency with respect to the solids loading is a balance between the significance of the above-mentioned effects.

The modification of cyclone efficiency with particle mass loading is not reported in the existing four-way coupled CFD studies as in most of the studies [35,107,109], the addition of large particles to the cyclone, e.g., 2000 μm in diameter, is investigated, leading to an ideal separation efficiency.

**Figure 6.** The predicted time-averaged values of axial and tangential velocities as well as the resolved turbulent kinetic energy of the gas for a pilot-scale cyclone with a body diameter of 0.29 m and operating with different mass loadings of particles using two-way coupled LES at 0.75D (top) and 2D (bottom) below the cyclone roof [116]. Reproduced with permission from J.J. Derksen, Simulation of mass-loading effects in gas-solid cyclone separators; published by Elsevier, 2006.

#### Agglomeration of Particles

As mentioned earlier, an important mechanism in solids separation in cyclones is particle agglomeration [19,30]. To model the particle agglomeration in a Lagrangian frame, it is common to assume either volume-equivalent (i.e., the new diameter has the same volume as the agglomerated particles) [117,118] or inertia-equivalent (i.e., the new diameter is calculated to maintain the inertia) [108,119,120] agglomeration. The agglomeration of particles in cyclones is studied in the CFD–DEM simulations of Sgrott and Sommerfeld [108] with a stochastic approach [121] to describe the particle–particle collisions. It is assumed that the only driving factor in interparticle adhesion is the van der Waals force. It has been shown that using the inertia-equivalent approach, the predicted agglomeration rate and the separation efficiency are slightly higher and lower, respectively, compared to the volume-equivalent method (see Figure 7). Furthermore, it is shown that for an inlet mass loading ratio of 0.1 kg*s*/kg*g*, the difference between the predicted velocity profiles considering twoand four-way coupling methods is small. However, when comparing the predicted grade efficiencies of two- and four-way coupled simulations, as shown in Figure 7, for very small particles, i.e., less than 2 μm in diameter, neglecting the agglomeration may lead to an underprediction of separation efficiency. This is in agreement with the agglomeration concept of Mothes and Löffler [30] regarding particle separation in cyclones, discussed in Section 2.2. The agglomeration effect is also investigated through one-way coupled CFD simulations of gas–solid cyclones [122].

#### 4.3.2. E–E Simulations

E–E studies of gas–solid cyclone separators are very scarce, with the early investigations being allotted to the prediction of the change of cyclone hydrodynamics with the addition of particles [67,123]. In most of the existing E–E studies of gas–solid cyclones, the comparison with experimental data is carried out for cyclones operating with very low particle mass loadings [124–126], i.e., below 0.01 kg*s*/kg*g*, while the true advantage of using the E–E approach is for the high loading of particles. Costa et al. [126] implemented an E–E simulation of gas–solid flow in a cyclone with four-way coupling. At a solid mass loading of 0.009 kg*s*/kg*g*, the predicted grade efficiencies were improved by increasing the number of solid phases (corresponding to each size bin), but still overpredicted compared to the experimental data. Variations of this model are further utilized for optimization

studies of gas–solid cyclones operating with either a very low mass loading of particles [127,128] or moderate loadings [129], i.e., 0.1 kg*s*/kg*g*, and higher.

For a particle mass loading of around 1 kg*s*/kg*g*, mixture models of cyclones are available [130,131], based on modeling the velocity slip between the phases (simplest approach mentioned in Table 1). The pressure drop estimated by these models is underpredicted [130,131], especially at high inlet solids loadings, i.e., 0.4 kg*s*/kg*g*, which, according to the authors, is due to the limitation of the granular and mixture models for densely loaded cyclones. Furthermore, the predicted separation efficiencies, compared to the measurements, are not conclusive [130]. In both studies, however, a reduction of swirl by increasing the mass loading is reported.

**Figure 7.** Grade efficiency predicted by the CFD simulation of Sgrott and Sommerfeld [108] for a pilot-scale cyclone loaded with particles with a diameter of 0.5–60 microns and a mass loading of 0.1 kg*s*/kg*g*. 1 W-C, 2 W-C, and 4 W-C refer to one-way coupling, two-way coupling, and four-way coupling methods (using CFD–DEM without agglomeration), respectively. Sphere and history models are volume-equivalent and inertia-equivalent approaches for agglomeration, respectively. Reproduced with permission from O.L. Sgrott and M. Sommerfeld, Influence of inter-particle collisions and agglomeration on cyclone performance and collection efficiency; published by John Wiley and Sons, 2018.

#### *4.4. Summary*

In this section, the existing CFD simulation studies of reversed-flow cyclone separators operating at ambient temperature are discussed and categorized based on the gas–solid coupling regime of the operating cyclone. For the single-phase gas flow in cyclones, the existing CFD simulation studies are able to accurately capture the velocity and pressure fields as compared to the available experimental data. According to these studies, the RSTM is favorable over other turbulence models with respect to accuracy and computational overhead. However, some modified versions of first-order turbulence models, e.g., k- curvature correction (cc) model [66], can still be used in the CFD simulations, especially when the solids loading is high and the swirl is less intense.

For dilute loading of particles and using the one-way coupling method, it is possible to properly capture the separation efficiency of particles in the cyclone as long as the one-way coupling condition is valid for the particle-laden flow. LES has proven to be more accurate in the prediction of grade efficiency, especially for small-sized particles, compared to RSTM, due to an improved prediction of velocity fluctuations and turbulent dispersion.

Among the available gas–solid simulation approaches listed in Table 1, the unresolved E–L and hybrid models are used more frequently, for dilute/medium and dense loading of particles in cyclones, respectively. For densely loaded cyclones, a few CFD studies focus on the effect of solids loading on performance, applying either DDPM–KTGF or CFD–DEM approaches, while only a handful of them investigate the loading of small-sized particles, i.e., 0.1–200 μm in diameter. In the existing CFD studies

of densely loaded cyclones, comparison of the modified flow field with measurements is carried out very rarely, mainly due to challenges in the measurement methods for such systems. Furthermore, comparison/validation of the predicted pressure drops against measurements is usually carried out, rather than the separation efficiency.

#### **5. CFD Simulation Studies of Gas–Solid Cyclones at Elevated Temperatures**

In general, the available studies on the CFD simulation of cyclones operating at elevated temperatures are more scarce than the ones for cyclones operating at ambient temperature. Similar to the previous section, the CFD simulations of cyclone separators operating at elevated temperatures can be divided into three groups: single-phase flow, one-way coupled flow, and two/four-way coupled flow simulations. In this section, the three groups of studies are discussed with a focus on the process or physics that is investigated through the CFD simulation.

#### *5.1. Group I: Single-Phase Flow CFD Simulations*

Once the operating temperature increases, the physics of gas flow inside the cyclone become more complicated, mainly due to the heat transfer and subsequently variable gas density and viscosity in local regions, if the insulation is not perfect. In turn, since an extra transport equation for energy needs to be solved, the CFD calculation of gas flow becomes computationally more complex and demanding. However, in some of the existing CFD studies of cyclones operating at elevated temperatures, in order to avoid this complexity, the simulations are carried out at a fixed gas temperature, while the gas density and viscosity values are modified according to the operating temperature. Using this method, the predicted pressure drop is reported to be of acceptable accuracy compared to the measurements [131,132].

Similar to the CFD simulations of the cyclones at ambient temperature, accurate modeling of the turbulence is a must to achieve acceptable accuracy for high-temperature single-phase flow simulations of cyclones. However, as mentioned earlier, the increased viscosity as well as decreased density of the carrier gas at elevated temperatures leads to a reduction of swirl intensity and turbulence, indicating the possibility of using a wider range of available turbulence models while keeping the prediction accuracy reasonable. Gimbun et al. [133] have compared the predicted pressure drop of a pilot-scale cyclone with the experimental data of [39] using the RNG *k* −  and RSTM. At lower temperatures, i.e., 300–800 *K*, RSTM is reported to provide a more accurate prediction of the pressure drop compared to the RNG *k* − *-*, while in the temperature range of 800–1200 *K*, the predicted pressure drops using RNG *k* −  and RSTM are nearly the same and in agreement with the experimental data [133]. In the existing CFD studies of cyclones operating at elevated temperatures, both the *k* − *-* [134] and RSTM [132] are used as turbulence models.

The tendency for weakened swirl and the modification of the internal velocity field, as an outcome of an elevated operating temperature, is reported in some of the existing CFD studies. According to the CFD simulation results of Shi et al. [132], an overall increase in the axial velocity profiles, weakened reverse flow at the center of the inner vortex, and weakened tangential velocity peaks (weakened swirl), are reported as the operating temperature of the studied pilot-scale cyclone increases from 20 ◦C to 800 ◦C. They reported that as the cyclone swirl intensity is weakened at elevated temperatures, the pressure drop is also reduced, and this reduction is dominantly affected by the gas density rather than its viscosity when the volumetric flow rate to the cyclone is kept constant [132]. The reduction of pressure drop with the gas temperature is reported in many of the existing CFD studies of cyclones operating at elevated temperatures [41,132,134–136]. An approach for the comprehensive analysis of local energy loss in cyclones is through evaluating entropy generation in CFD simulation, which has been carried out for cyclones operating at ambient [137,138] and elevated [135] temperatures. Duan et al. [135] have numerically investigated entropy generation in a cyclone system operating in a temperature range of 297–1123 *K*. According to this study, turbulence dissipation and wall friction are the main contributors to the energy loss. The ratio of contribution of both parameters is

nearly unaffected by the gas temperature, whereas their magnitude decreases at higher operating temperatures due to the reduction in gas density. The maximum energy loss takes place in areas close to the vortex finder and the entrance of the dust bin.

#### *5.2. Group II: One-Way Coupled Gas–Solid Flow Simulations*

In the existing CFD studies of cyclones operating at elevated temperatures, the predicted separation efficiency is negatively affected by the temperature. Gimbun [41] has performed one-way coupled simulations of a pilot-scale cyclone using RSTM and a stochastic tracking of particles in a temperature range of 20–900 ◦C. For a fixed inlet gas velocity, the predicted grade separation efficiencies and particle cut-sizes were in proper agreement with the experimental data, and overall, the separation efficiency worsened with the temperature. The reduction in separation efficiency was stated to be a direct consequence of the weakened centrifugal effect. In the same line of research, Gimbun et al. [139] presented an accurate prediction of the particle cut-size for different cyclone types compared to experiments, at ambient and elevated temperatures. The reduction of separation efficiency is also reported in some of other studies (for example, see [134,136]).

One-way coupled heat transfer to solids particles fed in to a pilot-scale cyclone was studied by Mothilal and Pitchandi [140,141] using E–L simulation and RNG *k* −  turbulence models and at an inlet air temperature of 200 ◦C. The predicted pressure drop was compared with the available experimental data for the particle-free case, while for the gas–solid cases, no validation was presented. By increasing the gas velocity at the cyclone inlet, the centrifugal effect on particles was improved and the predicted particle hold-up in the system was increased, especially for larger mass loadings of particles. This leads to a greater available surface area (for heat transfer) as well as a higher residence time of particles in the cyclone system.

#### *5.3. Group III: Two- and Four-Way Coupled Gas–Solid Flow Simulations*

In this group of studies, both E–L and E–E and hybrid approaches are applied for the CFD simulation of cyclones operating at elevated temperatures; these are discussed here with a focus on the studies that have investigated the performance of industrial-scale cyclones.

#### 5.3.1. Cyclone Heat Exchangers

An important industrial application of cyclones is the solids preheating process used in, e.g., cement industry. The performance evaluation for full-scale cyclone heat exchangers is carried out in some of the existing studies and summarized in Table 3. Cristea and Conti [142,143] have applied a DDPM–KTGF approach to simulate the gas–solid flow in the preheater system (excluding the calciner) of a cement factory with an approximate height of 58 m. The gas–solid heat transfer as well as the calcination reaction are included in the CFD model [143]. As compared to the industrial measurements for the first-stage twin cyclones, the predicted pressure drop and separation efficiencies are in good agreement with the measurements, while the exit gas temperature is overpredicted. This discrepancy has been attributed to the complications that arise with the modeling of the high mass loading of particles in the upper-stage calciner [143]. Mikulˇci´c et al. [144] have conducted a two-way coupled CFD simulation study of an industrial-scale cyclone belonging to the preheater system of a dry kiln process. The simulations aimed to assess the gas–solid heat transfer as well as the progress of the calcination process in the studied cyclone. The predicted pressure drop was in fair agreement with the measurements, while the outlet gas temperature was slightly overpredicted, and the authors did not provide any reasoning for this mismatch. In this CFD study, low gas temperature regions in the upper part of the cyclone and close to the walls were observed, due to the endothermic calcination reaction. Wasilewski [145] conducted a two-way coupled CFD simulation of a full-scale cyclone preheater and reported a reduction of the separation efficiency with increased temperatures and decreased particle loading. However, as compared to the measurements, an underprediction of the separation efficiencies is reported, which is claimed to be due to the limitation of the CFD model in taking into account the

particle–particle interactions, i.e., collision and agglomeration. Presented in Figure 8 is their predicted trajectories and heat-transfer rate to laden particles of different diameters. The smaller-sized particles have a higher tendency to travel through the central regions of the cyclone body and escape from the vortex finder. In addition, they have a higher rate of heat exchange with the gas, especially at regions close to the inlet, compared to the larger-sized particles. In contrast, larger-sized particles have a higher tendency to travel in spiral clouds along the cyclone walls and toward the particle exit chamber.

**Figure 8.** Predicted trajectories of particles of different diameters in an industrial-scale cyclone with a body diameter of 3.45 m. The trajectories are colored with particle temperature. The inlet gas and particle temperatures are 634 *K* and 611 *K*, respectively, and inlet solids mass loading is around 0.8 kg*s*/kg*g* [145]. Reproduced with permission from M. Wasilewski, Analysis of the effects of temperature and the share of solid and gas phases on the process of separation in a cyclone suspension preheater; published by Elsevier, 2016.

**Table 3.** Selected CFD studies of industrial-scale cyclones operating at elevated temperatures, in chronological order.


There are also a few E–E simulation studies of cyclones at elevated temperatures. In the study of Vegini et al. [124], 2D axisymmetric simulations of full-scale cyclone separators connected in series and operating at elevated temperature were carried out. The selected turbulence model is a combination of the standard *k* −  and Prandtl's longitudinal mixing model, which is claimed to be accurate enough for capturing the anisotropic effect of Reynolds stress in swirling turbulent flows. The numerical model was validated against the available grade efficiencies and pressure drops at ambient temperature and a dust load of around 0.004 kg*s*/kg*g*. In elevated temperature conditions, the predicted pressure drop values were generally overpredicted compared to the available long-term full-scale measurements. This overprediction was claimed to be due to the presence of false air in the system, and as a consequence, a reduction of swirl. The presence of false air is not considered in the simulations.

#### 5.3.2. CFD Simulation of Cyclones as a Part of a Bigger System

There are a number studies that have carried out CFD simulations of elevated-temperature cyclones as a part of a bigger process system, e.g., a circulating fluidized bed (CFB) boiler. An important phenomenon that is explored in this type of study is the fluctuation and maldistribution of gas and solid flows at the inlet of parallel cyclones, which can, in turn, affect cyclone performance. The subject is investigated both experimentally [147–150] and numerically [151,152]. Zhang et al. [151] have reported temporal synchronized fluctuations of solid flux at the inlet of two parallel cyclones based on four-way coupled E–E simulations of a 150 *MW* CFB boiler operating at 917 ◦C. Jiang et al. [152] have carried out hybrid E–L and E–L (MP–PIC) CFD simulations of a scaled down circulating fluidized bed comprised of six parallel cyclones. They reported a non-uniformity in the average solid concentration downstream of the cyclone diplegs, based on both simulations and experimental data. However, some of the predicted solid concentrations were somewhat different from those of the measurements, with the highest relative error being around 28%.

#### *5.4. Summary*

The existing CFD studies of cyclones operating at elevated temperatures were discussed in this section, with a focus on the ones investigating cyclone performance parameters, i.e., pressure drop, separation efficiency, and gas–solid heat transfer. The gas-flow simulation studies reveal that CFD simulations are capable of capturing the reduced swirl intensity and reduced pressure drop as a result of elevated temperature. Furthermore, in some of the studies, the first-order turbulence models are satisfactorily used to calculate the gas flow field. The reduction in separation efficiency of dilutely loaded cyclones, i.e., one-way coupled studies, due to the increase in temperature, is captured. A number of CFD studies of industrial-scale cyclones are also described. However, it is common for the full-scale measurements of these systems not to exist or to be limited, e.g., to only the pressure drop of the cyclone for a specific operating condition, indicating the lack of complete validation of the accuracy of current CFD studies.

#### **6. Outlook**

In this review, the capability of computational fluid dynamics in predicting the gas–solid flow field and performance of cyclone separators, operating at dilute–high loading of particles and ambient and elevated temperatures, is discussed. Summaries of ambient- and elevated-temperature CFD studies are provided in Sections 4.4 and 5.4. In general, it can be stated that the CFD simulation can be used as a powerful tool to successfully predict the flow field and performance of cyclone separators, i.e., pressure drop and separation efficiency, with the E–L and hybrid methods being the most frequently used approaches. However, as the mass loading and the operating temperature increase, the validation of the internal flow field, which is important for a better understanding of the separation process, becomes challenging due to the lack of experimental data. For these cyclones, the validation is so far limited to the comparison of the pressure drop and seldom the measured separation efficiency.

The drag models that are usually applied in two- and four-way coupled simulations are well-known homogeneous models, e.g., Wen–Yu [113] and Gidaspow [50], while heterogeneous models, e.g., the energy-minimization-multi-scale (EMMS) approach [153], are rarely used (for example, see [111]). In the authors' opinion, heterogeneous drag models can be more suitable for CFD simulation of heavily loaded cyclones (e.g., inlet mass loading on the order of 1–100 kg*s*/kg*g*), as it is likely to have cluster formation in the internal regions, apart from regions close to the cyclone walls, especially for the elevated operating temperatures with weakened centrifugal effect.

**Author Contributions:** Conceptualization, M.N. and H.W.; investigation, M.N.; writing–original draft preparation, M.N.; writing–review and editing, M.N.; supervision, H.W., B.L., Y.T., W.W., K.D.-J.; project administration, H.W.; funding acquisition, H.W., K.D.-J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project is supported by the Sino-Danish Centre for Education and Research; and the ProBu - Process Technology for Sustainable Building Materials Production - project (grant number: 8055-00014B) funded by Innovation Fund Denmark, FLSmidth A/S, Rockwool International A/S and Technical University of Denmark.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Comparison of Surface Tension Models for the Volume of Fluid Method**

#### **Kurian J. Vachaparambil \* and Kristian Etienne Einarsrud**

Department of Materials Science and Engineering, Norwegian University of Science and Technology (NTNU), 7491 Trondheim, Norway

**\*** Correspondence: kurian.j.vachaparambil@ntnu.no; Tel.: +47-45517425

Received: 24 July 2019; Accepted: 13 August 2019; Published: 15 August 2019

**Abstract:** With the increasing use of Computational Fluid Dynamics to investigate multiphase flow scenarios, modelling surface tension effects has been a topic of active research. A well known associated problem is the generation of spurious velocities (or currents), arising due to inaccuracies in calculations of the surface tension force. These spurious currents cause nonphysical flows which can adversely affect the predictive capability of these simulations. In this paper, we implement the Continuum Surface Force (CSF), Smoothed CSF and Sharp Surface Force (SSF) models in OpenFOAM. The models were validated for various multiphase flow scenarios for Capillary numbers of 10<sup>−</sup>3–10. All the surface tension models provide reasonable agreement with benchmarking data for rising bubble simulations. Both CSF and SSF models successfully predicted the capillary rise between two parallel plates, but Smoothed CSF could not provide reliable results. The evolution of spurious current were studied for millimetre-sized stationary bubbles. The results shows that SSF and CSF models generate the least and most spurious currents, respectively. We also show that maximum time step, mesh resolution and the under-relaxation factor used in the simulations affect the magnitude of spurious currents.

**Keywords:** surface tension modelling; VOF; rising bubbles; capillary rise

#### **1. Introduction**

For a comprehensive understanding of flow physics in multiphase systems, which is ubiquitous in both nature and technological processes, consideration of surface tension is important. For instance, the break down of a fluid jet into droplets is used to form droplets in inkjets [1] and lab-on-chip devices [2] while the thinning and breakdown dynamics of non-Newtonian fluid filaments is critical in its application in jetting [3,4]. Flow scenarios such as underground water flows [5], oil recovery [6] and paper-based microfluidics [7] are examples of flow through porous media where dominance of surface tension may produce a capillary rise. The detachment diameter of the bubble [8,9] and shape of rising bubble [10] during bubble evolution in champagne, boiling and electrochemical gas evolution is dependent on surface tension, as is the droplet size produced during atomisation of fuels [11], spraying [12,13] and growth of a bubble in confined geometries [14]. The effect of surface tension is also important in events such as nucleation of bubbles [15,16] and droplets [17].

Due to the importance and complex nature of multiphase flows, numerical simulations, especially computational fluid dynamics (CFD), are commonly used to study and understand these processes. The CFD strategies used to model multiphase flows can broadly be divided into four categories: Euler–Euler (EE), Euler–Lagrange (EL), interface tracking and capturing methods. The EE approach assumes that phases are interpenetrating, which is efficient when modelling large-scale industrial processes [18,19], while EL tracks the dispersed phases individually, which can be computationally expensive [20,21]. As both EE and EL approaches do not resolve the complete interactions between

the phases, they require so called "closure laws" (see [18–21]). Interface tracking methods, such as the moving mesh method, use a separate boundary-fitted moving mesh for each phase [22]. Although interface tracking methods are quite accurate, they are typically used to model bubble or droplets with mild-moderate deformations [22,23] but to handle complex interface deformations these methods require a global or local re-meshing [24]. Interface capturing methods use a fixed grid with functions to capture the interface such as the Volume Of Fluid method (VOF) [25], level-set [26] and diffuse interface methods [27]. Other methods available in the literature employ a hybrid interface tracking-capturing approach, such as the immersed boundary [28] and front tracking method [29]. Due to its ability to conserve mass (both level-set [30] and phase-field [31] models have difficulties in conserving mass), robustness and ability to produce reasonably sharp interfaces VOF is very popular in multiphase simulations [32–57] and implemented in both commercial (ANSYS Fluent® and Flow-3D®) and open source (OpenFOAM®) CFD packages.

Due to the popularity of open source CFD packages, this paper predominantly delves into the VOF formulation and reported development in interFoam, which is the VOF-based solver available in OpenFOAM®. In the VOF method, a scalar function representing the volume fraction of phases in the computational cells is advected. The advection of the volume fraction equation is done using specific discretisation schemes, such as the interface compression method [58], to prevent the excessive smearing of the interface thickness. Apart from interface compression method, recent work has explored reconstruction of interface using techniques such as the isoAdvector method [59,60] and piece wise linear interface calculation (PLIC) algorithm [61]. Although the VOF approach in theory produces a sharp interface, the "real" VOF, which is implemented in solvers such as interFoam, produces a non-sharp interface, which stretches over a few computational cells. This non-sharp nature of the volume fraction leads to errors in the calculated curvature which generates spurious currents that is quantified in the work by Harvie et al. [62], appearing as vortices around the interface (see [63,64]). The presence of these spurious currents introduces non-physical velocities near the interface, which can increase the interfacial mass transfer while modelling condensation [32] and evaporation [57] scenarios and adversely effects the accuracy of simulations. In the literature, spurious currents in VOF methods can be reduced by the following approaches:




To analyse the force balance (described in [65]), Deshpande et al. [63] evaluated interFoam and showed that there is no imbalance in the surface tension and pressure forces due to inconsistent discretisation. However, the iterative process, which is used to solve pressure equation, introduces an imbalance which is related to the user defined tolerance level of the solution [63]. An overview of literature that provides an improved estimate of the interface curvature and surface tension modelling approaches is provided in Table 1. The improved representation of the interface (which aids in accurate calculation of the interface curvature) is provided bycoupled level-set VOF (CLSVOF) method, height functions and interface reconstruction algorithms (like piecewise linear interface calculation (PLIC), parabolic reconstruction of surface tension (PROST) and efficient least squares volume-of-fluid interface reconstruction (ELIVRA) algorithms), whereas the other methods discussed in Table 1 provide alternative approaches to model surface tension. To ensure that spurious currents do not grow over time, a stability condition, proposed by Brackbill et al. [66], for explicit treatment of surface tension is

$$
\Delta t < \sqrt{\frac{\rho\_{\text{avg}} (\Delta x)^3}{2\pi \sigma}},
\tag{1}
$$

where Δ*x*, *σ* and *ρavg* are grid spacing, surface tension and average of density of both phases, respectively. As proposed by Galusinski and Vigneaux [79], a comprehensive constraint on the time step must consider the effect of both density and viscosity which can be expressed as

$$
\Delta t \le \frac{1}{2} \left( \mathbb{C}\_2 \tau\_\mu + \sqrt{(\mathbb{C}\_2 \tau\_\mu)^2 + 4\mathbb{C}\_1 \tau\_\rho^2} \right),
\tag{2}
$$

where *τμ* and *τρ* are time scales which are defined as *<sup>μ</sup>avg*Δ*x*/*<sup>σ</sup>* and *ρavg*(Δ*x*)3/*σ*, where *μavg* is the average dynamic viscosity between the phases, respectively. An evaluation of interFoam, by Deshpande et al. [63], proposed that time step should satisfy

$$
\Delta t \le \max \left( C\_2 \tau\_{\mu \nu} 10 C\_1 \tau\_{\rho} \right) \tag{3}
$$

along with the time step constraint discussed in Equation (2). Deshpande et al. [63] also calculated the values of *C*<sup>1</sup> and *C*<sup>2</sup> for interFoam to be equal to 0.01 and 10, respectively. Further details of the numerical methods used to model surface tension is available in the recent review work by Popinet [80].

In the literature, comparison between surface tension models is often done for a specific of flow phenomenon and at times a static scenario is used to quantify the spurious currents. Some examples of flow phenomena used to compare surface tension models are rising bubbles whose diameters are in the order of few millimetres [33–35], translating and rotating bubbles [64], oscillating droplets or bubbles [34], stagnant bubbles or droplets [34,35,39,64], Rayleigh–Taylor instability [37,38], Taylor bubbles [64], falling films [41], droplet splashing [38,39], capillary rise [42] and bubble evolution [37,40]. These typically compare the CSF model with height functions [33,34,64], PROST [37], PLIC [42], CLSVOF and its variants [37–40,64], FSF and SSF [42], and CSS [35,41] models. Although the flow scenarios that are used to compare surface tension models are diverse, they can be broadly categorised based on the dominance of surface tension in the flow using the Capillary number (*Ca*), which is defined as the ratio of viscous to surface tension forces in the system. Flow phenomena such as capillary rise and stationary bubbles are examples of low values of *Ca* whereas flows with larger values of *Ca* include rising bubbles and falling films.

During processes such as gas evolution during electrochemical reactions and boiling, nucleated bubbles grow by mass transfer across the interface [15,16] or coalescence [8], but once the bubble detaches it may deform as it rises up and/or interacts with other bubbles [53]. Other complex processes, such as splashing, involve droplet spreading on a surface which is accompanied by formation of secondary smaller droplets at the rim [81]. To reliably model these processes, surface tension models must be able to accurately handle flow scenarios with both small and large capillary numbers. In literature, the work by Hoang et al. [56] implemented the Smoothed CSF approach to model the steady motion of bubbles in a straight two-dimensional channel, the formation of bubbles in twoand three-dimensional T-junctions, and the breakup of droplets in three-dimensional T-junctions. A study by Heyns and Oxtoby [68] implemented a selection of surface tension modelling approaches (e.g., the CSF, a variant of Smoothed CSF and a force-balanced higher-resolution artificial compressive formulation) to model a stationary bubble. To the best of the authors' knowledge, a recent study by Yamamoto et al. [36] is the only one of its kind where different surface tension models (i.e., S-CLSVOF, density scaled S-CLSVOF and CSF) are compared based on a variety of processes with various capillary numbers (e.g., rising bubbles, capillary rise, capillary wave and thermocapillary flows).

In this study, we implemented three different surface tension models, namely CSF [66], Smoothed CSF [67] and SSF [76], in interFoam available in OpenFOAM 6. To investigate the capability of the surface tension models to handle various flow scenarios, we used two benchmark cases:


These two benchmark cases were selected due to their relevance in a variety of processes. To compare the spurious currents generated by the surface tension models, a stationary bubble was simulated. For practical applications, the time step constraint can substantially increase the computational time, thus the temporal development of the spurious currents with the surface tension models were also examined. Furthermore, the evolution of spurious currents with mesh resolution and under-relaxation factor used for the simulations was also investigated. In the interest of knowledge dissemination, the solvers and the test cases (implemented in OpenFOAM 6) discussed in the paper are available in the Supplementary Materials.

#### **2. Numerical Model**

#### *2.1. Governing Equations*

The VOF approach (developed by Hirt and Nichols [25]) denotes the individual phases using a scalar function called volume fraction, represented as

$$a\_1(\vec{x}, t) = \begin{cases} 0 \text{ (within Phase 2 or gas)}\\ 0 < a\_1 < 1 \text{ (at the interface)}\\ 1 \text{ (within Phase 1 or liquid)}, \end{cases} \tag{4}$$

where *α*<sup>1</sup> is the volume fraction of liquid. The fluid properties such as density (*ρ*) and viscosity (μ) in a control volume are calculated as

$$
\chi = \chi\_1 \mathfrak{a}\_1 + \chi\_2 \mathfrak{a}\_2 \text{ where } \chi \in [\mathfrak{p}, \mathfrak{p}], \tag{5}
$$

where *χ*<sup>1</sup> and *χ*<sup>2</sup> represent the fluid property of liquid and gas phase, respectively.

Considering the fluids to be incompressible, isothermal and immiscible, the VOF approach solves a single set of continuity and momentum equation for the whole domain. The continuity equation is written as

$$
\nabla \cdot \vec{\mathcal{U}} = 0,\tag{6}
$$

where *U*is the fluid velocity. The momentum equation is

$$\frac{\partial \rho \vec{\mathcal{U}}}{\partial t} + \nabla \cdot (\rho \vec{\mathcal{U}} \vec{\mathcal{U}}) = -\nabla p + \rho \vec{\mathcal{g}} + \nabla \cdot \mu \left(\nabla \vec{\mathcal{U}} + \nabla \vec{\mathcal{U}}^{T}\right) + \vec{F\_{st}} \tag{7}$$

where last term represents the surface tension forces, the second last term is the viscous term, *g* is the acceleration due to gravity and *p* is the pressure. Advection of the volume fraction of liquid (*α*1) is implemented in interFoam as

$$
\frac{
\partial \alpha\_1}{
\partial t} + \nabla \cdot (\alpha\_1 \vec{\mathcal{U}}) + \nabla \cdot (\alpha\_1 (1 - \alpha\_1) \vec{\mathcal{U}}\_r) = 0,\tag{8}
$$

where the third term is an artificial compression term used to sharpen the interface [58,61]. The artificial compression term uses a relative velocity (*Ur*) defined as

$$\mathcal{U}\_r = \mathbb{C}\_a \left| \frac{\Phi}{|S\_f|} \right| \vec{n}\_{f'} \tag{9}$$

where *φ*, *Sf* , *C<sup>α</sup>* and *n<sup>f</sup>* are the velocity flux, face surface area, an adjustable compression factor and unit normal vector to the interface, respectively. In the literature, *Cα* can be set between 0 and 4, where *C<sup>α</sup>* equal to zero and one correspond to the case of no and moderate compression, respectively [56]. The increase in the value of *Cα* sharpens the interface but increases the magnitude of spurious currents (see [51,56]). To model practical flow scenarios using interFoam, the value of *C<sup>α</sup>* is generally set to unity [32,63].

#### *2.2. Surface Tension Models*

This section introduces the three surface tension models: CSF, Smoothed CSF and SSF.

#### 2.2.1. The Continuum Surface Force (CSF) Model

Proposed by Brackbill et al. [66], the CSF model provides a volumetric representation of surface tension, represented as

$$F\_{st}^{\*} = \sigma \kappa \nabla \alpha\_{1\prime} \tag{10}$$

where *σ* is the surface tension and *κ* is the curvature, defined in Equation (13). The unit normal vector at the interface is calculated as

$$
\hat{m} = \frac{\nabla \alpha\_1}{|\nabla \alpha\_1| + \delta'} \tag{11}
$$

where *δ* is a small non-zero term to ensure that the denominator does not become zero. *δ* is calculated as 10<sup>−</sup>8/ <sup>∑</sup>*<sup>N</sup> Vi N* 1/3 , where *N* is the number of computational cells and ∑*<sup>N</sup> Vi* provides the sum of the volumes of individual cells (represented by *i*). Once *n*ˆ is calculated, it is corrected to account for wall adhesion through

$$
\hbar = \hbar\_w \cos \theta + \hat{l}\_w \cos \theta \tag{12}
$$

where *θ* is the contact angle of the gas–liquid interface at the walls (measured in the liquid phase), and *n*ˆ*<sup>w</sup>* and ˆ*tw* are unit vectors that are normal and tangential to the wall, respectively [82]. The curvature of the interface is then calculated as

$$
\kappa = -\nabla \cdot \mathbf{\hat{n}}.\tag{13}
$$

#### 2.2.2. The Smoothed CSF Model

The Smoothed CSF model (by Ubbink [67]) proposed modifying CSF by modifying the calculation of curvature of interface by using a smoothed volume fraction of liquid (*α*1).

The smoothed volume fraction field is calculated using a smoother proposed by Lafaurie et al. [78], which has been implemented in the literature [32,56] and is represented as

$$\widetilde{\alpha\_1} = \frac{\sum\_{f=1}^{N} < \alpha\_1 >\_{c \to f} S\_f}{\sum\_{f=1}^{N} S\_f} \, \tag{14}$$

where the indices *<sup>c</sup>* and *<sup>f</sup>* are the cell and face centre indices, respectively. <sup>&</sup>lt; *<sup>α</sup>*<sup>1</sup> <sup>&</sup>gt;*c*−→*<sup>f</sup>* represents the interpolation of *α*<sup>1</sup> from cell to face centre. The smoothening of volume fraction, done using Equation (14), is applied twice to obtain a smooth volume fraction field, which is used in Equation (15). Implementation of Equation (14) in interFoam is done using the subroutine developed in the work by [56]. Based on the smoothed volume fraction field, the unit normal to the interface is calculated as

$$
\widetilde{\boldsymbol{\hat{n}}} = \frac{\nabla \widetilde{\boldsymbol{\alpha}}\_1}{|\nabla \widetilde{\boldsymbol{\alpha}}\_1| + \boldsymbol{\delta}'} \tag{15}
$$

which is then corrected for wall adhesion (based on Equation (12)). The curvature of the interface is then calculated as

$$
\widetilde{\kappa} = -\nabla \cdot \dot{\mathfrak{h}}.\tag{16}
$$

The surface tension can be represented using the modified curvature (*<sup>κ</sup>* in Equation (16)), which can be represented as

$$
\hat{F}\_{st}^{\*} = \sigma \hat{\kappa} \nabla \alpha\_1. \tag{17}
$$

2.2.3. The Sharp Surface Force (SSF) Model

In the SSF model, proposed by Raeini et al. [76], smoothened and sharpened volume fraction fields are used to calculate curvature and gradient of of volume fraction.

The smoothened volume fraction (*αs*) is calculated based on interpolating the cell-centred values of *α*<sup>1</sup> to the cell faces using a three consecutive smoothening steps described using Equations (18a)–(18c)

$$\mathfrak{a}\_{s1} = \mathcal{C} \left\langle <\mathfrak{a}\_1 >\_{\mathfrak{c} \to f} \right\rangle\_{f \to \mathfrak{c}} + \left(1 - \mathcal{C}\right) \mathfrak{a}\_{1\prime} \tag{18a}$$

$$\mathfrak{a}\_{s2} = \mathcal{C} \left\langle < \mathfrak{a}\_{s1} >\_{\mathfrak{c} \to f} \right\rangle\_{f \to \mathfrak{c}} + \left(1 - \mathcal{C}\right) \mathfrak{a}\_{s1} \tag{18b}$$

$$\mathfrak{a}\_{\mathfrak{s}} = \mathcal{C} \left\langle <\mathfrak{a}\_{\mathfrak{s}2} >\_{\mathfrak{c} \to f} \right\rangle\_{f \xrightarrow{} \mathfrak{c}} + \left(1 - \mathcal{C}\right) \mathfrak{a}\_{\mathfrak{s}2},\tag{18c}$$

where C is set equal to 0.5. The unit normal to the interface is then calculated as

$$
\hat{n}\_s = \frac{\nabla a\_s}{|\nabla a\_s| + \delta'} \tag{19}
$$

which is then corrected for wall adhesion (based on Equation (12)). The curvature (*κs*) is calculated using Equation (19) as

$$
\kappa\_\sf s = -\nabla \cdot \mathring{n}\_\sf s.\tag{20}
$$

The interface curvature is smoothed by using a three step procedure, which can be broadly summarised into Equations (21a), (21c), and (21d). The first step involves smoothening the curvature calculated in Equation (20) as

$$\kappa\_{f1} = \left(2\sqrt{a\_c(1-a\_c)}\right)\kappa\_s + \left(1-2\sqrt{a\_c(1-a\_c)}\right)\kappa\_s^\* \tag{21a}$$

where *α<sup>c</sup>* is defined as min(1,max(*α*1,0)) and

$$\kappa\_s^\* = \frac{\left\langle \_{\varepsilon \to f} \right\rangle\_{f \to \varepsilon}}{\left\langle \_{\varepsilon \to f} \right\rangle\_{f \to \varepsilon}}, \\ w = \sqrt{a\_\varepsilon (1 - a\_\varepsilon) + 10^{-3}}. \tag{21b}$$

The second step further smoothens the curvature (calculated in Equation (21a)) as

$$\kappa\_{f2} = \left(2\sqrt{a\_c(1-a\_c)}\right)\kappa\_s + \left(1 - 2\sqrt{a\_c(1-a\_c)}\right)\kappa\_{s2}^\*, \text{ where } \kappa\_{s2}^\* = \frac{\left\langle \_{\varepsilon \to f} \right\rangle\_{f \to c}}{\left\langle \_{\varepsilon \to f} \right\rangle\_{f \to c}}.\tag{21c}$$

The final step calculates the the final curvature as

$$\kappa\_{final} = \frac{\_{c \to f}}{\_{c \to f}}.\tag{21d}$$

The surface tension is then given as

$$
\vec{F}\_{\text{sl}} = \sigma \kappa\_{\text{final}} \nabla \mathfrak{a}\_{\text{sl}\prime} \tag{22}
$$

where *αsh* is a sharpened volume fraction of liquid defined in Equation (23).

$$\alpha\_{sh} = \frac{1}{1 - \mathbb{C}\_{sh}} \left[ \min \left( \max \left( a\_{1\prime} \frac{\mathbb{C}\_{sh}}{2} \right), 1 - \frac{\mathbb{C}\_{sh}}{2} \right) - \frac{\mathbb{C}\_{sh}}{2} \right], \tag{23}$$

where *Csh* is a sharpening coefficient. A value of *Csh*=0 reduces *αsh* to *α*1, whereas *Csh*=1 provides sharp representation of the interface (which is numerically unstable). We used *Csh*=0.98 for static cases and *Csh*=0.5 for dynamic cases.

#### **3. Solver Settings**

To simplify the treatment of pressure boundary condition and density change across the interface, interFoam uses *prgh* which is defined as *p* − *ρg* · *x*, where *ρg* · *x* is the hydrostatic component of pressure [58]. The volume fraction evolution equation (Equation (8)) is solved using the Multidimensional Universal Limiter with Explicit Solution (MULES) algorithm, which preserves the boundedness of volume fraction [61,63]. Once volume fraction is solved, the continuity equation (Equation (6)) and momentum equation (Equation (7)) are solved using the Pressure Implicit with Splitting of Operator (PISO) algorithm [83]. In PISO, a predicted velocity is updated using a pressure correction procedure to advance velocity and pressure fields in time [58,63]. To understand the implementation and solution algorithm of the governing equations (Equations (6)–(8)) in interFoam, please refer to the work by Rusche [58] or Deshpande et al. [63]. The discretisation schemes, solvers and others parameters used to solve the governing equations for all the simulations discussed in this paper are presented in Tables 2–4, respectively. Under-relaxation factors, if set to less than unity, cause damping of the solution, which can lead to longer computational time for the solution reach to a steady state value. In flow scenarios where there is no steady state solution, using an under-relaxation factor can lead to erroneous results due to under-prediction of the flow variables. We used an under-relaxation factor in the solver equal to one for dynamic cases and 0.9 for static cases. The effect of using an under relaxation factor of one on static cases is also investigated.


**Table 2.** Discretisation schemes.

**Table 3.** Solvers used for the discretised equation.


**Table 4.** Other parameters used in solving the discretised equations.


#### **4. Validation: Benchmark Test Cases**

#### *4.1. Two Dimensional Rising Bubbles*

Due to the computational overhead of modelling a three-dimensional rising bubble, we model the buoyancy driven motion of a single bubble as proposed by Hysing et al. [54], Klostermann et al. [55]. The work by Hysing et al. [54] reported benchmarking data such as the bubble shape, rising velocity and circularity for two cases. These benchmarking data are produced based on numerical simulations using codes such as TP2D, FreeLIFE and MoonNMD [54]. In the work by Klostermann et al. [55], the benchmark proposed by Hysing et al. [54] was used to evaluate the VOF solver in OpenFOAM® (i.e., interFoam) for various meshes.

The computational domain used for the simulation is a rectangle of dimensions 1 m × 2 m where the bubble of diameter 0.5 m was initialised such that the centre of the bubble is at a distance of 0.5 m from the bottom and side walls. As mesh convergence could not be achieved perfectly in previous works [36,55], we used a uniform grid 160 × 320 for the simulations, corresponding to the fine mesh used in [54]. The pressure boundary conditions used in the simulations were zero gradient on the side and bottom walls, and a Dirchlet condition (equal to zero) at the top wall. The volume fraction of fluid used a zero gradient boundary condition on all walls. The velocity boundary conditions used for the simulations were no slip on top and bottom walls, but slip condition was implemented for the side walls. The fluid properties associated with the test cases, which are abbreviated as TC1 and TC2, are tabulated in Table 5. The maximum Courant number used by the solver was set equal to 0.01 and maximum time step permitted was based on Equations (2) and (3). The test cases are distinguished based on Reynolds (*Re*), Eötvös (*Eo*) and Capillary (*Ca*) numbers, which are defined as

$$Re = \frac{\mathcal{U}\_{\mathcal{S}}L}{\nu\_1}, Bo = \frac{\rho\_1 \mathcal{U}\_{\mathcal{S}}^2 L}{\sigma}, Ca = \frac{Eo}{Re} \tag{24}$$

with *L* and *Ug* being the characteristic length scale (equal to 0.5 m) and characteristic velocity (defined as |*g*|*L*), respectively. The bubble shape was obtained at *α*<sup>1</sup> = 0.5 and rising velocity was calculated based on bubble volume averaged vertical component of the velocity vector [54,55]. For validation, we used the the data reported by Klostermann et al. [55] and Hysing et al. [54] (for the predictions by the FreeLIFE solver, which is referred to as 'Benchmark' in this paper) for a uniform grid of 160 × 320.

**Table 5.** Physical parameters used for the rising bubble simulations (see [54]).


<sup>\*</sup>*g* is the reduced gravity as described in [54].

The first test case, TC1, corresponds to the case where surface tension effects are dominant [55]. The temporal evolution of the bubble as predicted by the various surface tension models is compared in Figure 1. Due to the stronger surface tension effects, the interface deforms into an ellipsoidal bubble (see Figure 2). The bubble shape (at *t* = 3 s) predicted by CSF model provides a slightly better agreement to the benchmark data compared to the other surface tension models. The surface tension models also tend to underpredict the position of the bubble at *t* = 3 s. This underprediction could be attributed to the lower rising velocity (see Figure 3), which has also been reported in previous studies using OpenFOAM [36,54,55]. Although bubble shape and rising velocity provide an overview of the capability of the surface tension models, the quantification of the errors associated with the models was based on the maximum rising velocity (*Vmax*) and the time at which the *Vmax* occurred (tabulated in Table 6). The benchmarking data show that SSF model provides a better agreement to the data reported by Hysing et al. [54] (absolute error is less than 2%) and Klostermann et al. [55] (absolute error is nearly 1.5%) in comparison to the other models.

**Figure 1.** Temporal evolution of the bubble for TC1: (**a**) *t* = 0.5 s; (**b**) *t* = 1.5 s; and (**c**) *t* = 2.5 s.

**Figure 2.** Validation Bubble shape for TC1 at *t* = 3 s: (**a**) bubble morphology; and (**b**) detailed.


**Figure 3.** Validation Bubble rising velocity for TC1: (**a**) temporal changes of bubble rising velocity; and (**b**) detailed.

The other test case, TC2, corresponds to a case where the surface tension effects are lower [55]. This results in larger deformation of interface as the bubble evolves (see Figure 4) and eventually forms a skirted bubble that has thin filaments that breaks down into smaller droplets (see Figure 5). Comparing the surface tension models to the benchmark for final bubble shape shows that the models agree quite well (see Figure 5) but there is a difference between the models with respect to the prediction of the skirted part of the bubble (see Figure 5b). Figure 6 shows that the surface tension models in comparison to the benchmark data under-predicts the rise velocity. Comparing with the benchmark, the SSF model provides the closest agreement for *Vmax*<sup>1</sup> (absolute error is nearly 3.5% [54] and less than 0.1% for [55]) and *t*(*Vmax*1) (absolute error is nearly 3–3.5% for both [54,55]) (see Table 7). On the other hand, CSF model agrees with the benchmarking data for *Vmax*<sup>2</sup> (absolute error is nearly 5.7% [54] and 0% for [55]) and *t*(*Vmax*2) (absolute error is nearly 0.6% for both [54,55]) (see Table 7).

**Table 7.** Benchmark quantities for TC2.


**Figure 5.** Validation Bubble shape for TC2 at *t* = 3 s: (**a**) bubble morphology; and (**b**) detailed.

**Figure 6.** Validation Bubble rising velocity for TC2: (**a**) temporal changes of bubble rising velocity; and (**b**) detailed.

In the previous work by Klostermann et al. [55], the spurious currents were reported to be the reason for the error between the benchmark ([54]) and their simulations (for both TC1 and TC2). Thus, the differences in the predictions, for the rising bubble simulations, between the three surface tension models considered in this paper and their departure from the benchmark can also be attributed to spurious currents generated by these models (which is discussed below). For TC2, the larger variation between the surface tension models after the first peak in the transient evolution of the rise velocity (see Figure 6) can be attributed to the differences in the shapes of filament or satellite droplets (based on the work of Yamamoto et al. [36]). Interestingly, there are also some differences in the predictions by the CSF model (for both TC1 and TC2) and the data reported by Klostermann et al. [55], which could be attributed to the difference in the solver settings (e.g., the discretisation schemes, linear solvers and number of iterations) and/or the variations within the different versions of OpenFOAM. The influence of the discretisation schemes on the predicting the flow variables has been previously investigated in [88,89] but further investigation into the effects of other solver settings (e.g., the choice of linear solver and number of iterations) on the solution is required to quantify its effect. As OpenFOAM gets updated, some of the functionalities and/or the algorithms are modified, for example, the artificial interface compression term used in advection of *α*<sup>1</sup> (defined in Equation (9)) is computed differently in the older versions of the software (see [55]). To the best knowledge of the authors, no study has reported a comparison of the performance of various versions of OpenFOAM for specific flow scenarios. These settings, especially discretisation schemes and interface compression algorithms, would effect the generation and evolution of spurious currents, which could be the potential source of the discrepancy between our simulations and the data reported in literature.

#### *4.2. Two-Dimensional Capillary Rise*

The rise of liquid through a narrow tube or between two parallel plates, which occurs as a consequence of the wetting of the walls by the liquid, is known as capillary rise. As the liquid rises, it reaches a point of equilibrium when the vertical component of the force exerted by surface tension is balanced by the gravitational force acting on the risen liquid column. This equilibrium point (for liquid rising between two vertical parallel plates) is denoted using a height (*hb*), which can be analytically calculated as

$$h\_b = \frac{2\sigma \cos \theta}{\Delta \rho |\vec{\varphi}| a'} \tag{25}$$

where Δ*ρ* is the difference between densities of liquid and gas, and *a* is the distance between the plates [90].

To study capillary rise, we used a rectangular domain of dimensions 1 mm × 20 mm, where *a* (defined in Equation (25)) is equal to 1 mm, with a uniform mesh of 20 × 400. This mesh resolution provided the most accurate prediction of capillary rise for the same computational domain while using CSF model in the previous work by Yamamoto et al. [36]. The boundary conditions for velocity field imposes a no slip boundary condition for the walls and pressure based condition (applied to both inlet and outlet) that computes inlet velocity based on the patch-face normal component of the internal-cell velocity and outflow using the zero gradient condition. The volume fraction field uses a zero gradient condition at walls (with a contact angle of 45◦) and outlet, along with a Dirchlet condition (equal to one) at inlet. The boundary condition for pressure uses a Dirichlet condition (equal to zero) at inlet and outlet whereas the walls use a Neumann boundary condition. The materials properties used for the simulations are described in Table 8. The initial volume fraction of liquid in the domain is set such that the liquid–gas interface is at a height of 8 mm from the bottom surface. The maximum time step (which satisfies both Equations (2) and (3)) and maximum Courant number were set equal to 3.5 μs and 0.1, respectively.


**Table 8.** Physical parameters used for the capillary rise simulations.

Once the interface position stabilised (see Figure 7), the capillary height *hb*,*calc* was calculated approximately from the volume fraction field as

$$
\eta\_{b,calc} = \frac{\int\_S \alpha\_1 dS}{a},
\tag{26}
$$

where the numerator is the area occupied by the liquid in the computational domain [36]. The capillary rise height calculated from the simulations is compared to the analytically derived *hb* (which was determined to be 9.9 mm using Equation (25)) in Table 9.

**Table 9.** Errors associated with the surface tension models on prediction of capillary rise.


**Figure 7.** Evolution of the water column during capillary rise.

Table 9 shows that SSF model provides a better prediction of the capillary rise height compared to CSF model. A previous work by Yamamoto et al. [36] reported an error of 7.7% for a capillary rise model using the CSF model. Interestingly, the Smoothed CSF model could not provide a reliable capillary rise prediction due to the oscillation of the water column (see Figure 7). This discrepancy can be explained based on the evolution of the spurious currents (*Usc* defined in Equation (27)), which are plotted in Figure 8. The magnitude of spurious currents (*Usc*) generated in the simulations was computed at each time step as

$$\mathcal{U}\_{\mathfrak{K}} = \max(|\vec{\mathcal{U}}|). \tag{27}$$

The periodic growth and decay of the spurious currents in the Smoothed CSF model (see Figure 8) results in the unrealistic motion of the interface whereas the CSF model which has much larger magnitude of spurious currents is much more periodic (see Figure 8), which reduces the net motion of the liquid–gas interface. Compared to CSF and Smoothed CSF models, the spurious current evolution in the SSF model is lowered by nearly two orders of magnitude (see Figure 8).

**Figure 8.** Evolution of spurious currents during the capillary rise simulations. It is worth pointing out that the figure is plotted using data extracted at every 500th point from the dataset obtained from simulations in order to reduce the rendering time of the image but care has been taken to showcase the larger temporal variations of *Usc*.

#### **5. Analysis: Spurious Current**

To study the spurious currents generated during the simulations, we simulated a stationary bubble where the effect of gravity was neglected. A bubble of diameter 2*R* was set at the centre of a square domain of dimensions 4*R* × 4*R*. The properties of the two phases and other physical parameters used for the simulations described in this section are tabulated in Table 10. For these simulations, the boundaries were assigned the Dirichlet condition, equal to 101325 Pa, for pressure and zero gradient condition for both *α*<sup>1</sup> and *U*- . The simulations were run until an end time of 0.05 s to ensure that initial transients (if any) were eliminated with maximum time step calculated based on Equations (2) and (3) along with maximum Courant number of 0.1.

**Table 10.** Physical parameters used for the simulations in the analysis of spurious current.


The accuracy of the surface tension models was calculated based on the following parameters: Laplace pressure, magnitude of spurious currents and mass imbalance. For a two-dimensional bubble, the Laplace pressure can be calculated using the Young–Laplace equation as

$$
\Delta p\_c' = \frac{\sigma}{R}.\tag{28}
$$

The Laplace pressure inside the bubble was calculated from the simulation as

$$
\Delta p\_{\varepsilon} = \frac{\int\_{V} a\_{2} p dV}{\int\_{V} a\_{2} dV} - p\_{0\prime} \tag{29}
$$

where *p*<sup>0</sup> is the operating pressure (which was equal to 101325 Pa). The mean error associated with the Laplace pressure calculated by the various surface tension models was determined as

$$\overline{E}(\Delta p\_c) = \frac{\overline{\Delta p\_c} - \Delta p\_c'}{\Delta p\_c'},\tag{30}$$

where the overbar represents the time averaged variables.

#### *5.1. Stagnant Bubble of Few Millimetres*

In this test case, we modelled a bubble with a radius of 2.5 mm using fluid properties described in Table 10 and under-relaxation factor of 0.9. The computations were performed using a uniform structured grid. The total number of mesh elements and maximum time step (which satisfies both Equations (2) and (3)) used in the simulations are described in Table 11.

**Table 11.** Details of mesh and the associated maximum time step calculated based on Equations (2) and (3) used for stationary bubble simulations.


\**R*/*δx* is the ratio of the radius of the bubble and the cell size.

To understand how spurious currents occur with various surface tension models, *Usc* is plotted at *t* = 0.05 s for the grid described by M3 in Figure 9. In the surface tension models considered in this study, the spurious currents occur around the interface but their magnitudes are much larger in the bubble than outside. To quantify the spurious currents from the simulations, the magnitude of spurious currents and capillary pressure are tabulated in Table 12. The spurious currents generated by the surface tension models tends to reduce with finer meshes for both SSF and Smoothed CSF. On the other hand, the increase in spurious current for CSF can be explained based on the dependence on the mesh size (Δ*x*) is given by

*C*Δ*<sup>x</sup>* ∼ *σ ρ*Δ*x* , (31)

where *C*Δ*<sup>x</sup>* is the magnitude of the spurious velocities (studied for CSF model [63,66]). Equation (31) indicates that smaller mesh sizes result in larger values of spurious currents for CSF model. As shown in Table 12, the Laplace pressure predicted by the surface tension models does not perfectly match Δ*p c* but both Smoothed CSF and SSF provides a better prediction in comparison to CSF.

(**c**) SSF

**Figure 9.** Comparison of spurious current generated by surface tension models at *t* = 0.05 s using M3 mesh. The gas–liquid interface in the domain is represented using a contour (in white) that is plotted at *α*<sup>1</sup> = 0.5.



#### *5.2. Effect of Time Step*

The two time step constraints were from Brackbill et al. [66] (Equation (1)) and Deshpande et al. [63] (Equations (2) and (3)). To study the effect of time step constraint, the simulations used a bubble of 2.5 mm with the M3 mesh (see Table 11) and fluid properties described in Table 10 using an under-relaxation factor of 0.9. The maximum time steps (Δ*t*) used for the simulations are 25 μs (based on [66]) and 6 μs (based on [63]).

The temporal evolution of *Usc* is compared for the surface tension models in Figure 10. Using the time step dictated by Deshpande et al. [63], the spurious currents generated by the CSF model are reduced by less than half in comparison to when time step constraint proposed by Brackbill et al. [66] was used. The other models show an absolute difference in the mean spurious current of nearly 7% and 6%, respectively, for the time step constraints (see Table 13).

**Figure 10.** Evolution of spurious currents for various surface tension models.

**Table 13.** Comparison of spurious currents for the time stepping constraints based on M3 mesh and surface tension models (while using an under-relaxation factor of 0.9).


#### *5.3. Effect of Under-Relaxation Factor*

To understand the effect of under-relaxation factor, we considered a case which used an under-relaxation factor of unity for modelling the stationary bubble of 2.5 mm with M3 mesh. Table 14 provides a summary of the spurious current and the Laplace pressure in the bubble. Comparison of the results from under-relaxation factor of 0.9 (see Table 12) and 1 (see Table 14) shows that spurious currents generated by Smoothed CSF model is substantially larger when using a larger under-relaxation factor (nearly twice). The SSF model provides the least amount of spurious currents for both the under-relaxation factors and the CSF model generates larger spurious currents with larger mesh density (as described by Equation (31)). It is also worth pointing out that the evolution of spurious currents for the time step constraints provide marginally higher spurious currents for CSF model (0.1% using the time step constraint by Equation (1)) but the Smoothed CSF and SSF models show a spurious current reduction by nearly 10% and 11%, respectively (see Table 15). Based on the evolution of spurious currents based on time step constraint, the SSF model generates the least spurious current when compared to Smoothed CSF and CSF models.

**Table 14.** Comparison of spurious currents based on mesh and surface tension models (using no under-relaxation and time step dictated by Deshpande et al. [63]).




#### **6. Conclusions**

In the study, we successfully implemented CSF, Smoothed CSF and SSF models in OpenFOAM and compared them based on their ability to simulate a two-dimensional stationary bubble, rising bubbles and capillary rise. The flow scenarios modelled corresponds to a variety of capillary numbers (in the order of 10−3, 0.1 and 1), which is relevant in various industrial processes. The numerical simulations show that:

• For a stationary bubble with a 2.5 mm radius, CSF and SSF models generate the most and least amount of spurious currents, respectively. For the finest mesh used, Smoothed CSF and SSF models reduce spurious currents by nearly one-tenth and one-twentieth of the CSF model (when no under-relaxation factor is used), respectively. When using a lower under-relaxation factor (for the finest mesh), Smoothed CSF and SSF models reduce the spurious currents by approximately one-fourth of the CSF model.


Although the surface tension models considered in this study did not eliminate spurious currents entirely, the comparison provides insights into the limitations of these models. Based on the simulations done in this study, the SSF model seems to provide a versatile surface tension formulation that generates small spurious currents and provides a more accurate representation of various processes in comparison to the standard CSF model.

#### **Supplementary Materials:** The following are available online at https://www.mdpi.com/2227-9717/7/8/542/s1.

**Author Contributions:** K.J.V. built the models, ran the simulations, post-processed and analysed the results, and wrote the manuscript. K.E.E. contributed contributed in supervising, reviewing of the results and revising manuscript.

**Funding:** We would also like to thank the Department of Material Science and Engineering, NTNU, for funding this research.

**Acknowledgments:** The authors would like to thank the OpenFOAM community (both developers and contributors). We would also like to thank the reviewers whose comments improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funding organisation had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
