Roughness on liquid-infused surfaces induced by capillary waves
UUnder consideration for publication in J. Fluid Mech.
Roughness on liquid-infused surfaces inducedby capillary waves
Johan Sundin † , St´ephane Zaleski and Shervin Bagheri Linn´e FLOW Centre, Dept. Engineering Mechanics, KTH, Stockholm, Sweden Sorbonne Universit´e & CNRS, UMR 7190, Institut Jean Le Rond d’Alembert, Paris, France(Received xx; revised xx; accepted xx)
Liquid-infused surfaces (LIS) are a promising technique for reducing friction,fouling and icing in both laminar and turbulent flows. Previous work has demon-strated that these surfaces are susceptible to shear-driven drainage. Here, wereport a different failure mode using direct numerical simulations of a turbulentchannel flow with liquid-infused longitudinal grooves. When the liquid-liquidsurface tension is small and/or grooves are wide, we observe traveling-waveperturbations on the interface with amplitudes larger than the viscous sublayer ofthe turbulent flow. These capillary waves induce a roughness effect that increasesdrag. The generation mechanism of these waves is explained using the theorydeveloped by Miles for gravity waves. Energy is transferred from the turbulentflow to the LIS provided that there is a negative curvature of the mean flow atthe critical layer. Given the groove width, the Weber number and an estimate ofthe friction Reynolds number, we provide relations to determine whether a LISbehaves as a smooth or rough surface in a turbulent flow.
1. Introduction
A protective and functional surface coating can be created by lubricating atextured surface with an appropriate liquid. The drag reducing properties ofthese liquid-infused surfaces (LIS) have been explored recently both numerically(Fu et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b J. Sundin, S. Zaleski and S. Bagheri contact angle at the liquid-liquid-solid interface. In particular, these surfaces mayexperience shear-driven drainage of the infused liquid, but this can be mitigated,for example with chemical patterning (Wexler et al. et al. τ ≈ a + ≈ − w + ) and a Weber number (We + ), bothnormalised with the viscous length scale. This design map is obtained from thecritical-layer theory and provides a means to design LIS that can be predicted toachieve a balance between performance (large w + ) and stability (smooth domain).
2. Numerical methods and configuration
We consider a fully developed turbulent open channel flow. The flow domain,shown schematically in fig. 1b, has the size ( L x , L y , L z ) = (6 . h, h + k, . h ), where x, y and z correspond to the streamwise, wall-normal and spanwise directions,respectively, and h is the half-channel height. The streamwise-aligned grooves atthe bottom wall have a height k , a width w and square cross-section, k = w .The fluid-solid ratio is set to 0 .
5. The infused and external fluids have the samedensity ρ i = ρ ∞ = ρ , but different viscosities ( µ i and µ ∞ ). We have used groovesof width w +0 = 18. Throughout this manuscript, +0 refers to normalisation usingthe friction velocity of a regular smooth wall (nominal wall units).We impose a constant mass flow rate through a uniform pressure gradientover 0 < y < h , where y = 0 corresponds to the crest of the texture so thatRe b = ρhU b /µ = 2820, giving Re τ = ρhu τ /µ ≈ b and Re τ areReynolds numbers based on bulk velocity U b and friction velocity u τ , respectively.Our simulations allow for a moving liquid-liquid-solid contact line with adynamic contact angle different from the static value, which is θ = 45 ° withrespect to the infused liquid. Our method also allows for interface deformation, oughness on liquid-infused surfaces induced by capillary waves We + − − − w + a ) η λ x λ z w c ( c ) µ i µ ∞ h w . h . hxy z free slip (symmetry) B.C. w ( b ) Figure 1: ( a ) Design map for LIS, spanned by We + and w + . Smooth and rough regions areseparated using eqs. (3.12), (3.14) and (3.15). Values from simulations with w +0 = 18 areincluded, with symbols referring to µ i /µ ∞ = 1 ( ◦ ), µ i /µ ∞ = 0 . (cid:3) ) and µ i /µ ∞ = 2 ( ♦ ) andcolours to We = 100 (blue), We = 150 (red) and We = 200 (yellow). Shown are also valueswith wider grooves, w +0 = 36, with We = 50 (green circle) and We = 100 (blue circle). Theasymptotic relations (3.16) and (3.17) are shown with (— — —) and (— · —), respectively. ( b )Sketch of the channel configuration and ( c ) of a wave on a groove. The infused liquid is shownin green and the solids in gray. which is typically quantified by the Weber number, defined as We = ρU h/γ ,where γ is the surface tension. We have simulated LIS for We = 100, 150 and200 and viscosity ratios µ i /µ ∞ = 0 .
5, 1 and 2. The Weber number in wall unitsis We + = ρu τ δ ν /γ = µ ∞ u τ /γ , where δ ν = µ ∞ / ( ρu τ ) is the viscous length scale.The numerical configuration described above corresponds to an infused liquidconsisting of some alkane (dynamic viscosities similar to water (Van Buren &Smits 2017)), a water channel with h = 0 . U b = 1 m/s. This results inRe b = 5000, which is close to our simulations. A typical surface tension γ = 50mN/m then results in We = 100.The code used for the simulations is based on the PArallel, Robust, InterfaceSimulator (PARIS), which employs a VOF method for the multi-phase description(Aniszewski et al. et al. N x , N y , N z ) = (256 , , h/U b , and statistics were collectedover a time of 500 h/U b .
3. Results
Dependency of drag on Weber number
Figure 2(a) shows an instantaneous snapshot of the liquid-liquid interface, viewedfrom top, for We = 100 and µ i = µ ∞ . There are oscillations on the interface, J. Sundin, S. Zaleski and S. Bagheri ( b )( a ) Figure 2: Top view of interfaces (green) and surface (gray) at one instance for µ i /µ ∞ = 1 andWe = 100 ( a ) and We = 200 ( b ). The flow is from left to right. The complete domain is shown. due to the finite surface tension, but they remain small. The deformation of theinterface increases with We, however. For We = 200, significantly larger wavesdevelop on the interface as shown in fig. 2(b).The consequences of the waves on the overlying flow can be quantified byDR = c f − c f c f , (3.1)where c f = 2 τ w / ( ρU b ) is the friction coefficient and τ w is the total stress at thecrest plane of the surface (computed from the pressure gradient). Here, c f is thecoefficient of a regular smooth wall at y = 0. For We = 100 and We = 200, weobtained DR = 0 .
09 and DR = − .
04, respectively. In other words, the capillarywaves observed at We = 200 increase frictional drag compared to a smooth andhomogeneous surface and have therefore induced failure of the LIS.Figure 3a shows DR for We = { , , } and viscosity ratios µ i /µ ∞ = { . , , } as a function of the apparent slip length, b +0 . The slip length b is thedistance at which the mean velocity would be zero if linearly extrapolated at thecrests of the surface. It is largely unaffected by changes in We, but it increases withdecreasing viscosity ratio. In fact, the slip lengths extracted from our numericalsimulations with w + ≈
18 are well approximated by the slip lengths obtained bysolving the Stokes equations for a single groove exposed to unit shear (Sch¨onecker et al. w + ≈
9) in turbulent flows by Fu et al. (2017); Arenas et al. (2019).When the interface is perfectly flat (We = 0), the drag reduction can be relatedto slip length as DR ≈ b +0 b +0 + U +0 b . (3.2)This relation – shown in 3a (black line) – can be obtained by neglecting changes inthe Reynolds shear stress above a smooth wall (Rastegari & Akhavan 2015). Weobserve from fig. 3a that for We = 100 (blue) and We = 150 (red), there is a dragreduction (DR >
0) for all three viscosity ratios. Moreover, the deviations from(3.2) are small, confirming that the drag reduction mechanism is indeed slippage.In contrast, the deviations from (3.2) are significant for We = 200 (yellow), wherewe observe a drag increase (DR <
0) for µ i /µ ∞ = 1 and 2 and a DR close tozero for µ i /µ ∞ = 0 .
5. The corresponding mean velocity and velocity fluctuationsreflect the increase in drag, and these are described in the SM (sec. S2).The waves formed on the interface at We = 200 are sufficiently large to causeroughness effects. Interface height profiles at different times are shown in fig. 4b,together with amplitudes of a wave in its initial stage (fig. 4c, yellow). Wave oughness on liquid-infused surfaces induced by capillary waves b +0 D R -0.1-0.0500.050.10.150.2( a ) µ i / µ ∞ − b / ( w ) b ) Figure 3: Drag reduction as a function of the slip length ( a ) and slip length normalised by thepitch as a function of viscosity ratio ( b ). Symbols and colours are the same as in fig. 1a. In ( b ),the points for different We are almost on top of each other. Analytical relations are shown byblack lines, —— Rastegari & Akhavan (2015) ( a ) and —— Sch¨onecker et al. (2014) ( b ). amplitudes of a + > a + <
1) (fig. 4a) and show a significantlysmaller growth rate (fig. 4c, blue). The exponential growth rate (fig. 4c, dashed)is an indication of a linear instability. In the next section, we provide evidence ofa critical-layer instability (Miles 1957), where energy is transferred to the waveperturbation from the turbulent flow.3.2.
Conditions for phase speed and growth rate of capillary waves
We assume a small perturbation on the liquid-liquid interface of the form η = Ae ik x ( x − ct ) cos( k z z ) , (3.3)where z = 0 is located in the centre of the groove. As illustrated in fig. 1c, η isthe height of the interface, A is the initial wave amplitude, c is a complex wavespeed and t is time. Moreoever, k x = 2 π/λ x and k z = 2 π/λ z are streamwiseand spanwise wavenumbers, respectively. The spanwise wavelength can have amaximum value of λ z = 2 w , due to the finite width of the grooves. This valuecan be seen to dominate in the snapshots of fig. 2, as most waves only have onecrest or one trough in the spanwise direction. We also observe from fig. 2 thatstreamwise wavelengths are generally similar to, or larger than, λ z . This three-dimensionality implies that both spanwise and streamwise curvature contributeto the capillary pressure of a wave, ∆p cap = p +0 − p − = γ (cid:18) ∂ ∂x + ∂ ∂z (cid:19) η = − γk η. (3.4)Here, k = p k x + k z and p +0 ( p − ) is the pressure above (below) the interface.Next, we consider a wall-normal velocity disturbance v ( x, y, z, t ) on the tur-bulent mean flow U ( y ) with the same wave-form as η . If we neglect viscous andnonlinear effects, the amplitude ˆ v ( y ) is governed by the Rayleigh equation (SM,sec. S3 A), 1 k ˆ v − (cid:20) U − c ) 1 k U (cid:21) ˆ v = 0 , (3.5) J. Sundin, S. Zaleski and S. Bagheri y + -101( a ) x/h y + -505( b ) tU b /h a + We = 100We = 200exp. fi t ( c ) Figure 4: Instantaneous interface heights in the center-line of a groove for µ i /µ ∞ = 1; ( a )We = 100 for 1 x b ) We = 200 for 3 x
6. The 5 profiles are separated by ∆t = 0 . h/U b . Note the difference in vertical scale. In ( c ), the wave amplitude developing at x/h = 4 . a ) and at x/h = 3 . b ) are shown. The phase speeds in ( a ) and ( b ) can beestimated as c + ≈
14 and c + ≈
10, respectively. where denotes a derivative with respect to y . The velocity perturbation mustvanish at infinity and satisfy the kinematic condition at the interface, v/ ( U − c ) = ik x η . The equation for the pressure amplitude, ˆ p ( y ), corresponding to eq. (3.5) isˆ pρ = − i k x k [( U − c )ˆ v − U ˆ v ] . (3.6)Our aim is to find an approximate solution to the equations (3.4-3.6) in order todetermine the phase speed < ( c ) (real part) and growth rate = ( k x c ) (imaginarypart) of the interface perturbation (3.3).Miles (1957) formulated a similar set of equations for describing wind-inducedwater waves, where gravity – instead of capillary – balances fluid pressure. Hesuggested the following approximate solution for v , v = ik x η ( U − c ) e − ky , y > . (3.7)This expression, which satisfies the boundary conditions at y → ∞ , impliesthat 1 /k is the relevant length scale over which v decreases. The assumptionof exponential decay can also be used in the grooves v = − ik x cηe ky , y < . (3.8)Here, we have assumed that the grooves are sufficiently deep such the velocityperturbation is nearly zero at the bottom of the groove. With a depth w , kw >k z w > π/ (2 w ) w = π , and since e − π (cid:28)
1, the assumption is valid for ourconfiguration. We have also neglected U and its derivative inside the groove.Inserting (3.8) into (3.6) results in the following expression for the pressureimmediately below the interface ( y → − ), p − = ρ k x k c η. (3.9)Similarly, by inserting (3.7) into (3.6) the pressure just above the interface is p +0 = ( α + iβ ) ρU k x k η, (3.10)where α and β are real constants and U is an arbitrary reference velocity. It isshown in SM (sec. S3 C) that the parameter α can be decomposed into two parts, oughness on liquid-infused surfaces induced by capillary waves λ x /w = 2 π / ( k x w ) c + w We = 100We = 150We = 200We = 200, k z = 0 ( a ) ky c − − − β − − − ( b ) Figure 5: ( a ) The free phase speed c +0 w for We = { , , } when k z = π/w (i.e. λ z /w = 2).For We = 200, the phase speed of a two-dimensional wave ( k z = 0) is also shown. ( b ) The growthrate coefficient β versus ky c , showing a fast decrease in β for ky c ≈ α = α + α , where α corresponds to eq. (3.9) and α incorporate the remainingcontributions from the slip velocity and the shear. Using this decomposition andinserting (3.9) and (3.10) into (3.4), we obtain (SM, sec. S3 D), c = c w (cid:18) α + iβ ) U c w + . . . (cid:19) . (3.11)Here c w = p γk / (2 ρk x ) is the free phase speed, i.e. the speed of a capillary wavewithout forcing from the overlying flow.The free phase speed is shown in fig. 5a as a function of λ x /w for We = { , , } . We note that two-dimensional capillary waves ( k z = 0) have aphase speed that monotonically decreases with λ x . However, for LIS, there is aminimum phase speed due to the finite spanwise wavelength. This minimum isapproximately (for analytical expression see SM, sec. S3 D) c + w, min ≈ r π We + w + . (3.12)For We = 200 (and µ i /µ ∞ = 1)¸ c + w, min = 7 .
0, which is slightly lower than thephase speed of the wave shown in fig. 4b. One may use c + w, min as a lower boundof the actual phase speed of capillary waves on LIS.We now turn our attention to the imaginary part of (3.11) to approximate thegrowth rate of the instability. As shown by Miles (1957) – and repeated in the SM(sec S3 E) – one may integrate the Rayleigh equation and evaluate the pressureequation (3.6) at the interface to find, β = − π (cid:12)(cid:12)(cid:12)(cid:12) v c k x ηU (cid:12)(cid:12)(cid:12)(cid:12) k U c U c , (3.13)where the subscript c denotes values at the position of the critical layer y c . Inorder for an infinitesimal wave to have a positive growth rate, i.e. β >
0, a firstrequirement is that U c <
0, i.e. negative curvature at the critical layer. This issatisfied if the critical layer is outside of the viscous sublayer.A second requirement for β > v c in eq. (3.13) is non-zero at the criticallayer. The approximate solution of v in eq. (3.7) implies however that v is zeroat the critical layer. As shown in SM (sec. 3 E), one may transform the conditionfor positive growth rate to an integral form to estimate v in the vicinity of the J. Sundin, S. Zaleski and S. Bagheri critical layer. By further assuming a logarithmic mean velocity profile and settingthe reference velocity to U = u τ /κ (where κ is the von-Karman constant), onemay evaluate the expression for β , and obtain what is shown in fig. 5b.We observe from fig. 5b that when ky c > /
2, then β < .
02, which results invery slow growing waves, whereas when ky c < /
2, we have β > .
8, resultingin a factor 40 or more faster growth. The gray region in fig. 5b marks the range1 / < ky c < / k > k z and the upper limit of λ z is 2 w , we may formulate bounds for theposition of the critical layer. When y + c . w + π , (3.14)the growth rate can be expected to be significant, in contrast to when y + c & w + π , (3.15)for which the growth is negligible.Equations (3.12), (3.14) and (3.15) provide relationships between We + , w + and c + w, min , y + c, max that can be confirmed by our simulations. Equation (3.12) states thata large We + and/or w + give a small phase speed. This is observed qualitatively byfollowing the travelling waves on the interface in fig. 4(a,b). More quantitatively,the space-time correlations of the interface height for We = 100 and We = 200give c + = 15 . c + = 10 .
5, respectively (SM fig. S8).Compared to We = 100, the lower phase speed for We = 200 results in a lowerposition of the critical layer. When the height of the critical layer approachesthe interface and satisfies eq. (3.14), the growth rate coefficient β of the waves(eq. 3.13) is large. This is confirmed by our simulations, where we observe infig. 4c that both the growth rate and interface amplitudes are larger for We = 200compared to We = 100.3.3. Implications for the design of LIS
The conditions (3.12) (3.14) and (3.15) can be used as design criteria for LIS.One may expect a high-performing LIS by choosing a groove width and a surfacetension of the infused liquid such that – for relevant friction Reynolds numbers– the design falls within the smooth region of fig. 1a. This region is defined by(We+ , w + ), where y + c > . w + /π and thus from eq. (3.15) very small growthrates of capillary waves are predicted. Conversely, the rough region in fig. 1ashows (We+ , w + ), where y + c . w + /π and therefore waves will amplify rapidly.Here, we may expect either very low-performing LIS or even drag-increasing LISdue to roughness effects. In between the smooth and rough domains, we showin fig. 1a a transitional region (gray), which corresponds to 0 . w + /π y + c . w + /π . Here, we cannot predict if the resulting waves induce roughness effectsusing our analytical approach. It should be mentioned that the boundaries of thetransitional region in fig. 1a are determined in three steps; given w + , determine y + c from (3.14) (lower boundary) or (3.15) (upper boundary); given y + c determine c + (phase speed) from U + ( y + c ) = c + , where U + ( y ) is a turbulent mean profile ofa smooth wall and finally; given c + assume c + ≈ c +0 w, min and determine We + from(3.12) (or the exact coefficient of (3.12) given in the SM).Figure 1a also shows scaling laws between smooth and rough regions. For small oughness on liquid-infused surfaces induced by capillary waves w + (and thus y + c ) we may assume that the critical layer velocity is U + c ≈ y + c .This is acceptable right above the viscous sublayer where the mean flow hassome curvature. Then eq. (3.12) gives that the height of the lowest critical layeris y + c ≈ q π/ (We + w + ). By assuming that y + c ∼ w + /π , we obtain w + ∼ (We + ) − / , (3.16)which is shown with dashed solid line in fig. 1a. It is observed that this asymptoticrelation represent a reasonable scaling law for w + . w + , away from the viscous sublayer, we assume U + = 1 /κ log( y + ) + B , where B is a constant. This gives a nonlinear relation1 / ( p We + w + ) ∼ /κ log( w + /π ) . (3.17)This curve is shown in fig. 1a with a dashed-dot line, and provides a scaling ofthe neutral curve for w + & w + increases, there needs to be rapiddecrease of We + to remain in the smooth region. For example, for w + ≈ + ≈ · − , which corresponds to We ≈
10. This is relevant fordrag reduction, since the width of the grooves should be maximised for a givensurface tension to optimize DR (see fig. 3), but without entering the rough zonein fig. 1a. Note that for a fixed geometry, increasing the flow speed, and thereby u τ , increases both We + and w + , so that the design needs to be made for thelargest flow speed the surface is exposed to.Finally, in fig. 1a, the values of our numerical simulations are shown withsymbols. In addition to the simulations at w + = 18 we also show points (circles)for larger grooves of width w +0 = 36 (using µ i /µ ∞ = 1). For these grooves, therewas a drag reduction by 17% for We = 50 (green circle), whereas for We = 100(blue circle), the drag reduction was lowered to 2% and we observed large waves.This implies that the growth rate rapidly increases between these two cases asthey fall in the transitional zone in fig. 1a.
4. Conclusions
We have explored the behaviour of LIS in a turbulent channel flow with squarelongitudinal grooves for Re τ ≈ w + and the Weber number We + , as illustrated in fig. 1a. It shouldalso be noted that these non-dimensional parameters depend on the flow speed.Using an analytical analysis, we have provided scaling laws and design criteria forrobust drag-reducing LIS. Specifically, the relations show how to achieve a balancebetween large groove widths (enhancing drag-reduction) and high surface tensionof the infused liquid (enhancing stability) for different flow speeds. Declaration of interests.
The authors report no conflict of interest. J. Sundin, S. Zaleski and S. Bagheri
REFERENCESAniszewski, W., Arrufat, T., Crialesi-Esposito, M., Dabiri, S., Fuster, D., Ling,Y., Lu, J., Malan, L., Pal, S., Scardovelli, R. & others arXiv preprint arXiv:2012.11744 . Arenas, I., Garc´ıa, E., Fu, M. K., Orlandi, P., Hultmark, M. & Leonardi, S.
J. Fluid Mech. , 500–525.
Arrufat, T., Crialesi-Esposito, M., Fuster, D., Ling, Y., Malan, L., Pal, S.,Scardovelli, R., Tryggvason, G. & Zaleski, S.
Comput. Fluids ,104785.
Boomkamp, P. A. M. & Miesen, R. H. M.
Int. J. Multiph. Flow , 67–88. Cartagena, E. J. G., Arenas, I., Bernardini, M. & Leonardi, S.
Flow Turbul. Combust. (4), 945–960.
Epstein, A. K., Wong, T., Belisle, R. A., Boggs, E. M. & Aizenberg, J.
Proc. Natl.Acad. Sci. (33), 13182–13187.
Fu, M. K., Arenas, I., Leonardi, S. & Hultmark, M.
J. Fluid Mech. , 688.
Fu, M. K., Chen, T., Arnold, C. B. & Hultmark, M.
Exp. Fluids (6), 100. Kim, P., Wong, T., Alvarenga, J., Kreder, M. J., Adorno-Martinez, W. E. &Aizenberg, J.
ACS nano (8), 6569–6577. Legendre, D. & Maglio, M.
Comput. Fluids , 2–13.
Ling, H., Katz, J., Fu, M. & Hultmark, M.
Phys. Rev. Fluids (12),124005. Miles, J. W.
J. Fluid Mech. (2),185–204. Rastegari, A. & Akhavan, R.
J. Fluid Mech. . Sch¨onecker, C., Baier, T. & Hardt, S.
J. Fluid Mech. , 168–195.
Seo, J., Garc´ıa-Mayoral, R. & Mani, A.
J. FluidMech. , 45–85.
Van Buren, T. & Smits, A. J.
J. Fluid Mech. , 448–456.
Wang, P., Lu, Z. & Zhang, D.
Corros. Sci. ,159–166. Wexler, J. S., Jacobi, I. & Stone, H. A.
Phys. Rev. Lett. (16), 168301.
Weymouth, G. D. & Yue, Dick K.-P.
J. Comput. Phys. (8), 2853–2865.
Wong, T., Kang, S. H., Tang, S. K. Y., Smythe, E. J., Hatton, B. D., Grinthal, A.& Aizenberg, J.
Nature (7365), 443–447. upplementary Material: Roughness on liquid-infused surfaces induced by capillarywaves
Johan Sundin , St´ephane Zaleski and Shervin Bagheri Linn´e FLOW Centre, Dept. Engineering Mechanics,Royal Institute of Technology, SE-100 44 Stockholm, Sweden Sorbonne Universit´e & CNRS, UMR 7190, Institut Jean Le Rond d’Alembert, F-75005 Paris, France
S1. DETAILS OF NUMERICAL METHODS
In this section, we describe the schemes of the numerical code that was used to solve the two-component flowover and inside the textured surfaces. For momentum convection, a second-order central scheme was used, whereasa third-order Runge-Kutta scheme was used for time integration [1]. The equation for the pressure was solved witha fast Fourier transform (FFT) solver (FFTW) and handling domain decomposition with the 2DECOMP&FFTlibrary. To describe solids, we used a grid-aligned immersed-boundary method [2]. With this method, the edge ofthe solid is located at the edge of the cells. For a staggered grid, this means that the top of the solid ridges are atthe location of the wall-normal velocity nodes.To impose a contact angle at the liquid-liquid-solid contact line, we specified the height function of the first ghostlayer [3]. We used a dynamic model of the contact angle based on hydrodynamic theory [4] – adapted to VOF[5] – , together with a no-slip velocity condition. As shown in [5], the dynamic contact angle also improves gridindependence. The method is explained further in sec. S1 A. It should be noted that the ridge corners’ ghost cellswere set to always be interface cells, and that the contact angle was imposed when the interface moved to adjacentcells. Depinning occured for the higher Weber numbers and also for θ = 90 ° .The turbulent flow was validated by comparing the mean flow and tje r.m.s. velocities for a full channel withsmooth walls to data from ref. [6]. The results are shown in fig. S1, having for the mean velocity a deviation of 0 . y + < y + ! U + a ) y + u +rms w +rms v +rms ( b ) y + ! h u v i + c ) FIG. S1: Simulation of a turbulent channel flow with smooth walls at Re τ ≈ a ), r.m.s. ve-locities ( b ), and Reynolds shear stress ( c ). The spanwise velocity r.m.s. value is represented by w +rms for convenience,but elsewhere w refers to the groove width. —— PARIS; — · — PARIS open channel; — — — Lee & Moser (2015)[6]The description of the streamwise liquid filled grooves was validated against the analytical expressions ofSch¨onecker et al. [7] for the slip length in a laminar flow. These expressions have recently also been used as areference for turbulent data of LIS simulations [8]. For these tests, we use a computational box corresponding to aunit cell of the surface, with height three times the height of the groove. On the top boundary, a constant shearwas applied in the streamwise direction. A similar setup was recently used to investigate the robustness of LISwith spanwise grooves [9]. Having 50 cells in the grooves, the errors were less than 4% for µ i /µ ∞ = 0 . , µ i /µ ∞ = 1, We = 100 and θ = 45 ° , where another grid with( N x , N y , N z ) = (384 , , a r X i v : . [ phy s i c s . f l u - dyn ] F e b i = ! ! b = ( w ) Sch B onecker et al. (2014)PARIS, N y = 150 FIG. S2: Comparison of slip length with the expressions of Sch¨onecker et al. (2014). N y refers to the total numberof grid cells in the wall-normal direction, with N y / y + U + a ) y + ! U + b ) y + u + r m s u +rms w +rms v +rms ( c ) FIG. S3: Grid refinement study: —— regular grid; — — — refined grid. The spanwise velocity r.m.s. value isrepresented by w +rms for convenience. A. Dynamic contact angle
Following ref. [5], the contact angle can be found from the equation G ∗ ( θ num ) = G ∗ ( θ stat ) + Ca ln (cid:18) ∆ / λ (cid:19) , (S1)where θ num is the numerically imposed contact angle, θ stat is the static angle, Ca = U cl √ µ i µ ∞ /γ , with contactline velocity U cl , λ is a microscopic length scale corresponding to an effective slip length of the contact line, heretaken to be constant, and ∆ is the wall-normal cell height. The angle θ stat is in this manuscript denoted only θ forsimplicity. The contact line velocity is measured half a cell above the wall. The function G ∗ ( θ ) is a monotonicallyincreasing function. It is defined as G ∗ ( θ ) = √ qG ( θ ), where q = µ ∞ /µ i and G ( θ ) is the original function derivedby ref. [10], with the extra factor for increased symmetry, G ∗ ( θ ) = Z θ f ∗− ( θ , q )d θ , (S2)with f ∗− ( θ, q ) = q . ( θ − sin θ )[( π − θ ) + sin θ cos θ ] + q − . (( π − θ ) − sin θ )[ θ − sin θ cos θ ]2 sin θ [(2 − q − q − ) sin θ + qθ + q − ( π − θ ) + 2 θ ( π − θ )] . (S3)When written in this form, it is apparent that f ∗ ( θ, q ) = f ∗ ( π − θ, q − ), which is a necessary requirement, sinceeq. (S1) should be independent of which of the two components we consider. This is identical to the hydrodynamictheory of an apparent angle model suggested in ref. [10], where θ num is the apparent angle. Legendre and Maglio [5]3 t= = R = R i n i t ! a ) t= = R = R i n i t ! -0.3-0.25-0.2-0.15-0.1-0.050( b ) FIG. S4: Radius at wall ( y = 0) of droplet spreading from 90 ° to 60 ° ( a ), and contracting from 90 ° to 120 ° ( b ), forthree different wall-normal cell sizes: —— N y = 32; — — — N y = 64; — · — N y = 128. The viscosity ratios are0 . λ should be on the order of the the physical slip length (typically on the order of 1 nm), andmay or may not be used in combination with a numerically applied slip length. We use a no-slip condition, with anumerical value of λ = 2 · − h , that with h ≈ . f ∗− ( θ, N ) was integrated with the trapezoidal rule. For numerical reasons, the valuehas been limited to 30 ° ≤ θ num ≤ ° .We evaluated the grid-independence of the contact line model by looking at droplet spreading for three differentcell sizes. A three dimensional half-sphere droplet of radius R init = 0 . L x , L y , L z ) =(2 , , x - and z -directions, and no-slip and shear-free boundary conditionsin the negative and positive y -direction, respectively. The initial radius sets the length scale of the problem. Thecontact angle was initially 90 ° , when the static contact was set to θ stat = 60 ° , so that the droplet started to spread.The droplet had a density ρ = 1, viscosity µ = 0 .
25 and surface tension γ = 7 .
5. Using the density and the viscositywe can define a time scale τ = ρR /µ . The surrounding fluid had the same density, and the viscosity ratio waschanged by varying the viscosity of the surrounding fluid. Viscosity ratios of both 0 . N x = 64 and N z = 64, and thenumber of cells of the wall-normal direction was varied between N y = 32, N y = 64 and N y = 128. The microscopiclength scale was set to λ = 2 · − . Resulting spreading radii are shown in fig. S4a. The differences between thecurves are a few percentage, considered small enough for this study. Corresponding droplet contraction tests, withcontraction from 90 ° to 120 ° , are shown in fig. S4b. S2. TURBULENCE STATISTICS
In this section, we show the turbulent statistics corresponding to the simulations We = { , , } andviscosity ratios µ i /µ ∞ = { . , , } . Mean velocity profiles are plotted in fig. S5a in wall units and fig. S5b in outerunits. As can be seen from the plot in wall units, the centerline velocity is increased by the LIS for low We. Thisis related to the drag reduction achieved at these We, whereas the opposite occurs for the cases of drag increase[11]. A reduction of the centerline velocity is followed by a decrease in the streamwise r.m.s. component, and anincrease of the wall-normal and the spanwise fluctuations, shown in fig. S5c. Over all, drag increase results inan increased isotropy of the velocity fluctuations. Increased wall-normal velocity fluctuations at the surface has astrong connection to the increase in drag [12, 13]. It implies stronger flow ejections, caused by the roughness thewaves impose.The pressure r.m.s. profiles are shown in fig. S5d. Outside the grooves, the profiles for We = 100 and We = 150are similar to the smooth wall profile, but are higher inside the grooves. For We = 200, when there are large waves,the pressure fluctuations are increased also outside the grooves. This could be due to the increased roughness, butalso due to formation of droplets.In the streamwise r.m.s. component, fig. S5c, there is a peak at y = 0, but not for the spanwise nor for thewall-normal. This can be related to the dispersive stresses from the solid structures, which are the r.m.s. of theroughness-coherent velocity, u RC . Due to the symmetry in the streamwise direction, the roughness-coherent velocityis the velocity field averaged in the streamwise direction and time. A peak is then created because the streamwisevelocity component on average is larger over the interface at the grooves than over the ridges. The peak height differsbetween the viscosity ratios, but not so much between the different Weber numbers. This reflects the behaviour ofthe mean velocity at the interface, fig. S5a for small y + , as well as the slip length, fig. 3b. It was seen that both4the spanwise and wall-normal roughness-coherent components were close to zero. The streamwise part is shown infig. S6 for We = 100. y + ! U + a ) y=h -0.2 0 0.2 0.4 0.6 0.8 1 U = U b b ) u +rms w +rms v +rms y + u + r m s c ) y + p + r m s d ) FIG. S5: Mean velocity profiles in wall units ( a ), and outer units ( b ) together with r.m.s. velocities ( c ) andr.m.s. pressure ( d ). The grey area represent the grooves. Colours represent different We as in fig. 3 (We = 100-blue, We = 150-red and We = 200-yellow), and the different lines represent different viscosity ratios, µ i /µ ∞ = 0 . µ i /µ ∞ = 1 (——) and µ i /µ ∞ = 2 (— · —). Also shown are statistics for a smooth wall ( · · · · · · ). Thespanwise velocity r.m.s. value is represented by w +rms for convenience, but elsewhere w refers to the groove width. . y + -20 -10 0 10 20 30 u + R C ; r m s FIG. S6: Streamwise dispersive stress for We = 100, θ = 45 ° and µ i /µ ∞ = 0 . µ i /µ ∞ = 1 (——) and µ i /µ ∞ = 2 (— · —).5 S3. THE MILES INSTABILITY ON GROOVESA. Governing equations
In this section, we look at the governing equations for the Miles instability. We start from the equations for ainfinitesimal disturbance on a baseflow U = U ( y ) [14], ∂u∂t + U ∂u∂x + vU = − ρ ∂p∂x , (S4) ∂v∂t + U ∂v∂x = − ρ ∂p∂y , (S5) ∂w∂t + U ∂w∂x = − ρ ∂p∂z , (S6)and ∂u∂x + ∂v∂y + ∂w∂z = 0 , (S7)where denotes a derivative in the y -direction. The spanwise velocity component is here represented by w forconvenience, which elsewhere refer to the groove width. Eliminating p gives Rayleigh’s equation for v (inviscidOrr-Sommerfeld equation), (cid:20)(cid:18) ∂∂t + U ∂∂x (cid:19) ∇ − U ∂∂x (cid:21) v = 0 , (S8)with p given by 1 ρ (cid:20) ∂ ∂x + ∂ ∂z (cid:21) p = ∂ v∂t∂y + U ∂ v∂x∂y − U ∂v∂x . (S9)We further assume an interface shape of a wave travelling in the streamwise direction, η = Ae ik x ( x − ct ) cos( k z z ) = Ae k x c i t e ik x ( x − c r t ) cos( k z z ) = a ( t ) e ik x ( x − c r t ) cos( k z z ) , (S10)where η is the location of the interface, A is the wave amplitude, k x is the streamwise wavenumber, k z is thespanwise wavenumber and c = c r + ic i . Henceforth, we’ll use the notation k = p k x + k z . The curvature of theinterface gives rise to a pressure difference over the surface of (going from negative to positive y ),∆ p cap = γ (cid:18) ∂ ∂x + ∂ ∂z (cid:19) η = − γk η. (S11)The perturbation must die out at infinity, v → y → ∞ . (S12)In addition to this, fluid parcels on the interface must remain on the interface, (so that the interface remains astreamline) [15], ∂η∂t + U ∂η∂x = v = ⇒ vU − c = ik x η on y = η ≈ . (S13)Returning now to Rayleigh’s equation, eq. (S8), with ( x, z, t )-dependence of v as η ,( U − c ) ∂ v∂y − [ k ( U − c ) + U ] v = 0 . (S14)From eq. (S9), with ( x, z, t )-dependence of p as η , pρ = − i k x k (cid:20) ( U − c ) ∂v∂y − U v (cid:21) . (S15)6 B. Non-dimensionalisation
In order to facilitate the derivations, we introduce the transform, ξ = ky, U − c = U ˜ u ( ξ ) and v = ik x U ˜ v ( ξ ) η, (S16)where U is an arbitrary reference velocity. From the kinematics of the interface, eq. (S13), it follows that˜ v = ˜ u , (S17)where a subscript 0 denote the location of the interface, y ≈
0. The Rayleigh equation becomes˜ v ˜ u − [˜ u + ˜ u ]˜ v = 0 . (S18)The equation for the pressure, eq. (S15), results in p = ρU k x k η (˜ u ˜ v − ˜ u ˜ v ) . (S19)The pressure just above the interface can be written as p +0 = ( α + iβ ) ρU k x k η. (S20)By comparing these two expressions, α + iβ = (˜ u ˜ v − ˜ u ˜ v ) = ˜ u (˜ v − ˜ u ) , (S21)where the last equality comes from eq. (S17). Following ref. [16], ˜ u is assumed to be approximately real and thus β must come from the imaginary part of ˜ u ˜ v , whereas α is the remaining part. C. Derivation of α We here derive an expression for α , which is the real part of eq. (S21), α = < (˜ u ˜ v ) − ˜ u ˜ u . (S22)We use an approximate solution suggested by Miles [16], ˜ v = ˜ u ( ξ ) e − ξ (satisfying the boundary conditions), insteadof solving the full equation. This approximate solution gives α = − ˜ u − ˜ u ˜ u = 1 U (cid:18) − ( U − c ) − U k ( U − c ) (cid:19) . (S23)Now, α can be decomposed into two parts, α = α + α , where α = − c /U and α is the remaining part. Weassume that U can be neglected in comparison to c . The inverse of the wave number has an upper bound,1 k = 1 p k x + k z < k z = λ z π < w π , (S24)since k x > λ z must be smaller than two groove widths. Therefore, U k = 1 kh Re τ Re b U b = 11 kh U b < π wh U b , (S25)for Re τ = 180. For this Reynolds number and groove width w + = 18, the limiting value is 5 . v = 0 at ξ = − kw. (S26)It can also be noticed that kw > k z w ≤ π/ (2 w ) w = π , and e − π = 0 . (cid:28)
1. Hence, if the same assumptionof exponential decrease of ˜ v is made inside the groove, it is reasonable. We therefore assume ˜ v = ˜ ue ξ for y < p − = ρ k x k c ηe ky , y < , (S27)as given by eq. (S9) or eq. (S19). This is the pressure corresponding to α .7 D. Phase speed and minimum phase speed
The difference between the pressure above (eq. S20) and below (eq. S27) the interface must be balanced by thecapillary pressure (eq. S11),∆ p = p − p − ( y = 0) = ( α + iβ ) ρU k x k η − ρ k x k c η = ∆ p cap = − γk η. (S28)Solving for c and including the expression for α ,2 c = γρ k k x + ( α + iβ ) U = ⇒ c = c w (cid:18) α + iβ ) U c w + . . . (cid:19) , (S29)where c w = s γ ρ k k x is the unforced phase speed neglecting any contributions from α and β and the dots denotehigher order terms.It is possible to find the k x for which c w assumes a minimum value by evaluating d c w / d k x = 0. The minimumvalue is, for k z = 2 π/ (2 w ), c w, min = s γρ / k z = s γρ / πw = ⇒ c + w, min = s + w + / π ≈ r π We + w + . (S30) E. Derivation of β In this section, we find an expression for β . Returning to eq. (S18), dividing by ˜ u , multiplying by the complexconjugate of ˜ v , namely ˜ v ∗ , integrating from ξ = ξ to ξ = ∞ , integrating ˜ v ˜ v ∗ by parts, and using the boundaryconditions above (eqs. S12 and S17), Z ∞ ξ {| ˜ v | + [1 + ˜ u / ˜ u ] | ˜ v | } d ξ = [˜ v ∗ ˜ v ] ∞ ξ = − ˜ u ˜ v . (S31)Hence, the imaginary part of eq. (S21) is β = = (˜ u ˜ v ) = −= (cid:18)Z ∞ ξ (˜ u / ˜ u ) | ˜ v | d ξ (cid:19) . (S32)By calculus of residues [16], this expression can be evaluated to β = − π | ˜ v c | ˜ u c ˜ u c = − π (cid:12)(cid:12)(cid:12)(cid:12) v c k x ηU (cid:12)(cid:12)(cid:12)(cid:12) k U c U c , (S33)where the subscript c denotes values at the location where ˜ u = 0, ξ = ξ c . To find an expression for ˜ v c , eq. (S18)can be re-written as (˜ u ˜ v − ˜ u ˜ v ) = ˜ u ˜ v. (S34)Integrating between ξ = ξ c and ξ = ∞ gives, ˜ v c = 1˜ u c Z ∞ ξ c ˜ u ˜ v d ξ. (S35)We now use the same approximate solution of ˜ v as above. This approximation gives β = − π ˜ u c ˜ u c (cid:20)Z ∞ ξ c e − ξ ˜ u d ξ (cid:21) . (S36)The solution to the full equation has been given by Conte & Miles [17]. Assuming a logarithmic profile, U = u τ /κ log( y/z ) = ⇒ ˜ u = log( ξ/ξ c ). We have here set the reference velocity, U = u τ /κ , where κ is the von K´arm´anconstant. With eq. (S36), β = πξ c (cid:20)Z ∞ ξ c e − ξ log( ξ/ξ c ) d ξ (cid:21) . (S37)8 F. Condition for wave growth
To the first order, the imaginary wave speed corresponds to c i = c w β U c w , (S38)giving wave growth when non-zero. The relation for β , eq. (S37), can be computed and is plotted in fig. S7. Thereis a sharp decrease in β at about ξ c ≈
1. This means that in order to have amplitude growth, ξ c = ky c . ⇒ y c . k < wπ , (S39)where we use the upper bound of 1 /k as given by eq. (S24). c ! ! ! - a ) c ! ! ! - ! ! ! ( b ) FIG. S7: Plot of β verses ξ c = ky c , with ( a ) linear and ( b ) logarithmic vertical axis. S4. SPACE-TIME CORRELATIONS
Space-time correlations are shown in fig. S8 for µ i /µ ∞ = 1 and We = 100 and 200. c + = 15 . r + t r + x a ) c + = 10 . r + t r + x b ) -1-0.500.51 FIG. S8: Space-time correlation of interface height in the middle section of a groove, for µ i /µ ∞ = 1 and ( a )We = 100 and ( b ) We = 200. Streamwise displacement is denoted r + x and time difference r + t . The slope of theregion with highest correlation corresponds to the dominating phase speed, marked by a black line. [1] P. Costa, “A fft-based finite-difference solver for massively-parallel direct numerical simulations of turbulent flows,” Comput. Math. with Appl. , vol. 76, no. 8, pp. 1853–1862, 2018.[2] W. P. Breugem and B. J. Boersma, “Direct numerical simulations of turbulent flow over a permeable wall using a directand a continuum approach,”
Phys. fluids , vol. 17, no. 2, p. 025103, 2005. [3] S. Afkhami and M. Bussmann, “Height functions for applying contact angles to 2d vof simulations,” Int. J. Numer.Methods Fluids , vol. 57, no. 4, pp. 453–472, 2008.[4] T. D. Blake, “The physics of moving wetting lines,”
J. Colloid Interface Sci. , vol. 299, no. 1, pp. 1–13, 2006.[5] D. Legendre and M. Maglio, “Comparison between numerical models for the simulation of moving contact lines,”
Comput.Fluids , vol. 113, pp. 2–13, 2015.[6] M. Lee and R. D. Moser, “Direct numerical simulation of turbulent channel flow up to re τ ≈ J. Fluid Mech. ,vol. 774, pp. 395–415, 2015.[7] C. Sch¨onecker, T. Baier, and S. Hardt, “Influence of the enclosed fluid on the flow over a microstructured surface in thecassie state,”
J. Fluid Mech. , vol. 740, pp. 168–195, 2014.[8] M. K. Fu, I. Arenas, S. Leonardi, and M. Hultmark, “Liquid-infused surfaces as a passive method of turbulent dragreduction,”
J. Fluid Mech. , vol. 824, p. 688, 2017.[9] Z. Ge, H. Holmgren, M. Kronbichler, L. Brandt, and G. Kreiss, “Effective slip over partially filled microcavities and itspossible failure,”
Phys. Rev. Fluids , vol. 3, no. 5, p. 054201, 2018.[10] R. G. Cox, “The dynamics of the spreading of liquids on a solid surface. part 1. viscous flow,”
J. Fluid Mech. , vol. 168,pp. 169–194, 1986.[11] J. Jim´enez, “Turbulent flows over rough walls,”
Annu. Rev. Fluid Mech. , vol. 36, pp. 173–196, 2004.[12] P. Orlandi, S. Leonardi, R. Tuzi, and R. A. Antonia, “Direct numerical simulation of turbulent channel flow with wallvelocity disturbances,”
Phys. Fluids , vol. 15, no. 12, pp. 3587–3601, 2003.[13] P. Orlandi and S. Leonardi, “Dns of turbulent channel flows with two-and three-dimensional roughness,”
J. Turbul. ,vol. 7, p. N73, 2006.[14] P. J. Schmid and D. S. Henningson,
Stability and transition in shear flows , vol. 142. Springer Science & Business Media,2012.[15] D. J. Acheson, “Elementary fluid dynamics,” 1991.[16] J. W. Miles, “On the generation of surface waves by shear flows,”
J. Fluid Mech. , vol. 3, no. 2, pp. 185–204, 1957.[17] S. D. Conte and J. W. Miles, “On the numerical integration of the Orr-Sommerfeld equation,”