The shape of an axisymmetric meniscus in a static liquid pool: effective implementation of the Euler transformation

Nastaran Naghshineh nxncad@rit.edu School of Mathematics and Statistics, Rochester Institute of Technology, Rochester, NY, 14623, USA Department of Mathematics and Sciences, Rochester Institute of Technology-Dubai, Dubai, 341055, UAE W. Cade Reinberger School of Mathematics and Statistics, Rochester Institute of Technology, Rochester, NY, 14623, USA Nathaniel S. Barlow School of Mathematics and Statistics, Rochester Institute of Technology, Rochester, NY, 14623, USA Mohamed A. Samaha Department of Mechanical and Industrial Engineering, Rochester Institute of Technology-Dubai, Dubai, 341055, UAE School of Mathematics and Statistics, Rochester Institute of Technology, Rochester, NY, 14623, USA Steven J. Weinstein School of Mathematics and Statistics, Rochester Institute of Technology, Rochester, NY, 14623, USA Department of Chemical Engineering, Rochester Institute of Technology, Rochester, NY, 14623, USA

Abstract

We examine the classical problem of the height of a static liquid interface that forms on the outside of a solid vertical cylinder in an unbounded stagnant pool exposed to air. Gravitational and surface tension forces compete to affect the interface shape as characterized by the Bond number. Here, we provide a convergent power series solution for interface shapes that rise above or fall below the horizontal pool as a function of contact angle and Bond number. We find that the power series solution expressed in terms of the radial distance from the wall is divergent, and thus rewrite the divergent series as a new power series expressed as powers of an Euler transformed variable; this series is modified to match the large distance asymptotic behavior of the meniscus. The Euler transformation maps non-physical singularities to locations that do not restrict series convergence in the physical domain, while the asymptotic modification increases the rate of convergence of the series overall. We demonstrate that when the divergent series coefficients are used to implement the Euler transformation, finite precision errors are incurred, even for a relatively small number of terms. To avoid such errors, the independent variable in the governing differential equation is changed to that of the Euler transform, and the power series is developed directly without using the divergent series. The resulting power series solution is validated by comparison with a numerical solution of the interface shape and the small Bond number matched asymptotic solution for the height of the interface along the cylinder developed by Lo (1983, J. Fluid Mech., 132, p.65-78). The convergent power series expansion has the ability to exceed the accuracy of the matched asymptotic solution for any Bond number given enough terms, and the recursive nature of the solution makes it straightforward to implement.

pacs:

I Background and Formulation

The shape of a static interface (meniscus) formed outside a vertical cylinder is partially submerged in an infinite horizontal pool has been studied since the early 1900s. This problem was first examined in the literature by [8] who obtained solutions for cylinders having large radii. Since then, the problem has been solved numerically for wire coating and other various applications [18, 12]. [13] and later [14] used the method of matched asymptotic expansions [17] to predict explicit expressions for the height of the meniscus along the cylindrical wall in the limit of small Bond number, $B$ , as defined in (1d) below. More recent extensions of this configuration is found in coating of optical fiber sensors and microfluids [4, 16] in liquid pools of finite extent.

The problem examined here is configured as follows. An infinite horizontal pool of liquid with density, $\rho$ , is subject to gravity, $g$ , and is in contact with air of negligible density. A cylindrical wall of radius, $R$ , is placed vertically in the pool and the liquid intersects the wall with a contact angle, $\theta\in(0,\pi)$ (measured through the liquid). The air–liquid interface has surface tension, $\sigma$ , and its height above (or below) the static pool, located at $z=0$ , is parameterized as $z=H(r)$ where $r$ is the radial distance from the axis of symmetry of the cylinder. For this configuration, the Young-Laplace equation couples with the hydrostatic pressure field to yield the following dimensionless equation and boundary conditions for $\bar{r}\in[1,\infty)$ :

	${d\over d\bar{r}}\left[{\bar{r}{d\bar{H}\over d\bar{r}}\over\left[1+({d\bar{H}\over d\bar{r}})^{2}\right]^{1/2}}\right]-B\bar{r}\bar{H}=0,$	(1a)
	${d\bar{H}\over d\bar{r}}=-\cot\theta~{}\textrm{at}\ \bar{r}=1,$	(1b)
	$\bar{H}\rightarrow 0~{}\textrm{as}\ \bar{r}\rightarrow\infty,$	(1c)
where
	$\bar{r}={r\over R}~{},~{}\bar{H}={H\over R}~{},~{}B={\rho gR^{2}\over\sigma}.$	(1d)

In (1d), the overbars denote dimensionless variables, and $B$ is the Bond number, which provides a ratio of characteristic scales for gravitational stress, $\rho gR$ , and surface tension stress, $\sigma/R$ , that compete to deform the interface.

The nonlinear system (1d), may be solved approximately on a finite domain by shooting, collocation, or other numerical methods. The intention of this work is to provide, for the first time, an analytical solution for the interface shape governed by (1d) over the entire semi-infinite domain, which is expressed here as a power series. Although the standard power series solution to (1d) diverges (as shown herein), we combine two resummation techniques to overcome this difficulty and to construct a single convergent expansion. Recently, [15, NonNewtonian] constructed convergent series solutions for the boundary layer along a moving wall (the Sakiadis Boundary Layer and its non-Newtonian counterpart) by judiciously choosing an independent variable motivated by its asymptotic behavior at large distances from the wall. In this work, we employ a similar approach to construct a convergent series solution that is asymptotically consistent at both ends of the domain. We combine this asymptotic consistency with the well-known Euler transformation [17] to arrive at our convergent power series solution. The relatively simple structure of this problem enables us to examine the efficacy of asymptotically consistent power series resummation methods for solving nonlinear ODEs introduced by [2].

We validate our power series solution by comparing its results with a numerical solution, as well as the $B\rightarrow 0$ matched asymptotic solution for the height of the meniscus at the cylinder by [14]. As a general note here, asymptotic series are often divergent and are thus limited by optimal truncation constraints. Additionally, the effort the effort to determine higher order corrections can be significant [17]. On the other hand, asymptotic methods are particularly efficient at describing limiting regimes with just a few terms and are a staple of fluid mechanics and other disciplines. The advantage of the power series method used here is that no overlap region is required as in matched asymptotics, and, by virtue of being a convergent power series series, our solution may be refined to within machine precision.

As stated above, the Euler transformation is embedded in our power series result. It is known that when the coefficients of the original divergent power series are used to construct the Euler coefficients, finite precision errors inherent to the original series coefficients become magnified [Scraton]—this restricts solution accuracy beyond a certain number of series terms. However, we show in this paper that the computational issue may be circumvented by transforming the radial variable in the system (1d) to the Euler transform variable and obtaining the power series solution to the transformed system directly. This result is important as this deficiency in the Euler transformation has been observed in other problems of mathematical physics (see, for example, Boyd1999 [5]).

The paper is organized as follows. In section II, the power series solution to the system (1d) is obtained and is found to be divergent. The series is transformed via resummation in section III such that it matches the large distance asymptotic behaviour and the influence of convergence-limiting singularities are mapped outside of the transformed domain via an Euler transformation; this results in a convergent representation of the solution as a power series in the Euler transformation variable. Section IV provides a comparison of the convergent power series and the numerical solution. Section V provides an algorithm, which uses the convergent series representation and Newton’s method to predict the height of the interface at the cylinder wall. In section VI, the results are compared with the first and second order matched asymptotic solutions of [14]. Concluding remarks of this study are provided in section VII. In Appendix A, the first and second order matched asymptotic solutions developed by [14] are rewritten for reference in the current notation, and in doing so, we correct two typos in that work. Formulae used to manipulate power series are provided in Appendix B. Appendix C provides a useful solution, valid for all Bond numbers, for configurations where the value of the contact angle, $\theta$ , lies near $\pi/2$ . Appendices D and E provide the details for determining the Euler coefficients by means of a variable substitution in the original governing ODE in system (1d).

II Divergent Power Series Solution

We now consider the solution of (1d) via power series. To begin, we make the transformation

	$\tilde{r}=\bar{r}-1,$	(2a)
such that $\tilde{r}=0$ corresponds to the location where the meniscus meets cylinder wall, and write $\bar{H}(\bar{r})=\bar{H}(\tilde{r}+1)=\bar{h}(\tilde{r})$ in equation system (1d) to obtain
	$(\tilde{r}+1)\bar{h}^{\prime\prime}=B(\tilde{r}+1)\bar{h}\left[(\bar{h}^{\prime})^{2}+1\right]^{3/2}-(\bar{h}^{\prime})^{3}-\bar{h}^{\prime},$	(2b)
	$\bar{h}^{\prime}=-\cot\theta~{}\textrm{at}\ \tilde{r}=0,$	(2c)
	$\bar{h}\rightarrow 0~{}\textrm{as}\ \tilde{r}\rightarrow\infty.$	(2d)

In (2d), note that we have utilized primes to denote derivatives of $\bar{h}$ with respect to $\tilde{r}$ , and $\tilde{r}\in[0,\infty)$ .

The standard power series solution of (2d) is found by first assuming the form

	$\bar{h}(\tilde{r})=\sum_{n=0}^{\infty}a_{n}~{}\tilde{r}^{n},~{}~{}\|\tilde{r}\|<\tilde{r}_{s}(B,\theta),$	(3a)
where $\tilde{r}_{s}(B,\theta)$ is the yet undetermined radius of convergence as function of $B$ and $\theta$ . Equation (3a) is then substituted into (2b) and—after applying JCP Miller’s formula and Cauchy’s product rule (see Appendix B)—terms of like powers are equated to obtain
	$a_{n+2}=\displaystyle{Bf_{n}-c_{n}-(n+1)a_{n+1}-n(n+1)a_{n+1}\over(n+1)(n+2)},$	(3b)
	$f_{n>0}=e_{n-1}+e_{n},~{}f_{0}=e_{0},~{}c_{n}=\displaystyle\sum_{j=0}^{n}(n-j+1)b_{j}a_{n-j+1},$	(3c)
	$e_{n}=\displaystyle\sum_{j=0}^{n}(a_{j})(d_{n-j}),~{}d_{0}=\displaystyle\tilde{b}_{0}^{3/2},~{}d_{n>0}={1\over n\tilde{b}_{0}}\sum_{j=1}^{n}\left({5\over 2}j-n\right)\tilde{b}_{j}d_{n-j},$	(3d)
	$\tilde{b}_{0}=1+b_{0},~{}\tilde{b}_{n>0}=b_{n>0},~{}b_{n}=\displaystyle\sum_{j=0}^{n}(j+1)(n-j+1)a_{j+1}a_{n-j+1},$	(3e)
with
	$a_{0}=\bar{h}(0),~{}\textrm{and}\ a_{1}=-\cot(\theta).$	(3f)

Note that in (3f), the value of $a_{0}$ is not specified.

At this point, we use the value of $a_{0}$ generated by a numerical solution to examine the behavior of the power series solution; in section V to follow, an algorithm to predict $a_{0}$ using the power series itself is provided. The numerical solution used here is a Chebyshev spectral method (the Chebfun package) [7] that implements piecewise polynomial interpolation to solve the nonlinear differential equation as a boundary value problem on a finite domain length $\mathcal{L}$ ; the length $\mathcal{L}$ is chosen so that doubling its size affects the finite domain solution below an infinity norm (taken between the $\mathcal{L}$ and $2\mathcal{L}$ solutions) of $O(10^{-13})$ . In figure 1, the dashed curves show the $N$ -term truncation of series (3f) for $B=0.1$ and $\theta=\pi/4$ . The power series solution (3f) agrees with the numerical solution of (2d) for small $\tilde{r}$ but ultimately diverges at a finite radius of convergence (indicated as a solid vertical line in figure 1); the radius of convergence $\tilde{r}_{s}\approx 0.26$ , which is predicted using the Domb-Sykes plot in figure 2 (the limit as $n\to\infty$ provides the radius of convergence); the Domb-Sykes plot itself is the numerical implementation of the ratio test [17]. Similarly divergent results are obtained for other values of $\theta$ and $B$ .

Refer to caption — Figure 1: The solution to (2d) is shown for Bond number $B=0.1$ and a contact angle of $\theta=\pi/4$ . The $N$ -term truncations of the divergent series solution (3f) (dashed curves) and the convergent series solution (8j) (solid curve) are compared against the numerical solution with $\mathcal{L}=60$ ( $\bullet$ ’s). The solid vertical line shows the radius of convergence, $\tilde{r}_{s}\approx 0.26$ , computed from the Domb-Sykes plot.

III Convergent Resummation

III.1 Euler Transformation

Our goal is to recast the divergent series (3f) as another power series that converges uniformly over the entire semi-infinite domain. This requires a variable transformation that maps the convergence-limiting singularity to a new location to increase the radius of convergence. To do so, we note that the coefficients of the divergent series expansion (3f) alternate in sign. This indicates (by a corollary to Pringsheim’s theorem) that the singularity responsible for divergence lies along the negative real axis [Hinch]. As such, it is likely that application of the Euler transformation will induce convergence.

The Euler transformation is defined as

	$\bar{h}(\tilde{r})=g(\eta(\tilde{r})),$	(4a)
with
	$\displaystyle\eta={\tilde{r}\over\tilde{r}+S},$	(4b)

where $S>0$ is a parameter corresponding to the approximation of the radius of convergence of the diverging series (3f). Given a power series in $\tilde{r}$ with a negative real convergence limiting singularity located at $\tilde{r}=-\tilde{r}_{s}$ , a power series in $\eta$ ,

g=\sum_{n=0}^{\infty}\hat{A}_{n}\eta^{n},

(5)

converges pointwise in $\eta$ for any positive real $S$ with $S<2\tilde{r}_{s}$ (see Appendix B.6). The Euler transformation (4b) maps the domain of the new expansion variable to be between $\eta=0$ and $\eta=1$ , and maps the singularity such that the transformed disc of convergence contains the entire physical domain, thus creating a convergent series. As shown in appendix B.5, the power series coefficients of the transformed function $g$ about $\eta=0$ can be written explicitly in terms of the original Taylor coefficients $a_{n}$ in (3f) about $\tilde{r}=0$ as

\hat{A}_{n>0}=\sum_{m=1}^{n}\binom{n-1}{m-1}a_{m}S^{m},~{}\hat{A}_{0}=a_{0}.

(6)

The above relation is attractive in that it may be applied to any divergent alternating series with a finite radius of convergence, whether it originates as a solution to an ODE or not. Thus, for purposes of generality, the Euler summation is often presented using an identity equivalent to (6) [Hardy]. However, the implementation of Euler summation via (6) is computationally unstable to finite precision errors [Scraton]. This is shown in figure 3, where coefficients computed using (6) deviate from their true value beyond some finite $n$ value (here $n\approx 50$ ). This error originates through the recurrence relation in (3f) for the divergent series coefficients, since it incorporates the sums and differences of expressions involving exponentially growing series coefficients (see Figure 3). Since these coefficients are used in (6) to compute $\hat{A}_{n}$ , the Euler summation coefficients incur this error.

To avoid finite precision error in computing the Euler transformation (4b), then, one must avoid using the original diverging series coefficients. This is accomplished by applying the variable transformation (4b) directly to the original ODE (2b) to obtain an ODE for $g(\eta)$ . The Taylor coefficients $\hat{A}_{n}$ are then found as the coefficients to the power series solution in $\eta$ to the transformed ODE. Upon applying the Euler transformation (4b), the ODE (2d) becomes

	$g^{\prime\prime}={2\over 1-\eta}g^{\prime}+{((1-\eta)^{4}g^{\prime 2}+S^{2})^{3/2}Bg\over(1-\eta)^{4}}-{(1-\eta)^{3}g^{\prime 3}\over S(\eta S+1-\eta)}-{Sg^{\prime}\over(1-\eta)(\eta S+1-\eta)},$		(7a)
	$g^{\prime}=-S\cot\theta~{}\textrm{at}\ \eta=0,$		(7b)
	$g\to 0~{}\textrm{as}\ \eta\to 1.$		(7c)

In (7c), the primes denote derivatives of $g$ with respect to $\eta$ , and $\eta\in[0,1)$ . We next assume a solution to (7a) of the form

	$g(\eta)=\sum_{n=0}^{\infty}\hat{A}_{n}\eta^{n},~{}~{}\|\eta\|<\eta_{s}(B,\theta).$	(8a)
Using JCP Miller’s formula (B.1b) and Cauchy’s product formula (B.2), a recurrence can be obtained (see Appendix D) for the Taylor coefficients $\hat{A}_{n}$ as
	$\hat{A}_{n+2}={S^{2}\over(n+1)(n+2)}\sum_{k=0}^{n}\binom{k+3}{3}w_{n-k},$	(8b)
where
	$w_{n}=f_{n}+\sum_{k=0}^{n}(q_{k}-t_{k}-u_{k})r_{n-k}~{}~{},~{}~{}f_{n}={2\over S^{2}}\sum_{k=0}^{n}(-1)^{k}\binom{3}{k}(n-k+1)\hat{A}_{n-k+1},$	(8c)
	$q_{n}=B\sum_{k=0}^{n}\ell_{k}p_{n-k}~{}~{},~{}~{}\bar{\ell}_{n}=S+(1-S)\delta_{n,0}~{}~{},~{}~{}\ell_{n}=\sum_{k=0}^{n}\bar{\ell}_{k}\hat{A}_{n-k}$	(8d)
	$u_{n}={1\over S}\sum_{k=0}^{n}(-1)^{k}\binom{2}{k}(n-k+1)\hat{A}_{n-k+1}$	(8e)
	$d_{n}=\sum_{k=0}^{n}u_{k}u_{n-k}~{}~{},~{}~{}t_{n}=\sum_{k=0}^{n}d_{k}u_{n-k}~{}~{},~{}~{}\bar{d}_{n}=d_{n}+\delta_{n,0}$	(8f)
	$r_{n}=(1-S)^{n}-(1-S)^{n-1}(1-\delta_{n,0}),$	(8g)
	$p_{n>0}={1\over n\bar{d}_{0}}\sum_{k=1}^{n}\left({5\over 2}k-n\right)\bar{d}_{k}p_{n-k}~{}~{},~{}~{}p_{0}=\bar{d}_{0}^{3/2},$	(8h)
with
	$\hat{A}_{0}=\bar{h}(\tilde{r}=0),~{}\textrm{and}~{}~{}\hat{A}_{1}=\cot(\theta),$	(8i)
where the Kronecker notation is used such that $\delta_{n,0}$ is 0 when $n\neq 0$ and 1 when $n=0$ . According to equation (8a) and the Euler transformation (4b), the interface shape can be written in terms of physical domain coordinates as
	$\bar{h}(\tilde{r})=\sum_{n=0}^{\infty}\hat{A}_{n}\left({\tilde{r}\over\tilde{r}+S}\right)^{n}.$	(8j)

The result (8j) is the desired uniformly convergent representation of the solution via the Euler transformation. Returning to figure 3, a comparison is provided between the coefficients computed by transforming divergent coefficients according to (6) and those directly from the transformed differential equation according to (8j). We observe that for $n<50$ , coefficients obtained are identical. However, for larger $n$ values, we see that the coefficients calculated from (6) begin to grow due to aforementioned finite precision error, and deviate from those obtained by (8j). We note here that the solution (8j) converges slowly, and that is evidenced by the slow decline in the magnitude of $\hat{A}_{n}$ as $n$ increases in figure 3.

III.2 Prefactor Inclusion

Convergence of the power series solution (8j) may be improved by incorporating the asymptotic behavior of the solution for $\bar{h}(\tilde{r})$ as $\tilde{r}\to\infty$ . We note that the power series solution (3f) implicitly matches the correct asymptotic behavior as $\tilde{r}\rightarrow 0$ , and thus efficiently represents the solution near the cylinder wall. To guide the form of a convergent resummation of (3f), the asymptotic behavior of $\bar{h}$ as $\tilde{r}\rightarrow\infty$ is first determined. The height constraint (2d) necessarily implies that $\bar{h}^{\prime}\to 0$ as $\tilde{r}\to\infty$ . Thus, using the method of dominant balance on (2b) we neglect terms that are quadratic or higher in $\bar{h}^{\prime}$ and obtain

(\tilde{r}+1)^{2}\bar{h}^{\prime\prime}+(\tilde{r}+1)\bar{h}^{\prime}-B(\tilde{r}+1)^{2}\bar{h}=0~{}\text{as}\ \tilde{r}\to\infty.

(9)

The solution of (9) is expressed in terms of the modified Bessel function (of the second kind) of zeroth order as

\bar{h}(\tilde{r})\sim D~{}K_{0}\left(\sqrt{B}(\tilde{r}+1)\right)~{}\text{as}\ \tilde{r}\to\infty,

(10)

where $D$ is an undetermined constant that can only be found through consideration of finite $\tilde{r}$ behavior. For $\theta$ near $\pi/2$ , note that the asymptotic result (10) is an excellent approximation for $\bar{h}$ for all $\tilde{r}$ , and the boundary condition (2c) may be applied to determine $D$ (see Appendix C). For smaller angles, the power series (8j) may be improved by combining the modified Bessel function (10) with the Euler transformation (4b) as

	$\bar{h}(\tilde{r})=K_{0}\left(\sqrt{B}(\tilde{r}+1)\right)\sum_{n=0}^{\infty}\bar{A}_{n}~{}\eta^{n}.$		(11a)
In (11a), $\eta$ is the Euler transformation defined in (4b), and the estimate for the radius of convergence, $S$ , is determined in what follows. Note that as $\tilde{r}\rightarrow\infty$ , the series in (11a) approaches a constant, preserving the asymptotic Bessel function behavior (10) in this limit.

To solve for the Euler coefficients, $\bar{A}_{n}$ , in (11a), we use Cauchy’s product formula (see Appendix B.2), and the formula for the Taylor coefficients of $K_{0}(\sqrt{B}(\tilde{r}+1))^{-1}$ in an expansion about $\eta=0$ , relating $r$ and $\eta$ through (4b). This is done using the defining ODE for the modified Bessel function, and applying the same techniques used to compute the coefficients $\hat{A}_{n}$ in Appendix D; details are provided in Appendix E to obtain the coefficients

\bar{A}_{n}=\sum_{k=0}^{n}\zeta_{k}\hat{A}_{n-k},

(11b)

where

\zeta_{n}={-1\over\bar{\zeta}_{0}}\sum_{k=1}^{n}\bar{\zeta}_{k}\zeta_{n-k}~{}~{},~{}~{}\zeta_{0}=\bar{\zeta}_{0}^{-1}~{}~{},~{}~{}\bar{\zeta}_{0}=K_{0}(\sqrt{B})~{}~{},~{}~{}\bar{\zeta}_{1}=-S\sqrt{B}K_{1}(\sqrt{B}),

(11c)

\bar{\zeta}_{n+2}={1\over(n+1)(n+2)}\left[\sum_{k=0}^{n}\left\{(k+1)\bar{\zeta}_{k+1}(2-\rho_{n-k})+BS^{2}\binom{k+3}{3}\bar{\zeta}_{n-k}\right\}\right]

(11d)

\rho_{n}=1-(1-S)^{n+1}.

(11e)

In writing (11a), we use the fact that the Bessel function prefactor does not introduce new singularities that impact convergence, since the actual function that is Euler-transformed and then subsequently expanded in a power series is $[\bar{h}(\tilde{r})][K_{0}(\sqrt{B}(\tilde{r}+1))^{-1}]$ . It is well known [Parnes] that the zeros of $K_{0}$ occur in the left-half plane, which means any value of $S$ that successfully Euler sums $\bar{h}$ must also do so for the function $[\bar{h}(\tilde{r})][K_{0}(\sqrt{B}(\tilde{r}+1))^{-1}]$ .

To utilize (8j), the estimate for the radius of convergence, $S$ , is needed (recall that Appendix B.6 indicates that S does not need to be exact). Note that the Domb-Sykes plot could be used to establish its value exactly for given parameter values $\theta$ and $B$ . However, to simplify implementation, a systematic estimate may be obtained by considering the asymptotic limits of small and large Bond numbers. When the Bond number is small ( $B\rightarrow 0$ ), both sides of (2b) can be integrated once, and is thus reduced to a $1^{\mathrm{st}}$ -order ordinary differential equation as

\displaystyle\bar{h}^{\prime}={-\cos\theta\over\sqrt{(\tilde{r}+1)^{2}-\cos^{2}\theta}},

which has a singularity at the location where the denominator is zero and the slope is infinite. The singularity lies on the negative $\tilde{r}$ axis, whose precise location depends on $\theta$ . The location of infinite slope imposes a radius of convergence on the physical domain from the closest root to the $\tilde{r}=0$ location, and the radius of convergence imparted is given by

\tilde{r}_{s}|_{B\rightarrow 0}=|1-|\cos\theta||,~{}~{}B\ll 1.

(12)

To find $\tilde{r}_{s}$ of (3f) as $B\rightarrow\infty$ , one can interpret the problem dimensionally as that of a cylinder with a large radius in a fluid with fixed physical properties. In such a limit, the cylinder surface may be viewed as flat over a large portion of the domain. The meniscus of a flat wall satisfies an autonomous ODE and has an analytical solution [3]. The convergence limiting singularity corresponds to the location where the slope of the interface solution is infinite even if outside the physical domain. That singularity is imparted because all flat wall interface solution lie along the same master curve that is translated to meet conditions at the wall. This location is given by [15] and is rewritten below in terms of current notation and non-dimensionalized variables as

\tilde{r}_{s}|_{B\rightarrow\infty}={\left|\cosh^{-1}\sqrt{2}-\cosh^{-1}\sqrt{{2\over 1-\sin\theta}}+\sqrt{2+2\sin\theta}-\sqrt{2}\right|\over\sqrt{B}},~{}~{}B\gg 1.

(13)

We have surveyed the parameter ranges $\theta\in\left(0,\pi\right)$ and $B\in[0.001,100]$ , and find that the ratio of $\tilde{r}_{s}$ found numerically (using a numerical $a_{0}$ to initiate the recursion in (3f)) to the estimated $\tilde{r}_{s}$ found using (12) and (13) always lies between 0 and 2 when chosen judiciously as follows—this satisfies the constraints for a convergent Euler summation (see Appendix B.6). We have found that the maximum error in the radius of convergence in using either (12) or (13) is obtained when the radii predictions are equal. Setting (12) equal to (13) leads to a master curve to solve for a critical value of $B$ , denoted by $B_{crit}$ , as a function of contact angle $\theta$ as

\displaystyle B_{crit}=\left[{\cosh^{-1}\sqrt{2}-\cosh^{-1}\sqrt{{2\over 1-\sin\theta}}+\sqrt{2+2\sin\theta}-\sqrt{2}\over 1-|\cos\theta|}\right]^{2},

(14)

which is plotted in figure 4. In this figure, errors in the radius of convergence (compared with the numerical solution) become smaller as one moves away from the plotted curve. We can use this fact, then, to determine the estimate for the radius of convergence, $S$ , in the Euler transformation (4b) which is used in the convergent series solutions (8j) and (11e). Consistent with the derivation, (13) is used to determine $S$ for $B>B_{crit}$ , and (12) is used to predict $S$ for $B<B_{crit}$ . For any combination of $B$ and $\theta$ values that fall along the solid line in figure 4, $\tilde{r}_{s}$ (i.e., $S$ ) obtained from (12) and (13) are identical, but (12) is preferred for simplicity. In the results and figures that follow, we use (12) and (13) along with figure 4 to obtain $S$ analytically.

IV Convergent Resummation Results

In this section, we examine the performance of the series (11e) by direct comparison with numerical results. We begin by noting that the system (1d) has a reflection property that the solutions for $\theta<\pi/2$ are identical in shape to those for that $\theta>\pi/2$ provided $\mid\theta-\pi/2\mid$ are the same. The only difference is that the solutions are mirror images, so solutions rising above the horizontal pool as $r\to\infty$ for $\theta<\pi/2$ fall below that same horizontal height. Thus, it suffices to only consider solutions for $\theta<\pi/2$ in what follows. Figure 5a shows the absolute error (the absolute difference) between $N$ -term truncations of the convergent series solution (11e) and the numerical solution for $B=0.1$ and $\theta=\pi/4$ . The value of $\bar{h}(\tilde{r}=0)$ from the numerical solution is used as the value of $a_{0}$ in constructing figure 5a (as well as figures 6a and 7a in this section). Thus, the lowest absolute error (the absolute value of the difference between the convergent series solution (11e) and the numerical solution) in figure 5a is observed at the wall. Note that the indicated errors are continually reduced for all $\tilde{r}$ as $N$ is increased beyond those values shown in the figures, but the trends indicated are maintained. Thus, desired accuracy can be achieved with sufficiently large $N$ .

Figure 6a shows the maximum absolute error (maximum of the absolute value of the difference over $\tilde{r}\in[0,\mathcal{L}]$ ) between the power series representation (11e) and the numerical solution for $B=0.1$ and various $\theta$ values, while figure 7a shows the same for $\theta=\pi/4$ and various $B$ values. We observe that the error between the numerical and series solution decreases as contact angles approach $\pi/2$ for the same number of terms. Similarly, for the same number of terms the power series solution (11e) becomes more accurate as the Bond number increases. Again, increased accuracy may be obtained with increasing $N$ beyond those shown in the figures. Moreover, we note that the series takes more terms to converge as $\theta\to 0$ , and less terms as $\theta\to\pi/2$ , for all Bond numbers.

V Prediction of the Meniscus Height at the Cylinder Wall

Up to this point in the paper, all discussed power series results have utilized the numerically predicted value of the height at the wall as an input. We now provide the methodology to determine the height of the wall, $a_{0}=\bar{h}(\tilde{r}=0)$ using the power series itself, which precludes the need for numerical inputs altogether. The recurrence for the coefficients provides data about the structure of the solution as a function of the chosen value of $a_{0}$ . While there are many conditions we impose on the coefficients of (11e) to yield an equation for $a_{0}$ , the simplest is to set the last $\bar{A}_{n}$ coefficient to zero. This approach has been used previously in [2]. The rationale for this choice is based on the fact that, in a convergent series, the last coefficient, $\bar{A}_{N}$ must approach zero as $N\to\infty$ ; thus the assumption that $\bar{A}_{N}=0$ is self consistent in this limit. From a purely mathematical perspective, setting $\bar{A}_{N}=0$ provides one equation in one unknown $a_{0}$ that enables solution. We can compute the $n$ -th coefficient $\bar{A}_{n}(a_{0})$ as a function of the chosen value of $a_{0}$ . Since the underlying recurrence is analytic, we can differentiate it to compute derivatives with respect to $a_{0}$ as well. As such, we can use Newton’s iteration to obtain the desired result by choosing $a_{0}$ to enforce the condition $\bar{A}_{N}=0$ .

The Newton’s iteration takes the form

	$a_{0,j+1}=a_{0,j}-{\bar{A}_{N}(a_{0,j})\over\bar{A}^{\prime}_{N}(a_{0,j})},$	(15a)
where $\bar{A}^{\prime}_{N}(a_{0,j})$ is the derivative of $\bar{A}_{N}$ with respect to $a_{0}$ . To utilize (15a), we differentiate the coefficients in (11e) with respect to $a_{0}$ , and obtain
	$\bar{A}^{\prime}_{n}=\sum_{k=0}^{n}\zeta_{k}\hat{A}^{\prime}_{n-k},$	(15b)
	$\hat{A}^{\prime}_{n+2}={S^{2}\over(n+1)(n+2)}\sum_{k=0}^{n}\binom{k+3}{3}w^{\prime}_{n-k},$	(15c)
	$w^{\prime}_{n}=f^{\prime}_{n}+\sum_{k=0}^{n}r_{n-k}\left(q^{\prime}_{k}-t^{\prime}_{k}-u^{\prime}_{k}\right)~{}~{},~{}~{}f^{\prime}_{n}={2\over S^{2}}\sum_{k=0}^{n}(-1)^{k}\binom{3}{k}(n-k+1)\hat{A}^{\prime}_{n-k+1},$	(15d)
	$q^{\prime}_{n}=B\sum_{k=0}^{n}\left\{\ell^{\prime}_{k}p_{n-k}+\ell_{k}p^{\prime}_{n-k}\right\}~{}~{},~{}~{}\ell^{\prime}_{n}=\sum_{k=0}^{n}\bar{\ell}_{k}\hat{A}^{\prime}_{n-k},$	(15e)
	$u^{\prime}_{n}={1\over S}\sum_{k=0}^{n}(-1)^{k}\binom{2}{k}(n-k+1)\hat{A}^{\prime}_{n-k+1},$	(15f)
	$d^{\prime}_{n}=\sum_{k=0}^{n}\left\{u^{\prime}_{k}u_{n-k}+u_{k}u^{\prime}_{n-k}\right\}~{}~{},~{}~{}t^{\prime}_{n}=\sum_{k=0}^{n}\left\{d^{\prime}_{k}u_{n-k}+d_{k}u^{\prime}_{n-k}\right\}~{}~{},~{}~{}\bar{d}^{\prime}_{n}=d^{\prime}_{n},$	(15g)
	$p^{\prime}_{n}={1\over n\bar{d}_{0}}\sum_{k=1}^{n}\left\{\left({5\over 2}k-n\right)\left(\bar{d}^{\prime}_{k}p_{n-k}+\bar{d}_{k}p^{\prime}_{n-k}\right)\right\}-{1\over n\bar{d}_{0}^{2}}\bar{d}^{\prime}_{0}\sum_{k=1}^{n}\left({5\over 2}k-n\right)\bar{d}_{k}p_{n-k},$	(15h)
	$p^{\prime}_{0}={3\over 2}\bar{d}_{0}^{1/2}\bar{d}^{\prime}_{0},$	(15i)
with
	$\hat{A}^{\prime}_{0}=1,~{}\textrm{and}~{}\hat{A}^{\prime}_{1}=0.$	(15j)

Using an arbitrary initial guess, (15j) is used to predict $a_{0}$ , which approaches the correct value as $N\to\infty$ , as evidenced by the error plots in figures 5b, 6b, and 7b.

Apart from comparing directly with numerical results, our predictions of the height of the interface at the wall ( $\bar{h}=a_{0}$ ) and resulting interface shape agree well with an analytical expression derived for small interfacial slopes (see Appendix C) over the whole radial domain. This occurs when $\theta$ is near $\pi/2$ . Under such circumstances, the constant $D$ in (10) can be determined based on constraints at the wall. According to (11e), this means that $\bar{A}_{0}=D$ for small slopes. If one utilizes the algorithm (15j) to determine the constant $A_{0}$ by setting $A_{1}=0$ , we obtain the same value of the constant $D$ obtained the small slope approximation. This provides additional insight into how the series (11e) with the modified Bessel prefactor accelerates convergence. For small slopes at the wall, one term in the expansion provides an excellent prediction of the interface shape, while higher order terms in $\eta$ contribute only incremental improvement. Far from the wall, the same comment applies—as the small slope approximation (10) is always valid, the Bessel function prefactor in (11e) assures that only few terms in the expansion in $\eta$ (or equivalently $\tilde{r}$ ) are needed to capture the solution behavior.

VI Performance of Power Series and Matched Asymptotic solutions

In this section, we compare predictions of the height of the interface at the wall, $a_{0}$ from the power series solution (11e) using the algorithm of section V, with those of the matched asymptotic solution by [14] in the limit $B\rightarrow 0$ . Figure 8 compares the solutions generated by the $N$ -term truncation of the convergent power series solution (11e), the first order matched asymptotic solution of [14], and the second order matched asymptotic solution of [14] (see (A.1b) and (A.2) in Appendix A, respectively), with the numerical solution for $B=$ 100, 10, 1, 0.1, 0.01, and 0.001 and $\theta=\pi/4$ , at the wall of the cylinder ( $\tilde{r}=0$ ).

As shown in figure 8, the accuracy of the asymptotic solution improves as $B$ approaches zero, as expected from asymptotic theory. That said, the series solution (11e) is effective over a broad range of Bond numbers, and in fact performs with high accuracy even in the range where asymptotic analysis is justified; that is, provided a sufficient number of terms in the series are included.

In the case of small Bond numbers where Lo’s solution has sufficient accuracy, one could choose to use the value of $a_{0}$ , the height of the liquid meniscus at the wall, from the method of matched asymptotic solutions [14], shown in (A.1b) and (A.2), and use the convergent series solution (11e) to generate the full shape of the interface. The accuracy of the matched asymptotic solution might be improved by considering a third order match, but as is evident from the work of [14], the generation of higher order matches is arduous. Additionally, asymptotic series are subject to optimal truncation errors, so there is a limit to the improvement that may be achieved by adding such corrections. By contrast, the power series approach yields a relatively simple recursion that can be used for solutions of desired accuracy by simply including additional terms in the series.

VII Conclusion

In this work, the asymptotically motivated Euler transformation was utilized to transform a divergent power series solution to a convergent one, where the influence of singularities were mapped outside the physical domain. Furthermore, we have obtained—for the first time—an exact solution for the shape of the meniscus formed outside a partially submerged vertical cylinder in an infinite horizontal static pool. The efficacy of the power series solution has been systematically tested against numerical and previously-derived matched asymptotic solutions. Moreover, the complexity of the matched asymptotic solution precludes additional higher-order corrections from being generated practically. The power series solution, generated recursively, is attractive for its ease of use for all values of $B$ , $\theta$ , and over the entire physical domain.

Appendix A Matched Asymptotic Solution

In this section, we re-visit the first and the second order matched asymptotic solutions implemented by [14] and convert their notation to ours. In doing so, we correct two typos in Eqns 3.18 and 3.21 of Lo’s paper [14], shown respectively with underlines in equations (A.2e) and (A.2f) below. Here we focus on Lo’s predictions of the height of the meniscus along the cylinder to benchmark against the power series solution. Lo uses the superscripts (1) and (2) to represent the first and second order matched asymptotic approximations, respectively, and denotes the well-known Euler’s constant as $\gamma=0.5772$ (rounded to four decimal places). Lo defines the parameter $\epsilon={R\over l_{c}}$ , where $R$ is the cylinder radius and $l_{c}=\sqrt{{\sigma\over\rho g}}$ is the capillary length. In the notation used in our paper, this means that the Bond number and $\epsilon$ are related as $B=\epsilon^{2}$ .

The lowest (first) order matched asymptotic solution for the height of the fluid along the wall, used in figure 8, is

	$\bar{h}^{(1)}(\tilde{r}=0)=\sin\phi\left\{\ln{4\over\epsilon(1+\cos\phi)}-\gamma\right\},$		(A.1a)
	$\phi={\pi\over 2}-\theta.$		(A.1b)

The second order matched solution (next order match) for the height of the interface along the cylindrical wall indicated in figure 8 is

	$\bar{h}^{(2)}(\tilde{r}=0)=\underbrace{z_{1}\ln\epsilon+z_{2}}_{\bar{h}^{(1)}~{}\textrm{in}\ (\ref{eq:FirstH0})}+z_{3}\epsilon^{2}\ln^{2}\epsilon+z_{4}\epsilon^{2}\ln\epsilon+z_{5}\epsilon^{2}.$		(A.2a)
The constants used in (A.2a) are defined as
	$z_{1}=-\sin\phi~{}~{},~{}~{}z_{2}=\sin\phi\left[\ln 4-\ln(1+\cos\phi)-\gamma\right]~{}~{},~{}~{}z_{3}={1\over 2}\sin\phi(1-\sin^{2}\phi),$		(A.2b)
	$z_{4}=\sin\phi\left\{\cos^{2}\phi\left[\gamma+{1\over 4}+\ln({1\over 4}(1+\cos\phi))+\right]+{1\over 4}-\cos\phi\right\},$		(A.2c)

	$\displaystyle z_{5}=$	$\displaystyle~{}{1\over 4}\sin\phi\cos\phi(\ln 4-\gamma)-{D\over\cos\phi}+{1\over 4}\sin\phi$	(A.2d)
		$\displaystyle+D(1-\ln 4{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\underline{+}}\gamma)+{1\over 4}\sin^{3}\phi\left({1\over 2}-\gamma^{3}-\ln^{2}2+\gamma\ln 4-\ln 2\right)$
		$\displaystyle+\left[\ln(1+\cos\phi)\right]\left[{1\over 4}\sin^{3}\phi\left(\ln 4-\gamma+{1\over\cos\phi}\right)+D-{1\over 4}\sin\phi\cos\phi\right]$
		$\displaystyle-{1\over 4}(\sin^{3}\phi)\ln^{2}(1+\cos\phi),$	(A.2e)
where
	$D={1\over 4}\mathcal{C}(2-\mathcal{C}^{2})\ln[1+\sqrt{1-\mathcal{C}^{2}}]-(1-\mathcal{C}^{2}){1\over 2}\mathcal{C}(\ln 4-\gamma)-{1\over 4}\mathcal{C}\sqrt{1-{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\underline{\mathcal{C}^{2}}}},$		(A.2f)
and
	$\mathcal{C}=\sin\phi.$		(A.2g)

Appendix B Useful formulae for manipulating series

B.1 Raising a Series to a Power

The following relation is JCP Miller’s formula for raising a series to a power [11]:

	$\left(\sum_{n=0}^{\infty}a_{n}x^{n}\right)^{\gamma}=\sum_{n=0}^{\infty}b_{n}x^{n},$	(B.1a)
where
	$b_{n>0}={1\over n~{}a_{0}}\sum_{j=1}^{n}(j\gamma-n+j)a_{j}b_{n-j},~{}b_{0}=\left(a_{0}\right)^{\gamma}.$	(B.1b)

B.2 Product of Two Series

The following relation is the well-known Cauchy product of two series [6]:

\sum_{n=0}^{\infty}a_{n}x^{n}\sum_{n=0}^{\infty}b_{n}x^{n}=\sum_{n=0}^{\infty}\left(\sum_{j=0}^{n}a_{j}b_{n-j}\right)x^{n},~{}a_{0}\neq 0.

(B.2)

B.3 Generalized Product Rule

The generalized product (Leibniz’s) rule applied to two functions is given as [10]

\displaystyle{d^{n}\left(uv\right)\over dy^{n}}=\sum_{k=0}^{n}{{n}\choose{k}}~{}{d^{n-k}u\over dy^{n-k}}~{}{d^{k}v\over dy^{k}}

(B.3)

and is important in determining the coefficients of an Eulerized series.

B.4 Chain Rule

The chain rule applied to a composite function $F(y(\eta))=f(\eta)$ is given, for first and second derivatives, as

	${df\over d\eta}={dF\over dy}{dy\over d\eta},$		(B.4a)
	${d^{2}f\over d\eta^{2}}={d^{2}F\over dy^{2}}\left({dy\over d\eta}\right)^{2}+{dF\over dy}{d^{2}y\over d\eta^{2}}.$		(B.4b)

The inversion of (B.4b) is given by

	${dF\over dy}={df\over d\eta}\left({dy\over d\eta}\right)^{-1},$		(B.5a)
	${d^{2}F\over dy^{2}}={d^{2}f\over d\eta^{2}}\left({dy\over d\eta}\right)^{-2}-~{}{df\over d\eta}{d^{2}y\over d\eta^{2}}\left({dy\over d\eta}\right)^{-3},$		(B.5b)

which will be useful in manipulating the Bessel series in terms of the Euler transformed variable in appendix E.

B.5 Euler Transformation

Given a power series

f=\sum_{n=0}^{\infty}a_{n}x^{n},

(B.6)

the Euler Sum is defined as follows

	$\displaystyle f=\sum_{n=0}^{\infty}b_{n}\left({x\over x+S}\right)^{n},$	(B.7a)
where
	$b_{n>0}=\sum_{m=1}^{n}{{n-1}\choose{m-1}}a_{m}S^{m},~{}b_{0}=a_{0}.$	(B.7b)

Note that (B.7a) is a power series in terms of the expansion variable

y\equiv x/(x+S),

(B.8)

such that the physical domain $x\in[0,\infty)$ maps to the transformed physical domain $y\in[0,1]$ . The derivation of (B.7b) is provided below.

First, equating the right-hand-sides of (B.6) and (B.7a) and using definition (B.8) to write in terms of $y$ , we have

\sum_{m=0}^{\infty}a_{m}\left({yS\over 1-y}\right)^{m}=\sum_{m=0}^{\infty}b_{m}y^{m},

(B.9)

where a deliberate switch has been made to using the index $m$ instead of $n$ so that we may extract $b_{n}$ at a specific value within the above (i.e., at $m=n$ ) as

b_{n}={1\over n!}\left\{{d^{n}\over dy^{n}}\left[\sum_{m=0}^{\infty}a_{m}\left({yS\over 1-y}\right)^{m}\right]\right\}_{y=0}.

(B.10)

Note that the above expression comes from applying the definition of a coefficient of a Taylor expansion about $y=0$ , which is precisely what the right-hand side of (B.9) is. Our task remains to show that the infinite series given by (B.10) may be rewritten as the finite series (B.7b). Interchanging the differentiation and summation operators and rewriting the inside as a product, (B.10) becomes

b_{n}={1\over n!}\left\{\sum_{m=0}^{\infty}a_{m}S^{m}{d^{n}\over dy^{n}}\left[y^{m}\left(1-y\right)^{-m}\right]\right\}_{y=0}.

(B.11)

In order to apply the generalized product rule to the $\left[y^{m}\left(1-y\right)^{-m}\right]$ term above, we first recognize that individual generalized derivatives of the two terms are

	${d^{j}(y^{m})\over dy^{j}}=\left[\prod_{\ell=0}^{j-1}\left(m-\ell\right)\right]y^{m-j}$	(B.12a)
and
	${d^{j}\left[\left(1-y\right)^{-m}\right]\over dy^{j}}=\left[\prod_{\ell=0}^{j-1}\left(m+\ell\right)\right]\left(1-y\right)^{-m-j}.$	(B.12b)

Substituting (B.12b) into the generalized product rule (B.3) to evaluate the derivative in (B.11), we obtain the expression

	$\displaystyle b_{n}=$	$\displaystyle{1\over n!}\left\{\sum_{m=0}^{\infty}a_{m}S^{m}\right.$
		$\displaystyle\left.\sum_{k=0}^{n}\left[{{n}\choose{k}}\left[\prod_{\ell=0}^{n-k-1}\!\!\!\!\!\left(m-\ell\right)\right]y^{m-n+k}\left[\prod_{\ell=0}^{k-1}\!\left(m+\ell\right)\right]\left(1-y\right)^{-m-k}\right]\right\}_{y=0}.$

Finally, we evaluate the above expression at $y=0$ and note that the only nonzero contributions will come from $y^{0}$ terms in the inner sum, which correspond with an index of $k=n-m$ . Since $k$ is a positive index, this tells us that $m$ cannot be greater than $n$ , thus placing an upper limit on the outer $m$ -sum. Setting $y=0$ , only retaining the $k=n-m$ term in the inner sum, and truncating the outer sum to $n$ terms, leads to

	$\displaystyle b_{n}=$	$\displaystyle\displaystyle~{}{1\over n!}\sum_{m=0}^{n}a_{m}S^{m}{{n}\choose{n-m}}~{}~{}~{}~{}~{}~{}\displaystyle\prod_{\ell=0}^{m-1}\left(m-\ell\right)~{}~{}~{}~{}~{}~{}\displaystyle\prod_{\ell=0}^{n-m-1}\left(m+\ell\right)$
	$\displaystyle=$	$\displaystyle\displaystyle~{}{1\over n!}\sum_{m=0}^{n}a_{m}S^{m}{n!\over m!\left(n-m\right)!}~{}~{}~{}~{}~{}~{}~{}~{}~{}m!~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\displaystyle{\left(n-1\right)!\over\left(m-1\right)!}$
	$\displaystyle=$	$\displaystyle\displaystyle~{}\sum_{m=0}^{n}a_{m}S^{m}{\left(n-1\right)!\over\left(m-1\right)!\left(n-m\right)!}$
	$\displaystyle=$	$\displaystyle\displaystyle~{}\sum_{m=0}^{n}{{n-1}\choose{m-1}}a_{m}S^{m},$

which is equivalent to (B.7b).

B.6 Estimation of the Radius of Convergence in an Euler Transformation

Typically, for nonlinear ODEs, the precise value of the radius of convergence of the power series solution is unknown, but can be estimated (via root test, ratio test, etc.) to within some tolerance. In such cases, we may estimate the true radius of convergence, $x_{s}$ , with an approximation, $S$ , as

S\equiv\epsilon x_{s},

(B.13)

where $\epsilon\in(0,2)$ . We are afforded this flexibility from Taylor’s theorem, which tells us that, through the mapping from $x$ to $y$ given by (B.8), a Taylor expansion about $y=0$ will converge over the full physical domain $y\in[0,1]$ as long as singularities lie outside of the circle $|y|=1$ in the complex $y$ plane. In other words, assuming the closest singularity is due to $x=-x_{s}$ , the following must hold for convergence of the corresponding Euler series:

\left|y\left(x=-x_{s}\right)\right|>1.

(B.14)

We know this to be true in the standard Euler series since then the left side of the inequality above becomes $\pm\infty$ , from the standard Euler gauge function (B.8). However, if we do not know the exact value of $x_{s}$ and substitute the Euler transformation

y\equiv x/(x+\epsilon x_{s}),

(B.15)

into (B.14), we obtain

\left|{1\over 1-\epsilon}\right|>1,

which tells us that $\epsilon\in(0,2)$ is a valid range to assure convergence of the transformed series. This means that one may use any arbitrarily small value for $S>0$ in (B.7b), as long as it is less than twice the value of the (unknown) exact $x_{s}$ (the modulus of the actual closest singularity). That said, as $\epsilon\to 1$ , the convergence rate of the Eulerized series is expected to improve, since—in this limit—the influence of the singularity is pushed further away from the mapped physical domain.

Appendix C Interface Solution for Small Slope

Equation (10) provides the leading order $\tilde{r}\to\infty$ asymptotic behavior of the interface solution given as:

\overline{h}(\tilde{r})\sim D{{K}_{0}}\left(\sqrt{B}\left(\tilde{r}+1\right)\right),

(C.1)

where $D$ is an unknown constant and ${{K}_{0}}\left(z\right)$ is the modified Bessel function of zeroth order. Note that this solution is valid for $d\overline{h}/d\tilde{r}\ll 1$ since $\overline{h}\to 0$ as $\tilde{r}\to\infty$ according to (2d). In a configuration for which $d\overline{h}/d\tilde{r}\ll 1$ for all $\tilde{r}\in\left(0,\infty\right)$ , (C.1) provides an excellent approximation to the interface solution with the constant, $D$ , determined as follows.

The boundary condition (2c) is rewritten here for convenience as:

{d\bar{h}\over d\tilde{r}}=-\cot\theta~{}\textrm{at}\ \tilde{r}=0.

(C.2)

Noting that:

{d{{K}_{0}}(z)\over dz}=-{{K}_{1}}(z),

where ${{K}_{1}}(z)$ is the modified Bessel function of first order, we differentiate (C.1) to obtain:

{d\bar{h}\over d\tilde{r}}\sim-D\sqrt{B}{{K}_{0}}\left(\sqrt{B}\left(\tilde{r}+1\right)\right)~{}\textrm{for}~{}~{}d\overline{h}/d\tilde{r}\ll 1.

(C.3)

Provided that $\cot\theta\ll 1$ , (C.3) is valid for all $\tilde{r}\in\left(0,\infty\right)$ , and thus the constraint (C.2) may be applied to determine $D$ as:

D={\cot\theta\over\sqrt{B}{{K}_{1}}\left(\sqrt{B}\right)}.

(C.4)

Finally, we substitute (C.4) into (C.1) to yield:

	$\overline{h}(\tilde{r})\sim{\cot\theta\over\sqrt{B}{{K}_{1}}\left(\sqrt{B}\right)}{{K}_{0}}\left(\sqrt{B}\left(\tilde{r}+1\right)\right)~{}\textrm{for}~{}~{}\cot\theta\ll 1,~{}\textrm{where}~{}~{}\tilde{r}\in\left(0,\infty\right).$	(C.5a)
In accordance with the restriction that $\cot\theta\ll 1$ , (C.5a) is valid for nearly horizontal interfaces for which values of $\theta$ lie near $\pi/2$ , regardless of the value of the Bond number, $B$ . The height of the interface, then at the cylindrical wall $(\tilde{r}=0)$ , is thus given as:
	$\overline{h}(0)\sim{\cot\theta\over\sqrt{B}}{{{K}_{0}}\left(\sqrt{B}\right)\over{{K}_{1}}\left(\sqrt{B}\right)},~{}\textrm{for}~{}~{}\cot\theta\ll 1.$	(C.5b)

The result (C.5b) serves as a useful analytical check on numerical and power series predictions for nearly horizontal interfaces.

Appendix D Computation of the Euler coefficients in (8j) via the transformed differential equation (7a)

After making the variable transformation (4b) in our governing ODE (2b), we obtain the transformed ODE (7a), rewritten here in a slightly different form as

g^{\prime\prime}={S^{2}\over(1-\eta)^{4}}\left\{{\displaystyle 2(1-\eta)^{3}\over\displaystyle S^{2}}g^{\prime}+{\displaystyle\left({\eta S\over 1-\eta}+1\right)\left({(1-\eta)^{4}\over S^{2}}g^{\prime 2}+1\right)^{3/2}Bg-{(1-\eta)^{6}\over S^{3}}g^{\prime 3}-{(1-\eta)^{2}\over S}g^{\prime}\over\displaystyle{\eta S\over 1-\eta}+1}\right\}.

(D.1)

Our goal is to obtain a power series solution of (D.1) in the form (8a), rewritten here for convenience as

g(\eta)=\sum_{n=0}^{\infty}\hat{A}_{n}\eta^{n}.

(D.2)

To do so, we note that, by direct differentiation, the first $n$ terms of the Taylor series of $g^{\prime}$ depend only on the first $n+1$ terms of $g$ . By Cauchy’s Product rule and JCP Miller’s Formula, the same is also true of powers and products of these terms. As such, the RHS of (D.1) will have an $n$ -th Taylor coefficient that can be expressed as a recursive expression involving the first $n+1$ Taylor coefficients of $g$ . Meanwhile, the LHS will give the $n+2$ Taylor coefficient of $g$ (up to a factor of $(n+1)(n+2)$ ), allowing one to recursive compute $\hat{A}_{n+2}$ from the previous Taylor coefficients $\{\hat{A}_{0},\hat{A}_{1},...,\hat{A}_{n+1}\}$ . More explicitly, we note that by the binomial theorem we have

{S^{2}\over(1-\eta)^{4}}=S^{2}\sum_{n=0}^{\infty}\binom{n+3}{3}\eta^{n}.

(D.3)

Hence by Cauchy’s product rule and direct differentiation if we write

{\displaystyle 2(1-\eta)^{3}\over\displaystyle S^{2}}g^{\prime}+{\displaystyle\left({\eta S\over 1-\eta}+1\right)\left({(1-\eta)^{4}\over S^{2}}g^{\prime 2}+1\right)^{3/2}Bg-{(1-\eta)^{6}\over S^{3}}g^{\prime 3}-{(1-\eta)^{2}\over S}g^{\prime}\over\displaystyle{\eta S\over 1-\eta}+1}=\sum_{n=0}^{\infty}w_{n}\eta^{n},

(D.4)

then the $w_{n}$ will have a recursive expression in terms of $\hat{A}_{n}$ only up to $\hat{A}_{n+1}$ , and equation (D.1) becomes

	$\displaystyle\sum_{n=0}^{\infty}(n+1)(n+2)\hat{A}_{n+2}$	$\displaystyle=S^{2}\sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}\binom{k+3}{3}w_{n-k}\right]\eta^{n},$		(D.5)
	$\displaystyle\hat{A}_{n+2}$	$\displaystyle={S^{2}\over(n+1)(n+2)}\sum_{k=0}^{n}\binom{k+3}{3}w_{n-k}.$		(D.6)

Continuing to work backwards, we can write by Cauchy’s Product Rule that

\displaystyle w_{n}

\displaystyle=f_{n}+\sum_{k=0}^{n}(q_{k}-t_{k}-u_{k})r_{n-k}

(D.7)

where

$\displaystyle{2(1-\eta)^{3}\over S^{2}}g^{\prime}$	$\displaystyle=\sum_{n=0}^{\infty}f_{n}\eta^{n}$	(D.8)
$\displaystyle\left({\eta S\over 1-\eta}+1\right)\left({(1-\eta)^{4}\over S^{2}}g^{\prime 2}+1\right)^{3/2}Bg$	$\displaystyle=\sum_{n=0}^{\infty}q_{n}\eta^{n}$	(D.9)
$\displaystyle{(1-\eta)^{6}\over S^{3}}g^{\prime 3}$	$\displaystyle=\sum_{n=0}^{\infty}t_{n}\eta^{n}$	(D.10)
$\displaystyle{(1-\eta)^{2}\over S}g^{\prime}$	$\displaystyle=\sum_{n=0}^{\infty}u_{n}\eta^{n}$	(D.11)
$\displaystyle\left({\eta S\over 1-\eta}+1\right)^{-1}$	$\displaystyle=\sum_{n=0}^{\infty}r_{n}\eta^{n}.$	(D.12)

We notice first that by term-by-term differentiation, Cauchy’s product rule, and the binomial formula we have the following formulae for $f_{n}$ and $u_{n}$ :

	$\displaystyle f_{n}$	$\displaystyle={2\over S^{2}}\sum_{k=0}^{n}(-1)^{k}\binom{3}{k}(n-k+1)\hat{A}_{n-k+1},$		(D.13)
	$\displaystyle u_{n}$	$\displaystyle={1\over S}\sum_{k=0}^{n}(-1)^{k}\binom{2}{k}(n-k+1)\hat{A}_{n-k+1}.$		(D.14)

Furthermore, by the formula for a geometric series, we have directly that

	$\displaystyle\left({\eta S\over 1-\eta}+1\right)^{-1}$	$\displaystyle={(1-\eta)\over 1-(1-S)\eta}$		(D.15)
	$\displaystyle r_{n}$	$\displaystyle=(1-S)^{n}-(1-\delta_{n,0})(1-S)^{n-1}.$		(D.16)

where the Kronecker notation is used such that $\delta_{n,0}$ is 0 when $n\neq 0$ and 1 when $n=0$ . The computation of $t_{n}$ in (D.10) could be obtained directly using JCP Miller’s formula on $g^{\prime 3}$ followed by Cauchy’s product rule; instead, we first build the formula for the term involving $g^{\prime 2}$ in (D.4), and then extract the result for $t_{n}$ from that. That is, we define

{(1-\eta)^{4}\over S^{2}}g^{\prime 2}=\sum_{n=0}^{\infty}d_{n}\eta^{n},

(D.17)

and by Cauchy’s product rule have

\displaystyle d_{n}

\displaystyle=\sum_{k=0}^{n}u_{k}u_{n-k}.

(D.18)

The expression for $t_{n}$ can be obtained by taking the Cauchy product of (D.17) and (D.11) to yield

\displaystyle t_{n}

\displaystyle=\sum_{k=0}^{n}d_{k}u_{n-k}.

(D.19)

From here it remains only for us to compute $q_{n}$ . We begin by writing

{(1-\eta)^{4}\over S^{2}}g^{\prime 2}+1=\sum_{n=0}^{\infty}\bar{d}_{n}\eta^{n},

(D.20)

then we have

\bar{d}_{n}=d_{n}+\delta_{n,0}.

(D.21)

Proceeding in the same fashion if we write

\left({(1-\eta)^{4}\over S^{2}}g^{\prime 2}+1\right)^{3/2}=\sum_{n=0}^{\infty}p_{n}\eta^{n},

(D.22)

then by JCP Miller’s formula we have

p_{0}=\bar{d}_{0}^{3/2}\hphantom{\sum\sum}p_{n>1}={1\over n\bar{d}_{0}}\sum_{k=1}^{n}\left({5\over 2}k-n\right)\bar{d}_{k}p_{n-k}.

(D.23)

Writing via geometric series that

	$\displaystyle{\eta S\over 1-\eta}+1$	$\displaystyle=\sum_{n=0}^{\infty}\bar{\ell}_{n}\eta^{n}$
	$\displaystyle\bar{\ell}_{n}$	$\displaystyle=S+\delta_{n,0}(1-S),$

then, if we write

\left({\eta S\over 1-\eta}+1\right)g=\sum_{n=0}^{\infty}\ell_{n}\eta^{n},

(D.24)

Cauchy’s product rule gives

\ell_{n}=\sum_{k=0}^{n}\bar{\ell}_{k}\hat{A}_{n-k}.

(D.25)

One final application of Cauchy’s Product Rule completes the recurrence by writing

q_{n}=B\sum_{k=0}^{n}\ell_{k}p_{n-k}.

(D.26)

In summary, the preceding coefficients $f_{n}$ , $q_{n}$ , $t_{n}$ , $u_{n}$ , $\ell_{n}$ , $d_{n}$ , $r_{n}$ , and $p_{n}$ are all utilized in the power series solution (8j) in the main text.

Appendix E Computation of the Euler coefficients in (11e) after inclusion of the modified Bessel function prefactor

Here, we re-sum to accelerate the convergence of the power series solution, making direct use of the previously derived series (8j) written in terms of the Euler transformed variable. The re-summation is enabled because the Euler coefficients in (8j) are computationally stable. This is in direct contrast to the original attempt to use the divergent series coefficients (3f) to construct the Euler coefficients using (6)–which leads to computational errors. Here, we rewrite for reference (8a) as

	$g(\eta)=\sum_{n=0}^{\infty}\hat{A}_{n}\eta^{n},$	(E.1a)
where the coefficients $\hat{A}_{n}$ are given by (8j) and
	$\eta(\tilde{r})={\tilde{r}\over\tilde{r}+S}.$	(E.1b)

We accelerate convergence of this summation by the inclusion of the leading asymptotic behavior in terms of the modified Bessel function of the second kind of zeroth order in accordance with (10) and (11a), written for reference here as

\bar{h}(\tilde{r})=K_{0}\left(\sqrt{B}(\tilde{r}+1)\right)\sum_{n=0}^{\infty}\bar{A}_{n}\left({\tilde{r}\over\tilde{r}+S}\right)^{n}.

(E.2)

By direct rearrangement of (E.1b), we have $\tilde{r}={S\eta\over 1-\eta}$ , and noting that $h(r)=h(r(\eta))=g(\eta)$ , we have

	$g(\eta)=K_{0}\left(\sqrt{B}\left({\eta S\over 1-\eta}+1\right)\right)\sum_{n=0}^{\infty}\bar{A}_{n}\eta^{n},$	(E.3a)
or equivalently:
	$g(\eta)\left[K_{0}\left(\sqrt{B}\left({\eta S\over 1-\eta}+1\right)\right)\right]^{-1}=\sum_{n=0}^{\infty}\bar{A}_{n}\eta^{n}.$	(E.3b)

Since $g(\eta)$ is known in terms of $\hat{A}_{n}$ from (E.1b), we can obtained a relationship between $\bar{A}_{n}$ (the expansion with the prefactor) and $\hat{A}_{n}$ (the original Euler expansion) as follows. We first define the argument of the Bessel function solution in (10) as

\sqrt{B}(\tilde{r}+1)=\sqrt{B}\left({\eta S\over 1-\eta}+1\right)\equiv y(\eta),

(E.4)

and express the modified Bessel function portion of (E.3b) as

K_{0}(y)=K_{0}\left(y(\eta)\right)=f(\eta)=\sum_{n=0}^{\infty}\bar{\zeta}_{n}\eta^{n}.

(E.5)

Using this expansion form, we employ JCP Miller’s formula (B.1b) to yield:

	$\left[K_{0}\left(\sqrt{B}\left({\eta S\over 1-\eta}+1\right)\right)\right]^{-1}=\sum_{n=0}^{\infty}\zeta_{n}\eta^{n},$		(E.6a)
	$\zeta_{n}={-1\over\bar{\zeta}_{0}}\sum_{k=1}^{n}\bar{\zeta}_{k}\zeta_{n-k}\hphantom{\sum\sum}\zeta_{0}=\bar{\zeta}_{0}^{-1}.$		(E.6b)

Subsequently, by inserting (E.6a) into (E.3b) and using Cauchy’s product rule (B.2), we obtain

\bar{A}_{n}=\sum_{k=0}^{n}\zeta_{k}\hat{A}_{n-k}.

(E.7)

The result (E.7) provides the necessary coefficients for the resummation–at this point, however, we have not yet determined the coefficients $\bar{\zeta}_{n}$ in (E.5), which we now obtain.

To proceed, we begin with the modified Bessel function governing equation (9) in the main text, rewritten in terms of $y$ as

{d^{2}K_{0}\over dy^{2}}+{1\over y}{dK_{0}\over dy}-K_{0}=0.

(E.8)

Our goal is to determine the Taylor expansion of $K_{0}(y(\eta))$ about $\eta=0$ , which is precisely the expression (E.5) that we require. Employing the chain rule formulae given by (B.5b) (letting $F=K_{0}$ ), substituting this result into (E.8), and solving for ${d^{2}f\over d\eta^{2}}$ leads to the expression

{d^{2}f\over d\eta^{2}}=\left({\hskip 5.0pt{d^{2}y\over d\eta^{2}}\hskip 5.0pt\over\hskip 5.0pt{dy\over d\eta}\hskip 5.0pt}~{}-~{}{\hskip 5.0pt{dy\over d\eta}\hskip 5.0pt\over\hskip 5.0pty\hskip 5.0pt}\right){df\over d\eta}+\left({dy\over d\eta}\right)^{2}f.

(E.9)

From here we will want to find the power series solution to (E.9) in the form (E.5). For the particular $y$ given by (E.4), we have the algebraic identities

$\displaystyle{\hskip 5.0pt{dy\over d\eta}\hskip 5.0pt\over\hskip 5.0pty\hskip 5.0pt}$	$\displaystyle={1\over 1-\eta}-{1-S\over 1-(1-S)\eta},$	(E.10)
$\displaystyle{\hskip 5.0pt{d^{2}y\over d\eta^{2}}\hskip 5.0pt\over\hskip 5.0pt{dy\over d\eta}\hskip 5.0pt}$	$\displaystyle={2\over 1-\eta},$	(E.11)
$\displaystyle\left({dy\over d\eta}\right)^{2}$	$\displaystyle={BS^{2}\over(1-\eta)^{4}}.$	(E.12)

Using the Binomial theorem, all of these expressions can be Taylor expanded explicitly as

$\displaystyle{\hskip 5.0pt{dy\over d\eta}\hskip 5.0pt\over\hskip 5.0pty\hskip 5.0pt}$	$\displaystyle={1\over 1-\eta}-{1-S\over 1-(1-S)\eta}=\sum_{n=0}^{\infty}\rho_{n}\eta^{n}=\sum_{n=0}^{\infty}(1-(1-S)^{n+1})\eta^{n},~{}\rho_{n}=1-(1-S)^{n+1},$	(E.13)
$\displaystyle{\hskip 5.0pt{d^{2}y\over d\eta^{2}}\hskip 5.0pt\over\hskip 5.0pt{dy\over d\eta}\hskip 5.0pt}$	$\displaystyle={2\over 1-\eta}=2\sum_{n=0}^{\infty}\eta^{n},$	(E.14)
$\displaystyle\left({dy\over d\eta}\right)^{2}$	$\displaystyle={BS^{2}\over(1-\eta)^{4}}=BS^{2}\sum_{n=0}^{\infty}\binom{n+3}{3}\eta^{n}.$	(E.15)

Furthermore, term-by-term differentiation of (E.5) gives

	$\displaystyle{df\over d\eta}$	$\displaystyle=\sum_{n=0}^{\infty}(n+1)\bar{\zeta}_{n+1}\eta^{n},$		(E.16)
	$\displaystyle{d^{2}f\over d\eta^{2}}$	$\displaystyle=\sum_{n=0}^{\infty}(n+1)(n+2)\bar{\zeta}_{n+2}\eta^{n}.$		(E.17)

We substitute the expansion (E.5) and results (E.13)-(E.17) into the ODE (E.9), and after making use of Cauchy’s product rule, we obtain

\sum_{n=0}^{\infty}(n+1)(n+2)\bar{\zeta}_{n+2}\eta^{n}=\sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}(k+1)\bar{\zeta}_{k+1}(2-\rho_{n-k})\right]\eta^{n}+BS^{2}\sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}\binom{k+3}{3}\bar{\zeta}_{n-k}\right]\eta^{n}.

(E.18)

Comparing like terms gives the recurrence relation

	$\bar{\zeta}_{n+2}={1\over(n+1)(n+2)}\sum_{k=0}^{n}\left\{(k+1)\bar{\zeta}_{k+1}(2-\rho_{n-k})+\binom{k+3}{3}\bar{\zeta}_{n-k}\right\},$	(E.19a)
where we have the initial conditions by direct substitution that
	$\bar{\zeta}_{0}=K_{0}\left(\sqrt{B}\right)~{}~{}\text{ and }~{}~{}\bar{\zeta}_{1}=-S\sqrt{B}K_{1}\left(\sqrt{B}\right).$	(E.19b)

The result (E.19b) can be inserted into (E.6b) and (E.7) which completely defines $\bar{A}_{n}$ in (E.2) and (11a), and the resummation is thus complete.

References

[1] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover, 1972, p. 376.
[2] N. S. Barlow, C. R. Stanton, N. Hill, S. J. Weinstein, and A. G. Cio, On the summation of divergent, truncated, and underspecified power series via asymptotic approximants, Q. J. Mech. Appl. Math., 70 (2017), pp. 21–48.
[3] G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge, 1967, ch. 1: The physical properties of fluids.
[4] L. Binetti, I. Del Villar, K. P. Dissanayake, A. Stankiewicz, T. Sun, K. T. V. Grattan, and L. S. M. Alwis, Monitoring of the critical meniscus of very low liquid volumes using an optical fiber sensor, IEEE Sensors Journal, 20 (2020).
[5] J. P. Boyd, The blasius function in the complex plane, Exper. Math., 8 (1999), pp. 381–394.
[6] R. V. Churchill, Complex Variables, McGraw-Hill, 1948, ch. VI: Power series.
[7] T. A. Driscoll, N. Hale, and L. Trefethen, Chebfun Guide, Patnuty Publications, Oxford, 2014.
[8] A. Ferguson, On the shape of the capillary surface formed by the external contact of a liquid with a cylinder of large radius, Philosophical Magazine and Journal of Science, 24 (1912), pp. 837–844.
[9] P. Flajolet and R. Sedgewick, Analytic Combinatorics, Cambridge University Press, 2009, ch. IV.3: Singularities and exponential growth of coefficients.
[10] I. S. Gradshteyn and I. M. Ryzhik, Table of integrals, series, and products: seventh edition, Academic Press, 2007.
[11] P. Henrici, Automatic computations with power series, JACM, 3 (1956), pp. 10–15.
[12] C. Huh and L. E. Scriven, Shapes of axisymmetric fluid interfaces of unbounded extent, J. of Colloid and Interface Science, 30 (1968), pp. 323–337.
[13] D. F. James, The meniscus on the outside of a small circular cylinder, J. Fluid Mech., 63 (1973), pp. 657–664.
[14] L. L. Lo, The meniscus on a needle - a lesson in matching, J. Fluid Mech., 132 (1983), pp. 65–78.
[15] N. Naghshineh, W. C. Reinberger, N. S. Barlow, M. A. Samaha, and S. J. Weinstein, On the use of asymptotically motivated gauge functions to obtain convergent series solutions to nonlinear ode. 2206.01536, 2022.
[16] Y. Tang and S. Cheng, The meniscus on the outside of a circular cylinder: From microscopic to macroscopic scales, Journal of Colloid and Interface Science, 533 (2019), pp. 401–408.
[17] M. Van Dyke, Perturbation Methods in Fluid Mechanics, Academic, 1964.
[18] D. D. White and J. A. Tallmadge, Static menisci on the outside of cylinders, J. Fluid Mech., 23 (1965), pp. 325–335.