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

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.


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 http://www.amses-journal.com/content/2/1/15 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][8][9][10][11][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 wellposed 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 http://www.amses-journal.com/content/2/1/15 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 ⊂ R 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, 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, 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: 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 http://www.amses-journal.com/content/2/1/15 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.
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, 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: −μ u + ∇p = 0 in , complemented by the boundary conditions, 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 http://www.amses-journal.com/content/2/1/15 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, 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).

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, We define C + ( ) 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: We can use a lift functionφ ∈ Z io such that φ = ϕ +φ with ϕ ∈ Z 0 . In this case, the weak form of the problem is given by: We now consider the non-dimensionalized stationary Stokes equations (with μ = 1) with slip boundaries, n · (σ · n) = 0 on io . (14f) http://www.amses-journal.com/content/2/1/15 We look for the velocity and pressure fields in the function spaces, Eq. (14c) requires that φ lie in the space H 2 ( ). Note that if ∂ is Lipschitz and convex, or C 1 -continuous, then the elliptic regularity theorem guarantees that φ ∈ H 2 ( ). 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 We have the following proposition.
. Then, the variational problem, is well-posed, and 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), Consider an H 1 ( ) extension of the function w, denoted byw. From corollary B. 53, on page 488 in [22], we know that for each w ∈ H 1 2 (∂ ), there exists an extension w * ∈ H 1 ( ), such that, From Proposition 12 in (F.J. Sayas: Weak normal derivatives, normal and tangential traces, and tangential differential operators on Lipschitz boundaries, unpublished), we have, Using the fact that the gradient of the curl is the zero vector, and the vector cross product identity, we have, We then easily have the following bound on |f (w)|, Therefore, f and an application of the Lax-Milgram theorem gives the well-posedness of the variational problem, and the bound g We can now define the operator : Z → X, 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, Therefore, is a bounded operator. Using this lift operator, we can write u = w+ (φ) = w + (ϕ) + (φ), where w ∈ X 0 , and write a weak form for Eq. (14), Combining Eq. (13) and Eq. (24) together, we get the coupled variational statement: 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, 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 : 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, 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 φ ∈ H 2 ( ). Correspondingly, to obtain a strong form corresponding to Eq (28), we require that the adjoint velocity w * ∈ H 2 ( ) 2 and the adjoint pressure p * ∈ H 1 ( ). We also have the adjoint stress tensor, Now integrating by parts we obtain, The boundary term in the formulation of the adjoint problem becomes where we have used integration by parts for the tangential derivative term along w [23], where ∇ w · v denotes the surface divergence of the vector v. Replacing these terms in the adjoint formulation and setting, We can write Eq. (31) as, The strong form of the adjoint system then reads: with three boundary conditions on io : and three boundary conditions on w : (37b) http://www.amses-journal.com/content/2/1/15 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, The boundary conditions on w for the strong form of the adjoint system would then read: 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 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 http://www.amses-journal.com/content/2/1/15 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: As was done in Eq. (26), we can recast the bilinear form above in more compact notation as, 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, where σ is the stress tensor corresponding to the penalized solution, The equivalent strong form for finite non-zero is, with the boundary conditions on io : and the boundary conditions on w : 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), Using integration by parts for the term involving the tangential derivative along w [23], i.e. and upon integrating by parts the higher-order terms and combining integrals with same test functions, one obtains: where σ * is the penalized adjoint stress tensor. The equivalent strong form for finite nonzero is, with the three boundary conditions on io : and the three boundary conditions on w : n · (σ * · n) + 1 (u * · n) = 0 on w , (51a) 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. http://www.amses-journal.com/content/2/1/15

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 while the penalized adjoint solution (φ * , u * ) satisfies, t · (σ * · n) + 1 (u * · t) = 0 on w , (53d) n · (σ * · n) + 1 (u * · n) = 0 on w , (53e) To formally interpret Eq. (53f), we can substitute Eq. (53d) into Eq. (53f) as follows, We can thus derive the following boundary conditions for the adjoint potential 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 http://www.amses-journal.com/content/2/1/15 be easily computed by taking the transpose of the stiffness matrix associated with the forward problem.

Results and discussion
We now consider the application of the new EOF formulation on specific microfluidic examples. First, we simulate a flow in a straight microchannel driven purely by electroosmosis. The objective here is to highlight the convergence and stability properties of the adjoint solution. We then showcase an adjoint-based adaptive strategy for mesh refinement on a T-shaped microchannel flow and adjoint-based parameter sensitivity analyses. We discuss the improvement of the convergence rates with respect to quantities of interest and their sensitivities when using adjoint-based techniques. Simulations are performed using the adjoint capabilities added to the libMesh Finite Element library [28]. For both applications, second-order Lagrange elements are employed for the potential and velocity approximations. Linear Lagrange elements are selected to approximate the pressure field in order to satisfy the inf-sup condition. Initial meshes in all the experiments dealing with the straight and T-channel domains consist of structured meshes of bi-quadratic quadrilateral elements. Numerical errors to generate the convergence plots are estimated in this work using so-called overkilled reference solutions of the two problems. These are obtained on a uniform mesh of 428,676 degrees of freedom for the straight channel problem and a combined adaptive-uniform mesh with 288,160 degrees of freedom for the T-channel problems. Numerical solutions are calculated using an ILU preconditioned GMRES iterative method for both problems. The linear algebra library PETSc is accessed through libMesh to obtain these solutions. The penalty parameter was set to the constant value of 10 −8 for all the numerical experiments.

Electroosmotic flow in a straight channel
Numerical experiments are performed here in the case of an electroosmotic flow in a straight channel. The channel has unit width and the length is five times the width. Since the objective of these simulations is to illustrate the numerical properties of the adjoint solution obtained by using Eq. (39), we set arbitrary values of the model parameters rather than choosing values representative of an actual flow. The fluid viscosity μ, electroosmotic slip parameter κ, and fluid density ρ are all taken to be unity. Constant potentials φ i = 8 and φ o = 0 are prescribed at the inlet in and outlet out boundaries, respectively.
The electric conductivity of the fluid is chosen as σ c = 1 + x (note that this particular form of the conductivity is chosen for no other reason than better illustrate the properties of the computed adjoint). The quantity of interest is defined here in terms of the bounded linear functional: where α = (1, 1). Such a bounded functional ensures that any oscillations observed in the numerical results solely arise from the definition of the bilinear form in the adjoint problem. We consider the formulation of the adjoint problem as given in Eq. (46). After computing the forward solution using the numerical set-up as above, we obtain the adjoint potential φ * as seen in Figure 2. By observing the norm of the computed adjoint http://www.amses-journal.com/content/2/1/15 solutions with successive mesh refinements, it was numerically verified that the adjoint potential and velocities were all in H 1 ( ). We also studied the convergence rates for the approximate primal and adjoint potential and x-component of the velocity. Recall that the potential and velocity fields are both approximated using second-order Lagrange elements so that one would expect first-order convergence rates with respect to the number of degrees of freedom in the H 1 norm. However, one can observe in Figure 3 that the primal velocity and the adjoint potential converge at a slower than optimal rate while the primal potential and adjoint velocity converge at the optimal rate. This is due to the tangential 'slip' coupling given by Eq. (8) for the forward problem and the Neumann conditions Eq. (55a) for the adjoint problem. Essentially, we can say that the forward Stokes problem and the adjoint potential problem have non-accurate data, leading to higher errors in the computation of their solutions.

Consequences of coupling both normal and tangential components As we saw in
Section "Ill-posedness of the adjoint problem for the naive formulation", coupling both the normal and tangential components of the velocity to the potential leads to an ill-posed adjoint problem. This ill-posedness is inherited by the penalized formulation as well. On directly enforcing the constraint Eq. (7c) on the wall boundary rather than splitting the two velocity components as in Eq. (8), one observes spurious oscillations in the numerical adjoint potential field φ * , as shown in Figure 4. One clearly observes in Figure 4(a) the http://www.amses-journal.com/content/2/1/15 presence of closed contour lines along the top and bottom wall boundaries. This result is confirmed in Figure 4(b), which shows the solution φ * along the top boundary. For uniform meshes, the instabilities manifest not as oscillations, but a blow-up of the adjoint potential solution in the H 1 ( ) norm.

Electroosmotic flow in a T-channel
Crossing T-and H-channels are commonly utilized in microfluidics. Applications typically involve mixing of two chemical species [8], purification [16], or fluid identification [12]. However, numerical modeling of electroosmotic flows with slip boundary conditions in such geometrical configurations poses distinctive challenges due to the presence of corner singularities [11]. One immediate consequence is the observation of reduced convergence rates in the approximation of the global solution. A possible remedy http://www.amses-journal.com/content/2/1/15 is to use adaptive finite element methods to help restore the optimal convergence properties of such singular problems [29]. Likewise, adaptive methods can also improve the convergence behavior of the adjoint solution and potentially restore the optimal rates that one may expect when estimating linear QoIs.
We consider below a T-channel geometry. The two upper ends of the T-channel, i,l and i,r , correspond to the left and right inlets, respectively, at which a high potential φ i is prescribed, while the bottom end of the channel o , the flow outlet, is set to the ground potential φ o = 0. The flow is assumed here to be purely electrically driven, in which case Dirichlet pressure boundary conditions p = 0 are considered at the inlet and outlet boundaries. The flow parameters used for the numerical experiments are provided in Table 1.
In the numerical experiments below, we consider the following quantity of interest: (57) http://www.amses-journal.com/content/2/1/15 We also estimate the sensitivity of the QoI with respect to the parameters φ i , φ o , and λ, evaluated in terms of the first derivatives dQ/dφ i , dQ/dφ o and dQ/dλ. We used ten adaptive refinement steps followed by two uniform refinements (for a total of 288,160 dofs) to calculate the reference values of these quantities. These values are reported in Table 2 and were used as exact values to compute numerical errors.
The adaptive strategy for mesh refinement with respect to the QoI is described in Algorithm 1, which has been implemented in libMesh. We show in Figure 5 the horizontal and vertical components of the primal velocity u. We note that the vertical component of the velocity, shown in Figure 5(b), is close to zero near the inlets, but then undergoes a stiff acceleration around the corners. Likewise, we observe in Figure 5(a) the rapid deceleration of the horizontal velocity near the corners. This clearly induces a singular behavior of the solution at the two corners. Note also that the solution is symmetric about the centerline of the vertical channel, as expected, given that the inlet potentials at stations i,r and i,l are equal.

Algorithm 1
Compute the finite element solution to Eq. (39) that either reaches a prescribed mesh size h min or is obtained after a given number of adaptive steps n max using an adaptive meshing strategy based on the dual approach with respect to the QoI Eq. (57).
1: Start step counter n step 2: Compute the finite element solution u h to the problem using a uniform mesh M start of resolution h elem = h start 3: Compute an a posteriori error indicatorẽ h for the QoI based on the adjoint residual error indicator (see [27]) and flag elements to be refined 4: if h elem ≤ h min OR n step > n max then 5: Go to step 11 6   Next, we show the adjoint solutions computed using the adaptive procedure described in Algorithm 1. The vertical velocity, displayed in Figure 6(a), exhibits a parabolic profile that reaches the maximum value along the centerline of the vertical channel and vanishes on its boundaries. Therefore, the presence of corners is solely responsible for the singular behavior in the velocity field. The adjoint potential solution is shown in Figure 6(a). The http://www.amses-journal.com/content/2/1/15 potential is of course singular at the corners due to the geometrical discontinuity and to the fact that the coupling boundary condition, although almost zero everywhere along the boundaries, becomes non-zero near the corners, since ∇ w · λt · (σ * · n)t may not be zero there. This should imply extensive refinement near the corners, as confirmed by the adapted mesh shown in Figure 7. http://www.amses-journal.com/content/2/1/15 Figure 7 Adaptive mesh obtained using adjoint-based error estimates. Note that the elements get refined almost exclusively near the corners due to the singularities in the primal velocity and adjoint potential.
We also used an adjoint method to compute parameter sensitivities for the given QoI to the parameters φ i , φ o , and λ. The advantage of using an adjoint method for sensitivity analysis is that the sensitivity to all three parameters could be found with a single adjoint solve. This is considerably more efficient than using a finite difference or a forward sensitivity method. In addition, we can also combine the adjoint-based mesh refinement and sensitivity analysis for further improvements in the convergence of the sensitivities.
Convergence plots are shown in Figure 8. In particular, the relative error in the quantity of interest estimated using uniform refinement and adjoint-based adaptive refinement is shown in Figure 8(a) against the total number of degrees of freedom (dofs). Relative errors in the estimated sensitivities of the QoI with respect to parameters are displayed in Figure 8(b). We note that the adaptive refinement strategy offers much improved error reduction than uniform refinement for both the estimation of the quantity of interest and its sensitivity derivatives.
In fact, on account of the geometric corner singularities present in the problem, we obtain an inferior convergence rate on using uniform refinement. However, with the adaptive method we obtain a rate of 1.5 (vs dofs) for the QoI, which can be said to be semi-optimal. We had observed earlier that there is a loss of one order in the convergence rate for the forward velocity and adjoint potential for the straight channel problem where there are no corner singularities. We recall that with second-order Lagrange Finite Elements this would result in a convergence rate of 1.5 (N 1 × N 1 2 ) for a linear QoI.

Conclusions
We have presented an analysis of an electroosmotic flow model with slip boundary conditions and its adjoint. The slip boundary conditions require the evaluation of potential http://www.amses-journal.com/content/2/1/15 derivatives on the boundary, which increases the regularity requirements on the potential. We emphasize that a naive enforcement of the standard slip boundary condition leads to an ill-posed adjoint problem (see Section "Ill-posedness of the adjoint problem for http://www.amses-journal.com/content/2/1/15 the naive formulation"). This leads to instabilities in the computed adjoint, illustrated by numerical experiments in Section "Results and discussion". A well-posed adjoint problem can be obtained by modifying the slip boundary condition (u + λ∇φ = 0), i.e. specifying the normal velocity at the wall independently of the potential (u·n = 0, u·t+λ∇φ ·t = 0).
We further proposed a penalty formulation of the forward problem that requires no extra regularity for the potential, and leads to a well-posed, asymptotically consistent adjoint formulation as well. The penalty boundary conditions lead to a weak enforcement of the boundary coupling, allowing us to easily compute the adjoint problem using the adjoint capabilities of libMesh.
Finally, we presented numerical experiments for a simple straight channel microflow and a more challenging T-channel flow. The convergence results for the straight channel problem indicate that the primal velocity and the adjoint potential converge at suboptimal rates due to the nature of the coupling between the potential and the velocity. For the T-channel, we presented QoI computation and QoI adjoint sensitivity results for a practical engineering QoI. We observed a loss of convergence order due to the singularities in the T-channel geometry, and substantial improvements in the rate on using an adjoint-based adaptive method. However, the fully optimal convergence rate for the QoI could not be achieved, possibly due to the convergence properties of the adjoint potential.
Future work will involve the application of more sophisticated adjoint-based error estimators to further improve the convergence properties of the approximate solutions and obtain reliable a posteriori error estimates for complex applications.

Appendix A: Well-posedness of coupled formulation
Let Z and X denote two reflexive Banach spaces. Consider the following abstract linear one-way coupled problem: where A(·, ·), B(·, ·), and C(·, ·) and continuous bilinear forms, and F(·) and G(·) are continuous linear forms.
If A(·, ·) and B(·, ·) satisfy the inf-sup conditions on Z × Z and on X × X, respectively, then the above problem can be solved sequentially: First solve (58a) for φ ∈ Z. Then, since C(φ, ·) is a continuous linear form on X, Eq. (58b) can be solved for u ∈ X.
The fact that one-way coupled problems (with bounded coupling) are well-posed is probably well known. We now provide an extension of this result for the aggregated bilinear form A(φ, ψ) + B(u, v) + C(φ, v).
Choosing which upon invoking (59b) yields φ = 0. Next, we derive the a priori estimates. From inf-sup stability of A(·, ·) and (58a) we obtain: Similarly, using inf-sup stability of B(·, ·), (58b), and continuity of C(·, ·) we obtain: Corollary. The weak formulation of the slip-BC EOF model, see (25), is well-posed. We have, Moreover, the following a priori estimates hold: where γ B is the inf-sup constant for the bilinear form B((w, p), (v, q)), λ as required by Proposition 1 and c( ) is as given by Eq. (19).
Corollary. The weak adjoint problem given by Eq. (28) is well-posed. Moreover, the following a priori estimates hold: