Numerical Analysis of Density-Driven Reactive Flows in Hele-Shaw Cell Geometry

−In this paper, a two-dimensional numerical simulation of the unsteady state of a two non-isothermal immiscible liquids layer system filling a reactor formed by two closely spaced parallel glass sheets, which is called an Hele-Shaw cell, vertically oriented, with an expected neutralization reaction between an acid and a base in the lower layer, under the action of gravity, is studied. Attention is given on the general behavior of the complete temporal pattern evolution (velocity, temperature, and concentration profiles) and the identification of the exothermic reaction’s role in giving birth to chemo-hydrodynamic patterns that occur because of concentration gradients. The effects of gravity and changes in initial acid and base concentrations on the formed patterns were studied. The mathematical model governing the phenomenon was solved numerically by the CFD software package COMSOL Multiphysics, with the finite element method and a comparison with the experimental data was made. The results show that this numerical tool is promising for the understanding of the reactive instabilities happening when two immiscible fluids come into contact. Keywords-liquid-liquid reactive system; Hele-Shaw cell; diffusion; chemo hydrodynamic patterns


I. INTRODUCTION
Modeling of the coupling between chemical reactions and hydrodynamics is one of the most challenging and trending topics in fluid mechanics. This interest arises from the large impact of the reaction on mass transfer rate and also from the fact that the coupling of chemical kinetics with other transport processes is far from being fully understood [1]. Simple instabilities such as Rayleigh-Benard [2], Marangoni [3], and Taylor-Benard [4], have been well investigated, but the coupling of reactions and hydrodynamics is quite recent and incomplete, despite the fact that their mechanics occur in a very wide range of fields such as pollution decontamination, earth mantle dynamics, geological formations and CO 2 sequestration in soils. When two superposed immiscible non-reactive fluids are put in contact in a gravity field, only mass transfer is produced and causes thermal and solutal gradients [5]. Those gradients are responsible of the changes in density, viscosity, heat transfer, and surface tension of the fluids [6] and trigger multiple types of simple dynamics such as Rayleigh-Benard, Rayleigh-Taylor, diffusive layer convection, double diffusive fingering and Marangoni instabilities [7]. However, in the case of reactive systems, chemical reactions take place at or on the interface, interplay with those instabilities, alter patterns and give birth to complex dynamics [8]. In many cases, these processes are accompanied by heat release, and in turn hydrodynamic instabilities, driven by these gradients, can drastically alter the chemical reaction rate [8]. In such a situation, interactions between heat transfer and mass transfer, on one hand, and between the reaction itself and the diffusion of the reactants on the other, are inevitable [9]. Systems of two reacting fluids have been studied thoroughly, theoretically and numerically. The first work done in this field was that of Quincke in 1888 who reported periodic twitching of a bubble of air in water when a thin jet of alcohol was sprayed towards it [6]. To obtain an overview of all works done in this field the interested reader is referred to recent researches [6,7,9,11,12]. The main objective of this paper is the numerical detection of the temporal development of the flow structure (velocity, temperature, and concentration profiles) that occurs because of the coupling between diffusion and neutralization between a base and an acid putted in contact in a Hele-Shaw cell. An effort was given in detecting the effect of gravity and to exhibit changes in patterns obtained when the initial values of acid and base concentrations varied. This study investigates the chemohydrodynamic flows in an Hele-Shaw cell without considering the Marangoni effect, the Soret effect or the effects associated to double diffusion instability (since the acid and the base are assumed to have the same diffusivities in water). Finite element method was used by utilizing Comsol Multiphysics. The aim is to obtain an accurate numerical solution to the stated problem and to add a contribution to the experimental and numerical works done in this field.

A. Physical Domain
Our study is based on the experimental work of [10] in a Hele-Shaw cell. Figure 1(a) illustrates the 2D model geometry of height H and width L. The two vertical plates are separated by a Teflon foil, forming a gap of width d=1mm, too small compared to H and L. For that reason and for the sake of simplicity and to reduce computational cost, the fluid flow may be considered as quasi-two dimensional (Figure 1(b)). The cell is operated in the vertical position, filled by two superposed immiscible liquids separated by a supposed plane and nondeformable interface located at z=H/2. The upper fluid is an organic solution (with density ρ 02 ) of propionic acid dissolved in isobutyl alcohol. The lower one is an aqueous solution (with density ρ 01 ) of sodium hydroxide.

B. Governing Equations
The fluids are assumed to be Newtonian and incompressible and undergo steady laminar flow. The base is insoluble in the organic medium of the upper layer in a way that the fluids react via the second order mechanism A+B→C, only in the lower layer. The cell walls are supposed to be rigid, adiabatic and impermeable. The thermo physical properties are assumed to be constant (because of the low concentrations of the two species), except for the density in the buoyancy term. A linear law is considered: Subscript 1 refers to the base and 2 refers to the acid. T is the temperature, ߩ ଵ and ߩ ଶ are the reference densities at reference temperature T ref , Mathematically the motion inside a vertically orientated Hele-Shaw cell is analogous to a two-dimensional flow inside a porous medium [13]. Usually the Darcy-Boussinesq system of equations is used, so the governing equations for transient analyses include the continuity equation (3), the Navier-Stokes equations (4), the energy equation (5), and the convectiondiffusion equation (6). Under assumption of laminar and incompressible flow, invoking the Boussinesq approximation in modeling the buoyancy force, and using the Hele-Shaw approximation in velocity term, the governing equations in each layer, can be written as: In the above equations, U ሬ ሬԦ ሺu, vሻ is two-dimensional

C. Initial Conditions
The concentrations of the acid and the base are

D. Boundary Conditions
At the impermeable but frictionless side walls (x=0, x=L, y=0 and y=H): u=v=0 and ൌ 0 . For 0x ‫ܮ‬ and y=0.071m (the interface), we impose continuity of these variables: . k 1 and k 2 denote the heat conductivities, D 1 and D 2 are the diffusivity coefficients, and ߤ ଵ , ߤ ଶ are the dynamic viscosities of fluids. Our aim is the detection of chemo-hydrodynamics patterns, their rate of growth, and the effects of gravity and initial acid and base concentrations on the structure of the flows formed. These results are compared to the previous experimental work in [10] and the two-dimensional numerical analysis recently derived in [7], showing that the agreement between the experimental results and numerical predictions is obtained.

III. NUMERICAL METHOD AND VALIDATION
The above assumptions and the (3)-(6) along with the boundary conditions were solved by the finite element CFD package Comsol Multiphysics version 4.3. Calculations were performed on a non-uniform grid with 8638 elements. A refinement of mesh at the interface and at the zones near the interface was done (Figure 2(b)). The numerical procedure was validated by a comparison with the experimental data provided by [14] concerning the experience of Marangoni convection in a cavity filled with acetone at a height h and subject to a horizontal temperature gradient ∆T for the case of low Marangoni and Rayleigh numbers. Figure 2 shows similarities and a closer behavior with the measurements, and provides confidence to the accuracy of CFD analysis, so we believe it can be safely used to understand the dynamics of instabilities occurring between two superposed reactive fluids.

IV.
RESULTS AND DISCUSSION Chemical reactions can change the density of a solution either by modifying the total volume or by releasing or consuming heat. These two effects provide a solutal and/or thermal contribution to the total density change [15]. Here, we focus on density changes induced by gradients of concentrations, which means that the dynamics here are mainly controlled by solutal effects, because it has been proven that the temperature difference over time and space is insignificant [10]. Ιndeed, heat diffuses faster than mass, the density gradients due to heat are quantitatively negligible [10]. Below we will present the numerical simulation obtained results. This code uses the finite element method. In the first instance we will show results related to the temporal evolution of patterns formed, keeping the Damkohler number Da (defined as the ratio of the hydrodynamic to the chemical time scale) and the initial concentrations ratio constant. Then we present the effects of gravity and the change in initial acid and base concentrations (C 0A , C 0B ) on the structure of the resulting flow.

E. Temporal Evolution of Chemo Hydrodynamic Patterns
In order to observe and to determine the evolution of chemo-hydrodynamic patterns formed in time, a series of numerical simulations at different times after onset of instability has been conducted for Da=1, and the ratio of initial concentrations C 0A /C 0B =1 (C 0A =C 0B =1mol/l) kept unchanged. Figure 3 shows a set of typical simulation streamlines at different times. As seen in [16], soon after contact, the acid enters into the basic solution ( Figure 4) and after a few seconds (≈2s) a downward moving front is obtained (Figure 3). Mass transfer and chemical reaction are taking place simultaneously, leading to instabilities in the system. In the zone slightly above the contact line, in earlier times and for a few seconds, the structure of the flow is not well determined (Figure 3(a)-3(c)). In the lower phase, small roll cells are observed confined between the interface and the reaction front (≈13mm) as seen in Figure 4. In the upper layer, convection arises because the fluid layers on the top of the Hele-Shaw cell are rich in acid and the fluid layers just above the interface are characterized by a depletion of acid concentration [18], because of its continuous diffusion across the horizontal contact line towards the lower layer ( Figure 5). Also, the fluid is warmer ( Figure 6) because of the heat released by the reaction between the acid and the base just underneath the interface. This unstable stratification leads to a concentration gradient in the vertical (y) direction, and diffusion will carry the solute from the bottom to the top, in a process called diffusive layer convection mechanism [7,10,15]. In the lower layer, as Da=1, the hydrodynamic and chemical effects are in equality. The fresh acid diffuses continuously through the interface and a part of it will be swiftly consumed by the neutralization reaction which takes place as soon as the acid comes in contact with the base near the interface, releasing heat and causing a warmer zone near the interface as observed in Figure 6, despite its weak effect [18], leads to destabilizing transient concentration gradients (Figure 7), that in turn triggers Rayleigh-Benard instability [10,12]. This allows the generation of regularly organized cellular structure forms at some distance from the interface penetrating into the bulk. The cells are confined between the liquid-liquid interface and the reaction front [8]. This front carries the other part of acid in the bulk of the lower layer and will react with the base and form irregular fingers sinking away from the reaction zone (Figure 3).  The blue color indicates a low acid concentration in the lower layer. This is obvious as long as the reaction consumes www.etasr.com Samia & Kadja: Numerical Analysis of Density-Driven Reactive Flows in Hele-Shaw Cell Geometry the majority of the acid diffused through the interface. In Figure 8, (0.071, 0.071) is a point at the interface, (0.071, 0.075) is a point at a distance of 4mm above the interface, and (0.071, 0.065) is a point at a distance of 6mm below the interface. The acid concentration profiles at those points give the proof that the reaction consumes most of the acid diffusing in the lower phase.  Figure 9 shows the velocity magnitude plotted as a function of time at four points in the lower phase: (071, 0.070), (0.071, 0.069), (0.071, 0.068), and (0.071, 0.065) in 1mm, 2mm, 3mm, and 6mm distance from the interface respectively. The velocity profiles are checked at those points because they represent the region between the interface and the reaction front. The profiles give a very clear idea about the dynamics in this area. As it is illustrated in Figure 4, there are two distinct regions: the region just underneath the interface where strong dynamics are noticed and the region of the bulk phase where dynamics are weak. This is consistent with the profiles in Figure 9, where it is clear that the velocity reaches its maximum just below the interface and weakens in further distances.  Figure 10 shows that gravity force has an important role in accelerating the occurrence of instabilities below the interface and also accelerates the appearance of the established state (t=2s). This is obvious as long as it can induce or affect buoyancy-driven instabilities in the cell [19].

I. Effect of Initial Acid and Base Concentrations
The effect of the initial acid and base concentrations on the patterns formed in the gravity field is inspected because it is the main factor changed in the experiments. Two important cases are considered: The case where the acid is in excess (when the initial concentration of acid is higher than that of the base) and the case where the base is in excess (the opposite). Figure 11 shows that the behavior of instability in the zone underneath the interface changes in a small amount qualitatively from case 1 to case 2. Considering case 1, the more the ratio (C A0 /C B0 ) increases, the less hydrodynamics occur near the interface, and the number of roll cells is reduced to two larger ones. Because the acid is in excess, it does not meet enough base to react (C 0B =0.5mol/l) and diffuses in the bulk of the lower layer. This acid reacts with the base which allows the appearance of strong dynamics in the form of two sinking irregular fingers ( Figure  11(c)). This result is in close agreement with the experimental results in [9]. Considering case 2, where the base is in excess, we see strong dynamics, in the form of small roll vortices confined near the interface. Because of the low acid concentration, when it diffuses to the lower layer, it will be consumed quickly near the interface by the neutralization reaction. The behavior of the fluid in Figure 12 is identical to that obtained in Figure 11. Maximum of velocity in case 2 is reached in a point in the lower phase, 6mm far from the interface.
V. CONCLUSION In this work, the coupling between a chemical reaction and hydrodynamic instabilities occurring in a Hele-Shaw cell containing two vertically superposed immiscible fluids has been numerically simulated, using a simplified mathematical model and taking into account the fact of the very closely spaced plates by considering the Hele-Shaw approximation. The goal was to detect the patterns formed and their temporal development. The results allow us to deduce that: • The diffusion of the acid from the upper layer into the lower layer is the main factor behind all dynamics. In fact, diffusion is the engine of the chemical reaction between the acid and the base. Also it causes the diffusive layer convection in the upper layer and the fingers appearing in the lower layer.
• Convective patterns appear quite quickly as soon as the acid comes in contact with the base, after less than one minute.
• The gravity force also plays a major role in accelerating the instabilities once the reaction takes place.
• The increase of the initial acid to base concentration ratio which means an excess in acid does not favorite the formation of chemo convective patterns and therefore this breaks the reaction front, however strong dynamics occur in the bulk of the lower layer. In contrast to this, a decrease in this ratio, which means an excess in base allows us to see strong dynamics underneath the interface and elongates the reaction zone.