VALIDITY OF THE REYNOLDS EQUATION FOR MISCIBLE FLUIDS IN MICROCHANNELS

In this paper, we consider asymptotic models for miscible flows in microchannels. The characteristics of the flows in microfluidics imply that usually the Hele-Shaw approximation is valid. We present asymptotic models in the Hele-Shaw regime for flows of miscible fluids in a channel in the case where the bottom and the top of the channels have been modified in two different ways. The first case concerns a flat bottom with slip boundary conditions obtained by chemical patterning. The second one is a non-flat bottom with a non-slipping surface. We derive in both cases 2.5D and 2D asymptotic models. We prove global well-posedness of the 2D model. We also prove that both approaches are asymptotically equivalent in the Hele-Shaw regime and we present direct 3D simulations showing that for passive mixing strategy, the Hele-Shaw approximation is not valid anymore.


MATHIEU COLIN, THIERRY COLIN AND JULIEN DAMBRINE
The aim of this paper is to study asymptotic models and their validity for mixing in microchannels.Many approaches have been considered in order to enhance fluid mixing.Active mixing seems to be the most efficient approach but very difficult to handle for practical applications since one needs to build sophisticated devices.Devices for passive mixing are easier to obtain but the efficiency of such process is not clear.In this paper, we propose to investigate numerically some passive mixing strategies in the Hele-Shaw limit.We consider fluid flows that are driven only by the pressure gradient and no external constraint on the system are taken into account.Two kinds of mixing strategies are considered.Both involve modification of the surface of the top and the bottom of the microchannels.The first one uses a flat bottom channel with slip boundary conditions obtained from chemical patterning ( [13], [14]); the second one uses a non-flat bottom channel with a non-slipping surface ( [21], [22]).It is almost impossible to compare these mixing strategies by making direct 3D numerical computations.We propose here another approach relying on the Hele-Shaw approximation.It is an asymptotic process in which the flow is finally governed by the Reynolds equation (that is similar to Darcy's law for flows in porous media), see for example [1].We show that, in this limit, it is possible to compare the efficiency of the mixing for both strategies: in this particular regime, they are equivalent in the sense that, given a non-flat bottom, one can find a patterning that leads to an equivalent asymptotic model.Nevertheless, we also present some numerical experiments proving that this result cannot be extrapolated to general flows in microchannels (that is, outside the range of validity of the Hele-Shaw approximation).More precisely, we show by a direct 3D computation that in the case of a non-flat bottom, the shapes of the Hele-Shaw solution and the 3D flow are very different and therefore, the asymptotic regimes are no longer valid.Note that we are not in the configuration where the Hele-Shaw approximation is used to study Saffman-Taylor fingering like in [18].Indeed, in our case, the fluids are both injected at the inlet of the channel as shown on Figure (1).Moreover, the microfluidic flows considered in this paper are characterized by: i) the jump of viscosity coefficient, ii) the effect of confinement, iii) the influence of the ratio of the flow-rates (pressure driven flows), iv) the absence of gravity effects, v) the effective mixing properties.Therefore, some phenomena like instability due to gravity (see [7]) or shape interfaces and viscous fingering (see [15]) are not present in our case.Moreover, we are not concerned by lubricant effects due to roughness as studied in [4] where the coupling between the size of the roughness and the thickness of the channel is investigated.One can see [9] for rigorous results concerning rough boundaries.
1.2.The basic equations.Usually, the flow dynamics in channels are governed by the Navier-Stokes equations.Nevertheless, the small dimensions involved in microfluidic applications lead to consider flows at low Reynolds numbers.As a consequence, the inertia terms in the Navier-Stokes model can be neglected in a first approximation.In this context, the velocity of the fluid is therefore given by the Stokes equations: where U denotes the velocity of the fluid flow, P the pressure, D(U ) is the strain rate tensor given by D(U ) = 1 2 (∇U + ∇U t ).In order to describe a mixing of two fluids, we introduce an order parameter ϕ (see for example [11]), which describes locally the volume fraction of each fluid in the mixture.It satisfies ϕ ∈ [0, 1], ϕ = 1 in the fluid 1, ϕ = 0 in the fluid 2. Equipped with this parameter, it is possible to write a one fluid formulation for the mixing by using ( 1.1).The viscosity η is given as a function of η 1 ,η 2 (each fluid viscosity) and depends on ϕ.There are many different ways to describe η(ϕ) (linear, cubic, tangential, non-monotonic [19]).In this paper we consider the linear approach: η(ϕ) = η 1 ϕ + η 2 (1 − ϕ).The evolution of ϕ is given by a convection-diffusion equation: ) where D(ϕ) is the diffusion coefficient from one fluid into the other one.
In order to close the system given by ( 1.1) and ( 1.2), we need to introduce suitable boundary conditions.Let us first present the geometry of the channel we consider here: . Schematic view of a microchannel and its characterisic dimensions.
h is the height of the channel, and for a non-constant height, we denote by h the maximum height and δh the amplitude of the variations.l denotes the width of the channel and L its length (see Figure 2).For the Stokes equations ( 1.1), we impose a Dirichlet boundary condition U = U 0 (y, z) at the inlet of the channel and U = 0 on the lateral walls; at the outlet we force the pressure to be P = 0. On the top and the bottom of the microchannel, we use various boundary conditions based upon the different approaches studied for mixing.We consider two types of channels.One is a flat bottom channel with slip boundary conditions and the other one is a non-flat bottom with a non-slipping surface (see Section 2 for a complete description).For the convection-diffusion equation, we consider a Dirichlet boundary condition at the inlet ϕ = ϕ 0 (y, z) and ∂ϕ ∂ν = 0 everywhere else where ν denotes the outward normal of the domain.
One of our motivations is to introduce a reliable model in order to perform quick numerical simulations, using finite volume method.To this end, since the full 3D Stokes system is expensive in terms of computations, we present reduced models, based on the Hele-Shaw approximation.It consists in introducing two different scales in space by assuming that the variations in the vertical direction are much more rapid than that in the horizontal plane.Therefore, the integration of the Stokes equations ( 1.1) with respect to z (height) leads, after transformations, to a 2-D equation on the pressure P. The velocity is recovered by a Darcy law.The convection-diffusion Equation ( 1.2) is writing respectively in 3-D and 2-D providing two different models: the first one is called the 2.5 D model while the second one is a complete 2-D model.In this paper, the validity of the 2D system is discussed for various situations.We consider first viscous displacements involving simple geometry with no surface modification on the microchannel.The second case concerns different mixing strategies involving boundary conditions already studied in literature.
1.3.Outline of the paper and main results.In Section 2, we perform explicetly the Hele-Shaw approximation in the case of a flat channel with chemical patterning and in the case of a non-flat bottom.We first introduce 2.5D models.In both cases, the limit system is given as follows.The pressure P (t, x, y) satisfies a 2D elliptic equation ∇ xy • K 2D j (ϕ)∇ xy P (x, y) = 0, with j = n or p (the index n stands for the case of a non-flat bottom while p stands for the case of a flat bottom with chemical patterning), see ( 2.25) or ( 2.33), where K 2D n (ϕ) or K 2D p (ϕ) are functions of (t, x, y) constructed from ϕ(t, x, y) (see Section 2 for the formulas).From this pressure field, one first build the horizontal components u(t, x, y, z) and v(t, x, y, z) of the velocity by for j = n, p, see ( 2.26) or ( 2.34) where K 3D j (ϕ) are functions of (t, x, y, z).The vertical component of the velocity is recovered finally through the incompressibility condition by

27) and ( 2.35).
In both cases, the order parameter ϕ solves the 3D convection-diffusion equation In Section 3, we propose some numerical computations.We first describe the numerical scheme that is used and then we perform some simulations in order to explain co-flow displacements and viscous fingering.Using the coefficients K 2D n and K 2D p , we show in what sense the two mixing strategies are equivalent and we explain how one can construct a patterning that leads to the same mixing than a given flat bottom in the Hele-Shaw regime.Then, in order to investigate the validity of the Hele-Shaw approximation, we compare numerically the 2.5D model with a full 3D computation in the case of a staggered herringbone micromixer.We show that, in this case, the Hele-Shaw approximation is not valid, therefore 3D computations are needed.
The appendix is devoted to the construction of a global, weak solution to the 2D model and to local in time strong solutions.The difficulty relies on the particular boundary conditions considered in this paper and for example, one has to prove that for any regular solution of the system, the velocity at the outlet of the channel is directed toward the exterior of the channel, which is not surprising from the physical point of view (see Lemma 5.2).
Notations: As usual, for a bounded domain Ω of R d , we denote by L p (Ω) the Lebesgue space where We define the Sobolev space H m (Ω) as follows The corresponding norm is denoted || • || H m .Let X be a Banach space.We denote by L p 0, T ; X (1 ≤ p ≤ +∞) the space of functions m : (0, T ) −→ X such that m is measurable and Different positive constants might be denoted by the same letter C.

2.
Obtention of the different models.The Hele-Shaw approximation relies on the hypothesis that the characteristic vertical length is smaller than the horizontal ones.The aim of this section is to construct some models derived from this approximation starting from ( 1.1)-( 1.2).Our strategy is somehow classical and we can refer to [11] for general theory of mixture as well as for the stability of an interface in a Hele-Shaw cell.Here the context is quite different: the gravity can be neglected and we are in a situation where the interface is stable.Our derivation of the models is classical, but we need the precise values of the coefficients and for the sake of completeness, we give the computations.In order to take into account that the velocity variations in the vertical direction are very rapid compared to the horizontal ones, we introduce two different space scales, one for the vertical direction z and the other one for the horizontal variables (x, y).Recalling that a microfluidic channel is a simple rectangular domain, we assume that h ≪ ℓ ≪ L (see Figure 2).It is then natural to perform the change of variable where ε is a small parameter describing the ratio between the height h and the others dimensions (L, ℓ) of the channel.We plug ( 2.3) into ( 1.1) and ( 1.2).The incompressibility condition ( 1.1) on the velocity U = (u, v, w) gives Equation ( 2.4) implies that the third component of the velocity w has to be of order ε.We hence need to perform the following change of unknown: Hence, the stationary Stokes system ( 1.1) can be written as In view of ( 2.3) and ( 2.5), the system ( 2.6) becomes Then Equations ( 2.10) and ( 2.11) can be rewritten as: According to ( 2.12), P is independent of z.Then, integrating along the z direction, we get: where A is a scalar function of (x, y) generated by the integration along z.In the general case, η depends on the local composition ϕ(t, x, y, z) of the mixing and Equation ( 2.13) becomes Integrating again Equation ( 2.14) along the z direction, we obtain where B is a function of (x, y) generated by the integration along z.Equation ( 2.15) can be understood as a Darcy law.As already mentioned, in this paper, for practical applications, we consider two different cases.The first one deals with a non-flat bottom with no-slip boundary conditions.The second one deals with a flat bottom with a chemical patterning of the surface inducing slip for the flow.We now derive the complete system in both cases.
• Non-flat bottom with no-slip boundary conditions.In this case, the boundary conditions for the velocity read This provides B(x, y) = 0, by taking z = 0 in Equation ( 2.15).Taking now z = h, we obtain (2.17) Denoting by we rewrite Equation ( 2.15) as Rewriting the incompressibility condition ( 1.1), and taking ( 2.15) into account, we obtain (2.20) By integrating ( 2.20) from 0 to h we get which provides denoting by the following equation on P The vertical velocity w(x, y, z) is recovered by integrating ( 2.20 By performing the change of variables ( 2.3) in the 3D convection-diffusion Equation ( 1.2), one gets: The complete model then reads (2.27) • Slip patterned flat channel.We consider the case of a flat channel with a chemical patterning on the top and the bottom.The classical slip condition reads , where U τ represents the tangential velocity, and L a characteristic length for the slip phenomenon.In our case, we define the top and bottom slip length patterns L t (x, y) and L b (x, y).The tangential velocity at the top and the bottom of a flat channel is the velocity along the (x, y) direction and the normal derivative is the derivative along the z direction.The slip conditions are: (2.28) Taking respectively z = 0 and z = h in Equation ( 2.14) and ( 2.15), we have . (2.30) From ( 2.29) and ( 2.30), one obtains (2.31) the complete model reads 2. Simplified 2D models.In this section, we introduced 2D simplified models as follows.As ε → 0, Equation ( 2.32) leads to ∂ z ϕ = 0 at the first order, and therefore, ϕ is independent of z.Averaging ( 2.32) over [0, h] and using ∂ z ϕ = 0 at z = 0 and z = h gives We now explain how one can compute this averaged velocity in both cases.
• 2D model for non-flat bottom.Since ϕ is independent of z, we obtain from ( 2.17) (2.36) In this case, we need to compute From ( 2.24)-( 2.27), we get the following 2D model • 2D model for a slip patterned flat channel.Let us rewrite ( 2.31) assuming that ϕ and therefore η is independent of z In the same way as above, we obtain which leads to the following model In addition, if we choose the same slip length patterns on the top and bottom of the channel (i.e.L t = L b = L) we get much simpler expressions for the permeability coefficients K 2D p (x, y) and K 3D p (x, y): Of course we recover in this case the formulas given in [8] for example.
3. Numerical methods.The aim of this section is to present an efficient numerical scheme in 2D and 3D for the different systems introduced in Section 2. First recall the general form of the system we deal with Here U represents the fluid velocity, this quantity being respectively 2D as in Section 2.2, U = (u, v) or 3D as in Section 2.1, U = (u, v, w).X represents the space variable ( (x, y) in 2D or (x, y, z) in 3D).In the Darcy law ( 3.44), the function F depends on the model we consider and more particularly on the geometry of the channel (slip patterns or non flat bottom).We use a very simple discretization of the problem that we present rapidly below.
3.1.Semi-discretization in time.Introducing t n = n∆t, we denote ϕ n = ϕ(t n , X), U n = U (t n , X), P n = P (t n , X).Since Equations ( 3.43) and ( 3.44) are stationary, their time discretizations are straightforward and reads In Equation ( 3.45), we use a classical finite difference scheme The transport part is treated by an explicit scheme while a semi-implicit θ-scheme is used for the diffusion part: Then, the semi-discretized system in time can be rewritten as: In this section, we consider a regular cartesian mesh, which is well-adapted when one deals with microfluidic flows, x i = i∆x, y j = j∆y in 2D and z k = k∆z in 3D, where i = 1, ..., N x , j = 1, ..., N y , k = 1, ..., N z .N x , N y , N z represent the mesh resolution.In general, for any variable V , we denote: ,k and we use a staggered grid.More precisely, the velocity nodes are placed in the following way (see Figure 3) whereas the pressure is computed on the nodes of the mesh P i,j,k = P (x i , y j , z k ).We now present briefly the essential point of our spatial discretization.We choose a finite volume formulation.The unit cell is given by ].The equation ( 3.48) include a transport term U n • ∇ which is approached with a Weighted Essentially Non-Oscillatory scheme (WENO 5) which reduces numerical diffusion for smooth enough solutions.This transport scheme is stable under the CFL condition (see [16]): We still have to give the discretization of the elliptic parts of System ( 3.46)-( 3.48), which can be written in the following form ( 3.48).Recall that the finite volume methods consists in computing the circulation of K(x, y)∇ xy P • n around each unit cell Ω ijk of the mesh, we obtain from Equation ( 3.46), denoting Γ ijk the boundary of Ω ijk Pointing out that both P and ϕ are calculated in the middle of each unit cell, ∇P is naturally evaluated on each part of Γ ijk while the coefficient K 2D n (ϕ) is computed by interpolation.We propose here to use the harmonic mean value on each lattices of the boundary, keeping in mind that this process conserves the flux continuity.Then introducing: we then obtain from ( 3.49): The same procedure is applied on Equation ( 3.48).The discretization of the other terms is classical and we refer to [5] for more details.
3.3.Viscous effects.The aim of this section is to investigate the validity of the models introduced in Section 2 with two classical test case dealing with viscous effects: the interface displacement and the viscous fingering.

3.3.1.
Co-flow displacements.The firts test case concerns the mixing of two fluids injected in a Y-shaped microchannel (Figure 4).Since we consider only the horizontal displacements, we choose the 2D model ( 2.37)-( 2.39) equipped with a linear viscosity law.We consider a flat geometry, that is h(x, y) is constant, and a no-slip condition at the top and the bottom of the channel.The diffusion coefficient is also supposed to be constant.It is usually accepted that the displacement of the interface is influenced by the viscosity contrast defined in the following way Of course, since one deals with miscible fluids, it is not easy to characterize the position of the interface since the notion of interface is not clear by itself!For that reason, in order to measure the interface displacement along the cross section direction we use a first order moment including the first order derivative with respect to y of the local composition parameter ϕ y∂ y ϕ(x, y, t)dy.
In the following simulations, we expect the most viscous fluid (the one with the biggest hydraulic resistance) to be pushed by the less viscous fluid in the cross direction of the channel.We also expect that the walls have an influence on this displacement.
For a low viscosity contrasts (λ ≪ 1) this influence should be neglectible, and hence the interface displacement M 1 should depend linearly on λ.On the other hand, if we choose a larger viscosity contrast, the walls effects should be the most important . From the physical point of view, one can say that at the end of the channel, where the viscous displacement process is over, the flow-rates in each fluid are equal, i.e.
where u 1 and u 2 are the horizontal components of the velocity in the fluids 1 and 2. l 1 and l 2 represent the width occupied by the fluid 1 and 2 in the channel.These lengths can be expressed as a function of the displacement of the interface M 1 : Fig 5).Thanks to the Darcy law ( 2.15), we have (3.52) Using ( 3.51) and ( 3.52) we get At the end of the displacement process, ∂ x P is the same in each fluid.We hence have Using the expressions of l 1 and l 2 and ( 3.53) we finally obtain , and according to ( 3.54), the relative displacement M 1 l of a fluid due to a viscous effect should be predicted by . This assertion is checked in the following numerical simulations.Before presenting our simulations, we ensure the discretization error to be small enough for the comparisons we wish to do.In order to find the mesh resolution which would fit this criteria, we need to perform a convergence calculus.Therefore, we compare the results given on the quantity M 1 (L, t) by simulations with the physical conditions described above, for different mesh resolutions (Figure 6).The convergence on M 1 is rather quick.This is due to the fact that M 1 , which only gives an information on the position of the middle of the interface, is not affected by any numerical diffusion.As we can see on Figure 6, an eligible resolution in this case would be 100 * 80.As we can see on Figure 7 (left), as the contrast λ between the viscosities η 1 and η 2 increases, the position of the interface M 1 (X, T ) gets closer to the borders.As foreseen, if the viscosity contrast is big enough (λ ≫ 1), the displacement tends to be independent of λ.Moreover, a linear behavior on M 1 for little viscosity contrasts seems to appear.It means that for little displacements, the interface is not influenced by presence of the walls.It seems that this behavior tends to disappear when λ > 1 (i.e.viscosity variations greater than 100%).We can see on Figure 7 (right) that simulations for different values of f (η 1 , η 2 ) (squares) nicely fits the theoretical displacement as described in ( 3.54) (black line) which validates the model in this case.Differences seem to occur as |f (η 1 , η 2 )| grows.It can be related to the wall effect that can be noticed on Figure 6 (left) for t ∼ 0.4.
All these results can hence lead us to think that the model used here works well for the co-flow interdiffusion case, see [6] for comparison with the experiments.

Viscous fingering.
In this section we study an other example of viscous displacement that involves the z direction: the viscous fingering.The experiment is the following (see Figure 8): a fluid 1 with a viscosity η 1 is injected in a channel containing a fluid 2 with a viscosity η 2 .In this study, we are interested in the profile of the interface in the ( x, z) plane.The flow along the z direction is a Poiseuille flow which tends to form finger shaped interfaces influenced by the viscosities η 1 and η 2 .The suitable 2.5D model fitted to this experiment is composed by Equations ( 2.24)-( 2.27) with a flat geometry and no-slip boundary conditions.The channel geometry chosen here is L = 2000µm, W = 200µm, h = 80µm.The profiles shown in Figure 9 correspond to two different viscosity gaps (η 1 > η 2 and η 1 < η 2 ).Those profiles seem to be physically relevant, however, further validation from the experimental point of view is needed in this case that can be done using advanced imaging techniques such as confocal microscopy, as shown in [10].

3.4.1.
Equivalence of the two models.In this section, we investigate the influence of the boundary conditions in the process of mixing fluids.We consider two types of boundary conditions: the slip ones and the no-slip ones.For simplicity, we assume that the viscosity η is constant and we limit our study to the 2D models presented in Section 2.2.The main goal is to make a connection between the two approaches used for mixing.Observe first that choosing h(x, y) = constant in ( 2.37)-( 2.39) and L t = L b = 0 in ( 2.40)-( 2.42) leads to two strictly equivalent models.The problem discussed here is the following: given a profile h(x, y), can we choose the slip lengths L t and L b so that the two 2D models are very close.
Recall first that (3.55) In this particular case, we supposed that the slip lengths L are the same at the top and bottom of the channel for K 2D p and K 3D p .We denote by h the height of the channel in the slip pattern model.The variable length in the no slip case shall be expressed as a function of h.We set h(x, y) = h + h(x, y).We want to find a couple of h(x, y) and L(x, y) that would make the results on the velocities U p and U n given by the two different models as close as possible.
Assuming that h(x, y) is known, let us try to find L(x, y) such that By multiplying the second equation by h and identifying the left hand side of each equation, we get which implies h(x, y) = 0, which is not relevant for the application we have in mind.However, there exists a choice of parameters h(x, y) and L(x, y) such that, under the constraint h = 0, the results obtained by the two models are very close.Indeed, assuming h(x, y) is known, let us minimize the norm of the following functional .
From ( 3.55), we have In each point of the domain, let us find L(x, y) that minimizes ||F x,y ( L)||, where ||.|| denotes the euclidian norm.By denoting a(s) = ||F x,y (s)||, we have Since ( h s − c 2 ) 2 = 0 for any acceptable s as seen above, it is possible to compute the derivative of a The optimum s * satisfies da ds (s * ) = 0, That gives (3.56) From ( 3.56), the choice of L that minimizes the error between the coefficients of the two models in every point (x, y) of the domain is hence It is then possible to state the following proposition.
Proposition 1.Let h be the height of the channel.For a given function h(x, y) = h + h(x, y), there exists a slip length L(x, y) which minimizes the error between the coefficients of the two models ( 2.37)-( 2.39) and ( 2.40)-( 2.42).Furthermore, L is given by Given a channel geometry h(x, y) (see Figure 10) and the choice of L(x, y) described in ( 3.57), we will measure the error on the results given by both models through numerical simulations.The channels geometries are: • For the no-slip case.
The height variations will vary from 0 to 50% of the total height.

Example of a herringbone micromixer geometry
• For the patterned case.
We use here the same geometry as described above, assuming h is constant.The slip length pattern L(x, y) is calculated as described in ( 3.57).As said before, we consider only one fluid in the channel.The viscosity and entry velocity are chosen such that Re = ρU l η ∼ 10 −9 .
The discretization steps are dx = 2.6µm and dy = 1µm which means 20 cells per ridge in the non-planar geometry case.We ensured the discretization steps to be small enough so that the numerical error is reasonable for the comparisons we wish to do.
We define the error between the results of the two models as Let us also introduce β = h h the normalized height variation.As we can see in Figure 11, the results on the velocity given by the two models are fairly close for the range of β we have chosen.At most, the error is around 20% for β = 0.5.In practice we never use such height variations (1/2 of total height).In literature, most of the time β ∼ 0.2 [21], which, here, generates an error lower than 10%.We can reasonably say that under the Hele-Shaw approximation, the two classical approaches for mixing enhancement studied are close.In order to study the validity of the 2.5D model for non planar geometry (Section 2), we compare numerical simulations performed on this model with predictions given by the Stokes model.The geometry used here is a ridgeshaped geometry, the height of the channel being given by where h max is the maximum height of the channel.First, we quantify the influence of h max on the error between the original Stokes model ( 1.1) and the 2.5D Reynolds model ( 2.25)-( 2.27) adapted to the non-flat bottom geometry.The error is computed on the velocities given by numerical simulations performed on both models on a N x × N y × N z = 400 × 100 × 50 grid.The following dimensions are taken for the channel length L = 2.5 mm, width l = 500µm, and h max varies between 50 µm and 200 µm.
The Figure 12 shows the expected convergence of the two models as h max tends to 0. The rate of convergence shown here is 1/2.The qualitative differences between the simulations performed on the two models are shown on Figure 13, by examining streamlines starting at the inlet of the channel.We notice that, for a sufficiently large h max , the streamlines computed from the Stokes model show a deviation of the flow from the main axis of the channel to the left wall.This deviation is not seen on the Reynolds model, even though the non-flat geometry is taken into account.This deviation phenomenon is a key mechanism in fluid mixing improvement by means of microchannel geometry modification (see [17] for a complete review on micromixers).As a consequence, the Reynolds model studied in this paper is not appropriate for the study of microfluidic mixers based on this approach.

Conclusion.
In this article, we have shown that the Hele-Shaw approximation (Section 2) applied to a 3D Stokes model can lead to various reduced models (Section 2.2).The main advantage of using such models is that numerical calculations are much quicker and simpler than for the full 3D Stokes model.Moreover, even though strong hypothesis are done on the microchannel geometry, these models allow us to take into account complex boundary conditions such as slippatterned surfaces or non-planar geometries.Those kind of boundary conditions have already been involved in various studies of enhanced fluid mixing in microchannels [21], [22] , [13], [14].For those reasons the reduced models developed here appear to be very interesting in general for the numerical study of flows in microchannels.
We also have shown that through this model, simulating a flow in a channel whith a complex geometry and simulating a flow in a patterned with hydrophilic/hydrophobic surfaces are almost equivalent (section (3.4.1)).We have seen by comparing calculations on the Hele-Shaw model with previous works [21], [22] that it doesn't describes correctly the flux deviation in the channel generated by a ridge-like shaped wall.These arguments lead us to think that chaotic mixing approaches involving a modification on the boundary conditions cannot be described using such an asymptotic reduction.A full 3D Stokes approach is needed for that. 5. Appendix: Existence results.In this appendix, we consider System ( 2.37)-( 2.39), that is we study the case of a non-flat bottom with no-slip boundary conditions in 2D.The aim is to prove the existence of weak and strong solutions to this system.The results are in the spirit of [2].Our aim is to show how one can handle physical boundary data using maximum principle technics and extension of functions in order to prove that this kind of models gives a physically reasonable behavior of the fluid at the outlet.Our proof is not very original and the same strategy can be found in a lot of papers in the literature.The proof of existence with periodic boundary conditions or on the whole space are straightforward.The interest relies here on the boundary conditions that correspond to the physical setting: the fluids are moving because of the pressure gradient.One speaks on pressure-driven flows (see [12]).5.1.Existence of weak solutions.In this section, we deal with weak solutions.The 2D bounded rectangular domain [0, L] × [0, ℓ] is denoted by Ω and its boundary is decomposed into Γ = Γ e ∪ Γ s ∪ Γ L where Γ e is the inlet, Γ s is the outlet and Γ L represents the lateral walls.The outward normals of the domain are denoted by ν.We first recall the boundary conditions associated with System ( 2.37)-( 2.39) on Γ e , we prescribe the velocity u = u 0 (y) and the order parameter ϕ = ϕ 1 (y), on Γ s , we impose P = 0 and ∂ϕ ∂ν = 0, on Γ L , we assume v = 0 and ∂ϕ ∂ν = 0.
We also add the following initial condition on ϕ ϕ(0, x, y) = ϕ 0 (x, y), where ϕ 0 is a given regular function.We also recall that, in this configuration, the coefficients K 2D n and K 3D n read, dropping the index n for simplicity, where η(s) = η 1 s + η 2 (1 − s) for all s ∈ [0, 1].As a consequence, there exists two positive constants K 1 and K 2 such that where K ′ 2D and K ′ 3D denote the derivative of K 2D and K 3D .We now can state the following theorem.
Theorem 5.1.Assume that D(ϕ) is constant, ϕ 1 (y) ∈ H 2 (0, ℓ), ∂ϕ1 ∂ν = 0 on Γ L , u 0 (y) ∈ L 2 (0, ℓ), 0 ≥ 0 and ϕ 0 ∈ L 2 (Ω) such that 0 ≤ ϕ 0 ≤ 1.Then there exists a global weak solution (P, u, v, ϕ) to System ( 2.37)-( 2.39) such that ϕ(0, x, y) = ϕ 0 (x, y) and Proof.The proof is based on a Galerkin method.To this end, we first get rid of the inhomogeneous Dirichlet boundary condition ϕ 1 at the inlet by introducing the following change of unknow φ We now introduce the Galerkin basis associated with the eigenfunctions of the Laplacian operator on Ω with the boundary conditions satisfied by φ, namely We first reorganized the sequence e ij i,j∈N 2 to obtain e n N ∈N and introduce the finite dimensional space V n spanned by e k 0≤k≤n . The approximate solution P n , u n , v n , φ n is given by φ n ∈ V n and x,y P n (x, y) = 0, (5.63) An integration by parts gives, recalling that (5.68) Using ( 5.62), we have (5.69) We then need the following lemma.Proof.We begin with the positivity of P n in the domain Ω.We decompose P n in the following way P n = P + n + P − n where P + n (resp.P − n ) denotes the positive part (resp.the negative part) of P n .Multiply Equation ( 5.63) by P − n , we obtain Using ( 5.62) and Equation ( 5.64), we then derive Since K 2D (ϕ 1 (y)) ≤ 0, K 3D (ϕ 1 (y)) ≤ 0, and u 0 (y) ≥ 0, one can see that It follows that |∇P − n | = 0 and then P − n = 0.As a consequence, P n ≥ 0 in the domain Ω.Furthermore, since on Γ s , one has P n = 0, one deduces that on Γ s ∂P n ∂x ≤ 0.
Recalling that, according to equation ( 5.64), on Γ s , one has the result follows from the fact that K Remark 2. Lemma 5.2 states that for any smooth solution to ( 5.59)-( 5.61), the velocity is directed toward the exterior of the channel at the outlet.
Equipped with Lemma 5.2, we are also able to derive the following maximum principle for φ n .Lemma 5.3.Let φ n be the solution to Equation ( 5.66).Then φ n satisfies 0 ≤ φ n ≤ 1 on Ω.
Proof.We use the notation of Lemma 5.2 and decompose φ n = φ + n +φ − n .We multiply Equation ( 5.66) by φ − n and we integrate over Ω.A straightforward integration by parts gives Using the boundary conditions ( 5.62) on φ and Lemma 5.2, we derive Then ( 5.70) provides We thus obtain φ − n = 0 since at t = 0, φ n ≥ 0, which implies that φ n ≥ 0 on Ω.In order to conclude the proof of Lemma 5.3, it remains to prove that φ ≤ 1 on Ω which is done by performing the same computation with φ = 1 − φ.
We are now able to end the proof of Theorem 5.1.From ( 5.68), we derive, using Lemma 5.2, ( 5.69) and Cauchy-Schwarz inequality We multiply Equations ( 5.74) and ( 5.75) by respectively ∂ x φ and ∂ y φ and integrate over Ω the resulting equations to get Using ( 5.62) and the fact that since φ ∈ V n , one has ∂ 2 x φ = 0 on Γ e , we have (5.79) Plugging ( 5.78) and ( 5.79) into ( 5.76), we get where the constant C depends only on ϕ 1 and u 0 .Again, dealing with ( 5.77), the boundary conditions ( 5.62) provides where the constant C depends only on u 0 and ϕ 1 .It remains to control ∇u and ∇v in L 2 (Ω) and this is done in the following lemma.
Lemma 5.5.There exists a constant C depending only on ϕ 1 , u 0 such that Proof.We proceed in two steps.
• Step 1: Estimates on P .We first recall that multiplying Equation ( 5.59) by P , integrating by parts and using ( 5.58) and ( 5.62), we obtain an H 1 -bound on P which reads as follows where the constant C depends only on u 0 .We apply ∂ y on Equation ( 5.59) to get (5.85) Multiplying Equation ( 5.85) by ∂ y P leads to, after integration over Ω: Note that at the inlet, we don't have any information on ∂ x φ.To overcome this difficulty, we introduce an extension of P in the following sense.Denote where f is a smooth function satisfying f (0) = f (L) = f ′ (L) = f ′′ (L) = 0 and f ′ (0) = 1.Then, there exists a constant C > 0 depending on ϕ 1 , u 0 and f such that (5.93) The equation for Q is The energy estimate on Equation ( 5.94) reads Here, the boundary terms have to be computed separetly and carefully.The boundary conditions on Q are the following ones.On the lateral walls, we have ∂ y P = 0.This implies ∂ x ∂ y P = 0 since ∂ x is a tangential derivative on Γ L .We then need the compatibility condition at the boundary, which is satisfied since on lateral walls ∂ y φ = 0, ∂ y ϕ 1 = 0 and ∂ y u 0 = 0 by hypothesis (see Theorem 5.4).This yields, on Γ L , At the inlet Γ e , one can write using Equation ( 5.60) and ( 5.62), which provides u 0 (y)f (x) = 0, since f (0) = 0 and f ′ (0) = 1.At the outlet, P = 0 implies that ∂ y P = ∂ 2 y P = 0.By Equation ( 5.59), we deduce Since, on Γ s , ∂ x (ϕ 1 + φ) = 0 and K 2D (ϕ 1 + φ) = 0, we have Taking x = L in the last equality and since f (L) = f ′ (L) = f ′′ (L) = 0, one concludes that ∂ 2 x Q = 0 on Γ s .As a consequence, we deduce, using also ( 5.62), that (5.97) (5.99) Collecting ( 5.97), ( 5.98) and ( 5.99), we obtain from ( 5.95)

Figure 1 .
Figure 1.Example of a channel containing two fluids.

Figure 3 .
Figure 3. Placement of the unknowns on the staggered grid.

Figure 4 .
Figure 4. Schematic view of the experimental setup for the study of two fluids interdiffusion.

Figure 5 .
Figure 5. Representation of the different lengths l 1 , l 2 , M 1 involved in the study of the displacement of the interface (thick dashed line).

Figure 6 .
Figure 6.Left: Numerical convergence on the displacement M 1 as a function of time.Right: example of a map on the mixing composition ϕ(x, y) in the channel.

Figure 7 .
Figure 7. Left: displacement M 1 of the interface in the channel as a function of the viscosity contrast λ.Right: relative displacement M 1 l of the interface in the channel as a function of f (η 1 , η 2 ).

Figure 8 .
Figure 8. Schematic view of the experimental setup for the study of the viscous fingering

Figure 11 .
Figure 11.Error between the two models as a function of β

Figure 12 .
Figure 12.Error on the velocities given by the 2.5D Reynolds model ( 2.25)-( 2.27) and the Stokes model ( 1.1), given as a function of the normalized maximum height of the channel.