Stéphane Glockner

Antoine Michael

Diego Jost

Arnaud Erriguible
email: erriguible@enscbp.fr

A M D Jost

Advanced petascale simulations of the scaling up of mixing limited flow processes for materials synthesis

Introduction

The constant evolution of supercomputer performance over the past 20 years has significantly changed the scientific landscape and the way problems are approached. While algorithms and computational codes have long been ahead of the available computing power, the paradigm has recently been reversed as the available power of petaflop and soon exaflop supercomputers now requires a considerable effort in developing computational codes to take advantage of this power. Nevertheless, this opens up new horizons for better understanding, representing, and predicting physical phenomena as well as for the development of Smart Manufacturing processes, in addition to or even replacing in situ experiments. As a demonstration, we propose, for the first time, petascale simulations to answer the open-ended question of scaling up of continuous material synthesis processes. We highlight the value of using the latest supercomputers to study the scaling up without any loss of information, i.e. keeping the relevant scale of the physical description to design safe, controlled, and flexible processes. This offers a new path for the industrial development of continuous, sustainable, and scalable technology for the production of innovative materials. Today, microfluidics is considered one of the most promising technologies for generating innovative materials [START_REF] Marre | Synthesis of micro and nanostructures in microfluidic systems[END_REF] in a reproducible manner, allowing for efficient mixing, process optimization, and fast screening. This is especially true for materials synthesis for which fast and homogeneous mixing [START_REF] Bałdyga | Mixing effects on particle formation in supercritical fluids[END_REF] is the key to obtaining a narrow size distribution of particles.

However, the problem of moving towards industrialization almost always requires increasing the production and therefore the scale of the reactor. This generally yields to losses in production quality, as the advantages of hydrodynamics at small scales are lost.

The case study, chosen here to demonstrate the proposed approach, is the nanoparticle precipitation in microfluidic reactor by using supercritical carbon dioxide as antisolvent. Recent studies [START_REF] Jaouhari | Process intensification for the synthesis of ultra-small organic nanoparticles with supercritical CO2 in a microfluidic system[END_REF] have shown that a combination of the two types of intensification of supercritical fluids and microfluidic reactors (µSAS) can significantly improve the performance and reproducibility of antisolvent processes. Moreover, it has recently been proven that turbulent conditions can be achieved in high pressure microfluidic devices [START_REF] Zhang | Mixing intensification under turbulent conditions in a high pressure microreactor[END_REF], resulting in ultra-fast mixing times as small as 10 -4 -10 -5 s which are highly favorable for the synthesizing of organic nanoparticles. However, it should be noted that the mass flow rates used in the microreactors are much lower than those used for large reactors. As a result, even though the former offer excellent mixing performances and better control of the conditions the latter have the big advantage of larger production rates. Finally, continuous flow processes [START_REF] Hartman | Flow chemistry remains an opportunity for chemists and chemical engineers[END_REF] currently offer greater advantages compared to batch processes as the former can maintain constant flow operating conditions and achieve very efficient mixing. However, one of the main drawbacks of the microfluidic flow processes remains its limited adaptation to an industrial scale [START_REF] Berton | Scaling continuous API synthesis from milligram to kilogram: extending the enabling benefits of micro to the plant[END_REF][START_REF] Rossetti | Chemical reaction engineering, process design and scale-up issues at the frontier of synthesis: Flow chemistry[END_REF]. It is commonly accepted that numbering up [START_REF] Al-Rawashdeh | Phenylacetylene Hydrogenation over [Rh(NBD)(PPh3)2]BF4 Catalyst in a Numbered-Up Microchannels Reactor[END_REF][START_REF] Qiu | Numbering-up of capillary microreactors for homogeneous processes and its application in free radical polymerization[END_REF][START_REF] Zhang | Design and Scaling Up of Microchemical Systems: A Review[END_REF] is probably the best approach, compared to sizing up [START_REF] Kockmann | Scale-up concept of single-channel microreactors from process development to industrial production[END_REF][START_REF] Yamamoto | Scale-up or numbering-up of a micro plasma reactor for the carbon dioxide decomposition[END_REF], to maintain an excellent control of the operating conditions and efficient mixing. However, scaling up through numbering up requires thousands of reactors to achieve significant production thereby posing the challenging task of a uniform fluid distribution [START_REF] Hartman | Flow chemistry remains an opportunity for chemists and chemical engineers[END_REF][START_REF] Zhang | Design and Scaling Up of Microchemical Systems: A Review[END_REF][START_REF] Dong | Scale-up of Micro-and Milli-Reactors: An Overview of Strategies, Design Principles and Applications[END_REF].

One of the solutions could be to seek a compromise between sizing up and numbering up [START_REF] Enzmann | Same but different-Scale up and numbering up in electrobiotechnology and photobiotechnology[END_REF].

A first step of sizing up could be a decrease in the number of units of the system to guarantee a uniform distribution of the fluids in the system. As such, the first question that needs to be answered is: can a continuous micro-scale/milli-scale reactor operate at the same flow rates as a large scale reactor while maintaining an excellent mixing performance? To answer this question, it is imperative to propose an efficient strategy to determine the relevant size of the flow reactor that can strike a good balance between size and production ability. The classical approach, consisting in using the canonical similitude law, tends to fail in giving appropriate solutions for the continuous processes. This is inherent to the complexity of the dynamic system as many parameters interact in the continuous flow processes. Therefore, even if satisfactory results with partial or complete similarity methods can be found in literature [START_REF] Rehage | A new scale-up method for competitive chemical model reactions based on complete similarity[END_REF][START_REF] Bałdyga | Particle formation in supercritical fluids-Scale-up problem[END_REF], a simulation-based strategy appears to be a suitable solution to find a good candidate for the continuous processes. For this purpose, direct numerical simulations (DNS) -which do not include any extra modeling of turbulence since the finest dissipative scale are solved -are the most accurate way to capture the interactions between turbulence and species mixing. Although DNS has been mainly first restricted to canonical cases, the continuous development of supercomputers has expanded its applicability to the design of complex systems. DNS yields information on the flow and mixing characteristics unattainable through in situ experiments.

In this paper, it will be shown that high performance computing (HPC), through petascale numerical simulations, is a relevant and promising tool to study and analyze scale up scenarios for chemical engineering processes and to further demonstrate that continuous flow processes are competitive. The current case study is the supercritical microfluidic antisolvent process for which we propose -for the first time -high fidelity DNS simulations for micro-scale to milliscale resulting in an increase of up to 90 times the initial production flow rate. The scaling up strategy takes advantage of HPC and latest supercomputers to pursue simulations of larger reactors while maintaining the same accuracy of the physical phenomena description simulated in the micro reactors (smaller dissipative scales). As such, special care must be taken to establish the strategy used to determine the larger flow rates and reactor dimensions as well as the numerical strategy needed to run petascale simulations. Results confirm both the fact that it is possible to increase the size of the reactor whist maintaining excellent mixing conditions.

Materials and Methods

Scaling up strategy.

As previously introduced, we have chosen to consider the supercritical antisolvent microfluidic process for which the mixing performance is the key point yet to be demonstrated [START_REF] Jaouhari | Process intensification for the synthesis of ultra-small organic nanoparticles with supercritical CO2 in a microfluidic system[END_REF]. Therefore, in this article we focus on the DNS simulations of the hydrodynamic and mixing phenomena. The process consists of a simple tube-in-tube reactor as shown in figure 1. The solution, solvent, is injected in the inner tube, diameter Di, whereas the supercritical CO2, acting as antisolvent, is injected in the coflow between the diameters Do and Dext. The typical pressure and temperature of p = 100 bar and T = 40°C, respectively, are chosen and the reference case is that at microscale. The mean Reynolds number 𝑅𝑒 = 𝐷 𝑒𝑥𝑡 𝜌 ̅ 𝑢 ̅ 𝜇 ̅ is calculated inside the outer tube by considering the mean values of velocity 𝑢 ̅ [START_REF] Zhang | Investigating laminar mixing in high pressure microfluidic systems[END_REF], viscosity 𝜇̅ , and density 𝜌̅ . The mean properties are estimated according to the molar ratio of 91% between the solvent and the antisolvent. It should be noted that an important constraint of the current process is the pressure loss. Actual back pressure regulator technology is able to absorb pressure losses of a few tens of bar. The classical approach used to compute the pressure loss P in the inner tube, outer tube, and coflow zone is Δ𝑃 = It has been widely reported that the mixing is governed by the dissipation rate of the turbulent kinetic energy [START_REF] Bałdyga | Turbulent micromixing in chemical reactorsa review[END_REF]. Thanks to the Kolmogorov theory, Baldyga and Bourne [START_REF] Baldyga | A Fluid Mechanical Approach to Turbulent Mixing and Chemical Reaction Part Ii Micromixing in the Light of Turbulence Theory[END_REF] have established a direct relation between the dissipation rate and the mixing time. From a mixing point of view, the most important global parameter is the dissipation rate of the turbulent kinetic energy and, in the tube, it can be approximated through dimensional analysis as 𝜀 ≈ 𝑓 2 𝑢 3

̅̅̅̅ 𝐷 𝑒𝑥𝑡

with f being the Blasius friction factor. It should be pointed out that as the key point of the process is to obtain an excellent mixing performance, the scaling up strategy consists in conserving the global turbulent kinetic energy dissipation ratein the reactor by adjusting the inlet mass flow rates and dimensions of the tube. In literature this strategy is coined as "constant energy dissipation scale up" [START_REF] Plouffe | On the scale-up of microreactors for liquid-liquid reactions[END_REF] or "constant pressure drop approach" [START_REF] Dong | Scale-up of Micro-and Milli-Reactors: An Overview of Strategies, Design Principles and Applications[END_REF]. Henceforth, the global turbulent kinetic energy dissipation rate is referred to as the global dissipation rate.

The initial case is chosen at microscale with Di =100 m and Dext = 250 m (Fig. 1). The initial mass flow rates of the solvent (inner tube) and CO2 (coflow zone) are 400 l/min and 3000 l/min, respectively, leading to a solvent-antisolvent ratio of 91% mol. The global dissipation rate  is estimated at 235 W/kg leading to a Kolmogorov scale of about a micrometer. As the Schmidt number is about 2.6 the Batchelor scale 𝜆 𝐵 = 𝜆 𝐾 √𝑆𝑐 it is approximately in the same range as the Kolmogorov length scale. Therefore, to capture all the mixing scales and to yield accurate results the mesh size of the simulation is chosen to be on the order of the micrometer (see Table 2). As previously mentioned, to determine the optimal flow rates of the reactor needed to maintain a good mixing, we chose to consider the constant energy dissipation scale up method. We, therefore, imposed the different sizes of the reactor (see dimensions cases 1, 2, 3, 4 in table 1) and through a minimization algorithm the flow rates of the solvent and antisolvent (QEtOH and QCO2 in table 1) are optimized. The target of this optimization problem is the global dissipation rate  (table 1) which coincides to the one determined for the initial case 0. The problem is closed by adding two constraints: 1) the pressure drop in the annular/coflow zone must be less than 50 bar and 2) the ratio of the molar flow rate of solvent and antisolvent must be 91%, as for initial case 0. This strategy yields 5 cases bringing the production flow rate up to nearly 90 times the initial production flow rate. The main operating conditions of the simulation cases are presented in Table 1.

Case

Governing equations of the fluid mixing.

All the equations established in previous studies [START_REF] Zhang | Mixing intensification under turbulent conditions in a high pressure microreactor[END_REF][START_REF] Zhang | Investigating laminar mixing in high pressure microfluidic systems[END_REF][START_REF] Jost | Direct numerical simulations of fluids mixing above mixture critical point[END_REF], are briefly presented. The hydrodynamic Navier-Stokes equations are solved for totally miscible and isothermal fluids. It is generally considered that a supercritical fluid is compressible (especially when close to the critical point). In our case, the fluid is quite above the critical point so the isothermal compressibility is relatively low (between 10 -8 and 10 -9 Pa - 1 ). Furthermore, as the flow is isothermal the thermo-compressible effects are negligible. The comparison of simulations between an incompressible formulation with variable density and a full compressible formulation [START_REF] Sharma | Vibration-induced thermal instabilities in supercritical fluids in the absence of gravity[END_REF][START_REF] Amiroudine | A Lagrangian-Eulerian compressible model for the trans-critical path of near-critical fluids[END_REF] has shown that the results were very close. As the CPU time was optimized within the incompressible formulation, as described in the numerical methodology section, we considered this formulation to solve the hydrodynamics phenomena.

As such, the fluid mixture flow composed of fully miscible ethanol solvent and CO2 antisolvent, is governed by the Navier-Stokes equations:

𝛁. 𝐮 = 𝟎 (1) 𝜌 ( 𝜕𝐮 ∂𝑡 + 𝐮. ∇𝐮) = -∇𝑝 + ∇. (𝜇(∇𝐮 + ∇ T 𝐮)) (2)
where u is the fluid velocity. p is the pressure,  is the fluid density, and  is the fluid viscosity.

As the Froude number is much higher than one resulting from a high Reynolds number flow for all reactor sizes in the range from micro-to milli-, gravity is neglected in the model.

The Reynolds number of the simulation ranges from 5245 to 71494. One of the strengths of the current work is its ability to perform DNS up to the Kolmogorov scale, thereby not requiring any turbulence model and yielding more accurate results.

At the inflow boundary, a laminar Poiseuille flow profile is imposed for the injector whereas the streamwise velocity component ux has a constant profile in the coflow and is superimposed by a uniformly distributed fluctuations ranging between ±0.05 ux [START_REF] Jost | Direct numerical simulations of fluids mixing above mixture critical point[END_REF]. A zero-gradient boundary condition is applied at the outlet boundary (length of the pipe is long enough to justify this approximation) and no-slip boundary conditions at the remaining boundary conditions.

The species mass fraction of the ethanol solvent is governed by the following advection/diffusion equation

𝜌 𝜕𝑥 𝐸𝑡𝑂𝐻 𝜕𝑡 + ∇ ⋅ (𝜌𝑥 𝐸𝑡𝑂𝐻 𝒖 + 𝜌𝐷 𝐸𝑡𝑂𝐻 ∇𝑥 𝐸𝑡𝑂𝐻 ) = 0 ( 3
)
where DEtOH is the diffusion coefficient.

The species mass fraction of the antisolvent CO2 is given by

𝑥 𝐶𝑂 2 = 1 -𝑥 𝐸𝑡𝑂𝐻 (4)
The thermophysical properties of the CO2 and ethanol fluid mixture are computed thanks to the cubic Peng-Robinson equation of state coupled with the Van der Waals mixing rules [START_REF] Maeta | Measurements of vapor-liquid equilibrium in both binary carbon dioxide-ethanol and ternary carbon dioxide-ethanol-water systems with a newly developed flow-type apparatus[END_REF]:

𝑝 = 𝑅𝑇 𝑉 𝑚 -𝑏 𝑚 - 𝑎 𝑚 𝑉 𝑚 (𝑉 𝑚 +𝑏 𝑚 )+𝑏 𝑚 (𝑉 𝑚 +𝑏 𝑚 ) (5)
where Vm is the molar volume (Vm = ZmRT/p with Zm the compressibility factor) and am and bm are the attraction parameter and the co-volume parameter of the mixture, respectively:

𝑎 𝑚 = ∑ ∑ 𝑥 𝑖 𝑥 𝑗 𝑎 𝑖𝑗 𝑛 𝑗=1 𝑛 𝑖=1
and 𝑏 𝑚 = ∑ ∑ 𝑥 𝑖 𝑥 𝑗 𝑏 𝑖𝑗 𝑛 𝑗=1 𝑛 𝑖=1 [START_REF] Berton | Scaling continuous API synthesis from milligram to kilogram: extending the enabling benefits of micro to the plant[END_REF] with

𝑎 𝑖𝑗 = (1 -𝑘 𝑖𝑗 )(𝑎 𝑖𝑖 𝑎 𝑗𝑗 ) 0.5 𝑏 𝑖𝑗 = (1 -𝑙 𝑖𝑗 )
𝑏 𝑖𝑖 +𝑏 𝑗𝑗 2 [START_REF] Rossetti | Chemical reaction engineering, process design and scale-up issues at the frontier of synthesis: Flow chemistry[END_REF] and

𝑎 = 0.45724𝛼𝑅 2 𝑇 𝑐 2 𝑝 𝑐 𝑏 = 0.0778𝑅𝑇 𝑐 𝑝 𝑐 (8) 𝛼 = (1 + (0.37464 + 1.5422𝜔 -0.26992𝜔 2 ) (1 -√ 𝑇 𝑇 𝑐 )) 2 (9)
where the binary interaction parameters kij and lij are given by

𝑘 𝑖𝑗 = -0.4652𝑒 -3 𝑇 + 0.238 (10)
𝑙 𝑖𝑗 = -0.8116𝑒 -3 𝑇 + 0.2491 (11)
The viscosities of the pure fluids CO2 and ethanol are obtained from the NIST database and the viscosity of the fluid mixture is given by

𝑙𝑛𝜇 𝑚 = 𝑥 𝐸𝑡𝑂𝐻 𝑙𝑛𝜇 𝐸𝑡𝑂𝐻 + 𝑥 𝐶𝑂 2 𝑙𝑛𝜇 𝐶𝑂 2 (12)
Finally. the diffusion coefficient of ethanol in CO2 is given by Hayduck-Minhas correlation [START_REF] Zhang | Investigating laminar mixing in high pressure microfluidic systems[END_REF]:

𝐷 𝑚 = 1.33𝑒 -7 𝑇 1.47 𝑉 𝐸𝑡𝑂𝐻 -0.71 𝜇 𝐶𝑂 2 10.2 𝑉 𝐸𝑡𝑂𝐻 -0.791 [START_REF] Dong | Scale-up of Micro-and Milli-Reactors: An Overview of Strategies, Design Principles and Applications[END_REF] where VEtOH is the molar volume of pure ethanol computed by the Peng-Robinson equation of state.

Numerical modeling and computational strategy. To tackle the numerical challenge of DNS of a supercritical microfluidic antisolvent process in large tubes, the Cartesian grid-based finite volume CFD open-source code Notus [25] (developed in-house) is used. The code is based on the finite volume method on a Cartesian grid with explicit or implicit discretization of the governing equations. The velocity-pressure coupling is solved thanks to the Goda pressure correction method (time splitting) [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows[END_REF]. Notus has been widely verified and validated for single and multiphase flows (see for instance [START_REF] Jost | Direct numerical simulations of fluids mixing above mixture critical point[END_REF]25,[START_REF] Jost | Direct forcing immersed boundary methods: Improvements to the ghost-cell method[END_REF][START_REF] Coquerelle | A fourth-order accurate curvature computation in a level set framework for two-phase flows subjected to surface tension forces[END_REF]). Concerning the present case, similar simulations have been validated by micro-PIV [START_REF] Zhang | Investigating laminar mixing in high pressure microfluidic systems[END_REF] and experimental validation for simulations including nucleation and growth phenomena have also been performed [START_REF] Jaouhari | Process intensification for the synthesis of ultra-small organic nanoparticles with supercritical CO2 in a microfluidic system[END_REF][START_REF] Neurohr | Challenge of the supercritical antisolvent technique SAS to prepare cocrystal-pure powders of naproxen-nicotinamide[END_REF].

For the current study, the code has been adapted to modern petascale supercomputers such as the Joliot-Curie supercomputer. It is based on AMD Epyc2 dual processor nodes with each processor having 64 cores, resulting in a total of 293 376 cores. This opens new possibilities in numerical simulations. The upscaling of Notus from thousands of cores to more than a hundred thousand cores has been made possible thanks to fully explicit numerical schemes that considerably reduce the amount of time spent in the linear system solvers. This approach takes advantage of favorable physical fluid characteristics that do not significantly reduce the Courant-Friedrichs-Lewy condition and thus maintain a low total number of time iterations of the simulations. A second order total variation diminishing (TVD) centrally discretized scheme with flux-limiters is used, which ensures a good compromise between accuracy and CPU time.

This approach has been recently proven to be efficient for DNS of supercritical microfluidic antisolvent process in micro-tube [START_REF] Jost | Direct numerical simulations of fluids mixing above mixture critical point[END_REF]. The only linear system that needs to be solved is associated with the pressure variable for which the very efficient PFMG geometric multigrid solver of the Hypre library [START_REF] Baker | Scaling hypre's multigrid solvers to 100,000 cores[END_REF] is used. Moreover, for such a large number of cores some communication bottlenecks between MPI processes appear, and have been overcome thanks to a change in the parallel programming paradigm, from a full MPI to a hybrid OpenMP/MPI approach. Before running large scale simulations, an analysis of the performance of the Notus code on such an architecture and programming environment is performed. Table 2 shows that the average number of cells per processor is reduced to 40 3 while keeping a very good scalability, thereby allowing to perform simulations on 131072 cores for the largest case we propose. This point is crucial to attain reasonable total run times of around only 48 hours for a mesh size of more than 10 billion cells.

Care must also be taken with the file system Input/Ouput (I/O) and data visualization strategies.

Two kinds of I/O are pursued. The first one deals with the code's checkpoint/restart procedure which needs to write on the parallel file system all the necessary information to restart a simulation without any loss of information. This leads to huge files of up to 3.6 TB which, thanks to the ADIOS library [START_REF] Liu | Hello ADIOS: the challenges and lessons of developing leadership class I/O frameworks[END_REF], are efficiently written. On a Lustre parallel file system, aggregate I/O rates of around 50 GB/s have been measured. Moreover, for visualization purposes, I/O performance constraint as well as a data file format that support domain decomposition while reading a file with a parallel visualization software (VisIt for this study) must be satisfied. For this purpose. the HDF5 Pixie file format is chosen and yields I/O performances of around 10 GB/s. The total amount of data needed to realize the proposed video of the present article is around 100 TB.

Results

The main operating conditions of the simulation cases are presented in Table 1. To capture all the mixing scales and to yield accurate results the mesh size of the simulation is chosen to be on the order of the micrometer (see Table 3). As further demonstrated, the Kolmogorov and Batchelor scales are on the same order of magnitude which is particularly advantageous for the current numerical simulations. Indeed, by solving the flow characteristics at the Kolmogorov length scale we resolve all of the mixing scales. The simulations were performed on 512 to 131072 processor cores on the Joliot Curie supercomputer, ranked 38th on the Top 500 list of the largest supercomputers [START_REF]TOP 500[END_REF] with mesh sizes ranging from 55 million to 11 billion (see Table 3). This yields extremely accurate results especially in terms of the mixing performance of the reactors. As detailed in a comprehensive review [START_REF] Sreenivasan | Turbulent mixing: A perspective[END_REF], a slope of -1 is only valid for Sc>>1.

Nb of Cores

For Sc = O(1), such as in our case, it is shown that the slope is not determined (probably around -3, a value close to the one obtained in the present work).

Qualitative flow and mixing analysis.

As the current study is focused on the evaluation of the mixing performance of the tube-in-tube processes, only the analysis of concentration fields of the species is presented. Fig. 3 shows the evolutions of the ethanol mass fraction in the median planes of the reactors for each configuration. As a first observation, the mixing of the species is correctly achieved in the domains and seem to be very fast. It should be recalled that the constraints imposed when determining the conditions for the scaling, such as the ratio of flow rates and the overall dissipation rate in the system, lead to different inlet velocity ratios which, subsequently, result in different jet shapes. However, overall, there are two different configurations. The first one, visible in cases 0, 3, and 4 is characterized by a relatively long jet. It should be noted that the potential core zone may vary significantly (especially in case 3)

for this first configuration. This jet morphology appears in the case of moderate velocity ratios (3 and 4). The second configuration concerns high velocity ratios such as those used in cases 1 and 2. In this configuration, the internal solvent jet breaks up as soon as it leaves the injector. This is due to a very high velocity in the coflow (up to 9 m/s in case 2) resulting in extremely intense and fast mixing.

As discussed in our previous paper [START_REF] Jost | Direct numerical simulations of fluids mixing above mixture critical point[END_REF], the turbulent transition of these coaxial jets is due to the significant development of primary and secondary Kelvin Helmoltz instabilities resulting from strong shear stress between layers of different densities. DNS emphasizes the multitude of scales and highlights the primary and secondary KH instabilities. These are clearly visible.

for example, in case 4 in Fig. 3 and4. Even though a simple visual analysis of the mass fraction fields leads to the conclusion that the mixtures are very good, the performances of the mixtures are further analyzed in order to quantify the influence of the changes of scale and its relevance in terms of processes. Quantitative mixing analysis.

To estimate the performance of the mixing in the different cases, we have chosen to analyze mixing both from a global point of view and from a local point of view thanks to the intensity of segregation and to the local estimation of micromixing time, respectively.

Global evolution of the mixing performance are checked through the segregation intensity of mixing [START_REF] Zhang | Mixing intensification under turbulent conditions in a high pressure microreactor[END_REF][START_REF] Zhang | Investigating laminar mixing in high pressure microfluidic systems[END_REF][START_REF] Jost | Direct numerical simulations of fluids mixing above mixture critical point[END_REF] based on the length of the reactor. The intensity of segregation is defined by [START_REF] Danckwerts | The effect of incomplete mixing on homogeneous reactions[END_REF]:

𝑰 𝒎 = ∑ (𝒙 𝑬𝒕𝑶𝑯 -𝒙 𝑬𝒕𝑶𝑯 ̅̅̅̅̅̅̅̅) 𝟐 𝑵 𝟏 𝑵(𝒙 𝑬𝒕𝑶𝑯 ̅̅̅̅̅̅̅̅(𝒙 𝑬𝒕𝑶𝑯 -𝒙 𝑬𝒕𝑶𝑯 ̅̅̅̅̅̅̅̅)) (14)
A segregation value of 1 and 0 correspond to total segregation between the species and perfect mixing, respectively. This index is calculated for each abscissa x in the yz plane [START_REF] Zhang | Investigating laminar mixing in high pressure microfluidic systems[END_REF]. Figure 5 we propose the evolution of the segregation index as a function of the length for each reactor.

Fig. 5. Evolution of the segregation intensities for all cases in function of the length of the reactor

Although mixing is achieved in a very short distance regardless of the case considered, two different situations can be clearly distinguished. The first concerns cases 1 and 2 for which the mixing seems to be extremely rapid and complete and occurs in the same range than the initial case (case 0). The second one concerns the bigger reactors which follows the same trend, a very good mixing, but with a slight decrease of the mixing performance compared to the original one. Indeed, the intensity of segregation decreases at a slower rate and does not reach zero at the end of the reactor. The intensity is equal to 0.018 at the end of the reactor indicating a good, but not fully complete, mixing. This remark is consistent with the evolution of the mass fraction field (figure 3) for the two larger reactors for which some concentration gradients are still present at the outlet. It should be noted, however, that the increases in the dimensions of the reactor are accompanied by increases in the external diameters to maintain acceptable pressure drops in the system. Therefore, the confinement of the jets becomes less important and tends to decrease the mixing performance of the reactors.

To complete the mixing performance analysis we currently propose an estimation of the mixing time. For that, we analyzed the statistic distribution of the local mixing time. As demonstrated in literature [START_REF] Bałdyga | Turbulent micromixing in chemical reactorsa review[END_REF], mixing can occur at different scales, such as macro-, meso-and micro-scale.

Many studies have emphasized the direct influence of the micromixing in obtaining favorable conditions for the process. Baldyga and Bourne [START_REF] Baldyga | A Fluid Mechanical Approach to Turbulent Mixing and Chemical Reaction Part Ii Micromixing in the Light of Turbulence Theory[END_REF] first proposed a definition of the micromixing time. They deduced from Kolomogorov's theory that this crucial characteristic time is directly correlated to the kinematic viscosity and the turbulent kinetic energy dissipation rate:

𝑡 𝑚 = 17.24√ 𝜐 𝜀 (15)
In this equation, the value of the constant (17.24) is subject of discussion in literature. Recently [START_REF] Lemenand | Micro-mixing measurement by chemical probe in homogeneous and isotropic turbulence[END_REF], authors have suggested a value of 12.7, very close to the constant of the hydrodynamic lifetime of vortex [START_REF] Baldyga | A Fluid Mechanical Approach to Turbulent Mixing and Chemical Reaction Part Ii Micromixing in the Light of Turbulence Theory[END_REF]. As this difference does not affect the mixing analysis of the present study, we considered in the present work the original model from Baldyga.

Direct numerical simulations allow the accurate determination of the local evolution of the dissipation rate of the turbulent kinetic energy ε. The latter is defined by: [START_REF] Bałdyga | Particle formation in supercritical fluids-Scale-up problem[END_REF] with ux', vy', wz' the velocity fluctuations in each direction. Figure 6 shows the fields of the dissipation rate of the turbulent kinetic energy calculated by equation 16 in the median plane of the reactors for each case and allows one to locate the areas of intense mixing and therefore provide information on the mixing. The range of the TKE dissipation rate is similar for all simulations indicating that mixing occurs probably in the same temporal range. As already observed in figure 5, mixing occurs in a short length after the tip of the injector. It can be observed that for case 1 and 2, strong mixing zones are observed in the periphery of the jet due to a high velocity ratio between the injection and coflow zones. To estimate the characteristic mixing time in the reactor, TKE dissipation rate fields (Fig. 6) have been used to estimate the micromixing time according to the equation [START_REF] Rehage | A new scale-up method for competitive chemical model reactions based on complete similarity[END_REF]. As such, to compare the mixing performances of the reactors, the micromixing time of each cell in the domains is directly deduced from the dissipation rate of the turbulent kinetic energy (Fig. 6). It is also interesting to note that the distributions of the mixing time for cases 3 and 4 are significantly wider than the previous ones. This is probably due to the lower velocity ratios and reduced confinement. These observations support the conclusion given by the evolution of the intensity of segregation in the reactor. So, even if there are some differences in the distributions of the mixing time, the order of magnitude is well preserved and shows that the currently proposed strategy allows significant increases in the production flow rate while keeping very favorable mixing conditions. It is worth noting that the increases in the dimensions of the reactor are accompanied, in our simulations, by a slight but acceptable decrease of the mixing performance. This suggests that the scaling up strategy can only be pursued up to a certain dimension, after which it will be impossible to combine increased production and conservation of mixing performance. Let us note that in our study, the last case corresponds to the maximum flow rate achievable by the usual CO2 pump. Finally, increases of the Reynolds numbers due to the increase in reactor size is needed to maintain the mixing performance. Furthermore, the velocity ratio is also very important and a high value contributes to increases in the mixing performance as demonstrated for the cases 1 and 2.

Conclusion

In this paper we are interested in the open-ended problem of the scaling up of a continuous flow process. Towards this end, we propose direct numerical simulations from lab/academic scale to production/industrial scale. The strategy consists in increasing the dimensions of the reactor whilst maintaining the same accuracy of the physical phenomena description. This results in DNS from micro-scale reactors to milli-scale reactors with production rates of up to two orders of magnitude larger. From a chemical engineering point of view, these simulations show that it is possible to increase the volume of the continuous reactor (the diameter has been enlarged by a factor of 6) while maintaining an excellent mixing quality, which is very favorable for the material synthesis process. With a larger reactor, the production rate has increased by a factor 100 while conserving the excellent performance of mixing, thereby indicating that this first step of sizing up yields a decrease, by two orders of magnitude, in the number of units in the numbered-up platforms. This leads to an assurance of a better distribution of the fluid in the system and also a reduction in the manufacturing costs and complexities.

To the best of our knowledge, the simulation performed on a mesh size of 11 billion cells and 131072 processors is the largest realized in the field of continuous flow processes. It shows that petascale high fidelity numerical simulations is an emerging and efficient tool for the analysis of different scaling up scenarios or reactor designs in the field of chemical engineering that open new insights at industrial scale. Indeed, all relevant space and time scales are resolved which gives information on the flow and mixing characteristics unattainable through in situ experiment. Moreover, the amount of data produced and their considerable precision are of primary interest to develop data reduction and learning AI algorithms to go through real time simulation of Smart Manufacturing processes.

We believe that this work helps to follow the increasing computing power of the supercomputers since these available quasi unique supercomputers will be increasingly available in the coming decade. Since the development of supercomputers, each new subsequent generation brings its own set of demonstrations of their capacity to handle more complex problems, higher Reynolds number DNS, etc. This work is part of this logic and is a starting point for further new challenges for continuous flow processes for materials synthesis that include, non isothermal processes, reaction kinetics, and geometry changes (zig-zag channels, T injectors, etc.). Regarding the larger reactors (stirred tank reactor, tower reactor, and fluidized bed reactor), the capability to solve all mixing scales without any turbulence model depends on the operating conditions and the size of the reactor. If DNS is not feasible, Large Eddy Simulation can be used to capture turbulence phenomena in the reactor and yield accurate local and temporal analysis of the flow field and, thereby, resulting in a successful scale up process. It should be noted that LES in large reactors also requires very large simulations.

with

  𝑓 = 0.3164𝑅𝑒 -0.25 being the Blasius friction factor for turbulent flow in pipes, u the velocity, and D the diameter in the considered zone (coflow zone D = Dext-Do, inner tube D = Di, outer tube D = Dext).

Fig 1 .

 1 Fig 1. Diagram of the tube-in-tube reactor. The solution, solvent, is injected in the inner tube, diameter Di, whereas the supercritical CO2, acting as antisolvent, is injected in the coflow between the diameters Do and Dext. The typical pressure, temperature and solvent/antisolvent molar ratio are 100 bar, 40°C and 91%, respectively, are chosen.

.

  Respectively, are performed. In the former relations, k is the wave number, and ux' and xs' are the velocity and mass fraction fluctuations, respectively. The fluctuation of a variable is defined as the difference between the variable's instantaneous value and its time averaged. i.e. 𝑢 𝑥 = 𝑢 ̅ + 𝑢 𝑥 ′ , where 𝑢 ̅ is the time averaged value. The turbulence spectra are then obtained by performing a Fast Fourier Transform with a Hann window of the desired quantities along the centerline of the domain and are shown in Fig. 2. As expected, regardless of the case, the spectra decay with slopes of approximately -5/3, as predicted by the Kolmogorov theory. It is widely accepted in the turbulence community that a slope of -5/3 is a good indicator that a simulation has captured the smallest scales of the turbulence. Moreover, as already mentioned above, as the Batchelor scales are on the same order of magnitude as the Kolmogorov scales, all the scales of interest of the mixtures have been captured which is the sine qua non condition of the forthcoming quantitative analysis of the mixture A precise analysis and interpretation of the second part of the mass fraction turbulence spectra is still an open question.

Fig. 2 .

 2 Fig. 2. Turbulence spectra analysis. 1D turbulence spectra for velocity RMS (right) and mass fraction fluctuation (left). The turbulence spectra are then obtained by performing a Fast Fourier Transform with a Hann window of the desired quantities along the centerline of the domain and are shown in Fig. 2. As expected. regardless of the case. the spectra decay with slopes of approximately -5/3, as predicted by the Kolmogorov theory.

Fig. 3 .

 3 Fig. 3. Instantaneous mass fraction fields in the median planes of the reactors for the cases 0 to 4 (top to bottom). As a first observation. the mixing of the species is correctly achieved in the domains and seem to be very fast. Secondly. different jet shapes are observed. Cases 0, 3, and 4 are characterized by a relatively long jet. With high velocity ratios such as those used in cases 1 and 2 the internal solvent jetbreaks up as soon as it leaves the injector. This is due to a very high velocity in the coflow (up to 9 m/s in case 2) resulting in extremely intense and fast mixing. Full length domains are presented, true scales between the cases are also respected.

Fig 4 .

 4 Fig 4. Instantaneous concentration volume renderings for the cases 0 to 4 (top to bottom) -see additional videos 2-6. The turbulent transitions of the coaxial jets is due to the significant development of primary and secondary Kelvin Helmoltz instabilities resulting from strong shear stress between layers of different densities. DNS emphasizes the multitude of scales and highlights the

Figure 6 :

 6 Figure 6: Instantaneous turbulent kinetic energy dissipation rate for the cases 0-4.

  Figure 7 only shows the normalized volume distributions of the micromixing time in the median planes of each reactor. Due to the reactors' cylindrical geometry, the distribution in the median plan is statistically representative of the 3D configuration. The first important result is the mean value of the micromixing time obtained for the initial case (case 0). As expected, we found a mean value between 100 and 200 s, indicating a very fast mixing leading to very favorable conditions for material synthesis. The previous observation of the mass fraction fields of cases 1 and 2 suggests an extremely fast mixing and is confirmed by the narrow distributions of the micromixing time for which the average values are around 50 s. The high velocity ratios and the important confinement in the reactors seem to lead to very favorable mixing efficiency conditions. Cases 3 and 4 result in mean micromixing times that are slightly higher than the initial case at around 200 s. These are still excellent mean mixing times for the elaboration of materials.

Fig. 7 .

 7 Fig. 7. Normalized volume distributions of the micromixing time in the median planes of the reactors for all cases. For case 0 we found a mean value between 100 and 200 s. indicating a very fast mixing. Cases 1 and 2 suggests an extremely fast mixing and is confirmed by the

Table 1

 1

		R	QEtoh	Qco2	VEtoh	Vco2	𝒖 ̅	CO2%		Re	Max. p	Di	Do	Dext
			(l/min)	(l/min)	(m/s)	(m/s)	(m/s)	mol	(W/kg)		(bar/m)	(m)	(m)	(m)
	0	1	400	3000	0.87	3.35	1.47	91	235	5245	16	100	170	250
	1	9	3539	25542	1.23	6.08	2.18	91	235	18496	21	250	450	600
	2	31	12400	89495	1.07	9.32	2.75	91	235	38884	31	500	800	1000
	3	55	21867	157820	1.89	6.25	2.75	91	235	54420	5	500	800	1260
	4	87	34655	250117	3.00	5.61	3.33	91	235	71494	3	500	800	1520

. Operating conditions of the simulation cases. R: production ratio. QEtoh: volume flow rate of solvent. Qco2: volume flow rate of antisolvent. VEtoh: velocity in inner tube. Vco2: velocity in the coflow zone. Max. P: the maximum calculated linear pressure drop in the reactor (injector or coflow zone).

Table 2 . Weak scalability study with 40 3 cells per core for an OpenMP/MPI strategy.

 2 The full MPI approach results in communication bottlenecks for more than 32768 cores; the 4 OpenMP threads / process MPI approach reduces the CPU time of the Notus code and of the Hypre library by a factor 3.

process MPI 2 threads / process MPI 4 threads / process MPI 8 threads / process MPI

	Reactor 1 thread / Case Rector Length (mm) diameter(mm)		Mesh size	Approx. Number of cells (x10 9 )	Mesh step size (x10 -6 )	Number of processors
	0		Mesh size 3	Notus (s) 0.2	Hypre (s) 2000 x 167 x 167 Notus (s)	Hypre (s)	Notus (s) 0.055	Hypre (s)	Notus (s) 1.5	Hypre (s) 512
	1	128 512	8.2E+06 3.3E+07 4	6.1 6.1 0.6	3.6 4.0 4000 x 600 x 600 5.6 6.1	3.3 5.8	1.4	6.1 6.5	3.4 4.5	1	7.8 7.9	3.5 4.3 32768
	2048 2	1.3E+08 5	6.6 1	4.7 5000 x 1000 x 1000 6.4	4.3	5	6.6	4.2	1	8.5	5.9 32768
	8192	5.2E+08	7.2	5.6	6.9	5.1		7.3	5.1		9.1	5.5
	32768 3	2.1E+09 5	8.7 1.26	9.1 5000 x 1260 x 1260 8.4	7.1	8	8.5	6.6	1	9.9	7.2 65536
	131072 4	8.4E+09 5	41.6 1.52	34.1 5000 x 1520 x 1520 14.5	13.0	11.5	14.1	10.5	1	19.1	11.7 131072

Table 3 . Reactor dimensions, mesh size and number of processors used for cases 0 to 4. Verification of resolved finest scales

 3

. First, to verify that the energy dissipation scales, i.e. the Kolmogorov length scales, are resolved in all the simulations currently pursued, analysis of the 1D turbulence spectra Eu and Ex of the root mean squared (RMS) of velocity and of the mass fraction fluctuations. Respectively, expressed as 〈𝑢 𝑥 ′ 𝑢 𝑥 ′ 〉 = ∫ 𝐸 𝑢 (𝑘)𝑑𝑘

Acknowledgments

We acknowledge the French National Research Agency for its support (ANR-17-CE07-0029 -SUPERFON). the MCIA (Mésocentre de Calcul Intensif Aquitaine) and GENCI (DARI project number A0062A10815 and "Les Grands Challenges" grant on Joliot-Curie supercomputer) for the HPC resources.


S S Shelistov

E A Demekhin

V S Shelistov

G S Ganchenko

Electric-permittivity-based instability of two dielectric miscible liquids under DC field

Electric-permittivity-based instability of two dielectric miscible liquids under DC field

Sakirhoussein Amiroudine, E. A. A Demekhin, V.

Introduction

Micro-and nanoscale phenomena of moving conductive and non-conductive liquids with an interface between them are of particular practical interest. Multi-phase flows in micro-scale systems have recently found numerous applications in a wide range of fields like biochemical processing, such as lab-on-chip reactors [1][2][START_REF] Gao | Colloids Surfaces Physicochem[END_REF][4], mixers [5,6], DNA extraction [7], drug delivery [START_REF] Squires | Induced-charge electroosmosis[END_REF], oil extraction from porous rock formations [9,10]. Miniaturized bioanalytical systems attempt to incorporate many of the necessary components and functionality on the surface of a typical laboratory substrate. Micro Total Analysis Systems (µTAS) components include reaction chambers, pumps, flow sensors, micromixers, diluters, and preconcentrators. Analyses requiring rapid mixing include immunoassays, DNA hybridization, and general cell-molecule interaction. Application of these techniques requires mixing of reagents that have relatively low diffusion coefficients.

Rapid homogeneous mixing becomes increasingly important [11][12][13][14][15][16][17][18][19][20][21][22][23][24]. Such a mixing can be caused by either mechanical vibration through hydrodynamical instability [25][26][27] or by a special type of hydro-electro instability [28][29][START_REF] Melcher | Continuum Electromechanics[END_REF][START_REF] Saville | [END_REF][32][33]. Electroconvective mixing may result from many physical factors, such as the first kind electroosmotic flow [START_REF] Di Fraia | Modelling electro-osmotic flow in porous media: a review[END_REF], the second kind electroosmotic flow [START_REF] Rubinstein | Electro-osmotically induced convection at a permselective membrane[END_REF][START_REF] Liu | Shear electroconvective instability in electrodialysis channel under extreme depletion and its scaling laws[END_REF], buoyant convection [START_REF] Shi | Length-dependent instability of shear electroconvective flow: from electroconvective instability to Rayleigh-Bénard instability[END_REF], conductivity gra-dient [START_REF] Hoburg | Electrohydrodynamic mixing and instability induced by co-linear fields and conductivity gradients[END_REF], and electric permittivities [START_REF] Jalaal | Electrohydrodynamic (EHD) mixing of two miscible dielectric liquids[END_REF]. The slow dif-fusion mechanism of mixing is thus complemented by much faster advection mixing.

The electrohydrodynamic instability of liquid flows with conductivity gradients stems from pioneer works by Hoburg and Melcher [28] which are based on classi-cal models by Melcher and Taylor [29] (see also [START_REF] Melcher | Continuum Electromechanics[END_REF][START_REF] Saville | [END_REF]). Hoburg and Melcher [28] described the key mechanism of this instability as caused by charge accumulation at the perturbed interface, and made qualitative compar-isons of their theory to experiments, but they neglected the diffusion of electric conductivity. As a result, in their calculation, the threshold of instability was absent. This was improved later by Baygents and Baldessari [32] with a more realistic model. However, Baygents and Baldessari [32] wrongly assumed the principle of exchange of stability in their linear stability analysis. This was corrected subsequently by Chang et al. [40] and Sharan et al. [41]. Research on this electrohydrodynamic instability was continued by Santiago's team [16][17][18][19][20]33]. A flow in a long rectangular-cross-section microchannel with a conductity gradient orthogonal to the main flow direction and an external electric field was considered both experimentally and theoretically in [16].

a e-mail: sakir.amiroudine@u-bordeaux.fr (corresponding author)

It was found experimentally that such a system exhibits a critical electric field above which the flow is unstable. In the theoretical part, the previous model was generalized in [32]. The volume charge was assumed small enough to be neglected in the ion transport equations, but not in the equations of fluid motion, where the Coulomb body force, based on the residual charge, was taken into consideration.

The electrohydrodynamic instability of liquid flows with conductivity gradients stems from pioneer works by Hoburg and Melcher [28] which are based on classical models by Melcher and Taylor [29] (see also [START_REF] Melcher | Continuum Electromechanics[END_REF][START_REF] Saville | [END_REF]). Hoburg and Melcher [28] described the key mechanism of this instability as caused by charge accumulation at the perturbed interface, and made qualitative comparisons of their theory to experiments, but they neglected the diffusion of electric conductivity. As a result, in their calculation, the threshold of instability was absent. This was improved later by Baygents and Baldessari [32] with a more realistic model. However, Baygents and Baldessari [32] wrongly assumed the principle of exchange of stability in their linear stability analysis. This was corrected subsequently by Chang et al. [40] and Sharan et al. [41]. Research on this electrohydrodynamic instability was continued by Santiago's team [16][17][18][19][20]33]. A flow in a long rectangular-cross-section microchannel with a conductity gradient orthogonal to the main flow direction and an external electric field was considered both experimentally and theoretically in [16]. It was found experimentally that such a system exhibits a critical electric field above which the flow is unstable. In the theoretical part, the previous model was generalized in [32]. The volume charge was assumed small enough to be neglected in the ion transport equations, but not in the equations of fluid motion, where the Coulomb body force, based on the residual charge, was taken into consideration. A general linear stability analysis of miscible sys-tems is not trivial due to the fact that the base state of the problem is time-dependent and the correspond-ing linear system ceases to be of classical type. Specif-ically, the case considered in this paper differs, in two ways, from the usually considered above-mentioned cases of electrohydrodynamic instability: (a) the liquids are dielectrics rather than conductors and, hence, the Maxwell-Wagner stresses are created by the nonuni-formity of the dielectric permittivity; (b) the liquids are miscible and the boundary separating them spreads with time by diffusion. We consider two viscous miscible dielectric liquids in a2 ˜h gap between conducting impermeable walls, see Fig. 1. Notations with tilde are used for the dimen-sional variables, as opposed to their dimensionless counterparts without tilde. (x, y ˜) are the Cartesian coordinates, where x ˜ is directed along the channel and y ˜ is normal to it.

A potential difference ΔV ˜ is applied between the walls; any external forces (such as gravity), except the electrical one, are neglected. It is assumed that at the initial time, t ˜ = 0, the interface is a straight line y ˜ = 0 perturbed with small-amplitude natural "room perturbations" and the first liquid is located in the h < ỹ < 0 region, while the second liquid is in the region 0 < ỹ < h. These liquids are denoted by the indices of 1 and 2, respectively. At time t > 0, the interface is subject to mutual diffusion of these liquids and an electrohydrodynamic instability may arise. The following two-dimensional equations describe the mixture behavior,

ρ ∂w ∂ t + Ũ ∂w ∂ x + Ṽ ∂w ∂ ỹ = D ∂ ∂ x ρ ∂w ∂ x + ∂ ∂ ỹ ρ ∂w ∂ ỹ , ( 1
)
∂ Π ∂ x - ∂ Txx ∂ x - ∂ Txy ∂ ỹ = 1 2 ∂ ∂ x ε Ẽ2 x -Ẽ2 y + ∂ ∂ ỹ ε Ẽx Ẽy , ( 2
)
∂ Π ∂ ỹ - ∂ Txy ∂ x - ∂ Tyy ∂ ỹ = ∂ ∂ x ε Ẽx Ẽy - 1 2 ∂ ∂ ỹ ε Ẽ2 x -Ẽ2 y , ( 3
)
These factors make the problem more difficult, but also more interesting. Besides the practical part connected with the desire of reducing the mixing time, this problem raises some fundamental questions, namely, wether it is possible or not to obtain instability in miscible flows with nonuniformly distributed permittivity of these dielectric liquids. To our knowledge, the stability problem of the interface of two viscous miscible liquids with different electric permittivities has not yet been addressed in the past. The present paper considers a two-phase microflow of dielectric miscible viscous liquids with different elec-tric permittivities bounded by two flat solid walls in an external electric field, normal to the interface. The interface is initially expanding due to diffusion and this process can be described by a self-similar solution. The expanded solution takes into account the nonuniformity of the permittivity of the two liquids which, in principle, can trigger an instability. There is an obvious mathematical difficulty in the solution of the given stability problem: in the self-similar solution, the coefficients of the stability equations depend exponentially on time and an exponential growth rate of the perturbations is thus expected. The classical frozen coefficient approach [START_REF] Lin | The Theory of Hydrodynamic Instability[END_REF][START_REF] Drazin | Hydrodynamic Stability[END_REF] is not applicable in our case. We use sophis-ticated ad hoc method proposed by Shtemler [START_REF] Shtemler | [END_REF]to solve our stability problem. In order to accomplish our investigation and to prove the results of the linear analysis, we apply a numerical integration of the full prob-lem without any simplification. These analytical results are then complemented by direct numerical simulation solutions of the full non-linear problem, and an estimation of the mixing efficiency is proposed as a function of the separation amplitude.

∂ ρ Ũ ∂ x + ∂ ρ Ṽ ∂ ỹ = 0, (4)
∂ ∂ x ε Ẽx + ∂ ∂ ỹ ε Ẽy = 0, (5)

˜˜˜Φ

where w is the mass fraction of the liquid 1 in the mixture with 0 ≤ w ≤ 1; D ˜,ρ,μ and ε ˜ are the diffusion, the density, the dynamic viscosity and the dielectric permittivity of the mixture, respectively; ( ˜U, V ˜)are the components of the velocity vector, ˜Π is the pressure and T ˜ij is the viscous stress tensor. The first equation describes the mixing by diffusion and advection ( [45]), the next three equations are Navier-Stokes equations with Maxwell forces in the right-hand side and the continuity equation. The channel is micron-sized and, thus, the Reynolds number is small so that Stokes' approximation for the creeping flow can be applied. The last equation is the Maxwell equation for the electric field which should be completed by the relation, (E x , E y )=-∇ ˜Φ, where is the electric potential.

In this work, we restrict our analysis to the partic-ular case where the densities of the fluids are equal to each other (this hypothesis will be further verified with experimental data, see Sect. 4.3.4). The liquids are assumed to be Newtonian, but the viscosity μ ˜ of their mixture depends on the mass fraction w, as well as the dielectric permittivity ε ˜. Both dependecies are assumed to be linear functions of w,

μ = μ1 w + μ2 (1 -w) = μ2 + (μ 1 -μ2 )w, ε = ε1 w + ε2 (1 -w) = ε2 + (ε 1 -ε2 )w. ( 6
)
The system (1-5) can be rewrittten into the following form,

∂w ∂ t + Ũ ∂w ∂ x + Ṽ ∂w ∂ ỹ = D ∂ 2 w ∂ x2 + ∂ 2 w ∂ ỹ2 , ( 7
) ∂ ∂ x ε ∂ Φ ∂ x + ∂ ∂ ỹ ε ∂ Φ ∂ ỹ = 0, (8)
∂ Π ∂ x + 2 ∂ ∂ x μ ∂ Ũ ∂ x + ∂ ∂ ỹ μ ∂ Ũ ∂ ỹ + ∂ ∂ ỹ μ ∂ Ṽ ∂ x = - 1 2 ∂ ε ∂ x ⎡ ⎣ ∂ Φ ∂ x 2 + ∂ Φ ∂ ỹ 2 ⎤ ⎦ , ( 9
) ∂ Π ∂ ỹ + 2 ∂ ∂ ỹ μ ∂ Ṽ ∂ ỹ + ∂ ∂ x μ ∂ Ũ ∂ ỹ + ∂ ∂ x μ ∂ Ṽ ∂ x = - 1 2 ∂ ε ∂ ỹ ⎡ ⎣ ∂ Φ ∂ x 2 + ∂ Φ ∂ ỹ 2 ⎤ ⎦ , ( 10
)
∂ Ũ ∂ x + ∂ Ṽ ∂ ỹ = 0. ( 11
)
This system must be complemented by relations (6) and proper boundary conditions (BCs). The solid surfaces y ˜ = ± ˜h are assumed to be impermeable to the mixture, and the no-slip conditions for the velocity are applied at these surfaces. The potential difference ΔV ˜ is applied between the walls. These conditions can be written as follows,

ỹ = ± h : Ũ = Ṽ = 0, ∂w ∂ ỹ = 0; Φ ỹ= h -Φ ỹ=- h = Δ Ṽ . ( 12
)
We neglect the net flow along the x ˜-direction in order to consider only the impact of electric forces on mixing.

The spatial domain is assumed to be infinitely large in the x ˜-direction, and the boundedness of the solution as x ˜ →±∞ is imposed. Therefore, periodic conditions will be assumed with respect to x ˜ in the numerical simula-tions described in Sect. 4.2.

Adding initial conditions for the fraction w completes the system ( 7)- (12). At t ˜=0,

w =1 at -˜h< ỹ<0, w = 0 at 0 < ỹ < h ( 13
)
3 Simplified analytical solution

A one-dimensional self-similar solution

At the initial time t = 0, liquids 1 and 2 are separated by a sharp boundary at ỹ = 0. At subsequent moments, t > 0, mutual diffusion of liquids occurs, which leads to the formation of a thin diffusion layer δ( t). Note that at this period of time, the influence of walls can be neglected and the system has no characteristic size and its solution is thus self-similar.

As long as the diffusion layer is uniform, ∂/∂ x = 0 and its thickness δ is a good choice for the dynamic characteristic length, δ = 2 Dt . The spatial variable ỹ is normalized to δ, and the self-similar variable is introduced as, η = ỹ/ δ. The derivatives with respect to t and ỹ can be written as,

∂ ∂ t = - ỹ 4 t Dt ∂ ∂η , ∂ ∂ ỹ = 1 2 Dt ∂ ∂η . ( 14
)
Equation ( 7) turns into the following ordinary differential equation and BCs,

d 2 w dη 2 + 2η dw dη = 0, η → -∞ : w → 1, η → +∞ : w → 0. ( 15
)
The solution of Eq. ( 15) is a well-known error function,

w = 1 2 [1 -erf(η)] . ( 16
)
The permittivity and the viscosity in the mixing layer can readily be found from relations (6),

ε = ε2 -ε1 2 erf(η) + ε2 + ε1 2 , μ = μ2 -μ1 2 erf(η) + μ2 + μ1 2 . ( 17
)
By introducing dimensionless parameters β and γ as,

β = ε1 ε2 -1, γ = μ1 μ2 -1, (18)
we get the dimensionless solution of the permittivity and the viscosity,

ε = β 2 [1 -erf(η)] + 1, μ = γ 2 [1 -erf(η)] + 1. (19)
The electric potential Φ ˜ is completely determined from the permittivity field ε ˜ (Eq. 8) and the corresponding BC (Eq. 12),

∂ ∂ ỹ ε ∂ Φ ∂ ỹ = 0, Φ| ỹ=-h = 0, Φ| ỹ=+ h = Δ Ṽ . ( 20
)
This problem is readily solved,

Φ = C ỹ - h dỹ ε , ( 21
)
where C is a constant of integration, which is found from the BCs at ỹ = ± h,

Δ Ṽ = C + h - h dỹ ε = C J + h 1 ε1 + 1 ε2 , ( 22
)
in which,

J = 0 -h 1 ε - 1 ε1 dỹ + h 0 1 ε - 1 ε2 dỹ, (23)
We can reduce Eq. ( 23) to dimensionless form using ( 14), ( 18), (19),

Δ Ṽ = C 2 Dt ε2 J + H 2 + β 1 + β , ( 24
)
where H = h/ δ and,

J = 0 -H 1 ε - 1 1 + β dη + H 0 1 ε -1 dη. ( 25
)
Assuming that δ h, the finite limits of integration can be approximated by infinite ones,

J = 0 -∞ 1 ε - 1 1 + β dη + ∞ 0 1 ε -1 dη, (26)
which means that J depends only on β, but not on H. The unknown of integration C ˜ can then be deduced from Eq. ( 24),

C = ε2 Δ Ṽ 2 Dt 1 + β (2 + β)H + (1 + β)J , ( 27
)
with its following dimensionless form,

C = 1 + β (2 + β)H + (1 + β)J , ( 28
)
where

C = 2 √ Dt ε2Δ
Ṽ C. Gravity forces are absent and the problem can thus be taken as symetric and without loss of generality, the bottom and the top liquids can be swapped. Therefore, only positive β can be considered, J(β) is evaluated numerically and tabulated (see Table 1).

The hydrodynamics is decoupled from the onedimensional solution and the velocity components are vanishing (see Eqs. 9-10), but the Coulomb forces produce a nonzero pressure field from Eq. (10). The corresponding solutions for the velocity and pressure field can be written as follows,

Ũ = Ṽ = 0, Π = 1 4 ε ∂ Φ ∂ ỹ 2 . ( 29
)

Stability of the self-similar solution

The interface between the two liquids is always slightly disturbed by natural "room perturbations". These small disturbances can either decay or grow. In the latter case, the 1D solution is unstable and must evolve into some more complex state, than the one obtained in Sect. 3.1.

We consider the linear stability of the 1D self-similar solution and we slightly perturb it in the x ˜-direction. It is more convenient to use different reference values in x ˜ and y ˜-directions, l0 = 1/α 0 : length in x-direction; Here α0 is some characteristic wavenumber and t0 is some characteristic time (we used t0 = l2 0 / D, but the resulting dimensionless equations do not depend on its choice). New independent variables and their derivatives are introduced,

τ = t t0 , ξ = xα 0 , η= ỹ δ = ỹ 2 Dt , ∂ ∂t = ∂ ∂τ - 1 2τ η ∂ ∂η , ∂ ∂ x = α0 ∂ ∂ξ , ( 30
)
∂ ∂ ỹ = 1 δ ∂ ∂η . ( 31
)
In self-similar variables, Eq. ( 7) turns into the following one,

4τ ∂w ∂τ -2η ∂w ∂η + α 2 U ∂w ∂ξ + V ∂w ∂η = α 2 ∂ 2 w ∂ξ 2 + ∂ 2 w ∂η 2 . ( 32
)
The dimensionless dielectric permittivity ε and the dimensionless viscosity μ are linear functions of w and can be found from (6),

ε = 1 + βw(τ, ξ, η), μ= 1 + γw(τ, ξ, η). ( 33
)
The electric potential (Eq. 8) turns into,

α 2 ∂ ∂ξ ε ∂Φ ∂ξ + ∂ ∂η ε ∂Φ ∂η = 0. ( 34
)
The dimensionless parameter α = δ ˜α˜0 can have two different interpretations: (i) a dimensionless diffusion length (derived from δ ˜) or (ii) a dimensionless wavenumber (derived from α˜0). We will use the second choice. This wavenumber depends on time and for sufficiently large time it is assumed to be a slow parameter. The diffusion layer between the two liquids expands exponentially and, therefore, the stability of the self-similar solution cannot be investigated by classical methods as in Refs. [START_REF] Lin | The Theory of Hydrodynamic Instability[END_REF][START_REF] Drazin | Hydrodynamic Stability[END_REF]. The details can be found in Ref. [START_REF] Shtemler | [END_REF], and the corresponding method was successfully applied in [46] and [47]. Details on this will be given later.

Introducing the stream function Ψ from the relations,

U = ∂Ψ ∂η , V = - ∂Ψ ∂ξ , ( 35
)
and substituting them into Eqs. ( 9)-( 10), we get,

∂ 2 ∂η 2 μ ∂ 2 Ψ ∂η 2 + 4α 2 ∂ 2 ∂ξ∂η μ ∂ 2 Ψ ∂ξ∂η + α 4 ∂ 2 ∂ξ 2 μ ∂ 2 Ψ ∂ξ 2 -α 2 ∂ 2 ∂ξ 2 μ ∂ 2 Ψ ∂η 2 -α 2 ∂ 2 ∂η 2 μ ∂ 2 Ψ ∂ξ 2 = Λ ∂ε ∂ξ ∂Φ ∂η ∂ 2 Φ ∂η 2 + α 2 ∂Φ ∂ξ ∂ 2 Φ ∂ξ∂η -Λ ∂ε ∂η ∂Φ ∂η ∂ 2 Φ ∂ξ∂η + α 2 ∂Φ ∂ξ ∂ 2 Φ ∂ξ 2 . ( 36
)
Here, the parameter,

Λ = ε2Δ Ṽ 2 D μ2
is introduced. It characterises the coupling between hydrodynamics and electrostatics and for fixed physical properties, it corresponds to the square of the dimensionless potential drop.

Imposing small perturbations on the 1D self-similar solution, substituting them into Eqs. ( 32)- [START_REF] Liu | Shear electroconvective instability in electrodialysis channel under extreme depletion and its scaling laws[END_REF] and linearizing with respect to the perturbations turns this system into a linear system with the coefficients independent of the coordinate ξ. It allows seeking sinusoidal elementary solution with respect to ξ-coordinate but its dependence on time must be changed from exponen-tial to the power one. The reason comes from the fact that in this linear stability analysis, the coefficient of the time derivative depends itself on time [4τ ∂w in Eq. ∂τ (32)]. As a result, the linear perturbations will not grow exponentially with time, as we used to have in most of the classical linear stability hydrodynamic problems, but in form of power law in time,

w = w 0 (η) + ŵ(η)e iξ τ λ , Φ = Φ 0 (η) + Φ(η)e iξ τ λ , Ψ = -i Ψ(η)e iξ τ λ , ( 37
)
where subscript 0 refers to the 1D self-similar solu-tion and λ is the growth rate. The factor -i before the amplitude of the stream function is introduced in order to make all the system coefficients real. Upon substi-tution of the expressions (37)intoEqs. (32), [START_REF] Di Fraia | Modelling electro-osmotic flow in porous media: a review[END_REF]and [START_REF] Liu | Shear electroconvective instability in electrodialysis channel under extreme depletion and its scaling laws[END_REF], omitting the subscript '0' and after linearization, we get the following system of ODEs,

ŵ -2η ŵ + α 2 w Ψ = ŵ -α 2 ŵ, ( 38
) (ε Φ + εΦ ) -α 2 ε Φ = 0, ( 39
)
μ ΨIV -2α 2 Ψ + α 4 Ψ +γ w Ψ + 2w Ψ -2α 2 w Ψ + α 2 w Ψ = ΛΦ (εΦ -Φ ε ), ( 40
)
where the prime means the derivative with respect to η. The BCs are, ŵ = 0, Φ = 0, Ψ = 0, Ψ = 0 at η = H. (41) Here from the self-similar solution,

ε = 1 + β 2 [1 -erf(η)] + 1, μ = 1 + γ 2 [1 -erf(η)] + 1, Φ = C(H, β) 1 ε . ( 42
)
h responds to a new type of instability in miscible fluids.

According to the classification of Cross and Hohenberg [48], the instability described above is a short-wave instability: the long waves and the very short waves are stable, while the intermediate ones are unstable. This classification determines the behavior of the perturbations in time. The peculiarity of our case is that while the dimensional wavenumber α˜ and the channel width are fixed, their dimensionless counterparts α =αδ ˜ and H = ˜h/δ ˜ are slowly changing in time. The fact that H and α are time-dependent is justified both experimentally and in numerical simulations. In Fig. 3 the results are presented in the coordinates,

1 H = 2 h Dt , α = 2α 0 Dt ,
for different fixed parameters Λ and β, which are independent of time.

The ratio k = α/H -1 does not depend on time, but the position of the point on this k-line depends on time. Each straight line, emanating from the origin (0, 0), with the inclination k, α = k 1 H , characterizes the scenario of the time evolution of the imposed perturbation with dimensional and constant in time wavenumber α0 . The point (1/H, α) on each straight line moves away from the origin as a square root of time, √ t. For sufficiently large k > k * (which corresponds to short waves) the straight line with inclination k will be in the stable region. The corresponding perturbation on the entire line will decay. For k < k * (intermediate and long waves), the straight line will cross the unstable region, but at small times, the perturbation will always decay, reflecting the short-wave instability. For the marginal case, k = k * , the straight line only touches the instability region, but does not cross it. Hence, this results in a rather sophisticated behavior. Depending on the value of k, two different scenarios are possible: (a) stable, when the perturbation always decays; (b) unstable, when the perturbation originally decays, then grows and eventually decays again. Such behavior will be confirmed later in our direct numerical simulation.

The formulated spectral stability problem is described by four outer parameters, β, γ, H and Λ and one inner parameter, the wavenumber, α (C(H, β) is defined via J(β)(seeEq.28) and is presented in Table 1). Note that β characterizes the ratio of the permittivity of the two liquids, γ is the ratio of their viscosity, Λ is the normalized strength of the electric field, H and α are slowly changing parameters with time.

The eigenvalue problem for the growth rate λ was solved numerically by the shooting method, using λ as a parameter, conditions (41)atη = -H as the initial conditions for the Cauchy problem and the corresponding conditions at η =+H as the target.

The dependence of the growth rate λ on the wavenumber α for equal viscosity of liquids γ = 0 for differ-ent dimensionless channel widths H = ˜h/δ ˜ is presented in Fig. 2a. The influence of different viscosities on the dependence λ(α) is shown in Fig. 2b. Increasing the value of γ =μ 1 /μ ˜2 -1leads toasystemwhich is more and more stable. Note that for equal liquid permittivity β = 0 the flow is stable for any value of other parameters. This fact allows us calling this instability as the electric-permittivity-based instability which cor-

Direct numerical simulation

To make the analysis complete, the 1D self-similar analytic solution and its stability investigation are supplemented by the direct numerical simulation (DNS) of the full nonlinear problem.

Dimensionless equations

In order to render the system dimensionless, the following characteristic quantities are chosen in our direct numerical simulations, h : half-width of the channel; h2 / D : time; Ũ0 = D/ h : velocity; Π0 = D μ2 / h2 : pressure; ε2 : permittivity of the fluid 2; Δ Ṽ : drop of potential over the channel width.

Unlike Sect. 3, we select a time-independent characteristic width here, which will simplify the interpretation of the results. The problem is described by the following dimensionless equations,

∂w ∂t + U ∂w ∂x + V ∂w ∂y = ∂ 2 w ∂x 2 + ∂ 2 w ∂y 2 , ( 43
)
∂ ∂x ε ∂Φ ∂x + ∂ ∂y ε ∂Φ ∂y = 0, (44)
∂U ∂x + ∂V ∂y = 0. ( 47
)
The BCs are,

y = ±1 : U = V = 0, ∂w ∂y = 0; Φ| y=+1 -Φ| y=-1 = 1. ( 48
)
and initial conditions are as follows,

t = 0 : w = 1 at -1 < y < 0, w = 0 at 0 < y < 1, ( 49
)
where small perturbations in the x-direction on the interface are assumed to be superimposed on this BC. The problem is characterized by three dimensionless parameters: β, γ and Λ, all defined previously in Sect. 3. The latter can be rewritten as,

Λ ≡ ε2 Δ Ṽ 2 D μ2 = Δ Ṽ 2 D μ2 /ε 2 = Δ Ṽ 2 Δ Ṽ 2 0 ,
which represents the square of the dimensionless drop of the electric potential with the electroviscous potential

Δ Ṽ0 = D μ2 /ε 2 .
In order to distinguish between dimensionless wavenumbers defined in the self-similar analysis and the present DNS model, we denote the wavenumber here as k = αh in contrast to α = αδ in the self-similar basis, so that k = α 1/H , where H = h/ δ. Note that k is a dimensionless wave number in the present basis not only but has also another meaning: it represents the inclination as defined in Sect. 3 (see Fig. 3).

Besides, if we consider the dependence between α, k and t, then α 2 =α 2 4D ˜ t ˜ or α 2 =4k 2 t and we can then obtain the following useful relation,

t = α 2 4k 2 . ( 50
)
The problem is solved for the following physical properties (the same viscosity has been taken for both fluids for the sake of comparison with the analytical solution): μ1 = μ2 = 1 × 10 -3 kg/(m s), D ≈ 10 -9 m 2 /s, ε2 = 2 × 10 -10 C 2 s 2 /(kg m 3 ), Δ Ṽ0 = 0.097 V. β = 2.25 is used in computations. These properties loosely correspond to water (ε 1 ) and ethanol (ε 2 ). The potential drop is varied within Δ Ṽ = 0-20 V. Thus, the dimensionless parameter Λ changes from 0 to 42,500.

The case of different viscosity will be used in Sect. 4.3.4 for sake of comparison with experimental data.

Numerical methodology

The numerical solution is based on the finite-difference method proposed in [49]. The space derivatives are approximated by a two-point difference scheme of the second order on a staggered grid. Spatial discretization with fine resolution leads to stiff problems and requires implicit methods for time advancement. Fully implicit methods produce a set of nonlinear coupled equations for the problem variables on the new time level, and are usually prohibitively costly for long-term calculations of the multi-parameter problems. Semi-implicit methods, in which only a part of the operator is treated implicitly, constitute a reasonable compromise for this class of problems. The semi-implicit third-order Runge-Kutta method [50] is used for time integration. The functions Φ and Ψ are found by a direct solution of the sparse system of linear algebraic equations, directly following from the discretization of the equations.

The infinite spatial domain is modelled by a finite domain that has dimensionless length l (so that waves with wavenumbers starting from k min =2π/lcan be captured in simulations). The condition that the solution at x →∞is bounded, is changed to periodic boundary conditions. The length of the domain has to be taken large enough to make the solution independent of the domain size. The value l =4π was typically chosen in most calculations, whereas l =8π was used to verify the results.

Results

Linear stage of evolution

For subcritical values of Λ and β, all perturbations in x-direction decay and the 1D solution is observed. However, this solution does not fully correspond to the selfsimilar solution due to the influence of the walls that breaks down the validity of the latter. The dimension-2 D the diffusion is the only one mechanism to mix the two fluids for the 1D solution.

Small perturbations are superimposed on the interface between the liquids in the numerical simulations. Two kinds of such perturbations are considered:

(a) The initial disturbances which are natural from the viewpoint of the experiment. The so-called "room perturbations" determining the external low-amplitude and broadband white noise are imposed on the surface.

(b) Artificial forced monochromatic perturbations with a fixed wavenumber α. These artificial perturbations allow better understanding of the behavior of the system.

For the supercritical parameters, a special kind of electro-hydrodynamical instability takes place. The instability triggers an additional mechanism of mixing, by advection, which is more powerful and strongly reduces the time for mixing. In all subsequent calculations, the initial conditions are defined as follows for both monochromatic and room perturbations,

w 0 = 1 2 (1 -sign (y)) +A 0 (1 -|tanh (ay)|) (sign (y)+s (x)) ,
where s (x) is either a uniformly distributed random value in the range [-1,+1], or a monochromatic sinusoidal wave. We define the following function (which will be also used in Sect. 4.3.2),

A(t) = max (y) max (x) w(t, x, y) -min (x)

w(t, x, y) .

A 0 = A (0) is a small parameter specifying the amplitude of disturbances, and 'a' controls the initial mixing layer thickness. Numerical simulations were performed with a = 99. Let us first consider simulations for the monochromatic perturbations as initial conditions. For such conditions, the initial perturbation is set up by a sinusoidal perturbation in x-direction with one single wavenumber k. These perturbations are artificial, but they allow us to understand the nature of the instability. So we return to the interpretation of the results of our selfsimilar analysis in Fig. 3. Keeping the wave number less time of the order of t ˜ ∼ 0.1 ˜h˜ still provides perfect matching between self-similar and DNS solutions when representing the permittivity as a function of the y-direction, see Fig. 4. The dimensionless time of full mixing is estimated of the order t ∼ 0.5. Evidently, k constant, originally along the straight line in Fig. 3, the disturbance decays first and after some time it may grow. As it can be seen from Fig. 5a even for the unstable parameters there is a latent period of time when perturbations do not grow up to a certain time. These observations are in good correspondence with the linear analysis of the self-similar solutions (see Fig. 5a).

The second kind of perturbations is the natural "room disturbances". In our case, all other perturbations practically decay in comparison with the most dangerous one with the wavenumber k r , see Fig. 5a. Thus, any evolution results in an almost monochromatic perturbation with the wavenumber k r . It is close to the wavenumber observed in DNS (Fig. 5b).

The results of the numerical simulation of the nonlinear system ( 43)-( 49) with the natural white noise initial conditions are in good correspondence with the linear stability results. Both are presented in Fig. 6. The disturbances are fully developed, t ≈ 0.14

Nonlinear evolution

The linear filtering mechanism singles out only one wave number k r with minimal time to instability, see Fig. 5a. When the amplitude of the perturbation becomes sufficiently large, the nonlinear stage of evolu-tion takes place.

For monochromatic disturbances, the function A(t) (defined in Sect. 4.3.1) is two times the amplitude of those disturbances. It is easy to see that A(t) is exactly zero for any 1D solution (including fully mixed state), and the range of its possible values is [0, 1]. The maximal value, 1, is reached when both pure liquids 1 and 2 are present at some y = const. The temporal evolu-tion of the amplitude as a function of its initial value and of the voltage is shown in Fig. 7. After the initial stage, which, in accordance to the linear stability anal-ysis, consists of a characteristic decay and followed by a growth, the disturbances reach saturation and then decay again. This decay seemingly justifies the picture in Fig. 2b, where all perturbations decay at sufficiently large times. In reality, the mechanisms differ. On the one hand, for a long-time process the perturbations are always nonlinear: the interface between the two phases becomes too bent for linearization to be valid. On the other hand, the diffusion smooths out the difference in concentration of the two liquids. The final stage of the evolution corresponds to the uniformly distributed (fully mixed) fluids.

If the potential difference is large enough, Λ > Λ * (approximately 5000), the amplitude can reach unity, i.e. the system contains areas filled with pure liquids 1 and 2 along some y-cross-section. Calculations show that these areas extend in "stripes" to almost the full microchannel width (Fig. 8e), and in this case the diffusion process occurs along x (Fig. 8f). For Λ < Λ * , transverse diffusion develops for sufficiently long time (see Fig. 9d ) and the concentration gradient does not reach sufficient high values to form "stripes" (Fig. 9e ).

Mixing efficiency

In order to estimate mixing efficiency, we introduce the separation amplitude,

A sep (t) = 1 x max w - 1 2 dxdy, (51)
enough, Λ > Λ * . As mentioned earlier, the value of Λ * is connected with the appearance of the "stripes", so these "stripes" are responsible for mixing enhancement over low-voltage regimes. However, with increasing Λ, this enhancement becomes small. The dynamics of the solution in time sequentially goes through four stages: I. A self-similar expansion of the diffusion layer according to Eq. ( 16) occurs. The solution does not depend on the strength of the external electric field and is determined only by the diffusion of the two layers. If we take the characteristic time associated with diffusion t ˜D = ˜h2 /D ˜ , then the curves in Fig. 10 shrink into one line and Eq. (51)for where the integration is performed over the whole simulation domain. This amplitude also varies from 0 to 1, where '0' means complete mixing and '1' means complete separation.

Our calculations (Fig. 10) show that the mixing efficiency practically does not depend on Λ unless it is large self-similar solution ( 16) can be written as,

A sep (t) = erf 1 2 √ t + 2 t π exp - 1 4t -1 .
∂w ∂t + U ∂w ∂x + V ∂w ∂y = 1 Λ ∂ 2 w ∂x 2 + ∂ 2 w ∂y 2 . ( 53
)
For Λ →∞, the solution of the system of equations does not depend on Λ and all curves (see Fig. 10) shrink into one. The diffusion is negligible for this case. IV. Diffusion of the formed regular structures. As the liquids are mixing at stage III, the "stripes" of the fraction w are eventually formed, where w practically does not depend on the y-coordinate, see Fig. 8. Also, dw/dy = 0 at the horizontal boudaries, so the velocity field vanishes, U = V = 1. Mixing is thus carried out only by the diffusion, hence the diffusion time t ˜D = ˜h2 /D ˜ must be taken as the characteristic one. Curves 1, 2 and 3 turn into parallel straight lines in semilogarithmic coordinates, see inset of Fig. 10a. It is possible to show it with trivial case. If we consider the diffusion problem and keep only the main harmonics in the series along x with wavenumber k,

w 0 = 1 2 (cos(kx) + 1), (54)
then this problem can be solve analytically and A sep takes the following form,

A sep = 2 π exp -k 2 t . ( 55
)
One can see that the dotted line in Fig. 10 is parallel to the others with non-zero Λ. It means that curves differ only by a multiplication constant, which transform into an additive constant in the logarithmic scale.

Mechanism of the instability

In order to understand the physical mechanism of the instability shown in Fig. 7, we consider a simplified theory. Note that there two competitive time-dependent phenomena in this process are: (i) the mixing length when the diffusion layer is expanding in the ỹ-direction as 2 Dt and (ii) the amplitude of the unstable perturbation in the x-direction which is increasing in time.

The first one depends only on the diffusivity (D ˜)and does not depend on the voltage (ΔV ˜ ) and the second depends strongly on ΔV ˜ and with its increase, the growth rate also increases. For sufficiently large ΔV ˜ ,the first characteristic time is much slower than the second one and, qualitatively, the instability can be considered with some frozen constant δ ˜. This is well illustrated in Fig. 7 in which one can see that the mixing layer is practically not changing with time while the instability runs into the non-linear phase of its evolution. The second assumption of this simplified theory is that we restrict ourselves with the long-wave perturbations.

The instability is connected with the ratio of permittivity ε ˜1/ε ˜2 which is taken as very large. Hence, our assumptions are the following, Let us consider the permittivity dependence as a piecewise linear function,

(a) Δ Ṽ Δ Ṽ0 , Δ Ṽ0 = D μ2 /ε 2 , (b) ∂/∂ x ∂/∂ ỹ, (c) ε1 /ε 2 1.
ε = ⎧ ⎨ ⎩ ε1 for -h < ỹ < Δ1 , ε1+ ε2 2 + ε1-ε2 ã (ỹ -ã) for Δ1 < ỹ < Δ2 , ε2 for Δ2 < ỹ < h. (56)
Δ Δ where 1 = -˜δ/2+a ˜(x, t ˜)and 2 = ˜δ/2+ a ˜(x, t ˜)(see Fig. 11 for details). Under the above assumptions, Eq. ( 8) can be integrated twice and the constants of integration can readily been found from the BCs, Eq. ( 12),

∂ Φ/∂ ỹ ≡ Ẽ = ε2 Δ Ṽ h B ε , ( 57
)
where B is a dimensionless constant of integration which is given by,

B = 1 - 1 2 δ h + δ h ln ε1 /ε 2 ε1 /ε 2 - ε1 ε2 ã(x, t) h -1
. (58)

The terms connected with velocity V ˜ and (∂Φ ˜/∂x ˜)2 can be omitted and the Maxwell pressure ˜Πcan be found from Eq. (10). Moreover, by simple algebra, we can show that, ˜Π2 / ˜Π1 =ε 1 /ε ˜2 1. In accordance with the hypothesis (c), ˜Π1 can be neglected in comparison with ˜Π2 . Hence, the pressure field ˜Π2 can be calculated from Eq. (10)andEq. (51) with the above assumptions, Π2 = 1 2

ε2 2 Δ Ṽ 2 h2 B 2 . ( 59
)
The Maxwell pressure Π2 is directly proportional to the squared voltage and inversely proportional to the Ẽ Fig. 11 Simplified sketch for the mechanism of instability of Fig. 7 width of the channel ˜h (the smaller is the width of the channel, the larger is ˜Π2 ). Note also that ˜Π2 (x, t ˜)has a maximum which coincides the maximum of a ˜(x, t ˜), the same applies with the minimum. Thus if the layer of liquid II becomes thinner because of some localized disturbance (as shown in Fig. 11), the electric field in this region becomes intense which, in turn, creates a high Maxwell pressure spot, see Fig. 11. This high-pressure spot drives the liquid away from it, the layer nearby the spot becomes very thin and the Maxwell pressure increases again, providing positive feedback. This mechanism of instability is expected to be stabi-lized by viscosity at short length scales. The behavior described above is perfectly confirmed with DNS for initial stages of evolution at Λ > Λ * , see Fig. 12a. However, whenever the assumption (56) fails (Λ is small or t is large), the resulting distribution becomes different, see Fig. 12b,c.

Conclusion

Mixing of dielectric miscible viscous liquids with different permittivities has been investigated analytically and numerically. Since the mixing layer is expanding with time, its instability is not described by the usual exponential law, but rather by the power law. This layer has been proven to be stable with respect to shortand long-wave disturbances. Intermediate waves have been found to be unstable for a finite period of time. The influence of walls does not allow to see the end of instability in practice, but the initial period, when all perturbations decay, has been confirmed with DNS. We have found by numerical simulations that with a strong enough electric field the amplitude of the waves grows, so that their edges reach the walls. Due to the Maxwell pressure between the top of the wave and the wall, the liquids reorganize into alternating stripes. This scenario of distinct stripes has not, to our knowledge, been identified previously, and our simplified theory has explained it. Mixing is enhanced with higher applied electric field as expected. However, the "stripes" of separated liquids may slow down the mixing, so moderate fields lead to better mixing in the long-time process. Moreover, we have verified that our results are consistent with experiments from [START_REF] Drazin | Hydrodynamic Stability[END_REF] with transformer and silicone oils with respect to the mixing index and time evolution of instability. The obtained results can be used in microfluidic applications as, for example, the control over on-chip assays that require rapid mixing of fluids.

Fig. 1

 1 Fig. 1 Schematic of the mixing channel, taken as the computational domain

Fig. 2

 2 Fig. 2 Growth rate λ vs the wavenumber α for a Λ = 1500, β = 2.25 and γ = 0 and different values of the dimensionless half-width of the channel H and, b Λ = 1500, H = 2 and different values of dimensionless viscosity γ

Fig. 3

 3 Fig. 3 Marginal stability curves for β = 2.25 and γ = 0 in coordinates α-1/H. The wavenumber α and 1/H are changing along the straight line as it is shown by the red arrows; λ during this evolution successively transfers between stable, unstable and again stable regions

Fig. 4

 4 Fig. 4 Evolution of the 1D permittivity solution for different time instants with β =2.25. Solid line stands for the DNS and dashed line, for the analytical solution (Eq. 27)

Fig. 5 aFig. 6

 56 Fig. 5 a Solid line stands for the marginal stability of the self-similar solution, triangles (DNS simulations) stand for the time when monochromatic perturbations with the wavenumber k change their decay to growth. The filled circles stand for the wavenumber k r realized for the room dis-

  turbances as intitial conditions, β =2.25 and Λ = 1500; b Permittivity at the cross-section (y = 0) for Λ = 1500, β =2.25, ε mix =1+β/2 is the permittivity of the mixture.

Fig. 7

 7 Fig. 7 Development of A (t) as a function of a initial amplitude A0 (1-A0 =10 -5 , 2-10 -4 , 3-10 -3 ,4-10 -2 ,5-10 -1 ) and for β =2.25 and Λ = 37000, bΛ: 1-Λ = 37,000, 2-16,500, 3-6100, 4-4000, 5-2500, 6-1500, 7-1000, A0 =10 -5 . The marked points are depicted in Fig. 8

Fig. 8

 8 Fig. 8 Temporal development of the fraction w for fluid 1 for β =2.25, A0 =10 -5 . a t =0, b 0.0071, c 0.0092, d 0.0149, e 0.0328, f 0.2. Snapshots (b-f) correspond to the markers in Fig. 7. Left column: Λ = 16,500, center column: Λ = 32,000, right column: Λ = 65,500

Fig. 9

 9 Fig. 9 Temporal development of the fraction w for fluid 1: β =2.25, Λ = 4000, A0 =10 -5 .The snapshots are marked in Fig. 7

  If the external field is strong enough, then the selfsimilar solution loses its stability. The initial disturbances grow until the velocity components become large enough. III. At some time t = t * , the disturbances become large enough to affect the convective terms and, hence, A sep .Fort>t * it is convenient to switch to a different characteristic time, namely, t ˜c = μ ˜ 2 . ε ˜2ΔV With this time, the parameter Λ disappears from the right hand-side of Eqs. (45-46) but appears on the right-hand side of the diffusion-convection equation of the fraction w (Eq. 43),

Fig. 10

 10 Fig. 10 Separation amplitude dynamics for different Λ and β = 2.25as a function a of time, and b Λt. 1: Λ = 16,500, 2: 13,000, 3: 8100, 4: 0. Dashed line corresponds to the self-

Fig. 12

 12 Fig. 12 Fraction (w) and Maxwell pressure distribution for Λ = 32,000 (b, e) and Λ = 4000 (e ). The snapshots are marked according to Fig. 7.

Table 1

 1 The function J for different values of β

	β	0.10	0.20	0.40	0.70	1.0	2.0	5.0	10.0
	J	-0.0035	-0.0121	-0.0383	-0.0870	-0.138	-0.289	-0.580	-0.830
		δ = 2 Dt : length in ỹ-direction;
		Ũ0 = D/ l0 : velocity in x-direction;
		Ṽ0 = D δ/ l2 0 : velocity in ỹ-direction; Π0 = D μ2 / δ2 : pressure;
		Δ Ṽ : drop of potential.

Acknowledgements

The investigation of the self-similar solution and its linear stability analysis were supported by the Russian Science Foundation, project N 20-79-00044, and the direct numerical simulation was supported by "Projets de Recherche Conjoints" (PRCCNRS). The research is carried out using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University.