A CFD model for the coupling of multiphase, multicomponent and mass transfer physics for micro-scale simulations

International Journal of Heat and Mass Transfer 113 (2017) 922–934

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A CFD model for the coupling of multiphase, multicomponent and mass transfer physics for micro-scale simulations Gim Yau Soh a,⇑, Guan Heng Yeoh a,b, Victoria Timchenko a a b

School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, NSW 2052, Australia Australian Nuclear Science and Technology Organisation (ANSTO), Locked Bag 2001, Kirrawee DC, NSW 2232, Australia

a r t i c l e

i n f o

Article history: Received 9 January 2017 Received in revised form 3 April 2017 Accepted 1 June 2017

Keywords: Volume-of-fluid (VOF) Multiphase flow Multiphase and multicomponent coupling Mass transfer Micro-scale

a b s t r a c t This paper presents a CFD model that couples the physics of multiphase and multicomponent flow involving mass transfer for micro-scale simulations. The volume-of-fluid (VOF) model is utilised to capture the multiphase physics, and smoothing operations are applied in the surface tension force calculation to minimise spurious velocities. For multicomponent tracking of the species transport, the a-factor and expulsion operation methods are implemented to ensure that the species are tracked within the correct phase. Three different mass transfer models were adopted with modifications to the original models. The developed coupled multiphase-multicomponent model was assessed via a series of validation cases: a simple one-dimensional (1D) diffusion problem, horizontal vapour flow over a smooth stationary liquid where Sherwood number was used to measure mass transfer performance, and the dissolution of a twodimensional (2D) droplet. In all the test cases, the model was able to yield reasonably close results with the analytical solution or empirical correlations used in the validation, thus establishing its accuracy and reliability of the developed model. Simulations using the coupled multiphase-multicomponent model represent a cost-effective approach to obtaining insights into the flow physics for a variety of applications: dissolution of droplets in microchannels, transport of drugs in blood vessels that entails mass transfer and simulations of carbon dioxide enhanced oil recovery at the pore-scale. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Many processes and applications in the industry involve fluid systems that encompass various complex flow physics: different phases co-existing, mass transfer and multiple chemical species. These systems prevail across various scales, from large-scale bubble column reactors used in chemical, petrochemical, biochemical and metallurgical industries [1] to the micro-scale transport of droplets in microchannels for drug delivery applications [2]. Numerical simulations can provide effective means of studying the flow behaviour and obtaining crucial insights into the design of such devices. However this can only be made possible if the numerical model used is able to capture all the complex flow physics aforementioned. Much advancement has been made in the development of multiphase computational methodologies: level-set method [3,4], front-tracking method [5] and the Shan-Chen model for multiphase modelling using the lattice Boltzmann method [6,7]. These models have subsequently been used to investigate problems ⇑ Corresponding author. E-mail address: [email protected] (G.Y. Soh). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.06.001 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

involving mass transfer. Yang and Mao [8] adopted the level-set approach to simulate interphase mass transfer and validated their model using experimental results of a rising deformable droplet. Koynov et al. [9] applied the front-tracking method to investigate the mass transfer and flow dynamics of bubble swarms in order to study their impact on chemical reactions. The most widely used interface capturing model is he volumeof-fluid (VOF) method which was first developed by Hirts and Nichols [10]. The model utilises a VOF function a to track the volume fraction of the different phases, and a Piecewise-Linear Interface Calculation (PLIC) method [11,12] to approximate the multiphase interface. Many authors have successfully implemented mass transfer models for multiphase systems that accommodate the VOF model as the base model [13–15]. Furthermore, the VOF model has also been used to investigate simultaneous heat and mass transfer for a variety of problems, including bubble growth during boiling in a microchannel [16,17], heat and mass transfer in an inclined channel [18], nucleate boiling in microchannels [19,20] and evaporation [21]. In these numerical models, the issue of spurious velocities persists which are the non-physical velocities that arise due to discretisation errors in the surface tension force calculation. This problem is exacerbated at the

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

micro-scale as capillary forces dominate, which may compromise the accuracy of simulations at the micro-scale. Recently, innovative approaches have been developed to minimise the issue of spurious velocities in the VOF model [22–25] and have been successfully used to investigate micro-scale problems involving pore-scale [26,27] and microchannel [28] flows. In addition, the VOF model has also been adopted to simulate the physics of multiphase flow with multicomponent species tracking in the presence of mass transfer simultaneously. Bothe and Fleckenstein [29] implemented a mass transfer model that is calculated based on the product of the diffusion coefficient and species concentration gradient. The model was assessed for the simulation of oxygen transfer from rising three-dimensional (3D) bubbles and obtained good agreement with experimental data. Haelssig et al. [30] developed a model that computes the mass transfer rate based on the conservation of species from the interface jump in the presence of heat transfer, and validated their results with empirical Sherwood number correlations. Other models have also been proposed that account for the discontinuity in species concentration across the multiphase interface [31–34]. The practicality of these models have been demonstrated through their usage in studying a variety of applications, including the dissolution of CO2 droplets in methanol and ethanol in a microchannel [35], a membrane separation process [36] and determining the mass transfer and liquid hold-up in structured packing [37]. Much development have been performed on the interface capturing models based on VOF coupled with the tracking of multicomponent species, as evidenced by the aforementioned literature. Nonetheless, many of these models fail to address the issue of spurious velocities which will encumber the reliability of results for multiphase-multicomponent simulations at the microscale. Furthermore, most models employ a single species transport equation to track the flow of a species type in all phases. This is useful in providing the overall flow behaviour of the species but may prove difficult when the behaviour of the species in a specific phase are required to be locally predicted. For instance, for the problem of a dissolving CO2 droplet within silicone oil flowing in a microchannel, one might be interested in studying the flow physics of the dissolved CO2 within the silicone oil, and in this case having a species transport equation that solves specifically for the transport of dissolved CO2 within silicone oil will be more relevant. In this paper, a coupled multiphase-multicomponent model that is applicable for micro-scale problems is presented. In the next section, the development of the model with all the required physics is described in detail: multiphase governing equations, method to minimise spurious velocities, mass transfer and calculation of interfacial area, multicomponent and tracking of species within the correct phase. In Section 3, three validation cases are employed to test the accuracy and reliability of the coupled multiphasemulticomponent model.

2. Mathematical and numerical formulation 2.1. CFD governing equations and multiphase modelling The coupled model developed in this work requires the simultaneous simulation of both multiphase and multicomponent physics. For the multiphase modelling, we adopt the interFoam volume-of-fluid (VOF) solver in OpenFOAM 2.3. The VOF multiphase model was first proposed by Hirts and Nichols [10] and uses a VOF function a to track the volume fraction occupied by each phase in the individual mesh cells. The model is assumed to be Newtonian, incompressible and isothermal. The isothermal assumption is valid for dissolution mass transfer with minimal

923

temperature variation (e.g. dissolution of CO2 in silicone oil). The density q and viscosity l are calculated via the a-weighted average of the phases occupying each cell:

q ¼ a1 q1 þ a2 q2

ð1Þ

l ¼ a1 l1 þ a2 l2

ð2Þ

where the subscripts 1 and 2 refer to the different phases under consideration, with a1 ¼ a and a2 ¼ 1  a. When mass transfer is included, a source term is added to the right hand side of the atransport equation:

@a 1 Ja þ r  ðaVÞ ¼ @t q1 VOF

ð3Þ

where V is the velocity vector, J is the mass transfer flux rate coefficient (to be discussed in Section 2.5) and aVOF is the interfacial area density, which will be further elaborated in Section 2.4. The atransport equation is then solved together with the continuity and momentum equations to fully define the flow field:




rV ¼

1

q1



1

q2

 JaVOF

@ qV þ rðqV  VÞ ¼ rp þ qg þ r  lðrV þ rVT Þ þ Fr @t

ð4Þ

ð5Þ

where p, g, and Fr are the pressure, gravity, and surface tension force, respectively. The right hand side term in Eq. (4) accounts for the fluid expansion/contraction due to the mass transfer whilst maintaining mass conservation. The surface tension force Fr is modelled using Brackbill’s [38] model of the Continuum Surface Force (CSF):

Fr ¼ rkn

ð6Þ

In Eq. (6), r is the interfacial tension, n ¼ ra=jraj is the interface normal which is non-zero only at the interface, and k ¼ r  n is the interface curvature. From a physical point of view, the expression for curvature k captures the degree of change in the interface normal vector n along the surface of the multiphase interface, which describes how curved the interface is (e.g. a sharp change in normal vector n over a short distance implies that the interface is highly curved). 2.2. Smoothing operations to alleviate spurious velocities in multiphase modelling Although the VOF model is able to simulate multiphase problems requiring the tracking of the interface, nevertheless as with all multiphase numerical models it suffers from the issue of spurious velocities or parasitic currents. These are non-physical velocities that arise due to the numerical errors from the discretisation of the interface curvature, which is used to calculate the surface tension force [39]. This error in the surface tension then propagates into the momentum equation and necessitates the velocity terms to balance them, resulting in the spurious velocities. The problem is exacerbated at the micro-scale where capillary forces dominate [24] and must therefore be addressed to simulate accurate flow behaviour at the micro-scale. To minimise the spurious velocities, we adopt a series of smoothing operations when computing the interface curvature. Firstly, the a-field is smoothed using a recursive interpolation between the cell and face centres for each grid cell [22,23]. The smoothing operation works by interpolating the cell centre avalues into the face centre values, and then interpolating from the new face centre values back into the cell centre values, recursively. This smoothing operation is applied twice to the a-field:

924

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

aiþ1 ¼ C s ½½ai c!f f !c þ ð1  C s Þai

ð7Þ

where C s is the smoothing coefficient (set to be 0.5). The subscript ½ c!f denotes interpolating from cell centre to face centre and vice versa. Note that the a-smoothing is only applied to the a values used to calculate the interface curvature, and not the a-field itself. The curvature k is then calculated using the smoothed a values, after which it goes through further smoothing via a similar recursive interpolation [24]. The smoothing for curvature k reduces the spurious velocities because it diffuses the k-values that are present at regions away from the interface (where a is non-zero but is close pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi to 0 or 1), through the multiplication by a factor of 2 að1  aÞ, which reduces the numerical errors that result in spurious velocities [24]. Similar to the a smoothing, the smoothing for k is also applied twice:

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ks;iþ1 ¼ 2 að1  aÞk þ ð1  2 að1  aÞÞðks Þ; 

ks ¼

½½ks;i wc!f f !c ½½wc!f f !c

;

w¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi að1  aÞ þ 106

ks;0 ¼ k

ð8Þ ð9Þ

The effectiveness of the smoothing operations described has been rigorously demonstrated through various simulations of microchannel flow in our previous works [25,28]. 2.3. Multicomponent coupling To describe the formulation of multiphase-multicomponent coupling, the example of a fluid system comprised of CO2 gas and silicone oil liquid (represented by phases 1 and 2 respectively) is considered. Capturing the fluid behaviour of such a system requires the modelling of both its multiphase (CO2 gas coexisting with silicone oil liquid) and multicomponent (CO2 dissolved in the silicone oil) characteristics. From a numerical simulation point of view, the multiphase and multicomponent behaviour can be captured by individually solving the VOF and species transport equations respectively. However when simulating both the multiphase and multicomponent physics simultaneously, there are some considerations that need to be taken into account in order to develop a holistic numerical model: 1. Mass must be conserved for the species studied during mass transfer. For instance to capture the dissolution of CO2 gas into the silicone oil, mass transfer from the CO2 gas phase (tracked via VOF) needs to be added as a source term for the dissolved CO2 in the liquid phase (tracked via species transport). The rate of mass transfer is determined by interface jump conditions to such that there is chemical equilibrium between the species (to be discussed in Section 2.5). It is worth mentioning that interface jump conditions for energy must be considered as well if thermal effects are to be studied, but will not be further discussed as we assume isothermal conditions in this work. 2. This also means that the dissolved species must be tracked only within the phase it exists in. For example, the solution to the species transport equation for dissolved CO2 concentration must only have non-zero values within the liquid phase and not the gas phase. Otherwise, having a non-zero solution in the gas phase would imply the existence of dissolved CO2 within the CO2 gas phase itself, which is non-physical. To ensure that the dissolved species are tracked within the correct phase, the species transport equation is multiplied by a factor of a, which is the volume fraction of the correct phase. This ensures that the advection and diffusion terms are only non-zero when solved for the correct phase. In the incorrect phase, the advection

and diffusion terms will be zero since a ¼ 0 for the phase not under consideration. The species transport equations for the vapour and liquid phases multiplied by the a factor for the respective phases can be presented as:

@ ðaV qV xa;V Þ þ r  ðaV qV Vxa;V Þ @t ¼ r  ðaV qV Dab;V rxa;V Þ þ Sa;V @ ðaL qL xa;L Þ þ r  ðaL qL Vxa;L Þ ¼ r  ðaL qL Dab;L rxa;L Þ þ Sa;L @t

ð10Þ ð11Þ

where xa refers to the mass fraction of species a within the vapour (subscript V) and liquid (subscript L) phases. The mass fraction for species b can be obtained using the simple relation xb ¼ 1  xa for the respective phases. Dab is the diffusivity of species a in species b. The mass transfer source terms Sa;V and Sa;L are not multiplied by a as they are zero in both phases and non-zero only at the interface where mass transfer occurs. These source terms are further discussed in Section 2.5. Nonetheless, employing the a factor alone is insufficient, as cells containing the incorrect phase that are immediately adjacent to the interface cells will still experience some advection and diffusion of the incorrect species into it. This is because OpenFOAM uses the Gauss’s theorem to calculate the gradient terms, where the face centre values are taken into account. For the cells containing the incorrect phase that are adjacent to interface cells, one or more of the faces will have a non-zero face a value (Fig. 1). When it is used to calculate the advection and diffusion gradients in the species transport equation, it will result in the transport of the incorrect species into the incorrect phase (e.g. dissolved CO2 into CO2 gas phase). In order to circumvent this problem, an ‘‘expulsion” operation is implemented, whereby the species that is transported into the cells containing the incorrect phase are ‘‘expelled” out into the adjacent interface cells (Fig. 2). The rationale for adopting such an operation is that since the species is incorrectly transported into the cells containing the incorrect phase from the interface cells in the first place, the expulsion operation merely reverses this by expelling the species back into the adjacent interface cells which is where they came from. The implementation of this expulsion operation in the solution algorithm is discussed in Section 2.6, and its effects on the accuracy of the results are examined in

Fig. 1. Schematic showing a incorrect-phase cell (a ¼ 0) having a non-zero aface , which results in an incorrect species transfer into the cell.

925

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

a < 0 or a > 1, which violates the volume fraction function a. This will also result in inaccurate jraj values for the following timestep, causing the subsequent mass transfer to be inaccurate as well. In the authors’ previous work [43], an algorithm was developed to obtain an improved estimate to the interfacial area for mass transfer. The algorithm adopts a similar concept to the Piecewise-Linear Interface Calculation (PLIC) method [12,44] by considering a plane cutting the mesh cell. The position and orientation of the plane cut are first determined from the volume fraction a and the interface normal n. Then the cut surface is used to approximate the interfacial area for the interface mesh cell of interest, using the algorithm aforementioned. For a detailed description of the algorithm, readers are kindly directed to the authors’ previous work [43], which also discusses various validation cases that established the accuracy of the algorithm. 2.5. Mass transfer models Fig. 2. Schematic illustrating the expulsion operation of the species from the cell containing the incorrect phase cell back into the adjacent interface cells.

Section 3.3. Finally, the species mass concentration C a is obtained using the following relations:

C a;V ¼ aV qV xa;V C a;L ¼ aL qL xa;L

ð12Þ

2.4. Interfacial area for mass transfer, aVOF Since mass transfer occurs across the multiphase interface, it is imperative that the interfacial area is calculated correctly in order to capture accurate mass transfer rates in numerical simulations. A common method found in literature has been to determine the interfacial area aVOF for mass transfer through the jraj function [16,18,32,40]. Indeed, the jraj function does represent the multiphase interfacial area on the global level across the entire compuR tational domain ( jrajdV ¼ Aint where Aint is the total interfacial area) [21,41]. Nonetheless at the local level, the jraj function does not represent the interfacial area adequately [42,43] and thus cannot be used to calculate the local mass transfer rates. As shown in Fig. 3, the jraj function can have non–zero values for the mesh cells directly adjacent to the interface cells since the gradient of a for the adjacent cells is normally calculated using the a values of the neighbouring cells including the interface cells (where 0 < a < 1). This leads to mass transfer occurring at the adjacent cells where a ¼ 0 and a ¼ 1, and can thus lead to a scenario where

In any multiphase-multicomponent systems, the chemical species will always move towards equilibrium. The state of nonequilibrium is what causes the mass transfer, as it is the mass transfer that drives the system towards equilibrium. Many mass transfer models used in numerical simulations are formulated based on this premise, and the mass transfer rate is calculated by considering the rate required for the system to reach equilibrium. The different models each have their distinct formulation but the aim of all models is to ultimately work out the mass transfer flux rate coefficient J. In this work, three mass transfer models are adopted: (i) basic model with constant mass transfer coefficient, (ii) mass transfer model by Bothe and Fleckenstein [29], and (iii) mass transfer model by Haelssig et al. [30]. These models are used and assessed in the validation test cases of the multiphase-multicomponent model. 2.5.1. Constant mass transfer coefficient This model is straightforward and only involves prescribing a constant mass transfer flux rate coefficient J in the source term for mass transfer (JaVOF ). The J values can be obtained from experimentally measured values for the fluid system being studied. Whilst easy to implement, this model is disadvantageous in that it does not account for the change in mass transfer rate as the concentration of the species changes. This is because the experimentally measured J values are usually averaged over a time period and therefore only captures the overall rather than instantaneous mass transfer behaviour. The model is still suitable however to study systems where the mass transfer rate is low and the species concentration is dilute, or when the overall mass transfer behaviour is of more interest than the local mass transfer physics. 2.5.2. Bothe and Fleckenstein’s model In Bothe and Fleckenstein’s [29] model, the mass transfer flux rate coefficient is determined by the concentration gradient at the interface:

J ¼ Dab rC int

ð13Þ

To calculate the concentration gradient rC , Bothe and Fleckenstein implemented an algorithm that computes two concentration values in the acceptor phase in the normal direction and the gradient is approximated using these concentration values. A different approach is utilised in this present work. Firstly, the concentration gradient is computed in the x; y and z directions to int

Fig. 3. An example of the distribution of volume fraction a across an interface. jraj values are non-zero at non-interface cells [43].

int int int yield rC int is comx , rC y and rC z , then the final gradient rC puted using an interface normal-weighted sum of the component gradients:

926

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

int int rC int ¼ nx rC int x þ ny rC y þ nz rC z

ð14Þ

For the calculation of the component gradients, a second order backward differencing scheme is used, by taking two concentration points away from the interface (Fig. 4). The example for the x-component gradient calculation is as follows:

rC int x ¼

3C int  4C x;1 þ C 0x;2 2 Dx

ð15Þ

Eqs. (20) and (21). To aid the explanation of the calculation, Fig. 5 is included as a reference, which illustrates a 2D case but the same method applies for the z-direction component. The idea is to use the mass fraction and coordinates of the interface’s neighbouring cells in the respective phases to calculate the mass fraction gradients:

xa;Vjiþ1;j  xint a;V

rxa;V  n ¼

where Dx is the distance between the interface centroid and centroid of the immediate neighbour (Cell 1). The value for C 0x;2 is obtained by linearly interpolating between C x;1 and C x;2 using a dis-

r

þ

xint ¼

p kH

C int ¼ aqxint

ð16Þ ð17Þ

Once J is determined, the source terms for the species transport equations can be calculated via Sa;V ¼ JaVOF and Sa;L ¼ JaVOF . Bothe and Fleckenstein’s model only considers mass transfer for one of the species, and is useful when investigating problems where mass transfer predominantly involves only one species (e.g. CO2/silicone oil system, where mass transfer is primarily the dissolution of CO2 only). For more complex systems where more than one species in the system experiences mass transfer, Haelssig et al.’s [30] model is more appropriate and is covered in the next section. 2.5.3. Haelssig’s model The mass transfer model by Haelssig et al. [30] is formulated by considering the jump conditions to ensure species mass conservation and chemical equilibrium across the multiphase interface. To enforce the conditions in the model, a series of governing equations are setup. These equations are then solved simultaneously to obtain the mass transfer rates for all species in the system. Species mass conservation in both liquid and vapour phases are satisfied via the following equations:

_ a þ qV Dab;V rxa;V  n  xint _ _ m a;V ðma þ mb Þ ¼ 0

ð18Þ

_ _ _ a þ qL Dab;L rxa;L  n  xint m a;L ðma þ mb Þ ¼ 0

ð19Þ

þ

zk1  zc

!

yi;jþ1  yc

 ny ð20Þ

 nz

!

xa;Ljk1  xint a;L

xa;Vji;jþ1  xint a;V

!

zkþ1  zc

xi1;j  xc

int

The species concentration at the interface C is calculated using Henry’s law, where the Henry’s law constant kH is obtained using experimentally measured values for the fluid system that is being studied:

xa;Vjkþ1  xint a;V

xa;Lji1;j  xint a;L

rxa;L  n ¼

are also determined using the same approach.

 nx þ

xiþ1;j  xc

tance of Dx from Cell 1. The other component gradients rC int y and C int z

!

xa;Lji;j1  xint a;L

 nx þ

yi;j1  yc

!

!  ny ð21Þ

 nz

where xc , yc and zc are the coordinates for the interface centroid and is used to work out the distance between the centroid and the neighbouring cells for the gradient calculation. In Haelssig et al.’s original work the centroid coordinates (xc , yc , zc ) were not used, but instead the coordinate is first extrapolated along the tangent of the interface up to the point where the extrapolated point, the interface cell centroid and the neighbouring cell centroid are all collinear. The distance used in the gradient calculation is then taken as the distance between this extrapolated point and the neighbouring cell centroid (see Haelssig et al.’s Fig. 1 for reference). This approach was not adopted herein as the extrapolated point in some instances lie outside the interface cell, when the interface is either highly shallow or steep. Further, using the interface centroid to directly calculate the distance for the gradient calculation is a reasonable approximation and it also precludes additional computational effort to calculate the extrapolated point. The chemical equilibrium condition is enforced via the following equations:

xint a Ka



yint a

¼ ¼0

wint a;L =M a int wint a;L =M a þ wb;L =M b

!

K1 þ

wint a;V =M a

!

int wint a;V =M a þ wb;V =M b

ð22Þ

_ a and m _ b refer to the mass transfer rates for species a and b where m respectively, and is positive when the species is condensing and negative during evaporation. The terms rxa;V  n and rxa;L  n refer to the mass fraction gradients at the interface and are calculated via

Fig. 4. Schematic showing implementation of the second order backward difference to calculate the component gradients.

Fig. 5. Structured grid showing the interface cell (i; j), the interface centroid (xc ; yc ; zc ), and the neighbouring cells.

927

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

int xint b K b  yb ¼

wint b;L =M b int wint a;L =M a þ wb;L =M b

! K2 þ

wint b;V =M b int wint a;V =M a þ wb;V =M b

! ¼0 ð23Þ

int xint a;V þ xb;V  1 ¼ 0

ð24Þ

int xint a;L þ xb;L  1 ¼ 0

ð25Þ

where x and y are the liquid and vapour phase mole fractions, K is the equilibrium constant and M is the molar mass. The six governing equations (Eqs. (18), (19), and (22)–(25)) are solved _ 2 , which are then added _ 1 and m simultaneously to determine m as volumetric source terms in the conservation equations via the following:

_ a aVOF Sa;V ¼ m

ð26Þ

_ a aVOF Sa;L ¼ m

ð27Þ

_ aþm _b J¼m

ð28Þ

2.6. Solution algorithm and numerical settings For the solution algorithm, the Pressure Implicit with Splitting of Operator (PISO) scheme by Issa [45] is adopted for pressurevelocity coupling. To ensure boundedness for the VOF function a, the Multi-dimensional Limiter for Explicit Solution (MULES) algorithm by Rusche [46] is adopted. The transient terms are discretised using the first order implicit Euler scheme, while the spatial discretisation is implemented using a combination of second order upwind and central differencing schemes via the van Leer limiter [47]. The convergence criteria used is 106 . The coupled multiphase-multicomponent model involves complex physics with numerous governing equations as described in Sections 2.1–2.5. The order in which all the equations are solved for each time-step is detailed as follows: 1. Calculate interfacial area aVOF using the algorithm described in Section 2.4. 2. Calculate mass transfer flux rate coefficient J using one of the selected mass transfer models in Section 2.5. 3. Solve for a transport (Eq. (3)). a. Update fluid properties based on new a distribution (Eqs. (1), (2)). 4. Solve for species transport (Eqs. (10) and (11). a. Perform expulsion operation described in Section 2.3. 5. Perform smoothing operations (Eqs. (7)–(9)) and calculate surface tension force Fr using smoothed variables (Eq. (6)). 6. PISO loop to solve for pressure and velocity. 7. Proceed to next time-step.

3.1. Pure diffusion validation The first test case conducted is a simple 1D diffusion case to validate the species transport solver of the coupled multiphasemulticomponent model. The analytical solution for a semiinfinite domain with a fixed concentration C o at one boundary and initial concentration being zero throughout the domain is given as [48]:


C ¼ 1  C o erf

x pffiffiffiffiffiffi 2 Dt

 ð29Þ

where the erf expression in Eq. (29) is the error function. The 1D domain used in our simulation has a length of 1 m, and the boundary conditions applied are a fixed value of C o ¼ 1 at the left end and zero species flux on the right end. To determine the appropriate mesh density for the domain, concentration profiles are obtained using 50, 100, 200 and 400 cells as shown in Fig. 6. The results show that all mesh densities yielded close results with one another which demonstrates convergence. A mesh density of 200 cells is therefore used in subsequent simulations. The simulation concentration profiles obtained using our model are compared with the analytical solution for both transient (different timeframes, Fig. 7) and varying diffusion coefficient (Fig. 8) scenarios. For a higher timeframe (t = 0.09 s) and higher diffusion coefficient (D = 2 m2/s), the tail-ends start to diverge by a small amount as the profile approaches the right end. This divergence is caused by the zero species flux applied to the right end of the simulation domain, which induces a zero concentration gradient at the right end of the domain. As a result, the concentration profiles from the simulation ‘‘bend-up” slightly to satisfy the zero gradient condition. This problem can be resolved by using a longer simulation domain (e.g. L ¼ 2) so that the zero gradient condition is only applied at x ¼ 2, whilst the accuracy of the results at x = 1 is maintained. Nevertheless for the sake of consistency, a domain length of L ¼ 1 is maintained for the results in this paper. Furthermore, the majority of the results exhibit a close match between the simulation and analytical results. Therefore this validation case establishes the accuracy of the species transport solver of our model. In the next section, we consider a more complex validation case which assesses the accuracy of the model not just in terms of the species transport but also the multiphase mass transfer within a 2D context. 3.2. Validation with local Sherwood number The accuracy of the coupled multiphase-multicomponent model is evaluated by comparing its simulated Sherwood number

The accuracy of the coupled multiphase-multicomponent model is assessed via a series of test cases which are covered next in Section 3. 3. Results and discussion Three validation cases are employed to test the accuracy and reliability of the coupled multiphase-multicomponent model: a simple 1D diffusion problem, horizontal vapour flow over a smooth stationary liquid where Sherwood number was used to measure mass transfer performance, and the dissolution of a 2D droplet.

Fig. 6. Concentration profiles ðD ¼ 1 m2 =s; t ¼ 0:05 s).

obtained

using

different

mesh

densities

928

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

and Fleckenstein (Section 2.5.2) and Haelssig et al.’s (Section 2.5.3) models were tested and the resulting Sherwood numbers simulated are compared with the following empirical correlations [30]:

Correlation 1 : ShDh;V ¼ 7:54 þ

0:03ðReDh;V ScDh =LÞ 1 þ 0:016ðReDh;V ScDh =LÞ2=3

Correlation 2 : ShDh;V ¼ 1:85ðReDh;V ScDh =LÞ1=3

Fig. 7. Transient concentration profiles compared with the analytical solution (D ¼ 1 m2 =s).

ð30Þ ð31Þ

Mass transfer occurs as the liquid pentane evaporates into the CO2 gas that is flowing over the liquid-vapour interface. The transient simulations were conducted until the composition and velocity at the outlet observed negligible change, which is defined as the pseudo steady-state. Simulations with different Schmidt numbers arrive at steady-state at different times, but to maintain consistency the simulations were carried out over a real flow time of t = 0.3 s which is more than sufficient duration for steady-state to be achieved. To calculate the Sherwood number obtained from the numerical simulation, the same method as proposed by Haelssig et al. [30] is adopted. The Sherwood number for the vapour phase is only considered as it is where most of the mass transfer and species transport occurs. The local mass transfer coefficient is determined via:

K V;local ¼ 

Dab;V rxb;V  n ðxint b;V  xb;V Þ

ð32Þ

xint b;V is the equilibrium vapour pentane mass fraction at the interface and xb;V is the bulk pentane mass fraction in the vapour phase. The average mass transfer coefficient for the entire domain can be calculated as: KV ¼ Fig. 8. Concentration profiles for different diffusion coefficients (t ¼ 0:05 s).

with empirical correlations found in literature. To this end, the case of vapour flows horizontally over a smooth stationary liquid is investigated, which us similar to the setup used by Haelssig et al. [30] in their validation. This setup is a good test evaluation case because it involves the physics of multiphase (vapour and liquid), mass transfer (evaporation of liquid into vapour phase), and multicomponent (species of evaporated liquid within original vapour phase). The Sherwood number which is the ratio of convective mass transfer rate over diffusion rate is able to assess the accuracy of the model in terms of the mass transfer models and the multicomponent species tracking. The geometry used for the computational domain together with the boundary conditions applied are shown in Fig. 9. The fluids used are CO2 for the vapour and n-pentane for the liquid, and the fluid properties are summarised in Table 1. The viscosity of the liquid pentane is artificially increased to prevent shear induced waves from forming along the vapour-liquid interface. This modification is acceptable since it does not affect the multicomponent or mass transfer physics being studied, nor does it influence the value of the Sherwood number correlations. A structured grid was used for the domain and convergence was achieved using 60  600 cells. The vapour phase was initialised with pure CO2, while the liquid phase was initialised with pure pentane. A semi-parabolic velocity profile which represents the steady-state velocity was applied at the inlet, with an average velocity magnitude of 1 m/s. This setup was simulated for a Schmidt number (Sc ¼ lV =qV Dab;V ) range of 0:4 < Sc < 1:2, which corresponds to diffusion coefficients used of 0:393  105 < Dab;V < 1:18  105 m2 =s, which are within the diffusion coefficient range of typical gases [49]. Both the Bothe

1 L

Z

L

K V;local dx

ð33Þ

0

where L is the length of the vapour-liquid interface. Finally, the vapour phase Sherwood number can be computed:

ShDh;V ¼

K V Dh Dab;V

ð34Þ

Fig. 10 shows the parity plot of the vapour phase Sherwood number, comparing the simulation results obtained using Bothe and Fleckenstein’s and Haelssig et al.’s models to the predictions using Correlations 1 and 2. The plot indicates a close relationship between the simulations and correlations, with an average deviation of around 5% (Table 2). The good agreement in the Sherwood number results indicates accurate capturing of the mass transfer performance and transport of species, and therefore establishes the reliability of the coupled multiphase-multicomponent model. We also present the concentration profiles obtained from both mass transfer models across different Schmidt numbers (Fig. 11). It can be observed that the concentration at the free-stream increases as the Sc increases. In addition, the concentration gradient for Bothe and Fleckenstein’s model is slightly steeper near the interface, which explains why it yielded marginally higher Sherwood numbers compared to Haelssig et al.’s model (Fig. 10). To better understand the local mass transfer performance, the values for the local mass transfer coefficient K V;local are also plotted as shown in Fig. 12. The results show that the K V;local values are greater for higher Schmidt numbers. Furthermore, the K V;local values obtained using both the Haelssig et al. and Bothe and Fleckenstein’s models closely correlate to one another. An interesting behaviour observed is the wavy nature of the K V;local values as it approaches the right end of the domain (x=L ¼ 1). This is caused by the minor shear induced waves along the vapour-liquid interface, which causes the vaporised n-Pentane concentration profile to be slightly wavy along the x-direction (Fig. 13). The shear waves

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

929

Fig. 9. Geometry for computational domain and the boundary conditions applied for the Sherwood number validation case.

Table 1 Fluid properties of CO2 and n-pentane for the operation conditions of p ¼ 172 kPa and T ¼ 294 K [49].

qV ðkg=m Þ lV ðPa  sÞ 3

M ðg=molÞ K ¼ y=x [50]

CO2

n-Pentane

3:10

5:08

1:47  105 44:01 41:02

6:68  106 72:15 0:329

Fig. 11. Concentration profiles obtained using Haelssig et al.’s [30] and Bothe and Fleckenstein’s [29] models.

Fig. 10. Parity plot of the vapour phase Sherwood number obtained using Haelssig et al.’s [30] and Bothe and Fleckenstein’s [29] models, compared with predictions using Correlations 1 and 2 (Re ¼ 2540 and 0:4 < Sc < 1:2).

Table 2 Average deviations between the simulation and correlation Sherwood numbers. Comparison

Average deviation (%)

Haelssig with Correlation 1 Haelssig with Correlation 2 Bothe with Correlation 1 Bothe with Correlation 2

5.03 4.63 4.87 4.52

are minimal in the entrance region (x=L ¼ 0) but gradually builds up as it travels along the interface and is the strongest near the exit region (x=L ¼ 1), which explains the greater ‘‘waviness magnitude” of the K V;local values towards the right end of the domain. 3.3. Dissolution of a 2D droplet The dissolution of a CO2 droplet within silicone oil in a 2D scenario is considered for the validation of the coupled multiphase-multicomponent model. Similar to the Sherwood number investigation in the previous section, this problem involves the physics of multiphase (CO2 gas co-existing with liquid silicone oil), multicomponent (dissolved CO2 within silicone oil) and mass transfer (rate of CO2 dissolution). Furthermore, the size of the domain is in the micrometre scale, in which the smoothing

Fig. 12. Plot of the local mass transfer coefficient K V;local values taken along the vapour-liquid interface.

operations adopted (Section 2.2) will play an important role in ensuring the accuracy of the simulations. In light of these considerations, this problem is able to assess the efficacy of all the physics incorporated into the coupled multiphase-multicomponent model holistically. To create the CO2 droplet, we use a T-junction microchannel geometry as depicted in Fig. 14. Silicone oil and CO2 inlets are set at the main and side channels respectively, with velocities of V main ¼ 0:0104 m=s and V side ¼ 0:0156 m=s. The droplet is formed when the silicone oil and CO2 interacts immiscibly with one another at the T-junction. The fluid properties of CO2 and silicone oil used are shown in Table 3. The contact angle is obtained by measuring from the experimental images of CO2 droplets within silicone oil in a microchannel [51]. In our simulations, both the Bothe and Fleckenstein’s (Section 2.5.2) and constant mass transfer coefficient models

930

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

Fig. 13. (Left) Magnified contour plot of the vaporised n-Pentane concentration distribution, taken at the right end of the computational domain. (Right) Concentration profile for the vaporised n-Pentane taken at y=ymax ¼ 0:45.

(Section 2.5.1) are employed. A structured grid for the domain is used and tested for convergence by checking the simulated droplet lengths and prediction of mass transfer coefficient J (Table 4), using mesh densities of 36, 48, 60 and 72 cells per D (D ¼ 800 lm). Convergence was found to have been reached at 60 cells per D and is therefore deemed as the final mesh density used for subsequent simulations. For the purposes of validating the simulation results, we compare it to the analytical solution for the length of a 2D dissolving droplet over time, which is derived in the following. The length of the dissolving CO2 droplet as a function of time can be expressed as:

Ldroplet ðtÞ ¼ Linitial  Ldissolv ed ðtÞ ¼ Linitial 

V dissolv ed ðtÞ Dlz

ð35Þ

The term lz is to account for the thickness in the z-direction, and is included to aid in the comprehension of the derivation albeit the problem being 2D. Dlz is the cross-sectional area of the channel that is perpendicular to the droplet length direction. The relation Ldissolv ed ðtÞ ¼ V dissolv ed ðtÞ=Dlz is valid because during dissolution, the droplet will always maintain the shape of the bubble caps at the frontal (cap EFG) and rear (cap ABC) ends (Fig. 15). Therefore to estimate the droplet length loss to dissolution, we only need to consider the reduction of the droplet’s rectangular core ACDF, which can be calculated as V dissolv ed ðtÞ ¼ Ldissolv ed Dlz . Dissolution only occurs in the regions where the droplet is in contact with the silicone oil, which are at the frontal and rear caps. In our derivation, we assume a constant overall mass transfer flux rate coefficient J (kg/m2 s), and multiplying this with the area of the frontal and rear caps (Acap , m2) will yield us our dissolution mass transfer rate (kg/s). Acap can be presented as:

Acap ¼ 2  arc length  lz

ð36Þ

where the arc length refers to the length of the arcs ABC and DEF. To calculate the arc length of a 2D droplet that makes a contact angle h with the wall (Fig. 15), we use the following relation:

Table 3 Fluid properties of CO2 and silicone oil used in the droplet dissolution simulation [49].

q ðkg=m3 Þ m ðcSÞ

kH ðbarÞ [52] Dab;L ðm2 =sÞ [53] r ðmN=mÞ [54] J ðkg=m2 sÞ [51] h ð Þ [51]

CO2

Silicone oil

1.98 7.5 199

818 1.0 n/a

11  109 17.4 0.00315 25

arc length ¼ r arc harc ¼ r cap ðp  2hÞ ¼

Dðp  2hÞ 2 cos h

ð37Þ

Using this relation in our Acap (Eq. (36)), we can express our dissolution mass transfer as:

  Dðp  2hÞ ðp  2hÞ _ ¼ JAcap ¼ J 2  m  lz ¼ JDlz 2 cos h cos h

ð38Þ

To get V dissolv ed ðtÞ, we simply divide the mass transfer rate by the CO2 density and multiply by the time of interest t (valid since we used a constant J):

V dissolv ed ðtÞ ¼

JDlz ðp  2hÞ t qCO2 cos h

ð39Þ

Therefore, the final expression for the length of the dissolving CO2 droplet is:

Ldroplet ðtÞ ¼ Linitial 

Jðp  2hÞ t qCO2 cos h

ð40Þ

We can also obtain the volume of the dissolving CO2 droplet using the V dissolv ed ðtÞ term (Eq. (39)) directly:

V droplet ðtÞ ¼ V initial 

JDlz ðp  2hÞ t qCO2 cos h

Fig. 14. T-junction microchannel geometry used to create the CO2 droplet.

ð41Þ

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

931

Table 4 Droplet length and mass transfer coefficient obtained using different mesh densities.

Droplet length (lm) J ðkg=m2 sÞ

36 cells

48 cells

60 cells

72 cells

3480 0.00246

3670 0.00309

3930 0.00341

3960 0.00343

Fig. 15. Schematic of droplet geometry used to derive length of dissolving CO2 droplet.

The simulations were carried out for a real flow time of t = 0.25 s. Results obtained for the dissolving droplet length and volume over time are plotted in Fig. 16. Overall, there is a good match between both mass transfer models and the analytical solution. The results for Bothe and Fleckenstein’s model overestimate the dissolution of the droplet, with differences of 2:87% and 4:88% with the analytical droplet length and volume respectively at t max . This is due to the slightly higher prediction of the overall mass transfer coefficient by Bothe and Fleckenstein’s model (0.00341 kg/m2 s) compared to the experimental measurement (0.00315 kg/m2 s), which is 8:25% higher. Nonetheless, given that Bothe and Fleckenstein’s model predicted the mass transfer rate using variables that are entirely independent from the experimental values (diffusion coefficient and concentration gradient), it corroborates that the model is reliable in predicting mass transfer performance. Furthermore, the concentration gradient in the model is calculated in light of the expulsion operation, and given that the accuracy of the mass transfer prediction is maintained, this also shows that the effects of the expulsion operation on the mass transfer model performance is minimal. In Fig. 16 it can also be observed that the constant mass transfer coefficient model slightly overestimates the droplet length but

under-approximates the droplet volume. This occurs because the frontal and rear cap areas for a dynamic droplet is marginally higher than those of a stationary droplet, since the caps are more curved during droplet motion. As a result the droplet length measured is slightly higher for the constant mass transfer model, since the caps form a larger protrusion, but will have more volume lost to dissolution since the cap area where mass transfer occurs is higher. Nevertheless, the difference for both the droplet length and volume at tmax is only around 1%, which indicates that both results are still highly correlated. This demonstrates the fact that the stationary droplet assumption for deriving the analytical solution is a valid approximation, and that the algorithm presented in Section 2.4 is accurate in estimating the interfacial area aVOF . To confirm the observation where the larger protrusion of the dynamic droplet results in a longer length and higher dissolution rate, a simulation for the dissolving CO2 droplet is also conducted for a wall contact angle of h ¼ 10 . The results in Fig. 17 show that the droplet with h ¼ 10 has a longer initial droplet length, but due to its larger frontal and rear cap areas, the droplet length eventually decreases to be less than the h ¼ 25 droplet due to the higher dissolution rate. The same observation can also be made for the droplet volume result, where the h ¼ 10 droplet volume decreases at a greater rate than the h ¼ 25 scenario. Fig. 18presents the simulation images of the dissolving CO2 droplet as it travels along the microchannel and also the ensuing species transport in silicone oil due to the dissolved CO2. The results in the bottom half of Fig. 18 show that the dissolved CO2 species only remain within the silicone oil regions and not within the CO2 gas region (represented by the red droplets in the top half of Fig. 18). This implies that the CO2 species are being tracked within the correct liquid phase, which demonstrates the effectiveness of the a-factor and expulsion operation adopted in Section 2.3. As the droplet dissolves, the dissolved CO2 species at the front cap form a protruding cone-like structure, while the dissolved CO2 at the rear cap diffuses away from the centre towards the walls of the microchannel. This observation is similar to the dissolved

Fig. 16. Length (left) and volume (right) of the dissolving CO2 droplet over time obtained using Bothe and Fleckenstein’s mass transfer model and constant mass transfer coefficient J, compared with the derived analytical solution.

932

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

Fig. 17. Length (left) and volume (right) of the dissolving CO2 droplet over time obtained using the constant mass transfer coefficient J for contact angles of h ¼ 10 and h ¼ 25 , compared with the derived analytical solution using h ¼ 25 .

Fig. 18. Dissolving CO2 droplet over time (top) and evolution of dissolved CO2 species in silicone oil (bottom).

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

species results by Jia and Zhang [35], and is due to the vortices in the vicinity of the front and rear caps that carry the dissolved species along the vortex flow direction. The current simulation only captures a single dissolving droplet, and it would be interesting in the future to study the species transport behaviour as a result of multiple droplets within the same microchannel. 4. Conclusion In this paper, the formulation of a coupled multiphasemulticomponent model that is applicable for micro-scale simulations has been described. The multiphase physics is captured using VOF, with smoothing operations applied to minimise the issue of spurious velocities. The species transport equation is adopted for multicomponent tracking, with the a-factor and expulsion operation methods implemented to ensure that the species are tracked within the right phase. Three mass transfer models are adopted: constant mass transfer coefficient, Bothe and Fleckenstein [29] and Haelssig et al.’s [30] models, with some modifications to the original models to improve their performance. The coupled multiphase-multicomponent model is validated using three validation cases to assess the model’s accuracy. The first case was a simple 1D diffusion problem, where the model yielded matching results with the analytical solution for different flow times and diffusion coefficients. This test case establishes the accuracy of the species transport solver in the model. The next case was the simulation of CO2 vapour flowing horizontally over a smooth stationary pentane liquid. Mass transfer occurs as the pentane evaporates into the CO2 phase, and the Sherwood number is used to measure the mass transfer performance. Both the Bothe and Fleckenstein and Haelssig et al.’s models were tested and were found to obtain Sherwood numbers that are in good agreement with empirical correlations. Lastly, a dissolving CO2 droplet within silicone oil in a 2D microchannel was simulated. The analytical solution for a dissolving droplet length and volume over time was derived, and compared with simulations results of the Bothe and Fleckenstein and constant mass transfer coefficient models. Overall, both models yielded close results with the analytical solution. The mass transfer coefficient predicted using Bothe and Fleckenstein’s model was within 10% of the experimental value, which demonstrated the reliability of the model. Simulation images are presented for the transport of dissolved CO2 species as the droplet travels along the microchannel, where it was observed that the species were tracked within the correct phase, thus establishing the effectiveness of the a-factor and expulsion operations. The coupled multiphase-multicomponent model developed in this work can be utilised to study a variety of microscale applications: dissolution of droplets in microchannels, transport of drugs in blood vessels that entails mass transfer and simulations of carbon dioxide enhanced oil recovery at the pore-scale. The simulation of such problems using the model is able to provide insights and understanding of flow physics which can aid in the design of devices or determining the right operating conditions for any experimental studies, in a cost-effective manner. Conflict of interest The authors declare that there is no conflict of interest. Acknowledgments The authors would like to express their gratitude for the support provided by the Australian Government Research Training Program Scholarship.

933

References [1] N. Kantarci, F. Borak, K.O. Ulgen, Bubble column reactors, Process Biochem. 40 (7) (2005) 2263–2283. [2] Q. Xu, M. Hashimoto, T.T. Dang, T. Hoare, D.S. Kohane, G.M. Whitesides, R. Langer, D.G. Anderson, Preparation of monodisperse biodegradable polymer microparticles using a microfluidic flow-focusing device for controlled drug delivery, Small 5 (13) (2009) 1575–1581. [3] M. Sussman, E. Fatemi, An efficient, interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow, SIAM J. Sci. Comput. 20 (4) (1999) 1165–1191. [4] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159. [5] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of multiphase flow, J. Comput. Phys. 169 (2) (2001) 708–759. [6] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiphase phases and components, Phys. Rev. E 47 (3) (1993) 1815–1819. [7] X. Shan, H. Chen, Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E 49 (4) (1994) 2941–2948. [8] C. Yang, Z.-S. Mao, Numerical simulation of interphase mass transfer with the level set approach, Chem. Eng. Sci. 60 (10) (2005) 2643–2660. [9] A. Koynov, J.G. Khinast, G. Tryggvason, Mass transfer and chemical reactions in bubble swarms with dynamic interfaces, AIChE J. 51 (10) (2005) 2786–2800. [10] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1) (1981) 201–225. [11] D. Kothe, W. Rider, S. Mosso, J. Brock, J. Hochstein, Volume tracking of interfaces having surface tension in two and three dimensions, AIAA Pap. 96 (0859) (1996) 1–18. [12] D.L. Youngs, Time-dependent multi-material flow with large fluid distortion, Numer. Methods Fluid Dyn. 24 (1982) 273–285. [13] D. Bothe, M. Koebe, K. Wielage, H.-J. Warnecke. VOF-simulations of mass transfer from single bubbles and bubble chains rising in aqueous solutions, in: ASME/JSME 2003 4th Joint Fluids Summer Engineering Conference, American Society of Mechanical Engineers, 2003. [14] M.R. Davidson, M. Rudman, Volume-of-fluid calculation of heat or mass transfer across deforming interfaces in two-fluid flow, Numer. Heat Transf.: Part B: Fundam. 41 (3–4) (2002) 291–308. [15] M. Kroger, A. Alke, D. Bothe, H. Warnecke, A VOF-based approach for the simulation of reactive mass transfer from rising bubbles, Fortschritt BerichteVDI Reihe 3, Verfahrenstechnik 883 (2007) 290. [16] M. Magnini, B. Pulvirenti, J.R. Thome, Numerical investigation of hydrodynamics and heat transfer of elongated bubbles during flow boiling in a microchannel, Int. J. Heat Mass Transf. 59 (2013) 451–471. [17] M. Magnini, B. Pulvirenti, J.R. Thome, Numerical investigation of the influence of leading and sequential bubbles on slug flow boiling within a microchannel, Int. J. Therm. Sci. 71 (2013) 36–52. [18] R. Banerjee, A numerical study of combined heat and mass transfer in an inclined channel using the VOF multiphase model, Numer. Heat Transf., Part A: Appl. 52 (2) (2007) 163–183. [19] R. Zhuan, W. Wang, Simulation on nucleate boiling in micro-channel, Int. J. Heat Mass Transf. 53 (1–3) (2010) 502–512. [20] R. Zhuan, W. Wang, Flow pattern of boiling in micro-channel by numerical simulation, Int. J. Heat Mass Transf. 55 (5–6) (2012) 1741–1753. [21] S. Hardt, F. Wondra, Evaporation model for interfacial flows based on a continuum-field representation of the source terms, J. Comput. Phys. 227 (11) (2008) 5871–5895. [22] D.A. Hoang, V. van Steijn, L.M. Portela, M.T. Kreutzer, C.R. Kleijn, Benchmark numerical simulations of segmented two-phase flows in microchannels using the volume of fluid method, Comput. Fluids 86 (2013) 28–36. [23] P. Horgue, M. Prat, M. Quintard, A penalization technique applied to the ‘‘Volume-Of-Fluid” method: wettability condition on immersed boundaries, Comput. Fluids 100 (2014) 255–266. [24] A.Q. Raeini, M.J. Blunt, B. Bijeljic, Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method, J. Comput. Phys. 231 (17) (2012) 5653–5668. [25] G.Y. Soh, G.H. Yeoh, V. Timchenko, Improved volume-of-fluid (VOF) model for predictions of velocity fields and droplet lengths in microchannels, Flow Meas. Instrum. 51 (2016) 105–115. [26] A.Q. Raeini, B. Bijeljic, M.J. Blunt, Numerical modelling of sub-pore scale events in two-phase flow through porous media, Transp. Porous Media 101 (2) (2014) 191–213. [27] A.Q. Raeini, M.J. Blunt, B. Bijeljic, Direct simulations of two-phase flow on micro-CT images of porous media and upscaling of pore-scale forces, Adv. Water Resour. 74 (2014) 116–126. [28] G.Y. Soh, G.H. Yeoh, V. Timchenko, Numerical investigation on the velocity fields during droplet formation in a microfluidic T-junction, Chem. Eng. Sci. 139 (2016) 99–108. [29] D. Bothe, S. Fleckenstein, A volume-of-Fluid-based method for mass transfer processes at fluid particles, Chem. Eng. Sci. 101 (2013) 283–302. [30] J.B. Haelssig, A.Y. Tremblay, J. Thibault, S.G. Etemad, Direct numerical simulation of interphase heat and mass transfer in multicomponent vapour– liquid flows, Int. J. Heat Mass Transf. 53 (19–20) (2010) 3947–3960. [31] M.M. Francois, N.N. Carlson, The global embedded interface formulation for interfacial mass transfer within a volume tracking framework, Comput. Fluids 87 (2013) 102–114.

934

G.Y. Soh et al. / International Journal of Heat and Mass Transfer 113 (2017) 922–934

[32] Y. Haroun, D. Legendre, L. Raynal, Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film, Chem. Eng. Sci. 65 (10) (2010) 2896–2909. [33] H. Marschall, K. Hinterberger, C. Schüler, F. Habla, O. Hinrichsen, Numerical simulation of species transfer across fluid interfaces in free-surface flows using OpenFOAM, Chem. Eng. Sci. 78 (2012) 111–127. [34] A. Onea, M. Woerner, D.G. Cacuci, A qualitative computational study of mass transfer in upward bubble train flow through square and rectangular minichannels, Chem. Eng. Sci. 64 (7) (2009) 1416–1435. [35] H. Jia, P. Zhang, Investigation of the Taylor bubble under the effect of dissolution in microchannel, Chem. Eng. J. 285 (2016) 252–263. [36] J.B. Haelssig, A.Y. Tremblay, J. Thibault, A new hybrid membrane separation process for enhanced ethanol recovery: process description and numerical studies, Chem. Eng. Sci. 68 (1) (2012) 492–505. [37] Y. Haroun, L. Raynal, D. Legendre, Mass transfer and liquid hold-up determination in structured packing by CFD, Chem. Eng. Sci. 75 (2012) 342– 348. [38] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335–354. [39] Y. Renardy, M. Renardy, PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method, J. Comput. Phys. 183 (2) (2002) 400–421. [40] M.J. Nieves-Remacha, L. Yang, K.F. Jensen, OpenFOAM computational fluid dynamic simulations of two-phase flow and mass transfer in an advanced-flow reactor, Ind. Eng. Chem. Res. 54 (26) (2015) 6649–6659. [41] Z. Guo, D. Fletcher, B. Haynes, A review of computational modelling of flow boiling in microchannels, J. Comput. Multiph. Flows 6 (2) (2014) 79–110. [42] J. Schlottke, B. Weigand, Direct numerical simulation of evaporating droplets, J. Comput. Phys. 227 (10) (2008) 5215–5237.

[43] G.Y. Soh, G.H. Yeoh, V. Timchenko, An algorithm to calculate interfacial area for multiphase mass transfer through the volume-of-fluid method, Int. J. Heat Mass Transf. 100 (2016) 573–581. [44] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (2) (1998) 112–152. [45] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operatorsplitting, J. Comput. Phys. 62 (1) (1986) 40–65. [46] H. Rusche, Computational Fluid Dynamics of Dispersed Two-phase Flows at High Phase Fractions, Imperial College of Science, Technology and Medicine, London, England, 2002. [47] B. Van Leer, Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow, J. Comput. Phys. 23 (3) (1977) 263–275. [48] J. Crank, The Mathematics of Diffusion, Oxford University Press, 1979. [49] R.H. Perry, D.W. Green, Perry’s Chemical Engineers’ Handbook, McGraw-Hill Professional, 1999. [50] H. Cheng, M.E. Pozo de Fernández, J.A. Zollweg, W.B. Streett, Vapor-liquid equilibrium in the system carbon dioxide+ n-pentane from 252 to 458 K at pressures to 10 MPa, J. Chem. Eng. Data 34 (3) (1989) 319–323. [51] M. Sauzade, T. Cubaud, Initial microfluidic dissolution regime of CO2 bubbles in viscous oils, Phys. Rev. E 88 (5) (2013) 051001. [52] M. Kohrt, I. Ausner, G. Wozny, J.-U. Repke, Texture influence on liquid-side mass transfer, Chem. Eng. Res. Des. 89 (8) (2011) 1405–1413. [53] A.C.M. Kuo, Poly(dimethylsiloxane), in: J.E. Mark (Ed.), Polymer Data Handbook, Oxford University Press, Oxford, 1999, pp. 411–435. [54] R. Sun, Carbon Dioxide Multiphase Flows in Microfluidic Devices, The Graduate School, Stony Brook University, Stony Brook, NY, 2010.