3.1 Square array: general features

We begin our analysis by examining the square array case. The periodic computational domain of side $L$ contains in this case a single sphere, and the area fraction can be simply calculated as $\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x03C0}a^{2}/L^{2}$ . We are particularly interested in the area-averaged streamwise velocity $\langle u\rangle$ , defined as

and how this quantity changes in the direction normal to the monolayer. In (3.1) the integral is extended over the region $-L/2\leqslant x\leqslant L/2$ , $-L/2\leqslant y\leqslant L/2$ .

In figure 3(a), the area-averaged velocity is plotted as a function of the coordinate $z$ normal to the monolayer for a fixed gap size, $d=3a$ . The velocity profile is approximately linear in two limits: when $\unicode[STIX]{x1D719}_{s}\ll 1$ and when $\unicode[STIX]{x1D719}_{s}$ is close to the maximum packing fraction $\unicode[STIX]{x1D719}_{s,max}=\unicode[STIX]{x03C0}/4\simeq 0.78$ . For intermediate values of $\unicode[STIX]{x1D719}_{s}$ , the curvature of the velocity profile has a maximum. This trend can be understood by applying the averaging operator defined in (3.1) to the streamwise component of the fluid momentum equation in the Stokes flow limit. The resulting averaged equation reads

where $\unicode[STIX]{x0394}\overline{p}(z)$ is the difference between the pressure at the plane $x=-L/2$ , averaged over the line $-L/2\leqslant y\leqslant L/2$ , and the corresponding average pressure at the plane $x=L/2$ . Expression (3.2) shows that the curvature of the average velocity profile can be neglected when the streamwise pressure drop occurring over a distance $L$ is negligible.

According to the Stokes equation, the pressure disturbance induced by the particle determines the curvature of the velocity profile. When $\unicode[STIX]{x1D719}_{s}\ll 1$ , $\unicode[STIX]{x0394}\overline{p}\simeq 0$ because the pressure disturbance set up by each sphere evaluated at the boundaries of the computational domain is small. Indeed, we will shortly see that in the dilute limit the pressure disturbance induced by a particle looks like a pressure dipole, and in the dilute limit we therefore expect $\unicode[STIX]{x0394}\overline{p}=O(\unicode[STIX]{x1D707}aU/L^{2})$ (Batchelor Reference Batchelor2000), as for a single sphere in uniform flow. The pressure disturbance $\unicode[STIX]{x0394}\overline{p}$ decreases linearly with $\unicode[STIX]{x1D719}_{s}$ as $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ , since $\unicode[STIX]{x1D719}_{s}\propto 1/L^{2}$ . On the other hand, when the monolayer is near maximum packing, the monolayer behaves almost as a flat wall. In this case, $\unicode[STIX]{x0394}\overline{p}\simeq 0$ . For intermediate values of $\unicode[STIX]{x1D719}_{s}$ , the inter-particle distance is simultaneously sufficiently small for the pressure disturbance produced by each particle on the boundaries of the computational domain to be significant, and sufficiently large for each particle not to block significantly the flow velocity incident on the other particles. As a consequence, the curvature has a maximum for intermediate values of $\unicode[STIX]{x1D719}_{s}$ .

Figure 3.

(a) Normalised area-averaged streamwise velocity versus coordinate normal to the monolayer for $d=3a$ . (b,c) Normalised area-averaged streamwise velocity along lines passing through the centre of each sphere (b) or through the midpoint between adjacent spheres (c).

To characterise how the velocity profile changes in the plane of the monolayer, we show in figure 3(b,c) the profiles of $u$ along lines perpendicular to the plane of the monolayer and passing through the sphere centre, $(x=0,y=0)$ , and through the midpoint between two spheres, $(x=L/2,y=0)$ , respectively. Because of the fore–aft symmetry of the Stokes equation, the velocity disturbances generated by two adjacent spheres cancel out at the midpoint $(x=L/2,y=0)$ . As a consequence, the velocity profiles corresponding to the midpoint are linear for any value of $\unicode[STIX]{x1D719}_{s}$ . As $\unicode[STIX]{x1D719}_{s}$ is reduced, the slope of the velocity profiles decreases. This trend gives a larger slip velocity. The velocity profiles that correspond to the sphere centres (figure 3 b) display significant curvature for any $\unicode[STIX]{x1D719}_{s}$ , except perhaps at the highest value of the area fraction.

The error induced by approximating a liquid interface covered by a packed monolayer of particles as a no-slip surface is due to two sources. First of all, a packed monolayer contains free-slip surfaces even at maximum packing. Secondly, owing to the convexity of the spheres, the packed monolayer is not flat. A more accurate interpretation is considering the monolayer as a collection of bluff solid protuberances over the no-shear plane $z=0$ . Each protuberance will locally block the flow and therefore produce a significant pressure disturbance.

Figure 4.

Normalised pressure (a,c) and velocity (b,d) in the plane $y=0$ for $\unicode[STIX]{x1D719}_{s}=12.57\,\%$ (a,b) and $\unicode[STIX]{x1D719}_{s}=62.05\,\%$ (c,d). The gap size is $d=4a$ .

To illustrate the spatial extent and magnitude of the pressure disturbance set up by each sphere in the monolayer, we show in figure 4(a,c) isocontours of the normalised pressure in the plane parallel to the mean flow and perpendicular to the monolayer, for two selected values of the area fraction, $\unicode[STIX]{x1D719}_{s}=12.57\,\%$ and $\unicode[STIX]{x1D719}_{s}=62.05\,\%$ . For $\unicode[STIX]{x1D719}_{s}=12.57\,\%$ , the pressure distribution bears a signature of the fore–aft symmetric pressure dipole characteristic of uniform flow past a sphere (Batchelor Reference Batchelor2000). As the area fraction increases (figure 4 b), the two pressure ‘lobes’ seen in figure 4(a) move upwards, and occupy a smaller region near the top apex of the sphere. The streamwise pressure drop caused by the sphere is smaller for $\unicode[STIX]{x1D719}_{s}=62.05\,\%$ than for $\unicode[STIX]{x1D719}_{s}=12.57\,\%$ . This observation supports our suggestion that $\unicode[STIX]{x0394}\overline{p}$ decreases as $\unicode[STIX]{x1D719}_{s}$ approaches the maximum packing fraction.

Figure 4(b,d) shows the velocity fields corresponding to figure 4(a,c), respectively. Flow recirculation regions do not seem to occur between the spheres. Flow recirculation between the spheres was evidenced in a simulation of pressure-driven slip flow over a porous medium interface, where the simulations were carried out using the same numerical method as used here (Liu & Prosperetti Reference Liu and Prosperetti2011). Flow recirculation is instead not apparent in the simulation results reported by Danov et al. (Reference Danov, Aust, Durst and Lange1995), who examined a single sphere straddling a fluid interface for a range of contact angles. One could expect that, owing to the presence of a re-entrant ‘wedge region’ between the particle surface and the free interface, recirculation would occur for contact angles for which the sphere is mostly immersed in the liquid.

Figure 5.

Normalised slip velocity versus (a) area fraction and (b) gap size.

The slip velocity is plotted as a function of $\unicode[STIX]{x1D719}_{s}$ for different values of the gap size $d$ in figure 5(a), and as a function of $d$ for different values of $\unicode[STIX]{x1D719}_{s}$ in figure 5(b). Following other authors (Ybert et al. Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007; Ng & Wang Reference Ng and Wang2009; Liu & Prosperetti Reference Liu and Prosperetti2011), we define the slip velocity $\langle u\rangle _{s}$ as the value of $\langle u\rangle$ at $z=a$ . As expected $\langle u\rangle$ tends to a uniform velocity of magnitude $U$ as $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ . The rate at which this limit is approached as $\unicode[STIX]{x1D719}_{s}$ changes depends on the gap size. As a consequence, the value $\unicode[STIX]{x1D719}_{s}$ for which $\langle u\rangle$ is a significant fraction of $U$ becomes smaller as $d$ increases. For instance, when $d/a=20$ , $\langle u\rangle /U\simeq 0.5$ when $\unicode[STIX]{x1D719}_{s}\simeq 1.5\,\%$ . When $d/a=1.5$ , $\langle u\rangle /U\simeq 0.5$ when $\unicode[STIX]{x1D719}_{s}\simeq 30\,\%$ . Since smaller gaps are associated with larger velocity gradients, this dependence on $d$ indicates a dependence of $\langle u\rangle _{s}$ on the shear rate.

The dependence of $\langle u\rangle _{s}$ on $d$ is nonlinear, but approximates a linear relation in the limit $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ (figure 5 b). The following argument shows that in the dilute limit the slope of the $\langle u\rangle _{s}$ versus $d$ curve is proportional to the gap size and inversely proportional to the slip length. Approximating the flow past the monolayer as a Couette flow past a flat partial slip surface, the macroscopic shear rate is $\langle \dot{\unicode[STIX]{x1D6FE}}\rangle _{s}\simeq (U-\langle u\rangle _{s})/d$ . Using this value in the Navier slip boundary condition (1.1) gives $\langle u\rangle _{s}/U\simeq 1/[(d/\unicode[STIX]{x1D706})+1]$ . For $\unicode[STIX]{x1D719}_{s}\ll 1$ , $\unicode[STIX]{x1D706}$ becomes much larger than $d$ and therefore $\langle u\rangle _{s}/U=1-(d/\unicode[STIX]{x1D706})$ with a small $O(d/\unicode[STIX]{x1D706})$ error.

3.2 Scaling of slip length for square array

Figure 6.

(a) Normalised slip length versus area fraction and normalised gap size. The plots in panels (b) and (c) are projections of the surface plot of panel (a) onto the $\unicode[STIX]{x1D706}$ – $d$ and $\unicode[STIX]{x1D706}$ – $\unicode[STIX]{x1D719}_{s}$ planes, respectively. The inset in panel (c) shows the $\unicode[STIX]{x1D706}$ – $\unicode[STIX]{x1D719}_{s}$ relationship in log–log scale.

Figure 6 shows the slip length $\unicode[STIX]{x1D706}$ as calculated from the definition (1.1). Consistently with the calculation of the slip velocity, the bulk shear rate at the interface $\langle \dot{\unicode[STIX]{x1D6FE}}\rangle _{s}=\text{d}\langle u\rangle /\text{d}z$ is evaluated at $z=a$ . In the surface plot of figure 6(a), $\unicode[STIX]{x1D706}$ is plotted as a function of both $\unicode[STIX]{x1D719}_{s}$ and $d/a$ . Projections of figure 6 onto the $\unicode[STIX]{x1D706}$ – $d$ and $\unicode[STIX]{x1D706}$ – $\unicode[STIX]{x1D719}_{s}$ planes are shown in linear–linear plots in figure 6(b,c), respectively. The inset of figure 6(c) shows the $\unicode[STIX]{x1D706}$ – $\unicode[STIX]{x1D719}_{s}$ relation in log–log scale.

The slip length is seen to increase for increasing gap sizes, eventually saturating to an asymptotic value. A gap size $d=10a$ already gives a value of $\unicode[STIX]{x1D706}$ close to the asymptotic limit.

In figure 6(c), the values of $\unicode[STIX]{x1D706}$ for different values of $d$ are seen to lie practically on the same curve, suggesting similar scaling laws for different gap sizes. The inset shows that the relation $\unicode[STIX]{x1D706}$ – $\unicode[STIX]{x1D719}_{s}$ follows approximately a power law. To guess possible scaling exponents, we have examined the literature on flows over microstructured superhydrophobic surfaces. Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007) developed a comprehensive theory for the dependence of the slip length on the area fraction for superhydrophobic surfaces. The theory was compared against literature data. Data for unconfined shear flow past a superhydrophobic surface composed of vertical pillars of circular or square cross-section could be fitted with good accuracy by a correlation of the form $\unicode[STIX]{x1D706}/L=(A_{1}/\sqrt{\unicode[STIX]{x1D719}_{s}})-B_{1}$ , where $L$ is the distance between the pillars, and $A_{1}$ and $B_{1}$ are constants. Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007) proposed $A_{1}=0.325$ and $B_{1}=0.44$ (the analytical expressions proposed by Davis & Lauga (Reference Davis and Lauga2010) give coefficient values very close to those indicated by Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007)). In the dilute limit, i.e. in our case for $\unicode[STIX]{x1D719}_{s}\ll (A_{1}/B_{1})^{2}\simeq 0.54$ , the correlation above reduces to $\unicode[STIX]{x1D706}/L\sim 1/\sqrt{\unicode[STIX]{x1D719}_{s}}$ , which is equivalent to $\unicode[STIX]{x1D706}/a\sim 1/\unicode[STIX]{x1D719}_{s}$ . This scaling can be understood from the fact that in the dilute limit the ratio of the hydrodynamic force on each pillar to the slip velocity is expected to be practically independent of $\unicode[STIX]{x1D719}_{s}$ . The tangential stress on the monolayer due to the bulk flow is proportional to the ratio of the hydrodynamic force and $L^{2}$ . Since the tangential stress for a Newtonian fluid is proportional to the macroscopic shear rate, then (1.1) yields $\unicode[STIX]{x1D706}\propto a/\unicode[STIX]{x1D719}_{s}$ for $\unicode[STIX]{x1D719}_{s}\ll 1$ .

Figure 7.

(a) Normalised slip length versus $\unicode[STIX]{x1D719}_{s}^{-1}$ . The continuous line is (3.2), developed by Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007) for flow past a superhydrophobic surface. The dashed line is (3.3), where (3.2) has been rescaled by the ratio of the Stokes drag coefficient for a sphere to that for a thin disk. (b) Normalised slip length versus $(1-\unicode[STIX]{x1D719}_{s})^{2}/\sqrt{\unicode[STIX]{x1D719}_{s}}$ showing the dense limit $\unicode[STIX]{x1D719}_{s}=35\,\%{-}70\,\%$ for $d=4a$ .

The argument above is expected to hold independently of the specific geometry of the solid object in contact with the fluid interface. Therefore, it should be possible to apply the scaling proposed by Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007) to our case. This expectation is confirmed in figure 7(a), where $\unicode[STIX]{x1D706}/a$ is plotted as a function of $1/\unicode[STIX]{x1D719}_{s}$ for values of $\unicode[STIX]{x1D719}_{s}$ corresponding to a relatively dilute monolayer. A linear correlation fits the data remarkably well. From figure 7(a) it appears that a linear scaling holds for any value of $d$ . The log–log plot in the inset of figure 6(b) shows that the power-law exponent has a small dependence on $d$ for $d\leqslant 4a$ . However, differences between the exponents corresponding to different values of $d$ are approximately $10\,\%$ of the $d=20a$ case, and therefore not noticeable if the data are plotted as in figure 7(a).

While the functional form proposed by Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007) does fit our data well, the prefactors are different. The function $\unicode[STIX]{x1D706}/a=(A_{1}\sqrt{\unicode[STIX]{x03C0}}/\unicode[STIX]{x1D719}_{s})-(B_{1}\sqrt{\unicode[STIX]{x03C0}}/\unicode[STIX]{x1D719}_{s}^{1/2})$ , obtained using the definition $\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x03C0}a^{2}/L^{2}$ , is plotted as a continuous line in figure 7(a), using the coefficient values for $A_{1}$ and $B_{1}$ suggested by Ybert et al. This function is seen to overpredict the magnitude of the slip length computed in our simulation.

There is a simple explanation for the difference between our result and that of Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007). The flow configuration considered by those authors can be interpreted as shear flow past a monolayer of infinitely thin disks (representing the top surfaces of the pillars composing the superhydrophobic surface). In our case the spheres protrude significantly into the liquid. Several studies in the context of superhydrophobic surfaces have shown that a larger protrusion produces smaller slip lengths (Sbragaglia & Prosperetti Reference Sbragaglia and Prosperetti2007; Ng & Wang Reference Ng and Wang2009; Kumar, Datta & Kalyanasundaram Reference Kumar, Datta and Kalyanasundaram2016; Shelley et al. Reference Shelley, Smith, Hibbins, Sambles and Horsley2016). Therefore, a smaller slip length in the current case of a sphere monolayer is not unexpected.

We can attempt a simple correction to (3.2) to account for the finite protrusion of the spheres into the flow. For a given free-stream velocity, the hydrodynamic drag on an infinitely thin disk in a uniform flow is $9\unicode[STIX]{x03C0}/16$ times smaller than the drag on a sphere having the same radius. Assuming that tangential stress due to the bulk flow is approximately proportional to the slip velocity (we will confirm this hypothesis in § 3.3), it should be expected that, for sufficiently dilute systems, the slip length for flow past a monolayer of spheres embedded in a gas–liquid interface is a factor of $9\unicode[STIX]{x03C0}/16$ smaller than that predicted by (3.2). To verify this approximation, we plot the function

as a thick dashed line in figure 7(a). Here $A_{2}=(16\unicode[STIX]{x03C0}\sqrt{\unicode[STIX]{x03C0}}/9)A_{1}$ and $B_{2}=(16\unicode[STIX]{x03C0}\sqrt{\unicode[STIX]{x03C0}}/9)B_{1}$ . The values given by this corrected expression are remarkably close to our data, suggesting that the difference between the data of Ybert et al. (Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007) and ours can be mainly attributed to the larger drag produced by the protruding spheres on the flowing liquid as opposed to the flat top surfaces of the pillars.

Figure 8.

Tangential stress due to the bulk flow normalised by (a) the wall velocity $U$ or (b) the slip velocity $\langle u\rangle _{s}$ as a function of the solid area fraction. The inset in panel (b) shows a zoom for small values of $\unicode[STIX]{x1D719}_{s}$ of the effective friction coefficient $\unicode[STIX]{x1D70F}a/(\unicode[STIX]{x1D707}\langle u\rangle _{s}\unicode[STIX]{x1D719}_{s})$ .

The scaling law discussed above holds for a relatively dilute limit (the data points clearly discernible in figure 7(a) correspond to $\unicode[STIX]{x1D719}_{s}\leqslant 4\,\%$ , while these corresponding to the dense limits are clustered near the origin and are barely visible). In the dense limit, the slip velocity results from a fraction $1-\unicode[STIX]{x1D719}_{s}$ of the plane $z=a$ occupied by the gas–liquid interfaces. Since the slip velocity in these regions is of the order of $\langle \dot{\unicode[STIX]{x1D6FE}}\rangle _{s}L(1-\unicode[STIX]{x1D719}_{s})$ (Ybert et al. Reference Ybert, Barentin, Cottin-Bizonne, Joseph and Bocquet2007) and $L\propto \unicode[STIX]{x1D719}_{s}^{-1/2}$ , we expect $\unicode[STIX]{x1D706}\sim a(1-\unicode[STIX]{x1D719}_{s})^{2}/\sqrt{\unicode[STIX]{x1D719}_{s}}$ . Figure 7(b) shows that the simulation data follow this scaling law with very good accuracy.

3.3 Tangential stress due to the bulk flow: square array

In many situations of practical interest it would be useful to estimate the drag per unit area on the particles due to the motion of the bulk fluid (or ‘subphase’, as it is often called in the surfactants literature). We call the hydrodynamic drag on the particles per unit area the tangential stress due to the bulk flow, and denote this quantity by $\unicode[STIX]{x1D70F}$ . In the current paper, $\unicode[STIX]{x1D70F}$ for a square array is calculated as $\unicode[STIX]{x1D70F}=F_{x}/L^{2}$ , where $F_{x}$ is the $x$ component of the hydrodynamic force acting on each sphere in the monolayer.

In figure 8(a), $\unicode[STIX]{x1D70F}$ is normalised by using $U$ as velocity scale. With this normalisation, the normalised value of $\unicode[STIX]{x1D70F}$ has a relatively strong dependence on both $\unicode[STIX]{x1D719}_{s}$ and $d$ . However, we note that the wall velocity is characteristic of the fluid velocity ‘seen’ by the particles only in the extremely dilute limit for which the velocity profile is almost a plug flow. An improved parametrisation of the shear stress and pressure exerted on each particle by the bulk flow is expected to be given by the slip velocity $\langle u\rangle _{s}$ . In figure 8(b), $\unicode[STIX]{x1D70F}$ is normalised by $\langle u\rangle _{s}$ . Using the slip velocity in the normalisation causes the data to collapse onto a single curve for any value of  $d$ .

The fact that $\unicode[STIX]{x1D70F}$ is practically proportional to the slip velocity suggests the introduction of a non-dimensional friction coefficient $\unicode[STIX]{x1D70F}a/(\unicode[STIX]{x1D707}\langle u\rangle _{s}\unicode[STIX]{x1D719}_{s})$ having only a marginal dependence on $d$ and $\unicode[STIX]{x1D719}_{s}$ . A linear fit to the data in figure 8(b) suggests a non-dimensional friction coefficient of approximately $7$ . For small values of $\unicode[STIX]{x1D719}_{s}$ the friction coefficient is smaller than 7 and non-constant (see inset), suggesting that the relation between $\unicode[STIX]{x1D70F}/\langle u\rangle _{s}$ and $\unicode[STIX]{x1D719}_{s}$ is linear only in an approximate sense.

In the limit $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ and $d\rightarrow \infty$ , we expect $F_{x}\simeq 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}aU$ (i.e. half the Stokes drag on a fully immersed sphere). Since $\langle u\rangle _{s}$ tends to $U$ as $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ (figure 5 a), and by definition $\unicode[STIX]{x1D70F}=F/L^{2}=F\unicode[STIX]{x1D719}_{s}/(\unicode[STIX]{x03C0}a^{2})$ for a square array, one might expect a friction coefficient of approximately $3$ in the limits $\unicode[STIX]{x1D719}_{s}\rightarrow 0$ and $d\rightarrow \infty$ . For small values of $\unicode[STIX]{x1D719}_{s}$ the data do appear to converge to a friction coefficient of 3 as $d$ increases (see inset of figure 8 b). However, the convergence is slow and for the values that we simulated the friction coefficient is significantly larger than 3 even for the smallest values of $\unicode[STIX]{x1D719}_{s}$ we considered. This fact may be due to the strong dependence of $\langle u\rangle _{s}$ on the surface fraction (figure 5 a). One would need to simulate truly negligible values of $\unicode[STIX]{x1D719}_{s}$ , and therefore use extremely large computational domains, to recover the expected asymptotic limit.

Figure 9.

Normalised pressure (a) and viscous (b) contributions to $\unicode[STIX]{x1D70F}$ .

The expansion in spherical harmonics in Physalis enables one to easily compute the pressure and viscous components of the hydrodynamic force on each sphere with excellent accuracy (Zhang & Prosperetti Reference Zhang and Prosperetti2005; Bluemink et al. Reference Bluemink, Lohse, Prosperetti and Van Wijngaarden2010; Botto & Prosperetti Reference Botto and Prosperetti2012). Upon normalisation, these force components yield $\unicode[STIX]{x1D70F}_{p}$ and $\unicode[STIX]{x1D70F}_{v}$ , namely the pressure and viscous contributions to $\unicode[STIX]{x1D70F}$ , respectively. Figure 9(a,b) shows $\unicode[STIX]{x1D70F}_{p}$ and $\unicode[STIX]{x1D70F}_{v}$ as a function of $\unicode[STIX]{x1D719}_{s}$ for different values of $d$ . As expected from the analytical solution for Stokes flow past a sphere, which suggests that pressure and viscous contributions to the drag are comparable, for $\unicode[STIX]{x1D719}_{s}\ll 1$ the viscous and pressure contributions to $\unicode[STIX]{x1D70F}$ are comparable in magnitude. For intermediate and relatively large values of $\unicode[STIX]{x1D719}_{s}$ , viscous stresses produce the dominant contribution to $\unicode[STIX]{x1D70F}$ . For instance, for $d=7a$ and $\unicode[STIX]{x1D719}_{s}>0.5$ , $\unicode[STIX]{x1D70F}_{v}$ is approximately one order of magnitude larger than $\unicode[STIX]{x1D70F}_{p}$ . The contribution $\unicode[STIX]{x1D70F}_{p}$ is only significant for small gap sizes, particularly at the highest values of the area fraction. This result confirms the intuitive notion that the largest contribution to $\unicode[STIX]{x1D70F}$ for a moderately dense monolayer originates from the shear forces exerted by the fluid on the portion of the sphere surfaces where the streamwise fluid velocity is larger.

For completeness, we also report in figure 10 the hydrodynamic force on each particle in the square monolayer as a function of $d$ for different values of $\unicode[STIX]{x1D719}_{s}$ . In contrast to the previous figures in which the tangential stress due to the bulk flow was reported, in figure 10 the force is not normalised by $L^{2}$ , and therefore the dependence on $\unicode[STIX]{x1D719}_{s}$ is perhaps more clear. This graph could be useful to quantify how the drag on a surface-active particle embedded at the boundary of a thin film between a particle-laden interface and a wall changes as a function of the thickness of the film.

Figure 10.

Hydrodynamic force on each sphere in the square array as a function of the gap size.

3.4 Reticulated array

The slip length depends on the specific microstructure of the monolayer. Two categories of microstructures are particularly relevant to particle-laden interfaces: well-dispersed systems, in which the particles are not in contact when $\unicode[STIX]{x1D719}_{s}$ is smaller than the maximum packing fraction; and percolating networks, where the particles form connected chains that span the boundaries of the interface (see e.g. Aveyard et al. Reference Aveyard, Clint, Nees and Paunov2000). The square array case examined in the previous section belongs to the first category. The percolating network case is examined in the current section. To get initial insights into more realistic situations, in which the percolating networks are disordered and are characterised by a variety of scales, in this section we simulate a periodic mesh-like reticulated arrangement. One example of such an arrangement is illustrated in figure 11(b).

Figure 11.

Schematic of a square array (a) and a reticulated array (b) for the same value of the solid area fraction, $\unicode[STIX]{x1D719}_{s}=43.63\,\%$ .

Table 1.

Normalised slip length, $\unicode[STIX]{x1D706}/a$ , and normalised tangential stress due to the bulk flow, $\unicode[STIX]{x1D70F}a/\unicode[STIX]{x1D707}U$ , for reticulated arrays characterised by different values of the area fraction $\unicode[STIX]{x1D719}_{s}$ and mesh size $L$ , for $d/a=7$ . The last two columns report the values of $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D70F}$ for a square array having the same solid area fraction as the corresponding reticulated array.

The simulations for the reticulated array case are carried out by including in the doubly periodic simulation box a variable number of spheres. The spheres form two orthogonal rectilinear chains that intersect each other in the centre of the domain. One of the chains is oriented along the flow direction.

Because of the need to consider several particles, simulating reticulated arrays requires a significantly larger domain than for a square array, for a given area fraction. To limit the number of simulations while allowing the exploration of the relevant parameter space, the majority of the simulations presented in this section are carried out for a fixed gap size, $d=7a$ . Preliminary tests confirm that this gap size is sufficiently large for the results to reasonably approximate the unbounded case $d=\infty$ .

Table 1 summarises values of the slip length $\unicode[STIX]{x1D706}$ and tangential stress $\unicode[STIX]{x1D70F}$ due to the bulk flow for different values of $L$ (and therefore of $\unicode[STIX]{x1D719}_{s}$ ). The last two columns in table 1 report the values of $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D70F}$ for a square array having the same area fraction as the corresponding reticulated array. These two quantities are denoted by the symbols $\unicode[STIX]{x1D706}_{sq}$ and $\unicode[STIX]{x1D70F}_{sq}$ , respectively.

Comparison of $\unicode[STIX]{x1D706}$ with $\unicode[STIX]{x1D706}_{sq}$ shows that, for a given surface coverage and particle size, the reticulated array gives a significantly larger value of the slip length. The effect of the specific microstructure (square array versus reticulated array) on $\unicode[STIX]{x1D706}$ becomes more and more significant as $\unicode[STIX]{x1D719}_{s}$ decreases. The difference between $\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D70F}_{sq}$ is comparatively small.

Figure 12.

(a) Normalised slip velocity and (b) normalised bulk shear rate at the interface as a function of solid area fraction for $d/a=7$ , comparing square and reticulated arrays.

Why does the reticulated array give a larger slip length? The slip length is defined as the ratio of $\langle u\rangle _{s}$ to $\langle \dot{\unicode[STIX]{x1D6FE}}\rangle _{s}$ . A difference in slip length can therefore be attributed to: (i) differences in $\langle u\rangle _{s}$ ; (ii) differences in $\langle \dot{\unicode[STIX]{x1D6FE}}\rangle _{s}$ ; or (iii) the combined effect of both quantities. Figure 12(a) compares values of $\langle u\rangle _{s}$ obtained with the reticulated array to those obtained with the square array. The corresponding values of $\langle \dot{\unicode[STIX]{x1D6FE}}\rangle _{s}$ are shown in figure 12(b).

It can be seen that the values of $\langle u\rangle _{s}$ given by the two microstructures are comparable only for relatively large values of $\unicode[STIX]{x1D719}_{s}$ . For moderate and small values of $\unicode[STIX]{x1D719}_{s}$ , the reticulated array gives a significantly larger value of $\langle u\rangle _{s}$ than the square array. Over a similar range of values of $\unicode[STIX]{x1D719}_{s}$ , the differences in the values of $\langle \dot{\unicode[STIX]{x1D6FE}}\rangle _{s}$ corresponding to the two microstructures are relatively small (figure 12 b). The observed dependence of the slip length on the microstructure is thus mainly due to changes in the slip velocity, while changes in the bulk shear rate at the interface play a relatively marginal role.

Figure 13.

Contours of normalised streamwise velocity, $u/U$ , in the plane $z=a$ for (a) square array and (b) reticulated array. The solid area fraction in (a) and (b) is $\unicode[STIX]{x1D719}_{s}=43.63\,\%$ and $d/a=7$ .

A question arises as to why the slip velocity is larger in the case of the reticulated array. Comparison of figure 11(a) with figure 11(b) shows that, for a given area fraction, the reticulated array is characterised by larger connected regions not occupied by particles. Because the flow in these regions is relatively unobstructed by the particles, the slip velocity in these regions could be substantially larger than the slip velocity in the interstices between the particles in the square array case. To verify this hypothesis, we plotted isocontours of $u$ in the plane $z=a$ , comparing the reticulated and square arrays (figure 13). The characteristic velocities in the open areas bounded by chains in the reticulated array case are at least $50\,\%$ larger than those in the interstices between particles in the square array case. These relatively large velocities are spread over regions of linear size comparable to the mesh size $L$ . In contrast, in the square array the smaller slip velocities are spread over regions of characteristic linear dimension $\ell \ll L$ , where $\ell$ is the average inter-particle separation. We have $\unicode[STIX]{x1D719}_{s}=(L/a-1)\unicode[STIX]{x03C0}a^{2}/L^{2}$ for the reticulated array and $\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x03C0}a^{2}/\ell ^{2}$ for the square array. Therefore, the ratio $L^{2}/\ell ^{2}$ of the interfacial areas occupied by high-velocity fluid (reticulated array) to that occupied by low-velocity fluid increases with decreasing $\unicode[STIX]{x1D719}_{s}$ . Since the slip velocity is an area-averaged quantity, larger local velocities spread over larger areas will give a larger value of the slip velocity.

We also note that, despite significant differences in slip length between the reticulated array and the square array, the interfacial shear stresses are comparable in the two cases (cf. table 1).

Figure 14.

Percentage contribution to the tangential stress due to the bulk flow  $\unicode[STIX]{x1D70F}$ for each particle in a reticulated array, $\unicode[STIX]{x1D719}_{s}=20.8\,\%$ and $d=7a$ . The horizontal chain in the diagram is parallel to the flow direction.

Comparing the values of $\unicode[STIX]{x1D70F}$ , we note that, despite significant differences in slip length between the reticulated array and the square array, the values of the tangential stress due to the bulk flow are comparable in the two cases (cf. table 1). A significant difference between the square and reticulated arrays is that, while in the square array each particle will contribute equally to $\unicode[STIX]{x1D70F}$ , in the reticulated array different particles can be subject to different hydrodynamic forces, and therefore the local tangential stress will in general be non-uniform. To characterise the degree of non-uniformity in $\unicode[STIX]{x1D70F}$ , we show in figure 14 the fraction of the total value of $\unicode[STIX]{x1D70F}$ contributed by each particle. This fraction is calculated as the ratio of the hydrodynamic force on each particle to the sum of the hydrodynamic forces exerted by all the particles within the computational domain.

The chains of particles arranged perpendicular to the flow direction, or transverse chains, carry the largest contribution to $\unicode[STIX]{x1D70F}$ . In the case considered in figure 14, the transverse chain contributes more than $66\,\%$ of the total value of $\unicode[STIX]{x1D70F}$ . The hydrodynamic force acting on the transverse chain is strongly non-uniform: the largest contribution, $12\,\%$ , is associated with the particles located near the midpoint of the transverse chain (farthest away from the intersection point). The central particle, located at the intersection between the longitudinal and transverse chains, contributes only $5.3\,\%$ of the total value of $\unicode[STIX]{x1D70F}$ . In comparison, particles belonging to the longitudinal chain, oriented along the flow direction, are subject to a relatively uniform bulk tangential stress.

Figure 15.

Normalised slip length versus (a) normalised mesh size and (b) normalised gap size; reticulated array, $d=7a$ . The area fraction in panel (b) is $43.63\,\%$ . The dashed line in panel (a) is (3.4) for $A=-0.039$ and $B=0.03$ .

The natural choice of length scale to characterise the flow past the reticulated array is the mesh size $L$ . It is therefore expected that the slip length will scale as $\unicode[STIX]{x1D706}=Lf(\unicode[STIX]{x1D719}_{s})$ , where $f(\unicode[STIX]{x1D719}_{s})$ is a function having relatively weak dependence on $\unicode[STIX]{x1D719}_{s}$ (and therefore on $L$ ). In figure 15(a), the values of $\unicode[STIX]{x1D706}$ for the reticulated array are plotted as a function of $L$ . Over the range of values simulated, the slip length is seen to increase only slightly faster than linearly with $L$ , confirming a scaling of the type $\unicode[STIX]{x1D706}\sim Lf(\unicode[STIX]{x1D719}_{s})$ . Davis & Lauga (Reference Davis and Lauga2009) carried out an analysis of the slip length for a flat partial slip surface composed of a mesh-like distribution of thin solid regions, obtaining a semi-empirical relation of the form

where $A$ and $B$ are constants. The values obtained by Davis & Lauga (Reference Davis and Lauga2009) were $A=-0.107$ and $B=0.003$ . Fitting the simulation data (dashed line in figure 15 a), we obtain $A=-0.039$ and $B=0.03$ . As for the square array, the effect of the protuberance of the no-slip region in the liquid gives a larger value of $\unicode[STIX]{x1D706}$ for a given value of $\unicode[STIX]{x1D719}_{s}$ as compared to a flat partial slip surface. We have attempted a simple rescaling by a factor of $9\unicode[STIX]{x03C0}/16$ as done for a square array (see § 3.2) to account for the difference between the results of Davis & Lauga (Reference Davis and Lauga2009) and ours, but the results have not been as satisfying as for the square array case. This is due to the fact that, for a percolating network, hydrodynamic interactions between the particles are important even in the dilute limit. Therefore, the ratio of the hydrodynamic force on a single sphere to that of a single disk cannot represent a good rescaling factor.

The results above suggest that – in analogy with the case of superhydrophobic surfaces – different scalings for the slip length hold depending on whether the solid or the liquid are the continuous phase in the plane of the particle-laden interface. In the reticulated network, the solid is the continuous phase, and the free-slip interface is the dispersed phase. On the contrary, in the case of the square array, the solid is the dispersed phase. When the solid is the dispersed phase and the fluid is the continuous phase in the particle-laden interface, the particle radius is the characteristic length scale, and the slip length scales as $\unicode[STIX]{x1D706}\propto a/\unicode[STIX]{x1D719}_{s}$ (for sufficiently small $\unicode[STIX]{x1D719}_{s}$ ). When the solid is the continuous phase, the mesh size $L$ is the characteristic length, and the influence of the particle radius enters only into a weak – potentially logarithmic – dependence on the solid fraction.