Predicting Mass Transfer in Liquid–Liquid Extraction Columns

Abstract

In this work, the GEneralised Multifluid Modelling Approach (GEMMA) is applied to the simulation of liquid–liquid extraction in a Rotating Disc Column (RDC) and a Pulsed Sieve-plate Extraction Column (PSEC). A mass transfer modelling methodology is developed, in which the multiphase flows, droplet size distribution and dispersed phase holdup predicted with computational fluid dynamics are coupled to mass transfer correlations to predict the overall mass transfer. The numerical results for the stage-averaged dispersed phase holdup, Sauter mean droplet diameter and axial solute concentration in the RDC and PSEC agree with experimental observations. The proposed modelling method provides an accurate predictive tool for complex multiphase flows, such as those observed in intensified liquid–liquid extraction, and provides an alternative approach to column design using empirical correlations or pilot plant study.

1. Introduction

Multiphase flows are found in countless industrial applications, spanning from power generation and chemical processes to food production and biomedical applications. In the context of the nuclear industry, multiphase flows also play a crucial role in current and advanced reprocessing technologies foreseen for the next generation of nuclear power plants [1]. Therefore, predicting the local and global transient evolution of multiphase flows is recognised to be of paramount importance for the nuclear industry at large [2].

Computational fluid dynamics (CFD) has the potential to provide improved predictive capability of multiphase flows. This, however, is hindered by the complexity of the multiphase flows that are of interest to industrial applications. In particular, most multiphase flows of practical interest exhibit a broad range of interfacial scales, ranging from small dispersed phase elements (DPEs), such as particles, droplets or bubbles embedded in a continuous phase, to large interfaces observed in segregated free-surface flows. Off-the-shelf multiphase flow numerical models generally assume either small or large interfacial scales, resulting in the so-called interface-averaging and interface-resolving approaches [3,4]. The former method is mainly used for dispersed flows where the interfacial scales are smaller than the size of the numerical grid; scale separation is assumed, and the governing equations are conditionally-averaged, resulting in the so-called multifluid formulation. Due to the averaging operation, suitable closures are needed to account for the interfacial transfer of momentum, heat and mass. Conversely, the interface-resolving approach assumes that the mesh size is small enough to allow for an adequate resolution of the morphology of the interface, which generally applies to large segregated interfaces; this approach leads to interface-tracking and interface-capturing models. In the former, the interface position is tracked in a Lagrangian fashion, whilst in the latter, it is reconstructed from a known indicator function; the well-known Volume of Fluid (VoF) approach [5] is an example of an interface-capturing model.

From the discussion above, it is clear that most multiphase flows of practical interest exhibit a marked multiscale behaviour, whilst standard multiphase modelling approaches assume the presence of either “small” dispersed interfacial scales or “large” segregated interfacial scales. In the numerical simulation context, the terms “small” and “large” are defined with respect to the size of the numerical grid. Therefore, it is evident that there is a need for a generalised multiphase modelling approach capable of handling the presence of a broad range of interfacial scales in the same computational domain; a number of these approaches have been proposed, mainly following the idea of embedding some form of large interface resolution within a standard multifluid framework [6]. However, most of these approaches either rely on a priori regime maps based on the local volume fraction [6,7,8] or lack the capability of adapting to the local flow regime altogether [9].

To overcome these shortcomings, the GEneralised Multifluid Modelling Approach (GEMMA) for the simulation of multiscale multiphase flows has been recently developed [10]. The approach adapts its formulation to the local resolution of the interfacial scales. The GEMMA approach reduces to a standard multifluid formulation suitable for small/dispersed interfaces in the numerical cells where the interfacial scales are small compared to the mesh size. In the cells where the mesh size is fine enough to guarantee an acceptable resolution of the interfacial morphology, a novel multifluid formulation suitable for the simulation of large/segregated interfaces is introduced; the latter formulation aims at mimicking the behaviour of an interface resolving approach such as VoF within the multifluid framework. However, given that the model is based on a multifluid description, dedicated closures for interfacial momentum transfer and surface tension remain necessary to describe the underlying physics of interfacial momentum exchange in large interface regions. The model has been assessed against different fundamental test cases in [10], where it has been shown that it is as accurate as the VoF approach for cases characterised by large/segregated interfaces, whilst a standard multifluid behaviour is recovered in dispersed flows; in the same work, the authors demonstrated the capability of the approach to adapt its formulation locally in a prototypical multiscale flow, i.e., a water jet plunging in a quiescent pool.

Furthermore, it has been demonstrated that the GEMMA approach can accurately represent the complex multiphase hydrodynamics encountered in liquid–liquid extraction devices [11,12]. Liquid–liquid extraction is when a solute is transferred between two immiscible fluids, typically a polar aqueous phase and a non-polar organic phase, with separation occurring based on relative solubility or chemical reaction at the interface. Liquid–liquid extraction can be performed batch-wise, where phases are sequentially contacted and then separated; however, this is often impractical when performing the process at scale. Performing liquid–liquid extraction at scale has traditionally employed large columns where the heavy phase enters the top before flowing down and exiting via the base. The less-dense phase of the two phases enters the column at the bottom, where it travels upwards before leaving at the top. A large surface-area-to-volume ratio between the two phases is desirable to ensure a high mass transfer rate, necessitating a small dispersed phase diameter; however, droplets must be sufficiently large to prevent entrainment and subsequent column flooding. By controlling the amount of turbulence within the column, it is possible to optimise the mass transfer performance. In the case of a rotating disc column (RDC), this is achieved by rotating discs attached to a central shaft using a variable speed motor, while in a pulse sieve-plate extraction column (PSEC), the column fluid is cyclically pulsed upwards and downwards through sieve-plates, resulting in jetting and droplet formation.

Traditionally, the design of liquid–liquid extraction columns has relied upon utilising some of the many published empirical correlations; however, these often perform poorly, resulting in large, expensive over-designed columns. Providing suitable time and research facilities are available, a pilot plant can be used to determine extraction as a function of column height and cross-section, which is then scaled accordingly. However, scaleup is a complicated process that can result in expensive overdesigned columns. Due to the inherent uncertainty associated with the empirical or pilot plant design of extraction columns, developing and utilising a modelling and simulation approach to design is desirable.

To date, there have been several investigations into CFD informed mass transfer modelling in liquid–liquid extraction columns, with several studies looking at application within an RDC [13,14] and a PSEC [15]. However, these were done with models unable to distinguish between the flow regimes. Instead, results in [11,12] demonstrated how an approach such as GEMMA can make available the key hydrodynamic parameters needed to characterise mass transfer in these applications. Although implementing mass transfer directly with GEMMA would provide a detailed understanding of performance, doing so would be computationally expensive as it requires solving the relevant mass transfer calculations in each cell in the computational domain. The present work builds on previous findings [11,12] and presents a modelling framework to evaluate the global mass transfer performance in complex multiphase systems. The mass transfer modelling framework relies on input from GEMMA to evaluate the hydrodynamic parameters that influence mass transfer; successively, a surrogate reduced-order model of the system is created based on the hydrodynamic information obtained from the CFD model. This surrogate model is used to infer the global mass transfer performance of the system. This is significantly quicker as mass transfer performance is only calculated in a small number of compartments as opposed to every cell within the mesh, and, therefore, it requires much less computational resource.

The paper is organised as follows: the GEMMA modelling approach and the integral mass transfer modelling methodology are described in Section 2; their application to an RDC and a PSEC are reported in Section 3. Finally, conclusions and future work are outlined in Section 4.

2. Numerical Modelling

2.1. Hydrodynamic Modelling

The GEMMA approach has been implemented in the well-known open-source CFD code OpenFOAM v7.0 [16,17], and a detailed description can be found in [10]. A high-level overview of the approach is provided here.

The GEMMA approach has been built on top of the standard multifluid modelling framework suitable for small/dispersed interfaces given by the OpenFOAM native reactingMultiphaseEulerFoam solver. GEMMA introduces two different formulations within the multifluid framework; in each cell of the computational domain, one of the two approaches is selected based on the numerical grid’s local capability to resolve the interface’s morphology. The two formulations are:

(1)

A standard multifluid formulation, suitable for small/dispersed interfacial scales: this approach is used in the cells where the local mesh size is larger than the local interfacial scales, and, therefore, it is not possible to directly resolve the morphology of the interface.

(2)

An ad-hoc multifluid formulation suitable for large/segregated interfacial scales: this approach is used in the cells where the local mesh size is smaller than the local interfacial scales and the mesh resolution is fine enough to guarantee an adequate resolution of the interface. This formulation aims to provide a form of interface resolution, similar to interface-resolving approaches, in the context of the multifluid framework.

A large interface identifier $C_{\alpha }$ is introduced to identify which formulation is used in each cell. $C_{\alpha }$ is equal to zero in the cells where the dispersed formulation is employed and is instead equal to one in the cells where the large interface formulation is used. A flow diagram describing the switching logic is shown in Figure 1. A detailed description of the logic controlling the local $C_{\alpha }$ value and of the multifluid formulation for large/segregated interfaces is provided in [10].

The Interface Resolution Quality (IRQ) index is a numerical indication of interface resolution and is a function of the local mesh size and the local interface curvature, $\kappa$. A user-specified critical IRQ value is used to determine whether there is sufficient interface resolution for $C_{\alpha }$ to be activated. IRQ and local interface curvature are defined as:

$$IRQ=\frac{2}{\sqrt[3]{V_{cell}}\kappa }$$

(1)

$$\kappa =-\nabla .\left(\frac{\nabla \alpha }{\left|\nabla \alpha \right|}\right)$$

(2)

To ensure $C_{\alpha }$ is only enabled in cells where an interface exists, minimum and maximum values for the dispersed phase volume fraction are specified. Finally, $C_{\alpha }$ is enabled if the calculated diameter of the dispersed phase is larger than the current cell size multiplied by a user-defined value, $\Gamma$.

In the case of adiabatic flows without mass transfer, the volume-averaged continuity equation for phase $x$ is:

$$\frac{\delta {\alpha }_x}{dt}+\nabla .\left(\alpha U_x\right)+\nabla .\left(U_c{\alpha }_x\left(1-{\alpha }_x\right)\right)=0$$

(3)

where the ${\alpha }_x\left(1-{\alpha }_x\right)$ term ensures the included compressive velocity term is only active in the presence of a large interface to maintain sharpness by preventing diffusion. The compressive velocity term, $U_c$, is:

$$U_c=C_{\alpha }\left|U_c\right|\frac{\nabla \alpha }{\left|\nabla \alpha \right|}$$

(4)

The corresponding momentum conservation equation is:

$$\frac{\delta {\alpha }_x{U}_x}{\delta t}+\nabla .\left({\alpha }_x{U}_x{U}_x\right)=-\frac{{\alpha }_x\nabla p}{{\rho }_x}+\nabla .\left({\nu }_x{\alpha }_x\nabla U_x\right)+{\alpha }_{x}g+\frac{\left(F_x+F_{st,x}\right)}{{\rho }_x}$$

(5)

where the interfacial exchange is described via the momentum exchange force, $F_x$, and the surface tension force, $F_{st,x}$. The underpinning phenomena for these are different depending on if the fluid is dispersed or segregated, and these are formulated accordingly. The formulation for a generic force is as follows:

$$F=\left(1-\left(1-C_{\alpha }\right)f_x-\left(1-C_{\alpha }\right)f_z\right)F_{LI}+\left(1-C_{\alpha }\right)f_x{F}_{xy}+\left(1-C_{\alpha }\right)f_z{F}_{zx}$$

(6)

The surface tension force, $F_{st,x}$, is formulated as:

$$F_{st,k}={\alpha }_x{\displaystyle \sum }_{i=1}^{n_k}\left(C_{{\alpha }_{k,i}}{\sigma }_{k,i}{\alpha }_{smooth}\nabla \alpha \frac{2\rho }{\nabla {\rho }_{k,i}}\right)$$

(7)

The RDC computational domain has been modelled in 3D as a 60° axisymmetric wedge with turbulence in both phases modelled with the Large Eddy Simulation (LES) approach with the standard Smagorinsky [18] closure used for the subgrid stresses. As LES resolves the grid-level velocity field fluctuations in the unsteady flow field, it is possible to determine the turbulence kinetic energy (TKE) and turbulence energy dissipation rate. Due to the large size of the PSEC and the requirement for mesh refinement at the plate region, a fully 3D LES was not deemed computationally feasible, and instead, the column was modelled as a 2D slice. However, this precludes the use of LES, and, therefore, both phases are instead modelled with the Unsteady Reynolds-averaged Navier–Stokes (URANS) approach [19]. As URANS can only solve for the averaged flow field and fluctuations in the velocity field, the TKE and its dissipation rate are modelled using the mixture k-ε model.

2.2. Reduced Population Balance

A population balance approach can be used within GEMMA to evaluate the DPE size distribution when working in dispersed-interface mode; this feature is particularly important in cases involving mass transfer since the interfacial area available for the interfacial exchanges is directly related to the Sauter mean diameter of the dispersed phase in the regions of small/dispersed interfaces. The formulation of GEMMA is fully compatible with the MUSIG [20] multigroup inhomogeneous population balance embedded within the reactingMultiphaseEulerFoam solver in OpenFOAM, which allows for the evaluation of the DPE’s diameter distribution.

A reduced-order population balance has been implemented within GEMMA as a less computationally intensive alternative to MUSIG. The One Primary One Secondary Particle Method (OPOSPM) [21] is used as a reduced population balance approach. OPOSPM allows for the conservation of two low-order moments, and the selection of these moments is arbitrary; however, the total number and volume concentrations are the most commonly employed moments to guarantee the conservation of the total number and mass of the DPEs. Since the total number and volume concentration are conserved, the population density is represented by a single particle (assumed to have a spherical shape) whose size is characterised by the diameter:

$$d_{30}=\sqrt[3]{\frac{{\displaystyle \sum }{d}^3}{{\displaystyle \sum }{d}^0}}=\sqrt[3]{\frac{6{\alpha }_d}{\pi N_{drop}}}$$

(8)

where ${\alpha }_d$ and $N_{drop}$ are the volume fraction and the particle number density of the dispersed phase, respectively, which are related to the zeroth and third moment of the distribution $d^0$ and $d^3$. Since the volume fraction is already known from the resolution of the standard multifluid governing equations, the OPOSPM formulation only requires the solution of one additional conservation equation for $N_{drop}$, which is given by:

$$\frac{\partial \left({\rho }_d{N}_{drop}\right)}{\partial t}+\nabla ·\left({\rho }_d{U}_d{N}_{drop}\right)={\rho }_d S$$

(9)

where the source term, $S$, is:

$$S=\left(n_{daughter}-1\right)r_{breakage}{N}_{drop}-\frac{1}{2} r_{coalescence}{N}_{drop}^2$$

(10)

Within the source term, droplet formation due to breakage is described as a function of the mean number of daughter droplets, $n_{daughter}$, which is assumed to be 2, and the breakage rate, $r_{breakage}$, given by the break-up model of Martinez-Bazan et al. [22]:

$$r_{breakage}=\frac{\sqrt{8.2{\left(ϵd_{30}\right)}^{\frac{2}{3}}-\frac{12\sigma }{{\rho }_c{d}_{30}}}}{4d_{30}}$$

(11)

Finally, droplet coalescence rate, $r_{coalescence}$, is modelled using the coalescence model of Prince and Blanch [23], where the initial film thickness, $h_0$, and final film thickness, $h_f$, are assumed to be 10⁻⁴ and 10⁻⁸ m, respectively [23]:

$$r_{coalescence}=1.409{ϵ}^{\frac{1}{3}} d_{30}^{\frac{1}{3}}\exp \left(-{\left(\frac{d_{30}}{2}\right)}^{\frac{5}{6}} \frac{{\rho }_c^{\frac{1}{2}} {ϵ}^{\frac{1}{3}}}{4{\sigma }^{\frac{1}{2}}} \ln \left(\frac{h_0}{h_f}\right)\right)$$

(12)

The OPOSPM formulation described above allows for the evaluation of $d_{30}$; however, knowledge of the Sauter mean diameter $d_{32}$ also allows evaluation of the interfacial area density as:

$$a=\frac{6{\alpha }_d}{d_{32}}$$

(13)

In the context of liquid–liquid extraction, Wardle [24] followed the approach of estimating the ratio $d_{30}/d_{32}$ from known DPE size distributions and combined this with the $d_{30}$ evaluated with the reduced population balance to infer the Sauter mean diameter. For liquid–liquid dispersions in Annular Centrifugal Contactors (ACCs), it has been observed that the $d_{30}/d_{32}$ ratio is consistently in the range of 0.75–0.8. Consistent with [24], a value of 0.76 is used throughout the simulations presented in this work.

2.3. Mass Transfer Modelling Approach

The GEMMA approach makes it possible to obtain key information on the hydrodynamic behaviour of complex multiscale multiphase flows that are of interest to the nuclear industry, as demonstrated in [11] for the case of intensified liquid–liquid extraction in ACCs and in [12] for pulsed sieve-plate extraction columns (PSECs). This section describes a methodology to predict the global mass transfer rates in these flows based on the hydrodynamic information provided by a detailed CFD model based on the GEMMA approach. A continuously stirred tank reactor (CSTR) represents a reasonable approximation for a single compartment of an extraction column, such as a PSEC and an RDC [25,26]. These systems operate in counter-current flow, and therefore, the resolution of a system of one-dimensional differential equations is needed to characterise them along their entire active length and account for effects such as backflow and axial dispersion. The model requires the molar flow rates and solute concentration for both phases at the inlet as input parameters. The description provided below is aimed at liquid–liquid extraction applications, but the methodology applies to any multiphase flow presenting mass transfer phenomena driven by concentration gradients.

The approach relies on the two-film theory [27], shown in Figure 2, to evaluate the mass transfer rate ($\mathrm{MTR}$), which is based on the total fluid volume, inclusive of dispersed and continuous phases, and is expressed in mol m⁻³ s⁻¹. The MTR of a solute from the bulk aqueous to the interface and from the interface to the bulk organic is given by:

$$MTR=k_{c}a\left(C_c-C_{c.i}\right)=k_{d}a\left(C_{d.i}-C_d\right)$$

(14)

where $k_c$ and $k_d$ are the mass transfer coefficients expressed in ms⁻¹ of the continuous and dispersed phases, respectively, $a$ is the interfacial area per unit volume of reactor in m² m⁻³ and $C$ is the solute concentration in mol m⁻³, with the subscripts $c$ and $d$ denoting if the concentration refers to the aqueous or dispersed phase and the subscript $i$ indicating if the concentration refers to the interface.

Assuming there is negligible resistance to mass transfer at the interface, the equilibrium distribution $K_{eq}$ can be calculated using Equation (15) below. The equilibrium assumption is frequently valid and is justified in most absorption applications; on the other hand, interfacial chemical reactions can be the limiting step in extraction with chemical reaction. Despite this, the equilibrium assumption is used in the present modelling methodology, and the inclusion of finite-rate interfacial chemistry will be part of future model development:

$$K_{eq}=\frac{C_{d,i}}{C_{c,i}}$$

(15)

Rearrangement and substitution of Equations (14) and (15) give Equation (16) for the rate of interphase mass transfer:

$$MTR=\left(\frac{1}{\frac{K_{eq}}{k_c}+\frac{1}{k_d}}\right)a\left(K_{eq}{C}_c-C_d\right)$$

(16)

It can be noted that a first coupling between the mass transfer and the CFD model is represented by the interfacial area density required in Equation (16), which is evaluated as the product of the interfacial area density $a$ and the volume of the reactor $V$. The interfacial area density is evaluated through the CFD model via Equation (13). Different correlations are available in the literature to evaluate the phase-specific mass transfer coefficients $k_d$ and $k_c$, which are usually expressed as a function of the Sherwood number $Sh$; good reviews of such correlations are given in [25,26], and more details on this are presented in Section 3.

The correlations express the Sherwood number as a function of the physical properties of the fluids (e.g., density, viscosity, diffusion coefficient) as well as a function of the hydrodynamic conditions within the system (e.g., droplet Reynolds number, the diameter of the dispersed phase, relative velocity between the two phases). An accurate hydrodynamic model of the system is needed to evaluate the latter parameters; this represents another level of coupling between the CFD and the mass transfer modelling methodology.

The above set of equations is applied to the extraction of acetone in a water/toluene/acetone system in an RDC and a PSEC. These cases are then used to validate the mass transfer modelling methodology against the experimental measurements of [25]. A compartment modelling approach will be utilised, with the unit operations being modelled as a network of well-mixed sub-volumes, with each stage represented by a single CSTR [28]. Due to the low solute concentrations and simple extraction mechanism based on relative solubility, the main assumptions underlying the integral mass transfer modelling approach described above are that:

(1)

The hydrodynamic behaviour of the system is not significantly impacted by mass transfer.

(2)

The interfacial chemistry is infinitely fast, i.e., interfacial concentrations are assumed to be equilibrium concentrations.

(3)

The saturation of the solvent is neglected.

3. Results and Discussion

3.1. Simulation of a Rotating Disc Column

The RDC considered in this simulation is taken from Garthe [24]. The column has a diameter of 0.08 m, consists of 59 compartments and each compartment has a height of 0.05 $\mathrm{m}$, resulting in a total active length of 2.95 m. The continuous phase is water, the dispersed phase is toluene and the solute is acetone; the physical properties of these components are reported in [25]. The organic and aqueous phases have a volumetric flow rate of 1.33 × 10⁻⁵ and 1.11 × 10⁻⁵ m³s⁻¹, respectively. The stirrer speed was 6.67 s⁻¹.

A simplified CFD model of the RDC with the operating conditions listed above was created following the approach used by [29]; a schematic representation of the computational domain is shown in Figure 3. The domain is three-dimensional and axisymmetric and comprises 59 compartments. A cut-cell mesh having a bulk size of 0.002 m and a size at the walls of 0.0008 m is used to discretise the computational domain, resulting in a total of 117,615 quadrilateral cells. A velocity inlet boundary condition was used at the bottom for both phases, whilst a pressure outlet condition was used at the top. For the volume fraction, values proportional to the volumetric flow rate of each phase were imposed at the bottom, and an inletOulet condition was used at the top; likewise, a fixed value of $N_{drop}$, corresponding to an inlet Sauter diameter equal to 3 × 10⁻³ m for toluene, was specified at the bottom, and a zero gradient condition was used at the top. The stirrer wall rotates around the axis with a rotating speed of 6.67 s⁻¹, whilst the stator wall is stationary; both walls were treated as no-slip walls. Interphase momentum transfer of dispersion is modelled with the drag model of Schiller and Naumann [30] and the lift model of Legendre and Magnaudet [31], whilst interphase momentum transfer between large interfaces is not modelled. Both phases were assumed to be turbulent and simulated with a wall-modelled Large Eddy Simulation (LES) approach with the standard Smagorinsky [18] closure for the subgrid stresses and a wall function [32] to relax mesh requirements at the walls.

The GEMMA approach was used to simulate the multiphase flow within the column; given the relatively coarse mesh size compared to the droplet size expected within the system (known from the experiments of [25]), it was expected that the $C_{\alpha }$ identifier will be equal to zero everywhere in the domain, effectively reducing GEMMA to a standard dispersed multifluid formulation. Therefore, the logical switching of $C_{\alpha }$ was disabled. The organic phase droplet size was evaluated using the reduced OPOSPM population balance described above.

The compartment-averaged organic phase holdup and Sauter mean diameter obtained along the column are compared with the experimental measurements of [25] in Figure 4. Overall, the CFD model is in good agreement with the experiments, with an average error of 9.2% for the organic phase holdup and 6.5% for the Sauter mean diameter, suggesting that a good representation of the hydrodynamic behaviour of the system has been obtained.

The continuous phase inlet solute concentration is 0.962 × 10³ mol m⁻³ in the experimental set-up, and the dispersed phase solute concentration is 0.171 × 10³ mol m⁻³. The resulting experimental solute concentration profiles along the column are shown in Figure 4. In the main body of the column, situated between the inlets, the solute concentration profiles can be reasonably approximated with a linear trend and are shaded in grey. This section of the column, between 1.45 and 3.2 m, was modelled using 1 CSTR per disk, and, in total, 59 CSTR-in-series were used, with organic flowing up and the aqueous phase flowing down. Once this simple approach had been proved in principle, additional complexity can be added when building the reduced-order model used to represent the system within the mass transfer model. The section associated with the disk near the measuring point located at the height of 2.6 m was used to back-calculate the overall mass transfer coefficient on the organic phase side via Equation (16), using the experimental values available for droplet size, holdup, $C_d$ and $C_c$; the resulting experimental $k_{ov}$ value is equal to 1.99 × 10⁻⁵ ms⁻¹. This value is used below as a benchmark for the $k_{ov}$ values calculated with the hydrodynamic parameters obtained from the CFD results and using different correlations available in the literature for the phase-specific mass transfer coefficients $k_d$ and $k_c$. A review of correlations used for the prediction of multiphase mass transfer coefficients was published by Attarakiha et al. [26] and is summarised in

Table 1. Correlations used to evaluate the phase-specific mass transfer coefficients (summarised from [26]) employed in estimating $k_c$ and $k_d$ for: $d_{32}$ = 2.9 × 10⁻³ m, $U_s$ = 0.025 ms⁻¹, $\alpha =7\%$.

Phase	Reference	Correlation	Value (ms−1)
Continuous	Ranz and Marshall [33]	$k_c=\frac{D_c}{d_{32}}\left(2+0.6R{e}_{drop}^{\frac{1}{2}}S{c}_c^{\frac{1}{3}}\right)$	1.93 × 10−5
	Treybal [34]	$k_c=\frac{D_c}{d_{32}}\left(0.725R{e}_{drop}^{0.57}S{c}_c^{0.42}\left(1-{\alpha }_d\right)\right)$	5.09 × 10−5
	Heertjes et al. [35]	$k_c=0.83\sqrt{\frac{D_c{U}_r}{d_{32}}}$	8.05 × 10−5
	Kronig and Brink [36]	$k_c=\frac{D_c}{d_{32}}\left(0.6\sqrt{R{e}_{drop}S{c}_c}\right)$	5.82 × 10−5
Dispersed	Handlos and Baron [37]	$k_d=0.00375\frac{U_d}{1+\frac{{\mu }_d}{{\mu }_c}}$	1.38 × 10−4
	Laddha and Degaleesan [38]	$k_d=0.023\frac{U_r}{S{c}_d^{0.5}}$	3.25 × 10−5
	Pilhofer and Mowes [39]	$k_d=0.002\frac{U_r}{1+\frac{{\mu }_d}{{\mu }_c}}$	3.17 × 10−5

.

As shown in Table 1, overall, it was observed that all the correlations give comparable predictions for both phases; the only exception is the Handlos and Baron [37] correlation, which overestimates $k_d$ markedly with respect to the other correlations. It was decided to use the correlation of Treybal [34] for $k_c$ since this was developed to be used in swarms of droplets rather than single droplets; the correlation of Laddha and Degaleesan [36] has been used for $k_d$ since this was derived from penetration theory, as opposed to the empirical nature of the Pilhofer and Mewes [39] correlation.

The overall mass transfer coefficient evaluated using the selected correlations, informed by the hydrodynamic parameters obtained from the CFD simulation, is equal to 2.12 × 10⁻⁵ ms⁻¹, which corresponds to a relative error of 6.19% with respect to the reference experimental value of 1.99 × 10⁻⁵ ms⁻¹. Such a relatively low error is deemed satisfactory considering the simplifying assumptions taken in the evaluation of $k_{ov}$.

To validate the mass transfer modelling approach, the mid-section of the column shaded in grey in Figure 5 was modelled as a series of CSTRs, each representing a stage. This region was selected as it eliminates the requirement for modelling mass transfer in the column inlets and separators, which cannot be approximated as a series of CSTRs. The inlet boundary conditions for the first CSTR were given by the experimental values at $x$ = 1.5 m. The CSTRs-in-series model used the predicted $\alpha$, $d_{32}$ and $U_r$, averaged over the volume of the stage, to calculate the mass transfer coefficient using Equation (16). The predicted solute concentrations in the dispersed and continuous phases are shown in Figure 5, with good agreement with the experimental results. The error is calculated to be 1.1 %, which provides confidence in the capability of the proposed modelling approach to predict the overall mass transfer performance in addition to the hydrodynamics in intensified liquid–liquid extraction devices.

3.2. Simulation of a Pulsed Sieve-Plate Extraction Column

Using the data presented by Garthe in [25], a PSEC was simulated. The column has an internal diameter of 0.08 m and contains 26 plates spaced 0.1 m apart. Plates are 0.001 m thick, with 0.002 m holes and a fractional free area of 0.2. The column is operated at the same flow rates at the RDC using the same fluids, with a pulse amplitude of 0.008 m and a pulse frequency of 1.25 s⁻¹.

A simplified CFD model of the PSEC with the above operating conditions was created, a schematic representation of the computational domain is shown in Figure 6. The domain was modelled as a two-dimensional slice consisting of 176,809 predominantly hexahedral cells with refinement towards the column walls and plates. A velocity inlet condition was used for the aqueous inlet, aqueous outlet, organic inlet, and pulse leg. Flowrates were scaled based on the difference in column cross-sectional free area. The organic outlet was modelled using the pressure inletOulet condition. A $N_{drop}$ corresponding to an inlet Sauter diameter of 2 × 10⁻³ m for toluene was specified at the organic inlet, with a zero gradient condition used at the top. Interphase momentum transfer of dispersion is modelled with the drag model of Schiller and Naumann [30], whilst interphase momentum transfer between large interfaces is not modelled. Both phases were assumed to be turbulent and simulated with an Unsteady Reynolds-averaged Navier–Stokes (URANS) [19] approach, with turbulence modelled using the mixture k-ε model and a wall function [32] to relax mesh requirements at the walls.

As above, the GEMMA approach was used to simulate the multiphase flow within the column. Unlike the RDC case, mesh refinement and a coalesced organic phase below the plates result in large values of IRQ, and therefore, logical switching of $C_{\alpha }$ was enabled. A critical IRQ value of 16, a dispersed phase volume fraction of 0.01 to 0.99 and a $\Gamma$ of 4 was used to ensure $C_{\alpha }$ is equal to a value of 1 in the near-plate region whilst being equal to 0 in the inter-plate region. The organic phase droplet size is evaluated using the reduced OPOSPM population balance described above.

The stage-averaged organic phase holdup and Sauter mean diameter along the column are compared with the experimental measurements of [25] in Figure 7. An error of 1.6% for the organic phase holdup and 8.0% for the Sauter mean diameter was calculated. A greater degree of variability is observed in both the experimental and modelled results compared to the RDC. This is expected due to the dynamic flow conditions observed in a PSEC, with each pulse consisting of an initial jetting phase with droplet breakage followed by coalescence at the next plate. As with the RDC, the CFD predicted values reasonably agree with the experimental observations and, therefore, can be considered a good representation of the system’s hydrodynamics. Simulating the entire column in 3D with an LES approach would help give a more accurate hydrodynamic field; however, this is not considered to be computationally feasible at present. Likewise, URANS-based predictions based on a Reynolds stress turbulence closure are likely to provide more accurate solutions than those given above and should be explored in further work. Nevertheless, good agreement with available data was obtained using the URANS approach adopted.

Figure 8 presents a block flow diagram showing the schematic representation of the continuous and dispersed phase flows within the PSEC. $Q$ represents volumetric flowrate in mol s⁻¹, $C$ represents solute concentration, $n$ represents the stage number, $N$ represents the upper stage and subscripts $c$ and $d$ denote the phase. As above, axial solute concentrations for the aqueous and organic phases were calculated based on the assumption that each stage of the PSEC is equivalent to a CSTR. Pulsation of the continuous phase was accounted for via the backflow ratio $\beta ,$ representing the degree of recirculation of the continuous phase between adjacent compartments.

Three equations were used to give the continuous phase mass balance. For the aqueous phase inlet, where $n$ =$N$, Equation (17) was used. For all stages between the aqueous phase inlet and the organic phase inlet, where 0 < $n$ < $N$, Equation (18) was used. For the organic phase inlet, Equation (19) was used:

$$Q_c C_{c,in}+Q_c \beta C_{c,N-1}=Q_c{C}_{c,N}+Q_c\beta C_{c,N}+V_{N}MT{R}_N$$

(17)

$$Q_c{C}_{c,n+1}+Q_c\beta C_{c,n+1}+Q_c\beta C_{c,n-1}=Q_c{C}_{c,n}+Q_c\beta C_{c,n}+Q_c\beta C_{c,n}+V_{n}MT{R}_n$$

(18)

$$Q_c{C}_{c,1}+Q_c\beta C_{c,1}=Q_c{C}_{c,0}+Q_c\beta C_{c,0}+V_{n}MT{R}_0$$

(19)

The following equation then gives the dispersed phase mass balance:

$$Q_d{C}_{d,n}=Q_d{C}_{d,n-1}+V_{n}MT{R}_n$$

(20)

The mass transfer rate is calculated using Equation (16). Equations (17)–(20) were solved iteratively until the axial solute concentrations were invariant with subsequent iterations. The continuous phase inlet solute concentration is 0.922 × 10⁻³ mol m³ in the experimental set-up, and the dispersed phase solute concentration is 0.131 × 10⁻³ mol m³, and the resulting axial solute concentration profile is shown in Figure 9.

The axial concentration of the organic and aqueous phases in the column section was modelled, with good results obtained over the entire column length. The model correctly predicts that the dispersed phase solute concentration increases along the length of the column, although the calculated values are slightly greater than those measured experimentally. The increase in aqueous phase solute concentration along the length of the column was also modelled, although the model slightly underpredicts in the region of the aqueous phase inlet. The error was calculated to be 24.3 % which is much larger than the error calculated in the RDC case. This could be reduced by considering mass transfer in the upper separator and further compartmentalisation of the PSEC stage, with the near-plate and inter-plate regions modelled separately, with the allocation of the volumes associated with each performed dynamically to account for the time-varying nature of PSEC operation.

4. Conclusions

This work introduces a mass transfer modelling approach informed by hydrodynamic information obtained via detailed CFD simulations performed using GEMMA. The approach can be used to evaluate mass transfer driven by concentration gradients and relies on a reduced-order representation of the device of interest. The following key conclusions can be drawn from this manuscript:

• It has been shown that the proposed modelling methodology represents the mass transfer performance in liquid–liquid extraction columns and provides a pragmatic tool for evaluating complex multiphase flows.
• The hydrodynamic modelling is in good agreement with the experimental observations for stage-averaged dispersed phase holdup and Sauter mean droplet diameter in an RDC and a PSEC.
• Furthermore, the modelled axial continuous and dispersed phase solute concentrations in the mid-section of the RDC are in good agreement with experimental observations when modelled as a series of CSTRs.
• Finally, the modelled axial solute concentrations in both the continuous and dispersed phases of a PSEC agree with experimental observations when modelling the column as a series of CSTRs, with pulsation accounted for via backflow in the continuous phase.

It is recommended that future work should focus on the validation of the modelling approach over a wider range of operating conditions and chemical systems, further compartmentalisation of the fluid domain, relaxation of the chemical equilibrium assumption at the interface and the implementation of local mass transfer modelling capabilities directly within the GEMMA approach.