Multiphase unsteady Stokes flow over grooved surfaces: an analytical study

TL;DR

This paper presents an analytical study of multiphase unsteady Stokes flow over grooved surfaces, deriving solutions and validating them with numerical simulations, focusing on flow response and shear stress under oscillatory conditions.

Contribution

It introduces analytical expressions for multiphase unsteady Stokes flow over grooved textures, including effects of oscillation frequency and interface position, validated by numerical simulations.

Findings

Analytical solutions match numerical simulations well.

Flow transfer function depends on oscillation frequency and interface location.

Effective slip length and shear stress are characterized for grooved surfaces.

Abstract

This paper studies multiphase flow within grooved textures exposed to external unsteadiness. We derive analytical expressions for multiphase unsteady Stokes flow driven by oscillating streamwise/spanwise velocity in the presence of periodic grooves. Good agreement is obtained between the analytical solution and direct numerical simulations performed using the volume of fluid method. The effects of oscillation frequency, Reynolds number, and the multiphase interface location on the transfer function between the input signal external to the groove and output near the interface, are examined. The effective slip length and the shear stress over the grooved plane are described.

Tables

Table 1. Table 1 : Convergence of ℜ ⁡ ( B 0 + C 0 ) subscript 𝐵 0 subscript 𝐶 0 \Re\left(B_{0}+C_{0}\right) , a = 0.6 𝑎 0.6 a=0.6 , b = 1.25 𝑏 1.25 b=1.25 , h = 0.9 ​ b ℎ 0.9 𝑏 h=0.9b , H = 5 𝐻 5 H=5 , μ r = 0.02 subscript 𝜇 𝑟 0.02 \mu_{r}=0.02 .

$L^2\omega /\nu$\M	5	10	15	20	25	30	35	40	45
0.1875	0.0232	0.0234	0.0234	0.0234	0.0235	0.0235	0.0235	0.0235	0.0235
0.0625	0.0366	0.0369	0.0369	0.0369	0.0370	0.0370	0.0370	0.0370	0.0370

Table 2. Table 2 : Convergence of ℜ ⁡ ( i ​ L 2 ​ ω ν ​ ( E 0 − F 0 ) ) 𝑖 superscript 𝐿 2 𝜔 𝜈 subscript 𝐸 0 subscript 𝐹 0 \Re\left(\sqrt{\frac{iL^{2}\omega}{\nu}}(E_{0}-F_{0})\right) , a = 0.6 𝑎 0.6 a=0.6 , b = 1.25 𝑏 1.25 b=1.25 , h = 0.9 ​ b ℎ 0.9 𝑏 h=0.9b , H = 5 𝐻 5 H=5 , μ r = 0.02 subscript 𝜇 𝑟 0.02 \mu_{r}=0.02 .

$L^2\omega /\nu$\M	10	20	30	40	50	60	70
0.1875	0.0069	0.0068	0.0068	0.0068	0.0067	0.0067	0.0067
0.0625	0.0102	0.0101	0.0101	0.0100	0.0100	0.0100	0.0100

Table 3. Table 3 : Location of the forcing H 𝐻 H , the penetration rate 1 − h / b 1 ℎ 𝑏 1-h/b , the relevant non-dimensional number L 2 ​ ω ν superscript 𝐿 2 𝜔 𝜈 \frac{L^{2}\omega}{\nu} , the viscosity ratio μ r subscript 𝜇 𝑟 \mu_{r} for all cases solved from analytical solution, and the applicable R ​ e τ 𝑅 subscript 𝑒 𝜏 Re_{\tau} . For all cases, L = 18 / 3500 𝐿 18 3500 L=18/3500 is considered.

$H$	Penetration rate	$\frac{L^2\omega }{\nu }$	${\mu }_r$	$R{e}_{\tau }$
$5.40$	$0,0.15,0.5,0.9$	$0.04-4.23$	$0.02,0.37,30.00$	$180$
$2.43$	$0,0.15,0.5,0.9$	$0.04-4.23$	$0.02,0.37,30.00$	$400$
$1.65$	$0,0.15,0.5,0.9$	$0.04-4.23$	$0.02,0.37,30.00$	$590$

Equations

$$u$$

$$y$$

$$t$$

$$\frac{\partial u}{\partial t}=\frac{\nu }{L^2\omega }\left(\frac{{\partial }^{2}u}{\partial z^2}+\frac{{\partial }^{2}u}{\partial y^2}\right),$$

$$u=\Re (\hat{u}\exp (it))=\Re (Y(y)Z(z)\exp (it)).$$

$$i\hat{u}=\frac{\nu }{L^2\omega }\left(\frac{{\partial }^2\hat{u}}{\partial y^2}+\frac{{\partial }^2\hat{u}}{\partial z^2}\right).$$

$$\begin{array}{rl}{\hat{u}}_{\mathrm{I}},& 0⩽z<a,\\ 0,& a<z⩽1;\end{array}$$

$$\frac{\partial {\hat{u}}_{\mathrm{I}}(0,z)}{\partial y}=\frac{\partial {\hat{u}}_{\mathrm{II}}(0,z)}{\partial y};0⩽z<a$$

$${\hat{u}}_{\mathrm{I}}(-b+h,z)={\hat{u}}_{\mathrm{III}}(-b+h,z);0⩽z<a,$$

$$\frac{\partial {\hat{u}}_{\mathrm{I}}(-b+h,z)}{\partial y}={\mu }_r\frac{\partial {\hat{u}}_{\mathrm{III}}(-b+h,z)}{\partial y};0⩽z<a$$

$$\frac{\partial {\hat{u}}_{\mathrm{I}}(y,0)}{\partial z}$$

$${\hat{u}}_{\mathrm{I}}(y,a)$$

$$Y_n(y)=A_{1n}\exp \left(\sqrt{{\alpha }_n^2+\frac{i{L}^2\omega }{\nu }}y\right)+A_{2n}\exp \left(-\sqrt{{\alpha }_n^2+\frac{i{L}^2\omega }{\nu }}y\right).$$

$${\hat{u}}_{\mathrm{I}}(y,z)$$

$$+A_{2n}\exp \left[-\sqrt{{\alpha }_n^2+\frac{i{L}^2\omega }{\nu }}(y+2b-2h)\right]\}.$$

$$\frac{\partial {\hat{u}}_{\mathrm{II}}(y,0)}{\partial z}$$

$$\frac{\partial {\hat{u}}_{\mathrm{II}}(y,1)}{\partial z}$$

$$\hat{u}(H,z)$$

$$Y_n(y)=B_n\exp \left(\sqrt{{\gamma }_n^2+\frac{i{L}^2\omega }{\nu }}y\right)+C_n\exp \left(-\sqrt{{\gamma }_n^2+\frac{i{L}^2\omega }{\nu }}y\right).$$

$${\hat{u}}_{\mathrm{II}}(y,z)=B_0\exp \left[\sqrt{\frac{i{L}^2\omega }{\nu }}y\right]+C_0\exp \left[-\sqrt{\frac{i{L}^2\omega }{\nu }}y\right]$$

$$+\sum _{n=1}^{\infty }\cos ({\gamma }_{n}z)\left[B_n\exp \left(\sqrt{{\gamma }_n^2+\frac{i{L}^2\omega }{\nu }}y\right)+C_n\exp \left(-\sqrt{{\gamma }_n^2+\frac{i{L}^2\omega }{\nu }}y\right)\right].$$

$$B_0\exp \left[\sqrt{\frac{L^2\omega }{2\nu }}\left(1+i\right)H\right]+C_0\exp \left[-\sqrt{\frac{L^2\omega }{2\nu }}\left(1+i\right)H\right]=1.$$

$${\hat{u}}_{\mathrm{III}}(y,z)$$

$$-\exp \left(-\sqrt{{\alpha }_n^2+\frac{i{L}^2\omega }{\nu }}(y+b+h)\right)].$$

$$(A_{1m}+A_{2m})\exp \left[-\sqrt{{\alpha }_m^2+\frac{i{L}^2\omega }{\nu }}(b-h)\right]=D_m\left[1-\exp \left(-2\sqrt{{\alpha }_m^2+\frac{i{L}^2\omega }{\nu }}h\right)\right],$$

$$(A_{1m}-A_{2m})\exp \left[-\sqrt{{\alpha }_m^2+\frac{i{L}^2\omega }{\nu }}(b-h)\right]={\mu }_r{D}_m\left[1+\exp \left(-2\sqrt{{\alpha }_m^2+\frac{i{L}^2\omega }{\nu }}h\right)\right].$$

$$B_0+C_0$$

$$C_m$$

$$-\frac{(-1{)}^m}{{\alpha }_m}\sqrt{\frac{i{L}^2\omega }{\nu }}(B_0-C_0)-\sum _{n=1}^{\infty }{C}_n{L}_{nm}\sqrt{{\gamma }_n^2+\frac{i{L}^2\omega }{\nu }}$$

$$=\frac{a}{2}\sqrt{{\alpha }_m^2+\frac{i{L}^2\omega }{\nu }}\left[A_{1m}-A_{2m}\exp \left(-2\sqrt{{\alpha }_m^2+\frac{i{L}^2\omega }{\nu }}(b-h)\right)\right].$$

$$i{\nabla }^2\hat{\Psi }=\frac{\nu }{L^2\omega }{\nabla }^4\hat{\Psi },$$

$${\hat{\Psi }}_I(y,z)=\sum _{n=1}^{\infty }\cos ({\alpha }_{n}z)f_n(y)+\sum _{n=1}^{\infty }{C}_n\sin ({\lambda }_{n}y){\varphi }_n({\lambda }_n;z),$$

$$f_n(y)$$

$$+A_{3n}\exp ({\alpha }_{n}y)+A_{4n}\exp [-{\alpha }_n(y+2(b-h))],$$

$${\varphi }_n({\lambda }_n;z)$$

Full text

1 Introduction

The potential drag reducing features of superhydrophobic surface (SHS) (Quéré, 2008) have renewed interest in the multiphase flow near patterned surfaces. To characterise the SHS, one needs to resolve both its micro-structure and multiphase features, which can be expensive at high Reynolds numbers (\Rey). Therefore, it would be advantageous to predict analytically the flow field near SHS. Pioneered by Philip’s study (1972a) of flow over alternating no-slip and no-shear boundary conditions, several aspects of SHS were analytically studied by Lauga & Stone (2003); Sbragaglia & Prosperetti (2007); Bazant & Vinogradova (2008); Crowdy (2010); Schönecker & Hardt (2013); Luca et al. (2018) and many other researchers.

The typical size of the SHS is often within the viscous sublayer (see Golovin et al. (2016) and references therein), therefore the convective terms in the Navier-Stokes equations are negligible as compared to the viscous terms, reducing the governing equations for the fluid to the Stokes equations. Moreover, the inertia terms are comparable to the convective terms, which further simplifies the problem into a quasi-steady Stokes problem or steady problem if the boundary condition is time-dependent (Luchini et al., 1991). Most analytical work on such problems have agreed on this assumption and adopted the steady Stokes equations as the governing equation.

Although the shear-free interface boundary condition has been widely applied (Philip, 1972a, b; Ng & Wang, 2009; Sbragaglia & Prosperetti, 2007; Crowdy, 2010; Schnitzer, 2017), there has been an emergence of liquid-infused surface in recent years: Rosenberg et al. (2016); Fu et al. (2017); Van Buren & Smits (2017); Keiser et al. (2017), which requires a more accurate estimation of the interfacial shear. Belyaev & Vinogradova (2010); Schönecker & Hardt (2013); Nizkaya et al. (2014) have considered the finite shear by determining a finite slip length, placing the interface at the crest of the corrugated surface and treating it as flat. The interface, however, can deform due to capillary effects, as considered in Sbragaglia & Prosperetti (2007) and Crowdy (2016, 2017) using a zero-shear approximation. Li et al. (2017) applied an approximate boundary condition to prescribe finite shear and considered interface deformation by placing the interface below the crest of the corrugated surface. Their solution is steady and accurate when the gas layer is thin. In this paper, the penetration effect is included similarly but with both phases solved simultaneously, and the solution being unsteady.

The flow may either be shear-driven (Wang, 2003; Ng & Wang, 2009; Kamrin et al., 2010; Schönecker et al., 2014; Li et al., 2017), or pressure-driven (Sbragaglia & Prosperetti, 2007; Belyaev & Vinogradova, 2010; Ng et al., 2010). For shear-driven flow, a shear rate is usually imposed at an infinite height, hence the solution is valid when the channel height is much greater than the geometric features on the surface. Pressure-driven flow can be decomposed into a Poiseuille flow and a perturbation due to the SHS for Stokes flow (Teo & Khoo, 2009), which works well when the shear rate on the multiphase interface is determined heuristically (Teo & Khoo, 2009; Schönecker & Hardt, 2013). If the second fluid inside the microstructures is included in the description, it could experience the same pressure gradient (Crowdy, 2017) or be driven by the interfacial stress (Ng et al., 2010), considering that the net cross-sectional mass flux in the second fluid is zero (Maynes et al., 2007). The influence of different geometric parameters of the grooves and interface condition is discussed in Li et al. (2017). In this paper, the flow is shear-driven by an oscillating velocity in either streamwise or spanwise direction, with the grooves aligned in the streamwise direction.

We study the effect of freestream unsteadiness on the multiphase flow inside the SHS grooves. We obtain analytical solutions for unsteady Stokes flow over an array of longitudinal grooves with two different fluids of arbitrary density and viscosity. The interfacial stress is finite and obtained as part of the solution. The solution is applicable to both air- and liquid-infused surfaces. §2 and §3 construct and solve the model problem analytically for longitudinal and transverse flows respectively. In §4, the analytical solution is validated against VOF simulation results. A parametric study of frequency, Reynolds number, and multiphase interface location is performed in §5. Concluding remarks are made in §6.

2 Longitudinal flow

We consider the incompressible and viscous flow over an array of grooves driven by an oscillating velocity. A schematic of the problem is shown in figure 1. We introduce the non-dimensional variables as:

[TABLE]

where the variables with the asterisks represent dimensional variables, $u$ is the streamwise velocity, $U_{0}$ is the velocity magnitude, $L$ is the half period of the groove, and $\omega$ is the frequency of the oscillating velocity. The groove geometry at the lower wall, $y=0$ is defined by its depth $b$ and half-groove width $a$. The multiphase interface is located at $y=-(b-h)$ and does not change with time. In this section, the velocity is along the direction of the grooves. The oscillating velocity is $u(H,z,t)=\cos(t)$. The groove is periodic in the span and symmetric at its centre.

The dimensionless governing equation is:

[TABLE]

where $L^{2}\omega/\nu$ is the non-dimensional Womersley number. The domain of interest is divided into three rectangular regions. Regions I and III represent the two types of fluid inside the grooves respectively and region II represents the fluid between $y=0$ and the oscillating plane at $y=H$.

For oscillatory flow, we assume the solution is of the form:

[TABLE]

Substituting (5) into (4) we obtain:

[TABLE]

By separation of variables $\hat{u}(y,z)=Y(y)Z(z)$, we obtain $\frac{Y^{\prime\prime}}{Y}+\frac{Z^{\prime\prime}}{Z}=\frac{iL^{2}\omega}{\nu}$.

We obtain the general solutions as a series of eigensolutions in each region independently and determine the constant coefficients by matching the boundary conditions between the neighbouring regions. The separable solutions are independent of time, therefore the boundary conditions can be written as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\mu_{r}=\mu_{A}/\mu_{B}$ is the viscosity ratio between the external fluid $A$ and the internal fluid $B$ (blue shaded area in figure 1), when $\mu_{r}=0$, the interface is shear-free. Fluid $A$ can penetrate into the groove.

2.1 Region I

The boundary conditions are no-slip at the wall ($z=a$) and symmetry at the groove centre ($z=0$). In terms of $\hat{u}(y,z)=Y(y)Z(z)$,

[TABLE]

From (11) and (12), we get $Z_{n}(z)=\cos(\alpha_{n}z)$, where $\alpha_{n}=(n-\frac{1}{2})\frac{\pi}{a}$. This implies $\frac{Y^{\prime\prime}}{Y}=\frac{iL^{2}\omega}{\nu}+\alpha_{n}^{2}$. The solution for $Y$ is therefore:

[TABLE]

A separable solution for $\hat{u}_{\mathrm{I}}$ is:

[TABLE]

The coefficients $A_{1n}$ and $A_{2n}$ will be obtained by matching adjoining regions below.

2.2 Region II

The boundary conditions in region II are symmetric at the centre ($z=0$) and between grooves ($z=1$). Also, at the fixed height $y=H$, the slip velocity oscillates with time.

[TABLE]

The solution for $Z$ must satisfy equations (14) and (15), i.e.: $Z_{n}(z)=\cos(\gamma_{n}z)$, where $\gamma_{n}=n\pi$. Hence the solutions for $Y$ must satisfy: $\frac{Y^{\prime\prime}}{Y}=\frac{iL^{2}\omega}{\nu}+\gamma_{n}^{2}$, for any integer $n\geqslant 0$. Therefore the eigensolution for $Y$ is:

[TABLE]

Therefore, the general solution in region II is:

[TABLE]

The inhomogeneous behaviour of the flow in the spanwise direction caused by the groove is represented by $\cos(\gamma_{n}z)$, which should diminish as the distance to the groove increases; therefore the coefficients $B_{n}=0$. Integrating the boundary condition (16) from [math] to $1$ yields:

[TABLE]

2.3 Region III

The general solution for region III that satisfies the no-slip condition on the walls and symmetry is:

[TABLE]

2.4 Solving for coefficients in the general solutions

2.4.1 Between region I and region III

The boundary at the interface satisfies equations (9) and (10). Multiplying (9) and (10) by $\cos(\alpha_{m}z)$ and integrating from [math] to $a$ respectively, yields:

[TABLE]

and

[TABLE]

2.4.2 Between region I and region II

Next, we determine the unknown constants $C_{0}$, $C_{n}$ and $D_{n}$ by applying the boundary condition between the two regions. Integrating (7) from 0 to 1,

[TABLE]

Multiplying (7) by $\cos(\gamma_{m}z)$ and integrating from 0 to 1,

[TABLE]

where $L_{mn}=\int_{0}^{a}\cos(\alpha_{n}z)\cos(\gamma_{m}z)dz$. Multiplying (8) by $\cos(\alpha_{m}z)$ and integrating from [math] to $a$,

[TABLE]

Finally, $C_{m}$ is truncated to $M$ terms and $D_{n}$ is truncated to $N=\mathrm{Int}[aM]$ terms. By solving equations (19), and (21) to (25) as a system of equations in MATLAB, we can obtain all the unknown coefficients in the solution. Table 1 shows the convergence of $B_{0}+C_{0}$ with $M$. $B_{0}+C_{0}$ is accurate up to the fourth digit when $M$ is greater than $25$.

3 Transverse flow

For transverse flow, the freestream oscillating velocity is in the $z$ direction: $w(H,z,t)=\cos(t)$. A schematic of the problem is shown in figure 2. The dimensionless governing equation is:

[TABLE]

with the streamfunction ($\hat{\Psi}$) to replace $\hat{w}=\frac{\partial\hat{\Psi}}{\partial y}$ and $\hat{v}=-\frac{\partial\hat{\Psi}}{\partial z}$.

3.1 Region I

$\hat{\Psi}$ and its normal derivatives are zero at the walls and symmetric at the groove centre ($z=0$). A separable solution for $\hat{\Psi}_{\mathrm{I}}$ is:

[TABLE]

where $\lambda_{n}=n\pi/(b-h)$, and

[TABLE]

[TABLE]

At the side wall

[TABLE]

Multiplying this by $\sin(\lambda_{m}y)$ and integrating from $-(b-h)$ to [math] yields:

[TABLE]

where

[TABLE]

The location of the interface is at $y=-(b-h)$ such that $\frac{\partial\hat{\Psi}_{I}}{\partial z}=0$, yielding,

[TABLE]

3.2 Region II

The boundary conditions in region II are symmetric at the centre ($z=0$) and periodic between grooves ($z=1$). At the fixed height $y=H$, the slip velocity oscillates with time. The boundary condition at $y=H$ is:

[TABLE]

The general solution in region II is:

[TABLE]

Integrating the top boundary condition (33) from [math] to $1$ gives:

[TABLE]

3.3 Region III

The general solution for region III that satisfies $\hat{\Psi}=0$ on the walls and symmetry at the centre is:

[TABLE]

where $\beta_{n}=n\pi/h$,

[TABLE]

and $\phi_{n}(\beta_{n};z)$ is the same as (29) with different eigenvalues $\beta_{n}$. $\hat{\Psi}_{\mathrm{III}}$ satisfies:

[TABLE]

at the bottom wall, and

[TABLE]

at the side wall $z=a$.

Equation (38) yields:

[TABLE]

Multiplying equation (39) by $\cos(\alpha_{m}z)$ and integrating from [math] to $a$:

[TABLE]

where

[TABLE]

Multiplying equation (40) by $\sin(\beta_{m}(y+b-h))$ and integrating from $-b$ to $-(b-h)$ gives

[TABLE]

where

[TABLE]

The location of the interface is fixed at $y=-(b-h)$, yielding,

[TABLE]

3.4 Solving for coefficients in the general solutions

3.4.1 Between region I and region III

The boundary conditions between region I and region III are:

[TABLE]

[TABLE]

Multiplying (46) and (47) by $\cos(\alpha_{m}z)$ and integrating from [math] to $a$, respectively:

[TABLE]

and

[TABLE]

where

[TABLE]

3.4.2 Between region I and region II

The boundary condition between regions I and II:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Integrating equation (51) from [math] to $1$ gives:

[TABLE]

Multiplying by $\cos(\gamma_{m}z)$ and integrating yields:

[TABLE]

Integrating equation (52) from [math] to $1$ gives:

[TABLE]

Multiplying by $\cos(\gamma_{m}z)$ and integrate equation (52) from [math] to $1$ yields:

[TABLE]

where

[TABLE]

[TABLE]

Multiplying equation (53) and equation (54) by $\cos(\alpha_{m}z)$ and integrating from [math] to $a$ yields:

[TABLE]

and

[TABLE]

respectively.

Finally, $E_{m}$ and $F_{m}$ are truncated to $M$ terms, $A_{n}$ and $B_{n}$ are truncated to $N=\mathrm{Int}[aM]$ terms, $C_{n}$ and $D_{n}$ are truncated to $K=\mathrm{Int}[bM]$ terms and solved from a system of equations constructed with equations (31), (32), (35), (41), (42), (44), (45), (48), (49), (55) to (58), (61), (62) in MATLAB. Table 2 shows the convergence of averaged velocity magnitude $\Re\left(\sqrt{\frac{iL^{2}\omega}{\nu}}(E_{0}-F_{0})\right)$ at $y=0$ with $M$. $\Re\left(\sqrt{\frac{iL^{2}\omega}{\nu}}(E_{0}-F_{0})\right)$ is accurate up to the fourth digit when $M$ is greater than $30$.

4 Validation

The solution is first validated against DNS data using the numerical methodology of Mahesh et al. (2004), implemented with VOF given in Li et al. (2017). Then the lower bound of the solution, representing a steady flow result is validated using the analytical solution in Philip (1972b), Schönecker et al. (2014), and the simulation result in Fu et al. (2017).

4.1 Validation of the multiphase unsteady solution

The multiphase unsteady solution is validated by DNS with VOF. The results are expressed in the form of the transfer equation defined as:

[TABLE]

where the input is considered as the energy of the forcing, and the output is the energy at $y=0$.

The transfer function shows good agreement between VOF simulation with flat interface ($\sigma\to\infty$) and the analytical solution in figure 3. In all figures, $L^{2}\omega/\nu$ varies from $0.1$ to $1.0$, $L=1$, $\nu=1$, $H=2.5$, $h/b=1.0$, but $\mu_{r}=0.02$ in figure 3(a, b), and $\mu_{r}=0.37$ in figure 3(c, d), respectively.

4.2 Validation of the multiphase steady solution

Schönecker et al. (2014) considered the flow over an array of grooves being filled with a secondary immiscible fluid driven by a constant shear stress $\tau_{\infty}$ as $y\to\infty$ and derived analytical expressions for effective slip length, with respect to the viscosity ratio of the two fluids. Here, $b_{\mathrm{eff}}\equiv\overline{u}_{\mathrm{II}}/\overline{\partial u_{\mathrm{II}}/\partial y}$, where the overline ‘$\overline{()}$’ represents the averaged quantity across time and span, with both values evaluated at $y=0$. This solution can serve as a lower bound of the current unsteady solution where $\omega\to 0$ and $y\to\infty$. This asymptotic trend is confirmed by first comparing two sets of analytical solutions with $L^{2}\omega/\nu=10^{-4}$, $H=5$; and $L^{2}\omega/\nu=10^{-5}$, $H=10$, with the same $a=0.5$, $b=1.0$, $h=b$, and $\mu_{r}=0.02$. For the two sets of parameters $b_{\mathrm{eff\parallel}}=0.2101$, agreeing up to the $10th$ digit. $b_{\mathrm{eff\perp}}=0.1032$, agreeing up to the $10th$ digit as well. Therefore, $L^{2}\omega/\nu=10^{-4}$, and $H=5$ are used to compare to the steady flow solutions.

When a small amount of penetration is considered, the data shows good comparison to the simulation data of Fu et al. (2017) at $\mu_{r}<1$.

5 Results and discussion

5.1 Instantaneous flow field

When the grooves are in parallel to the flow direction, the streamwise velocity changes monotonically along the wall normal direction. Its flow pattern has been shown in multiple literature (to name a few: Schönecker & Hardt (2013); Li et al. (2017)). When the grooves are transverse to the flow direction, the flow field inside the grooves behaves similarly to a lid-driven cavity flow. Figure 5 shows the instantaneous streamlines inside the grooves of different depths $b$ with different penetration rates $1-h/b$. The external oscillatory boundary condition is all applied at the same $H=5.40$. The penetration rate represents the area fraction of fluid $A$ inside the groove to the area of the groove $(w(b-h)/wb)$. On one hand, as the groove becomes deeper, multiple vortices form inside, which is similar to the eddy structure in lid-driven cavities with increasing height (Shankar, 2007). On the other hand, fluid $B$ in the blue shaded area generates vortices separately from fluid $A$ (white area). Multiple vortices form within fluid $B$ when its depths gets deep as well.

5.2 Parametric study

In this section, the effect of frequency, domain height, penetration rate, and viscosity ratio are studied using the analytic solution. All the following cases have a coverage ratio of $a=0.875$ otherwise stated. The parameters for each case are listed in table 3. The frequency range is selected to be equivalent to the nondimensional frequency range inside a turbulent channel flow at the given $Re_{\tau}$ listed on the last column. The location of the forcing $H$ is equivalent to $y^{+}=5$ for the $Re_{\tau}$ given. The viscosity ratios are $0.02,\ 0.37,\ 30.00$ to represent air, heptane, and Dupont Krytox to water, respectively (Rosenberg et al., 2016).

The transfer function of the cases from table 3 are shown in figure 6(a-c). The transfer function $\mathcal{H}\left(\frac{iL^{2}\omega}{\nu}\right)$ is defined as the ratio between the energy at $y=0$ and the energy of the forcing at $y=H$:

[TABLE]

Overall, as the oscillation frequency increases, the transfer function decreases. As the representative Reynolds number increases, the transfer function increases. However, note that as Reynolds number increases, the fluctuations of the forcing would also be closer to the surface. The shaded areas correspond to the condition $\mathcal{H}(L^{2}\omega/\nu,\mu_{r}\geq 1)$ to distinguish between high viscosity fluid and low viscosity fluid inside the groove. Note that the trend of $\mathcal{H}$ with penetration rate is different between shaded and non-shaded regions. When $\mu_{r}>1$, as the penetration rate increases, $\mathcal{H}$ increases, and the penetration of the external fluid improves the transfer of forcing energy. Conversely, when $\mu_{r}<1$, as the penetration rate increases, $\mathcal{H}$ decreases, and the penetration of the external fluid suppresses the transfer of the forcing energy. Neither of these two conditions, however, cross the border of the $\mu_{r}=1$ curve. Also, air performs the best over the other two fluids, which is intuitively reasonable, considering that the interface does not break or drain based on the assumption made in this paper. The transfer function for transverse flow (figure 6 d-f) behaves the same way as in longitudinal flow.

The root-mean-square (rms) of the unsteady analytical effective slip length is determined by:

[TABLE]

for longitudinal flow; and

[TABLE]

for transverse flow (Luchini et al., 1991; Maali & Bhushan, 2012; Cottin-Bizonne et al., 2004), where $\int_{0}^{1}\frac{\partial u_{\mathrm{II}}}{\partial y}(\infty,z)\mathrm{d}z=1$. In viscous units, $b_{0}^{+}=b_{\mathrm{eff}}Lu_{\tau}/\nu$, where $u_{\tau}$ is the friction velocity.

The unsteady effective slip length is plotted in figure 7(a-c) for longitudinal flow and figure 7(d-f) for transverse flow. Zero penetration of the external fluid with air pockets provide the highest slip length. Transverse effective slip length is more sensitive to the penetration of the external fluid: a $15\%$ penetration would result in an effective slip length overlapping with the fully-wetted grooved cases.

Comparing the longitudinal effective slip length to the transverse values in figure 8 using the ratio $b_{\mathrm{eff\parallel}}/b_{\mathrm{eff\perp}}$, it is found that the effective slip length for longitudinal flows is more than $1.5$ times of that of transverse flows with the same coverage ratio. The ratio drops at high frequency. A $15\%$ penetration rate could cause more decrease of the slip length in transverse direction, resulting a higher ratio. This agrees with the prior observation. But $b_{\mathrm{eff\parallel}}/b_{\mathrm{eff\perp}}$ is not a monotonic function of the penetration rate. For the combination of zero penetration and air-water interface, $b_{\mathrm{eff\parallel}}/b_{\mathrm{eff\perp}}=2$ at low frequencies. This agrees with the finding in Lauga & Stone (2003) where no-shear strips in steady pipe flow were considered.

The shear stress at the plane $y=0$ is evaluated by

[TABLE]

for longitudinal flows, and

[TABLE]

for transverse flows. Figure 9 shows the shear stress varying across the frequency range. Again, the combination of zero penetration and air pocket produces the least shear stress for both longitudinal and transverse flows. The effect of penetration of the external fluid has a larger impact on transverse flow than longitudinal flow. The trend observed in shear stress is consistent with those found in transfer function and effective slip length. As the frequency increases, the shear stress decreases, because the flow near the grooved surface is approaching a stationary state.

6 Conclusion

Analytical solutions of the multiphase Stokes flow generated by oscillatory velocity at a finite-height over a grooved surface have been derived. In the limit of the steady state, the analytical results agree well with the steady solution in Schönecker et al. (2014) and Fu et al. (2017). Numerical simulations using the VOF methodology were performed to validate the unsteady multiphase analytical solution. The analytical solution compares well with the DNS data. A transfer function based on the energy inside the grooved surface and the energy of the oscillatory input velocity was obtained. The analytical solution was parameterized by $\omega L^{2}/\nu$, a representative Reynolds number $Re_{\tau}$, and location of the multiphase interface $h/b$. We see that: (i) large values of $\omega L^{2}/\nu$ lower the transfer function, (ii) increasing $Re_{\tau}$ decreases the forcing height and therefore increases the transfer function, and (iii) When $\mu_{r}<1$, penetration of the external fluid suppresses the transfer function; when $\mu_{r}>1$, penetration of the external fluid increases it.

The variation of effective slip length and shear stress over the grooved plane show a consistent trend as the transfer function. The grooves show greater impedance to transverse flow than longitudinal flow, as is reported by other researchers (Luchini et al., 1991; Choi et al., 1993). Such an impedance is further strengthened when the external fluid penetrates into the groove. The ratio between the rms value of the longitudinal effective slip length $b_{\mathrm{eff\parallel}}$ and that of the transverse effective slip length $b_{\mathrm{eff\perp}}$ is optimal when $\mu_{r}=0.02$, $1-h/b=0.15$, meaning that the surface is providing a large slip while impeding the cross-flow motion the least. When $\mu_{r}=0.02$, $1-h/b=0$, $b_{\mathrm{eff\parallel}}/b_{\mathrm{eff\perp}}$ is $2$, which agrees with the finding of Lauga & Stone (2003).

The dimensionless parameter $\omega L^{2}/\nu$ can be related to the $\Rey$ in a turbulent flow as: $\frac{\omega L^{2}}{\nu}=\omega^{+}Re_{\tau}^{2}\left(\frac{L}{\delta}\right)^{2}$, which is approximately the order of $O(10)$ for a moderate $\Rey$ considering $L/\delta\sim O(10^{-2})$. As $\Rey$ increases, at the same $\omega^{+}$, both streamwise and spanwise velocity fluctuation characterised as the oscillating velocity in the analytical solution, approach the surface, the energy of the fluctuation will influence the interface stronger in both directions.

Acknowledgement

This work was supported by the United States Office of Naval Research (ONR) MURI (Multidisciplinary University Research Initiatives) program under Grant N00014-12-1-0874 managed by Dr. Ki-Han Kim. Computing resources were provided by the Minnesota Supercomputing Institute (MSI).