{'fraction_non_alphanumeric': 0.04176342525399129, 'fraction_numerical': 0.01527576197387518, 'mean_word_length': 4.431809223492314, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical estimate of the viscous damping of capillary-gravity waves: A macroscopic depth-dependent slip-length model

We propose a numerical approach to regularize the contact line singularity appearing in the computation of viscous capillary-gravity waves with moving contact line in cylindrical containers. The linearized Navier-Stokes equations are complemented by a macroscopic Navier-like slip condition on the container side wall, with a depth-varying slip-length bridging a free-slip condition at the contact line to a no-slip condition further away from the meniscus. In accordance with suggestions from the literature, this characteristic penetration depth is chosen as the typical Stokes layer thickness pertaining to the eigenfrequency. Since the latter is unknown, the resulting nonlinear eigenvalue problem is first simplified using the inviscid eigenfrequency to estimate the Stokes layer thickness. The solution is then shown to provide consistent results when compared to the asymptotic approaches of the literature using the inviscid eigenmodes to estimate the different sources of dissipation with the exception of the neglected corner regions. This demonstrates that such numerical regularization scheme is suitable to retrieve physically meaningful results.


I. INTRODUCTION
Sloshing motion represents an archetypal damped oscillator in fluid mechanics. The eigenfrequencies of capillary-gravity waves in closed basins, e.g. in circular cylinder, were traditionally derived in the potential flow limit 1 , while the linear viscous dissipation at the free-surface, at the solid walls and in the bulk for low-viscosity fluids was typically accounted for by a boundary layer approximation 2,3 . These classic theoretical approaches are built on the assumption that the free liquid surface, η, intersects the lateral wall orthogonally, and the contact line can freely slip with constant zero slope, ∂ n η = 0 (free-end edge condition), (1) where ∂ n is the spatial derivative in the direction normal to the lateral wall. However, when the fluid viscosity is formally accounted for, e.g. in a numerical scheme, the boundary conditions on the tangential velocity components at the lateral solid wall require a special treatment, in order to avoid the contact line singularity.
With regard to circular cylinders, at the sidewall, r = R (r is the radial coordinate and R is the container radius), a sliplength model is typically adopted, thus assuming that the fluid speed relative to the solid wall is proportional to the viscous stress [4][5][6] u φ + l s ∂ r u φ − u φ /r = 0, u z + l s ∂ r u z = 0, (2) with u φ and u z the tangential azimuthal and axial velocity components, respectively. Such conditions are indeed required in order to regularize the stress-singularity arising in the viscous solution of a moving contact line over a solid substrate, for which the no-slip condition is incompatible 7,8 . In its original formulation, the slip-length l s appearing in (2) is assumed finite and constant in space along the whole sidewall. The numerical implementation of conditions (1)-(2) with constant l s has been proposed Viola & Gallaire (2018) 6 , where a non-dimensional value of l s = 5 × 10 −4 , corresponding to a dimensional value of l s R = 25 µm for a container radius R = 0.05 m, was used. The numerically estimated damping coefficient associated with the first non-axisymmetric mode was then compared with the analytical prediction by Case & Parkinson (1957) 2 (see Sec. V), who assume a perfect no-slip condition at the wall (the stress-singularity at the contact line is simply ignored in their asymptotic analysis). However, the numerical damping strongly underestimates the one predicted according to Case & Parkinson (1957) 2 , which instead agrees fairly well with their experimental measurements. The origin of such a disagreement can be tentatively ascribed to the fact that a slip-length value of 25 µm is significantly larger than the one observed experimentally 9 , i.e. 0.1 − 10 nm. Lower values of l s could not be employed in Ref. 6 due to the grid resolution needed to guarantee numerical convergence of the results. Nevertheless, there is actually a more subtle aspect intrinsic to the formulation and the numerical implementation of conditions (1)- (2). The dynamics of small-amplitude viscous capillary-gravity waves is essentially governed by the unsteady Stokes equations (see Sec. II), which describe the problem from a macroscopic perspective, where the microscopic (molecular) dynamics at play at the contact line is completely filtered out. The free-end edge contact line condition has indeed only meaning in a macroscopic framework. It follows that the slip-length l s must also assume the significance of a phenomenological macroscopic slip-length, which only serves to regularize the contact line stress-singularity and which is not intrinsically linked to the actual molecular sliplength. When conditions (1) and (2) are simultaneously imposed in a numerical scheme, a paradox concerning the tangential velocity components arises.
A straightforward way to highlight such a paradox consists in combining the free-end edge condition (1), ∂ r η = 0, with the linearized kinematic equation, ∂ t η = u z , differenti-ated with the respect to the radial direction, Relation (3) points out that, by approaching the contact line from the sidewall, one must impose a sidewall stress-free condition on the axial velocity component u z .
A generalization to cylindrical coordinates of the more rigorous reasoning, based on the vorticity argument, outlined by Miles (1990) 10 , suggests indeed that, at the contact line, the only sidewall boundary condition compatible with the freeend edge condition is a stress-free condition for both tangential velocity components, u z and u φ .
In contradistinction with the argument above, the attempt to impose a constant and, at the same time, very small (experimentally inspired) slip-length l s along the entire wall, will translate in nearly no-slip condition everywhere at the sidewall, which is not compatible with the macroscopic free-end contact line condition and which may therefore lead to a numerical unstable solution or to incorrect estimates of the damping coefficients. Obviously, a stress-free condition everywhere would represent a mathematically correct alternative, but its imposition would automatically filter out the corresponding viscous boundary layer, leading to a strong numerical underestimation of the viscous damping 2 .
With the aim to marry the simultaneous need for a sufficiently accurate description of the viscous sidewall boundary layer and a formally correct treatment of the contact line region, we propose in this work the adoption of a phenomenological macroscopic and depth-dependent slip-length model, which constitutes a practical, ad-hoc, boundary condition for numerical calculations.
A similar slip-length model was proposed by Ting & Perlin (1995) 11 (see their Sec. 4.6 and 5). Such a model was developed on the base of their experimental observations for a vertically oscillating plate. In their model, wall-slip occurs within a specific distance from the contact line, whereas the flow obeys the no-slip condition outside this slip region.
The manuscript is organized as follows. The flow configuration and the physical model governing the problem are given in Sec. II. The adopted depth-varying slip-length model, core of the present work, is introduced in Sec. III. Sec. IV is devoted to the description of the numerical method on which results commented in Sec. V are based on. Sensitivity to the free parameters of the slip-length model is explored in Sec. VI. Lastly, final conclusions and comments are outlined in Sec. VII.

II. CONFIGURATION AND GOVERNING EQUATIONS
The natural dynamics of small amplitude viscous capillarygravity waves within a cylindrical vessel is governed by the unsteady Stokes equations, filled to a depth h with a liquid of density ρ and kinematic viscosity ν = µ/ρ. The air-liquid surface tension is γ. The origin of the Cartesian coordinate system is fixed at the center of the free liquid surface at rest, while the bottom is placed at z = −h. The dashed-dotted line is the geometrical axis of symmetry. The interface elevation is denoted by η, whereas g is the gravity acceleration. The macroscopic contact angle is assumed constant and equal to π/2 during the entire motion, i.e. the free surface intersect orthogonally the sidewall. l s (z) denotes the depth-varying slip-length imposed at the lateral wall. which are obtained by linearizing the Navier-Stokes equation around the static configuration, u st = 0, p st = −z, η st = 0, and which are made non-dimensional by using the container's characteristic length R and the velocity √ gR (see Fig. 1). Consequently, the Reynolds number is defined as Re = gR 3 /ν. At the free surface, z = 0, the linearized dynamic and kinematic boundary conditions hold, where ∂ η κ (η) is the first order variation of the curvature associated with the small interface perturbation η (r, φ ,t), n = {0, 0, 1} T is the unit vector normal to the interface, u (r, φ , z,t) = u r , u φ , u z T is the velocity field and p (r, φ , z,t) is the pressure field. The azimuthal coordinate is φ . The Bond number is defined as Bo = ρgR 2 /γ. At the bottom wall at z = −h/R = −H, the no-slip and no-penetration condition apply, so that u = 0. At the contact line, z = 0 and r = 1, we impose the classic free-end edge condition, which prescribes a constant zero slope for the moving interface, At the lateral solid wall, r = 1, it is legit to impose the nopenetration condition, i.e. u r = 0, however, as anticipated in Sec. I, the tangential velocity components require an adhoc treatment, which will be discussed in detail in Sec. III. In the following we proceed by assuming that the boundary conditions for u φ and u z at r = 1 are linear. Hence, the linear homogeneous system can be written in the following matrix compact form (9) We note that the dynamic and kinematic boundary conditions (6a)-(6b) do not explicitly appear in (9), but they are imposed as conditions at the interface 12,13 . More precisely, in the numerical scheme, the kinematic condition governing the state variable η is implemented as an additional equation dynamically coupled with u and p in (9) (see Appendix B of Ref. 12), whereas the three stress components of the dynamic condition are enforced as standard boundary conditions in the corresponding components of the vectorial momentum equation. The solution can be then expanded in terms of normal modes in time and in the azimuthal direction as Substituting the normal form (10a)-(10b) in system (9), we obtain a generalized linear eigenvalue problem, with (q mn (r, z) ,η (r)) = û r ,û φ ,û z ,p,η T the full eigenmode and λ mn = −σ mn + i ω mn the associated eigenvalue. σ mn and ω mn denote, respectively, the corresponding damping coefficient and natural frequency. Here the indices (m, n) represent the number of nodal circles and nodal diameters, respectively, with m also commonly known as azimuthal wavenumber. Owing to the normal mode expansion (10a)-(10b), we note that the operator A m is complex and it depends on the azimuthal wavenumber, as φ -derivatives produce i m terms. Lastly, in order to regularize the problem at the container axis, r = 0, depending on the selected azimuthal wavenumber m, different regularity conditions must be imposed 6,14 :

III. A MODIFIED DEPTH-VARYING SLIP-LENGTH MODEL
In the spirit of Ref. 11, we propose here the following phenomenological macroscopic and depth-dependent slip-length model with l s (z) described by the exponential law where l cl is the slip-length value at the contact line, r = 1 and z = 0, while l δ is its value at a distance δ below the contact line, r = 1 and z = −δ , with δ representing the size of the slip region 11 . In principle, l cl , l δ and δ are all free parameters. However, keeping in mind that one aims at mimicking a stress-free condition in the vicinity of contact line and a noslip condition after a certain distance δ , the natural choice is l cl 1 (∼ 10 2 ÷10 4 ) and l δ 1 (∼ 10 −4 ÷10 −6 ). The range of values proposed in brackets is based on the sensitivity analysis of Sec. 6 and on the discussion outlined in App. A.
The definition of the slip region penetration depth, δ , is more challenging. Miles (1990) 10 postulated that the sliplength l s is a function of the position along the vertical wall and vanishes at a distance comparable to the thickness of the Stokes boundary layer (expressed in non-dimensional terms), which represents the characteristic viscous length scale. Thus, in absence of any further informations, a reasonable choice for δ appears to be δ = δ St . As the oscillation frequency ω mn is not known a priori, the major problem in imposing a variable slip-length depending on δ St and, therefore, on ω mn , is that it turns the generalized eigenvalue problem (11) in a nonlinear eigenvalue problem in λ mn = −σ mn + i ω mn , whose solution would required a fixedpoint iteration-like scheme. Nevertheless, if the Reynolds number is sufficiently large, i.e. Re 1, the viscous correction to ω mn is essentially negligible in norm, i.e. σ mn ω mn , so that one can simply assume ω mn in (15) to be equal to the classical value derived in the potential flow limit and satisfying the well-known dispersion relation for inviscid capillarygravity waves 1 In (16) the wavenumber k mn is given by the nth-root of the first derivative of the mth-order Bessel function of the first kind satisfying J m (k mn ) = 0. This implies that, as a first order approximation, one can define δ = δ St ≈ 2/ (ω inv mn Re), hence making problem (11) fully linear in λ mn .
We note that the value of Reynolds number explored in this work (see Sec. V) ranges from Re ≈ 1 000 to 100 000 ( 1), thus legitimizing the aforementioned approximation.

IV. NUMERICAL METHOD
The numerical scheme used in the eigenvalue calculation is a staggered Chebyshev-Chebyshev collocation method implemented in Matlab. The three velocity components are discretized using a Gauss-Lobatto-Chebyshev (GLC) grid, whereas the pressure is staggered on Gauss-Chebyshev (GC) grid. Accordingly, the momentum equation is collocated at the GLC nodes and the pressure is interpolated from the GC to the GLC grid, while the continuity equation is collocated at the GC nodes and the velocity components are interpolated from the GLC to the GC grid. This results in the classical P N -P N−2 formulation, which automatically suppress spurious pressure modes in the discretized problem. A twodimensional mapping is then used to map the computational space onto the physical space. Lastly, the partial derivatives in the computational space are mapped onto the derivatives in the physical space, which depend on the mapping function. For other details see Refs. 6, 15. The smallest length scale in the flow is represented by the thickness of the oscillating Stokes boundary layer, δ St = 2/ (ω inv mn Re). In order to numerically solve the macroscopic flow at all scales, the computational grid used to discretized the eigenvalue problem (11) must be smaller than such a length scale. The computationally most challenging case analyzed in this work (see Sec. V) corresponds to H = h/R = 3, Re = 99 045 and ω mn = 10 (nondimensional frequency of mode (m, n) = (3, 10)). For this case one has δ St = 0.0016. By using a grid made of 100 and 150 radial and axial GLC nodes, respectively, the size of the first radial grid element at the lateral wall is ∆r 1 = 0.1571 δ St , whereas that of the first axial grid element at the free-surface (same for the bottom wall) is ∆z 1 = 0.2081 δ St . Smaller ratio ∆r 1 /δ St and ∆z 1 /δ St hold in all other cases examined. This ensures that the viscous boundary layers are properly accounted for and that results are sufficiently well converged.

V. COMPARISON WITH THEORETICAL ESTIMATES
We neglect more complex phenomena acting at the contact line, e.g. hysteretic contact angle dynamics due to solid-like friction [16][17][18][19][20][21][22][23][24] , or at the free surface, e.g. surface contamination, etc.. Under these assumptions, capillary-gravity waves, over their motion, are damped by viscous dissipation occurring (i) at the oscillating Stokes layers at the solid bottom and sidewall, (ii) in the fluid bulk and (iii) at the free surface, the latter being typically negligible for ideal surface waves (clean surface). The damping coefficient associated with these three sources of dissipation was evaluated asymptotically in the limit of small kinematic viscosity ν (or large Re 1) by Case & Parkinson (1957) 2 and by Miles (1967) 3 (17), Henderson & Miles (1990) 25 proposed the following viscous correction to the inviscid natural frequencies In this section, we compare our numerical estimates with these theoretical predictions. The free slip-length model parameters are here set to l cl = 10 3 , l δ = 10 −5 and δ = δ St . A sensitivity analysis to variations of these parameters is outlined in Sec. VI. For convenience, the fluid properties are fixed hereinafter, i.e. pure water with ρ = 1000 kg/m 3  In Fig. 2 the numerically computed damping coefficient, σ mn , associated with the first ten modes for m = 0, 1, 2, 3 is compared with the corresponding theoretical prediction according to (17). The non-dimensional fluid depth was kept constant and equal to H = 3 (deep water regime).
The agreement between the two estimates is in general fairly good in the all range of Re, with a maximum deviation of approximately 12 %. It is important to note that such a deviation does not represent an error, but only a difference between theoretical and numerical predictions. Indeed, neither the asymptotic formula (17), which assumes full no-slip along the wall till the contact line, nor the present numerical model, based on the spatially dependent slip-length model (14), are accurate representations of the actual sidewall boundary condition. They only serve as a practical (theoretical or numerical) phenomenological modeling of the actual condition and their resulting prediction should be therefore as closest as possible. The agreement appears striking if compared with a deviation of over 100 % observed in Ref. 6, where a small and constant slip-length was numerically imposed. Fig. 3 shows a similar comparison, but in terms of natural frequencies, with the theoretical estimates computed using (18). In this case, the maximum deviation is about 1 % only, hence confirming that the viscous correction to the natural frequencies when Re 1, e.g. Re 1000, is only mild and, essentially, negligible.
Lastly, with focus on the first non-axisymmetric mode,  Fig. 2, but in terms of natural frequencies, ω mn . The associated Bo = ρgR 2 /γ is reported in the legend of panel (d). Black empty squares: theoretical prediction according to (18). Colored crosses: present numerical estimates.  A similar comparison is shown in Fig. 4, where Fig. 2 of Ref. 2 has been reproduced. For a carefully polished sidewall, the agreement between experiments and theoretical predictions is rather satisfactory, thus suggesting that the analytical approximation (17) is adequate under such controlled conditions. The present numerical predictions seem to agree fairly well with these estimates.

VI. SENSITIVITY TO THE SLIP-LENGTH FREE PARAMETERS
As the slip-length model (14) requires the selection of three user-defined parameter values, i.e. l cl , l δ and δ , in this section we explore the sensitivity of the damping coefficient, σ mn , to variations of these parameters. Without loss of generality, we consider here, as a test case, the first non-axisymmetric eigenmode, (m, n) = (1, 1), for pure water and in a container of radius R = 0.02 and H = h/R = 3, for which Re = 8 859 and Bo = 53.8. Similar considerations qualitatively were found to apply to all modes considered in Sec. V.
A visualization of the slip-length variation along the sidewall, as the size of the penetration region δ is varied, is provided in Fig. 5, where colored filled circles correspond to actual grid points in the numerical implementation.
In Fig. 6(a), the length of the slip region is varied in the range δ ∈ [δ St /2, H]. The value of l δ is fixed to 10 −5 , while three large values of l cl , spanning three orders of magnitude, are considered. A similar analysis is displayed in Fig. 6(b), where the value of l cl = 10 3 is kept fixed and three values of l δ are used.
As anticipated in Sec. V, both panels in Fig. 6 show that, as soon as sufficiently large and small end-values, respectively, of l s (z) are imposed at the contact line and at a distance δ , the damping coefficient calculation is not very sensitive to variation of l cl and l δ . The range of values of l cl and l δ used in this analysis is based on the reasoning described in App. A. The only important parameter appears to be the size of the slip region, δ . In Fig. 7 we show how a variation of δ between two cases with, respectively, δ = H and δ = δ St affects the velocity fields close to the wall and in the vicinity of the contact line region. In particular, one can see that when δ = H, the entire lateral wall behaves as stress-free (see Fig. 7(b)). On the contrary, when δ = δ St , the sidewall is essentially a noslip wall and, indeed, both the tangential velocity components go to zero in the Stokes boundary layer region as r → 1 (see Fig. 7(c)), except for a very small region in the neighborhood of the contact line, i.e. r → 1 and z → 0.
These considerations are consistent with the damping coefficient trend displayed in Fig. 6. When δ = H, the sidewall contribution to the damping coefficient is negligible, so that the total damping is generated in the bulk and in the solid bottom only (the latter is actually negligible in the deep water regime). By gradually reducing δ , the damping increases and eventually, for δ ≈ δ St , it seems to saturate to a nearly constant value which is relatively close to the one predicted by (17), with a deviation of approximately 7% for the present case (see also Sec. V). With regard to Fig. 6, we note that the typical penetration depth of mode (1, 1) is 2π/ (Hk 11 ) = 1.1375 ≈ 1, so that a variation of δ expressed as a fraction of the fluid depth H is approximately equivalent to express the same variation as a fraction of the penetration length of the wave.

VII. CONCLUSION
In this manuscript, a phenomenological macroscopic and depth-dependent slip-length model for practical numerical implementation was presented. This model, inspired by that proposed by Ting & Perlin (1995) 11 , consists of a rapid exponential variation of the slip-length value along the sidewall, hence enabling one to simultaneously regularize the stresssingularity at the contact line and to properly account for the viscous boundary layer along the wall outside the slip contact line region.
A thorough comparison with previous approximated theoretical estimates and a sensitivity analysis to the free model parameters was outlined, showing that the present model is suitable to predict with fair accuracy the damping coefficients and natural frequencies of confined viscous and smallamplitude capillary-gravity waves under controlled conditions, i.e. in absence of free surface contamination and com- plex contact angle dynamics. Although these conditions are only ideal and rarely met in real-life lab-scale experiments, the present slip-length boundary condition represents a functional sidewall condition for numerical implementation, which could allow one to revisit a wide range of existing experimental and numerically-based studies dealing with a moving contact line 6,12,13 . A straightforward generalization of the present model, would be the combination of the slip-length condition with the linear contact angle law proposed by Hocking (1987) 17 . For instance, when these two conditions are combined, the artificial arbitrariness on the value of l cl in Eq. (14) is removed (see Refs. 10 and 11) and rather replaced by a value in consistence with the proportionality constant linking the contact angle deviation to the interface velocity. Furthermore, the adoption of such a slip-length model offers room for further modeling, which could allow one to generalize it to more involved scenarios including, e.g., presence of a static meniscus, non-ideal wall wetting conditions, etc..
Lastly, different depth-dependent slip-length laws could be tested, e.g. a gaussian law, a raised cosine distribution 26 or other similar bell-like functions which assume a maximum value centered in z = 0 (contact line) and quickly decaying. Some of these directions are being pursued and will be reported elsewhere.

SUPPLEMENTARY MATERIAL
See the supplementary material for further details on the derivation of the corrective factor β discussed in App. A.
Results presented in Sec. V and VI have been computed by imposing a large l cl value, i.e. ∼ 10 2 ÷10 4 , and a small l δ value, i.e. ∼ 10 −4 ÷10 −6 , in order to mimic, respectively, a stress-free condition in a small slip region in the neighborhood of the contact line and a no-slip condition outside that region along the wall. These two ranges of values have been estimated as discussed in the following.
By considering a modified version of the Stokes second problem (in dimensional terms), where the no-slip condition is replaced by the following sliplength condition u (z = 0,t) −U cos (ωt) = l s ∂ z u (z = 0,t) , with U and ω the velocity and the frequency, respectively, of the oscillating moving plate, we can compute a corrective factor (due to the slip-length condition), 0 ≤ β ≤ 1, to the classic dissipation induced by the viscous forces on the plate. This coefficient reads, β = 1 + 2ξ 1 + 2ξ (1 + ξ ) , ξ = l s δ St = l s ω/ (2ν) ≥ 0. (A3) See supplementary material for its full derivation. Hence, the theoretical estimate by Case & Parkinson (1957) 2 for the sidewall contribution to the damping coefficient could be generalized to the case of a constant slip-length (l s ) wall condition by accounting for such a correction as σ mod mn, sidewall = β σ th mn, sidewall , When ξ → 0 (β = 1) and ξ → +∞ (β = 0) the two limiting case of stress-free and no-slip sidewall boundary conditions, respectively, are retrieved. Fig. 8 shows the modified damping coefficient as the value of the constant slip-length is varied. For large l s , i.e. 10 2 , the sidewall contribution to the damping is negligible and the wall behaves essentially as a stress-free wall. On the other hand, for small l s , i.e. 10 −4 , the original theoretical prediction by Case & Parkinson (1957) 2 is retrieved, meaning that the sidewall behaves as a fully no-slip wall. The analysis outlined in Sec. VI indeed confirmed that numerical results based on the depth-varying slip-length model (14) are not very sensitive to variations of l cl and l δ values, once the beforehand mentioned limiting conditions are fulfilled.
As a last comment, we note that the approach followed in the computation of β (see supplementary notes) neglects the effect of the container curvature, by analogy with the asymptotic analysis by Case & Parkinson (1957) 2 .  (17) and accounting for the sidewall correction (A4). Colored dotted lines correspond to σ th mn, bulk+bottom (vanishing sidewall contribution), whereas colored dashed lines correspond to σ th mn, total . Parameters: R = 0.02 m, H = 3 and pure water.