The paper considers two static problems from capillarity. The first one consists in the determination of the surface of a liquid in a capillary tube and the second in the computation of the shape of a sessile or pendent drop of liquid of a given volume, both configurations being considered in the gravitational field. From the mathematical point of view both problems are nonlinear two-point boundary value problems which include in their formulations additional geometrical unknowns, i.e., the so-called free boundary value problems. Computationally they are laborious problems since they involve non-normal Jacobian matrices. The resulting numerical difficulties are considerably reduced by treating these problems by a collocation method in conjunction with the continuation with respect to the parameters. The method is implemented by the MATLAB code bvp4c. The numerical evaluations of the free surfaces are in reasonable accordance with asymptotic estimations.
Authors
Călin-Ioan Gheorghiu
(Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy)
Keywords
capillarity; free boundary value problem; liquid gas interface; sessile; pendent drop; collocation; matlab; bvp4c.
References
See the expanding block below.
Cite this paper as
C.I. Gheorghiu, Accurate numerical solutions to stationary free surface problems from capillarity, Phys. Part. Nuclei Lett., 5 (2008), 286-289.
doi: 10.1134/s154747710803031x
[1] Finn R. Capillary Free Surfaces. Preprint No. 7/1981. Bonn: Inst. féur Angewante Mathematik Univ.,1981. Google Scholar
2. Gheorghiu C.-I., Mathematica. Studia Univ. Babes-Bolyai. 2005. V. L. P. 61.
3. Kierzenka J., Shampine L. F. // ACM Trans. Math. Softw. 2001. V. 27. P. 299.
4. Miele A., Aggarwal A. K., Tietze J. L. // J. Comp. Phys. 1974. V. 15. P. 117.
5. Trefethen L. N., Bau III D. Numerical Linear Algebra. SIAM, 1997.
Paper (preprint) in HTML form
S154747710803031X
Accurate Numerical Solutions to Stationary Free Surface Problems from Capillarity ^(1){ }^{\mathbf{1}}
C.-I. Gheorghiu"T. Popoviciu" Institute of Numerical Analysis, Cluj-Napoca, Romania
Abstract
The paper considers two static problems from capillarity. The first one consists in the determination of the surface of a liquid in a capillary tube and the second in the computation of the shape of a sessile or pendent drop of liquid of a given volume, both configurations being considered in the gravitational field. From a mathematical point of view, both problems are nonlinear two-point boundary value problems which include in their formulations additional geometrical unknowns, i.e., the so-called free boundary value problems. Computational, they are laborious problems owing to the fact that they involve nonnormal Jacobian matrices. The resulting numerical difficulties are considerably reduced by treating these problems by a collocation method in conjunction with the continuation with respect to the parameters. The method is implemented by the MATLAB code bvp4cb v p 4 c. The numerical evaluations of the free surfaces are in reasonable accordance with asymptotic estimations.
INTRODUCTION
In many cases, problems in viscous fluid mechanics can be described by boundary value problems for partial differential equations. In some cases, using the geometrical symmetries of the configurations, such problems can be reduced to two-point boundary value problems. When there are additional geometrical unknowns, we speak of free boundary value problems (FB for short). More specific, these geometrical unknowns can be some values of the unknowns on subsets of boundary points or even the coordinates of the boundary points themselves. In this paper we solve numerically two classes of static capillary FBs. These problems are defined as FBs on which we must take into account capillary forces, i.e., forces due to the intermolecular attractions, which have a nonvanishing resultant at the boundary of the liquid.
The first one is to determine the liquid-air interface (surface) whose mean curvature is proportional to its height above the reference plane, and which makes with a container bounding wall a prescribed constant angle. The second one is to find the shape of a sessile or pendent liquid drop of a given volume, which rests in the gravitational field. The actual computations were carried out for symmetric cases. The equations of these problems, as well as the analytic properties of the solutions (existence, uniqueness, and asymptotic behavior for Bond number approaching 0), were established in the monograph of Finn [1]. Both problems lead to nonlinear second order two-point free boundary value problems, i.e., the values of the solution at some bound-
ary points are unknown, but some values of the derivatives are supplied in order to close the problem. Their solutions depend on the nondimensional Bond number (Bo for short) which is a measure of the "size" of the configuration, i.e., it is the ratio of gravitational forces to the capillary (superficial) forces. The main computational difficulty with these problems consists in the fact that they involve nonnormal Jacobian matrices. A test problem with a Jacobian matrix in close form is considered in order to underline these aspects.
The problems are solved numerically making use of a collocation type method, based on cubic interpolation and introduced in the paper of Kierzenka and Shampine [3]. It is implemented by the MATLAB code bvp4cb v p 4 c and provides a smooth solution continuous along with its first derivative. In contrast with well-known shooting, this solution is approximated on the whole interval and the boundary conditions are taken into account at all times. As the resulting nonlinear system of algebraic equations is solved iteratively, the question of coming up with a sufficiently accurate initial guess of unknowns is particularly discussed.
1. A TEST PROBLEM
In their paper [4], Miele, Aggarwal, and Tietze consider the following two-point boundary value problem introduced by Troesch, namely,
is characterized by the eigenvalues lambda=+-ksqrt(cosh(kx))\lambda= \pm k \sqrt{\cosh (k x)}, which at the endpoints become lambda(0)=+-k\lambda(0)= \pm k and lambda(1)=+-ksqrt(cosh(k))\lambda(1)= \pm k \sqrt{\cosh (k)}.
For relatively low values of kk, i.e., 0 < k <= 40<k \leq 4, the above eigenvalues remain relatively small and the problem can be solved employing the usual finite difference or the shooting methods. On the other hand, for relatively large values of kk, i.e., 5 <= k <= 105 \leq k \leq 10, the above quoted authors make use of a special multipoint boundaryvalue technique in order to solve the problem. They observe that the problem is characterized by the first integral
and use that as a "tool" in order to check the accuracy of the numerical solution. In our numerical experiments performed on a Pentium IV, 3.0 G , PC, the constancy of zz is verified to a maximal error equal to 4.082997 e-064.082997 e-06 for k=5k=5, to an error which equals 7.651421e -03 for k=10k=10, and to an error equal to 6.602193 e+006.602193 e+00 corresponding to k=15k=15. Our numerical results are with one significant figure worse than those reported in [4, Table 2, p. 126], for k=5k=5 and 10, respectively. For kk larger than 16 , the accuracy in the conservation of zz is lost rapidly. In spite of this fact, even at these values of kk, our numerical solutions behave correctly in the left neighborhood of t=1t=1. They show a boundary layer type behavior in this region.
It is worth noting at this point that we used essentially the technique of continuation with respect to parameter kk in order to handle the problem for k >= 15k \geq 15. Without this technique for such values of kk, the collocation equations become unsolvable due to singular Jacobian.
Remark 1. In fact the above Jacobian matrix JJ is nonnormal, i.e., J^(**)J-JJ^(**)!=0J^{*} J-J J^{*} \neq 0, where *** stands for the conjugate transpose of a matrix (see Trefethen [5, p. 187]). The pseudospectrum of the matrix JJ, at t=1t=1, i.e., the spectrum of the randomly perturbed matrix with an arbitrary small quantity, is depicted in Fig. 1. In our previous paper [2], we introduced a scalar measure of nonnormality. It reads H(J):=sqrt(epsi(J^(**)J-JJ^(**)))//epsi(J)H(J):=\sqrt{\varepsilon\left(J^{*} J-J J^{*}\right)} / \varepsilon(J) and we noticed that 0 <= H(J) <= root(4)(2),00 \leq H(J) \leq \sqrt[4]{2}, 0 being attained just in case of normal (symmetric) matrices. On this scale our matrix JJ is situated extremely close to the upper bound. This scalar measure as well as the pseudospectrum lead us to the conclusion that the matrix JJ is a genuine nonnormal one.
Fig. 1. The pseudospectrum of matrix JJ for k=4k=4.
2. THE LIQUID-GAS INTERFACE IN A CAPILLARY TUBE
In this section we consider static capillary liquidgas interfaces of liquids partly filling a vessel. More exactly we examine the behavior of capillary surfaces in a tube with circular section of radius R=aR=a. We assume the following system of two coupled nonlinear differential equations and boundary conditions (see [1, p. 11])
{:(3){[div(Tu)=Bou","quad" in "Omega],[vTu=beta","quad" on "quad Sigma:=del Omega","]:}:}\left\{\begin{array}{l}
\operatorname{div}(\mathrm{Tu})=\mathrm{Bo} u, \quad \text { in } \Omega \tag{3}\\
v \mathrm{Tu}=\beta, \quad \text { on } \quad \Sigma:=\partial \Omega,
\end{array}\right.
where, u(x,y)u(x, y) denotes the height over the wetted section Omega subR^(2),T(u):=(grad u)/(sqrt(1+|grad u|^(2))),beta\Omega \subset \mathbb{R}^{2}, T(u):=\frac{\nabla u}{\sqrt{1+|\nabla u|^{2}}}, \beta is known, and vv is the outward normal to Sigma\Sigma.
Due to the symmetry with respect to the axis of the tube, on the half of the axial section, i.e., 0 <= r <= 1,r:=R//a0 \leq r \leq 1, r:= R / a, the problem reads
{:(4){[((ru_(r))/(sqrt(1+u_(r)^(2))))_(r)=" Bo "ru","quad0 < r < 1],[u(0)=p_(1)","quadu_(r)(0)=0],[u(1)=p_(2)","quad cos(delta)=beta","quad r=1]:}:}\left\{\begin{array}{l}
\left(\frac{r u_{r}}{\sqrt{1+u_{r}^{2}}}\right)_{r}=\text { Bo } r u, \quad 0<r<1 \tag{4}\\
u(0)=p_{1}, \quad u_{r}(0)=0 \\
u(1)=p_{2}, \quad \cos (\delta)=\beta, \quad r=1
\end{array}\right.
where Bo :=rho ga^(2)//sigma,rho:=\rho g a^{2} / \sigma, \rho and sigma\sigma being respectively the density and the surface tension of the liquid, and gg is the usual gravitational acceleration. The parameters p_(1)p_{1} and p_(2)p_{2} are unknown and signify the heights of the liquid on the axis of symmetry and on the wall of the tube, respectively. The boundary condition u_(r)(0)=0u_{r}(0)=0 appears
Fig. 2. The liquid-gas interface.
due to the symmetry and delta\delta is the angle between the liq-uid-gas interface u(r)u(r) and the wall of the tube (see Fig. 2). The problem is effectively integrated as a first order differential system, i.e.,
The second differential equation in (5) is singular at r=0r=0, but the singular coefficient arises from the coordinate system and we expect a smooth solution in the neighborhood of this point. We deal with this singularity using the technique suggested in Example 6 from [3, pp. 18-19]. It means that we let r longrightarrow0r \longrightarrow 0 in this equation and obtain the value that must be used for evaluating the differential equation when r=0r=0.
The existence of symmetric solutions with prescribed angle delta\delta has been proved and the uniqueness theorem shows that the symmetric solutions are the only ones (see [1] and the paper quoted there). For Bolongrightarrow0\mathrm{Bo} \longrightarrow 0, in the same monograph, the Laplace formula, which estimates p_(1)p_{1}, and estimations for the parameter p_(2)p_{2} are also available. These estimations were fairly useful in guessing the initial shape of the unknown u(r)u(r).
For Bo=0.075\mathrm{Bo}=0.075 and Bo=10\mathrm{Bo}=10, the dependence of u(r)u(r) on the angle delta\delta is depicted in Fig. 2.
Remark 2. From the computational point of view, it is suggestive to notice that the code bvp4c works such that the maximal residual belongs to the range 4.563 e-4.563 e- 007 and 7.354 e-0047.354 e-004.
Fig. 3. The shape of the sessile drops.
3. THE SESSILE AND PENDENT DROPS
We consider a liquid drop of volume VV resting on a horizontal plane (surface) in a vertical gravity field and restrict our attention to the case of the contact angle delta\delta such that 0 < delta <= pi//20<\delta \leq \pi / 2 (see Fig. 3). The case delta > pi//2\delta>\pi / 2 reduces to that one under the transformation u longrightarrow-uu \longrightarrow-u. Supposing that the drop is symmetric, the equations for the surface of the drop are
where lambda\lambda is a Lagrangian multiplier corresponding to the volume constraint. It is shown (see [1, pp. 59-60]) that, if delta=\delta= const, then all equilibrium surfaces are symmetric and there is a one-to-one correspondence between capillary surfaces furnished by (4) and sessile drops implied by (6). Our numerical results are reported in Fig. 3.
Remark 3. An upper bound for the existence of the solution with respect to Bo can be obtained. The principle of virtual work applied to the energy functional of the drop of volume VV yields a slightly modified form of (3), namely,
{:(7){[div(Tu)=Bou+lambdaquad" in "Omega],[vTu=betaquad" on "Sigma:=del Omega","]:}:}\left\{\begin{array}{l}
\operatorname{div}(\mathrm{Tu})=\mathrm{Bo} u+\lambda \quad \text { in } \Omega \tag{7}\\
v \mathrm{Tu}=\beta \quad \text { on } \Sigma:=\partial \Omega,
\end{array}\right.
where lambda\lambda is associated with the volume restriction. The volume VV of the drop can be expressed explicitly in terms of the boundary data and this property holds for general section Omega\Omega. We integrate the first equation in (3) and introduce the boundary condition. Thus, the volume is expressed as
{:(8)V=(1)/(Bo)(Sigma cos delta-lambda Omega)",":}\begin{equation*}
V=\frac{1}{\mathrm{Bo}}(\Sigma \cos \delta-\lambda \Omega), \tag{8}
\end{equation*}
where the symbols Sigma\Sigma and Omega\Omega denote, respectively, the measures of these sets. For symmetric solutions, (8) becomes
As the Lagrangian multiplier lambda\lambda is delivered by bvp4cb v p 4 c, this upper bound for Bo can be verified numerically. The results of our numerical experiments show that the quantity (2-lambda)/(2int_(0)^(1)u(r)dr)\frac{2-\lambda}{2 \int_{0}^{1} u(r) d r} overestimates the Bo number for
Bo <= 10\leq 10 and unfortunately underestimates that for Bo > 10. It remains an open problem to explain why, for large values of Bo, this inequality is not fulfilled.
CONCLUSIONS
The numerical experiments reported in this paper, as well as some others which are in progress and refer to FB problems from laminar boundary layer theory, entitle us to think that the main difficulty with this type of problem consists in the high nonnormality of their Jacobian matrices. The collocation method implemented by bvp4cb v p 4 c, due to its global character, overcomes to some extent these difficulties. As a general remark, we have to observe that the problems corresponding to microgravity conditions, i.e., Bond number closed to zero, are easier to be solved and the results are much more accurate than those obtained when this number is increased.
REFERENCES
R. Finn, "Capillary Free Surfaces," Preprint 7/1981 (Institut für Angewante Mathematik Univ. Bonn, 1981)
C.-I. Gheorghiu, Studia Univ. "Babes-Bolyai," Mathematica L,61\mathbf{L}, 61 (2005).
J. Kierzenka and L. F. Shampine, ACM Trans. Math. Softw. 27, 299 (2001).
A. Miele, A. K. Aggarwal and J. L. Tietze, J. Comput. Phys. 15, 117 (1974).
L. N. Trefethen and D. Bau III, Numerical Linear Algebra (SIAM, 1997).
^(1){ }^{1} The text was submitted by the author in English.