Skip to main content
  • Research Article
  • Open access
  • Published:

Adjoint-consistent formulations of slip models for coupled electroosmotic flow systems

Abstract

Background

Models based on the Helmholtz `slip' approximation are often used for the simulation of electroosmotic flows. The objectives of this paper are to construct adjoint-consistent formulations of such models, and to develop adjoint-based numerical tools for adaptive mesh refinement and parameter sensitivity analysis.

Methods

We show that the direct formulation of the `slip' model is adjoint inconsistent, and leads to an ill-posed adjoint problem. We propose a modified formulation of the coupled `slip' model, which is shown to be well-posed, and therefore automatically adjoint-consistent.

Results

Numerical examples are presented to illustrate the computation and use of the adjoint solution in two-dimensional microfluidics problems.

Conclusions

An adjoint-consistent formulation for Helmholtz `slip' models of electroosmotic flows has been proposed. This formulation provides adjoint solutions that can be reliably used for mesh refinement and sensitivity analysis.

Background

Introduction

Emerging micro- and nano-electromechanical systems (MEMS/NEMS) have a growing number of applications, ranging from lab-on-a-chip DNA analysis to micro-actuators [1]. By scaling down processes, these systems offer savings in space, cost, and energy for scientific and technological advancement [2]. However, the high manufacturing costs and complex architectures of these systems necessitate the use of numerical simulation tools for optimal design and precise control of their operation [3]. Microfluidic devices operate over various length scales and are best described using multiphysics modeling that involves hydrodynamics, electroosmosis, and chemical species transport models. The development of accurate and efficient computational simulators of these devices is therefore challenging and resource intensive.

Numerical simulations of such complex engineering systems are typically targeted towards the calculation of specific Quantities of Interest (QoI) associated with the systems. Accurate estimation of local QoIs can be achieved using goal-oriented error estimation and adaptive techniques based on the use of adjoint methods [4],[5]. Adjoint methods can also be used to improve the computational performance of parameter sensitivity analyses [6], especially for systems with a large number of parameters. However, the application of adjoint methods to such coupled flow systems is an open area of study.

The objective of the current research work is to apply an adjoint-based Adaptive Finite Element Method (AFEM) to microfluidics problems for mesh refinement and parameter sensitivity analysis. There has been a growing interest in the modeling and numerical simulation of microfluidic systems [7]-[12]. A key issue that one faces while modeling and simulating microfluidic systems is the application of `slip' boundary conditions. Prachittham et al. [13] presented a space-time adaptive finite element method applied to an electroosmotic flow using large aspect ratio elements. However, their work did not consider adjoint techniques and did not use the `slip' boundary coupling condition. On the other hand, van Brummelen et al. [14] and Estep et al. [15] have shown the importance of the treatment of boundary flux coupling for the use of adjoint-based techniques. To our best knowledge, no advances in the application of adjoint-based techniques to microfluidics applications, particularly those involving `slip' boundary coupling, have yet been published in the literature.

In this work, we investigate coupled systems arising in electroosmotic flow (EOF). We show that a naive formulation of the `slip' model leads to an ill-posed adjoint problem and adjoint inconsistency [14]. We provide numerical evidence that the corresponding adjoint solution exhibits spurious oscillations on a simple example dealing with a straight channel flow. Accordingly, we propose a modified variational formulation of the `slip' model, for which the adjoint problem is well-posed and can be computed and used in adaptive mesh refinement and parameter sensitivity analysis.

The paper is organized as follows: Section "Microfluidics modeling" describes the `slip' electroosmotic flow (EOF) model, with a brief discussion of the relevant physics and the applicability of such a model. Then, a variational formulation for a modified version of the `slip' model is presented and its adjoint is derived in Section "Variational formulation of the slip BC EOF model". An analysis of the ill-posed adjoint problem for the naive formulation is presented in Section "Ill-posedness of the adjoint problem for the naive formulation". Next, a variational formulation using penalty boundary conditions (henceforth called the penalty formulation) is proposed in Section "Penalty formulation of the slip BC EOF model". This new formulation is also employed to derive the corresponding adjoint problem and it is shown that the adjoint problem thus obtained is asymptotically consistent with the one derived in Section "Variational formulation of the slip BC EOF model". Numerical experiments are presented in Section "Results and discussion" that support the fact that the adjoint problem obtained using the penalty formulation is well-posed and that the adjoint solution is free of spurious artifacts. The adjoint obtained using the penalty formulation is used for goal-oriented adaptive mesh refinement and sensitivity analyses for a T-channel problem. Section "Conclusions" provides some concluding remarks, followed by a discussion of further work and future applications related to this class of coupled flow models. Finally, the appendix elaborates on the well-posedness of the modified formulation and presents relevant theoretical developments.

Microfluidics modeling

Microfluidics is the branch of fluid mechanics concerned with the understanding, modeling, and control of flows that occur on the micron scale, i.e. where the characteristic length (L) is of the order 10−6 m. Squires and Quake [16], and later Whitesides [17], have presented reviews of the physics and applications of microfluidics. Prominent among them are microflows driven by applied electric fields, through electroosmosis, electrophoresis, or both. In this paper, we consider electroosmotic flow devices, which find wide use in commercial and industrial applications [8],[18]. These devices utilize the properties of the electric double layer (also called the Debye layer) to drive a bulk fluid flow. More detailed descriptions of the electric double layer and electroosmotically driven microfluidic devices can be found in the microfluidics literature [1],[19],[20].

Consider a rectangular open domain Ω· d , with d = 2, and boundary Ω. The boundary Ω is composed of the channel wall Γ w and its inlet/outlet Γ io , such that ∂Ω= Γ w · Γ io . For simplicity, we consider a single species flow through the channel. The one-way coupled, steady-state EOF in a straight rectangular channel can be modeled with the following system of equations,

Δψ = K 2 sinh ( ψ ) in Ω ,
(1a)
Δ · ( σ c φ ) = 0 in Ω ,
(1b)
μΔ u + p = ρ e ( Ψ ) φ in Ω ,
(1c)
· u = 0 in Ω ,
(1d)

where Ψ is the non-dimensional electric potential associated with the double layer, K is a non-dimensional constant, called the Debye-Huckel parameter, φ is the applied potential, σ c is the fluid conductivity, and u and p are the flow velocity and pressure, respectively. The charge density is given by ρ e =−2z v n e sinh(Ψ) and μ is the viscosity of the fluid. The parameters that the charge density depends on are the ion valence z v , electron charge e, and the bulk concentration of ions n . The above equations are supplemented with the boundary conditions,

Ψ = Ψ 0 on Γ w ,
(2a)
n · Ψ = 0 on Γ io ,
(2b)
n · ( σ c Φ ) = 0 on Γ w ,
(2c)
φ = φ io on Γ io ,
(2d)
u = 0 on Γ w ,
(2e)
u · t = 0 on Γ io ,
(2f)
n · ( σ · n ) = 0 on Γ io ,
(2g)

where Ψ0 is the electric zeta potential on the walls Γ w of the channel, φ io is the external potential applied at the inlet and outlet of the channel Γ io , n is the unit outward normal vector to Ω, t is a unit tangential vector along Ω, and σ is the stress tensor:

σ=pI+μ u + ( u ) T .
(3)

The no-slip condition (Eq. (2e)) is applied at the channel walls for the flow velocity. The no-slip conditions, in combination with the body force term in Eq. (1c), which is significant only near the boundary, lead to sharp velocity boundary layer near the wall. This makes the system given by Eq. (1) and Eq. (2) multi-scale. Note that the last two boundary (Eq. (2f) and (2g)) conditions imply that the velocity is normal to Γ io and that the pressure vanishes on Γ io (this is in the case of planar boundaries, see [21]), i.e.

p=0on Γ io .
(4)

Eqs. 1 and (2) define the complete EOF model and constitute a challenging system of coupled multiscale equations. They are computationally expensive, especially for complex geometries, due to the presence of extremely sharp layers near the channel walls because of the electric double layer. Therefore, to reduce the complexity and the computational cost associated with the full model, the Helmholtz-Smoluchowski velocity approximation is introduced into the model. The approximation states that the body-force term in the Stokes equations can be replaced by an effective `slip velocity' on Γ w as,

u wall = є Ψ 0 μ (Φ)=λE,
(5)

where E denotes the applied electric field, and a new parameter λ=є Ψ0/μ has been introduced. Here, є is the permitivity (or dielectric constant) of the fluid medium. A detailed derivation of this approximation can be found in a standard reference [19]. The validity of this approximation has been verified for several examples through both experiments and numerical simulations [10].

The associated slip model of EOF can thus be written as:

· ( σ c φ ) = 0 in Ω ,
(6a)
μΔ u + p = 0 in Ω ,
(6b)
· u = 0 in Ω ,
(6c)

complemented by the boundary conditions,

σ c n φ = 0 on Γ w ,
(7a)
φ = φ io on Γ io ,
(7b)
u + λ φ = 0 on Γ w ,
(7c)
u · t = 0 on Γ io ,
(7d)
n · ( σ · n ) = 0 on Γ io .
(7e)

We see that the slip model spares one from solving the Poisson-Boltzmann equation, whose solution exhibits a thin layer near the wall. As a remark, the slip boundary approximation model given by Eq. (6) and (7) is widely used throughout the microfluidics research and development community for modeling and simulation [1]. The model is even included in the commercial Finite Element software package COMSOL Multiphysics [12].

The slip condition/coupling makes sense only in the direction tangent to the wall. In Eq. (7c), the no-flux boundary condition on the potential (Eq. (7a)) automatically enforces a no penetration boundary condition on the velocity. Thus, expressing the coupling condition as Eq. (5) is a matter of convenience from a notational standpoint. However, as shown in Section "Ill-posedness of the adjoint problem for the naive formulation", using Eq. (5) as such in the formulation of the coupled problem leads to an ill-posed adjoint problem. This ill-posedness will also be illustrated on numerical examples in Section "Results and discussion". We decouple one of the velocity components from the potential as follows,

u · t + λ t φ = 0 u · n + λ n φ = 0 u · t + λ t φ = 0 , u · n = 0 ,
(8)

where we have denoted the normal (φn) and tangential (φt) derivatives of φ by n φ and t φ, respectively. This coupling is equivalent to the one given by Eq. (7c). We now proceed to derive the weak formulation for Eq. (6) using the modified coupling given by Eq. (8).

Methods

Variational formulation of the slip BC EOF model

Variational formulation of primal problem

We derive here the weak formulation of the equations in Eq. (6) and then combine them to give the coupled problem. The potential φ satisfies,

· ( σ c φ ) = 0 in Ω ,
(9a)
φ = φ io on Γ io ,
(9b)
σ c n φ = 0 on Γ w ,
(9c)

We define C+(Ωs) as the space of continuous, strictly positive functions on Ω, and assume that the conductivity σ c ЄC+(Ω). This assumption allows us to simplify the analysis later on and can be easily relaxed. We now introduce the spaces of admissible trial and test functions:

Z = H 1 ( Ω ) ,
(10)
Z 0 = { v Є Z ; v = 0 on Γ io } ,
(11)
Z io = { v Є Z ; v = φ io on Γ io } .
(12)

We can use a lift function φ ~ Є Z io such that φ=φ+ φ ~ with φЄZ0. In this case, the weak form of the problem is given by:

Given σ c Є C + ( Ω ) and φ ~ Є Z io , find φ Є Z 0 such that Ω σ c φ . Ψ dx = Ω σ c φ ~ . ψ dx , ψ Є Z 0 .
(13)

We now consider the non-dimensionalized stationary Stokes equations (with μ=1) with slip boundaries,

Δ u + p = 0 in Ω ,
(14a)
. u = 0 in Ω ,
(14b)
u . t + λ t φ = 0 on Γ w ,
(14c)
u . n = 0 on Γ w ,
(14d)
u . t = 0 on Γ io ,
(14e)
n . ( σ . n ) = 0 on Γ io .
(14f)

We look for the velocity and pressure fields in the function spaces,

X = H 1 ( Ω ) 2 ,
(15a)
X 0 = v Є X ; v = 0 on Γ w , v . t = 0 on Γ io ,
(15b)
M = p Є L 2 ( Ω ) ; Ω p dx = 0 .
(15c)

Eq. (14c) requires that Φ lie in the space H2(Ω). Note that if Ω is Lipschitz and convex, or C1-continuous, then the elliptic regularity theorem guarantees that ΦH2(Ω). However, to derive a well-posed adjoint problem for the slip EOF model, we need to enforce the coupling boundary conditions in a weak sense. In fact, if we enforce the condition in a weak sense, the curl theorem shows that it is sufficient to have ΦZ. Consider the function space H 0 1 2 (∂Ω),

H 0 1 2 (∂Ω)={vЄ H 1 2 (∂Ω);v=0on Γ io }.
(16)

We have the following proposition.

Proposition 1.

Let Ω be Lipschitz. Further, let ΦZ, and w,g H 0 1 2 (∂Ω). Define b(g,w)= g , w H 1 2 , H 1 2 and f(w)= λ t Φ , w Find g Є H 0 1 2 ( ∂Ω ) , such that b ( g , w ) = f ( w ) w Є H 1 2 ( ∂Ω ) ,

(17)

is well-posed, and

g H 1 2 ( ∂Ω ) c(Ω)|λ|Φ H 1 ( Ω ) ,
(18)

where c(Ω) is a positive constant that depends only on the domain Ω.

Proof.

Since the bilinear form b(g,w) is defined using the H 1 2 (∂Ω) inner product, it is bounded and coercive. We can use the curl theorem and a special extension of w to show the boundedness of f(w),

f(w)= λ t Φ , w H 1 2 , H 1 2 = ∂Ω λ t Φwds.

Consider an H1(Ω) extension of the function w, denoted by w ~ . From corollary B. 53, on page 488 in [22], we know that for each wЄ H 1 2 (∂Ω), there exists an extension w*ЄH1(Ω), such that,

w * H 1 ( Ω ) c(Ω)w H 1 2 ( ∂Ω ) .
(19)

From Proposition 12 in (F.J. Sayas: Weak normal derivatives, normal and tangential traces, and tangential differential operators on Lipschitz boundaries, unpublished), we have,

f(w)= λ t Φ , w H 1 2 , H 1 2 =λ Ω Φ . × w * k × Φ . w * k dx.

Using the fact that the gradient of the curl is the zero vector, and the vector cross product identity, we have,

f(w)=λ Ω Φ × w * .k.
(20)

We then easily have the following bound on |f(w)|,

| f ( w ) | | λ | Φ H 1 ( Ω ) w * H 1 ( Ω ) , | λ | Φ H 1 ( Ω ) c ( Ω ) w H 1 2 ( ∂Ω ) .
(21)

Therefore, f H Ω 1 2 c(Ω)|λ|Φ H 1 ( Ω ) , and an application of the Lax-Milgram theorem gives the well-posedness of the variational problem, and the bound g H 1 2 ( ∂Ω ) c(Ω)|λ|Φ H 1 ( Ω ) . □

We can now define the operator :ZX,

( Φ ) = ( g , 0 ) on ∂Ω , where g satisfies Eq. (17) ,
(22a)
( Φ ) X c g H 1 2 ( ∂Ω ) .
(22b)

The bounded extension corollary B. 53, on page 488 in [22] implies that there is at least one possible lift for any given Φ. For the purposes of deriving the variational formulation, we can pick any admissible lift. It is trivial to show that the operator is linear. Proposition 1 further gives,

(Φ) X cg H 1 2 ( ∂Ω ) c(Ω)|λ|Φ H 1 ( Ω ) .
(23)

Therefore, is a bounded operator. Using this lift operator, we can write u=w+(Φ)=w+(φ)+( Φ ~ ), where wX0, and write a weak form for Eq. (14),

Given Φ Z , find ( w , p ) X 0 · M such that , Ω w . v p . v q . w dx = Ω ( φ ) . v q . ( φ ) dx ( v , q ) X 0 · M.
(24)

Combining Eq. (13) and Eq. (24) together, we get the coupled variational statement:

Given σ c C + ( Ω ) and φ ~ Z b , find ( φ , w , p ) Z 0 × X 0 × M such that , Ω σ c φ · ψ dx + Ω [ w + ( φ ) ] · v p · v q · [ w + ( φ ) ] dx = Ω σ c Φ ~ · ψ dx Ω ( Φ ~ ) · v q · ( φ ~ ) dx ( ψ , v , q ) Z 0 × X 0 × M ,
(25)

where we emphasize that the integrals involving the lift velocity, which depend on the solution φ, will be part of the bilinear form for this coupled problem. We can recast the bilinear form above in more compact notation as,

Given σ c C + ( Ω ) and φ ~ Z b , find U Z 0 · X 0 × M such that , A ( U , V ) = F ( V ) V Z 0 × X 0 × M ,
(26)

where U = (φ,w,p) and V = (ψ,v,q).

We have thus incorporated the coupling condition within our bilinear form and can prove that the coupled problem (26) is well-posed. We refer the interested reader to Appendix A and proceed with the derivation of the corresponding adjoint problem.

Adjoint problem

Given the primal weak form Eq. (26), we have the corresponding weak form for the adjoint problem associated with the quantity of interest Q: Z 0 × X 0 ×M

Given σ c C + ( Ω ) , find U * Z 0 × X 0 × M such that , A ( V , U * ) = Q ( V ) V Z 0 × X 0 × M ,
(27)

where U = (φ*,w*,p) is the adjoint solution and Q(U) is a linear functional that corresponds to the quantity of interest (QoI). The full weak form for the adjoint problem reads,

Ω σ c ψ · φ * dx + Ω ( v · w * q · w * p * · v ) dx + Ω ( ( ψ ) · w * · ( ψ ) p * ) dx = Q ( V ) V = ( ψ , v , q ) Z 0 × X 0 × M.
(28)

Note that the adjoint Stokes problem solution w* and p* are coupled to the adjoint potential φ* through the lift operator (.) acting on the test function ψ. The adjoint problem (28) is also well-posed, see Corollary 2 in Appendix A.

We recall that applying the forward coupling condition Eq. (14c) in a strong sense required that the forward potential φH2(Ω). Correspondingly, to obtain a strong form corresponding to Eq (28), we require that the adjoint velocity w*[H2(Ω)]2 and the adjoint pressure p*H1(Ω). We also have the adjoint stress tensor,

σ * = p * I+( u * + ( u * ) T ).
(29)

Now integrating by parts we obtain,

Ω · ( σ c φ ) Ψ dx + ∂Ω σ c n φ Ψ ds + Ω · σ · v + ( Ψ ) + ∂Ω v + ( Ψ ) · ( σ · n ) ds + Ω q · w dx = Q ( V ) .
(30)

The boundary term in the formulation of the adjoint problem becomes

∂Ω v + ( Ψ ) · ( σ · n ) ds = Γ w v + ( Ψ ) · ( σ · n ) ds + Γ io v + ( Ψ ) · ( σ · n ) ds = Γ w n · v + ( Ψ ) # 0 n · ( σ · n ) + t · v + ( Ψ ) t · ( σ · n ) # λ t Ψ t · ( σ · n ) ds + Γ io n · v + ( Ψ ) n · ( σ · n ) + t · v + ( Ψ ) # 0 t · ( σ · n ) ds = Γ w Γ w · λ t · ( σ · n ) t Ψ ds + Γ io n · v n · ( σ · n ) ds ,
(31)

where we have used integration by parts for the tangential derivative term along Γ w [23],

Γ w (λ t Ψ)(t·z)ds= Γ w Γ w · λ t · z t Ψds,
(32)

where Γ w ·v denotes the surface divergence of the vector v. Replacing these terms in the adjoint formulation and setting,

Q(U)= Ω k(x)u·αdx+ Γ io k s (s)u·nds.
(33)

We can write Eq. (31) as,

Ω Ψ · ( σ c φ ) dx + Γ w Ψ n · ( σ c φ ) + Γ w · λ t · ( σ · n ) t ds Ω v + ( Ψ ) · · σ + k ( x ) α + Γ io n · v + ( Ψ ) n · ( σ · n k s ( s ) ) ds Ω q · w dx = 0 ( Ψ , v , q ) Z 0 · X 0 · M.
(34)

The strong form of the adjoint system then reads:

· ( σ c φ ) = 0 in Ω ,
(35a)
w + p = k α in Ω ,
(35b)
w = 0 in Ω ,
(35c)

with three boundary conditions on Γ io :

φ = 0 on Γ io ,
(36a)
w · t = 0 on Γ io ,
(36b)
n · ( σ · n ) = k s on Γ io ,
(36c)

and three boundary conditions on Γ w :

σ c n φ + Γ w · λ t · ( σ · n ) t = 0 on Γ w ,
(37a)
w = 0 on Γ w .
(37b)

We readily observe that the adjoint Stokes problem can be solved first, independently of the adjoint potential problem, but that the latter does depend on the former through the Neumann coupling condition Eq. (37a). We also note that this coupling condition involves the tangential derivatives of the adjoint stress tensor on the boundary.

Ill-posedness of the adjoint problem for the naive formulation

If both the normal and tangential coupling components are considered in the formulation of the primal problem, the corresponding adjoint boundary integrals (see Eq. (31)) would read,

Ω∂ v + ( Ψ ) · ( σ · n ) ds = Γ w v + ( Ψ ) · ( σ · n ) ds + Γ io v + ( Ψ ) · ( σ · n ) ds = Γ w n · v + ( Ψ ) n · ( σ · n ) + t · v + ( Ψ ) t · ( σ · n ) # λ t Ψ t · ( σ · n ) ds + Γ io n · v + ( Ψ ) n · ( σ · n ) + t · v + ( Ψ ) # 0 t · ( σ · n ) ds = Γ w λ n Ψ n · ( σ · n ) + Γ w Γ w · λ t · ( σ · n ) t Ψ ds + Γ io n · v n · ( σ · n ) ds.

The boundary conditions on Γ w for the strong form of the adjoint system would then read:

n · ( σ c φ ) + Γ w · λ t · ( σ · n ) t = 0 on Γ w ,
(38a)
n · ( σ · n ) = 0 on Γ w ,
(38b)
w = 0 on Γ w .
(38c)

Thus the adjoint problem for the ill-posed primal problem contains four boundary conditions, despite there being only three variables. Therefore, the ill-posed primal problem leads to an adjoint inconsistent formulation, specifically in the boundary terms. Such an adjoint inconsistency can lead to oscillations in the discrete solutions of the adjoint problem [24].

Penalty formulation of the slip BC EOF model

Penalty formulation of the primal problem

The imposition of the boundary conditions for the primal and adjoint problems derived in the previous section can be extremely challenging, mainly due to the regularity requirements for the corresponding spaces in the interior and the difficulty of constructing appropriate Finite Element spaces. Therefore, we seek to impose the coupling constraint using a penalty method and relax the regularity requirements on the spaces containing the primal solution and the adjoint. The penalty method is an easy and robust approach for applying Dirichlet boundary conditions. See [25] and [26] for analysis of the method. As we shall see, the penalty method is a natural method for weak enforcement of the coupling. In addition, the penalty formulation gives us an adjoint that is asymptotically consistent with the one obtained using the lift technique. The variational formulation Eq. (25) with equivalent penalty boundary conditions can be given as:

Given σ c C + ( Ω ) and Φ io H 1 2 ( Γ io ) , find ( Φ , u , p ) Z · X · M such that , Ω σ c Φ · Ψ dx + 1 Γ io Φ Ψ ds + Ω u · v p · v q · u dx + 1 Γ w ( u · n ) ( v · n ) ds + 1 Γ w ( u · t ) ( v · t ) ds + 1 Γ io ( u · t ) ( v · t ) ds + 1 Γ w λ ( t Φ ) ( v · t ) ds = 1 Γ io Φ io Ψ ds ( Ψ , v , q ) Z · X · M.
(39)

As was done in Eq. (26), we can recast the bilinear form above in more compact notation as,

Given σ c C + ( Ω ) and Φ io H 1 2 ( Γ io ) , find U Z · X · M such that , A ( U , V ) = F ( V ) V Z · X · M ,
(40)

where U = (Φ ,u ,p ) and V = (Ψ,v,q). We now verify that the weak form in Eq. (39) is indeed consistent and converges to the BVP Eq. (9) and Eq. (14) in the limit as the penalty parameter tends to zero. Integrating Eq. (39) by parts, we obtain,

Ω · ( σ c Φ ) Ψ dx + Γ w σ c n Φ Ψ ds + Γ io σ c n Φ + 1 ( Φ Φ io ) Ψ ds + Ω u + p · v q · u ds + Γ w n · ( σ · n ) + 1 ( u · n ) ( v · n ) ds + Γ w t · ( σ · n ) + 1 ( λ t Φ ) + u · t ) ( v · t ) ds + Γ io ( n · ( σ · n ) ) ( v · n ) ds + Γ io t · ( σ · n ) + 1 ( u · t ) v · t ds = 0 ( Ψ , v , q ) Z × X × M ,
(41)

where σ is the stress tensor corresponding to the penalized solution,

σ = p I+( u + ( u ) T ).
(42)

The equivalent strong form for finite non-zero is,

· ( σ c Φ ) = 0 in Ω ,
(43a)
u + p = 0 in Ω ,
(43b)
· u = 0 in Ω ,
(43c)

with the boundary conditions on Γ io :

σ c n Φ + 1 ( Φ Φ io ) = 0 on Γ io ,
(44a)
n · ( σ · n ) = 0 on Γ io ,
(44b)
t · ( σ · n ) + 1 ( u t ) = 0 on Γ io ,
(44c)

and the boundary conditions on Γ w :

σ c n Φ Γ = 0 on σ w ,
(45a)
n · ( σ · n ) + 1 ( u · n ) = 0 on Γ w ,
(45b)
t · ( σ · n ) + 1 ( λ t Φ + u · t ) = 0 on Γ w .
(45c)

One can observe that the penalty method replaces the Dirichlet boundary conditions with a Robin condition. However, upon taking the limit →0 one formally recovers the original problems Eq. (9) and Eq. (14).

Adjoint problem associated with the penalty formulation

We now derive the adjoint problem corresponding to the weak form Eq. (39),

Find ( Φ , u , p ) Z × X × M such that , Ω σ c Φ · Ψ dx + 1 Γ io Φ Ψ ds + 1 Γ w λ t Ψ ( u · t ) ds + 1 Γ w ( u · n ) ( v · n ) ds + 1 Γ w ( u · t ) ( v · t ) ds + 1 Γ io ( u · t ) ( v · t ) ds + Ω ( u · v p · v q · u ) dx = Ω k ( x ) v · α dx + Ω∂ k s ( s ) v · n ds ( Ψ , v , q ) Z × X × M.
(46)

Using integration by parts for the term involving the tangential derivative along Γ w [23], i.e.

Γ w λ t Ψ ( u ·t)ds= Γ w Γ w · ( λ u · t ) t Ψ ds,
(47)

and upon integrating by parts the higher-order terms and combining integrals with same test functions, one obtains:

Ω · ( σ c Φ ) Ψ dx + Γ w σ c n Φ 1 Γ w · ( λ u · t ) t Ψ ds + Γ io σ c n Φ + 1 Φ Ψ ds + Ω u + p k ( x ) α · v dx + Ω · u q dx + Γ w n · ( σ · n ) + 1 ( u · n ) ( v · n ) ds + Γ w t · ( σ · n ) + 1 ( u · t ) ( v · t ) ds + Γ io n · ( σ · n ) k s ( s ) ( v · n ) ds + Γ io t · ( σ · n ) + 1 ( u · t ) ( v · t ) ds = 0 ( Ψ , v , q ) Z × X × M ,
(48)

where σ is the penalized adjoint stress tensor. The equivalent strong form for finite non-zero is,

u + p = k in Ω ,
(49a)
· u Ω = 0 in ,
(49b)
· ( σ c Φ ) = 0 in Ω ,
(49c)

with the three boundary conditions on Γ io :

n · ( σ σ · n ) = k s ( s ) on Γ io ,
(50a)
t · ( σ · n ) + 1 ( u · t ) = 0 on Γ io ,
(50b)
σ c n Φ + 1 Φ = 0 on Γ io ,
(50c)

and the three boundary conditions on Γ w :

n · ( σ · n ) + 1 ( u · n ) = 0 on Γ w ,
(51a)
t · ( σ · n ) + 1 ( u · t ) = 0 on Γ w ,
(51b)
σ c n Φ 1 Γ w · ( λ u · t ) t = 0 on Γ w .
(51c)

In the next section, we show that above problem is asymptotically consistent with the previous formulation of the adjoint problem, in the sense that we recover the adjoint corresponding to the strong problem, i.e. Eq. (35), Eq. (36), and Eq. (37) as tends to zero.

Consistency of the adjoint penalty problem

The main issue is to ensure that the adjoint solution u to the adjoint problem obtained from the penalized formulation does in fact converge to the adjoint solution u obtained from the primal formulation as the penalty parameter tends to zero, as illustrated in Figure 1. In this case, one has to show that the resulting boundary conditions associated with the penalized and non-penalized formulations of the adjoint problems are consistent.

Figure 1
figure 1

Consistency of the adjoint problems associated with the original and penalty formulations. The question here is whether the adjoint problem obtained from the penalty formulation converges to the adjoint problem derived from the original formulation in the limit when the penalty parameter tends to zero.

Recall that the non-penalized adjoint solution (Φ,u) for the problem of interest satisfies the following boundary conditions

Φ = 0 on Γ io ,
(52a)
n · ( σ · n ) = k on Γ io ,
(52b)
u · t = 0 on Γ io ,
(52c)
u = 0 on Γ w ,
(52d)
σ c n Φ + Γ w · λ t · ( σ · n ) t = 0 on Γ w ,
(52e)

while the penalized adjoint solution ( Φ , u ) satisfies,

σ c n Φ 1 Φ = 0 on Γ io ,
(53a)
n · ( σ · n ) = k on Γ io ,
(53b)
t · ( σ · n ) + 1 ( u · t ) = 0 on Γ io ,
(53c)
t · ( σ · n ) + 1 ( u · t ) = 0 on Γ w ,
(53d)
n · ( σ · n ) + 1 ( u · n ) = 0 on Γ w ,
(53e)
σ c n Φ 1 Γ w · ( λ u · t ) t = 0 on Γ w .
(53f)

To formally interpret Eq. (53f), we can substitute Eq. (53d) into Eq. (53f) as follows,

t · ( σ · n ) + 1 ( u · t ) = 0 λ ( u · t ) = λ t · ( σ · n ) , σ c n Φ = Γ w · λ u · t t σ c n Φ = Γ w · λ t · ( σ · n ) t .
(54)

We can thus derive the following boundary conditions for the adjoint potential

σ c n Φ + Γ w · λ t · ( σ · n ) t = 0 on Γ w ,
(55a)
Φ + σ c n Φ = 0 on Γ io ,
(55b)

which are the same as those for the non-penalized adjoint in the limit →0. Eq. 53d corresponds to a penalty representation of the tangential boundary flux. Further discussion of this representation is presented in [27]. We thus see that the penalized formulation of the electroosmotic flow problem is adjoint consistent in the limit as →0, and the use of the discrete penalized adjoint solution ( ( Φ ) h , ( u ) h , ( p ) h ) in adjoint-based error estimation and sensitivity analysis is justified. Thus, if the forward problem is computed numerically using a discrete representation of Eq. (39), the adjoint solution can also be easily computed by taking the transpose of the stiffness matrix associated with the forward problem.

References

  1. Karniadakis G, Beskok A, Aluru NR: Microflows and Nanoflows: Fundamentals and Simulation. Springer, New York; 2005.

    Google Scholar 

  2. Feynman RP: There's plenty of room at the bottom. Eng Sci 1960,23(5):22–36.

    Google Scholar 

  3. Knio OM, Ghanem RG, Matta A, Najm HN, Debusschere B, Le Maiǐtre OP: Quantitative Uncertainty Assessment and Numerical Simulation of Micro-Fluid Systems. Technical report, Johns Hopkins University, Baltimore, MD; 2005.

    Google Scholar 

  4. Becker R, Rannacher R: An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica 2003, 10: 1–102.

    MathSciNet  Google Scholar 

  5. Estep D, Carey V, Ginting V, Tavener S, Wildey T: A posteriori error analysis of multiscale operator decomposition methods for multiphysics models. J Phys: Conf Ser 2008, 125: 012075. IOP Publishing IOP Publishing

    Google Scholar 

  6. Ionescu-Bujor M, Cacuci DG: A comparative review of sensitivity and uncertainty analysis of large-scale systems. i: Deterministic methods. Nuclear Sci Eng 2004,147(3):189–203. 10.13182/NSE03-105CR

    Article  Google Scholar 

  7. Zhang Y, Wong TN, Yang C, Ooi KT: Electroosmotic flow in irregular shape microchannels. Int J Eng Sci 2005,43(19–20):1450–1463. 10.1016/j.ijengsci.2005.05.017

    Article  Google Scholar 

  8. Ren L, Sinton D, Li D: Numerical simulation of microfluidic injection processes in crossing microchannels. J Micromech Microeng 2003, 13: 739. 10.1088/0960-1317/13/5/329

    Article  Google Scholar 

  9. MacInnes JM, Du X, Allen RWK: Prediction of electrokinetic and pressure flow in a microchannel T-junction. Phys Fluids 2003, 15: 1992. 10.1063/1.1580479

    Article  Google Scholar 

  10. Hahm J, Balasubramanian A, Beskok A: Flow and species transport control in grooved microchannels using local electrokinetic forces. Phys Fluids 2007, 19: 013601. 10.1063/1.2432893

    Article  Google Scholar 

  11. Craven TJ, Rees JM, Zimmerman WB: On slip velocity boundary conditions for electroosmotic flow near sharp corners. Phys Fluids 2008, 20: 043603. 10.1063/1.2906344

    Article  Google Scholar 

  12. Zimmerman WB, Rees JM, Craven TJ: Rheometry of non-newtonian electrokinetic flow in a microchannel T-junction. Microfluidics Nanofluidics 2006,2(6):481–492. 10.1007/s10404-006-0089-4

    Article  Google Scholar 

  13. Prachittham V, Picasso M, Gijs MAM: Adaptive finite elements with large aspect ratio for mass transport in electroosmosis and pressure-driven microflows. Int J Numerical Methods Fluids 2010,63(9):1005–1030.

    MathSciNet  Google Scholar 

  14. van Brummelen EH, van der Zee KG, Garg VV, Prudhomme S: Flux evaluation in primal and dual boundary-coupled problems. J Appl Mech 2012,79(1):010904. American Society of Mechanical Engineers American Society of Mechanical Engineers 10.1115/1.4005187

    Article  Google Scholar 

  15. Estep D, Tavener S, Wildey T: A posteriori error estimation and adaptive mesh refinement for a multi-discretization operator decomposition approach to fluid-solid heat transfer. J Comput Phys 2010, 229: 4143–4158. 10.1016/j.jcp.2010.02.003

    Article  MathSciNet  Google Scholar 

  16. Squires TM, Quake SR: Microfluidics: Fluid physics at the nanoliter scale. Rev Modern Phys 2005,77(3):977–1026. 10.1103/RevModPhys.77.977

    Article  Google Scholar 

  17. Whitesides GM: The origins and the future of microfluidics. Nature 2006,442(7101):368–373. 10.1038/nature05058

    Article  Google Scholar 

  18. Chen C, Lin H, Lele S, Santiago J: Convective and absolute electrokinetic instability with conductivity gradients. J Fluid Mech 2005, 524: 263–303. 10.1017/S0022112004002381

    Article  Google Scholar 

  19. Masliyah JH, Bhattacharjee S: Electrokinetic and Colloid Transport Phenomena. Wiley-Interscience, New York; 2006.

    Book  Google Scholar 

  20. Barz D: Comprehensive model of electrokinetic flow and migration in microchannels with conductivity gradients. Microfluidics and Nanofluidics 2009,7(2):249–265. Springer Springer 10.1007/s10404-008-0382-5

    Article  Google Scholar 

  21. Barth WL, Carey GF: On a boundary condition for pressure-driven laminar flow of incompressible fluids. Int J Numerical Methods Fluids 2007,54(11):1313–1325. 10.1002/fld.1427

    Article  MathSciNet  Google Scholar 

  22. Ern A, Guermond JL: Theory and Practice of Finite Elements, volume 159 of Applied Mathematical Sciences. Springer, New York; 2004.

    Book  Google Scholar 

  23. Delfour M, Zolésio J-P (2001) Shapes and Geometries: Analysis, Differential Calculus, and Optimization. SIAM, New York Vol. 4 of SIAM Series on Advances in Design and Control. Delfour M, Zolésio J-P (2001) Shapes and Geometries: Analysis, Differential Calculus, and Optimization. SIAM, New York Vol. 4 of SIAM Series on Advances in Design and Control.

  24. Hartmann R: Adjoint consistency analysis of discontinuous galerkin discretizations. SIAM J Numer Anal 2007,45(6):2671–2696. 10.1137/060665117

    Article  MathSciNet  Google Scholar 

  25. Babuška I: The finite element method with penalty. Math Comp 1973,27(122):221–228. 10.1090/S0025-5718-1973-0351118-5

    Article  MathSciNet  Google Scholar 

  26. Utku M, Carey GF: Boundary penalty techniques. Comput Methods Appl Mech Eng 1982,30(1):103–118. 10.1016/0045-7825(82)90057-3

    Article  MathSciNet  Google Scholar 

  27. Garg VV (2012) Coupled flow systems, adjoint techniques and uncertainty quantification. PhD Thesis, The University of Texas at Austin. Garg VV (2012) Coupled flow systems, adjoint techniques and uncertainty quantification. PhD Thesis, The University of Texas at Austin.

  28. Kirk B, Peterson J, Stogner R, Carey G: libMesh: a C++ library for parallel adaptive mesh refinement/coarsening simulations. Eng Comput 2006,22(3):237–254. 10.1007/s00366-006-0049-3

    Article  Google Scholar 

  29. Demkowicz L: Computing with Hp-adaptive Finite Elements: One and Two Dimensional Elliptic and Maxwell Problems. CRC Press, New York; 2006.

    Book  Google Scholar 

  30. Van der Zee K (2009) Goal-adaptive discretization of fluid-structure interaction Dissertation. Delft Institute of Technology, ISBN: 9789079488544 TU Delft institutional repository. Van der Zee K (2009) Goal-adaptive discretization of fluid-structure interaction Dissertation. Delft Institute of Technology, ISBN: 9789079488544 TU Delft institutional repository.

Download references

Acknowledgements

Vikram Garg is grateful for the support of the Bruton fellowship and the University Continuing fellowship from The University of Texas at Austin. Kris van der Zee is grateful for the support of this work by the 2010 NWO Innovational Research Incentives Scheme (IRIS) Grant 639.031.033. Serge Prudhomme is sponsored by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada. Serge Prudhomme is also a participant of the KAUST SRI center for Uncertainty Quantification in Computational Science and Engineering.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Vikram V Garg.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

The author's have presented an adjoint-consistency analysis of the widely used Helmholtz `slip' electroosmotic flow models. It was shown that a direct formulation of such models leads to an ill-posed adjoint problem. A modified formulation is proposed and it is shown that this formulation is well-posed and adjoint-consistent. Numerical experiments show that the adjoint obtained from the modified formulation can be used for goal-oriented adaptive mesh refinement and adjoint-based sensitivity analysis. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0), which permits use, duplication, adaptation, distribution, and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Garg, V.V., Prudhomme, S., van der Zee, K.G. et al. Adjoint-consistent formulations of slip models for coupled electroosmotic flow systems. Adv. Model. and Simul. in Eng. Sci. 1, 15 (2014). https://doi.org/10.1186/s40323-014-0015-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40323-014-0015-3

Keywords