# Effect of interfacial Maxwell stress on time periodic electro-osmotic flow in a thin liquid film with a flat interface

Manik Mayur, Sakir Amiroudine, Didier Lasseux, Suman Chakraborty

► To cite this version:

Manik Mayur, Sakir Amiroudine, Didier Lasseux, Suman Chakraborty. Effect of interfacial Maxwell stress on time periodic electro-osmotic flow in a thin liquid film with a flat interface. Electrophoresis, 2014, Special Issue: Fundamentals 2014, 35 (5), pp.670-680. 10.1002/elps.201300236. hal-03830429

HAL Id: hal-03830429

https://hal.science/hal-03830429

Submitted on 26 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effect of interfacial Maxwell stress on time periodic electro-osmotic flow in a thin liquid film with a flat interface

Electro-osmotic flows (EOF) have seen remarkable applications in lab-on-a-chip based microdevices owing to their lack of moving components, durability, and nondispersive nature of the flow profiles under specifically designed conditions. However, such flows may typically suffer from classical Faradaic artifacts like electrolysis of the solvent, which affects the flow rate control. Such a problem has been seen to be overcome by employing time periodic EOFs. Electric field induced transport of a conductive liquid is another nontrivial problem that requires careful study of interfacial dynamics in response to such an oscillatory flow actuation. The present study highlights the role of electric field generated Maxwell stress and free surface potential along with the electric double layer thickness and forcing frequency, toward influencing the interfacial transport and fluid flow in free-surface electro-osmosis under a periodically varying external electric field, in a semi-analytical formalism. Our results reveal interesting regimes over which the pertinent interfacial phenomena as well as bulk transport characteristics may be favorably tuned by employing time varying electrical fields.

**Keywords:**

Electro-osmosis / Maxwell stress / Thin film / Time-periodic

Manik Mayur<fnref n="1" />

Sakir Amiroudine<fnref n="1" />

Didier Lasseux<fnref n="1" />

Suman Chakraborty<fnref n="2" />

## 1 Introduction

Flow actuation in microchannels due to an externally applied electric field has found remarkable applications in lab-on-a-chip based microfluidics devices and systems [1–7]. In many situations, such flows are effectively actuated by interactions between a wall-adhering charged layer (also known as the electrical double layer) and an externally applied electrical field, resulting in so-called EOFs, which have an enduring advantage over classical fluidic pumps in a sense that these do not require any moving components. EOFs in microchannels with axially invariant interfacial potential typically have plug like (uniform) velocity profile which reduces the possibility of species dispersion and hence may turn out to be effective for the transport of biological species in microdevices, provided there is no sample overheating due to Joule effects.

The most commonly form of used constant electrical field to invoke EOF has some inherent problems due to electro-chemical reactions at the electrodes. Such problems include formation of hydrogen and oxygen bubbles because of the

hydrolysis of water which leads to fluctuations in the flow rate in microfluidic devices and sometimes eventually stops the flow altogether [8]. The other problem associated with electro-chemical reactions at the electrodes is the change of pH at electrode reservoirs [9]. Such a phenomenon creates a pH gradient in the channel affecting the substrate zeta potential and eventually the mobility of ions in the liquid. Some previous studies have shown that the use of time periodic EOF is instrumental in diminishing the Faradaic reactions at the electrodes and to achieve effective spatio-temporal control over liquid flow and mixing [10–13]. There have been a multitude of studies discussing various novel applications of AC-EOF including flow actuation [14–19], energy storage [20–22] etc.

Although EOF has been proven to be the preferred mechanism for liquid transport in microchannels, one of the biggest shortcomings of its usage in mainstream applications comes from its basic requirement of the concerned liquid to be electrically conductive. In most of the microdevices, the use of electrical field for fluid transport is limited by the non-conductive nature of a large group of fluids. For example, although several biological liquids (e.g. blood serum, insulin, etc.) show some degree of ionic conductivity, their low degree of ionic solvation renders them unresponsive to EOF. However, there have been some studies involving the transport of nonconductive liquid with the help of an immiscible conductive liquid through shear transfer at the fluid–fluid interface [7, 23–26]. The mentioned studies have been performed in a two-layer EOF, while considering classical interface

<footnote n="1">¹Université Bordeaux 1, Institut de Mécanique et d'Ingénierie –I2M, Pessac, France</footnote>

<footnote n="2">²Department of Mechanical Engineering, Indian Institute of Technology, Kharagpur, India</footnote>
matching conditions such as continuity of the velocity and hydrodynamic shear stress at the fluid–fluid interface.

The presence of a fluid–fluid interface in an EOF necessitates the characterization of the parameters responsible for instability of such systems, for which one has to ascertain the basic or unperturbed state solution of the system. This is characterized by several intricacies, the modeling of which is not trivial. Those intricacies are attributed to the existence of a free surface on which appropriate considerations on Maxwell stress need to be invoked [27]. This is in sharp contrast with flows in between rigid boundaries for which one imposes a velocity boundary condition (typically no slip at the walls) irrespective of any consideration on Maxwell stress. Moreover, in case of a symmetric flow between two rigid boundaries, the hydrodynamic stress and the Maxwell stress individually become necessarily zero at the centerline (because of the centerline symmetry). However, for a free surface flow, the total stress (hydrodynamic and Maxwell) at the air–liquid interface needs to be relaxed, which suggests that in the presence of a nonzero Maxwell stress, an equal and opposite hydrodynamic stress will be present contrary to the classical hydrodynamic stress-free notion of free-surface dynamics.

The present work is dedicated to bring out contrasting features of time periodic EOF in a thin liquid layer exposed to a low viscosity and low permittivity gas phase and the corresponding features of time periodic EOF between two rigid boundaries. In order to bring out a common comparative basis between the two cases, an unperturbed interface for the former case is considered here, which technically represents the basic state of a general stability problem. In an effort to bring out the interesting interplay between the characteristic length scale of the liquid layer and the characteristic electric double layer (EDL) length scale (Debye length), various orders of the liquid layer thickness to Debye length ratios are addressed. For further generalizations, the possibilities of addressing high zeta potentials by invoking the nonlinear Poisson–Boltzmann equation (PBE) are considered. The analytical expressions for the electric potential and velocity field are derived bringing out important implications of Maxwell stresses at the free surface. The present study opens new perspectives for further studies on instabilities associated with a two-layer EOF in the presence of a time periodic electric field.

## 2 Mathematical modeling

When an aqueous electrolyte comes into contact with a substrate like silica, glass, polymers, and some other chemically active substrates, there is a possibility of a series of ionic exchanges such as protonation, deprotonation, adsorption, and some chemical reactions at the solid–liquid interface. Such ionic exchanges, after attaining equilibrium, leave the surface charged [28].

The system under study consists of a thin electrolyte film spread over a rigid solid substrate exposed to an inert gaseous atmosphere (see Fig. 1). The film thickness is denoted by d. The dynamics of such a film is studied under the effect of a longitudinal oscillating electric field, $E_{app} = E_0 \text{Sin}(\omega t)$,

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematics diagram of the time periodic EOF system under study.</figcaption></figure>

where $E_0$ is the amplitude and $\omega$ is the frequency. The electrolyte concentration in this study is considered to be small enough in magnitude (~0.1–10 mM) in order to neglect the liquid property changes due to Joules heating [29] even in the case of an applied electric field of large amplitude. The low electrolyte concentration also avoids complexities in flow modeling by reducing the nonlinear dependence of electrophoretic mobility of ions on the sparse space charge distribution [30–36]. The solid substrate zeta potential is considered to be $\zeta_b$, which is a function of the substrate–fluid interaction, ionic concentration and pH of the solution [37]. The liquid surface exposed to a gaseous environment develops a charge, which is a function of various parameters like ionic concentration, pH of the solution and valence of the ions involved [38–40]. The associated zeta potential ($\zeta_{Interface}$) has been found to vary over a wide range in the reported literature [41–43].

### 2.1 Electric potential field

The thin electrolyte film under investigation contains a space charge distribution, which is electrostatically bounded by the EDL at the solid substrate and an interfacial potential at the free surface that is a function of electrochemical parameters of the system.

The transport of the charged species can be expressed by the Poisson-Nernst-Planck equation:

$$ \frac{\partial C_i}{\partial t} + \vec{\nabla} \cdot (\vec{u} C_i) = \vec{\nabla} \cdot (D_i \vec{\nabla} C_i) + \frac{z_i e}{k_B \theta} \vec{\nabla} \cdot (D_i C_i \vec{\nabla} \Phi), \quad (1) $$

where $C_i$ is the concentration of the i-th species, $\Phi = \Phi_{app} + \Phi_{sc}$ is the electrostatic potential field in the fluid, which is a sum of the potential field due to the space charge distribution ($\Phi_{sc}$) and the externally applied potential field ($\Phi_{app} = -\int E_{app} dx$), $z_i$ is the valence of the i-th ionic species, $D_i$ is the diffusivity of the i-th ionic species, $k_B$ is the Boltzmann constant, $\theta$ is the liquid temperature, $e$ is the electronic charge, and $\vec{u}$ is the fluid velocity. The forcing electric field frequency is considered to be smaller than the double layer relaxation frequency, hence the associated space charge distribution is assumed to be quasi-steady, independent of
the time-varying applied potential gradient ($\phi_{app}$). Also, in flow geometries having characteristic length scales of the order of micro and nanometers, the faster charge relaxation due to ionic diffusivity over small distances (thin films) and a weaker advective transport in a viscosity dominated regime, the diffusive relaxation of the free charge distribution dominates over their advective transport. Moreover, upon assuming a laminar flow situation ($\vec{u} = (u, 0, 0)$) and a nonreactive system connected to an infinite reservoir of the involved electrolytes with channel walls impermeable to ionic exchange, the concentration gradient in the longitudinal direction can be considered negligible ($\partial C_i / \partial x \sim 0$). The quasi-steady nature of the space charge distribution along with isotropic liquid properties reduces the Nernst-Planck equation to:

$$0 = D_i \nabla^2 C_i + \frac{z_i e D_i}{k_B \theta} \vec{\nabla} \cdot (C_i \vec{\nabla} \phi_{sc}), \quad (2)$$

which results into a Boltzmann distribution of the charged species as:

$$C_i = C_{0,i} \exp \left( -\frac{z_i e \phi_{sc}}{k_B \theta} \right), \quad (3)$$

where $C_{0,i}$ is the bulk ionic concentration of the i-th species. Using mean field approximations, the free charge density $\rho_e$ can be expressed in terms of volume averaged ion concentrations as:

$$\rho_e = \sum_i e z_i C_i. \quad (4)$$

This free charge density $\rho_e$ is related to the electric displacement field vector ($\vec{D}$) by the differential form of Gauss's Law as:

$$\vec{\nabla} \cdot \vec{D} = \vec{\nabla} \cdot (\vec{D}_{app} + \vec{D}_{sc}) = \rho_e, \quad (5)$$

where $\vec{D}_{app}$ is the electric displacement vector due to the applied electric field and $\vec{D}_{sc}$ is the electric displacement vector due to the space charge distribution. For a linear and isotropic dielectric medium under a time-dependent electric field, the electric displacement field is related to the electric field, $\vec{E}$ as:

$$\vec{D} = \varepsilon(\omega) \vec{E}(\omega), \quad (6)$$

where, $\varepsilon(\omega)$ is the complex permittivity of the system, $\vec{E}(\omega) = \vec{E}_{app}(\omega) + \vec{E}_{sc}$ is the electric field in the system and $\omega$ is the frequency of the applied field. The relaxation time of ionic species in a homogeneous dielectric media is given by the Maxwell-Wagner-O'Konski relaxation time [44] and the corresponding frequency can be written as:

$$\omega_{MWO} = \frac{1}{\tau_{MWO}} \approx \frac{2 D_i}{\lambda_D^2}, \quad (7)$$

where, $\lambda_D$ is the Debye length (which is discussed further). So in the limit of $\omega < \omega_{MWO}$, which sets the upper limit on the frequency in this study, charges relax to the equilibrium distribution faster than the time-dependent external perturbation and hence polarization effects are negligible and the permittivity can be assumed constant as $\varepsilon(\omega) \approx \varepsilon$. Also, since $\vec{E}_{app}$ is the time-varying homogeneous applied electric

field, $\vec{\nabla} \cdot \vec{D}_{app} = \varepsilon \vec{\nabla} \cdot \vec{E}_{app} = 0$. Hence, the electric displacement field generated by the space charge distribution can be written as:

$$\vec{\nabla} \cdot \vec{D}_{sc} = \varepsilon \vec{\nabla} \cdot \vec{E}_{sc} = \rho_e, \quad (8)$$

From electrostatics it is known that $\vec{E}_{sc} = -\vec{\nabla} \phi_{sc}$. Hence, the Poisson equation for the space charge potential is:

$$-\vec{\nabla} \cdot (\varepsilon \vec{\nabla} \phi_{sc}) = \sum_i e z_i C_i. \quad (9)$$

From Eqs. (3) and (9), the PBE can be written as:

$$-\vec{\nabla} \cdot (\varepsilon \vec{\nabla} \phi_{sc}) = \sum_i z_i e C_{0,i} \exp \left( -\frac{z_i e \phi_{sc}}{k_B \theta} \right). \quad (10)$$

For a symmetric binary electrolyte ($|z_+| = |z_-| = z$) and a solvent with constant permittivity, Eq. (10) reduces to:

$$\varepsilon \nabla^2 \phi_{sc} = 2 \varepsilon e C_0 \sinh \left( \frac{z e \phi_{sc}}{k_B \theta} \right), \quad (11)$$

where $C_0$ is the neutral bulk ionic concentration of the solution. In order to isolate the electric effects as the dominating forcing mechanism, the system under study is considered to have large lateral extents, which results into negligible x-gradients as compared to the y-gradients. Upon non-dimensionalizing the PBE using $\Phi_{sc} = \phi_{sc}/\zeta_b$ and $Y = \gamma/d$, Eq. (11) leads to:

$$\frac{d^2 \Phi_{sc}}{d Y^2} = \beta \sinh (\alpha \Phi_{sc}), \quad (12)$$

where $\beta = d^2/(\lambda_D^2 \alpha)$, $\alpha = (z e \zeta_b)/(k_B \theta)$ is the ionic energy parameter which measures the relative strength of the electrostatic energy of ions with respect to the thermal energy of ions and $\lambda_D = \sqrt{(\varepsilon k_B \theta)/(2z^2 e^2 C_0)}$ is the Debye length. The associated boundary conditions are:

$$\Phi_{sc}(0) = 1, \quad \Phi_{sc}(1) = Z_R, \quad (13)$$

where, $Z_R = \zeta_{Interface}/\zeta_b$. Interfacial potential value depends upon electrochemistry of the problem. The charging of an air-electrolyte interface is a complex phenomenon and its mechanism it not clearly understood so far. However, the widely accepted mechanism contributing to its surface charge is adsorption of ions at the interface. This surface charge is neutralized by a diffused distribution of oppositely charged ions (counter-ions), thus creating a double layer structure. The only difference between a solid-liquid and liquid-gas double layer is that in the case of solid-liquid interface, the Stern layer is bound to the solid surface and hence does not move when an external electric field is applied. However, at the free surface, it is free to move. Commonly such electrolytic systems have low substrate zeta potential corresponding to $\alpha < 1$, which for a monovalent symmetric electrolyte corresponds to $\zeta_b < 25$ mV at 25°C. In that case, Eq. (12) can be linearized as (also known as the Debye-Hückel linearization):

$$\frac{d^2 \Phi_{sc}}{d Y^2} = \frac{\Phi_{sc}}{D e^2}, \quad (14)$$
where $De = \lambda_D/d$ is the Debye number that represents the relative extent of the EDL with respect to the liquid film thickness. Under this formalism, the boundary value problem (BVP) for $\lambda_{sc}$ can be solved analytically to obtain the following closed form solution:

$$ \Phi_{sc}(Y) = \frac{1}{\sinh\left(\frac{1}{De}\right)} \left\{ Z_R \sinh\left(\frac{Y}{De}\right) + \sinh\left(\frac{1-Y}{De}\right) \right\}. \quad (15) $$

However, for larger wall zeta-potential systems ($\alpha \ge 1$), Eqs. (12) and (13) have to be solved. This system of equations is a BVP having two Dirichlet boundary conditions along with a nonlinear second-order ordinary differential equation, for which no tractable analytical solution is available. Hence, to obtain the resulting potential distribution one has to resort to numerical tools. It is also observed that for thin EDLs, $\beta$ can be larger than 1, making the differential equation very stiff with large variations over short distances (for example, within EDL). To solve such a highly stiff BVP, the Automated Continuation with Deferred Correction method [45] is used here. This method is based on implicit Runge-Kutta scheme on a Lobatto grid. A Lobatto grid includes boundary points in the integration interval facilitating the use of an adaptive mesh. Such an adaptive mesh is extremely useful for systems having large variations over small distances like a boundary layer or an EDL.

The electric potential due to the externally oscillating applied electric field ($\Phi_{app}$) can be written in the dimensionless form as:

$$ \Phi_{app}(X) = - \int \frac{E_0 d}{\zeta_b} \sin\Theta \, dX = - \frac{X}{E_R} \sin\Theta, \quad (16) $$

where, $\Phi_{app} = \Phi_{app}/\zeta_b$, $X = x/d$, $\Theta = \omega t$ are the dimensionless parameters and $E_R = \zeta_b/E_0 d$ is the ratio of the relative strength of electric field due to space charge distribution and the applied electric field. The net electric potential in the system can be written as a sum of the potential field due to the space charge distribution ($\Phi_{sc}$) and the externally applied electric field ($\Phi_{app}$). While, the space charge potential distribution for large wall zeta potential systems ($\alpha \ge 1$) has to be obtained numerically, for low wall zeta potential systems ($\alpha < 1$), the net potential field is obtained analytically and can be written as:

$$ \begin{aligned} \Phi(X, Y) = \Phi_{app}(X) + \Phi_{sc}(Y) \\ = -\frac{X}{E_R}\sin\Theta + \frac{1}{\sinh\left(\frac{1}{De}\right)} \left\{ Z_R \sinh\left(\frac{Y}{De}\right) \right. \\ \left. \qquad + \sinh\left(\frac{1-Y}{De}\right) \right\}, \end{aligned} \quad (17) $$

where $\Phi = \phi/\zeta_b$ is the dimensionless total electric potential of the system.

## 2.2 Velocity field

The oscillating electric field along with the space charge distribution induces a time-dependent Maxwell stress ($\Sigma^M$) in the liquid, which is represented as:

$$ \vec{\Sigma}^M = - \frac{\varepsilon |\vec{E}|^2}{2} \vec{I} + \varepsilon \vec{E} \otimes \vec{E} \quad (18) $$

Along with the hydrodynamic stress tensor ($\vec{\Sigma}^H$) for a Newtonian fluid, the total stress tensor ($\vec{\Sigma}^T$) can be expressed as [27]:

$$ \begin{aligned} \vec{\Sigma}^T = \vec{\Sigma}^H + \vec{\Sigma}^M = & - \left( p + \frac{\varepsilon |\vec{E}|^2}{2} \right) \vec{I} + \mu (\vec{\nabla} \vec{u} + \vec{\nabla} \vec{u}^T) \\ & + \varepsilon \vec{E} \otimes \vec{E} \end{aligned} \quad (19) $$

where $\vec{u}$ is the liquid velocity, $p$ is the hydrostatic pressure in the liquid, $\mu$ is the dynamic viscosity of the solvent, $\vec{I}$ is the unit tensor and $\vec{E} = -\vec{\nabla}\phi$ is the total electric field in the system.

Considering an incompressible flow, the conservation of mass and momentum results into the following equations, respectively:

$$ \vec{\nabla} \cdot \vec{u} = 0, \quad (20) $$

$$ \begin{aligned} \rho \left( \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \vec{\nabla}) \vec{u} \right) &= \vec{\nabla} \cdot \vec{\Sigma}^T = -\vec{\nabla} p + \mu \nabla^2 \vec{u} + \vec{\nabla} \cdot \vec{\Sigma}^M \\ &= -\vec{\nabla} p + \mu \nabla^2 \vec{u} + \varepsilon \vec{\nabla} \phi \nabla^2 \phi. \end{aligned} \quad (21) $$

While considering the continuity of shear and normal stresses at the planar free surface the boundary conditions at the free surface can be written respectively as:

$$ \left[ \vec{t} \cdot \vec{\Sigma}^T \cdot \vec{n} \right]_{y=d} = 0, \quad (22) $$

$$ \left[ \vec{n} \cdot \vec{\Sigma}^T \cdot \vec{n} \right]_{y=d} = 0, \quad (23) $$

where $\vec{t}$ and $\vec{n}$ are the tangent and normal vectors at the free surface, respectively.

At the solid wall, a nonpolar (e.g. PDMS, silanized glass etc.) substrate might induce a slip velocity on the polar solvent due to the ensuing hydrophobicity [46]. However, in the mentioned work, it has been stated that the roughness of the nonpolar substrate is generally greater than the threshold roughness of the induced slip. This results into a small slip length, which for all general purposes does not affect the flow characteristics of the polar liquids over such nonpolar substrates. Hence, considering a no slip condition at the solid wall gives:

$$ \vec{u}(0, t) = \vec{0}. \quad (24) $$
Equation (22) expresses the absence of shear stress at the interface while Eq. (23) expresses the absence of normal stress owing to a flat interface. The velocity profile is considered to be fully developed in the x-direction. Upon nondimension-alizing the flow velocity as, $U = u/U_{ref}$, where $U_{ref} = v/d$ is the momentum diffusion velocity scale. Without an external pressure gradient, the dimensionless momentum conserva-tion Eq. (21) reduces to:

$$ 
\begin{aligned}
Wo^2 \frac{\partial U}{\partial \Theta} &= \frac{\partial^2 U}{\partial Y^2} + \gamma_R E_R \frac{d\Phi_{app}}{dX} \frac{d^2\Phi_{sc}}{dY^2} \\
&= \frac{\partial^2 U}{\partial Y^2} - \gamma_R \sin\Theta \frac{d^2\Phi_{sc}}{dY^2},
\end{aligned} 
\quad (25) $$

$$ 0 = -\frac{dP}{dY} + \gamma_R E_R \frac{\partial \Phi}{\partial Y} \frac{\partial^2 \Phi}{\partial Y^2} = -\frac{dP}{dY} + \gamma_R E_R \frac{d\Phi_{sc}}{dY} \frac{d^2\Phi_{sc}}{dY^2}, \quad (26) $$

where $Wo = \sqrt{\omega d^2/\nu}$, is the Womersley number expressing the relative strength of temporal inertial force over the viscous dissipation force, $\gamma_R = (\epsilon E_0 \zeta_b d)/(\rho v^2) = -U_{HS}/U_{ref}$ is the relative strength of the electrical body forces over the viscous dissipation force and is henceforth referred to as electroviscous ratio and $U_{HS} = -\epsilon \zeta_b E_0/\mu$ is the Helmholtz-Smoluchowsky velocity. The boundary conditions in dimensionless form (see Eqs. (22)-(24)) are, respectively:

$$ \left. \frac{\partial U}{\partial Y} \right|_{Y=1} - \gamma_R \sin\Theta \left. \frac{d\Phi_{sc}}{dY} \right|_{Y=1} = 0, \quad (27) $$

$$ 
\begin{aligned}
P(1, \Theta) &= -\frac{\gamma_R E_R}{2} \left( \left( \frac{\partial \Phi}{\partial X} \right)^2 - \left( \frac{\partial \Phi}{\partial Y} \right)^2 \right)_{Y=1} \\
&= -\frac{\gamma_R E_R}{2} \left( \frac{\sin^2 \Theta}{E_R^2} - \left( \frac{d\Phi_{sc}}{dY} \right)^2 \right)_{Y=1},
\end{aligned} 
\quad (28) $$

$$ U(0, \Theta) = 0. \quad (29) $$

The solution of the system of Eqs. (25), (27), and (29) can be obtained by decomposing the velocity into time-dependent and space-dependent functions as:

$$ U(Y, \Theta) = F(Y) G(\Theta) = \operatorname{Im} (F(Y) e^{i\Theta}). \quad (30) $$

Upon substituting Eq. (30) in Eqs. (25), (27) and (29), an ordinary differential equation in $F(Y)$ is obtained as:

$$ i Wo^2 F = \frac{d^2 F}{dY^2} - \gamma_R \frac{d^2 \Phi_{sc}}{dY^2}. \quad (31) $$

The corresponding boundary conditions are:

$$ \left. \frac{dF}{dY} \right|_{Y=1} - \gamma_R \left. \frac{d\Phi_{sc}}{dY} \right|_{Y=1} = 0, \quad F(0) = 0. \quad (32) $$

For a generalized potential distribution, the resulting velocity profile can be obtained as:

$$ U(Y, \Theta) = \gamma_R \operatorname{Im} \left[ \begin{aligned} & \Phi_{sc}(Y) e^{i\Theta} - \frac{\cosh(\sqrt{i}(1-Y)W_0)}{\cosh(\sqrt{i}W_0)} \Phi_{sc}(0) e^{i\Theta} \\ & - \frac{\sinh(\sqrt{i}YW_0)}{\cosh(\sqrt{i}W_0)} e^{i\Theta} \sqrt{i} W_0 \\ & \times \int_0^1 \cosh(\sqrt{i}W_0(1-\eta)) \Phi_{sc}(\eta) d\eta \\ & + \frac{\sqrt{i}W_0}{2} e^{i\Theta} \left\{ e^{\sqrt{i}YW_0} \int_0^Y e^{-\sqrt{i}\eta W_0} \Phi_{sc}(\eta) d\eta \right. \\ & \left. - e^{-\sqrt{i}YW_0} \int_0^Y e^{\sqrt{i}\eta W_0} \Phi_{sc}(\eta) d\eta \right\} \end{aligned} \right], \quad (33) $$

which for a space charge potential obtained by Debye-Hückel linearization can be simplified analytically as:

$$ U(Y, \Theta) = U_a(Y) \sin(\Theta + \Delta(Y)), \quad (34) $$

where, $U_a(Y)$ is the amplitude of the velocity and $\Delta(Y)$ is the phase of the velocity. Detailed expressions of the amplitude and the phase are given in the Appendix.

To obtain the velocity field for a generalized potential one needs to perform the involved integrations in Eq. (33) numerically. In this study, the QUADPACK [47] numerical integration package was used, which effectively handles the integration of complex functions with the help of an adaptive automatic integration algorithm using Gauss-Kronrod rule [47]. The pressure distribution in the system can be obtained by solving the Eq. (26) along with the boundary condition from Eq. (28) as:

$$ P(Y, \Theta) = \frac{\gamma_R}{4E_R} \cos 2\Theta - \frac{\gamma_R}{4E_R} + \frac{\gamma_R E_R}{2De^2 \sinh^2\left(\frac{1}{De}\right)} \\ \times \left\{ Z_R \cosh\left(\frac{Y}{De}\right) - \cosh\left(\frac{1-Y}{De}\right) \right\}^2. \quad (35) $$

It is observed from Eqs. (25) and (26) that, in the absence of an externally applied pressure gradient, the resulting pressure and velocity distributions are decoupled and are parametrically dependent on the applied potential bias and the space charge potential.

### 3 Results and discussions

It can be seen from Eq. (27) that the interfacial shear jump condition is no more a classical hydrodynamic stress free boundary condition. The Maxwell stress, by the virtue of Coulombic force on the space charge distribution and the
free surface potential, contributes to the interfacial dynamics as well. To demonstrate the effects of Maxwell stress on EOFs, a comparison of flow characteristics between EOF under a Hele-Shaw configuration (parallel plate flow) [11, 13] and the free surface EOF has been carried out. Although being two completely different physical systems the mentioned problems share a similar physical modeling in terms of governing equation and boundary conditions. The most general modeling of a parallel plate EOF system involves Stokes equations with a Helmholtz-Smoluchowsky slip velocity at the walls and a symmetry condition applied on the velocity at the channel centerline. Such a system has been very successful in modeling microscale EOF owing to the thin EDL length scales as compared to the channel length scales. However, for systems where the EDL thickness becomes comparable to the channel length scales, one can add an electrical body force term in the Stokes equations, the involved electric charge distribution being obtained from the electrical potential distribution which results as a solution of the PBE [11, 13]. In the second case, the wall boundary condition is a no slip condition along with the previously used symmetry condition at the channel centerline. Similarly, the modeling of a thin free surface flow involves a Stokes equation with no slip condition at the wall and a stress-free condition at the free surface. In the absence of various surface phenomena such as surface tension variation, phase change, and mass transfer, the stress-free condition at the free surface ($\partial u/\partial y = 0$) is mathematically the same as a symmetry condition at the channel center ($\partial u/\partial y = 0$) in a parallel plate EOF [11, 13]. Hence, the commonly studied parallel plate EOF can be put forth as a case without the Maxwell stress, which in the simplified form adds a term involving transverse gradient of electric potential at the free surface (see Eq. (29)). As a consequence, it is seen that, in the presence of very thin EDLs ($De = 0.01$) and no interfacial potential ($Z_R = 0$), the free surface system behaves as a classical half channel EOF system where the electrical effects are localized close to the wall and the far field electric potential is zero. A validation of the present free surface model developed in the previous sections is presented through a comparison with a Hele-Shaw EOF, where the half channel velocity profile was taken from the seminal work of Dutta and Beskok [13] (see Fig. 2). Their usage of Debye length $\lambda_D$ as the characteristic length scale as compared to the film thickness $d$ used in this study leads to different dimensionless groups that can be mathematically adjusted without losing any physical details for comparing the results. Their various parameters can be recovered from the present parameters as:

$$ \kappa = Wo \times De; \eta = De; U_{HS} = -U_{\text{ref}}(\text{for } \gamma_R = 1), \quad (36) $$

where $\kappa$ is the dimensionless frequency, $\eta$ is the dimensionless Debye length, and $U_{HS}$ is the Helmholtz-Smoluchowski slip velocity used as the reference velocity by Dutta and Beskok [13]. The excellent agreement between the velocity profiles obtained from the present model and the one extracted from Dutta and Beskok [13] as reported in Fig. 2, justifies the success of non-Maxwell stress model in thin EDL systems.

<figure><img src="image_3.png" /><figcaption>Figure 2. Velocity profiles comparison with Dutta and Beskok [13], where $\kappa = Wo \times De$ is the dimensionless frequency used by Dutta and Beskok. Other fixed parameters are $De = 0.01$, $Z_R = 0$, $\gamma_R = 1$, $\Theta = \pi/2$.</figcaption></figure>

For cases where the EDL extent is comparable to the characteristic dimensions of the flow, which in the present case is the film thickness $d$, the effect of Maxwell stress modifying the shear stress balance criterion at the surface boundary is more pronounced. This leads to changes on the transverse velocity gradients as a function of the surface potential and potential field present in the bulk. In the absence of transverse velocity component, the distribution of such a gradient also corresponds to the absolute value of vorticity ($\omega_{abs}$) field defined as $\omega_{abs} = |\vec{\nabla} \times \vec{U}| = -\partial U/\partial Y$, which can account for vorticity induced flow instabilities in the system. To study the effect of Maxwell stress model on vorticity distribution, a comparison between free-surface EOF models with and without Maxwell stress was performed (see Fig. 3). The vorticity profile of non-Maxwell stress model was taken from Chakraborty

<figure><img src="image_4.png" /><figcaption>Figure 3. Vorticity profile comparison of Maxwell stress and non-Maxwell stress models [11]. The values of other fixed parameters are $De = 0.5$, $\alpha = 0.4278$, $Z_R = 0$, $\gamma_R = 1$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Velocity profile with different values of the Debye number ($De$) and ionic energy parameter ($\alpha$) for fixed values of Womersley number $Wo = 10$, electroviscous ratio $\gamma_R = 1$, $Z_R = 0$ and time, $\Theta = \pi/2$.</figcaption></figure>

and Srivastava [11], whose dimensionless parameters relate to the present parameters as mentioned in Eq. (36). It is seen from the figure that the neglecting Maxwell stress can lead to significant underestimation of the vorticity strengths in the bulk. Also, for thick EDLs the strong presence of velocity gradients renders Helmholtz-Smoluchowsky velocity a non-ideal slip condition at the liquid-solid interface. In order to strengthen this argument, a further discussion on velocity scales is provided with different values of $\gamma_R$. The resulting velocity gradients lead to dispersion in thick EDL flows over the time and with the help of this study one can identify the regimes of applied frequency and electrolyte characteristics to avoid the dispersion effects in a time periodic EOF.

As the ionic energy parameter, $\alpha$, defines the relative strength of electrostatic energy of ions over their thermal

energy, the changes in the velocity distribution due to the relative interplay between the mentioned competing energies is significant. To demonstrate the effect of $\alpha$, four cases have been considered with different values of the EDL thickness ($De$) and ionic energy parameter ($\alpha$) (see Fig. 4). It is seen that for the lower values of ionic energy parameter ($\alpha = 0.1$) the difference in the velocity distribution is significant for higher EDL thickness ($De = 0.1$). However, in the case of thin EDL, the difference is not noticeable. This also explains the reason why the Debye-Hückel linearization of PBE, which is valid for low values of $\alpha$, successfully models the thin EDL potential distribution even over a large range of $\alpha$. The above observation can also be explained by the fact that, when the thermal motion of ions dominates ($\alpha < 1$), the electric actuation is diffused strongly in the EDL and hence the local velocity magnitude is reduced. This diffusion of the electrical actuation is enhanced if the EDL is thicker as the presence of a greater amount of space charge distribution aides in more effective distribution of velocity in the liquid bulk.

The electroviscous ratio, $\gamma_R$, see Eq. (25) is directly proportional to the amplitude of the applied oscillating electric field and hence, increment in which amplifies the flow velocity profile for a given fluid. This is clearly highlighted in Fig. 5 for the values $De = 0.01$ and $De = 0.1$. It can also be seen that the maximum value of the velocity obtained for a thin EDL case ($De = 0.01$) is greater than in the case with thicker EDL ($De = 0.1$), which can be attributed to the fact that in the case of thick EDLs, the electrical energy is spent mobilizing a greater space charge distribution than in the case of a thin EDL. This leads to a lower maximum kinetic energy of the liquid within the EDL. Moreover, $\gamma_R$ also represents a ratio of the classically used Helmholtz-Smoluchowsky velocity ($U_{HS}$) as a reference velocity scale to the current velocity scale. It is interesting to note that $U_{HS}$ can be used as a reference velocity scale ($\gamma_R = 1$) when (i) EDL is thin (see Fig. 5),

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Velocity distribution over the film thickness with different values of electro-osmotic number, $\gamma_R$ for $\alpha = 5$, $Z_R = 0$, $Wo = 10$ at $\Theta = \pi/2$ with (A) $De = 0.01$ and (B) $De = 0.1$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Velocity distribution over the film thickness with different values of Wo and for fixed values of α = 5, Z<sub>R</sub> = 0 and γ<sub>R</sub> = 1 at Θ = π/2 with (A) De = 0.01 and (B) De = 0.1.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

(ii) forcing frequency is smaller than the viscous relaxation frequency (see Fig. 6). In remaining cases, U<sub>HS</sub> overestimates the characteristic velocity sometimes by an order of magnitude. Also, one can observe from Fig. 5 that the velocity oscillations propagate as a damped wave into the bulk. The extent of this bulk penetration of the damped velocity oscillations seem to be unaffected by the value of γ<sub>R</sub>. Taking a cue from the classical Stokes second problem, the typical penetration depth of momentum diffusion of oscillatory flows in a liquid can be estimated by the Stokes penetration depth δ<sub>s</sub>. Typically δ<sub>s</sub> which is a function of the forcing frequency as δ<sub>s</sub>/d ~ √v/(ωd<sup>2</sup>) = 1/Wo, where Wo is the Womersley number, which also represents the relative strength of temporal inertia over viscous force (see Eq. (25)). This inverse relationship between the Stokes penetration depth and Womersley number, Wo is also observed in Fig. 6 which shows the velocity profiles as a function of Wo.

It is observed that decreasing Wo increases the Stokes penetration depth and in turn enhances the sharp changes in the velocity gradients near the wall. This also results in the stronger diffusion of the near-wall vorticity in the liquid bulk. However, upon increasing Wo the viscous time lag of momentum propagation in the liquid also increases, which leads to inflection points in the velocity profiles. Presence of such inflection points leads to extrema in the vorticity distribution yielding possible fluid mixing and instability [48]. It is also observed that the magnitude and extent of vorticity propagation is a strong function of the EDL thickness. This can be explained in terms of electrical energy transfer over the extent of space charge distribution in the liquid. The thinner will be the space charge distribution, the stronger will be the velocity gradients.

Interfacial polarity, Z<sub>R</sub>, in the context of free surface flows has been demonstrated to be an important parameter

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Velocity profile with fixed values of De = 0.1, α = 5, γ<sub>R</sub> = 1 and Wo = 10 with (A) different values of interfacial polarity, Z<sub>R</sub> at Θ = π/2 and (B) at different times (Θ) and locations (Y).</figcaption></figure>

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Vorticity distribution in the fluid over one complete time cycle ($\Theta \in (0.2\pi)$). The values of fixed parameters are: $\gamma_R = 1$, $\alpha = 5$, $Z_R = 1$ and (A) $De = 0.1$ and $Wo = 10$; (B) $De = 0.5$ and $Wo = 10$; (C) $De = 0.1$ and $Wo = 50$; (D) $De = 0.5$ and $Wo = 50$.</figcaption></figure>

determining the base state velocity profile presented in Fig. 7. One can see that the nature of velocity profile is a strong function of the interfacial polarity. The presence of strong velocity gradients is a marker of high dispersive fields present near the interface and at the wall. It is seen that the interfacial dispersion can be reduced by reducing the interfacial potential or charge distribution. Such an effect can be achieved by controlling the solvent pH and electrolytic concentrations or using nonelectrolytic surfactants. The strong velocity gradients near the wall can be relaxed by using nonwettable surfaces that enhance velocity slip.

Another aspect of interfacial potential is the phase lag (see Eq. (34)) between velocities at different transverse locations in the flow (see Fig. 7B). It is observed that upon reversing the polarity of the interface induces a phase lag of $\pi$ in the interfacial velocity. However, the near-wall ($Y = 0.1$) and midfilm ($Y = 0.5$) velocities seem to be insensitive to the changes in the interfacial polarity.

Upon observing the vorticity distribution in the fluid over a complete time cycle, one can get an idea of vorticity fluctuation and diffusion in the fluid bulk as a function of time (See Fig. 8). It is observed that the vorticity transport in the bulk is controlled by the Womersley ($Wo$) and Debye ($De$) numbers. As high $Wo$ reduces the rate of viscous diffusion, the fluid bulk is not affected by the time-dependent fluctuations in the vorticity profile even when the EDL is thick. However, for lower $Wo$, the vorticity fluctuation is spread over the entire extent of the space charge distribution.

## 4 Concluding remarks

The present work focuses on the role of time-dependent Maxwell stress in free surface EOFs. It was highlighted that the Maxwell stress-generated dynamics introduces significant velocity gradients at the free surface as compared to the studies where interfacial Maxwell stress was not considered. Moreover, it was observed that for thick EDLs the non-Maxwell stress model significantly underestimates the velocity and vorticity distributions in the fluid film. Such gradients are instrumental in flow instability of time-periodic EOF [48]. The bulk penetration of near-wall oscillations is inversely proportional to the forcing frequency and at lower frequencies we observe a significant suppression of velocity dispersion that asymptotes to the classical dispersion free velocity in a DC EOF. By the variation of interfacial polarity one can not only control the interfacial velocity and the rate of shear transfer but also establish a significant phase difference in fluid velocity at different transverse locations in the fluid, which enhances the dispersion effects and can be effectively used in controlled species transport in microfluidic devices.

*The authors have declared no conflict of interest.*

## 5 References

[1] Stone, H. A., Stroock, A. D., Ajdari, A., *Annu. Rev. Fluid Mech.* 2004, **36**, 381–411.
[2] Sounart, T. L., Baygents, J. C., *J. Fluid Mech.* 2007, 576, 139-172.

[3] Xu, Z., Miao, J., Wang, N., Wen, W., Sheng, P., *Phys. Rev. E* 2011, 83, 066303.

[4] Squires, T. M., Quake, S. R., *Rev. Mod. Phys.* 2005, 77, 977-1026.

[5] Bazant, M., Thornton, K., Ajdari, A., *Phys. Rev. E* 2004, 70, 021506.

[6] Sheng, P., Wen, W., *Annu. Rev. Fluid Mech.* 2012, 44, 143-174.

[7] Lee, J. S. H., Li, D., *Microfluid. Nanofluidics* 2006, 2, 361-365.

[8] Schaeper, J. P., Sepaniak, M. J., *Electrophoresis* 2000, 21, 1421-1429.

[9] Persat, A., Suss, M. E., Santiago, J. G., *Lab Chip* 2009, 9, 2454-2469.

[10] Shin, S. M., Kang, I. S., Cho, Y.-K., *J. Micromech. Microeng.* 2005, 15, 455-462.

[11] Chakraborty, S., Srivastava, A. K., *Langmuir* 2007, 23, 12421-12428.

[12] Chakraborty, S., Ray, S., *Phys. Fluids* 2008, 20, 083602.

[13] Dutta, P., Beskok, A., *Anal. Chem.* 2001, 73, 5097-5102.

[14] Ramos, A., Morgan, H., Green, N. G., Gonzalez, A., Castellanos, A., *J. App. Phys.* 2005, 97, 084906.

[15] Green, N. G., Ramos, A., Gonzalez, A., Morgan, H., Castellanos, A., *Phys. Rev. E* 2000, 61, 4011-4018.

[16] Gonzalez, A., Ramos, A., Green, N., Castellanos, A., Morgan, H., *Phys. Rev. E* 2000, 61, 4019-28.

[17] Ramos, A., Morgan, H., Green, N., Castellanos, A., *J. Colloid Interface Sci.* 1999, 217, 420-422.

[18] Bose, N., Das, T., Chakraborty, D., Maiti, T. K., Chakraborty, S., *Lab Chip* 2012, 12, 69-73.

[19] Chakraborty, J., Ray, S., Chakraborty, S., *Electrophoresis* 2012, 33, 419-425.

[20] Takami, N., Ohsaki, T., Hasebe, H., Yamamoto, M., *J. Electrochem. Soc.* 2002, 149, A9-A12.

[21] Jang, J. H., Yoon, S., Ka, B. H., Jung, Y.-H., Oh, S. M., *J. Electrochem. Soc.* 2005, 152, A1418-A1422.

[22] Kötz, R., *Electrochim. Acta* 2000, 45, 2483-2498.

[23] Gao, Y., Wong, T. N., Yang, C., Ooi, K. T., *J. Colloid Interface Sci.* 2005, 284, 306-314.

[24] Haiwang, L., Wong, T. N., Nguyen, N.-T., *Int. J. Heat Mass Transf.* 2010, 53, 772-785.

[25] Lee, J. S. H., Barbulovic-Nad, I., Wu, Z., Xuan, X., Li, D., *J. App. Phys.* 2006, 99, 054905.

[26] Pascall, A. J., Squires, T. M., *J. Fluid Mech.* 2011, 684, 163-191.

[27] Mayur, M., Amiroudine, S., Lasseux, D., *Phys. Rev. E* 2012, 85, 046301.

[28] Israelachvili, J., *Intermolecular and Surface Forces, Academic Press, San Diego, 2011.*

[29] Cetin, B., Li, D., *Electrophoresis* 2008, 29, 994-1005.

[30] Wei, D., Patey, G. N., *J. Chem. Phys.* 1991, 94, 6795.

[31] Lorenz, C. D., Travesset, A., *Phys. Rev. E* 2007, 75, 061202.

[32] Borukhov, I., Andelman, D., Orland, H., *Europhys. Lett.* 1995, 32, 499-504.

[33] Yossifon, G., Mushenheim, P., Chang, Y.-C., Chang, H., *Phys. Rev. E* 2009, 79, 046305.

[34] Fedorov, M. V., Kornyshev, A. A., *Electrochim. Acta* 2008, 53, 6835-6840.

[35] Dufreche, J., Bernard, O., Turq, P., *J. Mol. Liq.* 2005, 118, 189-194.

[36] Song, J., Kim, M. W., *J. Phys. Chem. B* 2011, 115, 1856-1862.

[37] Kirby, B. J., Hasselbrink, E. F., *Electrophoresis* 2004, 25, 203-213.

[38] Manciu, M., Ruckenstein, E., *Coll. Surf. A* 2012, 400, 27-35.

[39] Li, C., Somasundaran, P., *J. Colloid Interface Sci.* 1991, 146, 215-218.

[40] Gray-Weale, A., Beattie, J. K., *Phys. Chem. Chem. Phys.* 2009, 11, 10994-11005.

[41] Yang, C., Dabros, T., Li, D., Czarnecki, J., Masliyah, J. H., *J. Colloid Interface Sci.* 2001, 243, 128-135.

[42] Graciaa, A., Morel, G., Saulner, P., Lachaise, J., Schechter, R. S., *J. Colloid Interface Sci.* 1995, 172, 131-136.

[43] Choi, W., Sharma, A., Qian, S., Lim, G., Joo, S. W., *J. Colloid Interface Sci.* 2011, 357, 521-526.

[44] Delgado, A. V, Gonzalez-Caballero, F., Hunter, R. J., Koopal, L. K., Lyklema, J., *J. Colloid Interface Sci.* 2007, 309, 194-224.

[45] Cash, J. R., Moore, G., Wright, R. W., *ACM Trans. Math. Soft.* 2001, 27, 245-266.

[46] Lin, C.-H., Ferguson, G. S., Chaudhury, M. K., *Langmuir* 2013, 29, 7793-7801.

[47] Favati, P., Lotti, G., Romani, F., *ACM Trans. Math. Soft.* 1991, 17, 218-232.

[48] Davis, S. H., *Annu. Rev. Fluid Mech.* 1976, 8, 57-74.

Appendix

Amplitude and phase of the velocity field
(see Eq. (34))

From the Eq. (34), the velocity amplitude, $U_a(Y)$, can be
written as:

$$
U_a (Y) = \sqrt{A^2 (Y) + B^2 (Y)} \quad (A.1)
$$

and the phase difference, $\Delta(Y)$, can be written as:

$$
\Delta(Y) = \arctan \left( \frac{A(Y)}{B(Y)} \right) \qquad (A.2)
$$

where:

$$
\begin{equation}
\begin{aligned}
A(Y) ={}& De^2 Wo^2 M_1 \left\{ Z_R \sinh\left(\frac{Y}{De}\right) + \sinh\left(\frac{1-Y}{De}\right) \right\} \\
& + N_1 \sin\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \sinh\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \\
& - N_2 \cos\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \cosh\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \\
& + N_3 \sin\left(\frac{Wo}{\sqrt{2}}Y\right) \cosh\left(\frac{Wo}{\sqrt{2}}Y\right) \\
& - N_4 \cos\left(\frac{Wo}{\sqrt{2}}Y\right) \sinh\left(\frac{Wo}{\sqrt{2}}Y\right),
\end{aligned}
\tag{A.3}
\end{equation}
$$
$$
\begin{align}
B(Y) ={}& M_1 \left\{ Z_R \sinh\left(\frac{Y}{De}\right) \sinh\left(\frac{1-Y}{De}\right) \right\} \nonumber \\
& + N_1 \cos\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \cosh\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \nonumber \\
& + N_2 \sin\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \sinh\left(\frac{Wo}{\sqrt{2}}(1-Y)\right) \nonumber \\
& + N_3 \cos\left(\frac{Wo}{\sqrt{2}}Y\right) \sinh\left(\frac{Wo}{\sqrt{2}}Y\right) \nonumber \\
& + N_4 \sin\left(\frac{Wo}{\sqrt{2}}Y\right) \cosh\left(\frac{Wo}{\sqrt{2}}Y\right), \tag{A.4}
\end{align}
$$

and:

$$
\begin{align}
M_1 &= \frac{\gamma_R}{(De^4 Wo^4 + 1) \sinh\left(\frac{1}{De}\right)}, \tag{A.5} \\
N_1 &= -2M_2 \left( De^2 Wo^2 \sin\left(\frac{Wo}{\sqrt{2}}\right) \sinh\left(\frac{Wo}{\sqrt{2}}\right) \right. \nonumber \\
    &\qquad \left. + \cos\left(\frac{Wo}{\sqrt{2}}\right) \cosh\left(\frac{Wo}{\sqrt{2}}\right) \right), \tag{A.6}
\end{align}
$$

$$
\begin{align}
N_1 = & -2 M_2 \left( De^2 Wo^2 \sin\left(\frac{Wo}{\sqrt{2}}\right) \sinh\left(\frac{Wo}{\sqrt{2}}\right) \right. \nonumber \\
      & \qquad \left. + \cos\left(\frac{Wo}{\sqrt{2}}\right) \cosh\left(\frac{Wo}{\sqrt{2}}\right) \right), \tag{A.6} \\
M_3 = & \frac{\sqrt{2} De Wo \left(Z_R \cosh\left(\frac{1}{De}\right) - 1\right)}{\sinh\left(\frac{1}{De}\right)}, \tag{A.11}
\end{align}
$$