Stability Analysis of Reactive Multiphase Slug Flows in Microchannels

Conducting multiphase reactions in micro-reactors is a promising strategy for intensifying chemical and biochemical processes. A major unresolved challenge is to exploit the considerable benefits offered by micro-scale operation for industrial scale throughputs by numbering-up whilst retaining the underlying advantageous flow characteristics of the single channel system in multiple parallel channels. Fabrication and installation tolerances in the individual micro-channels result in different pressure losses and, thus, a fluid maldistribution. In this work, an additional source of maldistribution, namely the flow multiplicities, which can arise in a multiphase reactive or extractive flow in otherwise identical micro-channels, was investigated. A detailed experimental and theoretical analysis of the flow stability with and without reaction for both gas-liquid and liquid-liquid slug flow has been developed. The model has been validated using the extraction of acetic acid from n-heptane with the ionic liquid 1-Ethyl-3-methylimidazolium ethyl sulfate. The results clearly demonstrate that the coupling between flow structure, the extent of reaction/extraction and pressure drop can result in multiple operating states, thus, necessitating an active measurement and control concept to ensure uniform behavior and optimal performance.


Introduction
An effective exploitation of the advantages of microscale operation for industrial production processes requires a reliable and robust numbering-up procedure.For multiphase systems, one must not only ensure uniform flow rates of each phase but also similar flow structures [1] and residence time distributions [2,3] in the individual microchannels, if the optimal performance is to be achieved [4].Furthermore, in reactive systems, multiplicities can arise, leading to drastic differences in the behavior between otherwise identical microchannels [5].The problem is compounded by the greater significance of fabrication tolerances at the microscale, since even slight discrepancies can give rise to considerable changes in volumetric flow rates, which scale with the fourth power of the channel diameter [6,7].
Passive distribution strategies thus necessitate high-precision machining, usually entail large upstream pressure drops and are unable to adapt to changing operating conditions, such as fouling [8].An active approach, in which the flow in each microchannel is monitored and regulated using simple non-invasive probes and actuators overcome these shortcomings [9].
Even so, the multiplicities in reactive flows still represent a particularly severe challenge, since they can not only diminish performance but may also lead to hazardous runaway situations or reactor failure.In polymerizations, for instance, the reaction and hydrodynamics are strongly coupled through the viscosity of the reaction medium.Should, for whatever reason, the flow in one channel be less than that in the others, the resulting increase in the residence time will yield higher conversion and thus viscosities, increasing the resistance to flow and lowering the flow rate still further.If no remedial action is taken, this self-reinforcing mechanism amplifies any flow divergence and can ultimately cause the channel concerned to plug.
Of course, such multiplicities are by no means confined to microscale systems.In multiphase flow in microchannels, on the other hand, certain reactive hydrodynamic multiplicities can arise which are unlikely or even impossible in macroscale equipment, due to the intensified interaction in the biphasic flow.In an alternating gas-liquid slug flow with absorption, for example, the gas uptake reduces the size of the gas slugs, thus, changing the pressure drop characteristics of the flow, raising the channel residence time and thus enhancing the absorption still further.Conversely, reactions that generate gas or decrease viscosity exert a stabilizing influence on the uniformity of the flow distribution over parallel microchannels.
The objective of the work presented here was to elucidate the relevance of such microscale-specific phenomena in the parallelization of multiphase flows in microchannels and to establish if it might be feasible to utilize self-regulating behavior in order to ameliorate irregularities in flow distribution.

Experimental Extraction
A liquid-liquid extraction was used to validate the results obtained after adapting the gas-liquid model presented in Section 3. The chemical, analytics and the experimental set up are described in the following sections.

Chemicals and Analytics
The chemicals n-heptane and acetic acid were supplied from VWR prolab with a purity of 99.7% w/w and 99% w/w.The ionic liquid 1-ethyl-3-methylimidazolium ethylsulfate [EMIM] [EtSO 4 ] was obtained from Io-Li-Tec (Ionic Liquids Technology, Tuscsaloosa, AL, USA) with a purity of 99% w/w.The most significant properties of the chemicals used are shown in Table 1.The viscosities were measured by means of a capillary viscosimeter (KPG Ubbelohde capillary viscosimeter, Schott AG, Mainz, Germany) and the concentration of acetic acid in n-heptane was analyzed by an Agilent 7890A gas chromatograph using an autosampler Agilent 7693 and a flame-ionization detector (FID) from Agilent Technologies (Agilent Technologies, Böblingen, Germany).A capillary column HP-05 (Agilent Technologies, Böblingen, Germany) with a length of 30 m and an inner diameter of 320 µm was employed and helium used as carrier gas with a flow rate of 1.22 mL/min at a temperature of 250 °C.

Experimental Setup
The experimental set-up is shown in Figure 1.n-Heptane containing 10% w/w acetic acid and EMIM EtSO 4 were introduced by two syringe pumps (LA-100 Landgraf Laborsysteme HLL GmbH, Langenhagen, Germany) into a "Mixer-Valve" where the slug flow of the immiscible phases was generated.This special valve enables the slug length to be manipulated whilst keeping all other parameters, for example flow rates and the physical properties of components, constant.The outlet of the mixer was connected to a Polytetrafluoroethylene-capillary (PTFE-capillary) with an inner diameter of 800 µm.To achieve well-defined extraction times a union T-separator was installed downstream of the contacting section to implement an instantaneous liquid-liquid phase separation by virtue of their different wetting properties.One of the separator outlets consists of a steel capillary, with an inner diameter of 750 µm, and the other comprises a PTFE-capillary with an inner diameter of 800 µm.The dispersed EMIM EtSO 4 phase, which does not wet the PTFE capillary wall, leaves via the steel outlet while the continuous n-heptane phase exits through the PTFE branch.The pressure drop measurement technique is depicted in Figure 1b.The pressure drop of interest is that in the reaction capillary denoted ΔP slugflow .To avoid the pressure sensor influencing the multiphase flow structure, the pressure sensors are placed upstream of the mixing point.The pressure drop ΔP 2 was measured without the separator and ΔP 1 following removal of the capillary from the mixer.The pressure drop of the mass transfer contacting section is then calculated as: In order to describe the mass transfer for a liquid-liquid system, the extraction efficiency was used: where C Eq is the equilibrium concentration at a phase ratio EMIM EtSO 4 :n-heptane of 0.2.This concentration was measured at room temperature in a 20 mL batch reactor.C in and C out are the inlet and outlet concentrations of acetic acid in the extract phase, i.e., EMIM EtSO 4 , respectively.The phase ratio (α) was defined via a snapshot with a camera and an optical sensor.For this case, the volume fraction of the film was considered negligible.The phase ratio can then be defined as: Observed flow structures and unit cells under different conditions are presented in Figure 2.

Gas-Liquid Model Equations
As shown in previous work, bubble length and bubble velocity play an important role for hydrodynamics and the mass transfer in the slug-flow regime.Even though some studies have focused on the possibility of the expansion of the bubble due to the inherent pressure drop [12,13], and in more recent studies [14,15] due to mass transfer from the liquid into the gas phase, none of them has considered the fact that, for very reactive systems, the effect of the consumption may well exceed that of the expansion of the bubble or the case where the mass transfer occurs from the gas to the liquid phase.The model presented in this work couples hydrodynamics with mass transfer and thus accounts for the change in the size of the bubble as a result of chemical absorption.
When the length of the bubble alters, several parameters are affected, starting with the superficial velocity.Since the volume of the unit cell comprising a single bubble and the associated slug is changing along the channel length, the two-phase superficial velocity will be a function of the position as given in Equation ( 4).This entails that all parameters that are velocity-dependent are also location-dependent, i.e., capillary and Reynolds numbers, film thickness and mass transfer coefficients.
In slug flow, it is known that bubbles travel faster than the nominal superficial velocity due to the presence of a thin liquid wall film surrounding the bubble.The true velocity can be approximated by a simple mass balance in round capillaries, as expressed in Equation ( 5) (a deviation of no more than 5% was reported by [3]).In this equation, δ represents the film thickness and can be obtained by the expression found by Aussillous and Quéré [16], describing the extent of significant films in micro-channels (Equation ( 6)).As can be seen, the film thickness is a function of the capillary number, = μ σ which in turn is a function of the velocity and therefore also of the location as stated earlier.
To describe the hydrodynamics, the pressure drop must be defined.In this work, the pressure drop is taken according to the unit cell based model presented by Warnier et al. [13].In this model the pressure drop across a unit cell (presented in Figure 3) is defined by the individual contributions of the liquid slug and the pressure drop generated by the presence of the bubble.The individual bubble pressure is described by the corrected Bretherton's approximation, and the liquid slug pressure drop is described by the fully developed Hagen-Poiseuille with a correction for the amount of liquid in the domed bubble cap region (ζ).
( ) ( ) ( ) The differential pressure drop can then be obtained as: However, the length of the bubble is not constant along the channel length.In the case of gas absorption into the liquid phase, for instance, it is necessary to couple the mass transfer with the pressure drop in order to determine how the bubble contracts along the channel length.
The change in the number of moles for the absorbed component with respect to the residence time can be given by: since the number of moles of the inert is not changing with time Equation ( 12) becomes: yielding: The residence time for the bubble is given by the bubble velocity yielding: The chemical system used for this study was the widely studied chemical absorption of CO 2 in NaOH [17].First a pseudo 1st order reaction was assumed, but since the reaction was found to be fast enough, the chemical absorption was then simply accounted for by the inclusion of the enhancement factor.Finally, the change in the molar flow due to chemical absorption can be given by: As for determining the mass transfer coefficient, several studies conducted on mass transfer in gas-liquid slug flow [18][19][20][21] were considered.In this work most of the available models were examined, with particular emphasis on the models of Vandu et al. [18] and Yue et al. [19], as they allow the introduction of changes in the parameters along the length of the channel.Vandu's model treats the cap and wall film contributions separately according to: Yue's model, on the other hand, considers the overall mass transfer across the entire interface: Finally, the link between the mass transfer and the pressure drop is given by the bubble length, since a relationship between bubble length and the total number of moles can be expressed as: the total number of moles is linked to the molar fraction by the relationship: yielding the influence of all parameters on the bubble length: It is important to notice that the expansion due to the pressure drop presented in previous works is included in this model as well, since the volume of the bubble is calculated by using the pressure at each location.

Simulations
For the simulations, as mentioned before, the absorption of carbon dioxide (CO 2 ) into sodium hydroxide (NaOH) was studied.For the modeling the resistance to mass transfer in the gas phase was neglected.The diffusion and Henry coefficients as well as the surface tension were taken for the system H 2 O/CO 2 .The system was considered to be isothermal at room temperature (298 K).The parameter values are summarized in Table 2.In order to validate the model, the results obtained were compared with those obtained by previous simulation and experimental studies.First of all, the effect of the mass transfer model on the hydrodynamics was evaluated.The parameters used in this study are presented in Table 3 and the results given in Figure 4.The results show that the choice of mass transfer model shows little or no influence on the hydrodynamics, only Vandu's Model yields slightly different behavior, due to the overestimation of the mass transfer coefficients and hence a faster contraction of the bubble.It can also be seen that the range of pressure drop per unit length is of the same order of magnitude as in previous studies.
The model was also used to calculate the mass transfer coefficients and good agreement was obtained between the simulations and the experimental values presented by Yue et al. [19] and Sobieszuk et al. [20] (Figure 5).Since this model allows for the change in velocity, a more realistic residence time is calculated for the bubble and hence the expansion behavior present in the experiments is also captured.For the final validation, the pressure drop model considered was compared to those previously presented in the literature (Figure 6).It can be seen that the unit cell based pressure drop model yields lower values for the pressure drop in comparison to Kreutzer's model [22].However, the results are in agreement with the values presented by Warnier et al. [13].[13] and the Kreutzer model is used as presented in [22].

Hydrodynamic Multiplicity
A hydrodynamic multiplicity has already been found in vertically operated monoliths by Kreutzer et al. [3].However, it was due to the hydrostatic pressure in comparison to the coupling of parameters presented in this study.
Once the model had been validated, the stability of the system was studied.Since the contribution to the pressure drop in the gas phase is small for long slugs, a short slug was taken as the basis for this study.The pressure drop was then calculated by varying the bubble length; the length of the reactor was taken as 10 cm and its diameter as 800 µm.For simplicity the calculations were first performed neglecting the film.
For short bubbles (Figure 7) the pressure rises as expected at a higher velocity.However, upon increasing the bubble length, this effect tends to become smaller indicating to the hydrodynamic instability sought.When bubble size was increased further, the multiplicity effect due to the increased pressure of the number of interfaces generated for higher conversions was found as illustrated in Figure 8.This effect can be explained as follows: at a low velocity, the bubble will spend more time inside the channel therefore the conversion will be high, resulting in an increase of the number of interfaces and the average viscosity of the unit cells.At a higher velocity, conversion will be lower leading to a system with a lower average viscosity a lower number of interfaces and hence an equal pressure drop to the case with lower initial velocity.This effect is depicted in Figure 9.

Film Effect
Since the results were obtained for the model without film, the next step was to calculate the complete model to check whether the discrepancy due to this simplification was significant.The results are presented in Figure 10.It can be seen that the inclusion of the film has a slight influence on the behavior, especially in the regions of higher superficial velocity.This can be explained by the wall film being thicker at higher velocities film since it depends on the capillary number.However, this effect was considered to be not large enough to affect the results qualitatively and therefore the subsequent calculations were performed without considering the film, since the computational effort was considerably smaller (with the film an iteration is required for each single step in order to establish the bubble velocity).

Interfacial Geometry
It has been found in some studies that the interfacial curvature factor used to calculate the interfacial pressure drop can change as a function of the velocity [23].The value of 7.16 normally used (as in Equation ( 6) was proposed initially by Bretherton for hemispherical caps.Jovanovic et al. found that these values can actually range between 3 and 10.Since the instabilities arise from the bubble pressure drop, as shown previously, this effect also had to be considered.To account for this phenomenon, the values from Jovanovic et al. were correlated using a simple linear regression (Equation ( 26)).The results are presented in Figure 11.It can be seen that the effect of the curvature exerts a noticeable effect on the pressure drop characteristics, making the region of the multiplicity more pronounced.The higher curvature for lower velocities increases the influence of the remaining caps at the end of the channel where the lowest velocities are found, and as explained earlier, the effect of these extra interfaces is the source of the multiplicity.

Effect of Viscosity Changes
In the Section 4.1.2,a hydrodynamic multiplicity was identified, albeit for systems with unrealistic bubble lengths.A factor playing a major role in this phenomenon was the average viscosity in the unit cell.It is, thus, intuitively obvious to imagine that if the average viscosity changes along the channel length, the hydrodynamic multiplicity could be present in a wider range of systems in comparison with the case presented in Section 4.1.2.One possible case where this behavior can occur is in the case where viscosity changes with the progress of a reaction/extraction/absorption.For such systems, a general formula for the viscosity as a function of the conversion could be represented as: where a would represent the net increase in viscosity, and b the mechanism.Using this behavior, multiplicities can be found even for systems with short bubbles.Figure 12 depicts the results for an initial bubble five times larger than the channel diameter and a value for a of 4 for various values of b.
Figure 13 presents the results for an initial bubble five times larger than the channel diameter and a value for b of 3 for different values of a.
As can be seen in both figures, any change increase viscosity leads to major changes in the pressure drop characteristics of the system and yielding instances of hydrodynamic multiplicity.
This leads to the question as to whether for systems with decreasing viscosity, the behavior would be the opposite, i.e., result in more stable systems with more uniform distribution.These results are presented in the same fashion as for increasing viscosity: first a fixed value of b of 3, and then for a fixed value of a of −0.5 in Figure 14.As it can be seen, no signs of instability were detected for this system, which exhibits self-regulating behavior.Experiments with a liquid-liquid extraction system in which the solvent viscosity diminishes with solute concentration were performed in an attempt to try and validate these results.4.1.6.Model Adaptation to the Liquid-Liquid Case Experiments were performed with a liquid-liquid system due to its relative simplicity in comparison to the gas-liquid case.The gas-liquid model was therefore adapted to an extractive aqueous-ionic liquid system.Conceptually the modifications were: • The pressure drop model was changed to that of a stagnant film model as presented by Jovanovic et al. [23].
• The extraction phenomenon (extraction of acetic acid from n-heptane into an ionic liquid) replaces the chemical absorption (absorption of CO 2 in NaOH).The equation representing the concentration change along the channel is the following: With this adaptation, it is implied that the unit cell length will remain unchanged for both of the two phases and therefore no change in velocity will occur.However, the coupling of the mass transfer presented in the gas-liquid system is still present.Additionally, the viscosity change presented in Section 4.1.5will be supplied by the change in viscosity of the ionic liquid.The experimental results are presented in the following section.
The hydrodynamic parameters for this system can be found in Table 1.

Experimental Validation: Liquid-Liquid System
The adapted model describing the coupling between fluid dynamics and the change in viscosity through mass transfer or reaction via the pressure drop and the residence time was validated through the extraction of acetic acid from n-heptane into the ionic liquid EMIM EtSO 4 .In this system, the viscosity of the ionic liquid changes considerably as a result of mass transfer, although the underlying behavior dictating the pressure drop interactions are the same whether mass transfer, reaction or gas slug contraction is involved.Two sets of experiments were carried out: one at a constant micro-channel length whilst varying the velocity and the other at constant velocity, varying the channel length.The pressure drop, the viscosity and the extraction efficiency were measured at each operating point as described in Section 2.2 and the results illustrated in Figure 15.

Validation
In order to validate the model for the liquid-liquid system, simulations were carried out for different configurations.The parameters used are listed in Table 4.  2)).The viscosity can now be described as: The model parameter a = −0.67 and b = 0.3 of the viscosity equation (Equation (30)) were taken from experimental results by fitting the dynamic viscosity of the ionic liquid at different levels of extraction efficiency (as presented in Figure 15) and performing a regression minimizing the least squares error (LSE), these parameter yielded a LSE of 0.001.The results from this procedure are presented in Figure 16.A comparison between experimental and model results shown in Figure 17 indicates the coupling between pressure drop, change in viscosity due to reaction or mass transfer and hydrodynamics is described in a correct way.It also can be seen that the pressure drop of the experiments fitted well with the prediction of the model.The mean relative error between model and experiment for various flow rates is about 13% and for changing capillary length 11% and therefore simulations results can also be taken as valid.

Conclusions
In this paper, a model has been developed to describe the behavior of a multiphase slug flow in circular microchannels.The model couples the hydrodynamics, reaction and mass transfer phenomena, and for gas-liquid systems, it accounts for both the contraction of the gas phase due to reactive absorption and the inherent change in volume linked to the pressure drop.
The model was first validated by comparing the results obtained with results from other researchers yielding results of the same order of magnitude.For the pressure drop it gave results consistent with those obtained by Warnier et al. [13], and to a lesser extent, Kreutzer et al. [22].For the mass transfer characteristics, the model predicted mass transfer coefficient values in close agreement with those from the most recent publications, like the experimental results of Sobieszuk et al. [20] and Yue et al. [19].Furthermore, the model was tested to get an idea of how large the effect of the film surrounding the bubble is on the hydrodynamics, showing only a slight deviation from the case without film.Finally, the effect of the slug geometry was studied, as it was found that this could play an important role for instabilities, and the work of Jovanovic et al. [23] has shown that it was not constant but a function of the velocity.
The model was used to analyze the stability of the flow for a gas-liquid system establishing that multiplicities can arise in particular cases, such as long bubbles, or for systems with changes in viscosity.These results led to the study of systems with decreasing viscosity, since this can enable one to achieve a better distribution in parallelized microchannels.Experimental results obtained for an aqueous-ionic liquid-system found a good agreement between the experimental and predicted values (11% and 13% error, depending on the parameter varied).
The results demonstrate that the coupling between flow structure, the extent of reaction/extraction and pressure drop can result in multiple operating states, thus, making an active measurement and control strategy to ensure uniform behavior and optimal performance essential.

Figure 2 .
Figure 2. Observed flow structure and unit cell for a superficial velocity of 0.044 m/s and (a) phase ratio 1; (b) phase ratio 2.

Figure 3 .
Figure 3. Unit cell representation.L b represents the length of the bubble, L s represents the slug length R b is the bubble diameter and ζ accounts for the liquid in the cap region.

Table 3 .Figure 4 .
Figure 4. Pressure drop for different mass transfer models.

Figure 6 .
Figure 6.Pressure drop models comparison.The unit cell model is based in Warnier et al.[13] and the Kreutzer model is used as presented in[22].

Figure 7 .
Figure 7. Pressure drop for small bubbles.

Figure 8 .
Figure 8. Pressure drop for long bubbles.

Figure 10 .
Figure 10.Film effect on pressure drop.

Figure 11 .
Figure 11.Effect of the curvature.

Figure 14 .
Figure 14.Pressure drop for systems with decreasing viscosity (a) Fixed a value (b) Fixed b value.

Figure 15 .
Figure 15.(a) Pressure drop, extraction efficiency and dynamic viscosity for different superficial velocities in slug flow regime.The capillary length was 100 cm, the phase ratio 0.2 and the number of unit cells 80; (b) Pressure drop, extraction efficiency and dynamic viscosity for different capillary lengths in slug flow regime.The superficial velocity was 0.022 m/s, the phase ratio 0.2 and the number of unit cells 80.

Figure 16 .
Figure 16.Measured values for the dynamic viscosity as a function of extraction efficiency and the fitting based on Equation (30).

Figure 17 .
Figure 17.Comparison between pressure drop calculated by the model and the experimental values obtained for different velocities and capillary lengths.(a) Pressure drop as a function of velocity (fixed capillary length 1 m); (b) Pressure drop as a function of capillary length (fixed superficial velocity 0.022 m/s).

Table 1 .
Physical properties of extractive system at 25 °C and 1 atm.

Table 2 .
Parameter values for simulations.

Table 4 .
Parameter values for liquid-liquid simulations.
Additionally, Equation (27) was modified to replace the conversion with the extraction efficiency h (defined in Equation (