MIXING CONTROL IN A CONTINUOUS-FLOW MICROREACTOR USING ELECTRO-OSMOTIC FLOW

In recent years, pharmaceutical production has been stimulating the gradual miniaturization of continuous-flow chemical reactors. This process eventually resulted in the emergence of a new generation of microreactors. The advantages of these new reactor types are the flexible production that allows us to quickly reconfigure the scheme, small reactant quantities used for the synthesis, the control of the main reaction parameters with high accuracy. Nevertheless, a decrease in the thickness of the channels where the species contact and react forces us to search for new non-mechanical mechanisms for mixing. This problem is relevant for the slow reaction occurring in a slot where diffusion alone cannot provide mixing at reasonable distances from the entrance. It is also true for the fast reaction that takes place in a frontal manner. In this work, we consider the efficiency of mixing the reactants induced by electro-osmotic flow in a Hele-Shaw configuration with non-uniform zeta potential distribution. As a test reaction, we take the neutralization reaction with simple albeit non-linear kinetics. The reaction occurs between two miscible solutions, which are initially separated in space and come into contact in a continuous-flow microreactor. The reaction proceeds frontally, which prevents the efficient mixing of the reactants due to diffusion. Using direct numerical simulations of 2D and 3D flows, we demonstrate that the zeta potential applied to boundaries can effectively control the mixing rate of fluids by lengthening the front of the reaction. This approach makes it possible to increase the yield of the reaction product. Mathematics Subject Classification. 76D55. Received September 30, 2020. Accepted July 31, 2021.


Introduction
Revolution in chemical technologies revolves around microreactors with the continuous flows of the reactant solutions and reaction products. The beginning of the 2000s, with the introduction of new technologies based on continuous-flow reactors, saw significant changes in fine organic synthesis. Since pharmaceutical manufacturing is looking for flexibility in synthesis readjustment more than for a larger volume of product output, miniature reactors with the specific size reaching dozens of microns, which turns these reactors into microfluidic devices, are being more and more developed [2,16,22,27,28]. Generally, continuous-flow microreactors have a small reactor zone and provide a high flow rate for a solution that maintains quick mixing of the reactants and efficient reaction behavior. formulate the problem and discuss all aspects of the proposed mathematical model. The details of the numerical method and parameter values to simulate the fluid flow in a closed reactor are given in Section 3. Here we evaluate the contribution of electrokinetic and diffusion processes on the reaction rate. We present a numerical solution to a two-dimensional problem formulated within the Hele-Shaw approximation and compare this solution with the results of direct numerical simulation of a three-dimensional problem. Section 4 looks at the control challenges for the continuous-flow microreactors with the EOF controller. Here we investigate the possibility of using dynamic EOF to obtain the optimal absolute and relative yield of the reaction product. Section 5 summarizes the results.

Mathematical formulation
Let us look at the chemical microreactor in the form of the Hele-Shaw rectangular cell (Fig. 1). It means that the fluid flow develops between two closely spaced parallel plates. We assume the spacing h to be sufficiently small to ensure that the viscous forces dominate. Figure 1 outlines a typical arrangement, where the flow can be generated either by a pressure gradient applied across the ends of the plates or by electro-osmotic flow applied to wide sidewalls of the cell. We assume that the Hele-Shaw cell constantly experiences an electric field directed along with the slot. In this case, reactant mixing is affected by the non-homogeneous distribution of the zeta potential of the electric double layer located at the solid-liquid interface. Feeding the electrodes or changes of the material (it's chemical composition) of the cell's wide borders introduce changes in zeta potentials ζ L , ζ U of the lower and upper interface, respectively. In what follows, for simplicity, we assume that the potentials are applied symmetrically: ζ L = ζ U = ζ. We will assume that the gap width between the Hele-Shaw plates is small in comparison with the characteristic inhomogeneity l of the zeta potential: The flow of the incompressible reacting fluid is governed by set of reaction-diffusion-convection equations coupling continuity equation (2.2) and Navier-Stokes equation (2.3) to evolution equation (2.4) for the the concentrations: where U = (U x , U y , U z ) is the three-component velocity field, C = (C 1 , . . . , C n ) is the vector of concentrations of chemical compounds entering into reactions, where n stands for the number of species. The term R(C) describes the reaction kinetics and depends on the concentrations C. Physical quantities characterizing the properties of a medium include the dynamic viscosity η, the solvent density ρ 0 , and the vector of the species diffusivities D = (D 1 , . . . , D n ). In what follows, we assume the diffusion coefficients of solutes D to be constant. We have previously shown that the effect of concentration-dependent diffusion can generate new types of instabilities (see, for example, Refs. [6][7][8]). This effect manifests itself best under the action of volumetric forces (inertial forces, gravity). In this work, however, we consider the system under zero gravity, and the primary mechanism, electric double layer, that can induce the fluid motion is of a surface nature. Therefore, we neglect the concentration dependence of diffusion in this paper. Under the theory of an electric double layer, electroosmotic rate and the tangential component of the applied electric field are connected by the Helmholtz-Smoluchowski equation: Here, = 0 r is the electric permeability of the liquid, 0 = 8.85 · 10 −12 C 2 /N · m 2 is the electrostatic constant, r is the relative dielectric permeability (for water r = 80.2), E is the external electric field directed along the Hele-Shaw cell's wide borders. Equation ( which is the ratio of surface conductivity σ s to bulk conductivity σ b . The resulting electric field and concentration fields in the case of moderate and high values of the Dukhin number cannot be considered homogeneous in the volume, which leads to a violation of the Helmholtz-Smoluchowski relation (2.5). Besides, in this case, one has to take into account the chemiosmotic mechanism [19]. In this work, we limit ourselves to considering the case of small Dukhin numbers, i.e. the case when the surface conductivity is negligible. The solute mass flux across the solid-liquid boundaries is assumed to be zero: In order to render the equations dimensionless, it is necessary to choose characteristic scales for variables. To adimensionalize the problem (2.2)-(2.5), (2.7), we follow [3] in choosing the following characteristic scales: Here, ζ 0 is the specific value of the zeta potential, E 0 is the reference value of the external field intensity, U * is the characteristic velocity defined by the Helmholtz-Smoluchowski equation (2.5) and U || = (U x , U y ) is the velocity field in the plane of the Hele-Shaw cell. The viscous scale for P * anticipates a lubrication-theory analysis presented below. The scale C * usually arises from the initial concentrations of species. In what follows, we keep the same symbols for dimensionless variables as for dimensional ones.
The dimensionless equations governing the velocity, pressure, and concentration fields in the bulk are the Navier-Stokes equations and the set of mass transfer equations for species, which read where ∇ || = (∂/∂x, ∂/∂y) stands for the two-dimensional vector differential operator. The dimensionless parameters Re, P e, and Da appearing in equations (2.9)-(2.12) are Reynolds, Péclet, and Damköhler numbers, respectively (see below for specific definition of these parameters).
We further assume that the parameter ε defined by (2.1) is small enough so that the fluid flow may be considered as quasi-two-dimensional, i.e. a Hele-Shaw approximation is applicable [5,32]. In addition, the inertial terms on the left-hand side of the Navier-Stokes equations (2.10) and (2.11) can be neglected due to the high hydrodynamic resistance of the wide plates of the Hele-Shaw cell.
We apply the lubrication theory [20] and seek solutions for U(x, y, z), P (x, y, z) and C(x, y, z) in the form of expansions in powers of ε: (2.14) As mentioned above, we will neglect convective terms of the Navier-Stokes equations in our analysis. To the order of approximation we will consider, this requires The order of the Péclet number define the order of terms in equation (2.12): diffusion plays a key role at P e ≤ ε 2 , convective terms prevail at P e ≥ 1. For generality, consider the case when both diffusion and convective terms are present in this equation (2.16) Here we assume that the size of the chemohydrodynamic structures is much larger than both the gap width h and the size l of the zeta potential inhomogeneity. This makes it possible to neglect the reaction term in the leading-order problem.
Then we obtain the leading-order problem in the following form: By integrating equations (2.17)-(2.20) and matching the results with boundary conditions (2.13) we obtain: The principal solution (2.22) for U 0 || was first obtained in [3]. Evolution equations in the Hele-Shaw approximation are then obtained by averaging with respect to the z-space direction perpendicular to the solid plates as: This procedure reduces the system geometry to a 2D domain spanned by the x and y-coordinates, in which the motion equations now read: are the fluid velocity and pressure averaged across the layer, respectively. In (2.26) and (2.27), we omitted the index of the del operator.
To define the concentration fields, we have to consider equation (2.12) for higher orders of the expansion. At two next orders, we have By successively solving equations (2.29), (2.30) with the boundary conditions (2.7) and then averaging across the Hele-Shaw gap, we obtain the following set of evolution equations for species: In what follows, we assume that two miscible liquids fill a Hele-Shaw cell. Let one part of the layer defined by y > 0 is filled with an aqueous solution of acid A, and the other part (y < 0) is filled with an aqueous solution of a base B. Then y = 0 determines the initial contact plane between the reacting species. Right after the process starts, the mixtures of acid and base diffuse into each other and react according to the scheme A + B → S with the formation of salt S at the rate k. Generally, this reaction is exothermic, but, for simplicity, we neglect this effect.
Let us reformulate the problem in terms of stream function defined as u = (∂ψ/∂y, −∂ψ/∂x). We can also write equations (2.31) explicitly. Then we obtain the following dimesionless set of governing equations: Equations (2.9)-(2.12) and (2.33)-(2.36) contain the following list of dimensionless parameters which are ratios of diffusion coefficients for the reactants and product, the Péclet number, the Damköhler number, and Reynolds number respectively. The Damköhler number is defined as the ratio of the chemical reaction rate to the advective mass transfer rate. It can be also considered as the ratio of the characteristic diffusion and reaction times: Similarly, the Péclet number is defined through the ratio of the characteristic diffusion and advection times: (2.39) 3. Batch reactor: contribution of diffusion and EOF to mixing

2D numerical solution
In this section, we evaluate the impact of diffusion and EOF on the reactants mixing in a closed chemical reactor. We consider the two-dimensional domain defined by −L/2 x L/2, −H/2 y H/2 (Fig. 1). Let us assume the walls of the reactor to be solid and impermeable to solutes: where n is the unit vector perpendicular to the boundary. Initially, the aqueous solutions of acid and base are separated in space and occupy the opposing halves of the cell: For simplicity, we assume in (3.2) that the initial concentrations of solutions are always equal. Table 1 presents the values of the dimensional parameters used in numerical simulations.
In what follows, we fix the amplitude of zeta potential ζ * = 0.1 V and the value of dimensionless gap ε = 0.1. The density of the external electric field is axially directed y and takes one of three values E * = 10 V/m, 10 2 V/m, 10 3 V/m. In Table 2, we evaluate the values for the dimensionless parameters for these three cases. Other parameters have been fixed to the following values: D b = 0.68, D s = 0.5. In this work, we use the tabular values of the diffusion coefficients for aqueous solutions of nitric acid HNO 3 , sodium hydroxide NaOH, and their salt NaNO 3 , being the standard and most common reactants for the neutralization reaction.
We have applied the finite-difference method to solve numerically the defined evolution problem (2.33)-(2.36) supplemented by boundary and initial conditions (3.1), (3.2) for the variables ψ, A, B, S. Differential equations were approximated by the standard finite difference equations with one-sided differences for the time derivatives and central differences for the spatial derivatives. Obtained finite-difference equations were solved with an explicit scheme. The main calculations were carried out on uniform grids with a step along the coordinate of 0.01. For example, a 500×500 node square cell was used for a 5×5 reactor.
The important measure to evaluate the role of EOF and diffusion in mixing is given by the mixing rate R(t) computed as the number of points where the concentration of salt S(t, x, z) is larger than an arbitrary threshold S 0 . In other words, we compute a time-dependent quantity where θ(x, y, t) is a Heaviside function, S 0 is the small threshold constant (it is set to be 0.001). The integral (3.3) is normalized by the area HL. This value characterizes the degree of mixing in the given miscible system since it shows how well the reaction product S has been redistributed over the two-dimensional domain.
As an example, consider the following steady distribution of the zeta potential: ζ = sin(2πx). In the case of a constant zeta potential, the velocity field, as a rule, quickly establishes and remains unchanged throughout the simulation time (Fig. 2). The concentration fields of the reactants and the reaction product continue to evolve under the fluid flow. Since the neutralization reaction is fast and frontal, the strategy for increasing the product yield is to lengthen the reaction front as much as possible. The velocity field caused by the spatially periodic zeta potential is directed perpendicular to the reaction front (Fig. 2), which leads to a gradual development of the deformation of the reaction front in the transverse direction. With an increase in the front length, the reaction involves more and more molecules, which leads to an intensification of the reaction-diffusion processes. Figure 3 shows the salt concentration field at t = 0.5, 4, 10 for P e = 30. The effectiveness of this approach is illustrated in Figure 4, which presents the variation of the spatial reaction rate R(t) with time for the pure diffusive mode only (E = 0) and the mode with electro-osmotic mixing. In the first case, we can see a slow involvement of molecules in the reaction due to their diffusion towards the reaction front. Since all possible mechanisms associated with the action of volumetric forces are disabled, one can observe the processes of pure reaction-diffusion. In the latter case, the state of complete mixing is reached already by the time t = 5. The reaction rate R(t) in the case of the switched-off zeta potential for the same time is five times less.
At small values Péclet number, diffusion plays the key role in reactants mixing, which maintains the reaction running. Figure 5 shows the salt distribution at times t = 0.5, 2 for P e = 0.3. Figures 4b and 4c give the comparison of the spatial reaction rate R between the pure diffusive mode (E = 0) and mode with moderate (P e = 3) and weak (P e = 0.3) electro-osmotic mixing, respectively. Diagrams clearly show that diffusion plays the main role in reactants mixing in the latter case.
Thus, we can state that the microreactor can successfully mix the reacting solutions under the diffusion and EOF competition. The governing equations define the ratio between the diffusive and electro-osmotic mixing via the Péclet number (2.39). The diffusion coefficients of acid and base dissolved in water are both of the order 10 −9 m 2 /s. Therefore, we can conclude that the diffusion is the primary mixing mechanism for P e < 1, competes with the EOF mixing for P e ∼ 1, and plays a minor role at larger values of the Péclet number. The Péclet number increases with the magnitude of the external electric field and zeta potential. We can adjust these two values and achieve efficient electrokinetic mixing at quite suitable sizes of the Hele-Shaw reactor. We found that a cell of 10 × 10 mm and a width of 1 mm is the most efficient among the parameters studied. Cells of smaller size decrease the Péclet number value and, thus, the diffusion impact grows; it turns into the primary mixing mechanism.

3D numerical solution
To verify the presented results obtained within 2D simulation, we performed a numerical analysis of the 3D evolution problem governed by equations (2.2)-(2.4) with the neutralization reaction kinetics. The 3D problem    10 × 10 × 0.2 mm. At the very beginning of the evolution, the opposite halves of the slot were filled with aqueous mixtures of reactants separated by the interface of y = 0 (Fig. 1). Table 1 provides the parameter values used in simulations. These values were selected according to the results of experimental studies [6][7][8][9]. The amplitude of the external electric field was fixed to 1000 V/m. Figure 3d-f shows the density plot of the salt concentration in the z = 0 plane developing under the action of static zeta potential. The evolution of the pattern is illustrated by three frames derived for three consecutive times t = 0.5, 4, 10. Comparison of the results of 2D (Fig. 3a-c) and 3D (Fig. 3d-f) simulations demonstrate a perfect agreement between the two approaches. Figure 6a shows the profile of the z-component of the velocity numerically calculated for the point x = 0.25, y = 0 and indicated by the open circles. We can compare this result with the same profile obtained analytically. Indeed, if we assume U x = 0, U z = 0 and ∂/∂y = 0, then the continuity equation (2.2) is satisfied identically and the Navier-Stokes equation (2.3) takes the form of the Laplace equation: Taking into account the boundary condition (2.5) and the given form of the zeta potential ς = ς * sin (k x x) we obtain: This solution is indicated in Figure 6a by the solid line. One can see that there is good agreement between numerical analysis and analytics. As the cell gap width decreases, we get an almost rectangular velocity profile. Figure 6b shows the variation in the reactant concentrations across the layer at the point x = 0.25, y = 0 obtained numerically. The concentrations are practically independent of z since the mass transfer of solutes across the layer is practically absent.

Controlling a continuous-flow microreactor using EOF mechanism
Let us look now at a continuous-flow chemical reactor, with the reactor zone being a 5 × 5 square Hele-Shaw cell. We assume that the aqueous solutions of acid and base have the same A 0 concentration and occupy the opposing halves of the cell at the very beginning. In this case, the reactants are fed into the cell with the same concentrations equal to A 0 = B 0 . At the same time, the product is removed through the opposite sidewall.
So, boundary conditions at the solid walls are: where Q is the liquid flux at the reactor entry. Boundary conditions for the output concentrations could be described with the Dankverts conditions [13]. The assumption is that a convective flow towards the exit channel is quite large, while the diffusive transfer in the same direction is small. Therefore, it is reasonable to say that the substance flow at the border perpendicular to the channel axis equals zero.
Boundary condition for the stream function is defined as follows: Overall reaction product output through the exit channel is the main reference value in the problem. This output can be evaluated by calculating an instantaneous solvent flux eliminating the salt from the reactor zone: To define a relative reaction product output, the instantaneous fluxes of the reactant's solutions through the exit channel are calculated: The process can be non-stationary, therefore it is reasonable to average the flow rates (4.4) and (4.5) over time: The same finite-difference approach was applied to solve numerically the evolutionary problem (2.33)-(2.36), First, let us refer to the results for the stationary distribution of zeta potential. Let us take, for example, ζ = sin (−k x x), where k x = 20π/L. Figure 7 gives the reaction product field at the time t = 5 calculated at Q = 0 (left) and Q = 5 (right). These findings show that the reactants are mixed slower in a continuous-flow reactor with the stationary distribution of zeta potential. In this case, zeta potential non-homogeneities move concerning the disturbances of the wavefront, which they create. Therefore, disturbances stop their growth. To enhance the EOF impact, one needs to provide the movement of zeta potential non-homogeneities at the same rate as the fluid flow. This case requires to replace a stationary potential ζ = sin (−k x x) with its wave analogue ζ = sin k x (vt − x), where v = Q/L is the phase rate, k x = 20π/L. Figure 8 shows A, B, S fields calculated for this potential at v = 5/7, Q = 25/7. Now, the figure clearly shows that the disturbances of the wavefront are accelerating. This process happens gradually with the fluid flow from the reactor entry to exit. Therefore, it shows a transition process with the time calculated as follows: Transition time t h is one of the problem parameters. It impacts the reactants' mixing rate and relative reaction product output. Figure 9 illustrates the time dependence of the product output and reactants for wave and stationary zeta potentials. Wave zeta potential provides a significantly larger reaction product output.
Numerical simulations show that transition time t h affects the absolute and relative reaction product output (Fig. 10). The absolute output diagram has a sloping extremum at the interval 3 < t h < 7. Areas t h with the optimal values for the absolute and relative output do not overlap. However, to achieve reasonable values for absolute and relative reaction product output, the parameters (t h , v, Q) can be chosen at the stage of designing  the flow reactor with EOF mixing. The working parameters to choose from can be narrowed if we refer to the condition for the specific time of the non-stationary electric field shown in [11]: to apply the Helmholtz-Smoluchowski equation to the liquid-solid body borders, the electric field should at least be quasistatic. This condition results in the requirement t h ≥ 10. Within this area of parameters, the joint action of EOF and diffusion give higher (more than 80 %) relative output and quite high absolute reaction product output.

Conclusion
In this paper, we have studied the control of the mixing process in a continuous-flow Hele-Shaw microreactor where a second-order chemical reaction occurs. To control the mixing process, we use an electroosmotic flow created by a combination of direct electric current and zeta potential. In this case, the EOF competes with diffusion. We found that the role of the first mechanism increases with a growth of the size of the microreactor, the voltage of the external electric field, and the intensity of the zeta potential. We found that 10×10 mm cell with 1 mm in gap width is the most reasonable one among the examined parameter values. Diffusion affects more in smaller cells (less than 10 −4 m): it turns into the principal mixing mechanism. EOF with a dynamic zeta potential is recommended for the flow microreactors to achieve high values for the absolute and relative reaction product output. To enhance the impact of EOF, we need to provide the movement of zeta potential non-homogeneities at the same rate as the fluid flow. The reaction front deformations created by the EOF gradually increase from the entrance to the exit from the channel. The characteristic growth rate of deformations defined by the ratio of the channel length and flow rate is a parameter that determines the most favorable operating mode of the microreactor.