2.1. Equations of motion

This section presents our plane Poiseuille flow of a Bingham fluid with an SH (or CP) lower wall, including the governing continuity and momentum balance equations, in a Cartesian coordinate system $(x,y,z)$ (see figure 1). Motivated by practical considerations, we consider a channel where the lower boundary is a groovy wall and the upper boundary has a no-slip condition; this is because, in practice, the construction of a channel with two symmetrically aligned patterned walls, typically on the micro- or nano-scale, is challenging (Schmieschek et al. Reference Schmieschek, Belyaev, Harting and Vinogradova2012) (we later discuss the extension of our models for the channels with two patterned walls in Appendix F). Considering a creeping flow motion, the dimensionless forms of the continuity and momentum balance equations are

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{U}} = 0, \end{gather}$$

(2.2)$$\begin{gather}- \boldsymbol{\nabla} P + \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\tau }} = 0, \end{gather}$$

where $t$ is the time, ${\boldsymbol {U} = U{\boldsymbol {e}_x} + V{\boldsymbol {e}_y} + W{\boldsymbol {e}_z}}$ is the dimensionless velocity vector, $P$ is the pressure and the deviatoric stress tensor is denoted by ${\boldsymbol {\tau }}$. The pressure gradient is in the $z'$ direction, which is at an angle $\theta$ with the $z$ axis (see figure 1). The dimensionless form of the equations of motion is obtained by considering the half-channel height ($\hat H$) as the characteristic length and the average velocity ($\hat U_{ave}$) as the velocity scale. In addition, the characteristic viscous stress, i.e. $\hat \mu _p ({\hat U_{ave}}/{\hat H})$ where $\hat \mu _p$ is the plastic viscosity, is used to obtain the dimensionless form of the pressure and stress terms.

The Bingham constitutive equation is considered to model the viscoplastic rheology, which in dimensionless form is presented as

(2.3)\begin{equation} \left\{\begin{array}{@{}ll} {\boldsymbol {\tau }} = \left( {1 + \dfrac{B}{{\dot \gamma }}} \right) {\boldsymbol {\dot \gamma }}, & \tau > B,\\ \boldsymbol{\dot \gamma} = 0, & \tau \le B, \end{array}\right. \end{equation}

where the strain rate tensor is shown by ${\boldsymbol {\dot \gamma } = \boldsymbol {\nabla } \boldsymbol {u} + {( {\boldsymbol {\nabla } \boldsymbol {u}} )^{\rm T}}}$, and the norms (magnitudes) of the stress and strain rate tensors are ${\tau = \sqrt {{{{\tau _{ij}}{\tau _{ij}}}/2}}}$ and ${\dot \gamma = \sqrt {{{{{\dot \gamma }_{ij}}{{\dot \gamma }_{ij}}}/2}}}$, respectively (here, $i$ and $j$ refer to the coordinate axes, i.e. $x$, $y$ and $z$). The ratio of the yield stress (${{{\hat \tau _0}}}$) to the characteristic viscous stress is represented by the Bingham number and defined as

(2.4)\begin{equation} B = \frac{{{\hat \tau _0}\hat H}}{{{{\hat \mu }_p}{{\hat U}_{ave}}}}. \end{equation}

Based on (2.3), formation of plug zones are expected for a viscoplastic flow when the applied stress is smaller than the yield stress. For example, formation of an unyielded plug zone in the channel centre is a characteristic of the creeping Poiseuille flow of viscoplastic materials. Based on (2.3), at the yield surfaces and inside the plug, the strain-rate tensor and its norm vanish.

2.2. No-slip base flow

Given the symmetry condition for the no-slip base flow, i.e. with the no-slip condition at both walls, the base flow pressure ($P_0^b$) and velocity ($U^b_0$) for the lower half of the channel can be easily derived as

(2.5)$$\begin{gather} P_0^b ={-} {\tau _w}z' + {\text{constant}}, \end{gather}$$

(2.6)$$\begin{gather}{U^b_0}(\kern0.7pt y) = \left\{\begin{array}{@{}ll} {C_1}y + {C_2}{y^2}, & 0 \le y \le {h_0}, \\ {C_3}, & {h_0} \le y \le 1, \end{array}\right. \end{gather}$$

where $C_1=\tau _w - B$, $C_2=-({\tau _w}/{2})$, $C_3={(\tau _w - B)^2}/{2\tau _w}$ and $h_0$ denotes the location of the lower yield surface. The wall shear stress, $\tau _w$, at $y=0$ corresponds to the largest positive root of the following equation:

(2.7)\begin{equation} 2\tau _w^3 - (3B + 6)\tau _w^2 + {B^3} = 0. \end{equation}

Subsequently, one can find the location of lower yield surface as

(2.8)\begin{equation} {{h_0} = 1 - \frac{B}{{{\tau _w}}}}. \end{equation}

Since the $U^b_0$ vector is in the $z'$ direction, the base no-slip velocity has two components in the $x$ and $z$ directions, i.e. $U_0=U^b_0 \sin \theta$ and $W_0=U^b_0 \cos \theta$, respectively.

2.3. Slip model

The linear Navier slip law is considered to account for the Bingham fluid slippage on the interface. Therefore, the following relation for the slip velocities in the $z$ ($w_s$) and $x$ ($u_s$) directions is derived for the Bingham fluid:

(2.9)\begin{equation} \{ {{w_s},{u_s}}\} = b{\{ {{\tau _{yz}},{\tau _{xy}}}\}_{y = 0}}, \end{equation}

where $b$ represents the dimensionless slip number, defined as $b = {{\hat b {\hat \mu }_p}/{\hat H}}$. Here, $\hat b$ is the dimensional slip number, correlating the dimensional values of the slip velocity and shear stress.

Based on the analyses developed for the Newtonian flows (Schönecker & Hardt Reference Schönecker and Hardt2013; Nizkaya, Asmolov & Vinogradova Reference Nizkaya, Asmolov and Vinogradova2014; Schönecker, Baier & Hardt Reference Schönecker, Baier and Hardt2014), in general, $\hat b$ can be related to the air viscosity ($\hat \mu _a$) and the groove aspect ratio ($A_r=\hat d/\hat \delta$). In addition, $\hat b$ (and, hence, $b$) has been proven to show a tensorial form for SH surfaces, i.e. with different values of the local slip length (number) in the longitudinal and transverse directions. For large values of the local slip length, where the flow condition on the liquid/air interface approaches the no-shear condition, the tensorial form of $\hat b$ (and, hence, $b$) becomes less important and can be effectively replaced by a scalar form, i.e. an equal value can be considered for the slip number in both the longitudinal and transverse directions (Schönecker & Hardt Reference Schönecker and Hardt2013; Nizkaya et al. Reference Nizkaya, Asmolov and Vinogradova2014; Schönecker et al. Reference Schönecker, Baier and Hardt2014). On the other hand, for the CP surfaces, the local slip length can be identical in the longitudinal and transverse conditions, rendering a scalar form of the slip number (Qian et al. Reference Qian, Wang and Sheng2005; Wang et al. Reference Wang, Qian and Sheng2008; Lee et al. Reference Lee, Charrault and Neto2014). Since the tensorial form of the slip number has not yet been explored for the viscoplastic flows (refer to the discussion in § 7), here we assume a scalar form of the slip number, i.e. a similar value of $\hat b$ (and, hence, $b$) for longitudinal, transverse and oblique configurations. However, assuming a tensorial slip number, i.e. any functionality of $\hat b$ (and, hence, $b$) with respect to the angle $\theta$, the models developed in this study would remain valid (see preliminary tensorial form analysis of the local slip number in Appendix D). We may expect that the assumed scalar form for the slip number can provide us with reasonably accurate predictions of the flow dynamics over the CP surfaces, while capturing leading physical trends for the flow over SH surfaces at large values of the slip number. In addition, the upcoming solutions may provide some insight about the viscoplastic flow dynamics over the liquid-infused (LI) surfaces; however, to provide accurate results for this case, the exact tensorial form of the slip number must be employed and the corresponding scalar form is not valid (Schönecker et al. Reference Schönecker, Baier and Hardt2014).

Before proceeding, it is worth mentioning that both the groove (stripe) period (width) (Sbragaglia & Prosperetti Reference Sbragaglia and Prosperetti2007; Hodes et al. Reference Hodes, Kirk, Karamanis and MacLachlan2017; Game, Hodes & Papageorgiou Reference Game, Hodes and Papageorgiou2019) and channel height (or pipe diameter) (Lauga & Stone Reference Lauga and Stone2003; Schnitzer & Yariv Reference Schnitzer and Yariv2017, Reference Schnitzer and Yariv2019; Kirk et al. Reference Kirk, Karamanis, Crowdy and Hodes2020) have been used within the literature as the characteristic lengths to make the slip length dimensionless. In this work, for definition of the slip number, i.e. $b={\hat b \hat \mu _p}/{\hat H}$, the half-channel height ($\hat H$) is used as the characteristic length. Although employing the groove (stripe) period ($\hat L$) as the characteristic length can be useful when interpreting the local slip dynamics on the patterned wall, our definition of the slip number is also physically relevant, as it provides an understanding regarding the overall patterned wall effects on the channel flow dynamics. In addition, the defined dimensionless slip number in our work can be simply converted to the slip number that is made dimensionless using the groove (stripe) period ($b_L$) as $b_L=b/\ell$. Therefore, one can also use such a conversion to be able to interpret the results, if required. In addition, the usage of $\hat H$ as the characteristic length is prevalent and advantageous in defining the Bingham number ($B$), which is a key parameter of our viscoplastic flow.

3.1. Perturbation equations

In this work, semi-analytical and explicit-form solutions are developed for the Poiseuille flow of Bingham fluids in channels with a patterned wall, considering the limiting cases of longitudinal ($\theta =0$) and transverse ($\theta =90^\circ$) grooves (stripes), as well as the general case of the oblique ($0 < \theta < 90^\circ$) flow configuration. These solutions are developed for the thick channel limit ($\ell ={\hat {L}}/{\hat {H}} \ll 1$), where the lower yield surface (located at $y=h$) remains flat (Rahmani & Taghavi Reference Rahmani and Taghavi2022). The solution can be derived by considering infinitesimal perturbations induced by the patterned wall (with respect to the no-slip flow) in the lower yielded zone ($-\ell /2 \le x \le \ell /2$ and $0 \le y < h$) and solving for the leading-order terms. Therefore, the total velocity vector ($\boldsymbol {U}$) is written as

(3.1)\begin{equation} \boldsymbol{U}=\boldsymbol{U}_0+\epsilon \boldsymbol{u}, \end{equation}

where $\boldsymbol {U}_0 = U_0{\boldsymbol {e}_x} + W_0{\boldsymbol {e}_z}$, i.e. $(U_0,0,W_0)$, represents the no-slip velocity vector and

(3.2)\begin{equation} \boldsymbol{u} = u{\boldsymbol {e}_x} + v{\boldsymbol {e}_y} + w{\boldsymbol {e}_z} \end{equation}

shows the perturbation velocity vector $(u,v,w)$. Here, $\epsilon =\kappa ^{-1}$ is the perturbation parameter where $\kappa =2{\rm \pi} /\ell$ is the patterned wall wavenumber. Based on this perturbation method, the effective viscosity of the Bingham fluid ($\mu$) and the stress components are expanded as follows:

(3.3)$$\begin{gather} \mu (\boldsymbol{U_0} + \epsilon \boldsymbol{u}) = \mu (\boldsymbol{U_0}) + \epsilon {\dot \xi _{ij}}(\boldsymbol{u})\frac{{\partial \mu }}{{\partial {{\dot \gamma }_{ij}}}} (\boldsymbol{U_0}) + O({\epsilon ^2}), \end{gather}$$

(3.4)$$\begin{gather}{\tau _{ij}} (\boldsymbol{U_0} + \epsilon \boldsymbol{u}) = {\tau _{ij}}(\boldsymbol{U_0}) + \epsilon {\sigma _{ij}}(\boldsymbol{u}) + O({\epsilon ^2}), \end{gather}$$

where $\dot \xi _{ij}$ and $\sigma _{ij}$ represent the component of the perturbation strain-rate and stress tensor, respectively. Henceforth, for presentation simplicity, we use $\mu _0=\mu (\boldsymbol {U_0})$.

3.1.1. Longitudinal configuration

For the longitudinal groove (stripe) configuration, the no-slip and perturbation velocity vectors reduce to ($0,0,W_0$) and ($0,0,w$), respectively. In addition, considering an infinitely long channel in both x and z directions, the velocity gradients in the $z$ direction are zero ($\partial \boldsymbol{U} /\partial z = 0$, a feature valid for any angle of $\theta$), only two non-zero perturbation stress components remain for the Bingham fluid as

(3.5)\begin{equation} \left.\begin{gathered} {\sigma _{yz}} = \frac{{\partial w}}{{\partial y}}, \\ {\sigma _{xz}} = \mu_0 \frac{{\partial w}}{{\partial x}}, \end{gathered}\right\} \end{equation}

where $\mu _0=1+{B}/{({\rm d} U^b_0/{{\rm d}y})}$. Therefore, the perturbed momentum balance equation shrinks to

(3.6)\begin{equation} \mu_0 \frac{{{\partial ^2}w}}{{\partial {x^2}}} + \frac{{{\partial ^2}w}}{{\partial {y^2}}} = 0, \end{equation}

in which the perturbation parameter is dropped.

In the present study, to facilitate the identification of various orders of perturbed momentum balance equations for a given flow configuration, a rescaling is introduced as

(3.7a–c)\begin{equation} \epsilon X=x,\quad \epsilon Y=y,\quad \epsilon \varPsi=\psi, \end{equation}

where $\psi$ is the perturbation stream function such that $u = {{\partial \psi }}/{{\partial y}}$ and $v = -{{\partial \psi }}/{{\partial x}}$. The introduced rescaling is motivated by the decay of perturbations in the $y$ direction, which renders the perturbation negligible beyond a distance comparable to the groove periodicity length. Consequently, the rescaled distance over which the perturbation field is considerable becomes of order $\epsilon ^0$. Thus, we may expand the viscosity over $Y=0$ to reach

(3.8)\begin{equation} \mu_0=1+\frac{B}{C_1}-\epsilon \left(\frac{2BC_2}{C_1^2}\right) Y, \end{equation}

in which coefficients $C_1$ and $C_2$ are given in § 2.2. Therefore, the rescaled form of (3.6) for the leading-order terms is obtained as

(3.9)\begin{equation} \left(1+\frac{B}{C_1}\right) \frac{{{\partial ^2}w}}{{\partial {X^2}}} + \frac{{{\partial ^2}w}}{{\partial {Y^2}}} = 0, \end{equation}

where $\epsilon$ drops from (3.9).

Following the convention in the field (Lauga & Stone Reference Lauga and Stone2003; Belyaev & Vinogradova Reference Belyaev and Vinogradova2010; Vinogradova & Belyaev Reference Vinogradova and Belyaev2011), we consider our creeping flow to be periodic in the $x$ direction with a period of $\ell$. Accordingly, in the rescaled system, the flow is periodic in the $X$ direction with a period of $2{\rm \pi}$. Therefore, the solution for $w$ can be written in the Fourier series form as

(3.10)\begin{equation} w (X,Y) = \sum_{n = 0 }^\infty {B_n}{\hat w_n } (Y){\cos ({nX})}, \end{equation}

where $B_n$ are unknown coefficients that will be calculated later with the use of appropriate patterned wall conditions. Here, the hat sign is used to define the Fourier coefficient (i.e. $\hat w_n$).

Substituting (3.10) into (3.9), one can find the following ordinary differential equation (ODE):

(3.11)\begin{equation} \frac{{{{\rm d}^2}\hat w_n}}{{{\rm d} {Y^2}}} - \left(1+\frac{B}{C_1}\right){n^2}\hat w_n = 0. \end{equation}

To close the boundary value problem, the following boundary conditions are considered ($n \ne 0$):

(3.12a,b)\begin{equation} \hat w_n(\kappa h) = 0,\quad \frac{{{\rm d} \hat w_n}}{{{\rm d} Y}}(\kappa h) = 0, \end{equation}

where the first and second conditions from left to right ensure zero perturbation velocity and zero strain-rate magnitude at the lower yield surface ($Y=\kappa h$), respectively. The zeroth term ($n=0$) has a linear distribution in $Y$ and is discussed later in § 3.2.

3.1.2. Transverse configuration

In the transverse configuration, the no-slip and perturbation velocity vectors reduce to $(U_0,0,0)$ and $(u,v,0)$, respectively. Perturbing the momentum balance equation, eliminating the pressure, using the definition of $\varPsi$, considering the rescaling, and expanding $\mu _0$, the following partial differential equation is eventually derived at the leading order:

(3.13)\begin{equation} \frac{{{\partial ^4}\varPsi }}{{\partial {X^4}}} + \frac{{{\partial ^4}\varPsi }}{{\partial {Y^4}}} + \left( {2 + \frac{{4B}}{{{C_1}}}} \right)\frac{{{\partial ^4}\varPsi }}{{\partial {X^2}\partial {Y^2}}} = 0, \end{equation}

where $\epsilon$ drops from (3.13). Considering the flow periodicity in $X$, one can write

(3.14)\begin{equation} \varPsi (X,Y) = \sum_{n = 0 }^\infty {A_n}{\hat \varPsi_n } (Y){\cos ({nX})}, \end{equation}

where $A_n$ are unknown coefficients and will be obtained later using the patterned wall condition. Here, the hat sign is used to define the Fourier coefficient (i.e. $\hat \varPsi _n$).

Substituting (3.14) into (3.13), we find the following ODE:

(3.15)\begin{equation} \frac{{{{\rm d}^4}\hat \varPsi_n }}{{{\rm d} {Y^4}}} - \left( {2 + \frac{{4B}}{{{C_1}}}} \right){n^2} \frac{{{{\rm d}^2}\hat \varPsi_n }}{{{\rm d} {Y^2}}} + {n^4}\hat \varPsi_n = 0. \end{equation}

To close the boundary value problem, we consider the following conditions ($n \ne 0$)

(3.16a–d)\begin{equation} \hat \varPsi_n (0) = 0,\quad \hat \varPsi_n (\kappa h) = 0,\quad \frac{{{\rm d} \hat \varPsi_n }}{{{\rm d} Y}}(\kappa h) = 0,\quad \frac{{{{\rm d}^2}\hat \varPsi_n }}{{{\rm d} {Y^2}}}(\kappa h) = 0, \end{equation}

where the first condition from left to right ensures the no-penetration assumption at the patterned wall (i.e. $V=0$ at $Y=0$). The zero perturbation velocity field at the lower yield surface ($Y=\kappa h$) is guaranteed using the second and third conditions (from left to right). Finally, the zero strain-rate magnitude at the lower yield surface is obtained using the fourth condition. Here, the zeroth term ($n=0$) would show a quadratic functionality of $Y$, as discussed in § 3.2.

3.1.3. Oblique configuration

In the oblique configuration, the no-slip and perturbation velocity vectors, i.e. $(U_0,0,W_0)$ and $(u,v,w)$, respectively, have their most general forms. Accordingly, the perturbation stress components can be obtained as given in Appendix A. Then, using the definition of $\varPsi$, adopting the rescaling introduced in (3.7a–c), expanding $\mu _0$ via (3.8) and performing considerable algebra, we derive a set of two coupled partial differential equations for $\varPsi$ and $w$ at the leading order:

(3.17) \begin{gather} {\left(1+\frac{B}{C_1}\right)}\left( {\frac{{{\partial ^4}\varPsi }}{{\partial {X^4}}} + \frac{{{\partial ^4}\varPsi }}{{\partial {Y^4}}} + 2\frac{{{\partial ^4}\varPsi }}{{\partial {X^2}\partial {Y^2}}}} \right) - {\frac{B}{C_1}}\left( {\frac{{{\partial ^4}\varPsi }}{{\partial {X^4}}} + \frac{{{\partial ^4}\varPsi }}{{\partial {Y^4}}} - 2\frac{{{\partial ^4}\varPsi }}{{\partial {X^2}\partial {Y^2}}}} \right){\sin ^2}\theta \nonumber\\ -{\frac{B}{C_1}}\left( {\frac{{{\partial ^3}w}}{{\partial {Y^3}}} - \frac{{{\partial ^3}w}}{{\partial {X^2}\partial Y}}} \right) \sin \theta \cos \theta = 0, \end{gather}

(3.18)\begin{gather} {\left(1+\frac{B}{C_1}\right)}\left( {\frac{{{\partial ^2}w}}{{\partial {X^2}}} + \frac{{{\partial^2}w}}{{\partial {Y^2}}}} \right) -{\frac{B}{C_1}} \frac{{{\partial ^2}w}}{{\partial {Y^2}}}{\cos^2}\theta -{\frac{B}{C_1}}\left( {\frac{{{\partial ^3}\varPsi }}{{\partial {Y^3}}} - \frac{{{\partial ^3}\varPsi }}{{\partial {X^2}\partial Y}}}\right) \sin \theta \cos \theta = 0, \end{gather}

where $\epsilon$ is dropped from the above equations. Equations (3.17) and (3.18) are unique to viscoplastic Bingham materials, significantly differing from the corresponding equations for Newtonian fluids. In the case of Newtonian fluids, we have $B=0$ (i.e. $\mu _0=1$; see (3.8)); consequently, (3.17) and (3.18) decouple, and become independent of $\theta$, i.e. reducing to the corresponding equations for the transverse and longitudinal flows, respectively. Thus, for Newtonian fluids, it suffices to solve the perturbed equation only for the longitudinal and transverse configurations, and the solution for the oblique configuration can be simply obtained via a linear vector transformation through the transform (or rotation) matrix (Bazant & Vinogradova Reference Bazant and Vinogradova2008; Vinogradova & Belyaev Reference Vinogradova and Belyaev2011). Therefore, it is clear that the nonlinear viscoplastic rheology sets our work apart from the Newtonian case, and it plays a crucial role in determining the slip dynamics on a patterned wall, i.e. a feature that we address in this study.

Considering the flow periodicity in the $X$ direction, one can substitute (3.10) and (3.14) into (3.17) and (3.18) to arrive at the following coupled ODEs:

(3.19)\begin{align} & {A_n}\left[ {{\left(1+\frac{B}{C_1}\right)}\left({\frac{{{{\rm d}^4}\hat \varPsi_n }}{{{\rm d} {Y^4}}} - 2{n^2} \frac{{{{\rm d}^2}\hat \varPsi_n }}{{{\rm d} {Y^2}}} + {n^4}\hat \varPsi_n}\right) -{\frac{B}{C_1}}\left({\frac{{{{\rm d}^4}\hat \varPsi_n }}{{{\rm d} {Y^4}}} + 2{n^2} \frac{{{{\rm d}^2}\hat \varPsi_n }}{{{\rm d} {Y^2}}} + {n^4}\hat \varPsi_n } \right){{\sin }^2}\theta } \right] \nonumber\\ &\quad -{B_n}{\frac{B}{C_1}}\left( {\frac{{{{\rm d}^3}\hat w_n}}{{{\rm d} {Y^3}}} + {n^2}\frac{{{\rm d} \hat w_n}}{{{\rm d} Y}}} \right)\sin \theta \cos \theta = 0, \end{align}

(3.20)\begin{align} & -{A_n}{\frac{B}{C_1}}\left( {\frac{{{{\rm d}^3}\hat \varPsi_n }}{{{\rm d} {Y^3}}} + {n^2}\frac{{{\rm d} \hat \varPsi_n }}{{{\rm d} Y}}} \right)\sin \theta \cos \theta \nonumber\\ &\quad + {B_n}{\left[ {{\left(1+\frac{B}{C_1}\right)}\left( {\frac{{{{\rm d}^2}\hat w_n}}{{{\rm d} {Y^2}}} - {n^2}\hat w_n} \right) -{\frac{B}{C_1}}\frac{{{{\rm d}^2}\hat w_n}}{{{\rm d} {Y^2}}}{{\cos }^2}\theta } \right] = 0.} \end{align}

The boundary conditions for these equations are the combination of those for the longitudinal and transverse configurations, i.e. (3.12a,b) and (3.16a–d). In addition, the zeroth term solutions ($n=0$) are treated similar to those of the longitudinal and transverse configurations and are discussed further in § 3.2.

3.2. Semi-analytical solution

In this section, we attempt to solve the ODEs (3.11), (3.15), (3.19) and (3.20). These ODEs have constant coefficients; thus, the solutions for $\hat \varPsi _n$ and $\hat w_n$ would be in the form ${\rm e}^{\varLambda _n Y}$. In the thick channel limit, since the perturbation field quickly decays in the $Y$ direction, the contribution of the terms having positive ${\varLambda _n}$, i.e. $\varLambda _n>0$, becomes negligible. In other words, one can simply show that the terms with positive ${\varLambda _n}$ belong to the higher orders of perturbations. Therefore, in the following sections, the solutions for $\hat \varPsi _n$ and $\hat w_n$ are obtained considering only the terms with negative ${\varLambda _n}$.

Before proceeding, let us consider the perturbation velocity field as $\boldsymbol u$ instead of $\epsilon \boldsymbol u$ in the following sections. This consideration leads to having simpler forms of equations, as the redundant $\epsilon$ symbol disappears.

3.2.1. Longitudinal configuration

For the longitudinal configuration, (3.11) is solved analytically, keeping the term with the negative ${\varLambda _n}$, as

(3.21a,b)\begin{equation} {\hat w_n}(Y) = {\exp({ - {\lambda ^\parallel } nY})},\quad {\lambda ^\parallel } = \sqrt {1 + \frac{B}{C_1}}, \end{equation}

where the symbol $\parallel$ represents the longitudinal flow, henceforth. The expression obtained for $\hat w_n$ would satisfy the conditions of (3.12a,b), since at $Y=\kappa h$, i.e. the lower yield surface, $\hat w_n$ and ${\rm d}\hat w_n/{\rm d} Y$ become negligible and belong to the higher orders of perturbation. After solving for $\hat w_n$, the solution for $w$ can be written in the following form:

(3.22)\begin{equation} w(X,Y) = {B_0}\left( {1 - \frac{Y}{\kappa h}} \right) + \sum_{n = 1}^\infty {{B_n}} \hat w_n (Y)\cos (nX), \end{equation}

where the zeroth term solution has a linear distribution in $Y$ and satisfies the zero perturbation velocity condition at the lower yield surface, i.e. $Y=\kappa h$. The remaining condition to be satisfied by the zeroth term solution is the zero strain-rate magnitude at the lower yield surface, which is satisfied through numerical iterations on $h$, as explained in § 4.1.1.

For the longitudinal configuration, the boundary conditions at the patterned wall is obtained as

(3.23)$$\begin{gather} {w}(X,0) - b{\left( {B + C_1 + \kappa \frac{{\partial w}}{{\partial Y}}} \right)_{Y = 0}} = 0,\quad 0 \leq X \leq {\rm \pi}\varphi, \end{gather}$$

(3.24)$$\begin{gather}{w}(X,0) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}, \end{gather}$$

where $\varphi ={\hat \delta }/{\hat L}$ is called the slip area fraction and represents the fraction of the patterned wall that undergoes the slip condition. Equations (3.23) and (3.24) together represent the patterned wall condition that has been already discussed.

After substituting the solution for the perturbation velocity field, i.e. (3.22), into the patterned wall condition, i.e. (3.23) and (3.24), the following dual trigonometric series problem emerges:

(3.25)$$\begin{gather} {B_0}\left( {1 + \frac{b}{h}} \right) + \sum_{n = 1}^\infty {{B_n}} \left[ {{{\hat w}_n}(0) - b\kappa{\frac{{\rm d}\hat w_n}{{\rm d} Y}(0)}} \right]\cos ({nX}) = b({B + {C_1}}),\quad 0 \leq X \leq {\rm \pi}\varphi, \end{gather}$$

(3.26)$$\begin{gather}{B_0} + \sum_{n = 1}^\infty {{B_n}} {{\hat w}_n}(0)\cos ({nX}) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}. \end{gather}$$

To calculate the unknown coefficients $B_n$ ($n=0,1,2,\ldots$), we first integrate (3.25) in $[0\ X]$ and then multiply it by ${\sin }(mX)$, where $m$ is a positive integer. Afterwards, we integrate the resulting equation in $[0\ {\rm \pi}\varphi ]$. Second, we multiply (3.26) by ${\cos }(mX)$ and then integrate it in $[{\rm \pi} \varphi \ {\rm \pi}]$. After the summation of the obtained terms, we establish a system of linear equations (truncated at the $N$th term):

(3.27)\begin{equation} \sum_{n = 0}^N {{P_{mn}^{{\parallel}}}{B_n}} = {M_m^{{\parallel}}}, \end{equation}

where the above coefficient and constant matrices are obtained as

(3.28)$$\begin{gather} {P_{m0}^{{\parallel}}} = \left( {1 + \frac{b}{h}} \right) I_3+ I_4, \end{gather}$$

(3.29)$$\begin{gather}{P_{mn}^{{\parallel}}} = \frac{1}{{n}}\left[ {{{\hat w}_n}(0) - b\kappa{\frac{{\rm d}\hat w_n}{{\rm d} Y}(0)}} \right] I_1 + {{\hat w}_n}(0) I_2,\quad n>0, \end{gather}$$

(3.30)$$\begin{gather}{M_m^{{\parallel}}} = b\left( {B + {C_1}} \right) I_3, \end{gather}$$

where $I_1$, $I_2$, $I_3$ and $I_4$ are the following integrals

(3.31)$$\begin{gather} I_1 = \int_0^{{\rm \pi} \varphi } {\sin } (nX)\sin (mX)\,{\rm d} X, \end{gather}$$

(3.32)$$\begin{gather}I_2 = \int_{{\rm \pi} \varphi }^{\rm \pi} {\cos } (nX)\cos (mX)\,{\rm d} X, \end{gather}$$

(3.33)$$\begin{gather}I_3 = \int_0^{{\rm \pi} \varphi } X \sin (mX)\,{\rm d} X, \end{gather}$$

(3.34)$$\begin{gather}I_4 = \int_{{\rm \pi} \varphi }^{\rm \pi} \cos (mX)\,{\rm d} X. \end{gather}$$

3.2.2. Transverse configuration

Solving the ODE (3.15) while keeping the terms with negative ${\varLambda _n}$, we find

(3.35a,b)\begin{align} {\hat \varPsi _n}\left( Y \right) = {{\exp({ - \lambda _1^ \bot nY})} - {\exp({ - \lambda _2^ \bot nY})}},\quad \left\{\begin{array}{@{}l} \lambda _1^ \bot = \dfrac{{\sqrt {{C_1}({2B + {C_1} + 2\sqrt {{B^2} + B{C_1}}})} }}{{{C_1}}},\\ \lambda _2^ \bot = \dfrac{{\sqrt {{C_1}({2B + {C_1} - 2\sqrt {{B^2} + B{C_1}}})} }}{{{C_1}}}, \end{array}\right. \end{align}

where the symbol $\bot$ represents the transverse flow henceforth. The solution form developed previously for $\hat \varPsi _n$ satisfies the no-penetration condition at the patterned wall, i.e. $\hat \varPsi _n(0)=0$, hence $v(X,0)=0$ (see (3.38) further below). At $Y=\kappa h$, the values of $\hat \varPsi _n$, ${\rm d}\hat \varPsi _n/{\rm d} Y$ and ${\rm d}^2\hat \varPsi _n/{\rm d} Y^2$ become negligible and they belong to higher orders of perturbations; thus, the conditions of zero perturbation velocity and zero strain-rate magnitude at the lower yield surface are satisfied (for $n>0$ in the leading order). The zeroth term solution (i.e. $n=0$) would have a quadratic form for $\varPsi$, leading to a linear distribution for $u$ (similar to $w$). Again, the no-penetration condition at $Y=0$ and the zero perturbation velocity at the lower yield surface ($Y=\kappa h$) is satisfied by the zeroth term solution. However, the condition of zero strain rate magnitude at the yield surface should be satisfied through an iterative approach on $h$, as described in § 4.1.1.

Having $\hat \varPsi _n$, one can now write the solution for $\varPsi$, $u$ and $v$ as

(3.36)$$\begin{gather} \varPsi (X,Y) = {A_0} \left( {Y - \frac{{{Y^2}}}{{2\kappa h}}} \right) + \sum_{n = 1}^\infty {{A_n}} \hat \varPsi_n (Y)\cos ({nX}), \end{gather}$$

(3.37)$$\begin{gather}u(X,Y) = {A_0} \left({1 - \frac{Y}{\kappa h}} \right) + \sum_{n = 1}^\infty {{A_n}} {\frac{{\rm d}\hat \varPsi_n}{{\rm d} Y}(Y)} \cos ({nX}), \end{gather}$$

(3.38)$$\begin{gather}v(X,Y) = \sum_{n = 1}^\infty {{A_n}} \hat \varPsi_n (Y) n \sin ({nX}). \end{gather}$$

The patterned wall condition for the transverse configuration holds

(3.39)$$\begin{gather} {u}(X,0) - b{\left({ B + C_1 + \kappa\frac{{\partial u}}{{\partial Y}}} \right)_{Y = 0}} = 0,\quad 0 \leq X \leq {\rm \pi}\varphi, \end{gather}$$

(3.40)$$\begin{gather}{u}(X,0) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}. \end{gather}$$

Substituting (3.37) into (3.39) and (3.40) leads to a dual trigonometric series problem:

(3.41)$$\begin{gather} {A_0}\left( {1 + \frac{b}{h}} \right) + \sum_{n = 1}^\infty {{A_n}} \left[ {{\frac{{\rm d}\hat \varPsi_n}{{\rm d} Y}(0)} - b\kappa{\frac{{\rm d}^2\hat \varPsi_n}{{\rm d} Y^2}(0)}} \right]\cos ({nX}) = b({B + {C_1}}),\quad 0 \leq X \leq {\rm \pi}\varphi, \end{gather}$$

(3.42)$$\begin{gather}{A_0} + \sum_{n = 1}^\infty {{A_n}} {\frac{{\rm d}\hat \varPsi_n}{{\rm d} Y}(0)}\cos ({nX}) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}, \end{gather}$$

where (3.41) and (3.41) can be solved for $A_n$, using the same method described for the longitudinal flow configuration. Thus, the following system is obtained

(3.43)\begin{equation} \sum_{n = 0}^N {{P_{mn}^{\bot}}{A_n}} = {M_m^{\bot}}, \end{equation}

where the coefficient and constant matrices for (3.43) are calculated as

(3.44)$$\begin{gather} {P_{m0}^{\bot}} = {P_{m0}^{{\parallel}}}, \end{gather}$$

(3.45)$$\begin{gather}{P_{mn}^{\bot}} = \frac{1}{{n}}\left[ {{\frac{{\rm d}\hat \varPsi_n}{{\rm d} Y}(0)} - b\kappa{\frac{{\rm d}^2\hat \varPsi_n}{{\rm d} Y^2}(0)}} \right] I_1 + {\frac{{\rm d}\hat \varPsi_n}{{\rm d} Y}(0)} I_2,\quad n>0, \end{gather}$$

(3.46)$$\begin{gather}{M_m^{\bot}} = {M_m^{{\parallel}}}. \end{gather}$$

3.2.3. Oblique configuration

Considering the solution for $\hat w_n$ and $\hat \varPsi _n$ to be in the form of ${\rm e}^{\varLambda _nY}$, the ODEs (3.19) and (3.20) lead to the following relation:

(3.47)\begin{equation} \left(\begin{array}{@{}l@{}} {\varOmega_{11}}\quad {\varOmega_{12}}\\ {\varOmega_{21}}\quad {\varOmega_{22}} \end{array} \right)\left(\begin{array}{@{}l@{}} {A_n}\\ {B_n} \end{array}\right) = 0, \end{equation}

where

(3.48)$$\begin{gather} \varOmega_{11} = {\left(1+\frac{B}{C_1}\right)} (\varLambda_n^2 - n^2)^2 - \frac{B}{C_1} (\varLambda_n^2 + n^2)^2 \sin^2 \theta, \end{gather}$$

(3.49)$$\begin{gather}\varOmega_{22} = {\left(1+\frac{B}{C_1}\right)} (\varLambda_n^2 - n^2) - \frac{B}{C_1} \varLambda_n^2 \cos^2 \theta, \end{gather}$$

(3.50)$$\begin{gather}\varOmega_{12} ={-}\frac{B}{C_1} (\varLambda_n^2 + n^2) \varLambda_n \sin \theta \cos \theta, \end{gather}$$

(3.51)$$\begin{gather}\varOmega_{21}=\varOmega_{12}. \end{gather}$$

In order to have non-zero solution for $A_n$ and $B_n$, the determinant of the coefficient matrix in (3.47) should be zero; thus, $\varOmega _{11}\varOmega _{22}-\varOmega _{12}\varOmega _{21}=0$. This condition leads to finding six solutions for ${\varLambda _n}$, with three being negative. Keeping the terms with negative ${\varLambda _n}$, the solution for $\hat w_n$ and $\hat \varPsi _n$ can be written as (for $n>0$)

(3.52)$$\begin{gather} {\hat w_n}(Y) = {\exp({ - {\lambda_1 ^\angle } nY})}+\varGamma_n^{(1)} {\exp({ - {\lambda_2 ^\angle } nY})}+\varGamma_n^{(2)} {{\rm exp}(-nY)}, \end{gather}$$

(3.53)$$\begin{gather}{\hat \varPsi _n}(Y) = \varDelta^{(1)}_n{{\exp({ - \lambda _1^ \angle nY})} - {\exp({ - \lambda _2^ \angle nY})}}+\varDelta^{(2)}_n {{\rm exp}(-nY)}, \end{gather}$$

where the symbol $\angle$ represents the oblique flow, henceforth, and

(3.54) $$\begin{gather} {\lambda _1^\angle }{ = \frac{{\sqrt {{C_1}\left( {2B \,{+}\, {C_1} \,{-}\, \dfrac{3}{2}B{{\cos }^2}\theta + 2\sqrt {{B^2} + B{C_1} + \dfrac{9}{{16}}{B^2}{{\cos }^4}\theta \,{-}\, B{{\cos }^2}\theta \left( {\dfrac{3}{2}B \,{+}\, {C_1}} \right)} } \,\right)} }}{{{C_1}}},} \end{gather}$$

(3.55) $$\begin{gather}{\lambda _2^\angle }{ = \frac{{\sqrt {{C_1}\left( {2B \,{+}\, {C_1} \,{-}\, \dfrac{3}{2}B{{\cos }^2}\theta - 2\sqrt {{B^2} + B{C_1} + \dfrac{9}{{16}}{B^2}{{\cos }^4}\theta \,{-}\, B{{\cos }^2}\theta \left({\dfrac{3}{2}B \,{+}\, {C_1}} \right)} } \,\right)} }}{{{C_1}}}.} \end{gather}$$

The solution for each $\hat w_n$ and $\hat \varPsi _n$ contains three terms; thus, for each solution, two coefficients are required in order to determine the contribution of each term to that solution, i.e. $\varGamma _n^{(1)}$ and $\varGamma _n^{(2)}$ for $\hat w_n$ and $\varDelta _n^{(1)}$ and $\varDelta _n^{(2)}$ for $\hat \varPsi _n$. Since the no-penetration condition should be satisfied at $Y=0$, one finds $\varDelta _n^{(2)}=1-\varDelta _n^{(1)}$. It is worth mentioning that the forms of solution written in (3.52) and (3.53) are comparable with those of the longitudinal and transverse configurations (3.21a,b) and (3.35a,b), as one might expect that when $\theta \to 0$, $\varGamma _n^{(1)} \to 0$ and $\varGamma _n^{(2)} \to 0$; however, when $\theta \to 90^\circ$, $\varDelta _n^{(1)} \to 1$ and, hence, $\varDelta _n^{(2)} \to 0$. One should also note that when $\theta \to 0$, $\lambda _1^\angle \to \lambda ^\parallel$ and once $\theta \to 90^\circ$, $\lambda _1^\angle \to \lambda _1^\bot$ and $\lambda _2^\angle \to \lambda _2^\bot$.

Substituting the solution terms for $\hat w_n$ and $\hat \varPsi _n$ with identical ${\varLambda _n}$, e.g. ${\exp ({- {\lambda _1 ^\angle } nY})}$ and $\varDelta ^{(1)}_n{\exp ({ - \lambda _1^ \angle nY})}$, into (3.20), one can obtain following relations:

(3.56)$$\begin{gather} {\varDelta^{(1)} _n} ={-} \left( {\frac{{{B_n}}}{{{A_n}}}} \right){\left( {\frac{{{\varOmega _{22}}}}{{{\varOmega _{12}}}}} \right)_{\varLambda_n ={-}\lambda _1^\angle n }}, \end{gather}$$

(3.57)$$\begin{gather}\varGamma _n^{(1)} = \left( {\frac{{{A_n}}}{{{B_n}}}} \right){\left( {\frac{{{\varOmega _{12}}}}{{{\varOmega _{22}}}}} \right)_{\varLambda_n ={-}\lambda _2^\angle n }}, \end{gather}$$

(3.58)$$\begin{gather}\varGamma _n^{(2)} = ({{\varDelta^{(1)} _n} - 1})\left( {\frac{{{A_n}}}{{{B_n}}}} \right) {\left({\frac{{{\varOmega _{12}}}}{{{\varOmega _{22}}}}} \right)_{\varLambda_n ={-} n}}. \end{gather}$$

The numerical procedure for calculating the above-mentioned coefficients is discussed at the end of this section (i.e. § 3.2.3).

The solutions provided in (3.52) and (3.53) are for $n>0$. Indeed, the solution form for $n=0$ of the oblique flow is identical to those of the longitudinal and transverse flow (i.e. the zero mode of the perturbation velocity is a linear function of $Y$). In fact, when considering $n=0$ in (3.52) and (3.53), the exponential forms vanish and the solutions for $n=0$ become similar to those obtained for the longitudinal and transverse flows (see (3.22) for $w_0$ and (3.36) for $\varPsi _0$).

In the oblique configuration, the patterned wall conditions are as follows:

(3.59) $$\begin{gather} w(X,0) - b{\left[ \begin{array}{@{}l@{}} ({B + C_1})\cos \theta \\ + \left( {1 + \dfrac{B}{{{C_1}}}{{\sin}^2}\theta } \right){\kappa }\dfrac{{\partial w}}{{\partial Y}}\\ - \dfrac{B}{C_1} \sin \theta \cos \theta {\kappa }\dfrac{{\partial u}}{{\partial Y}} \end{array}\right]_{Y = 0}} = 0,\quad 0 \le X \le {\rm \pi} \varphi, \end{gather}$$

(3.60) $$\begin{gather}{w}(X,0) = 0,\quad {\rm \pi} \varphi \leq X \leq {\rm \pi}. \end{gather}$$

(3.61) $$\begin{gather}u(X,0) - b{\left[ \begin{array}{@{}l@{}} ({B + C_1})\sin \theta \\ + \left( {1 + \dfrac{B}{{{C_1}}}{{\cos}^2}\theta } \right){\kappa } \dfrac{{\partial u}}{{\partial Y}}\\ - \dfrac{B}{C_1}\sin \theta \cos \theta {\kappa } \dfrac{{\partial w}}{{\partial Y}} \end{array}\right]_{Y = 0}} = 0,\quad 0 \le X \le {\rm \pi} \varphi, \end{gather}$$

(3.62)$$\begin{gather}{u}(X,0) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}. \end{gather}$$

These coupled patterned wall conditions for the oblique configuration highlight the strong effects of the nonlinear viscoplastic rheology, since new terms emerge in the shear stress components originated from the nonlinear viscosity.

Substituting the solutions in the form of (3.22) and (3.37) into the patterned wall conditions, i.e. (3.59), (3.60), (3.61) and (3.62), the following dual trigonometric series problems emerge:

(3.63)\begin{align} & {B_0}\left( {1 + \frac{{b( {{C_1} + B{{\sin }^2}\theta})}}{{{C_1}h}}} \right) - {A_0}\left( {\frac{{bB\sin \theta \cos \theta }}{{{C_1}h}}} \right) \nonumber\\ &\qquad + \sum_{n = 1}^\infty {{B_n}} \left[ {{{\hat w}_n}(0) - b\left( {1 + \frac{B}{{{C_1}}}{{\sin }^2}\theta } \right)\kappa \frac{{{\rm d}{{\hat w}_n}}}{{{\rm d} Y}}(0)} \right]\cos ({nX}) \nonumber\\ &\qquad + \sum_{n = 1}^\infty {{A_n}} b\frac{B}{{{C_1}}}\sin \theta \cos \theta \kappa \frac{{{{\rm d}^2}{{\hat \varPsi }_n}}}{{{\rm d}{Y^2}}}(0)\cos ({nX}) \nonumber\\ &\quad = b({B + {C_1}})\cos \theta,\quad 0 \leq X \leq {\rm \pi}\varphi, \end{align}

(3.64)\begin{equation} {B_0} + \sum_{n = 1}^\infty {{B_n}} {{\hat w}_n}(0)\cos ({nX}) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}, \end{equation}

(3.65)\begin{align} & {A_0}\left( {1 + \frac{{b( {{C_1} + B{{\cos }^2}\theta })}}{{{C_1}h}}} \right) - {B_0}\left( {\frac{{bB\sin \theta \cos \theta }}{{{C_1}h}}} \right) \nonumber\\ &\qquad + \sum_{n = 1}^\infty {{A_n}} \left[ {\frac{{{\rm d}{{\hat \varPsi }_n}}}{{{\rm d} Y}}(0) - b\left( {1 + \frac{B}{{{C_1}}}{{\cos }^2}\theta } \right)\kappa \frac{{{{\rm d}^2}{{\hat \varPsi }_n}}}{{{\rm d}{Y^2}}}(0)} \right]\cos ({nX}) \nonumber\\ &\qquad + \sum_{n = 1}^\infty {{B_n}} b\frac{B}{{{C_1}}}\sin \theta \cos \theta \kappa \frac{{{\rm d}{{\hat w}_n}}}{{{\rm d} Y}}(0)\cos \left( {nX} \right) \nonumber\\ &\quad = b({B + {C_1}})\sin \theta,\quad 0 \leq X \leq {\rm \pi}\varphi, \end{align}

(3.66)\begin{equation} {A_0} + \sum_{n = 1}^\infty {{A_n}} {\frac{{\rm d}\hat \varPsi_n}{{\rm d} Y}(0)}\cos ({nX}) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}. \end{equation}

Note that for a Newtonian fluid, as $B=0$, the dual series problems given above become decoupled and, since the solutions for $\hat w_n$ and $\hat \varPsi _n$ are identical to those of the longitudinal and transverse flow, respectively, we would simply have $B_n^\angle =B_n^\parallel \cos \theta$ and $A_n^\angle =A_n^\bot \sin \theta$ for $n=0,1,2,\ldots$. However, in a viscoplastic material, the dependence of $\hat w_n$ and $\hat \varPsi _n$ solutions on $\theta$ and the coupling between the ODEs and the patterned wall conditions make the problem different and more challenging.

Using the method developed for the longitudinal and transverse configurations, one can find the following coupled system of linear equations:

(3.67)$$\begin{gather} \sum_{n = 0}^N {{P_{mn}^{w\angle}}{B_n}} = {M_m^{w\angle}}- \sum_{n = 0}^N{{E_{mn}^{w\angle}}{A_n}}, \end{gather}$$

(3.68)$$\begin{gather}\sum_{n = 0}^N {{P_{mn}^{u\angle}}{A_n}} = {M_m^{u\angle}}- \sum_{n = 0}^N{{E_{mn}^{u\angle}}{B_n}}, \end{gather}$$

where

(3.69)$$\begin{gather} P_{m0}^{w\angle} = \left( {1 + \frac{b}{h}\left[ {1 + \frac{B}{{{C_1}}}{{\sin }^2}\theta } \right]} \right){I_3} +{I_4}, \end{gather}$$

(3.70)$$\begin{gather}E_{m0}^{w\angle} ={-} \frac{b}{h}\frac{B}{{{C_1}}}\sin \theta \cos \theta {I_3}, \end{gather}$$

(3.71)$$\begin{gather}P_{mn}^{w\angle} = \frac{1}{n}\left[ {{{\hat w}_n}(0) - b\left( {1 + \frac{B}{{{C_1}}}{{\sin }^2}\theta } \right)\kappa \frac{{{\rm d}{{\hat w}_n}}}{{{\rm d} Y}}(0)} \right]{I_1} + {\hat w_n}(0){I_2},\quad n > 0, \end{gather}$$

(3.72)$$\begin{gather}E_{mn}^{w\angle} = \frac{b}{n}\frac{B}{{{C_1}}}\sin \theta \cos \theta \kappa \frac{{{{\rm d}^2}{{\hat \varPsi }_n}}}{{{\rm d}{Y^2}}}(0){I_1},\quad n > 0, \end{gather}$$

(3.73)$$\begin{gather}{M_m^{w\angle}} = b({B + {C_1}}) \cos \theta I_3 \end{gather}$$

and

(3.74)$$\begin{gather} P_{m0}^{u\angle} = \left( {1 + \frac{b}{h}\left[ {1 + \frac{B}{{{C_1}}}{{\cos }^2}\theta } \right]}\right){I_3} + {I_4}, \end{gather}$$

(3.75)$$\begin{gather}{E_{m0}^{u\angle}} = {E_{m0}^{w\angle}}, \end{gather}$$

(3.76)$$\begin{gather}P_{mn}^{u\angle} = \frac{1}{n}\left[ {\frac{{{\rm d}{{\hat \varPsi }_n}}}{{{\rm d} Y}}(0) - b\left( {1 + \frac{B}{{{C_1}}}{{\cos }^2}\theta } \right)\kappa \frac{{{{\rm d}^2}{{\hat \varPsi }_n}}}{{{\rm d}{Y^2}}}(0)} \right]{I_1} + \frac{{{\rm d}{{\hat \varPsi }_n}}}{{{\rm d} Y}}(0){I_2},\quad n > 0, \end{gather}$$

(3.77)$$\begin{gather}E_{mn}^{u\angle} = \frac{b}{n}\frac{B}{{{C_1}}}\sin \theta \cos \theta \kappa \frac{{{\rm d}{{\hat w}_n}}}{{{\rm d} Y}}(0){I_1},\quad n > 0, \end{gather}$$

(3.78)$$\begin{gather}{M_m^{u\angle}} = b({B + {C_1}}) \sin \theta I_3. \end{gather}$$

Regarding the numerical procedure, in the first iteration, an initial guess for $\varDelta ^{(1)}_n$, $\varGamma ^{(1)}_n$ and $\varGamma ^{(2)}_n$ is considered, e.g. $\varDelta ^{(1,0)}_n=1$, $\varGamma ^{(1,0)}_n=0$ and $\varGamma ^{(2,0)}_n=0$. Then, one can find $\hat w_n$ and $\hat \varPsi _n$ based on (3.52) and (3.53). Therefore, in the first iteration, one can form the system of linear equations (3.67) and (3.68), while ignoring the rightmost term on the right-hand side of these linear systems (i.e. $E_{mn}^{w\angle } A_n$ and $E_{mn}^{u\angle } B_n$), and solve for $A_n$ and $B_n$. In the next iteration, the new values of $A_n$ and $B_n$ should be used to update the coefficients $\varDelta ^{(1)}_n$, $\varGamma ^{(1)}_n$ and $\varGamma ^{(2)}_n$ based on (3.56), (3.57) and (3.58), and, hence, to update $\hat w_n$ and $\hat \varPsi _n$ through (3.52) and (3.53). Then, the system of linear equations (3.67) and (3.68) can be solved, while using the new values of $A_n$ and $B_n$ in order to calculate the rightmost term on the right-hand side of these linear systems (i.e. $E_{mn}^{w\angle } A_n$ and $E_{mn}^{u\angle } B_n$), leading to finding updated values for $A_n$ and $B_n$. This iterative procedure continues until reaching converged values for $\varDelta ^{(1)}_n$, $\varGamma ^{(1)}_n$ and $\varGamma ^{(2)}_n$, which ensures the convergence of $A_n$ and $B_n$ as well (see Appendix C for an example of such a convergence).

3.3. Explicit-form solution

Let us attempt to derive an explicit-form solution for the perturbation velocity by means of analytical solutions to the dual trigonometric series problems, e.g. (3.25) and (3.26), emerging after applying the patterned wall conditions for the longitudinal, transverse and oblique flow configurations. To this end, we adopt a technique originally introduced by Sneddon (Reference Sneddon1966) and later extended to the study of Newtonian fluids in an SH channel by Belyaev & Vinogradova (Reference Belyaev and Vinogradova2010) to solve the resulting dual trigonometric series problem. As the first step, the dual series for the longitudinal and transverse configurations can be summarised in the following form:

(3.79)$$\begin{gather} {A^*_0}\left({1 + \frac{b}{h}}\right) + \sum_{n = 1}^\infty {{A_n^*}} [{1 + b\kappa {\lambda^*}n}]\cos (nX) = bS^*,\quad 0 \le X \le {\rm \pi}\varphi, \end{gather}$$

(3.80)$$\begin{gather}{A^*_0} + \sum_{n = 1}^\infty {{A_n^*}}\cos (nX) = 0,\quad {\rm \pi}\varphi \leq X \leq {\rm \pi}, \end{gather}$$

where

(3.81) $$\begin{gather} {\rm{for\ longitudinal}}:\left\{ \begin{array}{@{}l} A_0^* = {B_0},\\ A_n^* = {B_n},\\ {\lambda ^*} = {\lambda ^\parallel },\\ {S^*} = B + {C_1}, \end{array}\right. \end{gather}$$

(3.82) $$\begin{gather} {\rm{for\ transverse}} :\left\{ \begin{array}{@{}l} A_0^* = {A_0},\\ A_n^* = {A_n}({\lambda _2^ \bot - \lambda _1^ \bot})n,\\ {\lambda ^*} = \lambda _1^ \bot + \lambda _2^ \bot ,\\ {S^*} = B + {C_1}. \end{array}\right. \end{gather}$$

To obtain the dual series form (3.79) and (3.80) for the oblique configuration, further approximations are required. First, using the estimation (see Appendix B for more details)

(3.83)\begin{equation} \frac{B_n}{A_n} \approx{-}\alpha n \cot \theta,\quad \alpha > 0, \end{equation}

along with (3.56), (3.57) and (3.58) leads to the following approximations:

(3.84)$$\begin{gather} \varDelta _n^{(1)} \approx {\varDelta _1} = \alpha \frac{{({B + {C_1}}) ({{{({\lambda _1^\angle})}^2} - 1}) - B{{({\lambda _1^\angle})}^2} {{\cos }^2}\theta }}{{B({{{({\lambda _1^\angle})}^2} + 1})\lambda _1^\angle {{\sin }^2}\theta }}, \end{gather}$$

(3.85)$$\begin{gather}\varDelta _n^{(2)} \approx {\varDelta _2} = 1- {\varDelta _1}, \end{gather}$$

(3.86)$$\begin{gather}\varGamma _n^{(1)} \approx {\varGamma _1} ={-} \frac{1}{\alpha } \frac{{B({{{({\lambda _2^\angle})}^2} + 1})\lambda _2^\angle {{\sin }^2}\theta }}{{{(B + {C_1}})({{{({\lambda _2^\angle})}^2} - 1}) - B{{({\lambda _2^\angle})}^2}{{\cos }^2}\theta }}, \end{gather}$$

(3.87)$$\begin{gather}\varGamma _n^{(2)} \approx {\varGamma _2} = \frac{2}{\alpha}({{\varDelta _1}-1})\tan^2 \theta, \end{gather}$$

where, as should be expected, $\varGamma _{1}, \varGamma _{2} \to 0$ once $\theta \to 0$. Knowing the fact that $\varDelta _{1} \to 1$ once $\theta \to 90^\circ$, one can find

(3.88)\begin{equation} \alpha = \frac{{B({{{({\lambda _1^ \bot})}^2} + 1})\lambda _1^ \bot }}{{({B + {C_1}}) ({{{({\lambda _1^ \bot})}^2} - 1})}}. \end{equation}

These approximations enable us to obtain a reasonable estimate for the solution of $\hat w_n$ and $\hat \varPsi _n$ based on (3.52) and (3.53).

Using the estimation mentioned for ${B_n}/{A_n}$ and also estimating ${B_0}/{A_0} \approx \cot \theta$, which stems from a weak secondary flow stream in the $x'$ direction for the thick channel limit, we are able to make the patterned wall conditions (3.63), (3.64), (3.65) and (3.66) decoupled to find

(3.89) \begin{equation} {\rm{for\ oblique:}}\left\{\begin{array}{@{}l} {\rm{for}}\ w:\left\{\begin{array}{@{}l} A_0^* = {B_0},\quad A_n^* = {B_n} (1+\varGamma_1+\varGamma_2),\\ {\lambda ^*}={\lambda_w ^*} = \dfrac{{{\varUpsilon_1} \left({1 + \dfrac{B}{{{C_1}}}{{\sin }^2}\theta } \right) -{\varUpsilon_2} \dfrac{B}{{\alpha {C_1}}}{{\sin }^2}\theta }}{{1 + {\varGamma _1} + {\varGamma _2}}},\\ {S^*} = ({B + {C_1}})\cos \theta, \end{array}\right.\\ {\rm{for}}\ u:\left\{ \begin{array}{@{}l} A_0^* = {A_0},\quad A_n^* = {A_n}({ - {\varDelta _1} \lambda _1^\angle + \lambda _2^\angle - 1 + {\varDelta _1}})n,\\ {\lambda ^*} ={\lambda_u ^*}= \dfrac{{{\varUpsilon_2} \left({1 + \dfrac{B}{{{C_1}}}{{\cos }^2}\theta}\right) - {\varUpsilon_1} \dfrac{{\alpha B}}{{{C_1}}}{{\cos }^2}\theta}}{{{\varDelta _1} \lambda _1^\angle - \lambda _2^\angle + 1 - {\varDelta _1}}},\\ {S^*} = ({B + {C_1}})\sin \theta, \end{array}\right. \end{array}\right. \end{equation}

where

(3.90)\begin{equation} \left.\begin{gathered} {\varUpsilon _1} = ({\lambda _1^\angle + {\varGamma _1}{\mkern 1mu} \lambda _2^\angle + {\varGamma _2}}), \\ {\varUpsilon _2} = ({{\varDelta _1}{\mkern 1mu} {{({\lambda _1^\angle })}^2} - {{({{\mkern 1mu} \lambda _2^\angle })}^2} + 1 - {\varDelta _1}}). \end{gathered}\right\} \end{equation}

In (3.89), the terms corresponding to $w$-equations ($u$-equations) represent (3.63) and (3.64) ((3.65) and (3.66)) when written in the form of (3.79) and (3.80). One should note also that, $\theta \to 0$ ($\theta \to 90^\circ$), $A_n^*$, $\lambda ^*$ and $S^*$ for $w$-equations ($u$-equations) converge to those shown in (3.81) for the longitudinal ((3.82) for the transverse) flow.

Note that, in (3.79), the location of the lower yield surface in $X-Y$ coordinate, i.e. $Y=\kappa h$, is approximated with the corresponding value for the no-slip flow, i.e. $\kappa h \approx \kappa h_0$. As demonstrated by Rahmani & Taghavi (Reference Rahmani and Taghavi2022), with an increase in the slip number, before any plug formation at the liquid/air interface, $h$ only slightly deviates from $h_0$.

According to Sneddon (Reference Sneddon1966), to solve the dual series problem (3.79) and (3.80), one can write

(3.91)\begin{equation} {A^*_0} + \sum_{n = 1}^\infty {{A_n^*}} \cos (nX) = \cos \left(\frac{X}{2}\right)\int_X^{{\rm \pi} \varphi} {\frac{{\hbar(t)\,{\rm d} t}}{{\sqrt {\cos (X) - \cos (t)}}}},\quad 0 \le X \le {\rm \pi}\varphi \end{equation}

and, therefore, $A^*_0$ and $A^*_n$ are obtained as

(3.92)$$\begin{gather} {A^*_0} = \frac{1}{{\rm \pi} }\left( {\frac{{\rm \pi} }{{\sqrt 2 }} \int_0^{{\rm \pi} \varphi } {\hbar(t)\,{\rm d} t} }\right), \end{gather}$$

(3.93)$$\begin{gather}{A_n^*} = \frac{2}{{\rm \pi} }\left( {\frac{{\rm \pi} }{{2\sqrt 2 }}\int_0^{{\rm \pi} \varphi } {\hbar(t) [{{P_n}(\cos(t)) + {P_{n - 1}}(\cos (t))}]{\rm d} t} } \right), \end{gather}$$

where $P_n$ is the Legendre polynomial.

Integrating (3.79) in $[0\ X]$, while substituting (3.91) and (3.93), one can derive

(3.94)\begin{align} & b\kappa {\lambda ^*}\int_0^{{\rm \pi} \varphi } {\frac{{\hbar (t)}}{{\sqrt 2 }} \sum_{n = 1}^\infty {[{{P_n}({\cos (t)}) + {P_{n - 1}}({\cos (t)})}]\sin (nX)} }\,{\rm d} t \nonumber\\ &\quad = \left( {bS^* - \frac{{{A^*_0}b}}{h_0}} \right)X - \int_0^X {\cos \left(\frac{X}{2}\right) \int_X^{{\rm \pi} \varphi } {\frac{{\hbar (t)\,{\rm d} t}}{{\sqrt {\cos (X) - \cos (t)} }}}\,{\rm d} X}. \end{align}

The rightmost term on the right-hand side of (3.94) represents an integral of the slip velocity in $[0\ X]$, where $X \in [0\ {\rm \pi}\varphi ]$. The aforementioned term poses a challenge in obtaining an exact solution for the dual trigonometric series problem (Sneddon Reference Sneddon1966). In Belyaev & Vinogradova (Reference Belyaev and Vinogradova2010), this term was split into two integrals, and only one of them was considered in the solution, whereas the other was neglected (see Belyaev & Vinogradova (Reference Belyaev and Vinogradova2010) for details). Here, we instead make an approximation for this integral as follows:

(3.95)\begin{equation} \int_0^X {\cos \left(\frac{X}{2}\right)\int_X^{{\rm \pi} \varphi } {\frac{{\hbar (t)\,{\rm d} t}}{{\sqrt {\cos (X) - \cos (t)} }}} \,{\rm d} X} \approx \frac{{{X}}}{\varphi }A^*_0. \end{equation}

Such an approximation is based on the fact that the integral of the slip velocity in $[0\ X]$ may be linearly estimated based on the ${X}/{{\rm \pi} \varphi }$ fraction of the true integral, which is ${\rm \pi} A_0^*$, where $A_0^*$ represents the average slip velocity in $[0\ {\rm \pi}]$.

Substituting (3.95) in (3.94), according to Sneddon (Reference Sneddon1966), $\hbar (t)$ has the solution in the form of

(3.96)\begin{equation} \hbar (t) = \frac{2}{{\rm \pi} }\frac{{\rm{d}}}{{{\rm{d}}t}}\int_0^t {\frac{{\sin \left(\dfrac{\xi }{2}\right)}}{{\sqrt {\cos (\xi ) - \cos (t)} }} \left[ {\frac{{\left( {bS^* - {A^*_0}\left( {\dfrac{1}{\varphi } + \dfrac{b}{h_0}} \right)} \right)}}{{b\kappa {\lambda ^*}}}\xi } \right]{\rm d}\xi }. \end{equation}

Subsequently, using (3.92) and (3.96), $A^*_0$ can be obtained:

(3.97)\begin{equation} A_0^* = \frac{{2{S^*}\ln \left[{\sec \left({\dfrac{{{\rm \pi} \varphi}}{2}} \right)}\right]}}{{\kappa {\lambda ^*} + 2\left(\dfrac{1}{{b\varphi }} + \dfrac{1}{{{h_0}}}\right)\ln \left[{\sec \left({\dfrac{{{\rm \pi} \varphi}}{2}}\right)}\right]}}. \end{equation}

In addition, having $\hbar (t)$, $A_n^*$ (for $n\ne 0$) should be calculated using (3.93).

Considering the fact that $A_0^*$ represents the average slip velocity (in the $X$ or $z$ direction), (3.97) can be used to explicitly reveal the slip dynamics for any flow configuration type ($0 \le \theta \le 90^\circ$), thick channel geometry ($\ell \ll 1$, $0 < \varphi < 1$), liquid/air interface or slippery stripe dynamics (various $b$) and Bingham fluid rheology (various $B$).

4.1. Total velocity

The superimposition of the no-slip and perturbation velocities provides the total velocity profile. Since the perturbation field alters the position of the lower yield surface ($\kern0.7pt y=h$), impacting the velocity profiles, in general an iterative numerical scheme is needed to determine $h$, as discussed in § 4.1.1. Adhering to the no-slip boundary conditions at the upper wall and continuity of the velocity profile at $y=h$, we simply obtain the ensuing distribution for the total velocity:

(4.1)$$\begin{gather} U(x,y) = \left\{\begin{array}{@{}ll} {}[{C_1}y + {C_2}{y^2}] \sin \theta + u(x,y), & 0 \le y \le h,\\ {}[{C_1}h + {C_2}{h^2}] \sin \theta, & h \le y \le {h^u},\\ {}[{C_5}(2 - y) + {C_6}{(2 - y)^2}] \sin \theta, & {h^u} \le y \le 2, \end{array}\right. \end{gather}$$

(4.2)$$\begin{gather}V(x,y) = \left\{\begin{array}{@{}ll} v(x,y), & 0 \le y \le h,\\ 0, & h \le y \le {h^u},\\ 0, & {h^u} \le y \le 2, \end{array}\right. \end{gather}$$

(4.3)$$\begin{gather}W(x,y) = \left\{\begin{array}{@{}ll} {}[{C_1}y + {C_2}{y^2}] \cos \theta + w(x,y), & 0 \le y \le h,\\ {}[{C_1}h + {C_2}{h^2}] \cos \theta, & h \le y \le {h^u},\\ {}[{C_5}(2 - y) + {C_6}{(2 - y)^2}] \cos \theta, & {h^u} \le y \le 2, \end{array}\right. \end{gather}$$

where $U$, $V$ and $W$ are total velocity components in the $x$, $y$ and $z$ directions, respectively. We remind that $h^u$ is the location of the upper yield surface. Note that, after calculating $h$ (as explained in § 4.1.1), one can calculate $h^u$ and the unknown coefficients $C_5$ and $C_6$ (as discussed in § 4.1.2).

4.1.1. Finding lower yield surface location ($h$)

In the vicinity of the lower yield surface, i.e. far from the patterned wall, the total velocity vector is predominantly oriented in the $z'$ direction, with insignificant higher-order contributions in the $x'$ direction. Therefore, in each iteration of the numerical procedure explained in § 3.2.3, extra iterations are required to determine the point where the gradient of the total velocity with respect to the $y$ direction vanishes, leading to a vanishing magnitude of the total strain rate. Thus, one can write

(4.4)\begin{equation} \frac{{\rm d} \langle W' \rangle}{{\rm d}y} = \frac{{\rm d} U^b_0}{{\rm d}y} + \frac{{\rm d} \langle w^\angle \rangle}{{\rm d}y} \cos \theta + \frac{{\rm d} \langle u^\angle \rangle}{{\rm d}y} \sin \theta= 0,\quad \text{at }y=h, \end{equation}

leading to the following relation for $h$:

(4.5)\begin{equation} h = \frac{{ - {C_1} - \sqrt {C_1^2 + 8{C_2}({{B_0}\cos \theta + {A_0}\sin \theta})} }}{{4{C_2}}}. \end{equation}

The solution can be achieved via a nested iteration, explained as follows. At each iteration, once coefficients $A_0$ and $B_0$ are updated (in an inner loop), extra iterations are required to find a converged value of $h$ using (4.5) and the patterned wall condition, while keeping the former values of $\hat w_n$ and $\hat \varPsi _n$. After finding a converged $h$, $\hat w_n$ and $\hat \varPsi _n$, and, hence, $A_0$ and $B_0$ are updated through the next iteration. Our results (omitted for brevity) reveal that a converged value of $h$ is typically achieved after few iterations, with a relative error of less than $10^{-6}$.

4.1.2. Finding $C_5$, $C_6$ and $h^u$

The remaining unknowns $C_5$, $C_6$ and $h^u$ can be obtained using the following three conditions, ensuring the velocity profile continuity (4.6) and zero strain-rate magnitude (4.7) at the upper yield surface, while maintaining the fixed flow rate (4.8)

(4.6)$$\begin{gather} {C_5}(2 - {h^u}) + {C_6}{(2 - {h^u})^2} = {C_1}h + {C_2}{h^2}, \end{gather}$$

(4.7)$$\begin{gather}{C_5} + 2{C_6}(2 - {h^u}) = 0, \end{gather}$$

(4.8)$$\begin{gather}U_{ave} \sin \theta + W_{ave} \cos \theta =1, \end{gather}$$

where $U_{ave}$ and $W_{ave}$ are the cross-sectional averaged total velocities in the $x$ and $z$ directions, respectively, and they are calculated through

(4.9) \begin{equation} \left.\begin{gathered} \left[\begin{array}{@{}l@{}} \displaystyle\int_0^h {({{C_1}y + {C_2}{y^2}})\,{{\rm d}y}} + \int_h^{{h^u}} {({{C_1}h + {C_2}{h^2}})\,{{\rm d}y}} \\ + \displaystyle\int_{{h^u}}^2 {({{C_5}(2 - y) + {C_6}{{(2 - y)}^2}})\,{{\rm d}y}} \end{array}\right]\ell \cos \theta \\ \quad + \int_{ - \ell /2}^{\ell /2}{\int_0^h {w(x,y)\,{\rm d}\kern0.7pt x\,{\rm d} y}} = 2\ell {W_{ave}}, \\ \left[\begin{array}{@{}l@{}} \displaystyle\int_0^h {({{C_1}y + {C_2}{y^2}})\,{{\rm d}y}} + \int_h^{{h^u}} {({{C_1}h + {C_2}{h^2}})\,{{\rm d}y}} \\ + \displaystyle\int_{{h^u}}^2 {({{C_5}(2 - y) + {C_6}{{(2 - y)}^2}})\,{{\rm d}y}} \end{array}\right]\sin \theta \\ \quad + \int_0^h {u(x,y)\,{{\rm d}y}} = 2{U_{ave}}. \end{gathered}\right\} \end{equation}

One can simply show that

(4.10)\begin{equation} \int_{ - \ell /2}^{\ell /2} {\int_0^h {w(x,y)\,{\rm d}\kern0.7pt x\,{\rm d} y}} = B_0 \ell h/2, \end{equation}

and

(4.11)\begin{equation} \int_0^h {u(x,y)\,{{\rm d}y}} = A_0 h/2. \end{equation}

Thus, the solution of (4.6), (4.7) and (4.8) provides us with the values of $C_5$, $C_6$ and $h^u$.

4.2. Onset of SH wall plug formation

For viscoplastic flows over SH surfaces, it is expected that, at some critical value of the slip number (i.e. $b_{cr}$, which is relatively large), the no-shear condition becomes locally met on the liquid/air interface. On the other hand, assuming no formation of nano-bubbles on the slippery stripes of CP surfaces, herein, we may consider partial-shear condition on the slippery stripes; thus, we can assume $b < b_{cr}$ for CP surfaces (Lee et al. Reference Lee, Charrault and Neto2014). Therefore, on an SH wall, at $y=0$ and $x=0$ (i.e. the middle of the slip zone), we might have $\partial \{ {W,U} \}/\partial y = 0$ (due to the no-shear condition) and $\partial \{ {W,U} \}/\partial x = 0$ (due to the symmetry with respect to $x=0$ for a creeping flow). These lead to the zero strain rate magnitude and, hence, formation of an unyielded plug zone on the liquid/air interface, called the SH wall plug.

Before proceeding, to elaborate more on the fact that our model is capable of capturing the onset of no-shear condition, let us consider the transverse flow (i.e. $\theta =90^\circ$), for which the Navier slip law (2.9) reduces to

(4.12)\begin{equation} u_s=b \tau_{xy}|_{y=0}=b\left[\left(1+\frac{B}{\dot \gamma} \right) {\dot \gamma_{xy}}\right]_{y=0}. \end{equation}

Based on (4.12), the local slip length on the liquid/air interface is $\chi _l=b( 1+{B}/{\dot \gamma })_{y=0}$. An increase in the slip number ($b$) leads to a more slippery interface, thus, reducing the strain-rate magnitude ($\dot \gamma$) on the interface. Therefore, based on the above Navier slip law, an increase in $b$ and, simultaneously, a decrease in $\dot \gamma$ lead to growth in the local slip length ($\chi _l$). Theoretically, a further increase in $b$ leads to a zero strain-rate magnitude at $x=0$ and $y=0$ (i.e. the middle of the liquid/air interface); thus, an infinite local slip length ($\chi _l \to \infty$) is achieved representing the onset of the no-shear condition at the interface. Therefore, at the onset of the no-shear condition, the local slip length ($\chi _l$) becomes infinite and the slip number ($b$) has a finite value. In the other words, based on the perturbed form of the Navier slip law for the transverse flow (see (3.39)), an increase in $b$ eventually leads to the condition $\partial U/\partial y = \kappa \partial u/\partial Y+C_1=0$ at $x=0$ and $y=0$ (i.e. the no-shear condition on the liquid/air interface), causing the local shear stress to become $\tau _{xy}=B$, which describes the onset of SH wall plug formation. It is important to emphasise that our proposed solutions are strictly valid only up to the establishment of the SH wall plug. In fact, through the Navier slip law used in our work, we are able to apply the partial-shear ($b< b_{cr}$) and no-shear ($b \approx b_{cr}$) conditions at the liquid/air interface based on the developed mathematical solutions. To accurately determine the critical slip number ($b_{cr}$), which indicates the onset of SH wall plug formation, it is necessary to increase the slip number until the no-shear condition is attained at $x=0$ and $y=0$, using the semi-analytical solution approach. We elaborate more on the flow dynamics close to and at the no-shear condition in Appendix E.

To obtain an explicit relation for the onset of the SH wall plug formation, we exploit the nearly uniform distribution of the shear stress/strain at the liquid/air interface of the SH wall, with the exception of two peaks near the groove edges (Teo & Khoo Reference Teo and Khoo2009; Rahmani & Taghavi Reference Rahmani and Taghavi2022). This fact has been employed in Schönecker et al. (Reference Schönecker, Baier and Hardt2014) to derive analytical solutions for the flow of Newtonian fluids over micro-structured surfaces in the Cassie state, with an enclosed second fluid inside the cavities. By assuming a no-shear condition around the liquid/air interface at the critical slip condition, one can estimate $C_1 \cos \theta \approx \kappa {\partial w}/{\partial Y}$ and $C_1 \sin \theta \approx \kappa {\partial u}/{\partial Y}$ at $Y=0$ and $0\le X \le {\rm \pi}\varphi$. Using these estimations, integrating (3.59) and (3.61) over the interval $[0;{\rm \pi} \varphi ]$, summing the resulting terms, and performing some algebra, one can find the following relation for $b_{cr}$:

(4.13)\begin{equation} {b_{cr}} \approx \frac{{{{\langle w_s^\angle \rangle }_{cr}} + {{\langle u_s^\angle \rangle }_{cr}}}}{{B\varphi ({\sin \theta + \cos \theta})}}. \end{equation}

Since in (4.13), ${{\langle w_s^\angle \rangle }_{cr}}$ and ${{\langle u_s^\angle \rangle }_{cr}}$ are themselves functions of $b_{cr}$, the following equation could be solved to calculate $b_{cr}$:

(4.14) \begin{equation} \frac{{2{b_{cr}}{C_1}{\zeta _2}({{b_{cr}}\varphi \kappa {C_1} ({\lambda _u^*\cos \theta \,{+}\, \lambda _w^*\sin \theta }) \,{+}\, 2({\sin \theta + \cos \theta})({{b_{cr}}\varphi {\zeta _2} \,{+}\, {\zeta _1}})})}}{{B({\sin \theta \,{+}\, \cos \theta}) ({{b_{cr}}\varphi \kappa {C_1}\lambda _w^* \,{+}\, 2({{b_{cr}}\varphi {\zeta _2} \,{+}\, {\zeta _1}})})({{b_{cr}}\varphi \kappa {C_1}\lambda _u^* + 2({{b_{cr}}\varphi {\zeta _2} + {\zeta _1}})})}} \,{-}\, {b_{cr}} \,{=}\, 0, \end{equation}

where

(4.15)\begin{equation} \left.\begin{gathered} {\zeta _1} = {C_1}\ln \left[ {\sec \left( {\frac{{{\rm \pi} \varphi }}{2}} \right)} \right], \\ {\zeta _2} = ({B + {C_1}})\ln \left[ {\sec \left( {\frac{{{\rm \pi} \varphi }}{2}} \right)} \right]. \end{gathered}\right\} \end{equation}

For the longitudinal and transverse flow limits, the solution can be simplified to

(4.16)$$\begin{gather} b_{cr}^\parallel{=} \frac{{2{C_1}\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}{{\left( {\kappa {\lambda ^\parallel } + 2\left( {1 + \dfrac{B}{{{C_1}}}} \right)\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]} \right)B\varphi }}, \end{gather}$$

(4.17)$$\begin{gather}b_{cr}^ \bot = \frac{{2{C_1}\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}{{\left( {\kappa \left( {\lambda _1^ \bot + \lambda _2^ \bot } \right) + 2\left( {1 + \dfrac{B}{{{C_1}}}} \right)\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]} \right)B\varphi }}. \end{gather}$$

These relations can be also rewritten in the form

(4.18)\begin{equation} \frac{1}{{\{{b_{cr}^\parallel ,b_{cr}^ \bot}\}\varphi }} = \frac{{B\{{{\lambda ^\parallel},({\lambda _1^ \bot + \lambda _2^ \bot})}\}}}{{2{C_1}}} \left( {\frac{\kappa }{{\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}}\right) + \frac{B}{{{C_1}}}\left( {1 + \frac{B}{{{C_1}}}} \right). \end{equation}

Before proceeding, let us mention that the above-calculated $b_{cr}$ for the longitudinal and transverse flows is not affected by the assumption of the scalar slip number, since for each flow configuration we independently calculate the slip number that causes the no-shear condition. However, the scalar slip number assumption might affect the distribution of $b_{cr}$ vs $\theta$ for the oblique flows, i.e. the transition of $b_{cr}$ from $\theta =0$ to $\theta =90^\circ$; this feature is less important compared with the fact that the values of $b_{cr}$ for the longitudinal and transverse flows can provide us with accurate upper and lower bounds of the critical slip number.

4.3. Effective slip length

The effective slip length, which characterises the slip dynamics of a flow system, is defined as the average slip velocity over the average of total velocity gradient normal to the wall at $y=0$. In fact, it can represent an imaginary distance within the slippery wall where the average total velocity becomes zero, thereby, providing a physical insight into the average slip behaviour. The tensorial form of the effective slip length can be used to compute the slip velocities in different directions. Considering the longitudinal and transverse configurations, the tensor is expressed as

(4.19) \begin{equation} \left({\begin{array}{@{}c@{}} {\langle {w_s^\parallel }\rangle }\\ {\langle {u_s^ \bot }\rangle } \end{array}}\right) = \left( {\begin{array}{@{}cc@{}} {{\chi ^\parallel }} & 0\\ 0 & {{\chi ^ \bot }} \end{array}} \right)\left( {\begin{array}{@{}c@{}} {\langle {{\partial _{y_0}}{W^\parallel }}\rangle }\\ {\langle {{\partial _{y_0}}{U^ \bot }}\rangle} \end{array}}\right), \end{equation}

where $\langle \cdot \rangle$ is the averaging operator in the $x$ direction, from $-\ell /2$ to $\ell /2$ (i.e. one period of the patterned wall). In addition, $w_s=w(x,0)$ and $u_s=u(x,0)$ are the slip velocities and, as explained earlier, the symbols $\parallel$ and $\bot$ represent the longitudinal and transverse configurations, respectively. In addition, the symbol ${{\partial _{y_0}}}$ is the gradient operator in the $y$ direction and at $y=0$, which is used to mark the total velocity shear gradient at the SH wall. Accordingly, ${\chi ^\parallel }$ and ${\chi ^\bot }$ represent the effective slip length for the longitudinal and transverse configurations, respectively. From the definition and considering (3.22) and (3.37), one can write

(4.20)$$\begin{gather} \chi^\parallel{=} \frac{{{B_0}}}{{{C_1} - \dfrac{{{B_0}}}{h}}}, \end{gather}$$

(4.21)$$\begin{gather}\chi^\bot = \frac{{{A_0}}}{{{C_1} - \dfrac{{{A_0}}}{h}}}. \end{gather}$$

One can correlate the average values of the slip velocities and the total velocity gradients in the oblique flow configuration (presented by the symbol $\angle$ throughout the manuscript) with those of the longitudinal and transverse flow configurations as

(4.22)$$\begin{gather} \langle {w_s^\angle } \rangle = \langle {w_s^\parallel } \rangle {g_1}, \end{gather}$$

(4.23)$$\begin{gather}\langle {u_s^\angle }\rangle = \langle {u_s^ \bot } \rangle {f_1}, \end{gather}$$

and

(4.24)$$\begin{gather} {\langle {{\partial _{y_0}}{W^\angle }}\rangle } = {\langle {{\partial _{y_0}}{W^\parallel }}\rangle } {g_2}, \end{gather}$$

(4.25)$$\begin{gather}{\langle {{\partial _{y_0}}{U^\angle }}\rangle } = {\langle {{\partial _{y_0}}{U^\bot }}\rangle } {f_2}. \end{gather}$$

As discussed in § 3.2.3, one can simply show that, for the Newtonian fluids, $g_1=g_2=\cos \theta$ and $f_1=f_2=\sin \theta$.

One can then calculate the average of the slip velocities in the $z'$ and $x'$ directions as

(4.26)$$\begin{gather} \langle {w_s^{'\angle }}\rangle = \langle {u_s^\angle } \rangle \sin \theta + \langle {w_s^\angle } \rangle \cos \theta, \end{gather}$$

(4.27)$$\begin{gather}\langle {u_s^{'\angle }} \rangle = \langle {u_s^\angle } \rangle \cos \theta - \langle {w_s^\angle}\rangle \sin \theta. \end{gather}$$

Therefore, using (4.22), (4.23), (4.26) and (4.27), one can write

(4.28) \begin{equation} \left( {\begin{array}{@{}c@{}} {\langle {w_s^{'\angle }} \rangle }\\ {\langle {u_s^{' \angle }} \rangle } \end{array}} \right) = \left( {\begin{array}{@{}cc@{}} {{g_1}\cos \theta } & {{f_1}\sin \theta }\\ {- {g_1}\sin \theta } & {{f_1}\cos \theta } \end{array}} \right)\left( {\begin{array}{@{}c@{}} {\langle {w_s^\parallel } \rangle }\\ {\langle {u_s^ \bot } \rangle } \end{array}}\right). \end{equation}

Following the same procedure for the total velocity gradients leads to

(4.29) \begin{equation} \left({\begin{array}{@{}c@{}} {\langle {{\partial _{{y_0}}}{W^{'\angle }}}\rangle }\\ {\langle {{\partial _{{y_0}}}{U^{'\angle }}} \rangle } \end{array}}\right) = \left( {\begin{array}{@{}cc@{}} {{g_2}\cos \theta } & {{f_2}\sin \theta }\\ {- {g_2}\sin \theta } & {{f_2}\cos \theta } \end{array}}\right)\left( {\begin{array}{@{}c@{}} {\langle {{\partial _{{y_0}}}{W^\parallel }} \rangle }\\ {\langle {{\partial _{{y_0}}}{U^ \bot }} \rangle } \end{array}}\right). \end{equation}

Finally, considering (4.19), (4.28) and (4.29), the tensorial form of the effective slip length is obtained as

(4.30) \begin{equation} \left({\begin{array}{@{}c@{}} {\langle {w_s^{'\angle }}\rangle }\\ {\langle {u_s^{'\angle }}\rangle } \end{array}} \right) = \left( {\begin{array}{@{}cc@{}} {\chi _{z'z'}^\angle } & {\chi _{z'x'}^\angle }\\ {\chi _{x'z'}^\angle } & {\chi _{x'x'}^\angle } \end{array}} \right)\left( {\begin{array}{@{}c@{}} {\langle {{\partial _{{y_0}}}{W^{'\angle }}}\rangle }\\ {\langle {{\partial _{{y_0}}}{U^{'\angle }}}\rangle } \end{array}}\right), \end{equation}

where

(4.31)$$\begin{gather} \chi _{z'z'}^\angle = \frac{{{g_1}}}{{{g_2}}}{\chi ^\parallel }{\cos ^2}\theta + \frac{{{f_1}}}{{{f_2}}}{\chi ^ \bot }{\sin ^2}\theta, \end{gather}$$

(4.32)$$\begin{gather}\chi _{x'x'}^\angle = \frac{{{g_1}}}{{{g_2}}}{\chi ^\parallel }{\sin ^2}\theta + \frac{{{f_1}}}{{{f_2}}}{\chi ^ \bot }{\cos ^2}\theta, \end{gather}$$

(4.33)$$\begin{gather}\chi _{x'z'}^\angle ={-} \left( {\frac{{{g_1}}}{{{g_2}}}{\chi ^\parallel } - \frac{{{f_1}}}{{{f_2}}}{\chi ^ \bot }} \right)\sin \theta \cos \theta, \end{gather}$$

(4.34)$$\begin{gather}\chi _{z'x'}^\angle = \chi _{x'z'}^\angle. \end{gather}$$

Equations (4.20), (4.21), (4.30), (4.31), (4.32), (4.33) and (4.34) represent the exact form of the effective slip length tensor. However, using the explicit-form solution, one can also find

(4.35)$$\begin{gather} {\chi ^\parallel } = \frac{{2\left( {1 + \dfrac{B}{{{C_1}}}} \right)\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}{{\kappa {\lambda ^\parallel } + 2\left( {\dfrac{1}{{b\varphi }} - \dfrac{B}{{{C_1}}}\left( {1 + \dfrac{B}{{{C_1}}}} \right)} \right)\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}, \end{gather}$$

(4.36)$$\begin{gather}{\chi ^ \bot } = \frac{{2\left( {1 + \dfrac{B}{{{C_1}}}} \right)\ln \left[ {\sec \left({\dfrac{{{\rm \pi} \varphi }}{2}} \right)}\right]}}{{\kappa ({\lambda _1^ \bot + \lambda _2^ \bot }) + 2\left({\dfrac{1}{{b\varphi }} - \dfrac{B}{{{C_1}}}\left( {1 + \dfrac{B}{{{C_1}}}} \right)} \right)\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}, \end{gather}$$

and

(4.37)$$\begin{gather} {g_1} = \frac{{b\varphi {\kern 1pt} \kappa {C_1}{\kern 1pt} \lambda^\parallel{+} 2( {b\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) + {C_1}})\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}{{b\varphi {\kern 1pt} \kappa {\kern 1pt} {C_1}\lambda _w^* + 2({b\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) + {C_1}})\ln \left[ {\sec \left({\dfrac{{{\rm \pi} \varphi }}{2}} \right)}\right]{\kern 1pt} }}\,\cos \theta, \end{gather}$$

(4.38)$$\begin{gather}{g_2} = \frac{{b\varphi {\kern 1pt} \kappa {\kern 1pt} C_1^2\lambda _w^* - 2({bB\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) - C_1^2})\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)}\right]{\kern 1pt} }}{{b\varphi {\kern 1pt} \kappa {\kern 1pt} C_1^2{\lambda ^\parallel } - 2({bB\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) - C_1^2})\ln \left[ {\sec \left({\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}{g_1}, \end{gather}$$

(4.39)$$\begin{gather}{f_1} = \frac{{b\varphi {\kern 1pt} \kappa {C_1}{\kern 1pt} ({\lambda _1^ \bot + \lambda _2^ \bot }) + 2({b\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) + {C_1}})\ln \left[{\sec \left({\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}{{b\varphi {\kern 1pt} \kappa {\kern 1pt} {C_1}\lambda _u^* + 2({b\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) + {C_1}})\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]{\kern 1pt} }}\sin \theta, \end{gather}$$

(4.40)$$\begin{gather}{f_2} = \frac{{b\varphi {\kern 1pt} \kappa {\kern 1pt} C_1^2\lambda _u^* - 2( {bB\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) - C_1^2})\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]{\kern 1pt} }}{{b\varphi {\kern 1pt} \kappa {\kern 1pt} C_1^2({\lambda _1^ \bot + \lambda _2^ \bot}) - 2({bB\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) - C_1^2})\ln \left[{\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}\,{f_1}. \end{gather}$$

It is worth mentioning that, when $B \to 0$, $g_1,g_2 \to \cos \theta$ and $f_1,f_2 \to \sin \theta$, recovering the effective slip length tensor for the Newtonian fluids, as presented in (Bazant & Vinogradova Reference Bazant and Vinogradova2008; Vinogradova & Belyaev Reference Vinogradova and Belyaev2011) (with a negative sign difference in $\chi _{x'z'}$ due to the choice of the coordinate system).

At the critical slip number, i.e. the no-shear condition, one can find

(4.41)\begin{equation} {\left( {\frac{{{\chi ^\parallel }}}{{{\chi ^ \bot }}}} \right)_{cr}} = \frac{{\lambda _1^ \bot + \lambda _2^ \bot }}{{{\lambda ^\parallel }}}, \end{equation}

which reveals that, at the no-shear condition, the ratio of the largest effective slip length (i.e. $\chi ^\parallel$ for the longitudinal flow) to the smallest one (i.e. $\chi ^\bot$ for the transverse flow) is only a function of the Bingham number. Interestingly, this ratio remains almost equal to $2$ for a wide range of Bingham numbers, identical to the corresponding value for Newtonian fluids, as discussed in more detail in § 6.

4.4. Slip angle

Let us define the slip angle ($s$) as the angle between the average slip velocity vector and the $z$ axis. Thus, one can obtain

(4.42)\begin{equation} s = {\tan ^{ - 1}}\left( {\frac{{\langle {u_s^ \angle}\rangle }}{{\langle {w_s^\angle}\rangle }}}\right). \end{equation}

Using the developed explicit-form solution, we find

(4.43)\begin{equation} s = {\tan ^{ - 1}}\left({\frac{{b\varphi {\kern 1pt} \kappa {\kern 1pt} {C_1}\lambda _w^* + 2({b\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) + {C_1}})\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)}\right] {\kern 1pt} }}{{b\varphi {\kern 1pt} \kappa {\kern 1pt} {C_1}\lambda _u^* + 2({b\varphi {\kern 1pt} ({B + {\kern 1pt} {C_1}}) + {C_1}})\ln \left[ {\sec \left( {\dfrac{{{\rm \pi} \varphi }}{2}} \right)} \right]}}\tan \theta}\right). \end{equation}

Note that the slip angle must be always smaller than $\theta$, since the flow streamlines near the patterned wall are strongly deviated towards the $z$ axis, i.e. the groove (stripe) direction. Therefore, we are interested in quantifying such a difference, i.e. $\theta -s$, as a variable of this study vs our flow parameters. Such an analysis explores the direction where the flow slippage occurs on the patterned wall, thus, revealing the direction towards which the average shear force is applied to the patterned wall by the Bingham fluid. Such an analysis may be of interest for practical applications, such as passive mixing (Stroock et al. Reference Stroock, Dertinger, Ajdari, Mezic, Stone and Whitesides2002b; Vagner & Patlazhan Reference Vagner and Patlazhan2019) and particle fractionation and manipulation (Asmolov et al. Reference Asmolov, Dubov, Nizkaya, Kuehne and Vinogradova2015, Reference Asmolov, Dubov, Nizkaya, Harting and Vinogradova2018; Nizkaya et al. Reference Nizkaya, Asmolov, Harting and Vinogradova2020) in channels with patterned walls, where the flow streamline structure near the patterned wall is important.

4.5. Flow mixing

Following the literature (Vinogradova & Belyaev Reference Vinogradova and Belyaev2011), in this section, the ratio between the flow rates generated in the $x'$ and $z'$ directions is used to quantify the mixing capability for the oblique flow configuration. The study of this phenomenon is important as the effects of various flow parameters can be analysed on mixing efficiency, leading to finding the optimum parameter combination corresponding to the maximum mixing. In our Bingham flow, since mixing occurs in the lower yielded zone, we define a mixing index ($I_M$) as the ratio between the flow rates in the $x'$ and $z'$ directions in the lower yielded zone as

(4.44)\begin{equation} I_M = \left|{\frac{{{Q_{x'}}}}{{{Q_{z'}}}}}\right| = \frac{{|{{A_0}\cos \theta - {B_0}\sin \theta}|}}{{{C_1}h + \dfrac{2}{3}{C_2}{h^2} + {A_0}\sin \theta + {B_0}\cos \theta }}. \end{equation}

We can simply evaluate the maximum value of $I_M$, representing the optimum mixing condition, through finding the value of $\theta$ where ${\partial I_M}/{\partial \theta }=0$.