Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://nanofluidics.phys.msu.ru/paper/JCP2013.pdf Äàòà èçìåíåíèÿ: Fri Nov 8 16:49:17 2013 Äàòà èíäåêñèðîâàíèÿ: Thu Feb 27 20:12:07 2014 Êîäèðîâêà:

Effective slippage on superhydrophobic trapezoidal grooves
Jiajia Zhou, Evgeny S. Asmolov, Friederike Schmid, and Olga I. Vinogradova Citation: J. Chem. Phys. 139, 174708 (2013); doi: 10.1063/1.4827867 View online: http://dx.doi.org/10.1063/1.4827867 View Table of Contents: http://jcp.aip.org/resource/1/JCPSA6/v139/i17 Published by the AIP Publishing LLC.

Additional information on J. Chem. Phys.
Journal Homepage: http://jcp.aip.org/ Journal Information: http://jcp.aip.org/about/about_the_journal Top downloads: http://jcp.aip.org/features/most_downloaded Information for Authors: http://jcp.aip.org/authors


THE JOURNAL OF CHEMICAL PHYSICS 139, 174708 (2013)

Effective slippage on superhydrophobic trapezoidal grooves
Jiajia Zhou,1 Evgeny S. Asmolov,
1

2,3,4

Friederike Schmid,1 and Olga I. Vinogradova2,5,6

Institut fÝr Physik, Johannes Gutenberg-UniversitÄt Mainz, D55099 Mainz, Germany 2 A.N. Frumkin Institute of Physical Chemistry and Electrochemistry, Russian Academy of Science, 31 Leninsky Prospect, 119071 Moscow, Russia 3 Central Aero-Hydrodynamic Institute, 140180 Zhukovsky, Moscow Region, Russia 4 Institute of Mechanics, M.V. Lomonosov Moscow State University, 119991 Moscow, Russia 5 Department of Physics, M.V. Lomonosov Moscow State University, 119991 Moscow, Russia 6 DWI, RWTH Aachen, Forckenbeckstraúe 50, 52056 Aachen, Germany

(Received 19 August 2013; accepted 17 October 2013; published online 6 November 2013) We study the effective slippage on superhydrophobic grooves with trapezoidal cross-sections of various geometries (including the limiting cases of triangles and rectangular stripes), by using two complementary approaches. First, dissipative particle dynamics (DPD) simulations of a flow past such surfaces have been performed to validate an expression [E. S. Asmolov and O. I. Vinogradova, J. Fluid Mech. 706, 108 (2012)] that relates the eigenvalues of the effective slip-length tensor for one-dimensional textures. Second, we propose theoretical estimates for the effective slip length and calculate it numerically by solving the Stokes equation based on a collocation method. The comparison between the two approaches shows that they are in excellent agreement. Our results demonstrate that the effective slippage depends strongly on the area-averaged slip, the amplitude of the roughness, and on the fraction of solid in contact with the liquid. To interpret these results, we analyze flow singularities near slipping heterogeneities, and demonstrate that they inhibit the effective slip and enhance the anisotropy of the flow. Finally, we propose some guidelines to design optimal one-dimensional superhydrophobic surfaces, motivated by potential applications in microfluidics. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4827867]
I. INTRODUCTION

diagonalized by a rotation S = cos - sin sin cos . (2)

The design and fabrication of micro- and nanotextured surfaces has received much attention during the past decade. In case of a hydrophobic texture, a modified surface profile can lead to a very large contact angle, which induces exceptional wetting properties.1 A remarkable mobility of liquids on such superhydrophobic surfaces is observed in the Cassie state, where the textures are filled with gas, which renders them "self-cleaning" and causes droplets to roll (rather than slide) under gravity2 and rebound (rather than spread) upon impact.3, 4 Furthermore, patterned superhydrophobic materials are important in context of fluid dynamics and their superlubricating properties.5, 6 In particular, superhydrophobic heterogeneous surfaces in the Cassie state exhibit very low friction, and this drag reduction is associated with the large slippage of liquids.7­9 It is very difficult to quantify the flow past heterogeneous surfaces. However, analytical results can often be obtained using an effective slip boundary condition, beff , at the imaginary smooth homogeneous, but generally anisotropic surface.9, 10 For anisotropic textures, the effective slip generally depends on the direction of the flow and is a tensor, beff ef {bij f } represented by a symmetric, positive definite 2 â 2 matrix,11

Equation (1 ) allows us to calculate an effective slip in any direction given by an angle , provided the two eigenvalues of the slip-length tensor, beff ( = 0) and beff ( = /2), are known. The concept of an effective slip length tensor is general and can be applied for arbitrary channel thickness.12 It is a global characteristic of a channel,9 and the eigenvalues usually depend not only on the parameters of the heterogeneous surfaces, but also on the channel thickness. However, for a thick channel (compared to a texture period, L) they become a characteristic of a heterogeneous interface solely. An important type of such surfaces are highly anisotropic textures, where a (scalar) partial local slip b(y) varies in only one direction. Such surfaces can be created experimentally by making textured surfaces with one dimensional surface profiles, and in the Cassie state, the local slip length will always reflect the relief of the texture. This is the consequence of the "gas cushion model," which takes into account that the dissipation at the gas/liquid interface is dominated by the shearing of a continuous gas layer,13, 14 b (y ) e (y ), g (3)

beff = S

beff 0

0
beff

S- ,

(1 )

where g and are dynamic viscosities of a gas and a liquid, and e is the thickness of the gas layer. Equation (3) represents an upper bound for a local slip at the gas area, which is
© 2013 AIP Publishing LLC

0021-9606/2013/139(17)/174708/11/$30.00

139, 174708-1


174708-2

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

attained in the limit of a small fraction of the solid phase. At a relatively large solid fraction15 or if the flux in the "pockets" of superhydrophobic surfaces is equal to zero,16, 17 Eq. (3) could overestimate the local slip. However, even in all situations the local slip length profile, b(y), necessarily follows the relief of the texture provided it is shallow enough, e(y) L.17 Since in case of water Eq. (3) leads to b(y) 50e(y), it is expected that the "gas cushion model" could be safely applied up to b(y)/L = O(10) or so. The flow properties on highly anisotropic surfaces can be highly peculiar. In particular, off-diagonal effects may be obtained where a driving force (such as pressure gradient, shear rate, or electric field) in one direction induces a flow in a different direction, with a measurable and useful perpendicular component. This could be exploited for implementing transverse pumps, mixers, flow detectors, and more.18, 19 The main task is to find the connection between the eigenvalues of the slip-length tensor and the parameters of such one-dimensional surface textures. Up to now analytical results for the effective slip length have only been obtained for the simplest geometry of rectangular grooves.12, 20­24 We are unaware of any prior work that has obtained analytical results for other types of one-dimensional textures. However, some recently derived relations for the case of an anisotropic isolated surface (corresponding to the thick channel limit) allow one to simplify the analysis of more complex geometries, and to rationalize the limiting behavior of the effective slip length in certain regimes. We mention first that the transverse component of the slip-length tensor was predicted to be exactly half the longitudinal component that one would obtain if the local slip were multiplied with a factor two, 2b(y),25 beff 2b (y ) /L . (4) 2 This relation can greatly simplify the analysis since it indicates that the flow along any direction of the one-dimensional surface can be easily determined, once the longitudinal component of the effective slip tensor is found. Equation (4) has been recently verified for a texture with sinusoidal superhydrophobic grooves by using a Lattice Boltzmann (LB) method,26 and for weakly slipping stripes by using dissipative particle dynamics (DPD) simulations.24 Another useful relation addresses the particular case where the local slip is large, b(y)/L 1. In that case, the effective slip length tensor beff becomes isotropic and is given by26, 27
beff b (y ) /L =

beff

beff

1 b (y )

-1

,

(5)

where 1/b(y) is the mean inverse local slip. Finally, for weakly slipping patterns with b(y)/L 1, the flow again becomes isotropic, and its value is equal to the surface average of the local slip length,10, 21 beff
beff

to Eq. (6) become questionable. Higher order contributions to the effective slip become comparable to the first-order term given by Eq. (6) and can no longer be ignored.24 Computer simulations and numerical modeling can shed light on flow phenomena past one-dimensional surfaces. In efforts to better understand the connection between the parameters of the texture and the effective slip lengths, several groups have performed continuum and Molecular Dynamics simulations of flow past anisotropic surfaces. Most of these studies focused on calculating eigenvalues of the effective slip-length tensor for superhydrophobic stripes.28­30 To study various aspects of a flow past striped walls, recent works employed Dissipative particle dynamics (DPD)23, 24, 31 and Lattice Boltzmann (LB)12 methods. In the case of stripes with a piecewise constant local slip profile (i.e., the profile features steplike jumps at the boundaries of the stripes), the singularity in the slip profile induces singularities both in the pressure and velocity gradient, which generates an additional mechanism for dissipation even at small local slip.24, 25 Anisotropic onedimensional texture with continuously varying patterns can produce a very large effective slip. This has been confirmed in recent LB simulations studies of sinusoidal textures.26 Onedimensional continuous surfaces may also exhibit singularities (for example, the first derivative may be discontinuous). However, we are unaware of any previous work that has addressed the question of effective slip lengths past such textures. In this paper, we study the friction properties on trapezoidal superhydrophobic grooves [see Fig. 1(a)]. This texture is more general than rectangular grooves since it provides an additional geometrical structure parameter, i.e., the trapezoid base width. Also, the trapezoid texture can be extended to its extremes, such as a triangular patterns, where the liquid/solid contact area vanishes, and rectangular reliefs. Superhydrophobic trapezoidal textures were shown to give large contact angles, and to provide stable Cassie states.32 An obvious advantage of the trapezoidal surface relief is that the manufactured texture becomes mechanically more stable against bending compared to rectangular grooves. Furthermore, it can be naturally formed from dilute rectangular grooves as a result of elasto-capillarity.33 This type of surfaces has already been intensively used in slippage experiments.34 However, the impact of trapezoidal texture reliefs on transport and hydrodynamic slippage remains largely unknown. Our study is based on the use of DPD simulations, which we have put forward in our previous papers,23, 24 on the numerical solution of the Stokes equations (a collocation method), and on theoretical analysis. In particular, we validate Eq. (4) for trapezoidal textures with various parameters, we apply Eq. (5) to evaluate an effective slip for strongly slipping patterns, we discuss the singularities in velocities caused by inflection points in the local slip profiles, which in particular lead to deviations from Eq. (6), and we finally suggest the optimal design of trapezoidal patterns.
II. MODEL AND GOVERNING EQUATIONS

b (y ) = b .

(6) We consider creeping flow along a plane anisotropic superhydrophobic wall, using a Cartesian coordinate

Note that Eq. (6) is expected to be accurate for a continuous local slip profile only. If the local slip length exhibits step-like jumps at the edges of heterogeneities, the expansions leading


174708-3

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

(a)

z x y e

(b) b
b2 b1

L 1L 2L

shear flow.7, 9 However, it does not apply to a thin37 or an arbitrary12 channel situation, where the effective slip scales with the channel width. Note that the flow in thin channels with superhydrophobic trapezoidal textures has been considered in recent work.17 Dimensionless variables are defined by using L as a reference length scale, the inverse shear rate 1/G sufficiently far from the surface as a reference time. At small Reynolds number, Re = GL2 / , the dimensionless velocity and pressure fields, v and p, satisfy Stokes equations · v = 0, (10) p - v = 0. Due to Eq. (4), it is sufficient to consider the flow parallel to the pattern (x-axis). In this case, the velocity has only one component in the x-direction v = (v (y, z), 0, 0), and we seek a solution for the velocity profile of the form v = U + u, (11)

y

FIG. 1. Sketch of the trapezoidal superhydrophobic surface in the Cassie state (a), and, according to Eq. (3), of the corresponding local slip length (b).

system (x, y, z) (Fig. 1). As in most previous publications,7, 12, 15, 20, 29, 35­37 we model the superhydrophobic plate as a flat heterogeneous interface. In such a description, we neglect an additional mechanism for a dissipation connected with the meniscus curvature.15, 38, 39 Note however that such an ideal situation has been achieved in many recent experiments.40­42 The origin of coordinates is placed at the flat interface, characterized by a slip length b(y), spatially varying in one direction, and the texture varies over a period L. The local slip length is taken to be a piecewise linear function. It includes two regions with a constant slip: one where the slip length assumes its minimum value, b1 , which has area fraction of 1 , and another with area fraction 2 where the slip length is maximal. These regions are connected by two transition regions with equal surface fractions (1 - 1 - 2 )/2, where the local slip length varies linearly from b1 to b2 . If we rewrite the coordinate y in the form y = (n + t)L, where n is an integer number and |t | 1 , the local slip-length 2 profile can be expressed in terms of t: |t | 1 b1 , 2 1 1 2(b2 -b1 ) b(y ) = b1 + 1-1 -2 |t |- 2 , 2 < |t | 1-2 . (7) 2 1-2 b2 , < |t | 1 2 2 For 1 = 2 = 0, we get a triangular profile b(y ) = b1 + 2(b2 - b1 )|t |, |t | 1 . 2 (8)

where U = z is the undisturbed linear shear flow. The perturbation of the flow, u(y, z), which is caused by the presence of the texture and decays far from the interface. It follows from Stokes' equations that y p = z p = 0. Since the pressure is constant far from the superhydrophobic surface this implies that it also remains constant in the entire flow region. The velocity perturbation satisfies a dimensionless Laplace equation, u = 0. (12)

The boundary conditions at the wall and at infinity are defined in the usual way, z=0: u- b (y ) b (y ) z u = , L L (13)

z : z u = 0.

(14)

The general solution to Eq. (12), decaying at infinity, can be represented in terms of a cosine Fourier series as u= a0 + 2


a n exp(-2nz) cos(2ny ),
n=1

(15)

Note that regular trapezoidal and triangular profiles are continuous functions (albeit with a discontinuous first derivative). However, in the limit 1 + 2 = 1, the profile reduces to the discontinuous rectangular slip length profile (alternating stripe profile) studied in earlier work, b (y ) = b1 , b2 , |t |
1 2 1 2 1 2

where an are coefficients to be found by applying the appropriate boundary condition (13). It can be written in terms of the Fourier coefficients as a0 + 2


1 + 2n
n=1

b (y ) b (y ) (y ) a n cos (2ny ) = . (16) L L

< |t |

.

(9)

Our results apply to a single nel, where effective hydrodynamic at the scale of roughness, so that ciently far above the surface may

surface in a thick chanslip is determined by flow the velocity profile suffibe considered as a linear

Equation (16) was used to determine the longitudinal effective slip numerically. This was done using a method similar to Ref. 24. We truncate the sum in Eq. (16) at some cutoff number N (usually N = 1001) and evaluate it in the points yl = l/2(N - 1), where l are numbers varying from 0 to N - 1. The problem is then reduced to a linear system of equations, Aln a n L = bl , where Al0 = 1/2; Aln = [1 + 2nb(yl )/L] cos(2nyl ), n > 0 and bl = b(yl ). The


174708-4

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

system is solved using the IMSL routine LSARG. The transverse effective slip was then calculated using Eq. (4).
III. DISSIPATIVE PARTICLE DYNAMICS SIMULATION

with the velocity replaced by the relative particle velocity with respect to the wall. For immobile walls, the dissipative part can be written as FD = -L L (z)vi . i (22)

Mesoscale fluid simulations were performed using the dissipative particle dynamics (DPD) method.43­45 More specifically, we use a version of DPD without conservative interactions46 and combine that with a tunable-slip method47 to model patterned surfaces. The tunable-slip method allows one to implement arbitrary hydrodynamic boundary condition in coarse-grained simulations. In the following, we give a brief description of the model and introduce the implementation of planar surfaces with arbitrary one-dimensional slip profile. The basic DPD equations involve pair interactions between fluid particles. The pair forces are equal in magnitude but opposite in direction; thus the momentum is conserved for the whole system. This leads to the correct long-time hydrodynamic behavior (i.e., Navier-Stokes). The pair force consists of a dissipative contribution and a random component, FDPD = FD + FR . ij ij ij (17)

The parameter L characterizes the strength of the wall friction and can be used to tune the value of the slip length. For example, L = 0 results in a perfectly slippery wall, whereas a positive value of L leads to a finite slip length. The random force has the form FR = i 2kB TL L (z) i , (23) where each component of i is a random variable. Again, the magnitude of the random force is chosen such that the fluctuation-dissipation theorem is satisfied. For homogeneous surfaces, a useful analytical expression can be derived which establishes a relation between the slip length b and the simulation parameters (the wall friction parameter L and the cutoff rc ) to a good approximation.47 However, the accuracy of the analytic result is not satisfactory at small slip lengths. For the purpose of the present work, we need an accurate expression for the relation between the slip length and the tuning parameter L . Therefore, we have simulated plane Poiseuille and Couette flows with different L to obtain the slip length and the position of the hydrodynamic boundary. The simulation results were fitted to a truncated Laurent series up to second order in L and first order in 1/ L . The simulation data and the fitted curve are shown in Fig. 2. The fit yields the expression b= 1.498 2 - 0.3291 + 0.01457L - 0.001101L . L (24)

The dissipative part FD is proportional to the relative velocity ij between two particles vij = vi - vj , FD = - ij
DPD

^^ D (rij )(vij · rij )rij ,

(18)

with a local distance-dependent friction coefficient DPD D (rij ). The weight function D (rij ) is a monotonically decreasing function of rij and vanishes at a given cutoff, and DPD is a parameter that controls the strength of the dissipation. The random force FR has the form ij FR = ij 2kB T
DPD

^ D (rij )ij rij ,

(19)

where ij = ji is a symmetric random number with a zero mean and a unit variance. The amplitude of the stochastic contribution is related to the dissipative contribution by a fluctuation-dissipation theorem in order to ensure that the model reproduces the correct equilibrium distribution. The tunable-slip method47 introduces the wall interaction in a similar fashion. The force exerted by the wall on particle i is given by Fwall = FWCA + FD + FR . i i i i (20)

This formula relates the wall friction coefficient L , which is set in the simulation, to the slip length b. Also shown in Fig. 2 is the analytic result in dashed line, which agrees with the simulation data quite well at large slip length. One important observation is that both the slip length b and the friction coefficient L are local quantities. Thus, it is reasonable to assume that the relation (24) also holds for
350 300 250 200 b [] 150 100 50 5 4 3 2 1 0 -1 0 1 2 3 4 5 6 7 8 9 10

The first term is a repulsive interaction that prevents the fluid particles from penetrating the wall. It can be written in term of the gradient of a Weeks-Chandler-Andersen potential,48 FWCA = - V (z), i V (z ) = 4 0
12 z

-

6 z

+

1 4

z < 6 2 , z 6 2

0

(21)

-50 0.01 0.1 L [(m)
1/2

1 /]

10

where z is the distance between the fluid particle and the wall, assuming the wall lies in xy plane, and and set the energy and length scales. In the following, physical quantities will be reported in a model unit system of (length), m (mass), and (energy). The dissipative contribution is similar to Eq. (18),

FIG. 2. The relation between the slip length b and the wall friction coefficient L , for a wall interaction cutoff 2.0 . The inset shows an enlarged portion of the region where the slip length is close zero. The no-slip boundary to condition is implemented by using L = 5.26 m / . The solid lines are the fitting result Eq. (24), and the dashed curves show the analytical prediction.


174708-5

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

TA B L E I . Parameters used in the DPD simulations. Fluid density Friction coefficient for DPD interaction DPD Cutoff for DPD interaction No-slip boundary condition L (b = 0) Cutoff for wall interaction Dynamic viscosity Hydrodynamic boundary from the simulation wall Simulation box Pattern period, L 3.75 -3 5.0 m / 1.0 5.26 m / 2.0 1.35 ± 0.01 m / 1.06 ± 0.12 20 â 40 â 82 40

the pressure gradient (the body force in simulation). To obtain accurate simulation data at large b2 /L, one has to reduce the body force significantly, and the necessary simulation times increase prohibitively. Therefore, only simulation results for b2 /L = O(1) and smaller are presented for this texture.
2

IV. EIGENVALUES OF THE SLIP-LENGTH TENSOR

inhomogeneous surfaces, as long as the slip length varies on length scales much larger than the range of the weight functions, L (z) and D (r). For a patterned surface with L much larger than , the local friction coefficient L (y) can be calculated from the slip profile b(y) by inverting Eq. (24). The resulting profile L (y) is then used in simulations to set up the patterned surfaces. The simulations are carried out using the open source simulation package ESPResSo.49 Modifications have been made to incorporate the patterned surfaces. Table I summarizes the simulation parameters. In simulations, the patterned surface has a period of L = 40 . The simulation box is a rectangular cuboid of size 20 â 40 â 82 , periodic boundary conditions are assumed in the xy plane, and the two symmetric patterned surfaces are located at z = 0 and z = 82 .The separation between the two surfaces is about twice the pattern periodicity. We have checked that increasing the channel thickness does not change the results. The simulation starts with randomly distributed particles inside the channel, and flow is induced by applying a body force to all particles. Body forces used in this simulation were chosen small enough that the flow velocity near the wall was less than 0.1 / m, in order to ensure that the flow is strictly laminar and effects of finite Reynolds numbers do not affect the simulation results. The magnitude of the slip velocity also affects the spatial resolution in modeling the variation of the slip length. Since the trajectories of the fluid particles are integrated at discrete time steps, the friction experienced by the particle between time steps is determined from the early instant. Therefore, we used a small time step of 0.01 m/ . Due to the large system size, the relaxation time required for reaching a steady state was very large as well (over 106 time steps). The flow velocity is small in comparison with the thermal fluctuation; thus many time steps are necessary to gather sufficient statistics. In this work, velocity profiles were averaged over 4 â 105 time steps. A fit to the plane Poiseuille flow then gives the effective slip length beff .23 The error bars are obtained by six independent runs with different initialization. Finally, we note that in case of a triangular local slip, only a very small region near the no-slip point contributes to the shear stress at large b2 . This presents a complication for simulations since fluid particles with a finite translational velocity travel a small distance over the no-slip region before reaching the zero averaged velocity. This distance is related to the initial velocity of the fluid particles when they enter the no-slip region, which is in general proportional to b2 and

Here we present the results for the eigenvalues of the sliplength tensor obtained by simulations of different local slip profiles. We first verify Eq. (4), derived earlier by some of us25 to relate the longitudinal and transverse effective slip lengths. We use local slip profiles with 1 = 2 . In that special case, the profile b(y) can be characterized by three parameters: the area fraction 1 of the no-slip region, the surface-averaged slip length b = (b1 + b2 )/2, and the heterogeneity factor, , which is a measure for the relative variation of the slip length and defined as = b2 - b1 b2 - b1 = . 2b b2 + b1 (25)

A homogeneous surface then corresponds to = 0, and low values of corresponds to small b2 /b1 . Patterns involving noslip regions (i.e., b1 = 0) have = 1. Such a no-slip assump1, which is typical for tion can be justified provided b2 /b1 many situations. Indeed, for most smooth hydrophobic surfaces, the slip length is of the order of a few tens of nm,50­53 which is very small compared to the slip length values of up to tens or even hundreds of m slip lengths that may be obtained on the gas areas.

A. Fixed average slip

Let us first investigate the effect of on the effective slip at fixed b /L = 0.5. Figure 3(a) shows the simulation results for eigenvalues of the slip length tensor as a function of obtained for a trapezoidal local slip ( 1 = 0.25). Also included are numerical results (solid curves) calculated with Eqs. (4) and (16). The agreement between the simulation and numerical results is quite good for all . This demonstrates both the accuracy of our simulations, and the validity of Eq. (4). The data presented in Fig. 3(a) show that the longitudinal slip is larger than the transverse slip, i.e., the flow is anisotropic, and the anisotropy increases with . The effective slip assumes its maximum value at = 0, when the surface is homogeneous, and then decreases monotonously. This implies that the effective slip is found to be always smaller than the surfaceaveraged value, and the difference from the average value increases with the heterogeneity. A similar behavior was reported in prior work for other types of patterns.26, 54 In the case of large local slip, the effective slip length can be obtained from Eq. (5), which gives for a trapezoidal slip profile, Eq. (7): beff 1 b2 2 1 - 1 - 2 + + ln b1 b2 2 (b2 - b1 ) b1
-1

.

(26)


174708-6

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

(a)

0.60 0.50 0.40

0.30 0.20 0.10 0.00 0 0.2 0.4 0.6 0.8 1

To examine the difference between the different slip profiles, the results for a longitudinal effective slip for a texture with a trapezoidal local slip from Fig. 3(a) are reproduced in Fig. 3(b). Also included in Fig. 3(b) are numerical and simulation data for rectangular ( 1 = 0.5) and triangular ( 1 = 0) profiles with the same average slip. The data show that the triangular profile gives a larger, and the rectangular one a smaller effective slip than that resulting from the trapezoidal local slip. We note that 1 is very different for these slip-length profiles. This suggests that one of the key parameters determining the effective slip is the area fraction of solid in contact with the liquid. If it is very small, the effective slip length becomes large. Another important factor is the presence of flow singularities near slipping heterogeneities as we will discuss below.
B. Fixed heterogeneity factor

(b)
0.50

beff/L

0.40 beff/L

0.30

We now fix = 1 and consider flows that are maximally anisotropic for a given pattern, which corresponds to b1 = 0 (no-slip areas) and b2 = 2 b . It is instructive to focus first on the role of b2 . Figure 4 shows the eigenvalues of the slip-length tensor as a function of b2 /L as obtained from

0.20

0.10 0 0.2 0.4
FIG. 3. (a) Longitudinal (open symbols) and transverse (closed symbols) effective slip for a trapezoidal texture with an average slip b /L = 0.5. Solid curves show corresponding numerical solutions. Dash-dotted curve shows asymptotic solution in the limit of the large local slip [Eq. (26)]. (b) Longitudinal effective slip length simulated (from top to bottom) for triangular, trapezoidal, and rectangular profiles with the same average slip b /L = 0.5 (symbols) and corresponding numerical data.

(a)
0.30 0.6 0.8 1

beff/L

0.20

0.10

0.00 0.1

Figure 3(a) includes the theoretical curve expected for this (isotropic) case (dash-dotted line). We note that the slip length is always smaller than the average slip, b /L = 0.5. Perhaps the most interesting and important result here is that Eq. (26) remains surprisingly accurate even in the case of moderate, but not vanishing, local slip, provided is small enough. We remind the reader that Eq. (26) may not be used when = 1, i.e., for textures with no-slip regions, b1 = 0, where it would simply predicts the effective slip length to be zero. Qualitatively similar results were obtained for rectangular stripes and triangular local slip profiles. Therefore we do not show these results here. We only recall that in the case of stripes, the asymptotic result for large slip gives (see Refs. 26 and 28 for the original derivations), beff 1 2 + b1 b2
-1

1 b2/L

10

(b)
1.60

1.20 beff/L

0.80

0.40

0.00 0.1 1 b2/L 10

,

(27)

and for triangular profiles, one can derive beff 2(b2 - b1 ) . ln(b2 /b1 ) (28)

FIG. 4. Effective slip lengths as a function of b2 /L at = 1 for (a) trapezoidal and (b) triangular local slip. Symbols show simulation data for longitudinal (open) and transverse (closed) effective slip lengths. The solid curves are results of numerical solutions. Dash-dotted curves show the predictions of analytic formulae [Eqs. (29) and (30)] for stripes with the same fraction of no-slip areas. Dashed lines show the asymptotic values for stripes at large b2 /L [Eq. (31)].


174708-7

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

simulations (symbols). As before we include numerical results (solid curves) calculated with Eqs. (4) and (16). We remind that the "gas cushion model", Eq. (3), becomes very ap1, and the local slip profile could not proximate at b2 /L necessarily follow the relief of the texture. It is, however, instructive to include large values of b2 /L into our consideration below since this will provide us with some guidance. The eigenvalues of the slip-length tensor for a superhydrophobic trapezoidal texture are shown in Fig. 4(a). The agreement between the theoretical and simulation data is very good for all b2 /L, again confirming the accuracy of the simulations and the validity of Eq. (4). Qualitatively, the results are similar to those obtained earlier for a striped surfaces:21 At very small b2 /L, the effective slip is small and appears to be isotropic, in accordance to predictions of Eq. (6). It becomes anisotropic and increases for larger b2 /L, and saturates at b2 /L 1. In case of stripes with alternating regions of no-slip (b1 = 0) and partial-slip (b2 = 2 b ) areas the situation can be described by simple analytic formulae,21 beff L 2 1 2 L 2 + tan 1+ ln sec b2 2 2 ln sec

0.12 0.10 0.08 beff/L 0.06 0.04 0.02 0.00 0 0.1 0.2 1
FIG. 5. Effective slip lengths as functions of 1 for surfaces with = 1 and b /L = 0.125. The open and closed symbols are simulation results for the longitudinal and transverse effective slip lengths, respectively. Solid lines are numerical results. Dashed line shows predictions of Eq. (6).

0.3

0.4

0.5

,
2

(29)

beff L

1 2

2 2 L 2 + tan 1+ ln sec 2b2 2 2 ln sec

.
2

(30)

L, these expressions turn into those deIn the limit of b2 rived earlier for perfect slip stripes20, 55 beff 2b = eff L L 1 ln sec 2
2

.

(31)

For comparison, we plot the above expressions (using the value 1 = 0.25 of the trapezoidal texture) in Fig. 4(a) as dash-dotted [Eqs. (29) and (30)] and dashed [Eq. (31)] lines. One can see that with these parameters, the results for stripes practically coincides with those for trapezoidal slip at very small and very large b2 /L. At intermediate local slip, there is a small discrepancy suggesting that the equations for partial slip stripes overestimate the effective slip for trapezoidal texture. We will return to a discussion of these results later in the text. Figure 4(b) plots the simulation results for a triangular local slip profile, which is particularly interesting since 1 = 0. As explained above, simulation data could only be obtained for small values of b2 in this case. In contrast, our numerical approach allowed us to calculate the effective slip length over a much larger interval of b2 /L. These numerically computed slip lengths are also included in Fig. 4(b). In the regime where simulation data are available, the agreement with the simulation and numerical results is excellent, therefore validating Eq. (4). It is interesting to note that according to the numerical calculations, the effective slip lengths grow weakly with b2 /L at large b2 /L and likely do not saturate, in contrast to rectangular and trapezoidal textures.

Finally, we compare the effective slip of surfaces with the same and b /L, but different 1 . Figure 5 shows the corresponding curves for a maximal heterogeneity factor, ( = 1), and a relatively small b /L = 0.125. The results for larger b /L are qualitatively similar and not shown here. Only the values of the effective slip lengths become larger.) We note that by varying 1 from 0 to 0.5, we change the local slip profile from triangular ( 1 = 0) to rectangular (stripes, 1 = 0.5). Regular trapezoidal textures correspond to intermediate values of 1 . Figure 5 also shows the analytical prediction of Eq. (6) for weakly slipping surfaces, i.e., an effective slip , b which is isotropic and independent on 1 . The efbeff fective slip according to simulations and numerical results is always smaller than the first-order theoretical prediction, and the deviations from the area-average slip model increase with 1 . Furthermore, the data show that the longitudinal and transverse slip lengths are different from each other, i.e., the flow is anisotropic, even though it should be nearly isotropic according to the asymptotic theory, Eq. (6). For striped patterns, similar deviations from the asymptotic theory were shown to result from flow perturbations in the vicinity of jumps in the slip length profile.24 Here, we find that the friction and anisotropy of the flow past triangular and trapezoidal, weakly slipping surfaces is increased, compared to the expectation from the asymptotic theory, Eq. (6). This somewhat counter-intuitive result suggests that similar, although perhaps weaker, viscous dissipation and singularities of the shear stress can also appear if singularities in the slip length profile only appear in its first derivative. Below we will focus on this phenomenon.
V. FLOW SINGULARITIES NEAR SLIPPING HETEROGENEITIES

Motivated by the findings from Sec. IV, we now discuss the flow near slipping heterogeneity in more detail. We remind the reader that the shear stress is known to be singular near sharp corners of rectangular hydrophilic grooves, being proportional to r-1/3 for longitudinal configurations, and to r-0.455 for transverse configurations56 (r is the distance from the corner). In the case of a striped superhydrophobic surface,


174708-8

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

(a)

r

(b)

1

r



0.9 0.8 0.7 0.5 0.4 0.3 0.2 0.1 1 2 3 4 5 6 7 8 9 10
FIG. 7. Exponents = 1/2 (see text for definition) for the trapezoid [Eq. (34), solid line] and triangle slip profiles [Eq. (38), dashed line] As a reference, the exponent = 1/2 for striped surfaces is also shown as dotted line. Symbols are simulation data.

partial slip

no-slip

symmetric partial slip

FIG. 6. Polar coordinate systems for (a) trapezoidal and (b) triangular profiles of the local slip length. For a mathematical convenience the two coordinates differ in the choice of the origin of the angular variable .

the shear stress behaves as r-1/2 ,24, 25, 38 i.e., the singularity is stronger and creates a source of additional viscous dissipation that reduces effective slip.24 We are now in a position to prove that similar singularities also arise at the border between noslip and linear local slip regions. As above, we develop our theory for a longitudinal configuration only. The solution of Laplace equation, Eq. (12), can be constructed in the vicinity of the edge of slipping regions by using polar coordinates (r, ).24, 56 It is convenient to chose different coordinates for trapezoidal and triangular textures, and to consider these two cases separately. For a trapezoidal slip profile, we use the same coordinates as Refs. 24, 56 (illustrated in Fig. 6(a)). In the case of triangular textures, for symmetry reasons it is more convenient to use the coordinates shown in Fig. 6(b). We begin by studying the trapezoidal slip profile. For the trapezoid, the origin of coordinates is at (ye , ze ) = ( 1 /2, 0), with y = 1 /2 + r cos and z = r sin , 0 , so that the no-slip and the slip regions correspond to = 0 and = . The slip condition (13) becomes r 1, = : u - r z u = r, (32)

or equivalently as ( )
-1

||

0.6

tan( ) = -.

(37)

where is defined as = 2b2 , L ( 1 - 1 - 2 ) (33)

and can vary from 0 to , depending on the local slip, b2 /L, and on the value of 1 + 2 . In the case of stripes, one has 1 + 2 = 1, resulting in at finite b2 /L. The solution of Eq. (12) that satisfies the no-slip boundary condition at = 0 is given by u = cr sin , and the components of the velocity gradient are z u = c r
-1

(34)

The theoretical prediction of Eq. (37) for the behavior of as a function of is shown in Fig. 7 (solid curve). Note that decays monotonically from ( 0) 1 down to ( ) 1/2. In other words, for any finite , we have < 1, and the shear stress becomes singular. As shown in Fig. 7, these theoretical predictions are in excellent quantitative agreement with the simulation results, where we have measured the slip velocity along the surface and then obtained by using Eq. (34). Altogether the theoretical and numerical results confirm that the singularity is weaker for trapezoidal profiles than for stripes, such that the viscous dissipation is smaller. This explains qualitatively the results presented in Fig. 5, where the slip length deviates more strongly from the area-averaged slip length on striped than on trapezoidal textures. In the theoretical analysis of the triangular slip profile ( 1 = 0), it is convenient to shift the origin of the angular variable compared to the trapezoidal case (see Fig. 6), so that z = r cos , y = r sin , | | /2. The even solution of the Laplace equation (12) is u = cr cos ( ), (38)

cos[ (1 - )], sin[ (1 - )].

(35)

and the components of the velocity gradient are then z u = c r
-1

cos [ (1 + )], sin [ (1 + )].

(39)

y u = -c r

-1

(36) y u = -c r
-1

(40)

The velocity at the edge is finite provided > 0, but its gradient becomes singular when < 1. In this case, the two terms in the left-hand side of Eq. (32) are of the same order, r , but the right-hand side is smaller and can be neglected, since r r as < 1. As a result, Eq. (32) can be rewritten, taking into account (34)­(36) as sin( ) - cos[(1 - ) ] = 0,

The slip condition (13) at =± /2 can be rewritten as cr cos - cos 1 + = r, (41) 2 2 with = 2b2 /L. (42)


174708-9

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

Applying the condition < 1, we can neglect the righthand side to obtain = -1 . (43) tan 2 The numerical solution of this equation is included in Fig. 7 (dashed curve), and supported by simulation data. We see that again decays monotonically from ( 0) = 1, but at large the solution of Eq. (43) is = 2
1/2

b2 > L/2, the singularity becomes stronger than for stripes ( > 1/2).

VI. OPTIMIZATION OF THE EFFECTIVE SLIP

1,

i.e., different from expected for a trapezoidal profile, where would tend to 0.5. This implies that for given , the triangular slip profile always gives a stronger singularity compared to a trapezoid slip profile. At small , the singularity is weaker than for a striped surface. A striking conclusion from our analysis is that nonsmooth surface textures with continuous local slip can generate stronger singularities in flow and shear stress than surfaces with a discontinuous, piece-wise constant slip length profile. Indeed, we have > 1/2 for > 2(or b2 > L), i.e., the singularity for a triangular texture becomes stronger than for stripes and should significantly reduce the longitudinal effective slip. This is likely a reason why the effective slip in Fig. 4(b) grows only weakly at very large b2 /L. Coming back to Fig. 5,weremark that the results for triangular profile there correspond to 0.77, i.e., the singularity is weaker in this = 0.5 and case than for stripes with = 0.5. The singularities for trapezoidal profiles may be even weaker: The exponents vary from 0.87 to 0.5. Note however that since the no-slip area remains the key factor determining the effective slip, beff always decreases with 1 . To calculate the exponents for a transverse configuration, one can use the relation between the velocity fields in eigendirections:25 v= 1 ud ud + z 2 z p=- , w=- z ud , 2 y (44)

Based on the results obtained in Secs. IV­V,wenow discuss how to optimize the effective slippage, and what maximum slip lengths may actually be expected, taking into account the actual limitation in surface engineering. Very large effective slip lengths were recently found for sinusoidal profiles b(y) with a no-slip point.26 If we request that the local slip pattern includes a no-slip area, we also find here that the largest values of beff are obtained for triangular profiles, where the no-slip area reduces to just one point. However, real textures always have finite contact area with liquid. The effective slip length of a trapezoidal texture with large b2 /L is close to that of stripes [see Fig. 4(a)]. In this section, we therefore focus on the case of a trapezoidal profile with small 1 , where both b2 /L and the slope of the inclined regions is large. Such a texture should lead to large beff . To find one-dimensional textures which optimize the effective

(a)

25

20

b(y)/L

15

10

5

0 -0.4 -0.2 0 y/L 0.2 0.4

ud , y

(45)

where ud (y, z) = u[y, z, 2b(y)/L] is the velocity field for the longitudinal pattern with a twice larger local slip length. It follows from Eqs. (44) and (45) that the velocity gradients and the pressure in the transverse configuration for both trapezoidal and triangular textures have stronger singularities at the edges of slipping regions than those for the longitudinal configuration: v w , , pr z z
-1

(b)

1.4 1.2 1.0

beff/L

0.8 0.6 0.4 0.2 0.0 -2 10 1 10
-1

,

where = (2 ) < ( ). This leads to larger viscous dissipation. As a result the anisotropy of the slip-length tensor is finite even at small b2 /L and is approximately the same for triangular, trapezoidal, and stripe profiles (see Fig. 5). For example, the eigenvalues shown there for a triangular texture 0.77 and 0.64. Therefore, in the correspond to transverse configuration even for a shallow triangular texture,

FIG. 8. (a) Profiles of local slip length resulting in the same beff . (b) Effective slip lengths for stripe (solid lines and square symbols) and trapezoidal (dotted lines and circular symbols) local slip profiles. The top and bottom lines correspond to beff and beff , respectively.


174708-10

Zhou et al.

J. Chem. Phys. 139, 174708 (2013)

forward slip, we consider profiles of the following form |t | 1 0, 2 b (y ) 2|t |-1 1 = . (46) , 2 < |t | 51 2 91 L -1 1 , 51 < |t | 1 , 1 < 0.1 2 The bottom base of the trapezoid for such textures is wide, 10 1 [see Fig. 8(a)]. For comparison, we also consider stripes - with the same maximal local slip length, b2 /L = 1 1 : b (y ) = - L 1 1 , 0, |t |
1 2 1 2 1 2

lar Structures of New Generations," and by the DFG through SFB TR6 and SFB 985. The simulations were carried out using computational resources at the John von Neumann Institute for Computing (NIC JÝlich), the High Performance Computing Center Stuttgart (HLRS), and Mainz University (MOGON).
1 2 3 4

< |t |

.

(47)

5 6 7 8 9 10 11 12 13 14 15 16 17

Both profiles are characterized by only one parameter 1 . The simulation and numerical results for beff and beff are shown in Fig. 8(b). The effective slip lengths for the trapezoid and stripe profiles are nearly equal for 1 < 0.1. For example, beff for the stripe profile with 1 = 0.05 is only slightly larger than the corresponding value for the trapezoid profile with the same 1 and wide bottom base [see Fig. 8(a)]. Note that such narrow trapezoids with wide bottom base allowed researchers to achieve a very large effective slip in recent experiments.34 The same value of beff can be achieved with a triangular profile with bmax /L = 4. Thus the effect of the inclined regions is relatively small for large slopes, and the main contribution to the effective slip length stems from the no-slip region and the singularity of the shear stress.
VII. CONCLUDING REMARKS

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

In conclusion, we have presented simulation and numerical results and some asymptotic law analysis that allowed us to assess the frictional properties of superhydrophobic surfaces as a function of the texture geometric parameters. We have focused on one-dimensional superhydrophobic surfaces with trapezoidal grooves and systematically investigated the effect of several parameters, such as the maximum local slip and the heterogeneity factor. The slip properties on a trapezoidal surface were compared with two limiting cases: striped and triangular slip length profiles. Our simulation and numerical results suggested that the effective slip length is mainly determined by (i) the area fraction of the no-slip region and (ii) the singularity of the shear stress at the no-slip edge. The first factor is obvious; a smaller fraction of no-slip area results in less friction of the surface, which in turn leads to large slip length. The second effect is less intuitive, especially for a continuous slip profile. Here we demonstrated that the singularity develops for the trapezoidal and triangular local slip profiles, and in the latter case can be even stronger than for discontinuous slip of striped textures. The singularity is stronger for a transverse configuration in comparison to a longitudinal configuration, and as a result, it always contributes to the flow anisotropy. Our results can guide the design of superhydrophobic surfaces for micro/nanofluidics.
ACKNOWLEDGMENTS

42 43

This research was supported by the RAS through its priority program "Assembly and Investigation of Macromolecu-

D. Quere, Rep. Prog. Phys. 68, 2495 (2005). D. Richard and D. Quere, Europhys. Lett. 48, 286 (1999). D. Richard and D. Quere, Europhys. Lett. 50, 769 (2000). P. Tsai, R. C. A. van der Veen, M. van de Raa, and D. Lohse, Langmuir 26, 16090 (2010). O. I. Vinogradova and A. L. Dubov, Mendeleev Commun. 22, 229 (2012). G. McHale, M. I. Newton, and N. J. Schirtcliffe, Soft Matter 6, 714 (2010). L. Bocquet and J. L. Barrat, Soft Matter 3, 685 (2007). J. P. Rothstein, Annu. Rev. Fluid Mech. 42, 89 (2010). O. I. Vinogradova and A. V. Belyaev, J. Phys.: Condens. Matter 23, 184104 (2011). K. Kamrin, M. Bazant, and H. A. Stone, J. Fluid Mech. 658, 409 (2010). M. Z. Bazant and O. I. Vinogradova, J. Fluid Mech. 613, 125 (2008). S. Schmieschek, A. V. Belyaev, J. Harting, and O. I. Vinogradova, Phys. Rev. E 85, 016324 (2012). O. I. Vinogradova, Langmuir 11, 2213 (1995). D. Andrienko, B. DÝnweg, and O. I. Vinogradova, J. Chem. Phys. 119, 13106 (2003). C. Ybert, C. Barentin, C. Cottin-Bizonne, P. Joseph, and L. Bocquet, Phys. Fluids 19, 123601 (2007). D. Maynes, K. Jeffs, B. Woolford, and B. W. Webb, Phys. Fluids 19, 093603 (2007). T. V. Nizkaya, E. S. Asmolov, and O. I. Vinogradova, "Flow in channels with superhydrophobic trapezoidal textures," Soft Matter (published online). A. D. Stroock, S. K. Dertinger, G. M. Whitesides, and A. Ajdari, Anal. Chem. 74, 5306 (2002). A. Ajdari, Phys. Rev. E 65, 016301 (2001). E. Lauga and H. A. Stone, J. Fluid Mech. 489, 55 (2003). A. V. Belyaev and O. I. Vinogradova, J. Fluid Mech. 652, 489 (2010). C. O. Ng,H.C.W.Chu, andC.Y.Wang, Phys. Fluids 22, 102002 (2010). J. Zhou, A. V. Belyaev, F. Schmid, and O. I. Vinogradova, J. Chem. Phys. 136, 194706 (2012). E. S. Asmolov, J. Zhou, F. Schmid, and O. I. Vinogradova, Phys. Rev. E 88, 023004 (2013). E. S. Asmolov and O. I. Vinogradova, J. Fluid Mech. 706, 108 (2012). E. S. Asmolov, S. Schmieschek, J. Harting, and O. I. Vinogradova, Phys. Rev. E 87, 023005 (2013). S. C. Hendy and N. J. Lund, Phys. Rev. E 76, 066313 (2007). C. Cottin-Bizonne, C. Barentin, E. Charlaix, L. Bocquet, and J.-L. Barrat, Eur. Phys. J. E 15, 427 (2004). N. V. Priezjev, A. A. Darhuber, and S. M. Troian, Phys. Rev. E 71, 041608 (2005). N. V. Priezjev, J. Chem. Phys. 135, 204704 (2011). N. Tretyakov and M. MÝller, Soft Matter 9, 3613 (2013). W. Li, X. S. Cui, and G. P. Fang, Langmuir 26, 3194 (2010). T. Tanaka, M. Morigami, and N. Atoda, Jpn. J. Appl. Phys. 32, 6069 (1993). C.-H. Choi, U. Ulmanella, J. Kim, C.-M. Ho, and C.-J. Kim, Phys. Fluids 18, 087105 (2006). A. V. Belyaev and O. I. Vinogradova, Soft Matter 6, 4563 (2010). E. S. Asmolov, A. V. Belyaev, and O. I. Vinogradova, Phys. Rev. E 84, 026330 (2011). F. Feuillebois, M. Z. Bazant, and O. I. Vinogradova, Phys. Rev. Lett. 102, 026001 (2009). M. Sbragaglia and A. Prosperetti, Phys. Fluids 19, 043603 (2007). J. HyvÄluoma and J. Harting, Phys. Rev. Lett. 100, 246001 (2008). A. Steinberger, C. Cottin-Bizonne, P. Kleimann, and E. Charlaix, Nature Mater. 6, 665 (2007). E. Karatay, A. S. Haase, C. W. Visser, C. Sun, D. Lohse, P. A. Tsai, and R. G. H. Lammertink, Proc. Natl. Acad. Sci. U.S.A. 110, 8422 (2013). A. S. Haase, E. Karatay, P. A. Tsai, and R. G. H. Lammertink, Soft Matter 9, 8949 (2013). P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. 19, 155 (1992).


174708-11
44 45 46 47 48 49

Zhou et al.
50 51 52 53 54 55 56

J. Chem. Phys. 139, 174708 (2013) O. I. Vinogradova and G. E. Yakubov, Langmuir 19, 1227 (2003). O. I. Vinogradova, K. Koynov, A. Best, and F. Feuillebois, Phys. Rev. Lett. 102, 118302 (2009). C. Cottin-Bizonne, B. Cross, A. Steinberger, and E. Charlaix, Phys. Rev. Lett. 94, 056102 (2005). L. Joly, C. Ybert, and L. Bocquet, Phys. Rev. Lett. 96, 046101 (2006). O. I. Vinogradova, Int. J. Min. Process. 56, 31 (1999). J. R. Philip, J. Appl. Math. Phys. 23, 353 (1972). C. Y. Wang, Phys. Fluids 15, 1114 (2003).

P. Espaßol and P. Warren, Europhys. Lett. 30, 191 (1995). R. D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997). T. Soddemann, B. DÝnweg, and K. Kremer, Phys. Rev. E 68, 046702 (2003). J. Smiatek, M. Allen, and F. Schmid, Eur. Phys. J. E 26, 115 (2008). J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 54, 5237 (1971). H. Limbach, A. Arnold, B. Mann, and C. Holm, Comput. Phys. Commun. 174, 704 (2006).