# A consistent methodology to transport a passive scalar with the geometric Volume-of-Fluid method isoAdvector

Alexis Tourbier, Lionel Gamet, Philippe Béard, Typhène Michel, Joelle Aubin, Hrvoje Jasak

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alexis Tourbier, Lionel Gamet, Philippe Béard, Typhène Michel, Joelle Aubin, et al. A consistent methodology to transport a passive scalar with the geometric Volume-of-Fluid method isoAdvector. Journal of Computational Physics, 2024, 513, pp.113198. 10.1016/j.jcp.2024.113198. hal-04642305

HAL Id: hal-04642305

https://ifp.hal.science/hal-04642305v1

Submitted on 9 Jul 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ELSEVIER

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

# A consistent methodology to transport a passive scalar with the geometric Volume-of-Fluid method isoAdvector

Alexis Tourbier<sup>a,b,*</sup>, Lionel Gamet<sup>a</sup>, Philippe Béard<sup>a</sup>, Typhène Michel<sup>a</sup>,
Joelle Aubin<sup>b</sup>, Hrvoje Jasak<sup>c</sup>

## ARTICLE INFO

**Keywords:**

Transport of concentration

Face interpolation

Volume-of-Fluid

isoAdvector

OpenFOAM

## ABSTRACT

This paper presents a new methodology allowing the discretisation of phase specific transport equations within the Volume-of-Fluid interface capturing method framework. The method uses a sharp interface algorithm to compute the transport of species concentration. The interface reconstruction and advection is provided by the geometric advection scheme *isoAdvector* [1] implemented in the OpenFOAM® library. When discretising the transport equation for species, lack of consistency with the free surface advection scheme can lead to numerical errors, causing conservation or boundedness issues. This work addresses the issue of consistency in convective transport of species and is divided in two parts. First, a new interpolation procedure is used to compute face values from cell-centered values. Then, the diffusive operator of the transport equation is corrected. Finally, a set of test cases are presented to validate the transport equation's consistency with the free surface advection. Species transport across the interface is not part of the scope of this article, however, this methodology can further be used to study mass transfer at the gas-liquid interface using additional mass source terms that are not discussed here.

## 1. Introduction

Multiphase catalytic processes are prevalent in diverse industries including commodity and fine chemicals, as well as pharmaceuticals. There is a continual innovation and development for such processes in order to improve reaction performance but also to comply with evolving environmental and safety regulations [2]. Catalyst testing is often performed in batch or fed-batch reactors with volumes of several litres. However, temperature and mixing may be difficult to control due to the relatively large volume of the reactor. The reactors are also often operated at high pressure, which requires a further depressurization step. This step might be harmful for the products in case of secondary reactions [3]. One alternative to batch reactors is the choice of continuous flow capillary microreactors. These types of reactors require much smaller volumes of reagents than batch reactors and temperature control is easier. For these reasons, microfluidics – as a viable alternative to conventional reactors for catalytic tests – have attracted considerable attention over the last years [4–9]. To ensure micro capillaries can perform at least as well as the agitated batch and fed-batch reactors, proof is required to demonstrate adequate reaction performance.

<footnote><sup>a</sup> IFP Énergies nouvelles, Rond-point de l'échangeur de Solaize, 69360 Solaize, France</footnote>

<footnote><sup>b</sup> Université de Toulouse, INP, LGC (Laboratoire de Génie Chimique), Toulouse, 4 allée Emile Monso BP-84234, 31432 Toulouse Cedex 4, France</footnote>

<footnote><sup>c</sup> University of Cambridge, Department of Physics, The Cavendish Laboratory, 19 J.J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom</footnote>

<footnote>* Corresponding author at: IFP Énergies nouvelles, Rond-point de l'échangeur de Solaize, 69360 Solaize, France.
E-mail address: [alexis.tourbier@ifpen.fr](mailto:alexis.tourbier@ifpen.fr) (A. Tourbier).</footnote>
Computational Fluid Dynamics (CFD) is now being considered as a design tool and used in the aim of process optimization and intensification. Nevertheless, simulating reactive two-phase flows requires numerical tools rigorously tested and validated. Indeed, this type of flow requires the interface between the fluids to be resolved and also an additional transport equation for the concentration to track the species in the domain. Ultimately, this additional transport equation can have an influence on momentum conservation, for example due to volume changes [10]. In the present work, the Volume of Fluid (VOF) method is used. This method was first detailed in the work of Hirt and Nichols [11]. The coupling of the hydrodynamics solver with species transfer has been addressed in several studies, depending on the method used to solve the free-surface flow. For example, front tracking methods [12–16], arbitrary Lagrangian-Eulerian (ALE) [17–19], Level-Set (LS) [20,21] and also Volume of Fluid (VOF) [22–28].

In the OpenFOAM® library [29], two different VOF methods exist, the so-called algebraic and geometric VOF. While the algebraic VOF solves an equation for the transport of an indicator function [30], the geometric VOF uses a geometric representation of the interface and its intersection with cell faces to compute the phase specific fluxes [1,31]. Predicting the surface tension force is the main difficulty when using the VOF methods. The formulation of the surface tension force used in most solvers was described as the Continuous Surface Force (CSF) developed by Brackbill et al. [32]. However, errors in the estimation of interface curvature are responsible for the generation of non-physical (or spurious) velocities [33]. The geometric VOF solvers were designed to reduce spurious velocities. An additional interface reconstruction step allows an explicit location of the interface in the cell and thus a better prediction of the variable interpolations.

When dealing with species transport and chemical reactions, an advection-diffusion equation for each species must be added to the problem. To describe the transport of these species, two methods are possible. Either a unique concentration field is used to describe the species in the whole domain (this is the method proposed by Haroun et al. [24] and applied in [25,34] as the Continuous Species Transfer (CST)) or two concentration fields can be created, one for each phase. The corresponding concentration field should therefore be null in the other phase. This two-field approach was successfully applied in the work of Bothe and Fleckenstein [26].

In order for the concentration fields to be consistent with the interface, one must use the same interpolation weights to compute the face values of the field. However, in the geometric VOF method, the volume fraction field $\alpha$ is not explicitly interpolated onto the cell faces. Our work focuses on a new methodology being able to interpolate a concentration field onto the cell faces, in order to stay consistent with the interface. The most important point in a two-field approach is to avoid artificial mass transfer due to badly computed face fluxes.

A geometric VOF method is also available in the free software Basilisk (more information on the software can be found in [35,36]) and several authors used it to study scalar-transport related cases. Farsoiya and coworkers [37] studied mass transfer in the wake of a gas bubble in a turbulent flow, López-Herrera et al. [38] worked on electrohydrodynamics, and multicomponent droplet evaporation was simulated by the group of Cipriano [39]. All of these cases rely on the transport of a scalar that must be advected consistently with the VOF method provided by Basilisk. To do so, the Bell-Colella-Glaz second order upwind scheme [40] is used to compute the face value of the scalar from the cell-centered value. However, this scheme was not implemented in OpenFOAM® and the comparison of its performances with classical interpolation schemes available in OpenFOAM® is not part of the scope of the present study.

Consistency with the diffusion operator is also critical. The transport of the free surface is an advection-only equation and there is no diffusion of the interface. To be able to stay consistent with the interface, diffusion can only occur in the bulk of the phases and the diffusion coefficient must be corrected to stay consistent with the interface. If consistency with both advection and diffusion operator is not enforced, artificial mass transfer might arise, causing numerous numerical errors. These numerical errors will lead to further numerical issues when adding physical mass transfer on top of the transport equation, which is why consistency with the VOF equation is so important when dealing with multiple transport equations.

In this work, a new methodology for the interpolation of a scalar field from the cell centres to the cell faces is presented. These face values are used in the divergence operator to compute the flux of the scalar quantity. A correction of the diffusion coefficient in order to avoid the diffusive fluxes through the faces of cells containing the interface between phases is also detailed. Finally, test-cases evaluating the efficiency of the method are discussed. To our knowledge, no work in the literature focuses specifically on consistency of the advection and diffusion operators. The focus of this work is to give an example on how to start with the transport equation in order to ensure consistency with the interface. The source code and the test cases presented in this document are available in the git repository [41].

## 2. Modelling

### 2.1. The Volume of Fluid method

The VOF methods use an indicator function to represent the discontinuity between the two phases. This step function is defined as:

$$H(x, t) = \begin{cases} 0 & \text{in the gaseous phase} \\ 1 & \text{in the liquid phase} \end{cases}$$

The volume fraction $\alpha$ of the liquid phase for a cell $i$ can be derived as follows:

$$\alpha_i(t) = \frac{1}{V_i} \int_{\Omega_i} H(x, t) dV \qquad (1)$$
where $\Omega_i$ represents the cell $i$ and $V_i = \int_{\Omega_i} dV$ its volume. The volume fraction of the gaseous phase is implicitly defined as $1 - \alpha$. The resulting evolution equation for $H$ is obtained from the continuity equation:

$$ \frac{\partial H}{\partial t} + \nabla \cdot (\mathbf{U} H) = 0 \qquad (2) $$

with the condition that $\nabla \cdot U = 0$. Two categories of VOF methods exist. The algebraic VOF method, uses equation (2) to implicitly capture the interface position. In OpenFOAM®, it is used along with the Multidimensional Universal Limiter for Explicit Solution (MULES) to provide a sharper interface. The solver is named interFoam and was historically the first implemented in OpenFOAM®. More details about its formulation and performance can be found in the work of Deshpande and coworkers [30] and in the work of Bilger et al. [42].

The issue with algebraic VOF methods is the smearing of the interface. Geometric methods used in combination with a reconstruction algorithm provide a sharper interface representation [33]. For surface tension dominated flows (i.e. low capillary numbers flows), the prediction of the surface tension force offered by the interFoam solver is not sufficient and generates spurious currents [43,44]. Therefore, the geometric VOF method isoAdvector has been chosen in this work and will be presented in the following section.

## 2.2. The geometric VOF algorithm isoAdvector

Recently, a geometric VOF solver has been developed in OpenFOAM®. It differs from the algebraic VOF solver in the explicit reconstruction of the interface from the volume fraction field. With a sharper interface calculation, the curvature prediction, and therefore the computation of surface tension force, can potentially be improved. The drawback of geometric methods is the difficulty of implementation compared to algebraic method. isoAdvector is a numerical method developed for the advection of a passively advected surface. It was implemented in OpenFOAM® to be used in multiphase solvers to advect the interface between two incompressible, isothermal fluids. Full details on the algorithm are available in [1]. The present section will be used to give some background on the method, which will be necessary to explain the contribution of this work.

isoAdvector introduces two novel features. First, the reconstruction step uses the concept of isosurfaces to locate the interface in each cell. Then, the intersection between the interface and the face (the face-interface intersection line, FIIL) is updated with a vertex-interpolated velocity to allow the computation of the phase specific face fluxes. The starting point is the phase indicator function, which is computed in each cell based on the phase densities $\rho_1$ and $\rho_2$:

$$ H(\mathbf{x}, t) = \frac{\rho(\mathbf{x}, t) - \rho_2}{\rho_1 - \rho_2} \qquad (3) $$

The volume fraction is defined exactly as equation (1). After some manipulation and integration in time, the new time step volume fraction of cell $i$ is given by:

$$ \alpha_i(t + \Delta t) = \alpha_i(t) - \frac{1}{V_i} \sum_f s_{if} \Delta V_f(t, \Delta t) \qquad (4) $$

In equation (4), $s_{if}$ on face $f$ equals +1 or -1 to ensure that $s_{if} dS$ points out of cell $i$. $\Delta V_f$ is the volume of fluid transported across face $f$ during a time step. It is computed as:

$$ \Delta V_f(t, \Delta t) = \int_{t}^{t+\Delta t} \int_{A_f} H(\mathbf{x}, \tau) \mathbf{U}(\mathbf{x}, \tau) \cdot d\mathbf{S} \, d\tau \qquad (5) $$

with $A_f$ the surface of face $f$. To guarantee mass conservation, this “volumetric face flux” is used to update both cells sharing face $f$.

isoAdvector is originally distributed with a reconstruction algorithm called isoAlpha. This algorithm is based on the idea of “isosurfaces”, i.e. surfaces of same volume fraction. The isosurface 0.5 is often used to visualize the concept of interface between two phases. However, this isosurface does not necessarily cut a cell into two subsets of volumes corresponding to the volume fraction of the cell. Therefore, isoAlpha uses the volume fraction field information to find which isovalue best fits the liquid distribution in each cell. This description of isovalues implies that the interface is not connected between cells.

The computation of fluxes will now be explained. All the details are inspired from the papers of the authors [1,45]. The basic concepts are presented to understand the importance of the developments made in the current study. Fig. 1 illustrates a simple case where an interface translates in a cubic cell. To translate the interface, the cell-centered velocity field is interpolated onto the interface centroid. The previous time step velocity is used and considered constant.

Fig. 2 and 3 illustrate the time evolution of the interface. The initial area immersed in the liquid phase is given in Fig. 2 by $A_{im}(t)$. During a time step, the portion of the face that will be swept by the interface is given by the area of the polygon $A_{im}(\tau) - A_{im}(t)$. Using the letters defined in Fig. 3 together with the approximation of constant velocity during a time step, the position of the intersection at an intermediate time $\tau$ can be expressed as:

$$ B = A + (C - A) \frac{\tau - t}{(t + \Delta t) - t} \qquad E = D + (F - D) \frac{\tau - t}{(t + \Delta t) - t} \qquad (6) $$
<figure><img src="image_6.png" /><figcaption>Fig. 1. Advection of a free surface with isoAdvector during a time step $\Delta t$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 2. Liquid immersed part of the top face ($A_{im}$) during a time step $\Delta t$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 3. Illustration for the computation of the liquid immersed area.</figcaption></figure>

It is possible to obtain the area from Fig. 3 with the relation:

$$A_{im}(\tau) - A_{im}(t) = \frac{1}{2} |AE \times DB| \qquad (7)$$

After manipulation, it appears that the time dependent immersed area is a second-order polynomial:

$$A_{im}(\tau) - A_{im}(t) = a\tau^2 + b\tau + c \qquad (8)$$

Letters a, b and c are obtained after combining equations (6) and (7) and identifying the terms according to their order in the polynomial. When this second-degree polynomial is integrated over time, it gives a third-degree polynomial, which is quite easy to fit. For this reason, isoAdvector is computationally efficient compared with other geometric VOF solvers. The procedure is applicable to any grid topology, including polyhedral meshes. The volumetric face flux $\Delta V_f(t, \Delta t)$ is used to define $\alpha_f$, which is the part of the face that is immersed in the liquid. The definition of the liquid immersed part of face f is:

$$\alpha_f = \frac{\Delta V_f}{F_f \Delta t} \qquad (9)$$
The volumetric face flux $F_f$ results from the velocity-pressure solver. `isoAdvector` uses a converged velocity from the old time step. Species transport is then performed thanks to the face-area-averaged volume fraction after the interface advection. This face value $\alpha_f$ is the cornerstone of our developments as it will be used to enforce consistency with `isoAdvector`.

## 2.3. The piecewise linear interface construction

The original algorithm proposed by [1] in combination with `isoAdvector` to reconstruct the interface from the volume fraction field was called `isoAlpha`. The algorithm determines the best isovalue for the isosurface to match the actual amount of liquid inside the cell. The drawback of the method is the connectivity of the interface. Indeed, the interface is not necessarily connected between cells. Note that this is also the case for some other reconstruction schemes such as the Piecewise Linear Interface Construction (PLIC) (see [46] for a description of the method). A new implementation of PLIC is now available in combination with `isoAdvector` (see [31] for more details). The PLIC method is interesting because it allows a very precise curvature prediction. This new implementation is available as a separate git repository [47] and the solver is called `interFlow`. Benchmark tests have shown the robustness of `interFlow` for the resolution of surface tension dominated flows [43], reducing the intensity of parasitic currents by two orders of magnitude compared with MULES or `isoAdvector` `isoAlpha`. Our choice of using `interFlow` was thus motivated by its good performance in capillary flows. The remainder of this document considers only `interFlow` as the VOF solver. Nevertheless, other VOF solvers could be used provided that adequate face values are returned by the VOF solver.

## 2.4. The two-field approach

In chemical engineering, concentration or mass fractions are often used to describe chemical species. Concentration is useful when comparing multiple reactions to obtain reaction rates or yields. To represent the species distribution in the domain, a transport equation is required to follow the species concentrations in every control volume. As was mentioned in the introduction, there are two possibilities to describe the concentration field. The single-field approach has been used by Haroun and coworkers [24] to study chemical mass transfer at a gas-liquid interface. Marschall et al. [25] then developed the Continuous Species Transfer using a different computation for the diffusion coefficient in the interface region. Maes and Soulaine [34] also have successfully included volume change and derived the Compressive-CST (C-CST), which improves consistency with MULES of `interFoam`. To avoid the difficulty of the discontinuous concentration field at the interface, the two-field approach is interesting as it uses two separate fields to describe the concentration in the domain. It also avoids the complicated computation of the phase properties in the cells where both phases exist. In the two-field approach, the interface concentration is directly accessible through the linear extrapolation of the phase concentration to the interface. The two-field approach is also more consistent with a geometric VOF method, see [34]. Bothe and Fleckenstein [26] used the two-field approach for the transport of concentration. Later, the same authors included volume change to multi-component mass transfer [10]. Both single field and two-field approaches can yield good results. However, the most important thing to keep in mind is artificial mass transfer. Artificial mass transfer occurs when the concentration of a given species becomes non-physical in a phase. It is artificial because the fluxes responsible for the creation of non-zero concentration come from the discretisation of the transport equation.

The difficulty with modelling two-phase flows is making the transport of concentration with the interface advection resulting from the VOF solver consistent. Consider a simple case of a gas bubble containing a dilute species that is advected in a channel filled with liquid. The dilute species in the gas is not soluble in the liquid. Therefore, its concentration should remain zero everywhere in the liquid phase. If the concentration of the species becomes non-zero in the liquid phase, this is because there is a flux of species entering a cell filled with liquid.

The averaging tools needed to derive a transport equation consistently with a geometric VOF method such as `isoAdvector` are as follows. Consider a gas-liquid two-phase flow; in a two-field approach, the concentration of a species $C_k$ is split into a liquid concentration field and its gaseous counterpart:

$$C_{k,1}(\mathbf{x}, t) = \begin{cases} C_k(\mathbf{x}, t) & \text{in phase 1} \\ 0 & \text{in phase 2} \end{cases} \qquad (10)$$

$$C_{k,2}(\mathbf{x}, t) = \begin{cases} 0 & \text{in phase 1} \\ C_k(\mathbf{x}, t) & \text{in phase 2} \end{cases} \qquad (11)$$

Two averaged quantities can be expressed, namely the volume average and the phase average. The volume average is the quantity that is transported on the computational grid in a Finite Volume Method (FVM) framework. It is the concentration averaged over the volume $V$ of a cell:

$$\langle C_{k,1} \rangle_V = \frac{1}{V} \int_{V_1} C_{k,1} \, dV \qquad (12)$$

Note that since $C_{k,1}$ must be zero in the gas phase, integrating over $V_1$, which is the volume of liquid inside the cell, or over $V$, which is the total volume of the cell, are equivalent. A similar definition holds for phase 2. The phase average can be defined as:

$$\langle C_{k,1} \rangle_{V_1} = \frac{1}{V_1} \int_{V_1} C_{k,1} \, dV \qquad (13)$$
<figure><img src="image_9.png" /><figcaption>Fig. 4. Schematic of a one-dimensional channel used to test the advection performances of interpolation schemes.</figcaption></figure>

The phase average represents how the concentration is distributed within the phase. The two averages can be linked by the volume fraction:

$$ \langle C_{k,1} \rangle_V = \frac{1}{V_1} \int_{V_1} C_{k,1} \, dV = \frac{V_1}{V} \frac{1}{V_1} \int_{V_1} C_{k,1} \, dV = \alpha \langle C_{k,1} \rangle_{V_1} \quad (14) $$

Equation (14) is very important since both averaged quantities will be used in the following, depending on the situation.

## 2.5. Face interpolation

In this section, the basic interpolation tools used in section 3 are presented. Consider the transport equation for a scalar quantity $\phi$ with a diffusivity $D$ and a volumetric source term $S_v$:

$$ \frac{\partial \phi}{\partial t} + \nabla \cdot (U\phi) = \nabla \cdot (D\nabla\phi) + S_v \quad (15) $$

In FVM, this equation is integrated over a control volume and time. The divergence theorem is then used to convert the convection and viscous volume integrals into surface integrals on the contour of the control volume of unit normal $n$:

$$ \int_V \int_{\Delta t} \frac{\partial \phi}{\partial t} dt \, dV + \int_{\Delta t} \int_A U \cdot n \, \phi dt \, dS = \int_{\Delta t} \int_A D n \cdot \nabla \phi dt \, dS + \int_V \int_{\Delta t} S_v dt \, dV \quad (16) $$

The mid-point assumption stating that the integral is approximately equal to the mean value of the face variable multiplied by the surface of the face is used. This product can be approximated by the value at the centre of the face $f$ multiplied by the face area:

$$ \int_A \phi dA \approx \phi_f A_f \quad (17) $$

In the cell-centered OpenFOAM® FVM framework, the only known value of $\phi$ is at the cell centres. Therefore, a surface interpolation scheme is needed to provide the face value $\phi_f$ from the cell centered values of the two cells sharing the face. In most numerical codes, a weight $w$ is computed on each face depending on the selected scheme. This weight is then used to compute the face value $\phi_f$ with respect to the two cell centered values $\phi_p$ and $\phi_N$:

$$ \phi_f = \phi_N + w(\phi_p - \phi_N) \quad (18) $$

In the following of the text, different linear interpolation schemes are used: first order upwind differencing (UD), second order central differencing (CD) and Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL), which is also known as van Leer's scheme [48–52]. The interest of MUSCL is the bounding option that limits over- and under-shoots.

## 2.6. The importance of consistency for the transport of a passive scalar

Numerical diffusion is a well-known issue when trying to advect a sharp profile. This issue has motivated significant numerical developments, including the VOF method. When transporting a concentration field in two-phase flow, the same issue occurs for the concentration field. Indeed, the preservation of the species interface discontinuity imposed by thermodynamics requires the same kind of treatment as free surface flows. In an analogous manner, this “chemical interface”, where the fluid properties and concentrations undergo a jump, must coincide with the hydrodynamic gas-liquid interface resulting from the VOF solver. If it does not coincide, then artificial mass transfer has occurred because the concentration value is not what it should be. This issue can be seen on a simple test case of a bubble of gas flowing in a channel of liquid. Consider a dilute species inside the liquid phase that does not transfer through the interface and has uniform concentration. Since there is no concentration gradient in the liquid phase, the system corresponds to a purely advective case and the concentration should follow the interface. The set-up is illustrated in Fig. 4; the liquid is water and the gas air. To simulate this test-case, a 1D channel of 0.01 m is discretised with 1000 1D cells. A gas bubble is initialized and in the liquid the species has a concentration of $1 \text{ kg m}^{-3}$. The CFL number is set to 0.1 with an adaptive time step, the simulations last 0.7 s. The objective is to verify if the advection of the concentration field is consistent with the advection
<figure><img src="image_10.png" /><figcaption>Fig. 5. Comparison between the volume fraction field profile and nondimensionalized concentration field profiles obtained with different interpolation schemes. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

of the interface. An OpenFOAM® function object called scalarTransport is added in the controlDict to add a passive scalar transported by the two-phase flow. The diffusivity of this scalar is set to zero and several interpolation schemes used in the divergence term are compared. In Fig. 5, the concentration profiles obtained with these different interpolation schemes after 0.4 s are shown. The use of a first order UD scheme to interpolate the concentration onto the cell faces leads to an over-diffused concentration profile. This numerical diffusion spreads the initially sharp concentration profile, making it non-zero in the gas phase. Fig. 5 shows that the UD scheme in red fails to prevent artificial mass transfer. The use of higher order interpolation schemes to preserve the sharp profile has been proposed. However, second order schemes tend to be oscillatory, which might be an issue for bounded variables. In Fig. 5, the CD scheme, displayed as a blue curve, shows over- and under-shoots of the concentration profile. They result from the inaccurate face concentration prediction of the CD interpolation scheme. These over- and under-shoots are damaging for the simulation as they break thermodynamics of saturation or provide negative concentrations. The MUSCL scheme is interesting as it avoids over- and under-shoots. Furthermore, the profile obtained with the MUSCL scheme is sharper than with the UD scheme, thus, reducing artificial mass transfer. Nevertheless, none of these three schemes give fully satisfactory results, i.e. a bounded concentration profile, without artificial mass transfer. This work aims at presenting a new methodology that was successfully applied to the transport of concentration consistently with the geometric VOF method using isoAdvector.

## 3. Sharp interface consistent scalar transport equation

In section 2.6, it has been shown that discretising the transport equation consistently with the interface advection is crucial to avoid artificial mass transfer. In this section, a new methodology to discretise the transport of concentration will be described to overcome the issue of consistency with the advection of the free surface. Concentration transport is the first step to build a coupled reactive two-phase flow solver. If the transport step is not as accurate as possible, the prediction of reaction will inevitably be impacted.

### 3.1. General expression of a transport equation for concentration

The starting point is the following transport equation for the concentration of species k:

$$ \frac{\partial C_k}{\partial t} + \nabla \cdot (U C_k) = \nabla \cdot (D_k \nabla C_k) + S_v \quad (19) $$

This concentration field is split into two fields: two transport equations exist, one for each phase. Note that the phase concentrations are defined over the whole domain, but their values must be zero in the “other” phase.

### 3.2. Discretisation of the advective part of the transport equation

The first step is the discretisation of the advection operator. Only the advective transport of the concentration of species k in phase 1 is considered; as the same equation can be applied to the other phase:

$$ \frac{\partial C_{k,1}}{\partial t} + \nabla \cdot (U C_{k,1}) = 0 \quad (20) $$

In equation (20), the volume average concentration, as established by equation (12), is used after integration over the volume $V_i$ of cell i and a time step $\Delta t$:
<figure><img src="image_11.png" /><figcaption>Fig. 6. Schematic diagram of two cells cut by the interface with a volume fraction of 0.5.</figcaption></figure>

$$ \langle C_{k,1} \rangle_{V_i}^{n+1} = \langle C_{k,1} \rangle_{V_i}^{n} - \frac{\Delta f}{V_i} \times \sum_{f} (n \cdot U C_{k,1}^{n} A)_{f} \quad (21) $$

The main issue is to determine $(C_{k,1})_f$ in order to respect the consistency with isoAdvector. In the first implementation, relation (18) was applied to the VOF equation (2) to calculate the weights that would have produced the fluxes. The linear weight, $w_\alpha$, can be computed as:

$$ w_{\alpha} = \frac{\alpha_{f} - \alpha_{N}}{\alpha_{P} - \alpha_{N}} \quad (22) $$

This weight is then used in relation (18) to obtain the face value for the concentration field:

$$ C_{k,f} = C_{k,N} + w_{\alpha}(C_{k,P} - C_{k,N}) \quad (23) $$

This computation of the face weights is consistent with isoAdvector because it will enforce proportionality between the face flux of liquid and the face flux of concentration. However, this method is not applicable to a geometric VOF method such as isoAdvector. A simple example can illustrate cases where the face weight $w_\alpha$ is not defined.

In the idealized case depicted in Fig. 6, two neighbouring cells P and N with a volume fraction of 0.5 share a face f. A similar situation might occur in cases with a large curvature for instance. The portion of face f that is immersed in the liquid phase is 0.75. Therefore, $\alpha_f$ is not included between the values of cells P and N. For isoAdvector, this is not an issue since $\alpha_f$ comes from relation (9). The undefined weights can be replaced by some weights calculated with an UD scheme for example. UD is preferred because it is unconditionally bounded and stable. This solution has been implemented and tested. In practice, there were far too many faces where the undefined weights were replaced using UD leading to an over-diffused interface, with a lot of inconsistencies and artificial mass transfer. Using a hybrid formulation with an already existing interpolation scheme was not satisfactory either. The new methodology uses the phase average concentration. Since isoAdvector calculates the amount of liquid that flows through a face, it is consistent to use the phase average concentration from relation (13) to determine the amount of species k that flows through a face during a time step. The phase average is determined from the volume average, which is stored at a cell centre and from relation (14):

$$ \langle C_{k,1} \rangle_{V_1} = \frac{\langle C_{k,1} \rangle_{V}}{\alpha} \quad (24) $$

To obtain the face value, the cell phase average concentration is multiplied by $\alpha_f$, and gives the amount of species k in phase 1 that flows through face f:

$$ (C_{k,1})_f = \alpha_f \frac{\langle C_{k,1} \rangle_V}{\alpha} \quad (25) $$

Following this, the question of which cell centre value $(\langle C_{k,1} \rangle_{V_P}$ or $(\langle C_{k,1} \rangle_{V_N})$ should be used to compute the face value is raised. In line with the isoAdvector approach, the direction of the flow is followed so that the upwind cell is used. The basic algorithm is shown below in Fig. 7.

In this algorithm, alphaPhi is the flux returned by isoAdvector. The method can be applied to any solver (e.g. interFoam) as long as alphaPhi is accessible. The algorithm works in parallel and the basic communication function was modified to exchange data only when it is needed. In the bulk of each phase, a MUSCL scheme is applied to limit numerical diffusion, however, any higher order scheme could be used. The tolerance called “alphaTol” should be chosen accordingly to the tolerances of isoAdvector; these tolerances are used to determine whether the liquid volume in a cell should be considered as meaningful or not. If the liquid volume is not meaningful, it can be regarded as noise left by the advection process and the algorithm should not try to transport concentration in the cell. Since relation (25) is divided by the volume fraction, the scheme can become very oscillatory, especially when the volume fraction field is wiggly. In order to prevent these issues, it was necessary to work at low CFL numbers ($\le 0.1$) and to use a limiter, for which the development will be presented later.
scalar Cf = 0;
if (alphaf > alphaTol)
{
    if (alphaPhi > 0 && alphaOwn > alphaTol)
    {
        Cf = alphaf * COwn / alphaOwn;
    }
    if (alphaPhi < 0 && alphaNei > alphaTol)
    {
        Cf = alphaf * CNei / alphaNei;
    }
}
else
{
    Cf = 0;
}
return Cf;

<figure><img src="image_12.png" /><figcaption>Fig. 8. Simple 1D uniform structured stencil for the integration of the diffusion equation.</figcaption></figure>

### 3.3. Discretisation of the diffusive part of the transport equation

The second part of the transport equation consists in a diffusion equation:

$$ \frac{\partial C_{k,1}}{\partial t} = \nabla \cdot (D_{k,1} \nabla C_{k,1}) \qquad (26) $$

The diffusion equation must be discretised implicitly in time:

$$ \langle C_{k,1} \rangle_{V_i}^{n+1} = \langle C_{k,1} \rangle_{V_i}^{n} + \frac{\Delta t}{V_i} \times \sum_f (D_{k,1} \mathbf{n} \cdot \nabla C_{k,1}^{n+1} A)_f \qquad (27) $$

Two things must be considered carefully. Firstly, the diffusion coefficient must be defined consistently with the two-field approach. The liquid diffusivity of species k must be zero in the gas phase for example. Besides, the concentration gradient must be examined thoroughly. Indeed, the difference between the volume average and the phase average can lead to a non-zero face concentration gradient. When this gradient is computed based on the volume average concentration, it can be non-zero even though there is actually no difference in concentration in the phase. These circumstances can occur on a face separating a bulk cell from a cell cut by the interface. Fig. 8 shows an example of three 1D cells of a structured uniform mesh. The interface between the gas G and the liquid L is shown in red.

The gradients of volume average concentration on faces w and e are computed as:

$$ \mathbf{n}_w \cdot (\nabla C_{k,1})_{w}^{n+1} = \frac{C_{k,1,P}^{n+1} - C_{k,1,W}^{n+1}}{\Delta x} \qquad (28) $$

$$ \mathbf{n}_e \cdot (\nabla C_{k,1})_{e}^{n+1} = \frac{C_{k,1,E}^{n+1} - C_{k,1,P}^{n+1}}{\Delta x} \qquad (29) $$

The interface cell P is partially filled with liquid and therefore $C_{k,1,P}^{n+1} < C_{k,1,E}^{n+1}$ (since there is no mass transfer). However, $C_{E,L} = C_{P,L}$ because the species is perfectly mixed in the liquid phase. There is no concentration gradient in the liquid and therefore, no chemical diffusion should happen. If the diffusion equation is solved using the volume average concentration, the difference in volume fraction will create artificial diffusion. To tackle this problem, the phase average must be used to compute the diffusive fluxes. This is
<figure><img src="image_13.png" /><figcaption>Fig. 9. Flowchart of the global algorithm for the transport of concentration.</figcaption></figure>

consistent with the two-field approach since it is the concentration gradient within the phase that must promote chemical diffusion. The new equation can be solved to obtain the phase average concentration:

$$ \langle C_{k,1} \rangle_{V_{i,1}}^{n+1} = \langle C_{k,1} \rangle_{V_{i,1}}^{n} + \frac{\Delta t}{V_i} \times \sum_{f} (D_{k,1} \mathbf{n} \cdot \nabla C_{k,1}^{n+1} A)_{f} \quad (30) $$

$V_{i,1}$ corresponds to the volume of phase 1 inside cell $i$. The new face gradients in Fig. 8 are expressed in terms of phase average concentration:

$$ \mathbf{n}_{iW} \cdot (\nabla C_{k,1})_{iW}^{n} = \frac{1}{\Delta x} \left( \frac{C_{k,1,P}^{n}}{\alpha_{P}^{n}} - \frac{C_{k,1,W}^{n}}{\alpha_{W}^{n}} \right) \quad (31) $$

$$ \mathbf{n}_{e} \cdot (\nabla C_{k,1})_{e}^{n} = \frac{1}{\Delta x} \left( \frac{C_{k,1,E}^{n}}{\alpha_{E}^{n}} - \frac{C_{k,1,P}^{n}}{\alpha_{P}^{n}} \right) \quad (32) $$

Once the diffusion equation is solved and the new time step phase average concentration is obtained, the new time step volume average concentration is computed using relationship (14), together with the new time step volume fraction field.

### 3.4. Final description of the algorithm

The final implementation of the new algorithm uses both volume average and phase average concentrations, as it was presented in the previous section 3.3. This implementation also integrates limiters to avoid over- and under-shoots of the concentration field. The flowchart describing the algorithm is shown in Fig. 9.
The major improvement of this algorithm is the statement “bounded”, which is used to verify if the advection step has produced unboundednesses in the concentration field. Limiting the face concentration is not trivial because it can take any positive value depending on the system, the flow etc. To derive a limiter, the definition of the phase average concentration and the thermodynamic value of the saturation of the liquid phase $C_{k,1}^*$ are used. The phase average represents the concentration of a given species in the liquid or the gas phase of the cell. If there is only one species, it will be equivalent to density since the species is alone in its phase. On the other hand, the saturation is a thermodynamics limit which imposes the maximum concentration in a mixture. This depends on the other species in the mixture, but also on temperature and pressure. The ratio between the phase average concentration and the saturation limit is computed. This ratio represents the quantity of the species the phase can still accept. Finally, the creation of under- or over-shoots by the advection process can be evaluated. This is useful when dealing with cells that are close to saturation or almost empty of fluid. However, it is not possible to limit the concentration of a cell when the species is still being transferred and the phase is far from being saturated because there are no ways of limiting the concentration other than this physical limit. In the algorithm, the logical statement “bounded” uses the degree of filling of phase 1 to assess the “space” left for species k:

$$ \text{degree of filling} = \frac{\langle C_{k,1} \rangle_{V_1}}{C_{k,1}^*} \qquad (33) $$

The sequence of the “bounded” statement is described with the flowchart in Fig. 10.

The important part in the correction of the advection step is the evaluation of unboundedness. The bounded nature of the “degree of filling” is used in order to verify how the interpolation scheme produces over- and under-shoots. Since the degree of filling comes from the phase average concentration, the interpolation scheme is slightly modified, and the degree of filling value in the upwind cell is multiplied by the face values of the volume fraction field. The result gives information on how much the face is “saturated” by species k. The degree of filling of the face cannot be negative nor can it be greater than one. There are two options to bound the face values. The first is to clip values that are negative and greater than unity. However, this may create a conservation error. The second option is to redistribute the excess of concentration (i.e. values over the saturation concentration) in the neighbouring cells. In the current implementation, the clipping option was chosen for two reasons. Firstly, because isoAdvector itself offers a clipping option and it is compulsory to be as consistent as possible with isoAdvector. The other reason comes from the origin of the unboundednesses. The error occurs because the volume average concentration is divided by the volume fraction field to obtain the phase average concentration. When the cells are almost empty and the volume fraction is close to zero, any small error on the volume fraction has a large impact on the phase average concentration. Redistributing this error does not result in better advection since the error does not have any physical meaning in the first place.

# 4. Validation of the numerical method

This section will detail several test cases set up to demonstrate the performance of the methodology presented and its compatibility with isoAdvector.

## 4.1. Validation of the advective transport

The discretisation of the advection operator and its consistency with isoAdvector is tested first. The test case uses a dilute species with a uniform and constant concentration in the liquid phase equal to the saturation concentration $C_{k,1}^* = 1 \text{ kg m}^{-3}$. The concentration gradients are initially zero such that the case is a pure advection case. The volume average concentration field must be equal to $\alpha C_{k,1}^*$ at every time step. To evaluate this condition, the difference $\alpha - \frac{\langle C_{k,1} \rangle_{V_1}}{C_{k,1}^*}$ must be monitored. This test case uses the fact that the phase average concentration in the cell should stay equal to the saturation concentration at all times. The set up is the same as the one presented earlier in Fig. 4.

Fig. 11 compares the custom interpolation scheme with the MUSCL in a 1D context as presented in Fig. 5. The custom interpolation scheme provides a much more accurate concentration profile than the classical MUSCL scheme. The black and blue curves cannot be differentiated, hence proving that the custom scheme yields a concentration profile that is perfectly consistent with the volume fraction.

The second test case presented corresponds to the same problem as the first test case but for a Taylor flow in two dimensions. Increasing the number of dimensions renders the problem more challenging, due to more degrees of freedom and to the surface tension force calculation.

Fig. 12 compares the liquid phase concentration field obtained after advection with three different interpolation schemes in a 2D channel with 50 cells per diameter. The numerical diffusion induced by the UD scheme creates a lot of artificial mass transfer as it can be seen on the diffused concentration profile at the top of Fig. 12. Even with a second order MUSCL scheme, a lot of artificial mass transfer arises. Nevertheless, with the interpolation scheme developed in this work, the concentration field coincides strictly with the interface, as we will see.

In order to assess the performances of the new interpolation scheme, the $L_1$ norm is introduced to measure the difference between the volume fraction field and the ratio between the volume average concentration and the saturation concentration:

$$ L_1 = \sum_{i}^{n} \left| \alpha_i - \frac{\langle C_{k,1} \rangle_{V_i}}{C_{k,1}^*} \right| \qquad (34) $$
<figure><img src="image_14.png" /><figcaption>Fig. 10. Flowchart for the correction of the advection process.</figcaption></figure>

In equation (34), n is the total number of mesh cells. If the advection scheme is consistent with isoAdvector, then $\langle C_{k,1} \rangle_{V_i} = \alpha_i \langle C_{k,1} \rangle_{V_{i,1}}$. Since the relation $\langle C_{k,1} \rangle_{V_{i,1}} = C_{k,1}^*$ is imposed, the $L_1$ norm should be small. To compare multiple refinement levels, $L_1$ is divided by the total number of mesh cells to obtain what is called “mean error”:

$$ \text{mean error} = \frac{L_1}{n} \qquad (35) $$

Fig. 13 compares the mean error generated over time by three interpolation schemes on the liquid phase concentration field of a two-dimensional Taylor flow in a microchannel. The artificial mass transfer arising from numerical diffusion in Fig. 12 is quantified with the time evolution of the error. The error created by the UD scheme and the MUSCL scheme is around 0.01 which is seven orders of magnitude larger than the error created by the custom interpolation scheme developed in this work which is around $10^{-9}$. Fig. 14 shows the mean error resulting from advection for different mesh refinement levels. Even though the error does increase with time for the coarsest mesh, its maximum value remains very small compared with the saturation concentration (1 kg m$^{-3}$). Increasing the mesh resolution decreases the error as it can be seen by the red and blue curves. However, such a level of refinement
<figure><img src="image_15.png" /><figcaption>Fig. 11. Comparison between the volume fraction field profile and nondimensionalized concentration field profile obtained with the custom interpolation scheme.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Liquid phase concentration field in a 2D Taylor flow obtained with a UD scheme (top), a MUSCL scheme (middle) and the custom interpolation scheme (bottom) developed in this work.</figcaption></figure>

has a significant impact on numerical cost and a compromise should be made with respect to the users needs.

The final test case uses a 3D channel with a square cross section of 1 mm. The refinement is the same as the 2D case in the axial direction. The base mesh is made of 50 cells across the width and height of the channel. Fig. 15 compares the mean error for the 2D and the 3D cases. The 3D simulation gives lower error than the 2D case and this may be explained by a better interface advection by *isoAdvector* in 3D than in 2D.

The results obtained with the custom interpolation scheme show that it can be used to advect a concentration field with minimal artificial mass transfer. Increasing mesh refinement seems to enhance the results, which is consistent with the fact that it is based on *isoAdvector*. However, a higher level of refinement is at the expense of a larger numerical cost.

## 4.2. Validation of the diffusive transport

Since there is no molecular diffusion in the VOF equation (2), the diffusion coefficient in the species concentration equation must be corrected to make it consistent with a two-field approach. If the diffusion coefficient is not updated at each new time step, diffusive fluxes can make the phase average concentration become non-zero in cells that do not contain phase 1.

To test the method proposed in section 3.3, a structured 1D mesh is used. The channel length L = 0.01 m has 1000 cells. The gas-liquid interface is initially located at 0.2505 L and moves at 0.01 m s⁻¹. There are three species: the liquid is water, the gas is air and there
<figure><img src="image_17.png" /><figcaption>Fig. 13. Comparison of the error evolution with time for the three different interpolation schemes of the 2D case presented in Fig. 12.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Evolution of the mean error with time for three levels of refinement of a two-dimensional Taylor flow.</figcaption></figure>

is also a species diluted only in water. Water and air are treated by the VOF algorithm itself, while the dilute species is handled by a concentration equation as discussed in the former sections. The diffusion coefficient of the diluted species is set to $1 \times 10^{-6} \text{ m}^2 \text{ s}^{-1}$. This value is much larger than typical physical values for diffusion coefficients in liquids in order to amplify any possible issues. The diluted species in water will not cross the interface because no mass transfer occurs. Therefore, the diluted species is confined to the liquid phase. Its saturation concentration is arbitrarily set to $1 \text{ kg m}^{-3}$. The concentration of the diluted species is initialized as a step at the saturation concentration between 0.1 L and 0.2 L and zero everywhere else. Since the diffusion coefficient is large, the species will rapidly diffuse towards the interface and towards the inlet of the channel.

As explained in section 3.3, the phase average concentration and the volume average concentration differ because the volume of liquid is smaller than the volume of the cell, in a cell cut by the interface. If the classical FVM volume average variable is used to solve the diffusion equation, the concentration profile given in Fig. 16 is obtained. The volume average concentration, which is computed using the volume average concentration for diffusion and is indicated by the red curve, is compared with the volume average concentration computed using the phase average concentration, shown by the blue. The position of the interface is shown by the black curve indicated as “alpha”. The global shapes of the curves are very close, however, the profiles near the interface are quite different. The red profile displays a decrease before showing a spike at the interface that comes from the direction of the face gradient. Since the volume fraction in the interface cell is smaller than unity, the volume average concentration of the interface cell is smaller than the neighbouring cell volume average concentration. This is true before reaching steady state since the concentration is initially zero near the interface and there is a transition time before the diffusion from the step profile reaches the interface. Once steady state has been reached, however, the maximum concentration will be at the interface because the species cannot evaporate into air. At 0.5 s, the concentration profile is established such that the interface concentration is at the maximum. There is a volume
<figure><img src="image_19.png" /><figcaption>Fig. 15. Comparison of the mean error evolution between a two-dimensional and a three-dimensional Taylor flow.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. Non-dimensional volume average concentration profiles at 0.5 s.</figcaption></figure>

average concentration gradient near the interface because there is a difference of volume fractions (see relation (14)), and thus, the concentration profile shows a kind of wave near the interface. Apart from close to the interface region, the profile looks as expected. However without diffusion correction, the profile is badly predicted. This appears clearer when the phase average is plotted instead of the volume average, as shown in Fig. 17.

It can be seen in the inset of Figs. 16 and 17 that the slope of the blue line is not vertical. This is because the volume fraction field changes from 1 to 0 over three cells, even though the interface is located in only one cell. In Fig. 17, the red curve displays the same kind of decrease in concentration before the interface. Indeed, the volume fraction in this region equals one, which is why the phase average and volume average are equal. At the interface, the spike in the concentration profile is much larger than in the volume average concentration profile of Fig. 16. This is expected since the volume fraction is smaller than one in the interface cell. Nonetheless, a spike of concentration at the interface is not physical since diffusion is supposed to smoothen the concentration profile with time.

In both Figs. 16 and 17, the blue curve entitled “diffusion correction” displays a much more accurate concentration profile near the interface. These blue curves were obtained when the diffusion equation is solved in terms of phase average. Therefore, there is no gradient of phase average due to the change in volume fraction in the interface cell. This basic test case shows the importance of the phase average when solving for diffusion in a two-field approach. For a single field approach, the issue is not the same since the concentration field is unique, but not necessarily continuous.
<figure><img src="image_21.png" /><figcaption>Fig. 17. Non-dimensional phase average concentration profiles at 0.5 s.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 18. Gas bubble in a diagonal uniform flow with constant velocity $U = (1, 0.5)$.</figcaption></figure>

### 4.3. Disc in steady uniform two-dimensional flow

In part 4.1, test cases with relatively simple flow conditions and uniform structured meshes were used. The present section aims at demonstrating the capabilities of the proposed method in more realistic conditions. The same error evaluation based on equation (34) and the mean error will be used. The test case is the same as proposed by Roenby et al. [1] in their original work on *isoAdvector*. A gas bubble moves in a rectangular domain with a diagonal, uniform and constant velocity field as illustrated in Fig. 18. The concentration of this bubble is set identically to section 4.1.

An unstructured mesh made of triangular 2D cells provided by GMSH [53], based on a Delaunay standard algorithm, is compared to the uniform, structured, mesh generated with `blockMesh`. Three levels of refinement are used. For the square mesh, the edge length is uniform in all directions. For the triangle mesh, the same edge length is prescribed on the boundaries. The unstructured meshes display a larger cells count for the same edge length.

In Fig. 19 we compare three different edge lengths of 0.01 m, 0.02 m and 0.04 m for both square and triangular meshes. Apart from the square mesh with 0.04 m long edges, which gives a surprisingly low deviation, all cases yield results of a mean error between $1 \times 10^{-13}$ and $1 \times 10^{-9}$. As these simulations on a 2D triangular mesh give satisfactory results, the method proposed in this work can be used for various configurations. However, we need to keep in mind that the error displayed on Fig. 19 quantifies the artificial mass transfer by comparing the concentration with the volume fraction computed by *isoAdvector*. Therefore, a small deviation on Fig. 19 does not mean that the results are realistic, but only that the method does not introduce additional error in the simulation.
<figure><img src="image_23.png" /><figcaption>Fig. 19. Comparison of the mean error evolution for three levels of refinement between square and triangular 2D meshes.</figcaption></figure>

## 5. Conclusions

This work has been devoted to the implementation of a transport equation for the concentration of chemical species in a Volume-Of-Fluid framework. It was shown that existing interpolation schemes such as UD, CD or even MUSCL used to compute face values from cell centered values are not adapted to the geometric nature of the *isoAdvector* method and lead to over-diffused concentration profiles and artificial mass transfer. To overcome this issue and avoid inconsistencies, a new methodology to interpolate scalar values from cell centres to cell faces has been developed. This methodology uses the phase average concentration to compute the face values in an adequate way, which allows the concentration fields to be non-zero only in their corresponding phases. Numerical tools to test the consistency of the new methodology with *isoAdvector* have been described. Using the $L_1$ norm and the mean error, it was shown that the consistency error caused by advection is very low with the new methodology. It has also been shown that increasing mesh refinement provides better results, however, the additional numerical cost induced by refinement is not always necessary since the consistency error caused by advection is already very small.

Additionally, it was shown that the diffusive fluxes computed using the classical volume average variables can lead to non-physical results. Thus, a new formulation of concentration gradients based on the phase average concentration has been proposed. Comparison between corrected and non corrected fluxes has been provided with a test case and it was shown that the use of volume average concentration leads to a non-physical concentration profile in the interface region.

In future work, the transport equation for concentration will be extended to include additional source terms. Species transfer at the interface will be computed based on the local concentration at the interface. The implementation of chemical reaction with power law kinetics will also be investigated.

## CRediT authorship contribution statement

**Alexis Tourbier:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Lionel Gamet:** Writing – review & editing, Validation, Supervision, Resources, Project administration, Funding acquisition, Conceptualization. **Philippe Béard:** Writing – review & editing. **Typhène Michel:** Writing – review & editing. **Joelle Aubin:** Writing – review & editing, Validation, Investigation, Conceptualization. **Hrvoje Jasak:** Validation, Methodology, Investigation.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgements

The authors thank IFP Energies nouvelles for funding this research and providing HPC resources.
References

[1] J. Roenby, H. Bredmose, H. Jasak, A computational method for sharp interface advection, R. Soc. Open Sci. 3 (2016) 160405.

[2] M.N. Kashid, L. Kiwi-Minsker, Microstructured reactors for multiphase reactions: state of the art, Ind. Eng. Chem. Res. 48 (2009) 6465-6485.

[3] J. Pérez-Ramírez, R.J. Berger, G. Mul, F. Kapteijn, J.A. Moulijn, The six-flow reactor technology: a review on fast catalyst screening and kinetic studies, Catal. Today 60 (2000) 93-109.

[4] C. de Bellefon, R. Abdallah, T. Lamouille, N. Pestre, S. Caravieilhes, P. Grenouillet, High-throughput screening of molecular catalysts using automated liquid handling, injection, and microdevices, Chimia 56 (2002) 621.

[5] J. Keybl, K.F. Jensen, Microreactor system for high-pressure continuous flow homogeneous catalysis measurements, Ind. Eng. Chem. Res. 50 (2011) 11013-11022.

[6] C.P. Park, M.M. Van Wingerden, S.-Y. Han, D.-P. Kim, R.H. Grubbs, Low pressure ethenolysis of renewable methyl oleate in a microchemical system, Org. Lett. 13 (2011) 2398-2401, PMID: 21456591.

[7] L. Vanoye, M. Pablos, N. Smith, C. de Bellefon, A. Favre-Réguillon, Aerobic oxidation of aldehydes: selectivity improvement using sequential pulse experimentation in continuous flow microreactor, RSC Adv. 4 (2014) 57159-57163.

[8] L. Vanoye, M. Pablos, C. de Bellefon, A. Favre-Réguillon, Gas-liquid segmented flow microfluidics for screening copper/tempo-catalyzed aerobic oxidation of primary alcohols, Adv. Synth. Catal. 357 (2015) 739-746.

[9] M. Kamaleddine, C. Bonnin, T. Michel, L. Brunet-Errard, J. Aubin, L. Prat, Gas-liquid flow characterization and mass transfer study in a microreactor for oligomerization catalyst testing, Chem. Eng. Process. Proc. Intensif. 166 (2021) 108476.

[10] S. Fleckenstein, D. Bothe, A volume-of-fluid-based numerical method for multi-component mass transfer with local volume changes, J. Comput. Phys. 301 (2015) 35-58.

[11] C. Hirt, B. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201-225.

[12] J.G. Khinast, Impact of 2-D bubble dynamics on the selectivity of fast gas-liquid reactions, AIChE J. 47 (2001) 2304-2319.

[13] B. Aboulhasanzadeh, S. Thomas, M. Taeibi-Rahni, G. Tryggvason, Multiscale computations of mass transfer from buoyant bubbles, Chem. Eng. Sci. 75 (2012) 456-467.

[14] J.G. Khinast, A.A. Koynov, T.M. Leib, Reactive mass transfer at gas-liquid interfaces: impact of micro-scale fluid dynamics on yield and selectivity of liquid-phase cyclohexane oxidation, Chem. Eng. Sci. 58 (2003) 3961-3971.

[15] S. Radl, A. Koynov, G. Tryggvason, J.G. Khinast, DNS-based prediction of the selectivity of fast multiphase reactions: hydrogenation of nitroarenes, Chem. Eng. Sci. 63 (2008) 3279-3291.

[16] Z. Tukovic, H. Jasak, Simulation of free-rising bubble with soluble surfactant using moving mesh finite volume/area method, in: Proceedings of 6th International Conference on CFD in Oil & Gas, Metallurgical and Process Industries, 2008, no. CFD08-072.

[17] K. Bäumler, Simulation of single drops with variable interfacial tension, Ph.D. thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), 2014.

[18] C. Lehrenfeld, On a space-time extended finite element method for the solution of a class of two-phase mass transport problems, Ph.D. thesis, Universitätsbibliothek der RWTH Aachen, 2015.

[19] P.S. Weber, H. Marschall, D. Bothe, Highly accurate two-phase species transfer based on ALE interface tracking, Int. J. Heat Mass Transf. 104 (2017) 759-773.

[20] C. Yang, Z.-S. Mao, Numerical simulation of interphase mass transfer with the level set approach, Chem. Eng. Sci. 60 (2005) 2643-2660.

[21] K.B. Deshpande, W.B. Zimmerman, Simulations of mass transfer limited reaction in a moving droplet to study transport limited characteristics, Chem. Eng. Sci. 61 (2006) 6424-6441.

[22] VOF-Simulations of Mass Transfer from Single Bubbles and Bubble Chains Rising in Aqueous Solutions, Symposia, Parts A, B, and C of Fluids Engineering Division Summer Meeting, vol. 2, 2003.

[23] D. Bothe, M. Koebe, K. Wielage, J. Prüss, H.-J. Warnecke, Direct Numerical Simulation of Mass Transfer Between Rising Gas Bubbles and Water, Springer Berlin Heidelberg, Berlin, Heidelberg, 2004, pp. 159-174.

[24] Y. Haroun, D. Legendre, L. Raynal, Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film, Chem. Eng. Sci. 65 (2010) 2896-2909.

[25] 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.

[26] D. Bothe, S. Fleckenstein, A Volume-of-Fluid-based method for mass transfer processes at fluid particles, Chem. Eng. Sci. 101 (2013) 283-302.

[27] D. Deising, H. Marschall, D. Bothe, A unified single-field model framework for Volume-of-Fluid simulations of interfacial species transfer applied to bubbly flows, Chem. Eng. Sci. 139 (2016) 173-195.

[28] A. Weiner, D. Bothe, Advanced subgrid-scale modeling for convection-dominated species transport at fluid interfaces with application to mass transfer from rising bubbles, J. Comput. Phys. 347 (2017) 261-289.

[29] H.G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys. 12 (1998) 620-631.

[30] S.S. Deshpande, L. Anumolu, M.F. Trujillo, Evaluating the performance of the two-phase flow solver interfoam, Comput. Sci. Discov. 5 (2012) 014016.

[31] H. Scheufler, J. Roenby, Accurate and efficient surface reconstruction from volume fraction data on general meshes, J. Comput. Phys. 383 (2019) 1-23.

[32] J. Brackbill, D. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335-354.

[33] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annu. Rev. Fluid Mech. 31 (1999) 567-603.

[34] J. Maes, C. Soulaine, A unified single-field Volume-of-Fluid-based formulation for multi-component interfacial transfer with local volume changes, J. Comput. Phys. 402 (2020) 109024.

[35] S. Popinet, Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries, J. Comput. Phys. 190 (2003) 572-600.

[36] G. Gennari, R. Jefferson-Loveday, S.J. Pickering, A phase-change model for diffusion-driven mass transfer problems in incompressible two-phase flows, Chem. Eng. Sci. 259 (2022) 117791.

[37] P.K. Farsoiyi, S. Popinet, L. Deike, Bubble-mediated transfer of dilute gas in turbulence, J. Fluid Mech. 920 (2021) A34.

[38] J. López-Herrera, A. Ganán-Calvo, S. Popinet, M. Herrada, Electrokinetic effects in the breakup of electrified jets: a volume-of-fluid numerical study, Int. J. Multiph. Flow 71 (2015) 14-22.

[39] E. Cipriano, A. Essamade Saufi, A. Frassoldati, T. Faravelli, S. Popinet, A. Cuoci, Multicomponent droplet evaporation in a geometric Volume-of-Fluid framework, J. Comput. Phys. 507 (2024) 112955.

[40] J.B. Bell, P. Colella, H.M. Glaz, A second-order projection method for the incompressible Navier-Stokes equations, J. Comput. Phys. 85 (1989) 257-283.

[41] scalartransportflow, https://github.com/alexistourbier/scalarTransportFlow.git, 2024.

[42] C. Bilger, M. Aboukhedr, K. Vogiatzaki, R. Cant, Evaluation of two-phase flow solvers using level set and volume of fluid methods, J. Comput. Phys. 345 (2017) 665-686.

[43] L. Gamet, M. Scala, J. Roenby, H. Scheufler, J.-L. Pierson, Validation of volume-of-fluid OpenFOAM® isoAdvector solvers using single bubble benchmarks, Comput. Fluids 213 (2020) 104722.

[44] F. Jamshidi, H. Heimel, M. Hasert, X. Cai, O. Deutschmann, H. Marschall, M. Wörner, On suitability of phase-field and algebraic volume-of-fluid OpenFOAM solvers for gas-liquid microfluidic applications, Comput. Phys. Commun. 236 (2019) 72-85.
[45] J. Pedersen, B. Larsen, H. Bredmose, H. Jasak, A new Volume-of-Fluid method in OpenFOAM, in: M. Visonneau, P. Queutey, D. Le Touzé (Eds.), MARINE 2017 Computational Methods in Marine Engineering VII, 7th International Conference on Computational Methods in Marine Engineering, MARINE 2017, International Center for Numerical Methods in Engineering, Conference Date: 15-05-2017 Through 17-05-2017, 2017, pp. 266–278.

[46] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (1998) 112–152.

[47] *Twophaseflow*, https://github.com/DLR-RY/TwoPhaseFlow.git, 2022.

[48] B. Van Leer, Towards the ultimate conservative difference scheme I. The quest of monotonicity, in: Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics, Fundamental Numerical Techniques July 3–7, 1972 Universities of Paris VI and XI, in: General Lectures, vol. I, Springer, 1972, pp. 163–168.

[49] B. van Leer, Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme, J. Comput. Phys. 14 (1974) 361–370.

[50] B. Van Leer, Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow, J. Comput. Phys. 23 (1977) 263–275.

[51] B. Van Leer, Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, J. Comput. Phys. 23 (1977) 276–299.

[52] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method, J. Comput. Phys. 32 (1979) 101–136.

[53] C. Geuzaine, J.-F. Remacle, Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Methods Eng. 79 (2009) 1309–1331.