1. Introduction
The study of floating objects has garnered attention since the middle of the 20th century due to various practical applications. These include the use of colloidal particles to stabilise emulsions (Binks & Horozov Reference Binks and Horozov2006; Tavacoli et al. Reference Tavacoli, Katgert, Kim, Cates and Clegg2012), the ability of insects and spiders to move on water surfaces (Gao & Jiang Reference Gao and Jiang2004; Bush & Hu Reference Bush and Hu2006) and capillary-driven self-assembly involving particles floating at a liquid–air interface (Whitesides & Boncheva Reference Whitesides and Boncheva2002; Whitesides & Grzybowski Reference Whitesides and Grzybowski2002; Delens, Collard & Vandewalle Reference Delens, Collard and Vandewalle2023). Liquid–air interface deflections by partially wet ellipsoids are known to cause chaining (Loudet et al. Reference Loudet, Alsayed, Zhang and Yodh2005, Reference Loudet, Yodh and Pouligny2006; van Nierop, Stijnman & Hilgenfeldt Reference van Nierop, Stijnman and Hilgenfeldt2005; Lewandowski et al. Reference Lewandowski, Bernate, Searson and Stebe2008) and result in aggregates of particles with complex shapes and selectively functionalised faces, which can be used to create tuneable devices like lenses with adjustable focal lengths (Bucaro et al. Reference Bucaro, Kolodner, Taylor, Sidorenko, Aizenberg and Krupenkin2009).
Self-assembly processes of objects that float on the surface of a liquid are driven by local liquid interface deformation around the objects, which depends on their shape, buoyancy and wetting properties (Kralchevsky & Nagayama Reference Kralchevsky and Nagayama1994; Danov & Kralchevsky Reference Danov and Kralchevsky2010; Gat & Gharib Reference Gat and Gharib2013; Poty, Lumay & Vandewalle Reference Poty, Lumay and Vandewalle2014). Based on the calculations of the surface height in the neighbourhood of a single bubble at rest on a horizontal surface of a fluid, Nicolson (Reference Nicolson1949) examined the implications of the fluid–fluid interface distortions for bubble–bubble interactions in particle rafts. Chan, Henry & White (Reference Chan, Henry and White1981), Kralchevsky & Denkov (Reference Kralchevsky and Denkov2001) and Vella & Mahadevan (Reference Vella and Mahadevan2005) found that the resulting shape of the liquid surface elevation, which was previously studied by Nicolson (Nicolson Reference Nicolson1949), causes the attraction/repulsion between bubbles.
Interaction between two floating particles occurs when the distance between them is of the order of the capillary length. For similar particles with similar capillary lengths, this means that the regions of interface deformation caused by each particle overlap, lead to mutual interaction (Loudet, Yodh & Pouligny Reference Loudet, Yodh and Pouligny2006). Particles with similar wetting properties and shapes tend to attract one another and form floating aggregates (Kralchevsky & Nagayama Reference Kralchevsky and Nagayama2000). Lucassen (Reference Lucassen1992) showed theoretically that interparticle interaction will take place as soon as the regions with deformed surfaces around two adjoining particles begin to overlap, where this interaction was calculated for particles with wetting perimeters of a sinusoidal shape. Yao et al. (Reference Yao, Botto, Cavallaro, Bleier, Garbin and Stebe2013) studied particles at interfaces with contact-line undulations having a wavelength significantly smaller than the characteristic particle size. By using theory and experiment, they showed that identical microparticles with features in phase attract each other, and microparticles with different wavelengths, under certain conditions, repel each other in the near field, leading to a measurable equilibrium separation. Equilibrium and mutual attraction or repulsion of objects supported by surface tension were also studied by Mansfield, Sepangi & Eastwood (Reference Mansfield, Sepangi and Eastwood1997). These phenomena are commonly called the Cheerios effect and were highlighted by Vella & Mahadevan (Reference Vella and Mahadevan2005) as well, who discussed some examples of particles’ attraction and repulsion taking into account their weight and buoyancy. Capillary interactions between particles trapped at the interface have been widely studied, both experimentally and theoretically. Botto et al. (Reference Botto, Lewandowski, Cavallaro and Stebe2012) provides a comprehensive discussion of the physical mechanisms governing interface deformations around anisotropic particles and their resulting pairwise interactions. More recently, Forth et al. (Reference Forth, Kim, Xie, Liu, Helms and Russell2019) expanded this by highlighting how such interactions can be harnessed to design dynamic and reconfigurable assemblies at complex liquid interfaces, connecting fundamental physics to functional material architectures. These studies emphasise the central role of interface shape in determining particle interactions and configurations.
Singh & Joseph (Reference Singh and Joseph2005) proposed a direct numerical simulation method to describe the nonlinear mechanism of capillary attraction between floating particles. Based on level sets, their numerical scheme allowed them to simulate the evolution of single spheres and disks to their equilibrium depth and the evolution to clusters of two and four spheres and two disks under capillary attraction. The interaction between adjacent floating particles may also be explained by using an analogy between capillary and electrostatic forces. Attachment of particles is introduced by ‘capillary charges’ referred as poles (Kralchevsky & Nagayama Reference Kralchevsky and Nagayama2000; Kralchevsky & Denkov Reference Kralchevsky and Denkov2001; Danov et al. Reference Danov, Kralchevsky, Naydenov and Brenn2005; Danov & Kralchevsky Reference Danov and Kralchevsky2010), which characterise the magnitude of the interfacial deformation and could be both positive and negative. For multipoles, the sign and magnitude of the capillary force depend on the interplay between the curvature signs of the interfacial deformations in the vicinity of the interacting particles. Capillary multipoles have been experimentally investigated by Bowden et al. (Reference Bowden, Terfort, Carbeck and Whitesides1997, Reference Bowden, Choi, Grzybowski and Whitesides1999) and Wolfe et al. (Reference Wolfe, Snead, Mao, Bowden and Whitesides2003), who observed that particles located along an interfacial surface meniscus with the same sign attract one another, while those along a meniscus with opposite signs repel one another. The strength and directionality of the interactions can be tailored by manipulating the heights of the faces, the pattern of the hydrophobic/hydrophilic faces and the densities of the three interacting phases. In a recent study, Eatson et al. (Reference Eatson, Morgan, Horozov and Buzza2024) explored the impact of particle shape on both short-range and long-range interactions by using colloidal building blocks composed of polygonal plates with uniform surface chemistry and wavy edges. Utilising minimum energy calculations alongside Monte Carlo simulations, which defined the interaction potential through multipoles, they showcased how colloidal building blocks could be customised to form diverse two-dimensional structures, such as hexagonal, honeycomb, open Kagome and quasicrystalline lattices. In addition, Soligno, Dijkstra & van Roij (Reference Soligno, Dijkstra and van Roij2016) numerically investigated, by energy minimisation, the capillary deformations induced by adsorbed cubes at fluid–fluid interfaces for various Young’s contact angles. In particular, they showed that strongly directional capillary interactions drive the cubes to self-assemble into hexagonal or graphene-like honeycomb lattices.
The notion of poles was also used in the theoretical studies of Kralchevsky & Denkov (Reference Kralchevsky and Denkov2001), Danov et al. (Reference Danov, Kralchevsky, Naydenov and Brenn2005) and Danov & Kralchevsky (Reference Danov and Kralchevsky2010), who theoretically predicted the meniscus shape around the particle under the assumption that it can be expanded into Fourier series. By integrating the linearised Young–Laplace equation along the midplane between the two particles and by using superposition, they derived a general integral expression for the capillary force and interaction energy between capillary multipoles of arbitrary orders.
A number of studies focused on the case of semi-submerged cylinders resting at the liquid interface. Janssens, Chaurasia & Fried (Reference Janssens, Chaurasia and Fried2017) introduced a study utilising the Young–Laplace equation to explore the configuration of the liquid–gas interface surrounding a floating particle. They examined the forces and torques on a partially submerged circular cylinder, considering surface tension disparities. Zhang, Zhou & Zhu (Reference Zhang, Zhou and Zhu2018) evaluated how surface tension influences both the vertical and rotational stability of floating cylinders with differing cross-sectional shapes by altering the surface’s curvature radius. Their research concluded that a diminished radius of curvature enhances rotational stability. Tan, Zhang & Zhou (Reference Tan, Zhang and Zhou2022) applied bifurcation theory to analyse the stability of multiple menisci formed in small concave arcs of cylinders. Additionally, Zhang & Zhou (Reference Zhang and Zhou2023) combined the Young–Laplace equation and bifurcation theory to study the stability of a floating cylinder situated between two parallel vertical plates and partly immersed in a liquid bath. Zhang et al. (Reference Zhang, Zhang, Tan and Zhou2023) developed a novel methodology to forecast the stability of floating objects of any shape, anchored on the contact angle and geometrical features at the contact point, including object concavity and convexity.
The motivation for this study is to provide a solution to the Young–Laplace equation with varying boundary height for non-circular boundaries, assuming pinned conditions at the boundary, and to evaluate the correlation between this theoretical solution and experimental data. Although domain perturbation is a well-established asymptotic technique for solving differential equations in general domains, it has not previously been applied to the particle–liquid configurations examined here. To the best of our knowledge, all previous works have utilised multipoles as a solution method. Furthermore, while the assumption of pinned conditions at particle edges simplifies the theoretical analysis, experimental validation of this boundary condition remains uncertain. The complexity is further amplified for non-circular geometries, where significant deviations from circularity could render first-order solutions inadequate. This study demonstrates that, for weakly non-circular geometries, the domain perturbation technique is effectively applicable and yields good agreement with experimental results, confirming its efficacy as a versatile theoretical tool. However, for geometries with larger deviations from circularity, higher-order terms become necessary to accurately capture the observed behaviours.
The insights gained from these findings are vital to advancing our understanding of the collective behaviour and aggregation dynamics of multiple anisotropic floating particles and can be used to explain how different particle shapes affect their aggregation processes. In § 2, we present our solution, which employs the domain perturbation methodology, as well as the leading- and first-order solutions for the liquid–gas interface for a general particle shape. Then, in § 3, we discuss our experimental procedure and our methodology for analysing the data. Section 4, provides a comparison between the experimentally obtained liquid–gas interface and the analytical results for a set of anisotropic particle shapes. We conclude our findings in § 5. Appendix A provides additional details of the derivation of the leading- and first-order solutions for the second-order partial differential equation describing the interface shape. In Appendix B we express the solution in the dimensional form. Appendix C contains expressions for the deviations of equilateral polygons from circular arcs of the auxiliary circle in several specific cases. Finally, Appendix D discusses the error estimation in contact inclination angle at the particle boundary due to the liquid height uncertainty and Appendix E provides the formula for the second-order term in the asymptotic expansion of the height of the liquid-gas interface.
2. Analysis
We examine the steady-state interfacial deformations induced by anisotropic particles held at the liquid–gas interface. By applying domain perturbations, we formulate a general solution for the height of the liquid–gas interface that is pinned to the circumference of anisotropic particles. The assumption of pinned boundary conditions was explicitly used in several studies concentrating on particles with pinned contact lines, e.g. Yao et al. (Reference Yao, Botto, Cavallaro, Bleier, Garbin and Stebe2013, Reference Yao, Sharifi-Mood, Liu and Stebe2015) and Prakash, Perrin & Botto (Reference Prakash, Perrin and Botto2024). While the assumption of the pinned boundary condition is reasonable, it is typically taken as granted. As described below, our study shows that pinned boundary conditions provide good agreement with experimental measurements even for non-circular particles with sinusoidal boundary conditions.
Consider a particle on top of a liquid where the particle consists of a hollow cylinder attached to a thin surrounding surface, as seen in figure 1(a). A cylindrical coordinate system
$(r,\theta ,z)$
, as depicted in figure 1(b) is adopted, where the origin
$(r,\theta )=(0,0)$
is located at the particle’s centre and where
$z=0$
is the liquid height sufficiently far from the particle, as shown in figure 1(c). We denote the outer radius of the circular particle by
$r_{\kern-1pt o}$
. In the case of non-circular particles,
$r_{\kern-1pt o}$
represents the radius of the auxiliary circle inscribed in the particle’s boundary circumference. Additionally,
$d_r(\theta )\geqslant 0$
denotes the deviation of the outer circumference of the boundary from the circular arc of the inscribed circle with radius
$r_{\kern-1pt o}$
. We also assume that the thickness of the surface along the circumference boundary is very small relative to the thickness of the particle. Thus, we can use the pinned boundary conditions along the rim of the particle. Accordingly, the rim of the models that we fabricated for the experiments was kept as thin as possible (see § 3).
Figure 1. (a) A 3D printed particle held at the liquid–gas interface, where the outer particle’s height is prescribed by
$0.5\sin (8\theta )$
mm. (b) A sketch of the rotated particle and (c) a side view of a particle which is held fixed on a liquid, where the height of its boundary is prescribed by
$f(\theta )$
.
We denote the height of the liquid–gas interface above
$z=0$
outside of the particle as
$h(r,\theta )$
. Here, we shall find an asymptotic approximation for
$h(r,\theta )$
in the domain
$(r,\theta )\in (r_{\kern-1pt o}+d_r(\theta ),\infty )\times [0,2\pi )$
, while accounting for various particle rim boundary shapes. For
$\theta \in [0,2\pi )$
the outer boundary of the particle is determined by the functions
$d_r(\theta )$
and
$f(\theta )$
, where
$f(\theta )$
denotes the particle’s height at the contact with the liquid–gas interface. Both these functions determine the shape of the interface and will serve in defining the boundary conditions for our problem.
The interface geometry is described by the Young–Laplace equation (Landau & Lifshitz Reference Landau and Lifshitz1987),
$\Delta p = 2\gamma \kappa$
, where
$\kappa$
denotes the mean curvature of the surface,
$\Delta p=p_a-p(h)$
,
$\gamma$
is the surface tension of the liquid,
${p}_{a}$
is the atmospheric pressure and
${p}({h})$
denotes the pressure acting on the liquid–gas interface
$ {p}({h}(r,\theta ))={p}_a- {\rho }{g} {h}({r},{\theta }),$
with
$\rho$
and
$g$
denoting the fluid density and gravitational constants, respectively. For a surface
$(r,\theta ,h(r,\theta ))$
it is well known (see e.g. Kralchevsky & Nagayama Reference Kralchevsky and Nagayama1994; Rotman Reference Rotman2009) that the mean curvature is given by
$\kappa =0.5\boldsymbol{\nabla }\,{\cdot}\,[ \boldsymbol{\nabla }h/\sqrt {1+(\boldsymbol{\nabla }h)^2} ]$
. Under the assumption that
$h(r,\theta )$
is sufficiently small (
$\vert h(r,\theta ) \vert \ll 1$
) and with bounded derivatives, the Young–Laplace equation can be linearised yielding that
$\gamma \Delta h=\Delta p$
, which may be expressed in polar coordinates as
\begin{equation} {\gamma } \left ( \frac {1}{{r}^2}\frac {\partial ^2{h}}{\partial {\theta }^2}+\frac {1}{r}\frac {\partial {h}}{\partial {r}}+\frac {\partial ^2{h}}{\partial {r}^2} \right )={\rho }{g} {h}({r},{\theta }). \end{equation}
Using the aforementioned definitions we can write the boundary conditions for our problem as
\begin{equation} h\vert _{r=r_{\kern-1pt o}+d_r(\theta )}={f}({\theta }), \end{equation}
where, for simplicity, we assume that
$f(\theta )$
is a 2
$\pi$
periodic function and is sufficiently regular so it is possible to expand it into Fourier series in
$\theta \in (0,2\pi )$
, namely
${f}({\theta })=w_0/2+\sum _{n=1}^{\infty } (q_n\sin {(n {\theta })} +w_n\cos {(n {\theta })} )$
, where
$q_n$
and
$w_n$
denote the Fourier coefficients.
First, we render the problem dimensionless, by defining the following transformations:
\begin{equation} \begin{aligned} &H=\frac {{h}}{h_c}, \quad R=\frac {{r}}{r_{\kern-1pt o}}, \quad D_r=\frac {d_r}{\varepsilon r_{\kern-1pt o}}, \quad Q_n=\frac {q_n}{h_c}, \quad W_n=\frac {w_n}{h_c}, \quad n=0,1,\ldots , \end{aligned} \end{equation}
where we assume that the maximal deviation of the circumference of the outer particle’s boundary relative to the inscribed circle is sufficiently small, namely
$\varepsilon :=\max _{\theta \in [0,2\pi )}{(|d_r(\theta )|/r_{\kern-1pt o})}$
satisfies
$0\lt \varepsilon \ll 1$
. Note that
$h_c \neq 0$
denotes the characteristic liquid–gas interface height, whose value is unimportant, since both the equation and the boundary conditions are linear and thus
$h_c$
cancels out. Moreover, note that the definition of
$r_{\kern-1pt o}$
as the radius of the inscribed circle implies that we solve the dimensionless problem in the domain
$R\gt 1$
. As we shall show in the continuation of this section, our solution involves infinite series which contain Bessel functions, and thus the requirement that
$R\gt 1$
is necessary to ensure convergence of the corresponding series.
Substituting the transformations in (2.3) into (2.1)–(2.2), we get that the dimensionless version of the problem is given by
\begin{align} &\qquad\quad \frac {\partial ^2{H}}{\partial {\theta }^2}+R \left (\frac {\partial {H}}{\partial {R}}+R \frac {\partial ^2{H}}{\partial {R}^2}\right )-{B}_{\textit{o}} {{R}^2} H=0 , \end{align}
\begin{align} &H\big \vert _{{R}=1+ \varepsilon D_r({\theta })}=\frac {W_0}{2}+\sum _{n=1}^{\infty } \left (W_n\cos {(n {\theta })} +Q_n\sin {(n {\theta })} \right )\!, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\quad H\big \vert _{R\to \infty }=0, \end{align}
where
${B}_{{o}}=g r_{\kern-1pt o}^2 ({\rho }_{L} - {\rho }_{G})/\gamma$
is the Bond number, which reflects the ratio of the gravitational force to the surface tension force. The Bond number typically appears in the analysis of interfacial phenomena such as of bubbles and drops (de Gennes, Brochard-Wyart & Quere Reference de Gennes, Brochard-Wyart and Quere2003). Here,
${\rho }_{L}$
denotes the liquid density and
${\rho }_{G}$
denotes the gas density, where the underlying assumption is that
${\rho }_{G}\ll {\rho }_{L}$
, thus we define
${\rho }_{L}-{\rho }_{G}={\rho }$
.
Note that the application of the boundary conditions at the particle’s boundary given in (2.4b
) to obtain a closed-form analytic solution is not straightforward in a general case since the shape of particle’s boundary may be rather complicated. To overcome this problem, we use the domain perturbations method (Leal Reference Leal2007) which provides an approximate mean to solve this problem under the assumption that
$0\lt \varepsilon \ll 1$
, allowing us to replace the exact boundary condition with an approximate boundary condition that is asymptotically equivalent for
$0\lt \varepsilon \ll 1$
, but applied at the boundary
$R=1$
.
Note that by the definition of
$\varepsilon$
it follows that
$ |D_r(\theta )|\leqslant 1$
. Hence, we may expand
$H\vert _{R=1+\varepsilon D_r(\theta )}$
into a Taylor series in
$\varepsilon$
around
$\varepsilon =0$
, as follows:
\begin{equation} H\big \vert _{{R}=1+\varepsilon D_r({\theta })}=H\big \vert _{R=1}+ \varepsilon D_r(\theta )\frac {\partial H}{\partial R}\bigg \vert _{R=1}+O(\varepsilon ^2). \end{equation}
Moreover, expanding
$H$
in regular asymptotic series in
$\varepsilon$
, namely
\begin{equation} H(R,\theta )=H_0(R,\theta )+\varepsilon H_1(R,\theta )+O(\varepsilon ^2), \end{equation}
and using these expansions in (2.4), we get that the leading-order problem with respect to
$\varepsilon$
is given by
\begin{align} &\quad\,\,\, R^2 \frac {\partial ^2H_0}{\partial {R}^2}+R\frac {\partial H_0}{\partial {R}}+\frac {\partial ^2 H_0}{\partial {\theta ^2}}-{B}_{{o}} {{R}^2} H_0=0, \end{align}
\begin{align} &H_0\big \vert _{{R}=1}=\frac {W_0}{2}+\sum _{n=1}^{\infty } \left (W_n\cos {(n {\theta })} +Q_n\sin {(n {\theta })} \right )\!, \end{align}
\begin{align} &\qquad\qquad\qquad\qquad H_0\big \vert _{R\to \infty }=0. \end{align}
Using separation of variables, we obtain that the general solution to the problem in (2.7) is given by
\begin{equation} H_0(R,\theta )=\frac {K_0\big(\sqrt {{B}_{{o}}}R\big)A_0}{2}+\sum _{n=1}^{\infty }K_n\big(\sqrt {{B}_{{o}}}R\big)\left [A_n \cos {(n\theta )}+B_n\sin {(n\theta )}\right ]\!, \end{equation}
where
$K_n$
,
$n=0,1,2,\ldots$
, denotes the modified Bessel function of the second kind (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010) and
$A_n$
,
$B_n$
are constants which are determined by using the boundary conditions in (2.7b
), namely
\begin{equation} A_n=\frac {W_n}{K_n\big(\sqrt {{B}_{{o}}}\big)} \quad \text{and} \quad B_n=\frac {Q_n}{K_n\big(\sqrt {{B}_{{o}}}\big)}. \end{equation}
Next, similarly to the discussion above, we get at the first,
$O(\varepsilon )$
, order the following problem:
\begin{equation} \begin{aligned} &R^2 \frac {\partial ^2 H_1}{\partial R^2}+R \frac {\partial H_1}{\partial R}+\frac {\partial ^2 H_1}{\partial \theta ^2}-{B}_{{o}} R^2 H_1=0,\\ &H_1\big \vert _{R=1}=-D_r(\theta ) \frac {\partial H_0}{\partial R}\bigg \vert _{R=1},\\ &H_1\big \vert _{R\to \infty }=0. \end{aligned} \end{equation}
Using again the separation of variables and the solution for
$H_0$
given in (2.8), we obtain that the first-order solution is given by
\begin{equation} H_1(R,\theta )=\frac {K_0\big(\sqrt {{B}_{{o}}}R\big)\tilde {A}_0}{2}+\sum _{n=1}^{\infty }K_n\big(\sqrt {{B}_{{o}}}R\big)\big[\tilde {A}_n \cos {(n\theta )}+\tilde {B}_n\sin {(n\theta )}\big], \end{equation}
where the constants
$\tilde {A}_n$
and
$\tilde {B}_n$
, which are determined by using the boundary conditions in (2.10), are given by
\begin{equation} \begin{aligned} &\tilde {A}_n=\frac {\sqrt {{B}_{{o}}}}{2\pi K_n\big(\sqrt {{B}_{{o}}}\big)}\int _{0}^{2\pi }G(\theta )\cos {(n\theta )}{\rm d}\theta , \quad n=0,1,2,\ldots ,\\ &\tilde {B}_n=\frac {\sqrt {{B}_{{o}}}}{2\pi K_n\big(\sqrt {{B}_{{o}}}\big)}\int _{0}^{2\pi }G(\theta )\sin {(n\theta )}{\rm d}\theta , \quad n=1,2,\ldots , \end{aligned} \end{equation}
where for brevity we use the following notation:
\begin{equation} \begin{aligned} G(\theta )&=D_r(\theta )\left \{\frac {\left (K_{-1}\big(\sqrt {{B}_{{o}}}\big)+K_{1}\big(\sqrt {{B}_{{o}}}\big)\right )A_0}{2}\right . \\ &\quad \left . +\sum _{m=1}^{\infty }\big(K_{m-1}\big(\sqrt {{B}_{{o}}}\big)+K_{m+1}\big(\sqrt {{B}_{{o}}}\big)\big)\bigl [A_m \cos {(m\theta )}+B_m\sin {(m\theta )}\bigr ]\vphantom{\frac {\left (K_{-1}\big(\sqrt {{B}_{{o}}}\big)+K_{1}\big(\sqrt {{B}_{{o}}}\big)\right )A_0}{2}}\right \}\!. \end{aligned} \end{equation}
Thus, we obtained an approximate solution for the liquid–gas interface height outside of a particle held at the liquid–gas interface with a non-circular shape and anisotropic height at the pinned contact line. The dimensionless solution is provided by
$H(R,\theta )=H_0(R,\theta )+\varepsilon H_1(R,\theta )$
, where the leading- and first-order terms
$H_0(R,\theta )$
and
$H_1(R,\theta )$
are given in (2.8)–(2.9) and (2.11)–(2.13), respectively. Note that, by our construction, the solution is valid for
$R\gt 1$
, which follows from defining
$r_{\kern-1pt o}$
as the radius of the inscribed circle. The reason for this choice is the convergence issues. It is known (Olver et al. (Reference Olver, Lozier, Boisvert and Clark2010), see formula 10.30.2) that
$K_{\nu } (R)\to \varGamma (\nu )/(2(R/2)^{\nu })$
when
$\nu \neq 0$
,
$R\to 0$
. Therefore, the sum
$\sum _{\nu =1}^{\infty }K_{\nu }(\sqrt {{B}_{{o}}}R)/K_{\nu }(\sqrt {{B}_{{o}}})$
is not convergent for
$|R|\lt 1$
. On the other hand (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010, see formula 10.25.3),
$K_{\nu }(R)\sim (\pi /2R)^{1/2}e^{-R}$
for
$R\to \infty$
and fixed
$\nu$
, so that
$\sum _{\nu =1}^{\infty }K_{\nu }(\sqrt {{B}_{{o}}}R)/K_{\nu }(\sqrt {{B}_{{o}}})$
converges for
$|R|\gt 1$
. Thus, defining
$r_{\kern-1pt o}$
as the radius of the inscribed circle, we obtain that the dimensionless radius satisfies
$R\gt 1$
, so that our solution is indeed convergent.
The dimensional form of the solution is given in Appendix B.
3. Experiments
Measuring relatively large surface deformations on a transparent fluid surface poses unique experimental challenges. While interferometric techniques offer high accuracy for very small surface deformations, they are generally unsuitable for large-amplitude measurements. This limitation arises because interferometry relies on precise optical phase differences, and large height variations typically exceed the coherence length of conventional interferometric set-ups, leading to ambiguity and loss of measurable interference fringes (Hinsch Reference Hinsch1978). Alternative techniques such as structured-light profilometry or standard optical three-dimensional scanning methods depend on a diffusely reflecting, opaque surface to reconstruct three-dimensional profiles accurately. Consequently, their application is severely limited by the transparent and reflective nature of water surfaces, resulting in unreliable or non-existent signal detection (Zhang Reference Zhang2018). Due to these constraints, the present study utilised a fluorescence-based laser sheet scanning approach combined with stereo imaging, enabling accurate measurement of significant liquid surface elevations around particles held at the liquid–gas interface without the aforementioned limitations.
A schematic representation of the experimental set-up is shown in figure 2(a). The set-up was developed to quantitatively measure the liquid surface elevation around particles placed at the free water surface, using a laser sheet scanning method combined with stereo imaging from two high-resolution cameras (LaVision CX-25). The measured surface elevations were subsequently compared with the theoretical asymptotic solution. Particles used in the experiment were fabricated using a Formlabs Form 3 stereolithography (SLA) 3D printer utilising Formlabs’ Tough 1500 UV-cured resin. Each particle featured a thick, stable base and a thin peripheral rim, specifically designed to minimise geometric distortions caused by resin shrinkage during printing and the curing processes. After complete UV curing, the particles underwent visual inspection and manual dimensional verification to confirm conformity with the intended design. Although precise deviations were not systematically quantified, no significant discrepancies were detected. Given Formlabs’ documented high-resolution capabilities and extensive prior experience with the Form 3 printer, we considered the fabricated geometries sufficiently accurate for the purposes of this study. Each particle was carefully placed onto the free water surface using a rigid supporting base rod to eliminate any motion during measurements, thus ensuring the best possible reconstruction quality. The tank was filled with purified water doped with a dilute concentration of Rhodamine 6G fluorescent dye. Although the water level was adjusted carefully, experimental constraints limited precise control to within
$\pm 0.05$
mm of the desired water height. The analytical solution assumes that the far-field liquid level is zero. Hence, considering this uncertainty will introduce an effective deviation to the relative height of the model edge, where the fluid is pinned. The uncertainty in the water level will effectively introduce a change to inclination angle of the meniscus close to the contact line. Appendix D presents a derivation of the inclination angle and examines the effects of the water level uncertainty in the experiments. It shows that the maximum uncertainty of
$\pm 0.05$
mm in the measurement of the liquid–gas interface at far field results in the error of the order of
$\pm 1\,\%$
in the inclination angle of the meniscus. Because the solution is decaying to the flat surface boundary, the maximal error in the solution is bounded by the measurement uncertainty applied to the boundary. In addition, analysis of the measurement presented in figure 2(b) shows that overall shape of the measured surface is correct, i.e. flat in this case. Consequently, it is reasonable to assume that the main findings and comparisons are mostly unaffected by this relatively small experimental limitation. A Nd:YAG laser (Quantel Evergreen-200), emitting at a wavelength of 532 nm, was utilised to excite the Rhodamine dye, resulting in fluorescent emission at approximately 550 nm. The laser beam was transformed into a thin laser sheet (
$\sim$
0.1 mm thickness) using a cylindrical lens, illuminating a vertical cross-section of the water’s surface around the particle. Fluorescent emissions were recorded using two high-resolution LaVision CX-25 cameras, each equipped with a 100 mm macro lens achieving a magnification of 54.3 (pixels mm–1). To enhance image quality and signal-to-noise ratio, each camera lens was fitted with a 550 nm long-pass optical filter (Edmund Optics #84-757), effectively eliminating scattered laser light and isolating the dye fluorescence. Cameras were symmetrically positioned
$\pm 30^\circ$
horizontally from the particle’s central axis and inclined at approximately
$10^\circ$
above the horizontal plane. This dual-camera configuration significantly reduced shadowing artefacts caused by elevated portions of the particle rim, thereby ensuring clear visualisation of the illuminated intersection line on the water’s surface. The laser sheet and both cameras were mounted on a precision motorised translation stage, enabling incremental spatial scanning in steps of 0.1 mm across the particle with a positional accuracy of approximately 0.001 mm. To maintain sharp imaging throughout the scan, each camera employed a Scheimpflug adapter to compensate effectively for focal-depth variations occurring during stage translation. As the cameras were placed on one side of the tank, only the front-facing half of the particle was captured. The recorded images were analysed using an in-house image processing algorithm designed to extract surface elevation profiles from each cross-sectional image. These profiles were then compiled to reconstruct a detailed three-dimensional representation of the water surface surrounding the particle. To validate the reconstruction accuracy, we performed repeated scans of the undisturbed free surface in the absence of any model. An example of such a scan is shown in figure 2(b). Across the scanned region and between both cameras (each camera scan region is marked by a different colour), the standard deviation in the measured surface height was found to be approximately 0.02 mm. This value is significantly smaller than the typical interface deformations observed in our experiments, and thus confirms that the measurement system provides sufficient precision for reliable comparison with the theoretical predictions. An illustrative example of the reconstructed water surface is shown in figure 2(c), where the elevation scale has been amplified by a factor of ten for improved visualisation. Prior to analysis, the mean far-field water level was subtracted from the measurements, establishing the undisturbed free-surface reference at zero elevation. While there are occasionally small regions which cannot be detected in the raw images due to light reflections or some impurities at the water’s surface, leading to small void regions in the experiments, the technique has proven reliable and the results are reproducible. To ensure that the measurement system performs well we quantified the standard deviation of the measured water surface elevation in the undisturbed far-field region away from the particle model, which consistently yielded an estimated precision of approximately 0.02 mm.
Figure 2. (a) The experimental set-up. The particles are held at the water–air interface, where the water is painted with a low concentration of Rhodamine 6G. A Nd:Yag 532 nm laser is used to excite the Rhodamine 6G, causing it to fluoresce at a wavelength of 550 nm. An optical set-up is used to form a thin laser sheet (
${\sim}0.1$
mm) that illuminates a cross-section of the water’s interface across the particle. The images were acquired by a camera with a band-pass filter attached to its lens, with a cutoff wavelength of 550 nm. (b) A raw scan of the free surface without any model present. The figure also indicates the regions captured by each camera: red indicates the area recorded by camera 1, and blue corresponds to camera 2. The standard deviation of the free-surface measurements is 0.02 mm, although this variation is too small to be discernible in the figure. (c) A raw scan result obtained for a flat hexagon, which is placed 0.5 mm higher than the air–liquid interface. The displacement is magnified by a factor of 10 for visual representation. These results are obtained by a camera which is placed at an angle of 30
$^\circ$
in relation to the central plane, marked in grey, and at an angle of 10
$^\circ$
to the horizon. The far side of the model was captured by the second camera.
Figure 3. The 3D printed particles made of Tough 1500 V1, with a reference auxiliary diameter circle of 10 mm held at the liquid–gas interface. The outer particles’ height is prescribed by
$0.5\sin (8\theta )$
mm for (a) circle, (b) triangle, (c) square, (e) pentagon, (f) hexagon and by
$0.5\sin (4\theta )$
mm for (d) circle.
Figure 3 shows a series of SLA 3D-printed particles held at the liquid–gas interface, each featuring a central reference circle of 10 mm diameter. The images illustrate various particle geometries: (a,d) circles, (b) triangle, (c) square, (e) pentagon and (f) hexagon. The outer edges of the particles incorporate sinusoidal undulations defined by amplitudes of
$0.5\sin (8\theta )$
mm for the polygonal shapes (b,c,e,f) and
$0.5\sin (4\theta )$
mm for the circular shape shown in (d). Off-axis illumination was employed to enhance visualisation of the particle surface features and the associated deformation of the liquid–gas interface around the rims. The visible interface distortions around each particle highlight the interplay between particle geometry and surface tension effects, emphasising both the fidelity achievable through SLA printing and the subtle liquid interface deformation investigated in this study.
4. Results
This section presents results obtained for two cases, circular and non-circular particles, as well as a comparison between the theoretical and experimental results. In the case of the circular boundary of the particle, the solution to the Young–Laplace equation is obtained analytically. However, in the case of a non-circular boundary, we use the domain perturbation technique, which yields an asymptotic solution to the Young–Laplace equation. This asymptotic solution depends on the difference in absolute value between the circular and the perturbed boundary. For all results presented in this section, except § 4.4, we employ the sum of the leading- and the first-order terms in the asymptotic expansion of
$H(R,\theta )$
, namely
$H(R,\theta )\approx H_0(R,\theta )+\varepsilon H_1(R,\theta )$
, where both,
$H_0(R,\theta )$
and
$H_1(R,\theta )$
, contain 20 modes in the Fourier series. In § 4.4, we show an example of convergence relative to the number of terms in the asymptotic expansion. This is achieved by comparing the results obtained using
$H^{(0)}(R,\theta )=H_0(R,\theta )$
,
$H^{(1)}(R,\theta )=H_0(R,\theta )+\varepsilon H_1(R,\theta )$
and
$H^{(2)}(R,\theta )=H_0(R,\theta )+\varepsilon H_1(R,\theta )+\varepsilon ^2 H_2(R,\theta )$
, where
$H_0(R,\theta )$
,
$H_1(R,\theta )$
and
$H_2(R,\theta )$
contain 20 modes in the Fourier series. Note that the expressions for
$H_0(R,\theta )$
and
$H_1(R,\theta )$
are given in § 2 and the expression for
$H_2(R,\theta )$
is derived in Appendix E. Convergence details relative to the number of modes are given in figure 9 in Appendix A.
4.1. Analytic results for circular particles
Figure 4(a,b) shows two examples of the liquid–gas interface obtained by using the solution given in (2.8), calculated for a circular particle with sinusoidal boundary conditions. As expected according to the boundary conditions, the interface is sinusoidal in the rotating angle
$\theta$
along the particle’s boundary, and it is flat and vanishing far enough away from the particle. This is reasonable since the Bessel functions appearing in the leading-order solution, which is exact in the case of the circular particle, are decreasing functions of
$r$
. The boundary conditions of figures 4(a) and 4(b) have a proportional value to
$\sin (4\theta )$
and
$\sin (8\theta )$
mm, respectively. Based on these figures, when the frequency of the sinusoidal wave increases, not only does the density of the peaks increase, but they also vanish at a faster rate as
$r$
increases towards the far boundary. The dependency of the influence of the particle’s distance on its rim undulation frequency is further emphasised in figure 4(c), where we show the maximal liquid–gas interface height (at any
$r$
) versus
$r$
for different frequencies of the sinusoidal wave prescribed along the outer boundary of the particle. According to figure 4(c), as the rim undulation frequency increases the liquid–gas interface flattens closer to the particle.
Figure 4. Analytical results for the liquid–gas interface for particles placed by held at the liquid–gas interface with
$r_{\kern-1pt o}=20$
mm, where the outer particles’ height is prescribed by: (a)
$0.5\sin (4\theta )$
mm and (b)
$0.5\sin (8\theta )$
mm. (c) Maximum heights (at any
$r$
) of the liquid–gas interfaces versus
$r$
obtained by solving the Young–Laplace equation with different boundary conditions, given by
$h\vert _{r=20}=0.5\sin (n\theta )$
mm, where the uppermost and the lowermost curves correspond to
$n=2$
and
$n=42$
, respectively.
Figure 5. The perturbation functions
$\varepsilon D_r(\theta )$
relative to the unit circle for (a) a hexagon, (b) a pentagon, (c) a square and (d) a triangle. The sketches of the corresponding circumferences of the particles in our experiments and in the asymptotic solutions are shown in the insets.
Figure 6. The liquid–gas interface obtained by the asymptotic solution for (a) hexagonal, (b) pentagonal (c) square and (d) triangular particles, whose outer heights are prescribed by
$0.5\sin (8\theta )$
mm and where
$r_{\kern-1pt o}=5$
mm in all cases.
4.2. Asymptotic results for non-circular particles
Here, we show the liquid–gas interfaces for the following particle’s shapes: equilateral hexagon, pentagon, square and triangle. As presented in the insets in figure 5, we define an auxiliary inscribed circle and scale the problem by using its radius
$r_{\kern-1.2pt o}$
. To obtain the asymptotic solution by using the domain perturbation technique, we refer to the circumferences of these shapes as to a small deviation
$\varepsilon D_r(\theta )$
from a circle. Note that as the sides number of an equilateral particle decreases, the deviation from the circle increases. Figure 5 illustrates, however, that the largest deviation between the studied shapes, which is equal to one, occurs at the corners of the triangle. As a result, our method becomes inapplicable, regardless of the number of terms in the asymptotic expansion, at the vicinity of the triangle’s corners only.
More specifically, for each
$\theta$
,
$\varepsilon D_r(\theta )$
denotes the distance between the corresponding shape and the inscribed circle, expressed in radial coordinates. The graphs of
$\varepsilon D_r(\theta )$
, are shown in figure 5 and the corresponding expressions for
$\varepsilon D_r(\theta )$
are given in Appendix C. As expected, the maximum of
$\varepsilon D_r(\theta )$
increases as we move from a hexagon to a pentagon, from a pentagon to a square and from a square to a triangle, which reflects the fact that the circumference of the hexagon is the closest (among the discussed shapes) to those of the inscribed circle and the circumference of the triangle is the furthest from those of the inscribed circle.
Figure 6 presents results for the liquid–gas interface obtained using the asymptotic solution given in (2.6) including both the leading- and first-order solutions prescribed in (2.8) and (2.11), respectively. The asymptotic results were calculated for equilateral hexagonal, pentagonal, square and triangular particles with sinusoidal boundary conditions. The liquid–gas interface is pinned to the particle, and hence its boundary conditions are determined according to the particle’s boundary outer height. As the particle circumference approaches the circle, the liquid–gas interface near the particle’s boundary approaches the sinusoidal wave (for a given
$r$
, as a function of
$\theta$
), as the asymptotic approximation becomes more accurate. As we move along the particle’s boundaries towards the corners, where the deviation from a circle is maximal, the amplitude of the sinusoidal-like wave decreases, since the asymptotic approximation is less accurate there. When
$0\lt \varepsilon \lt 1$
, this inaccuracy may be improved by adding additional terms of higher orders to the asymptotic expansion in (2.6), (for an additional discussion, see § 4.4). Moreover, similar to the circular particles, the maximal height (in absolute value) of the liquid–gas interface decreases along the radial direction as we move from the particle’s boundary until the interface flattens.
4.3. Comparison between theoretical and experimental results
In figure 7 we show a comparison between experimental (upper half of all panels) and theoretical (lower half of all panels) results for the liquid–gas interfaces for the following particle shapes: an equilateral flat hexagonal particle which is above the liquid level by 0.5 mm, a circular particle where the outer particle’s height is prescribed by
$0.5\sin (8\theta )$
mm and an equilateral hexagonal particle where the outer particle’s height is prescribed by
$0.1+0.5\sin (8\theta )$
mm. The circumference radius of the inscribed circle is 20 mm for all particles. As can be seen, in all cases there is a quantitative agreement between theory and experiment, except probably in close vicinity of the hexagon’s corners. This inaccuracy in theoretical results might be easily improved by taking into account more terms in the asymptotic expansion of
$h(r,\theta )$
, as discussed in § 4.4.
Figure 7. A comparison between experimental (upper half of the panels a–c) and theoretical (lower half of the panels a–c) results for the liquid–gas interface in three cases: (a) a hexagonal particle held at the liquid–gas interface at the height of 0.5 mm, (b) a circular particle with the outer height of
$0.5\sin {(8\theta )}$
mm and (c) a hexagonal particle with the outer height of
$0.1+0.5\sin {(8\theta )}$
mm. In all panels the radius of the inscribed circle is
$r_{\kern-1pt o}=20$
mm and the height of the particles is relative to the liquid–gas interface sufficiently far from the particle. The colour bar refers to both cases (experimental and theoretical) for all panels.
Figure 8. An example of the convergence analysis relative to the number of terms in the asymptotic expansion for a square-shaped particle with an outer height of
$0.5\sin {(8\theta )}$
mm using (a) only the leading-order term
$H^{(0)}(R,\theta )$
, (b) the sum of the leading- and the first-order terms
$H^{(1)}(R,\theta )$
, (c) the sum of the leading-, the first- and the second-order terms
$H^{(2)}(R,\theta )$
and (d) a comparison between experiment (upper half of the panel) and
$H^{(2)}(R,\theta )$
(lower half of the panel). In all panels the radius of the inscribed circle is
$r_{\kern-1pt o}=20$
mm.
4.4. Limitations of the method
As mentioned in § 2, the asymptotic nature of the theoretical approach employed in this study implies the need to use small values of
$\varepsilon$
,
$0\lt \varepsilon \ll 1$
. Moreover, note that the error of the asymptotic expansion
$H^{(M)}(R,\theta ):=\sum _{m=0}^M \varepsilon ^m H_m$
, where
$M=0,1,2,3\,\ldots$
, is of order
$O(\varepsilon ^{M+1})$
, namely
$(H(R,\theta )-H^{(M)}(R,\theta ))=O(\varepsilon ^{M+1})$
. Thus, as
$\varepsilon$
increases towards 1, we get a larger deviation of the asymptotic expansion
$H^{(1)}(R,\theta )$
discussed above from
$H(R,\theta )$
, which implies that, in order to get satisfactory accuracy for relatively large
$0\lt \varepsilon \lt 1$
, higher-order terms in the asymptotic expansion are required.
To illustrate the effects of the expansion order on the accuracy of the solution, we performed an experiment with a square-shaped particle, where the radius of the inscribed circle is 20 mm and the particle’s outer height is given by
$0.5\sin {(8\theta )}$
mm. It is easy to verify that, for a general square-shaped particle,
$\varepsilon =1/\cos {(\pi /4)}-1\approx 0.414$
, which is visualised by the maximum values in figure 5(c). This value of
$\varepsilon$
is relatively large and, as shown in figure 8(a), results in a significant error in the leading-order term solution
$H^{(0)}(R,\theta )=H_0(R,\theta )$
given in (2.8)–(2.9). As shown in figure 8(b), this error considerably decreases when taking into account the sum of the leading- and first-order terms
$H^{(1)}(R,\theta )=H_0(R,\theta )+\varepsilon H_1(R,\theta )$
, where
$H_1(R,\theta )$
is given in (2.11)–(2.13). However, the decrease in error is not sufficient, since the solution
$H^{(1)}(R,\theta )$
misses the correct behaviour near the corners of the square, where
$\varepsilon D_r(\theta )$
is the largest and approaches
$0.414$
. Thus, figure 8(c) depicts the solution
$H^{(2)}(R,\theta )$
, which takes into account the sum of three terms: the leading-, the first- and the second-order terms, where the expression for
$H_2(R,\theta )$
is derived in Appendix E. Figure 8(d) visualises the agreement between the experimental measurement of the liquid–gas interface for a square-shaped particle (upper half of the panel) and the solution
$H^{(2)}(R,\theta )$
(lower half of the panel). According to this figure, it is obvious that including additional terms in the asymptotic expansion significantly improves the accuracy of the solution, and in particular, taking 3 terms in the expansion allows us to get the correct behaviour along the sides of a square, including its corners. Note that there is an excellent agreement between the asymptotic solution with 3 terms in the expansion and the experimental measurements.
In summary, figure 8 demonstrates that the theoretical method proposed in the current study can be applied to more complex geometries. However, if
$\varepsilon$
increases, additional orders in the asymptotic expansion are required. Specifically, while asymptotic solutions are typically employed for
$0\lt \varepsilon \ll 1$
, here, we show that, by employing higher-order terms in the expansion, we may obtain satisfactory results even for a relatively large
$\varepsilon$
,
$\varepsilon \approx 0.414$
. Note that, for the hexagon shown in figure 7, the value of
$\varepsilon$
was relatively small,
$\varepsilon \approx 0.15$
, and hence even the solution
$H^{(1)}(R,\theta )$
composed of the leading and first terms only agrees well with the experiments.
5. Conclusions
This study leverages the domain perturbation technique to analyse the liquid–gas interface deformation in the vicinity of an anisotropic particle held at the liquid–gas interface, where the liquid–gas interface is assumed to be pinned to the particle’s edges. We obtain a closed-form solution for particles with non-circular complex shapes and varying contact-line heights. This technique provides an appealing method for non-circular domains, complementing previously used techniques, such as multipoles, which were implemented for circular or ellipsoidal particles (Kralchevsky & Denkov Reference Kralchevsky and Denkov2001; Danov et al. Reference Danov, Kralchevsky, Naydenov and Brenn2005; Danov & Kralchevsky Reference Danov and Kralchevsky2010).
Our theoretical results were validated by experiments performed using high-resolution laser scanning that captured the deformation of the liquid–gas interface around 3D printed particles held at the liquid–gas interface. We tested circular, hexagonal and square-shaped particles with flat and sinusoidal rims which are placed at, or at a shifted height from, the liquid–gas interface. The overall comparison shown in figures 7 and 8 reveals that the surface shape around the particle and the consequent effects it has on its vicinity seem to be well correlated between experiments and theory. This provides us with confidence that the analysis method and the assumption of pinned boundary conditions employed in this study provide a good estimation of the effects of particles on the surrounding liquid–gas interface. Note that we decided to concentrate on polygonal shapes in order to highlight the robustness of our approach, which succeeds in yielding a good agreement between theory and experiment, even though polygons contain sharp edges, which generally may complicate the theoretical modelling. For the cases of polygons that have relatively large deviations from a circle, such as a square, taking into account additional terms in the asymptotic expansions is required to resolve the accuracy issues in the vicinity of the corners, as shown in figure 8.
Building on this work, future studies can generalise the present approach to systems of multiple floating anisotropic particles, enabling analysis of their collective capillary interactions and self-assembly dynamics. By combining the domain perturbation method with bipolar coordinate techniques, as discussed by Kralchevsky & Nagayama (Reference Kralchevsky and Nagayama2000), the model can be extended to derive explicit solutions for pairwise interface deformations and interaction forces induced by complex particle shapes. This will help to clarify how particle anisotropy and orientation dictate not only attractive or repulsive capillary forces but also significant capillary torques that can induce mutual rotations, an effect discussed by Koens et al. (Reference Koens, Wang, Sitti and Lauga2019), who highlighted the distinct roles of near-field and far-field effects in the self-assembly of magnetic capillary disks. Moreover, adopting superposition approximations similar to those employed by Singer (Reference Singer2004) and Beatus, Bar-Ziv & Tlusty (Reference Beatus, Bar-Ziv and Tlusty2012) could simplify the modelling complex particle self-assembly. Collectively, such advances will offer a more complete picture of anisotropic capillary interactions and inform the design of self-assembled interfacial particle networks, paving the way toward programmable interfacial materials with tailored assembly properties (Poty et al. Reference Poty, Lumay and Vandewalle2014).
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2025.10727.
Acknowledgements
J.T.J. is partially supported by a Technion Post-Doctoral fellowship. This work is partially supported by the Zuckerman STEM Leadership Program.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Approximate solution of the problem in (2.4)
We begin by solving the leading-order problem given in (2.7). We assume that the solution may be obtained by separation of variables, namely that there exist two functions,
$\mathcal{R}(R)$
and
$\varTheta (\theta )$
, so that the solution of (2.7) may be expressed as
$H_0(R,\theta )=\mathcal{R}(R)\varTheta (\theta )$
. Substituting this separation of variables into equation (2.7a
) and employing standard mathematical manipulations, we obtain the following equation:
\begin{equation} \frac {R^2}{\mathcal{R}}\frac {\partial ^2 \mathcal{R}}{\partial R^2}+\frac {R}{\mathcal{R}}\frac {\partial \mathcal{R}}{\partial R}-{B}_{{o}}R^2=-\frac {1}{\varTheta }\frac {\partial ^2 \varTheta }{\partial \theta ^2}. \end{equation}
Since the left-hand side of the equation in (A1) depends only on
$R$
and the right-hand side of this equation depends only on
$\theta$
, there exists a constant
$\lambda$
, so that
\begin{equation} \frac {R^2}{\mathcal{R}}\frac {\partial ^2 \mathcal{R}}{\partial R^2}+\frac {R}{\mathcal{R}}\frac {\partial \mathcal{R}}{\partial R}-{B}_{{o}}R^2=-\frac {1}{\varTheta }\frac {\partial ^2 \varTheta }{\partial \theta ^2}=\lambda . \end{equation}
Considering the equation for
$\varTheta$
in (A2) and noting that, according to the boundary condition in (2.7b
)
$\varTheta (\theta )$
is 2
$\pi$
periodic in
$\theta$
, we may conclude that the eigenvalues
$\lambda$
are necessarily of the form
$\lambda =n^2$
, for
$n=0,1,2,\ldots$
, and the corresponding eigenfunctions are given by
\begin{equation} \varTheta _n(\theta )=A_n \cos {(n\theta )}+B_n\sin {(n\theta )}, \quad n=0,1,2,\ldots . \end{equation}
Then, for these eigenvalues the resulting equations for
$\mathcal{R}$
are given by
\begin{equation} R^2\frac {\partial ^2 \mathcal{R}}{\partial R^2}+R\frac {\partial \mathcal{R}}{\partial R}-\big({B}_{{o}}R^2+n^2\big)\mathcal{R}=0, \quad n=0,1,2,\ldots . \end{equation}
This type of equation is known as a modified Bessel equation and its solutions are called modified Bessel functions (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010). Due to the fact that the modified Bessel function of the first kind,
$I_n(R)$
, is unbounded as
$R$
tends to infinity, and we seek a solution that satisfies (2.7c
), we may only use the modified Bessel function of the second kind,
$K_{n}(R)$
(which tends to 0 when
$R\to \infty$
). Hence, for the problem in (A4) we will use the eigenfunctions
\begin{equation} \mathcal{R}(R)=K_{n}\big(\sqrt {{B}_{{o}}}R\big), \quad n=0,1,2,\ldots , \end{equation}
so that the general solution to the equation in (2.7a
) is given by (2.8). Imposing on the solution in (2.8) the boundary conditions in (2.7b
), we obtain that the constants
$A_n$
and
$B_n$
are as prescribed in (2.9).
Note that since the equation for
$H_1$
as given in (2.10) is exactly of the same form as the equation for
$H_0$
, by using the same procedure as described above, we get the same general solution for
$H_1$
as for
$H_0$
, see (2.11). Imposing the boundary conditions in (2.10) on the general solution for
$H_1$
in (2.11) and substituting (2.8), we get that
\begin{equation} \begin{aligned} &\frac {K_0\big(\sqrt {{B}_{{o}}}\big)\tilde {A}_0}{2}+\sum _{n=1}^{\infty }K_n\big(\sqrt {{B}_{{o}}}\big)\bigl [\tilde {A}_n \cos {(n\theta )}+\tilde {B}_n\sin {(n\theta )}\bigr ]\\ &=-D_r(\theta )\!\left \{\frac {\partial K_0\left(\!\sqrt {{B}_{{o}}}R\right)}{\partial R}\bigg \vert _{R=1}\frac {A_0}{2}+\!\sum _{n=1}^{\infty }\frac {\partial K_n\left(\!\sqrt {{B}_{{o}}}R\right)}{\partial R}\bigg \vert _{R=1}\!\left [A_n \cos {(n\theta )}+B_n\sin {(n\theta )}\right ]\!\right \}\!, \end{aligned} \end{equation}
where
$A_n$
and
$B_n$
are given in (2.9).
In order to get explicit expressions for
$\tilde {A}_n$
and
$\tilde {B}_n$
, we use the following property of the Bessel functions
$K_n$
(Olver et al. Reference Olver, Lozier, Boisvert and Clark2010):
\begin{equation} \frac {\partial K_n\big(\sqrt {{B}_{{o}}}R\big)}{\partial R}=-\frac {\sqrt {{B}_{{o}}}}{2}\big(K_{n-1}\big(\sqrt {{B}_{{o}}}R\big)+K_{n+1}\big(\sqrt {{B}_{{o}}}R\big)\big). \end{equation}
Next, expanding the right-hand side of the equation in (A6) into Fourier series in the interval
$\theta \in (0,2\pi )$
, we get exactly the expressions which were given in (2.12)–(2.13).
In order to check Fourier series convergence relative to the number of modes in the approximate solution, containing both the leading- and the first-order terms, we calculated the error which we denote by Err and define as
\begin{equation} \text{Err}=h_c||H(r,\theta ;K=50)-H(r,\theta ;K)||_{\infty }, \end{equation}
where
$H(r,\theta )$
is as given in (2.6), (2.8) and (2.11), with the number of terms in the truncated series denoted by
$K$
and where
$||\boldsymbol{\cdot }||_{\infty }$
denotes the
$L^{\infty }-$
norm in the two-dimensional domain
$(r,\theta )\in (r_{\kern-1pt o},3.5r_{\kern-1pt o})\times [0,2\pi )$
. In figure 9, we show the error graphs for the results shown in figure 6, corresponding to the hexagon, pentagon, square and triangle versus
$K$
. According to figure 9, it can be seen that in all cases the Fourier series converges exponentially relative to the number of modes in the series.
Figure 9. The error, Err defined in (A8), arising from the truncation of a finite number of modes in the Fourier series of the asymptotic solution
$H(r,\theta )$
versus the number of modes,
$K$
. The full circles represent results obtained by simulation for the hexagon, pentagon, square and triangle, which correspond to the results shown in figure 6, and the curves, which were obtained by least square fitting to the simulation results, are to guide the eyes.
Appendix B. The solution in dimensional form
To obtain the dimensional solution, we first recall that the boundary conditions in dimensional notation are given by
\begin{equation} h\vert _{r=r_{\kern-1pt o}+d_r(\theta )}=\frac {w_0}{2}+\sum _{n=1}^{\infty } \left (q_n\sin {(n {\theta })} +w_n\cos {(n {\theta })} \right )\!. \end{equation}
Now, we substitute the variables transformation given in (2.3) back into the dimensionless solution given in (2.8)–(2.9) and (2.11)–(2.13) and use the definition of
${B}_{{o}}$
. Thus, we get that
$h(r,\theta )=h_0(r,\theta )+\tilde {h}_1(r,\theta )$
, where we denote for brevity
\begin{align} \tilde {h}_1(r,\theta )=\frac {\max _{\theta \in (0,2\pi )}|d_r(\theta)|}{r_{\kern-1pt o}}h_1(r,\theta ). \end{align}
As to the leading-order solution, we get for
$r\gt r_{\kern-1pt o}$
that
\begin{equation} h_0(r,\theta )=\frac {K_0\left(\sqrt {g \rho /\gamma }r\right)a_0}{2}+\sum _{n=1}^{\infty }K_n\big(\sqrt {g \rho /\gamma }r\big)\left [a_n \cos {(n\theta )}+b_n\sin {(n\theta )}\right ]\!, \end{equation}
where
$a_n$
and
$b_n$
, the dimensional versions of
$A_n$
and
$B_n$
, are given by
\begin{equation} a_n=\frac {w_n}{K_n\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)} \quad \text{and} \quad b_n=\frac {q_n}{K_n\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)}. \end{equation}
As to the first-order solution, we get for
$r\gt r_{\kern-1pt o}$
that
\begin{equation} \tilde {h}_1(r,\theta )=\frac {K_0\big(\sqrt {g \rho /\gamma }r\big)\tilde {a}_0}{2}+\sum _{n=1}^{\infty }K_n\big(\sqrt {g \rho /\gamma }r\big)\big[\tilde {a}_n \cos {(n\theta )}+\tilde {b}_n\sin {(n\theta )}\big], \end{equation}
where
$\tilde {a}_n$
and
$\tilde {b}_n$
, the dimensional versions of
$\tilde {A}_n$
and
$\tilde {B}_n$
, are given by
\begin{equation} \begin{aligned} &\tilde {a}_n=\frac {\sqrt {g \rho /\gamma }}{2\pi K_n\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)}\int _{0}^{2\pi }g(\theta )\cos {(n\theta )}{\rm d}\theta , \quad n=0,1,2,\ldots ,\\ &\tilde {b}_n=\frac {\sqrt {g \rho /\gamma }}{2\pi K_n\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)}\int _{0}^{2\pi }g(\theta )\sin {(n\theta )}{\rm d}\theta , \quad n=1,2,\ldots , \end{aligned} \end{equation}
and where
$g(\theta )$
, the dimensional version of
$G(\theta )$
, is given by
\begin{equation} \begin{aligned} g(\theta )&=d_r(\theta )\Bigl \{\frac {\left (K_{-1}\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)+K_{1}\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)\right )a_0}{2}\\ &\quad +\sum _{m=1}^{\infty }\big(K_{m-1}\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)+K_{m+1}\big(\sqrt {g \rho /\gamma }r_{\kern-1pt o}\big)\big)\bigl [a_m \cos {(m\theta )}+b_m\sin {(m\theta )}\bigr ]\Bigr \}. \end{aligned} \end{equation}
Appendix C. The expressions for
$\boldsymbol{\varepsilon D_r(\theta )}$
for different shapes
In this section we shall give the expressions for
$\varepsilon D_r(\theta )$
in four cases: equilateral hexagon, pentagon, square and triangle, which are shown in the insets of figure 5. In all cases the dimensionless radius of the inscribed circle is assumed to satisfy
$R=1$
.
Let us start our discussion from equilateral hexagon, which is shown in the inset of figure 5(a). It is easy to verify that the sides of the hexagon, which we denote by
$L_i$
,
$i=1,2,\ldots ,6$
(as indicated in the inset of figure 5
a), are given by
\begin{equation} \begin{aligned} &L_1:\quad y_1=\sqrt {3} x-2, \\ &L_2:\quad y_2=-\sqrt {3} x-2, \\ &L_3:\quad y_3=-\sqrt {3}x+2, \\ &L_4:\quad y_4=\sqrt {3}x+2, \\ &L_5:\quad y_5=1, &&x=\frac {y_5}{\tan {(\theta )}}=\cot {(\theta )},\\ &L_6:\quad y_6=-1, &&x=\frac {y_6}{\tan {(\theta )}}=-\cot {(\theta )}. \end{aligned} \end{equation}
Now, if we wish to express the sides
$L_i$
,
$i=1,2,3,4$
, in polar coordinates, namely in the form
$(r_i\cos {(\theta )},r_i\sin {(\theta )})$
, where
$\theta$
is within the corresponding ranges, which are denoted by
$I_i$
,
$i=1,2,3,4$
, respectively, and are given in (C7), we need only to find
$r_i$
. However, according to (C1), it follows that
\begin{equation} r_i(\theta )=\left \{\begin{aligned} &\frac {2}{\sqrt {3}\cos {(\theta )}-\sin {(\theta )}}, &&i=1,\\ &-\frac {2}{\sqrt {3}\cos {(\theta )}+\sin {(\theta )}}, &&i=2,\\ &\frac {2}{\sqrt {3}\cos {(\theta )}+\sin {(\theta )}}, &&i=3,\\ &-\frac {2}{\sqrt {3}\cos {(\theta )}-\sin {(\theta )}}, &&i=4. \end{aligned}\right . \end{equation}
Then, using the definition of
$\varepsilon D_r(\theta )$
, for sides of the form,
$(r_i\cos {(\theta )},r_i\sin {(\theta )})$
\begin{equation} \begin{aligned} \varepsilon D_r(\theta )&=\sqrt {\left (1-r_i(\theta )\right )^2\cos ^2{(\theta )}+\left (1-r_i(\theta )\right )^2\sin ^2{(\theta )}}=-(1-r_i(\theta )), \end{aligned} \end{equation}
with (C2), we readily find
$\varepsilon D_r(\theta )$
for
$i=1,2,3,4$
. The correct sign of the square root in (C3) is minus, because if one refers to
$D_r(\theta )$
as to a vector, it points from the circle towards the sides of the hexagon.
Next, expressing
$\varepsilon D_r(\theta )$
for
$L_5$
and
$L_6$
(given in (C1)), where the corresponding ranges of
$\theta$
are
$I_5$
and
$I_6$
, respectively, we get that
\begin{equation} \varepsilon D_r(\theta )=\sqrt {\left (\cos {(\theta )}-\cot {(\theta )}\right )^2+\left (\sin {(\theta )}-1\right )^2}, \end{equation}
and
\begin{equation} \varepsilon D_r(\theta )=\sqrt {\left (\cos {(\theta )}+\cot {(\theta )}\right )^2+\left (\sin {(\theta )}+1\right )^2}. \end{equation}
To summarise our results in (C2), (C4) and (C5) we may conclude that
$\varepsilon D_r(\theta )$
in this case is given by
\begin{equation} \varepsilon D_r(\theta )=\left \{\begin{aligned} &-1+\frac {2}{\sqrt {3}\cos {(\theta )}-\sin {(\theta )}}, &&\theta \in I_1,\\ &-1-\frac {2}{\sqrt {3}\cos {(\theta )}+\sin {(\theta )}}, &&\theta \in I_2,\\ &-1+\frac {2}{\sqrt {3}\cos {(\theta )}+\sin {(\theta )}}, &&\theta \in I_3,\\ &-1-\frac {2}{\sqrt {3}\cos {(\theta )}-\sin {(\theta )}}, &&\theta \in I_4,\\ &\sqrt {\left (\cos {(\theta )}-\cot {(\theta )}\right )^2+\left (\sin {(\theta )}-1\right )^2}, &&\theta \in I_5,\\ &\sqrt {\left (\cos {(\theta )}+\cot {(\theta )}\right )^2+\left (\sin {(\theta )}+1\right )^2}, &&\theta \in I_6, \end{aligned}\right . \end{equation}
where the various intervals are given by
\begin{equation} \begin{aligned} &I_1=\left [\frac {5\pi }{3},2\pi \right )\!,\\ &I_2=\left [\pi ,\frac {4\pi }{3}\right )\!,\\ &I_3=\left [0,\frac {\pi }{3}\right )\!,\\ &I_4=\left [\frac {2\pi }{3},\pi \right )\!,\\ &I_5=\left [\frac {\pi }{3},\frac {2\pi }{3}\right )\!,\\ &I_6=\left [\frac {4\pi }{3},\frac {5\pi }{3}\right )\!. \end{aligned} \end{equation}
The function given in (C6) is shown in figure 5(a).
Next, we examine the expression for
$\varepsilon D_r(\theta )$
in the case of the pentagon, which is shown in the inset of figure 5(b). Note that the length of the pentagon’s side is
\begin{align} A=2\tan {(\pi /5)}. \end{align}
Hence, it is easy to verify that the sides of the pentagon, which we denote by
$\bar {L}_i$
,
$i=1,2,3,4,5$
(as indicated in the inset of figure 5
b) are given by
\begin{equation} \begin{aligned} &\bar {L}_i:\quad y_i=a_i x-\left (1+\frac {A}{2\tan {(\pi /10)}}\right )\!, &&i=1,2,\\ &\bar {L}_i:\quad y_i=a_i x+\frac {1}{\sin {(3\pi /10)}}, &&i=3,4,\\ &\bar {L}_5:\quad y_5=-1, &&x=\frac {y_5}{\tan {(\theta )}}=-\cot {(\theta )}, \end{aligned} \end{equation}
where
\begin{equation} a_1=\tan {\left (\frac {2\pi }{5}\right )}, \quad a_2=-\tan {\left (\frac {2\pi }{5}\right )}, \quad a_3=\tan {\left (\frac {4\pi }{5}\right )}, \quad a_4=-\tan {\left (\frac {4\pi }{5}\right )}. \end{equation}
Now, if we wish to express the sides
$\bar {L}_i$
,
$i=1,2,3,4$
, in polar coordinates, namely in the form
$(r_i\cos {(\theta )},r_i\sin {(\theta )})$
, where
$\theta$
is within the corresponding ranges, which are denoted by
$I_i$
,
$i=1,2,3,4$
, respectively, and are given in (C16), we need only to find
$r_i$
. However, according to (C9), it follows that
\begin{equation} a_i r_i\cos {(\theta )} -\left (1+\frac {A}{2\tan {(\pi /10)}}\right )=r_i\sin {(\theta )}, \quad i=1,2, \end{equation}
and
\begin{equation} a_i r_i\cos {(\theta )}+\frac {1}{\sin {(3\pi /10)}}=r_i\sin {(\theta )}, \quad i=3,4. \end{equation}
Hence,
\begin{equation} r_i(\theta )=\left \{\begin{aligned} &-\frac {1+\dfrac {A}{2\tan {(\pi /10)}}}{\sin {(\theta )}-\tan {(2\pi /5)}\cos {(\theta )}}, &&i=1,\\ &-\frac {1+\dfrac {A}{2\tan {(\pi /10)}}}{\sin {(\theta )}+\tan {(2\pi /5)}\cos {(\theta )}}, &&i=2,\\ &\frac {1}{\sin {(3\pi /10)}\left [\sin {(\theta )-\tan {(4\pi /5)}\cos {(\theta )}}\right ]}, &&i=3,\\ &\frac {1}{\sin {(3\pi /10)}\left [\sin {(\theta )+\tan {(4\pi /5)}\cos {(\theta )}}\right ]}, &&i=4. \end{aligned}\right . \end{equation}
Then, using the formula in (C3) with (C13), we readily find
$\varepsilon D_r(\theta )$
for
$i=1,2,3,4$
.
Finally, expressing
$\varepsilon D_r(\theta )$
for
$\bar {L}_5$
(given in (C9)), we get that
\begin{equation} \varepsilon D_r(\theta )=\sqrt {\left (\cos {(\theta )}+\cot {(\theta )}\right )^2+\left (\sin {(\theta )}+1\right )^2}. \end{equation}
To summarise our results in (C13) and (C14), we may conclude that
$\varepsilon D_r(\theta )$
in this case is given by
\begin{equation} \varepsilon D_r(\theta )=\left \{\begin{aligned} &-1-\frac {1+\dfrac {A}{2\tan {(\pi /10)}}}{\sin {(\theta )}-\tan {(2\pi /5)}\cos {(\theta )}}, && \theta \in I_1, \\ &-1-\frac {1+\dfrac {A}{2\tan {(\pi /10)}}}{\sin {(\theta )}+\tan {(2\pi /5)}\cos {(\theta )}}, && \theta \in I_2,\\ &-1+\frac {1}{\sin {(3\pi /10)}\left [\sin {(\theta )-\tan {(4\pi /5)}\cos {(\theta )}}\right ]}, && \theta \in I_3,\\ &-1+\frac {1}{\sin {(3\pi /10)}\left [\sin {(\theta )+\tan {(4\pi /5)}\cos {(\theta )}}\right ]}, &&\theta \in I_4,\\ &\sqrt {\left (\cos {(\theta )}+\cot {(\theta )}\right )^2+\left (\sin {(\theta )}+1\right )^2}, && \theta \in I_5, \end{aligned}\right . \end{equation}
where the various intervals are given by
\begin{equation} \begin{aligned} &I_1=\left [0,\frac {\pi }{10}\right ]\cup \left [\frac {17\pi }{10},2\pi \right )\!,\\ &I_2=\left (\frac {9\pi }{10},\frac {13\pi }{10}\right ]\!,\\ &I_3=\left (\frac {\pi }{10},\frac {\pi }{2}\right ]\!,\\ &I_4=\left (\frac {\pi }{2},\frac {9\pi }{10}\right ]\!,\\ &I_5=\left (\frac {13\pi }{10},\frac {17\pi }{10}\right )\!. \end{aligned} \end{equation}
The function given in (C15) is shown in figure 5(b).
Now consider the expression for
$\varepsilon D_r(\theta )$
in the case of the square, which is shown in the inset of figure 5(c). It is easy to verify that the sides of the square, which we denote by
$\tilde {L}_i$
,
$i=1,2,3,4$
(as indicated in the inset of figure 5
c), are given by
\begin{equation} \begin{aligned} &\tilde {L}_1: \quad x=1, \quad &&y=x\tan {(\theta )}=\tan {(\theta )},\\ &\tilde {L}_2: \quad x=-1, \quad &&y=x\tan {(\theta )}=-\tan {(\theta )},\\ &\tilde {L}_3: \quad y=1, \quad &&x=\frac {y}{\tan {(\theta )}}=\cot {(\theta )},\\ &\tilde {L}_4: \quad y=-1, \quad &&x=\frac {y}{\tan {(\theta )}}=-\cot {(\theta )}. \end{aligned} \end{equation}
Hence, the distance
$\varepsilon D_r(\theta )$
between the circle and
$\tilde {L}_1$
is given by
\begin{equation} \varepsilon D_r(\theta )=\sqrt {\left (\cos {(\theta )}-1\right )^2+\left (\sin {(\theta )}-\tan {(\theta )}\right )^2}, \end{equation}
where the corresponding range of
$\theta$
is given in
$I_1$
, see (C20).
The rest of the sides of the square are treated similarly, thus overall resulting that for
$\theta \in [0,2\pi )$
,
$\varepsilon D_r(\theta )$
is given by
\begin{equation} \varepsilon D_r(\theta )=\left \{\begin{aligned} &\sqrt {\left (\cos {(\theta )}-1\right )^2+\left (\sin {(\theta )}-\tan {(\theta )}\right )^2}, &&\theta \in I_1,\\ &\sqrt {\left (\cos {(\theta )}+1\right )^2+\left (\sin {(\theta )}+\tan {(\theta )}\right )^2}, &&\theta \in I_2,\\ &\sqrt {\left (\cos {(\theta )}-\cot {(\theta )}\right )^2+\left (\sin {(\theta )}-1\right )^2}, &&\theta \in I_3,\\ &\sqrt {\left (\cos {(\theta )}+\cot {(\theta )}\right )^2+\left (\sin {(\theta )}+1\right )^2}, &&\theta \in I_4,\\ \end{aligned}\right . \end{equation}
where
\begin{equation} \begin{aligned} &I_1=[0,0.25\pi )\cup (1.75\pi ,2\pi ),\\ &I_2=[0.75\pi ,1.25\pi ),\\ &I_3=[0.25\pi ,0.75\pi ),\\ &I_4=[1.25\pi ,1.75\pi ]. \end{aligned} \end{equation}
The function given in (C19) is shown in figure 5(c).
Finally, we derive the expression for
$\varepsilon D_r(\theta )$
in the case of the triangle, which is shown in the inset of figure 5(d). It is easy to verify that the sides of the triangle, which we denote by
$\hat {L}_i$
,
$i=1,2,3$
(as indicated in the inset of figure 5
d), are given by
\begin{equation} \begin{aligned} &\hat {L}_i:\quad y_i=a_i x+\frac {1}{\sin {(\pi /6)}}, &&i=1,2,\\ &\hat {L}_3:\quad y_3=-1, &&x=\frac {y_3}{\tan {(\theta )}}=-\cot {(\theta )}. \end{aligned} \end{equation}
where
\begin{align} a_1=\tan {(2\pi /3)}=-\sqrt {3} \quad \text{and} \quad a_2=\tan {(\pi /3)}=\sqrt {3}. \end{align}
Now, if we wish to express the sides
$\hat {L}_i$
,
$i=1,2$
, in polar coordinates, namely in the form
$(r_i\cos {(\theta )},r_i\sin {(\theta )})$
, where
$\theta$
is within the corresponding ranges, which are denoted by
$I_i$
,
$i=1,2$
, respectively, and are given in (C27), we need only to find
$r_i$
. However, according to (C21), it follows that
\begin{align} a_i r_i \cos {(\theta )}+\frac {1}{\sin {(\pi /6)}}=r_i \sin {(\theta )}, \quad i=1,2. \end{align}
Hence, we may conclude that
\begin{equation} \begin{aligned} r_i(\theta )&=\frac {1}{\sin {(\pi /6)}\left (\sin {(\theta )}-a_i\cos {(\theta )}\right )}\\ &=\frac {1}{2\big(\sin {(\theta )}\pm \sqrt {3}\cos {(\theta )}\big)}, \quad &&i=1,2, \end{aligned} \end{equation}
where
$\sqrt {3}$
corresponds to
$i=1$
and
$-\sqrt {3}$
corresponds to
$i=2$
. Then, using the formula in (C3) with (C24), we readily find
$\varepsilon D_r(\theta )$
for
$i=1,2$
.
Finally, expressing
$\varepsilon D_r(\theta )$
for
$\hat {L}_3$
(given in (C21)), where the corresponding range of
$\theta$
is
$I_3$
(given in (C27)), we get that
\begin{equation} \varepsilon D_r(\theta )=\sqrt {\left (\cos {(\theta )}+\cot {(\theta )}\right )^2+\left (\sin {(\theta )}+1\right )^2}. \end{equation}
To summarise our results in (C24) and (C25), we may conclude that
$\varepsilon D_r(\theta )$
in this case is given by
\begin{equation} \varepsilon D_r(\theta )=\left \{\begin{aligned} &-1+\frac {1}{2\big(\sqrt {3} \cos (\theta )+\sin (\theta )\big)}, && \theta \in I_1, \\ &-1-\frac {1}{2\big(\sqrt {3} \cos (\theta )-\sin (\theta )\big)}, && \theta \in I_2,\\ &\sqrt {\left (\cos {(\theta )}+\cot {(\theta )}\right )^2+\left (\sin {(\theta )}+1\right )^2}, && \theta \in I_3, \end{aligned}\right . \end{equation}
where the various intervals are given by
\begin{equation} \begin{aligned} &I_1=\left [0,\frac {\pi }{2}\right ]\cup \left [\frac {11\pi }{6},2\pi \right )\!,\\ &I_2=\left (\frac {\pi }{2},\frac {7\pi }{6}\right ]\!,\\ &I_3=\left (\frac {7\pi }{6},\frac {11\pi }{6}\right )\!. \end{aligned} \end{equation}
Appendix D. Estimation of the error in contact inclination angle at the particle boundary due to the liquid height uncertainty
The inclination angle of the meniscus near the particle boundary is denoted by
$\alpha$
and may be found according to the formula
\begin{equation} \alpha =\text{arccot}\left (\frac {\partial h}{\partial r}\Bigl |_{r=r_{\kern-1pt o}+d_r(\theta )}\right )=\text{arccot}\left (\frac {h_c}{r_0}\,\frac {\partial H}{\partial R}\Bigl |_{R=1+\varepsilon D_r(\theta )}\right )\!, \end{equation}
where we used the scaling parameters which were defined in (2.3). Now, from the expansion in (2.6), it easily follows that
\begin{equation} \frac {\partial H}{\partial R}=\frac {\partial H_0}{\partial R}+\varepsilon \frac {\partial H_1}{\partial R}+O(\varepsilon ^2). \end{equation}
According to the solution for the leading-order problem, which was given in (2.8)–(2.9), and the formula for the derivative of the Bessel functions given in (A7), it follows that
\begin{equation} \begin{aligned} \frac {\partial H_0}{\partial R}=&-\frac {A_0\sqrt {{B}_{{o}}}}{4}\big(K_{-1}\big(\sqrt {{B}_{{o}}}R\big)+K_{1}\big(\sqrt {{B}_{{o}}}R\big)\big)\\ &-\frac {\sqrt {{B}_{{o}}}}{2}\sum _{n=1}^{\infty }\left [A_n \cos {(n\theta )}+B_n\sin {(n\theta )}\right ]\big(K_{n+1}\big(\sqrt {{B}_{{o}}}R\big)+K_{n-1}\big(\sqrt {{B}_{{o}}}R\big)\big). \end{aligned} \end{equation}
Note that, since, according to (2.11), the expression for
$H_1(R,\theta )$
is the same as for
$H_0(R,\theta )$
, except that
$A_0$
,
$A_n$
and
$B_n$
,
$n=1,2,3,\ldots$
, are replaced with
$\tilde {A}_0$
,
$\tilde {A}_n$
and
$\tilde {B}_n$
, respectively, also the expression for
$\partial H_1/\partial R$
is the same as for
$\partial H_0/\partial R$
, except that
$A_0$
,
$A_n$
and
$B_n$
,
$n=1,2,3,\ldots$
, are replaced with
$\tilde {A}_0$
,
$\tilde {A}_n$
and
$\tilde {B}_n$
, respectively.
In order to assess the effect of uncertainty of the water level in the experiments, which is
$\pm$
0.05 mm, on the inclination angle
$\alpha$
, we calculated
$\alpha$
by using the equations in (D1)–(D3) and the derivative
$\partial H_1/\partial R$
for 3 representative points, which we denote by
$P_1$
,
$P_2$
and
$P_3$
(see figure 10) located on the boundaries of the square and hexagonal-shaped particles. One point was chosen close to the centre of the particles’ edge, one is chosen closer to the corner and one was randomly chosen in between them. This calculation was repeated for three cases: (i) the centreline of the particle is the same as the liquid–gas interface height at
$r\to \infty$
, (ii) the centreline of the particle is above the liquid–gas interface height at
$r\to \infty$
by 0.05 mm, and (iii) the centreline of the particle is below the liquid–gas interface height at
$r\to \infty$
by 0.05 mm. Our results are summarised in table 1.
Figure 10. The liquid–gas interface shape for a square and a hexagonal-shaped particles, for which we calculated the inclination angles in table 1, where the points
$P_1$
,
$P_2$
and
$P_3$
are indicated.
Table 1. The inclination angle
$\alpha$
, calculated according to (D1), on 3 points
$P_1$
,
$P_2$
and
$P_3$
in three cases: case (i), case (ii), and case (iii) where in case (i) the centerline of the particle is the same as of the liquid-gas interface at infinity, in case (ii) the centerline is 0.05 mm above the liquid-gas interface at infinity, and in case (iii) the centerline is 0.05 mm below the liquid-gas interface at infinity.
According to the representative results shown in table 1, it can be seen that the uncertainty of
$\pm 0.05$
mm in the measurement of the liquid–gas interface at far field results in the error of the order of
$\pm 1\,\%$
in the inclination angle of the meniscus.
Appendix E. The second-order expansion
Here, we derive an expression for the second-order term in the expansion of
$H(R,\theta )$
in asymptotic series. Similarly to (2.5), we expand
$(H_0+\varepsilon H_1+\varepsilon ^2 H_2)\vert _{R=1+\varepsilon D_r(\theta )}$
into a Taylor series in
$\varepsilon$
around
$\varepsilon =0$
, as follows:
\begin{align} \big(H_0+\varepsilon H_1+\varepsilon ^2 H_2\big)\big \vert _{{R} =1+\varepsilon D_r({\theta })} &=H_0\big \vert _{R=1}+ \varepsilon \left (\!H_1+D_r(\theta )\frac {\partial H_0}{\partial R}\right )\bigg \vert _{R=1} \nonumber \\ & \quad +\varepsilon ^2 \!\left (\!H_2+D_r(\theta)\frac {\partial H_1}{\partial R}+0.5D^2_r(\theta)\frac {\partial ^2 H_0}{\partial R^2}\!\right )\!\bigg \vert _{R=1} \nonumber \\ & \quad + O(\varepsilon ^3). \end{align}
Next, similarly to the discussion for
$H_1(R,\theta )$
(see (2.10)), we get at the second,
$O(\varepsilon ^2)$
, order the following problem:
\begin{equation} \begin{aligned} &R^2 \frac {\partial ^2 H_2}{\partial R^2}+R \frac {\partial H_2}{\partial R}+\frac {\partial ^2 H_2}{\partial \theta ^2}-{B}_{{o}} R^2 H_2=0,\\ &H_2\big \vert _{R=1}=-\left (D_r(\theta )\frac {\partial H_1}{\partial R}+0.5D^2_r(\theta )\frac {\partial ^2 H_0}{\partial R^2 }\right )\bigg \vert _{R=1},\\ &H_2\big \vert _{R\to \infty }=0. \end{aligned} \end{equation}
Note that, since the equation for
$H_2$
as given in (E2) is exactly of the same form as the equation for
$H_0$
and
$H_1$
, by using the same procedure as described above, we get the same general solution for
$H_2$
as for
$H_0$
and
$H_1$
, see (2.8) and (2.11), namely
\begin{equation} H_2(R,\theta )=\frac {K_0\big(\sqrt {{B}_{{o}}}R\big)\hat {A}_0}{2}+\sum _{n=1}^{\infty }K_n\big(\sqrt {{B}_{{o}}}R\big)\left [\hat {A}_n \cos {(n\theta )}+\hat {B}_n\sin {(n\theta )}\right ]\!, \end{equation}
where the constants
$\hat {A}_n$
and
$\hat {B}_n$
are determined by using the boundary conditions in (E2).
In order to determine the constants
$\hat {A}_n$
and
$\hat {B}_n$
explicitly, analogously to the coefficients
$\tilde {A}_n$
and
$\tilde {B}_n$
appearing in the Fourier series of
$H_1(R,\theta )$
, see (2.12), we get now that
\begin{equation} \begin{aligned} &\hat {A}_n=\frac {\sqrt {{B}_{{o}}}}{2\pi K_n\big(\sqrt {{B}_{{o}}}\big)}\int _{0}^{2\pi }[\tilde {G}(\theta )-\hat {G}(\theta )]\cos {(n\theta )}{\rm d}\theta , \quad n=0,1,2,\ldots ,\\ &\hat {B}_n=\frac {\sqrt {{B}_{{o}}}}{2\pi K_n\big(\sqrt {{B}_{{o}}}\big)}\int _{0}^{2\pi }[\tilde {G}(\theta )-\hat {G}(\theta )]\sin {(n\theta )}{\rm d}\theta , \quad n=1,2,\ldots , \end{aligned} \end{equation}
where for brevity we use the following notation:
\begin{equation} \begin{aligned} \tilde {G}(\theta )&=D_r(\theta )\left \{\frac {\left (K_{-1}\big(\sqrt {{B}_{{o}}}\big)+K_{1}\big(\sqrt {{B}_{{o}}}\big)\right )\tilde {A}_0}{2}\right . \\ &\quad\left . +\,\sum _{m=1}^{\infty }\big(K_{m-1}\big(\sqrt {{B}_{{o}}}\big)+K_{m+1}\big(\sqrt {{B}_{{o}}}\big)\big)\bigl [\tilde {A}_m \cos {(m\theta )}+\tilde {B}_m\sin {(m\theta )}\bigr ]\!\right \}\!, \end{aligned} \end{equation}
and
\begin{equation} \begin{aligned} \hat {G}(\theta )&=\frac {D^2_r(\theta )\sqrt {{B}_{{o}}}}{4}\left \{\frac {\left (K_{-2}\big(\sqrt {{B}_{{o}}}\big)+2K_{0}\big(\sqrt {{B}_{{o}}}\big)+K_{2}\big(\sqrt {{B}_{{o}}}\big)\right )A_0}{2}\right . \\ &\quad +\sum _{m=1}^{\infty }\big(K_{m-2}\big(\sqrt {{B}_{{o}}}\big)+2K_{m}\big(\sqrt {{B}_{{o}}}\big)+K_{m+2}\big(\sqrt {{B}_{{o}}}\big)\big)\\ &\quad \left . \times \left [A_m \cos {(m\theta )}+B_m\sin {(m\theta )}\right ]\vphantom{\frac {\left (K_{-2}\big(\sqrt {{B}_{{o}}}\big)+2K_{0}\big(\sqrt {{B}_{{o}}}\big)+K_{2}\big(\sqrt {{B}_{{o}}}\big)\right )A_0}{2}}\right \}\!. \end{aligned} \end{equation}
Note that, in the derivation of the expression in (E6), we used (A7), according to which it easily follows that
\begin{equation} \frac {\partial ^2 K_n\big(\sqrt {{B}_{{o}}}R\big)}{\partial R^2}=\frac {{B}_{{o}}}{4}\big(K_{n-2}\big(\sqrt {{B}_{{o}}}R\big)+2K_{n}\big(\sqrt {{B}_{{o}}}R\big)+K_{n+2}\big(\sqrt {{B}_{{o}}}R\big)\big). \end{equation}
Open access
