A Shape Newton Scheme for Deforming Shells with Application to Capillary Bridges

We assume that $\Gamma=\partial\Omega$ is the surface of some volume of fluid $\Omega$. This capillary surface $\Gamma$ can be disjointly decomposed into $\Gamma_{\text{LA}}$, the liquid air interface, and $\Gamma_{\text{LS}_{i}}$, the interface between the liquid and solid body $i$. The resulting interface lines, or solid-gas-liquid triple phase lines are subsumed as boundaries $\partial\Gamma$. The geometrical setting is shown in Figure 1.

In particular, we use $\mathbf{n}$ to denote the outer normal of the volume $\Omega$ and $\bm{\mu}$ to denote co-normals of $\Gamma$, i.e., unit vectors that are orthogonal to $\partial\Gamma$ and to $\mathbf{n}$. In order to reduce the computational workload and to avoid additional rank deficits of the Hessian, we are in particular interested in problem formulations that are tailor made for shell meshes, i.e., grids of topological dimension two and geometrical dimension three, and which do not require volume integrals. Hence, the surface area of each interface and the total volume of fluid are then given by the following surface integrals

Furthermore, the gravitational force acting on the liquid is given by

where $\mathbf{g}$ is the gravitational vector and $\tilde{\mathbf{g}}(s):=\frac{1}{2}(g_{1}s_{1}^{2},g_{2}s_{2}^{2},g_{3}s_{3}^{2})^{T}$.

The formation of the capillary surface follows the minimization of the total interface energy, leading to the problem

where, $\beta_{i}>0$ is the non-dimensional relative adhesion coefficient of solid $i$ and $b\left\|\mathbf{g}\right\|\geq 0$ is the Bond number, the non-dimensional gravitational influence. Finally, $\Gamma_{\text{S}_{i}}$ denotes the surface of solid body $i$.

The interface condition $\Gamma_{\text{LS}_{i}}\subset\Gamma_{{\text{S}}_{i}}$ is not directly numerically tractable and can also not readily be addressed via shape calculus. Suppose the discretization of the subdomain $\Gamma_{\text{LS}_{i}}$ consists of the vertices $s_{i,j}$, $j=1,...,k_{i}$. We then replace the analytic subset condition with the constraint of minimal distance for each vertex $s_{i,j}$, namely

where $\operatorname{dist}(.,\Gamma_{\text{S}_{i}}):\mathbb{R}^{3}\rightarrow\mathbb{R}$ is some signed distance measure from any point in $\mathbb{R}^{3}$ to $\Gamma_{\text{S}_{i}}$, the surface of the solid body $i$. We follow the convention of the distance being negative if the point $s$ is in the interior of $\Gamma_{\text{S}_{i}}$.

To find the surface $\Gamma$ solving (3), we employ techniques from shape optimization, in particular shape calculus. Following the approach of perturbation of identity, a deformed surface $\Gamma_{\varepsilon}$ is defined via

where $\mathbf{V}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ is a suitably smooth vector field and $s_{\varepsilon}:=s+\varepsilon\mathbf{V}(s)$. Following [11, 44], the Eulerian-Semi Derivative of a shape functional

is then given by

provided the shape is of sufficient regularity, with (5) requiring the least. Here, ${\operatorname{div}_{\Gamma}\,}$ is the divergence in the tangent space of $\Gamma$ and $\kappa:={\operatorname{div}_{\Gamma}\,}\mathbf{n}$ is the additive curvature, with $\mathbf{n}$ being the outer normal to $\Omega$ and $\bm{\mu}$ being the normal to $\partial\Gamma$, which also fulfills $\langle\bm{\mu},\mathbf{n}\rangle=0$ and is also called the co-normal. Finally, $df[\mathbf{V}](s):=\frac{d}{d\varepsilon}_{|_{\varepsilon=0}}f(\varepsilon,s_{\varepsilon})$ is the material derivative and $f^{\prime}[\mathbf{V}](s):=\frac{\partial}{\partial\varepsilon}_{|_{\varepsilon=0}}f(\varepsilon,s)$ is the local derivative. Due to the chain rule, the identity

holds, if both local and material derivative exist. Here, $\operatorname{D}$ is the classical spatial Jacobian. Hence, for the spatial coordinates $x$ themselves, one arrives at

It is worth mentioning that, contrary to local derivatives, the material derivative and spatial differentiation do not commute, meaning

For the outer normal, one has in particular [41, 44]

where $\operatorname{D}_{\Gamma}\mathbf{V}=\operatorname{D}\mathbf{V}-\operatorname{D}\mathbf{V}\mathbf{n}\mathbf{n}^{T}$ is the tangential Jacobian of $\mathbf{V}$ with respect to the spatial coordinate. Hence, (6) is essentially the chain rule relating the material derivative to the local derivative, where the tangential part of the convective part has been merged with the divergence. Furthermore, Equation$~{}\eqref{eq:SurfDerive3}$ is created by integration by parts on surfaces, where the curvature term arises because $\mathbf{V}$ is not tangent to $\Gamma$.

We revisit the necessary optimality conditions of (3) from [19] using shape calculus, in comparison to a differential geometric perspective used in the former. The Lagrangian of (3) is given by

We now discuss the shape derivative of the Lagrangian term by term. Because the integrand is a constant with respect to shape perturbations, the shape derivative of most terms in the Lagrangian can immediately be found to be

where $\delta_{ij}$ is the Kroenecker delta. The distance function $\operatorname{dist}_{\Gamma_{\text{LS}_{i}}}(s)$ is independent of the deformation parameter $\varepsilon$. Hence, the shape derivative vanishes and the material derivative is given by transport only, i.e.,

where $\operatorname{D}$ is the classical spatial Jacobian and $\mathbf{d}_{i}(s)$ is the spatial gradient of the distance field with respect to $\Gamma_{\text{LS}_{i}}$ at point $s$. Picking the respective curvature free expressions from (10) with the least regularity requirements on $\Gamma$, one arrives at an alternative variational formulation of (YLE), well-suited for discretization via triangular meshes. Indeed, the following section shows how the respective high-regularity expressions from (10) create (YLE).

Omitting the subset constraint, one arrives at

where $\partial\Gamma_{i}$ is the contact line between the liquid to air interface and the contact surface on solid body $i$. In the absence of gravity, i.e. $b=0$, the above expression can only vanish, if the liquid air interface $\Gamma_{\text{LA}}$ is a surface of constant mean curvature, i.e., $\kappa=-\lambda_{\text{Vol}}$, which is indeed the case where (YLE) is spatially constant. Furthermore,

has to hold on $\partial\Gamma_{i}$, which is not possible for arbitrary perturbation fields $\mathbf{V}$. If $\Gamma$ is feasible, then $\bm{\mu}_{i}$ is in the tangent space of $\Gamma_{\text{S}_{i}}$. Hence, if one restricts the motion to tangent directions only, i.e. $\mathbf{V}=\alpha(s)\bm{\mu}_{i}(s)$, then the necessary optimality condition becomes

Because the co-normals need to satisfy $\left\langle{\bm{\mu}_{i}},{\bm{\mu}_{\text{A}}}\right\rangle=\beta_{i}$ and have unit length, the critical surface must form a contact angle of $\cos\beta_{i}$ and stationarity is only possible for tangent motions.

Our desired computational approach is a Newton-iteration to find the roots of $d\mathcal{L}(\Omega,\lambda_{\text{vol}},\lambda_{i,j})[\mathbf{V},d\lambda_{\text{vol}}[\mathbf{V}],d\lambda_{i,j}[\mathbf{V}]]$. Within the optimization context, this is also called sequential quadratic programming (SQP), because during each update, the second order Taylor-expansion of $\mathcal{L}$ is minimized. For the problem at hand, this means we require the second order shape derivative of functionals of the types

the latter arising due to our reformulation of volume integrals to surfaces in (1) and (2).

Indeed, it is possible to transform $J_{2}$ back into a volume integral, which would, however, result in two problems: First, using the volume formulation of the Hessian over $\Omega$ necessitates using a tetrahedral mesh, resulting in considerable numerical overhead. Furthermore, the Hessian matrix post discretization would then suffer from an additional rank deficit for vertices in the interior. Second, using a boundary formulation of the volumetric problem description requires the computation of the total curvature $\kappa$, which is an edge concentrated quantity for the shell meshes consisting of planar triangles our code is using. Hence, optimization and discretization would not commute in this situation.

A repeated application of (5) leads to

The above expression is more involved than those readily available in the literature, as the consideration of surface integrals here requires the computation of $d({\operatorname{div}_{\Gamma}\,}\mathbf{V})[\mathbf{W}]$, the material derivative of the tangent divergence, whereas the volume situation is usually considered otherwise. Material derivatives commute with algebraic operations, such that

Due to (8), one arrives at

In the last step, the remaining full Jacobian $\operatorname{D}\mathbf{V}$ can be replaced with the tangential Jacobian, because

Hence, taking the trace results in

To achieve symmetry, we assume $d\mathbf{V}[\mathbf{W}]=0$ here and in the following. This makes the shape Hessian derived via repeated differentiation align with the expressions obtained by enforcing a vector space structure. It is worth mentioning that the latter part of (14) does not appear when volume integrals are considered. Summarizing the above, the shape Hessian of $S(\Gamma)$ from (1) is given by

The shape Hessian of $F(\Gamma)$ and $G(\Gamma)$ in (1) and (2) is more involved, as the second material derivative $d^{2}\mathbf{n}[\mathbf{V},\mathbf{W}]$ of the normal is needed. For the first Eulerian derivative, Equation (5) is applied to the $J_{2}$ case above, leading to

and the symmetric repeated differentiation second order shape derivative is given by

The second material derivative of the normal can be explicitly computed using (9),

where we have again used the independency $d\mathbf{V}[\mathbf{W}]=0$ in the last step. Thus, when inserting (17) into (16), one arrives at an explicit expression for the second order shape derivative $d^{2}J_{2}[\mathbf{V},\mathbf{W}]$ in a general setting, which can now easily be adapted to the second derivative of $F(\Gamma)$ and $G(\Gamma)$ in (1) and (2), using the respective material derivatives and the independency $d\mathbf{V}[\mathbf{W}]=0$, namely

for $F(\Gamma)$ and

for $G(\Gamma)$. Finally, the second derivative of the $\Gamma$-independent distance-function is readily given by

where $\mathbf{H}$ is the spatial Hessian matrix of the distance field.

To numerically find the shape of the capillary bridge, we use finite elements on a shell mesh $\Gamma_{h}$ of topological dimension two consisting of planar triangles $T$ in $\mathbb{R}^{3}$ provided by the FEniCS framework [27]. The classical continuous Galerkin finite element space of order $r$ of $d$-dimensional vectors is given by

where $\mathcal{P}_{r}$ is the space of polynomials of order $r$. To include the respective scalar constraints of volume and distance, we augment this space to $Q:=\mathcal{CG}_{1}^{3}(\Gamma)\times\mathbb{R}\times\mathbb{R}^{n_{i}n_{j}}$ to include the adjoint multipliers. The Newton update for iteratively finding the capillary bridge is then interpretable as a variational problem, namely to find $q[\mathbf{W}]:=(\mathbf{W},d\lambda_{\text{vol}}[\mathbf{W}],d\lambda_{i,j}[\mathbf{W}])\in Q$, such that

Post discretization, this variational problem creates the typical Karush-Kuhn-Tucker (KKT) system. Finally, the state is updated via

We use the discrete shape differentiability capabilities of FEniCS [18] to numerically validate the $\mathbf{V}$-$\mathbf{W}$ block in (18), confirming that indeed for our first and second order shape derivatives the “discretize-then-optimize” and “optimize-then-discretize” approach commutes irrespective of mesh size. Including the off-diagonal blocks in this confirmation is unfortunately not conveniently possible, as FEniCS does not provide a means to compute the $\mathbf{V}$-$d\lambda[\mathbf{W}]$ off-diagonals discretely.

A straightforward implementation of (18) results in the challenge of globalizing the convergence, i.e., guaranteeing a positive definite Hessian matrix away from the optimum. To address this, we have also implemented an approximative Newton scheme, where the $\mathbf{V}$, $\mathbf{W}$ block in (18) is replaced by the $H^{1}(\Gamma)$ scalar product

where $\gamma>0$ is some smoothing parameter. This approach is sometimes also called a Sobolev gradient descent and the resulting KKT-system is positive definite.

There are several reasons for rank deficiencies of the proper KKT system for shapes, also at the optimum. Using a full 3D deformation $\mathbf{V},\mathbf{W}\in\mathcal{CG}_{1}^{3}(\Gamma_{h})$ is indeed too rich. As can be seen from (6) and (7), tangential motions are in the kernel of (18). Furthermore, the second order shape derivatives (12), (15) and (16) only involve spatial derivatives of $\mathbf{V}$ and $\mathbf{W}$, such that translations are also within the kernel. As such, one typically restricts the deformation unknown to $\mathbf{V}=\alpha\mathbf{N}$, where $\alpha\in\mathcal{CG}_{1}^{1}(\Gamma_{h})$ and $\mathbf{N}$ is some vertex average of the cell normal $\mathbf{n}$. However, as discussed in Section 2.3, the contact angle requirement (11) demands a tangential motion of the vertices constituting $\Gamma_{\text{LS}_{i}}$ and in particular $\partial\Gamma_{i}$. As such, we restrict the motion fields

where $\alpha_{1}$, $\alpha_{2}$ are $\mathcal{CG}_{1}^{1}(\Gamma_{h})$ functions, likewise for $\mathbf{W}$. Thus, we admit a motion in normal direction only for the vertices, except for those forming the solid-gas-liquid triple phase lines. There, each vertex is allowed to move in their respective two dimensional subspace spanned locally by the two edge-averaged co-normals $\bm{\mu}_{\text{LA}}$ and $\bm{\mu}_{\text{LS}}$. It is worth mentioning that as a consequence, once feasibility with respect to the distance constraint is achieved, there is no longer any possible motion of the interior vertices of $\Gamma_{\text{LS}}$, except for the contact line $\partial\Gamma:=\Gamma_{\text{LA}}\cap\Gamma_{\text{LS}}$, which can still move in the plane spanned locally by the vertex co-normals. Hence, if one switches to the Newton-solver too early, the mesh quality and point spacing of $\Gamma_{\text{LS}}$ is going to degrade. At present, we counter this problem by remeshing, which is necessary as we aim at being able to deal with large scale deformations from the initial guess, anyway. Incorporating equal vertex spacing into the optimization goal to remove the tangential kernel of the Hessian for vertices in the interior of $\Gamma_{\text{LS}}$ is part of future work.

Depending on the geometrical setup of the solid bodies $\Gamma_{\text{S}_{i}}$, translations can also be in the kernel of (18). One such case would be a capillary bridge between two planes. Full rank in this situation is achieved by activating additional centroid constraints. Surface and volume centroid are given by

While using the surface centroid might seem to be the approach of choice for a code operating on the shell alone, we rather choose the volume centroid, as the expression can be considerably simplified by exploiting the volume constraint in this setting. Thus, we constrain the $i$-th component of the centroid via

where $\text{vol}(\Omega^{*})$ is the fixed target volume, $\mathbf{f}_{2}^{i}(s):=\frac{1}{2}\left(s_{j}^{2}\delta_{ij}\right)_{j=1}^{3}$ and $\delta_{ij}$ is the Kronecker-Delta. The respective first and second order Eulerian derivatives are again given by (15), (16) and (17), where

More details on these constraints can also be found in [42].

In order to validate our second order optimization scheme, we consider the case of a spherical body with zero distance from a plane, as this situation is very well understood theoretically. In the following, we base our evaluation on the data available in [29]. In particular, there is a closed form solution of the total force available as

where $\theta_{i}$ are the wetting angles at the sphere and plane. Furthermore, $\psi$ is the so called filling angle. The available data is given in non-dimensionalized form as the ratio of contact force $F_{t}$ to radius $R$ times surface tension $\sigma$. However, there is no closed form solution for the curvature $H$ as a coupled system of non-linear equations has to be solved. Tabularized values are provided in [29] depending on the angles $\theta_{1}$ and $\psi$ for the curvature. Because our code requires the desired contact angle and volume as input rather than the angle $\theta_{1}$, we use a volume constraint of $0.165$ as per the tabularized values in [29], corresponding to the angles $\theta_{1}=\theta_{2}=40^{\circ}$. Seeing that curvature is not readily available for a discrete surface of planar triangles, we use the adjoint of the volume constraint for comparisons, utilizing the correspondence

as described in Section 2.4 and validate the resulting effective pressure difference $\Delta p$ as the corresponding quantity. Instead of using the close-form representation of the forces (21) for spheres, we follow the general expression for forces [47] and compute the force acting on obstacle $\Gamma_{\text{S}_{i}}$ via numerical quadrature of the following integrals

with $\mathbf{F}_{t,i}:=\mathbf{F}_{p,i}+\mathbf{F}_{s,i}$ being the resulting force acting on obstacle $i$. The vector $\bm{\mu}$ is again the co-normal on $\partial\Gamma_{i}$, e.g. a vector orthogonal to $\partial\Gamma_{i}$ and tangential to $\Gamma_{\text{LA}}$. In particular the integrands in $\mathbf{F}_{p}$ and $\mathbf{F}_{s}$ are constants per cell or edge for our mesh consisting of planar triangles. As such, the above integrals can readily be computed by summing the contribution of each triangle and edge individually.

The shapes of the two obstacles, the sphere and plane, enters our code via closed form descriptions for the signed distances in (4), in particular

where $z=0$ is twice the distance between the objects, i.e., $z=0$ here. First and second order spatial derivatives of these functions, which are needed for the shape derivative and shape Hessian, can easily be computed and the non-differentiability of the square root at zero is irrelevant, as the centroid of the sphere is outside of our computational mesh $\Gamma_{h}$.

Our initial geometry is a cylinder around the $x_{3}$-axis of radius $1.0$, spanning from $x_{3}=-0.1$ to $x_{3}=0.1$, consisting of $3,686$ triangles. Hence, we start infeasible with respect to all constraints. The top cap of the cylinders is subject to the touch condition based on $\operatorname{dist}(x,\Gamma_{\text{S}_{1}})$, while the bottom cap has to adhere $\operatorname{dist}(x,\Gamma_{\text{S}_{2}})=0$ as a constraint. We first conduct $5$ steps using the approximate Newton-Scheme, where the second order partial shape derivative in the KKT-System, i.e., the corresponding diagonal block, is approximated with the $H^{1}$-inner product (19) with smoothing $\gamma=1$. After $5$ steps, we remesh $\Gamma_{h}$ using the compounding of gmsh 3.0.6 [28, 35]. After the remeshing, we switch to the proper Newton iteration. As the geometric situation under consideration is rotationally invariant, we also have to restrict the motion of the vertices once we move to the Newton scheme as discussed in (20). As a consequence, interior vertices of the cap and bottom of the original cylinder, or their descendants after remeshing, can no longer slide over the obstacles once we move to the Newton scheme. However, we found that after the first initial $5$ gradient steps, where sliding is possible, the unknown surface $\Gamma$ already clings to the solid bodies in a well-behaved manner.

Initial and final mesh are shown in Figure 2 and the respective quantities of interest are listed in Tab. 1. The convergence of the optimization residual is shown in Figure 3.

Our methodology achieves very accurate results in comparison to the reference data in [29]. In particular, the relative error in the computed total force is less than $0.5\%$ for all meshes considered, in particular including the coarsest grid of only $13,728$ cells. We also observe excellent convergence speed of the Newton method, which is activated after one remeshing in iteration 5 for all grids. Observing the convergence history in Figure 3, one can see that the reduction of the residual with the Newton method is indeed very rapid, but not fully grid independent. A possible source for this slight mesh dependency can possibly found in the subset constraint, which is implemented discretely on a per-vertex level without considering a possible analytic or function space equivalent of the subset constraint.

We also use the same sphere-plate setting to study the non-axisymmetric situation. To this end, we apply gravity in the negative $x_{2}$-direction with varying bond numbers of $8.0$, $4.0$, $1.0$ and $0.5$. We did not obtain solutions for bond numbers $10.0$ and up, as these would lead to topological changes with the droplet splitting and detaching in part from the obstacle. The respective shapes are shown in Figure 4.

In particular, [29] states that the gravitational effect is negligible provided that the Bond number satisfies $b\ll 1$. For comparison, we have tabulated the computed force for different Bond numbers in Table 2.

We revisit the situation of a sphere of unit radius $1.0$ placed over a plane with a non-dimensional distance of $2.2$. The desired contact angles are $30^{\circ}$ on the sphere and $80^{\circ}$ on the plane. This setting was studied analytically in [39], where four unduloid bridges are given, all satisfying the same non-dimensional volume of $V=3.6$, a constant curvature and the desired touch and angle constraints with the obstacles. Hence, they are critical shapes of (3). However, as mentioned in [39], it is unclear which of these candidates constitutes a physical solution. The respective unduloid bridges do not admit a closed-form description of their geometry. However, the authors of [39] kindly provided us with a polygonal graph of two of their unduloids as a cut through the $z=0$ plane: The unduloid with a meniscus of $32^{\circ}$ was given as a plane graph of $363$ points and the unduloid with a meniscus angle of $47.2^{\circ}$ was kindly provided as a plane graph of $443$ points. We then use the compounding functionality of gmsh 3.0.6 [28, 35] to create a shell mesh of the resulting bodies of revolution with optimal point spacing independent of the original points provided in the graphs. These two shells are then used as the starting geometry in our program.

Due to floating-point inaccuracies, a saddle point cannot be expected to be a stable point for a numerical gradient descent scheme. Indeed, the unduloid shape for a meniscus angle of $32^{\circ}$ proved to be almost stationary initially. However, a rapid descent of the objective was eventually achieved with the geometry indicating the desire to split the single capillary bridge into two droplets. This ultimately terminates the program as topological changes are not possible, yet. As expected, the Newton-iteration was unusable for this geometry due to the ill-conditioning of the Hessian. Contrary to this saddle-point, a Newton-scheme is possible for the initial geometry of the unduloid of meniscus angle $47.2^{\circ}$, where the provided starting guess of the discretized analytical shape quickly converges to a residuum of $7.940805\cdot 10^{-9}$, numerically identifying this shape as a proper local minimum. The resulting forces on the sphere were computed to be $F_{p}=-2.140941$ and $F_{s}=4.481142$. The respective shapes are shown in Figure 5.