Arbitrary Lagrangian--Eulerian finite element method for curved and deforming surfaces. I. General theory and application to fluid interfaces

TL;DR

This paper introduces an ALE finite element method for curved, deforming surfaces, enabling stable simulations of fluid interfaces with complex geometries and flows, validated through numerical benchmarks and physical insights.

Contribution

It develops a general ALE finite element framework for curved surfaces, addressing numerical instabilities and applying it to fluid interface stability analysis.

Findings

ALE method effectively handles curved, deforming surfaces.

Surface tension projection stabilizes in-plane incompressibility.

Cylindrical fluid films are stable to non-axisymmetric perturbations.

Abstract

An arbitrary Lagrangian--Eulerian (ALE) finite element method for arbitrarily curved and deforming two-dimensional materials and interfaces is presented here. An ALE theory is developed by endowing the surface with a mesh whose in-plane velocity need not depend on the in-plane material velocity, and can be specified arbitrarily. A finite element implementation of the theory is formulated and applied to curved and deforming surfaces with in-plane incompressible flows. Numerical inf--sup instabilities associated with in-plane incompressibility are removed by locally projecting the surface tension onto a discontinuous space of piecewise linear functions. The general isoparametric finite element method, based on an arbitrary surface parametrization with curvilinear coordinates, is tested and validated against several numerical benchmarks. A new physical insight is obtained by applying the ALE developments to cylindrical fluid films, which are computationally and analytically found to be stable to non-axisymmetric perturbations, and unstable with respect to long-wavelength axisymmetric perturbations when their length exceeds their circumference. A Lagrangian scheme is attained as a special case of the ALE formulation. Though unable to model fluid films with sustained shear flows, the Lagrangian scheme is validated by reproducing the cylindrical instability. However, relative to the ALE results, the Lagrangian simulations are found to have spatially unresolved regions with few nodes, and thus larger errors.

Equations

$$\mathbf{x}({\theta }^{\alpha },t)=\mathbf{x}({\theta }^{\alpha }({\xi }^{\beta },t),t)=\hat{\mathbf{x}}({\xi }^{\beta },t),$$

$$\bm{v}\,=\,\dfrac{\partial\theta^{\alpha}}{\partial t}\Big{\rvert}_{\xi^{\beta}}\dfrac{\partial\bm{x}}{\partial\theta^{\alpha}}\,+\,\dfrac{\partial\bm{x}}{\partial t}\Big{\rvert}_{\theta^{\alpha}}\,=\,\dfrac{\partial\theta^{\alpha}}{\partial t}\Big{\rvert}_{\xi^{\beta}}\bm{a}_{\alpha}\,+\,\dfrac{\partial\bm{x}}{\partial t}\Big{\rvert}_{\theta^{\alpha}}~{}.$$

$$\mathbf{x}=\hat{\mathbf{x}}({\xi }^{\alpha },t)=\mathbf{x}({\theta }^{\alpha },t)=\check{\mathbf{x}}({\zeta }^{\alpha },t),$$

$$\check{\mathbf{x}}({\zeta }^{\alpha },t)=\check{\mathbf{x}}({\zeta }^{\alpha }({\xi }^{\beta },t),t)=\hat{\mathbf{x}}({\xi }^{\beta },t),$$

$$\bm{v}\,=\,\dfrac{\partial\zeta^{\alpha}}{\partial t}\Big{\rvert}_{\xi^{\beta}}\,\dfrac{\partial\check{\bm{x}}}{\partial\zeta^{\alpha}}\,+\,\dfrac{\partial\check{\bm{x}}}{\partial t}\Big{\rvert}_{\zeta^{\alpha}}\,=\,\dfrac{\partial\zeta^{\alpha}}{\partial t}\Big{\rvert}_{\xi^{\beta}}\,\bm{a}_{\check{\alpha}}\,+\,\dfrac{\partial\check{\bm{x}}}{\partial t}\Big{\rvert}_{\zeta^{\alpha}}~{},$$

$$\bm{v}^{\mathrm{m}}\,:=\,\dfrac{\partial\check{\bm{x}}}{\partial t}\Big{\rvert}_{\zeta^{\alpha}}\,=\,\bm{x}^{\prime}~{}.$$

$$c^{\check{\alpha}}\,:=\,\dfrac{\partial\zeta^{\alpha}}{\partial t}\Big{\rvert}_{\xi^{\beta}}\,=\,\dot{\zeta}^{\alpha}\hskip 20.0pt\mathrm{and}\hskip 23.0pt\bm{c}\,:=\,c^{\check{\alpha}}\bm{a}_{\check{\alpha}}~{},$$

$$\mathbf{v}=\mathbf{c}+{\mathbf{v}}^{\mathrm{m}}.$$

$$\dfrac{\mathrm{d}}{\mathrm{d}t}(\,~{}\,)\,=\,(\,~{}\,)^{\prime}\,+\,c^{\check{\alpha}}\,(\,~{}\,)_{,{\check{\alpha}}}\,=\,(\,~{}\,)^{\prime}\,+\,\big{(}\bm{v}-\bm{v}^{\mathrm{m}}\big{)}\cdot\bm{a}^{\check{\alpha}}\,(\,~{}\,)_{,{\check{\alpha}}}~{},$$

$$\dot{\bm{v}}\,=\,\bm{v}^{\prime}\,+\,c^{\check{\alpha}}\,\bm{v}_{,{\check{\alpha}}}\,=\,\bm{v}^{\prime}\,+\,\big{(}\bm{v}_{,{\check{\alpha}}}\otimes\bm{a}^{\check{\alpha}}\big{)}\,\big{(}\bm{v}-\bm{v}^{\mathrm{m}}\big{)}~{}.$$

$$\check{\mathbf{x}}({\zeta }^{\alpha },t)=\check{\mathbf{x}}({\zeta }^{\alpha },t_0)+{\int }_{t_0}^t{\mathbf{v}}^{\mathrm{m}}({\zeta }^{\alpha },t^{\prime })\mathrm{d}{t}^{\prime },$$

$${\mathbf{v}}^{\mathrm{m}}\cdot \mathbf{n}=\mathbf{v}\cdot \mathbf{n},$$

$${\mathbf{v}}^{\mathrm{m}}\cdot {\mathbf{a}}_{\check{\alpha }}=0,$$

$${\mathbf{v}}^{\mathrm{m}}=(\mathbf{n}\otimes \mathbf{n})\mathbf{v}.$$

$${\mathbf{a}}^{\check{\alpha }}\cdot {\mathbf{v}}_{,\check{\alpha }}=v_{;\check{\alpha }}^{\check{\alpha }}-2vH=0.$$

$$\rho \dot{\mathbf{v}}=\rho \mathbf{b}+{\mathbf{T}}_{;\check{\alpha }}^{\check{\alpha }},$$

$${\sigma }^{\check{\alpha }\check{\beta }}=\lambda a^{\check{\alpha }\check{\beta }}+{\pi }^{\check{\alpha }\check{\beta }}.$$

$$\pi^{{\check{\alpha}}{\check{\beta}}}\,=\,\zeta\,\dot{a}_{\mu\nu}\,a^{\alpha\mu}\,a^{\beta\nu}\,=\,\zeta\,\bm{v}_{,{\check{\gamma}}}\cdot\big{(}\bm{a}^{\check{\alpha}}\,a^{{\check{\gamma}}{\check{\beta}}}+\,\bm{a}^{\check{\beta}}\,a^{{\check{\gamma}}{\check{\alpha}}}\big{)}~{},$$

$$\rho\dot{\bm{v}}\cdot\bm{a}_{\check{\alpha}}=\rho b_{\check{\alpha}}+\lambda_{,{\check{\alpha}}}+\zeta\,\Big{(}a^{{\check{\beta}}{\check{\gamma}}}\,v_{{\check{\alpha}};{\check{\beta}}{\check{\gamma}}}+K\,v_{\check{\alpha}}+2\,v_{,{\check{\alpha}}}\,H-2\,v_{,{\check{\beta}}}\,b^{\check{\beta}}_{\check{\alpha}}-2\,v\,H_{,{\check{\alpha}}}\Big{)}~{},$$

$$\rho\dot{\bm{v}}\cdot\bm{n}=p+2\,\lambda\,H+\zeta\,\Big{(}2\,b^{{\check{\alpha}}{\check{\beta}}}\,v_{{\check{\alpha}};{\check{\beta}}}-4\,v\,\big{(}2H^{2}-K\big{)}\Big{)}~{}.$$

$$\mathbf{v}=\bar{\mathbf{v}}\text{on}{\partial }^{\text{v}}\mathcal{P}\text{and}\mathbf{T}=0,v=\bar{v}\text{on}{\partial }^{\text{t}}\mathcal{P},$$

$$\int_{\mathscr{P}}\big{(}\ldots\big{)}~{}\mathrm{d}a\,=\,\int_{\Omega}\big{(}\ldots\big{)}~{}J^{\mathrm{m}}~{}\mathrm{d}\Omega~{}.$$

$$\int_{\partial\!\mathscr{P}}\big{(}\ldots\big{)}~{}\mathrm{d}s\,=\,\int_{\Gamma}\big{(}\ldots\big{)}~{}J^{\mathrm{m}}_{\Gamma}~{}\mathrm{d}\Gamma~{},$$

$${\int }_{\Omega }\delta \mathbf{v}\cdot \rho \dot{\mathbf{v}}{J}^{\mathrm{m}}\mathrm{d}\Omega -{\int }_{\Omega }\delta \mathbf{v}\cdot \rho \mathbf{b}{J}^{\mathrm{m}}\mathrm{d}\Omega -{\int }_{\Omega }\delta \mathbf{v}\cdot {\mathbf{T}}_{;\check{\alpha }}^{\check{\alpha }}{J}^{\mathrm{m}}\mathrm{d}\Omega =0\forall \delta \mathbf{v}\in {\mathcal{V}}_0,$$

$$-{\int }_{\Omega }\delta \mathbf{v}\cdot {\mathbf{T}}_{;\check{\alpha }}^{\check{\alpha }}{J}^{\mathrm{m}}\mathrm{d}\Omega ={\int }_{\Omega }\delta {\mathbf{v}}_{,\check{\alpha }}\cdot {\mathbf{T}}^{\check{\alpha }}{J}^{\mathrm{m}}\mathrm{d}\Omega -{\int }_{\Gamma }\delta \mathbf{v}\cdot \mathbf{T}{J}_{\Gamma }^{\mathrm{m}}\mathrm{d}\Gamma ,$$

$$-{\int }_{\Omega }\delta \mathbf{v}\cdot {\mathbf{T}}_{;\check{\alpha }}^{\check{\alpha }}{J}^{\mathrm{m}}\mathrm{d}\Omega ={\int }_{\Omega }\delta {\mathbf{v}}_{,\check{\alpha }}\cdot {\mathbf{T}}^{\check{\alpha }}{J}^{\mathrm{m}}\mathrm{d}\Omega .$$

$$\begin{array}{rl}{\mathcal{G}}_{\mathrm{v}}:=& {\int }_{\Omega }\delta \mathbf{v}\cdot \rho \dot{\mathbf{v}}{J}^{\mathrm{m}}\mathrm{d}\Omega +{\int }_{\Omega }\delta {\mathbf{v}}_{,\check{\alpha }}\cdot {\pi }^{\check{\alpha }\check{\beta }}{\mathbf{a}}_{\check{\beta }}{J}^{\mathrm{m}}\mathrm{d}\Omega \\ & +{\int }_{\Omega }\delta {\mathbf{v}}_{,\check{\alpha }}\cdot {\mathbf{a}}^{\check{\alpha }}\lambda J^{\mathrm{m}}\mathrm{d}\Omega -{\int }_{\Omega }\delta \mathbf{v}\cdot \rho \mathbf{b}{J}^{\mathrm{m}}\mathrm{d}\Omega =0\forall \delta \mathbf{v}\in {\mathcal{V}}_0.\end{array}$$

$$\mathcal{G}_{\mathrm{m}}\,:=\,\alpha^{\mathrm{m}}\!\int_{\Omega}\delta\bm{v}^{\mathrm{m}}\cdot\Big{[}\bm{v}^{\mathrm{m}}-\big{(}\bm{n}\otimes\bm{n}\big{)}\bm{v}\Big{]}~{}J^{\mathrm{m}}~{}\mathrm{d}\Omega\,=0\hskip 15.0pt\forall\,\,\delta\bm{v}^{\mathrm{m}}\in\mathcal{V}~{}.$$

$$\mathcal{G}_{\lambda}\,:=\,\int_{\Omega}\delta\lambda\,\big{(}\bm{a}^{\check{\alpha}}\cdot\bm{v}_{,{\check{\alpha}}}\big{)}~{}J^{\mathrm{m}}~{}\mathrm{d}\Omega\,=0\hskip 15.0pt\forall\,\,\delta\lambda\in\Lambda~{}.$$

$${\int }_{\Omega }\delta \mathbf{v}\cdot \rho \dot{\mathbf{v}}{J}^{\mathrm{m}}\mathrm{d}\Omega +{\int }_{\Omega }\delta {\mathbf{v}}_{,\check{\alpha }}\cdot {\pi }^{\check{\alpha }\check{\beta }}{\mathbf{a}}_{\check{\beta }}{J}^{\mathrm{m}}\mathrm{d}\Omega +{\int }_{\Omega }\delta {\mathbf{v}}_{,\check{\alpha }}\cdot {\mathbf{a}}^{\check{\alpha }}\lambda J^{\mathrm{m}}\mathrm{d}\Omega -{\int }_{\Omega }\delta \mathbf{v}\cdot \rho \mathbf{b}{J}^{\mathrm{m}}\mathrm{d}\Omega$$

Full text

\Hy@raisedlink

${}^{\ddagger}\,$amaresh.sahu*@*berkeley.edu

\Hy@raisedlink    ${}^{\S}\,$yannick.omar*@*berkeley.edu

\Hy@raisedlink    ${}^{\flat}\,$sauer*@*aices.rwth-aachen.de

\Hy@raisedlink    ${}^{\dagger}\,$kranthi*@*berkeley.edu

1 Introduction

In this paper, and the subsequent manuscript in the series​***From now on, we refer to the present paper as “Part I” and the following one [1] as “Part II.” [1], we develop an arbitrary Lagrangian–Eulerian (ALE) theory for arbitrarily curved and deforming two-dimensional interfaces with in-plane fluidity. The theory is based on a surface discretization which is independent of the in-plane material flow, such that the surface mesh need not convect with the material. Consequently, two-dimensional materials with large in-plane flows on arbitrarily deforming surfaces can be modeled. In Part I, we develop our theory and use standard numerical techniques to devise an isoparametric ALE finite element method for incompressible fluid films. We then implement the finite element formulation, model the deformations and flows of such materials over time, and provide several numerical results for both flat and cylindrical geometries. In Part II, we hope to extend the finite element formulation to lipid membranes and study membrane behavior in several biologically relevant situations. As the equations governing single- and multi-component lipid membranes reduce to the fluid film equations in the limit where no elastic energy is stored in the membrane, such a separation is natural and allows us to present our results in a more accessible manner.

Two-dimensional fluids have played an increasingly important role in many engineering applications, in which they often arise at phase boundaries in multiphase systems [2]. For example, under the influence of gravity and capillary forces, foams will drain over time until the constituent bubbles burst [3]. Foam lifetime plays a key role in their viability for engineering applications, and there have consequently been many efforts to improve foam stability [2]. Similar efforts have been made to stabilize emulsions and colloidal dispersions, which again are of much industrial value [4, 5]. Surfactants are often used to stabilize vapor–liquid and liquid–liquid interfaces by lowering the local surface tension [6]; surface tension gradients can drive Marangoni flows [7, 8, 9] and in some cases have been shown to significantly affect material properties [10].

Two-dimensional materials with in-plane fluidity also play a fundamental role in biology. Biological membranes, which are interfaces composed of lipids and proteins, are in-plane fluid and out-of-plane elastic materials [11]. They make up the boundary of the cell as well as many of its internal organelles, including the nucleus, endoplasmic reticulum, and Golgi complex. Lipid membranes thus add structure and organization to the cell, and furthermore play an important role in many cellular processes. Endocytosis, for example, begins when proteins in the surrounding bulk fluid bind to the cell membrane’s constituent lipids and proteins at a specific location. The membrane forms an initially shallow invagination, which then develops into a mature bud and eventually pinches off into a membrane-bound vesicle that enters the cell [12]. The vesicular membrane contains lipids and proteins which were previously on the cell boundary, and furthermore the vesicle may enclose nutrients or other cargo. Endocytosis is thus a key process in transferring nutrients to the cell, regulating the expression of proteins on the cell surface, and cell homeostasis [13]. It involves nontrivial coupling between protein binding and unbinding reactions, in-plane lipid flow, and out-of-plane membrane shape changes. In particular, the in-plane flow and out-of-plane bending are coupled because lipid membranes are nearly area-incompressible [11] and therefore lipids are required to flow in-plane to accommodate any shape changes. In another example, lipid membranes can phase separate into liquid–ordered $(\mathrm{L}_{\mathrm{o}})$ and liquid–disordered $(\mathrm{L}_{\mathrm{d}})$ domains under physiological conditions [14]; the energetic penalty of the $\mathrm{L}_{\mathrm{o}}$–$\mathrm{L}_{\mathrm{d}}$ interface plays a major role in the fusion of HIV-containing vesicles with target immune cells [15]. This phenomenon demonstrates the value in understanding the coupling between elastic membrane shape changes and the thermodynamically irreversible processes of in-plane lipid flow and in-plane species diffusion.

The interfacial materials discussed thus far are of fundamental importance to engineering and biology. Consequently, there have been significant theoretical efforts to better understand their physics. The pioneering work of L.E. Scriven [16] and R. Aris [17] was crucial to our current understanding of interfacial flows. In particular, Scriven recognized it is prohibitively difficult to use standard Cartesian, cylindrical, and spherical coordinate systems to solve for fluid flows on arbitrarily curved surfaces, where even expressing the surface Laplacian of the velocity field at every point on the surface is nontrivial. Consequently, Scriven used a mathematically elegant differential geometric framework to naturally represent two-dimensional flows and their gradients on arbitrarily curved surfaces [16]. Subsequently, Aris worked to describe three-dimensional fluids using the machinery of differential geometry, and incorporated Scriven’s surface flow description within his general differential geometric perspective [17]. The powerful formalism developed by Scriven and Aris continue to be in widespread use today. An excellent review of the interfacial dynamics of fluid interfaces in multiphase systems is provided in Ref. [6], and for a wonderful perspective on interfaces in fluid mechanics see Ref. [18].

While the equations of motion characterizing two-dimensional interfacial flows on surfaces are now widespread and well-understood, theoretical developments for lipid membranes are in a less mature stage. A major complexity arises in modeling lipid membranes because they behave as in-plane fluids, out-of-plane elastic solids, and the surface on which dynamical equations are to be written is itself curved and deforming over time. Early membrane models were modifications of P.M. Naghdi’s seminal contributions to shell theory. However, while Naghdi used a balance law formulation [19], the first membrane models used variational methods and focused only on elastic membrane behavior. In particular, P. Canham [20] and W. Helfrich [21] proposed an elastic membrane bending energy in the early 1970’s; Helfrich also used variational methods to determine the Euler–Lagrange equations governing axisymmetric membrane shapes in the absence of in-plane flows. The Euler–Lagrange equations, which by construction include only thermodynamically reversible phenomena and thus do not contain viscous forces, were not extended to non-axisymmetric settings until 1999 [22]. However, by this time various other models which restricted membrane shapes to small deviations from flat planes [23, 24, 25, 26] and cylindrical or spherical shells [26, 27] had also emerged. Since then, variational methods encompassing different physical phenomena have continued to be developed [28, 29, 30, 31, 32, 33, 34, 35, 36]. In a parallel development, in-plane velocities were included in some models about simple geometries [26, 34, 37]. It was not until 2009, however, that the general equations governing a single-component, arbitrarily curved and deforming lipid membrane with in-plane viscosity were determined [38]—using a combination of variational methods to determine elastic contributions and the so-called Rayleigh dissipation potential to determine the viscous terms. Since then, variational methods have been extended, with viscous stresses sometimes included in an ad-hoc manner [39, 40, 41]. Membrane models have also recently been developed by building on the work of Naghdi [19] and using fundamental balance laws and associated constitutive equations [42, 43, 44, 45].

While such theoretical developments have had success in modeling certain membrane phenomena, they were difficult to extend to study how elastic out-of-plane membrane bending is coupled to different irreversible phenomena, such as in-plane lipid flow, in-plane phase transitions involving multiple components, and chemical reactions between membrane components and species in the surrounding bulk. Our recent work [46], inspired by the pioneering works of I. Prigogine [47], L. Onsager [48, 49], and S.R. de Groot & P. Mazur [50], developed the general theory of irreversible thermodynamics for arbitrarily curved lipid membranes, provided a formalism to determine the equations governing membrane dynamics, and presented comprehensive models for all of the irreversible phenomena described thus far.

Though the equations governing both fluid interfaces and lipid membranes are now determined, the equations are highly nonlinear and in general cannot be solved analytically. Our work entails developing an ALE theory for two-dimensional materials with in-plane fluidity. The theory involves a surface discretization whose in-plane velocity can be (i) zero, as in an Eulerian formulation, (ii) equal to the in-plane material velocity, as in a Lagrangian formulation, or (iii) specified arbitrarily. The flexibility of our ALE theory, as well as its similarities to bulk methods of the same name [51], explain our nomenclature. With the ALE theory, we numerically solve the equations of motion governing the aforementioned materials of interest. We split this effort into two pieces: in Part I we derive the general ALE theory, develop its finite element formulation, and apply it to two-dimensional fluid films. In Part II we extend the finite element formulation to lipid membranes, which elastically resist out-of-plane bending, and present results from our numerical simulations.

The challenges in theoretically modeling deforming interfaces with in-plane flow, and lipid membranes especially, extend to their numerical modeling as well: to model a material with arbitrarily large shape deformations and in-plane flows, standard techniques from fluid mechanics and solid mechanics are insufficient. Regarding fluid films, many previous studies simplified the problem by assuming the film was fixed in space. The resultant fixed-surface flow equations, derived by Scriven [16], have been solved using various methods: by modeling the fluid interface as a level set in $\mathbb{R}^{3}$ [52], with projection-based finite element methods [53], and through the discretization of exterior calculus operators [54]. On the other hand, a different study using level set methods [55] made considerable advances in numerically modeling bubble deformation and breakup in foams, however they separated the dynamics into different steps and in each step made simplifying assumptions. In addition, interfaces have been modeled as the boundary between bulk fluid domains, using both ALE and level set techniques—however, such works are either restricted to simple geometries [56] or do not include in-plane interfacial flow [57, 58, 59, 60]. ALE methods were also used to study the evolution of scalar fields on surfaces whose time evolution is known [61]. While the numerical methods discussed thus far have modeled different fluid film phenomena, they do not seem easily amenable to the study of general deforming fluid interfaces or lipid membranes, with the latter having their own constitutive behavior.

Just as in the case of fluid films, several lipid membrane studies assumed the membrane shape to be fixed, and under these conditions studied how surface flows are coupled to flows in the surrounding bulk in two cases: spherical surfaces with protein inclusions [62] and radial surfaces in a one-to-one correspondence with a sphere [63]. Alternatively, many of the studies modeling the deformation of lipid membranes [64, 65, 66, 67, 68, 69, 70, 71, 72] consider only the Euler–Lagrange equations, and thus predict membrane shapes without knowledge of the in-plane flow. However, as the in-plane flow and out-of-plane deformations are coupled through the in-plane viscosity [46], the predictions of the previously mentioned works are only physically relevant in the limit where velocities are negligible. Another approach has been to include in-plane fluid flow and limit the membrane to remain in one-to-one correspondence with a flat plane [42, 44]. While such an approach is theoretically sound, it is limited in its use as it cannot, for example, model the large shape deformations observed in endocytosis [12]. Several recent works have avoided the computational complexity of modeling the full lipid membrane equations by assuming only axisymmetric shapes [45, 73], however this turns out to be a poor assumption which in many cases yields incorrect results [74]. In our previous work we modeled the full non-axisymmetric membrane equations using a Lagrangian finite element method [75, 74], which is computationally valid yet can attain locally singular Jacobians and uninvertible matrices when there are moderate in-plane flows. Lagrangian methods are thus not suitable for the study of general fluid and lipid membrane phenomena.

The limitations of our Lagrangian finite element formulation and other computational techniques in modeling fluid interfaces and lipid membranes motivate our development of an ALE finite element method for curved and deforming surfaces. The following aspects are new in this work: we

develop an ALE theory, within a differential geometric setting, for general arbitrarily curved and deforming two-dimensional interfaces with in-plane flow, 2. 2.

apply the theory to two-dimensional fluid films and derive a corresponding isoparametric, fully implicit ALE finite element method, 3. 3.

prevent numerical inf–sup instabilities associated with the in-plane areal incompressibility by adapting the method of C.R. Dohrmann and P.B. Bochev [76] to curved surfaces, 4. 4.

numerically simulate an arbitrarily curved and deforming fluid film, from which we find a physical instability that is confirmed analytically with a linear stability analysis, and 5. 5.

demonstrate how our ALE formulation can be altered to yield a Lagrangian scheme as a special case.

As mentioned earlier, we limit our numerical calculations to fluid films in this manuscript, as the extensions of the theory and numerical methods to lipid membranes will be presented in Part II [1]. We note Ref. [77] describes a concurrent effort with a similar objective. In particular, Ref. [77] also derives a general ALE theory. However, their numerical implementation is based on a surface parametrization involving small membrane deformations, generalized to arbitrarily curved surfaces, with periodic updates of the reference surface as required. Additionally, Ref. [77] uses a Hodge decomposition of the membrane velocity field, while the present work employs isoparametric finite element methods.

Our paper is organized as follows: In Sec. 2, we present our ALE theory for general two-dimensional interfaces. We provide the equations of motion governing arbitrarily curved and deforming fluid films in Sec. 3 and develop the corresponding finite element formulation in Sec. 4. Numerical results of an ALE implementation are presented in Sec. 5; the modifications leading to a pure Lagrangian formulation and corresponding results are provided in Sec. 6. We end with conclusions and avenues for future work in Sec. 7. Several of the more detailed calculations regarding fluid films and our finite element implementation, as well as additional numerical benchmarks, are relegated to Appendices A–D. Relevant movies are provided in Appendix E, and important symbols are listed in Appendix F.

2 Arbitrary Lagrangian–Eulerian Theory

While the equations of motion governing arbitrarily curved and deforming fluid films and lipid membranes were presented in Refs. [16, 38, 42, 46], solving the resultant equations is nontrivial. These equations are highly nonlinear and cannot be solved analytically, yet many issues arise when trying to solve them numerically as well. We have described the shortcomings of our previous Lagrangian finite element formulation [74, 75], which is not appropriate for materials with in-plane flow. Namely, when the surface is discretized and the corresponding mesh travels in-plane with the material, mesh elements become highly distorted and attain nearly singular Jacobians (see Fig. 1(a)). In such a case, a simple example of vortex flows is not possible. We therefore require a new numerical method to solve the equations of motion.

In this section, we derive an ALE theory for arbitrarily curved and deforming surfaces, which provides equations more amenable to numerical solution. In particular, we seek a description of the material surface which can be easily discretized, with individual elements not undergoing large distortions when material flows in-plane. To this end, we introduce an ALE parametrization of the two-dimensional material, which endows the surface with a mesh that deforms in the normal direction with the material—yet whose in-plane motion can be specified arbitrarily and need not depend on the material flow. We note Ref. [61] introduced a mesh evolving in a similar manner, for the study of scalar fields on evolving surfaces whose time evolution is prescribed. Our ALE description, on the other hand, introduces three new unknowns corresponding to the three components of the mesh velocity. In what follows, we discuss various parametrizations of the surface, and in the ALE case provide the additional three equations required for our problem to be mathematically well-posed. We also describe how the equations governing fluid films and lipid membranes, which are based on an Eulerian surface parametrization, are modified when the ALE parametrization is employed.

2.1 Surface Geometry

As discussed in Refs. [42, 46], an arbitrarily curved and deforming patch of surface $\mathscr{P}$ can be parametrized by either the convected coordinates $\xi^{\alpha}$, which are attached to material points and are convected with the material, or the surface-fixed coordinates $\theta^{\alpha}$, which are defined such that a point of fixed $\theta^{\alpha}$ moves only normal to the surface. Schematics of the movement of points of constant $\xi^{\alpha}$ and $\theta^{\alpha}$ are provided in Figs. 1(a) and 1(b), respectively. Here and from now on, we prescribe Greek indices to span the set $\{1,2\}$, and use the Einstein summation convention in which Greek indices repeated in a subscript and superscript are summed over. As the material occupying a point of constant $\theta^{\alpha}$ will in general change in time, we formally write $\theta^{\alpha}=\theta^{\alpha}(\xi^{\beta},\,t)$ and express the position in terms of either surface-fixed or convected coordinates as

[TABLE]

where the ‘hat’ accent is used to denote the position expressed in terms of the $\xi^{\alpha}$ parametrization (see Fig. 2).

As discussed below, the $\theta^{\alpha}$ parametrization yields a surface description which is in-plane Eulerian and out-of-plane Lagrangian, and is most natural to theoretically model lipid membranes—which are in-plane fluids and out-of-plane elastic solids. Consequently, our theoretical developments [46] used the surface-fixed parametrization throughout. We use the same notation in this work, and review it here. Partial and covariant differentiation with respect to $\theta^{\alpha}$ are respectively denoted $(\,~{}\,)_{,\alpha}$ and $(\,~{}\,)_{;\alpha}$. The surface-fixed parametrization yields the in-plane basis vectors $\bm{a}_{\alpha}:=\bm{x}_{,\alpha}$ and unit normal $\bm{n}:=(\bm{a}_{1}\times\bm{a}_{2})/\lvert\bm{a}_{1}\times\bm{a}_{2}\rvert$. The metric and curvature components are respectively given by $a_{\alpha\beta}\,:=\,\bm{a}_{\alpha}\cdot\bm{a}_{\beta}$ and $b_{\alpha\beta}:=\bm{n}\cdot\bm{x}_{,\alpha\beta}$, with which the mean and Gaussian curvatures are respectively calculated as $H:=\tfrac{1}{2}a^{\alpha\beta}b_{\alpha\beta}$ and $K:=\det(b_{\alpha\beta})/\det(a_{\alpha\beta})$. A comprehensive geometric description of the surface can be found in Ref. [46] and the references provided therein.

2.2 Surface Kinematics

As a point of constant $\xi^{\alpha}$ follows a material point over time, the velocity $\bm{v}$ of a material point is defined as $\bm{v}:=[\partial\hat{\bm{x}}/\partial t]\rvert_{\xi^{\beta}}$, which upon substitution of Eq. (1) and application of the chain rule yields

[TABLE]

At this point, the parametrization $\theta^{\alpha}$ is defined such that the second term on the right-hand side of Eq. (2) lies entirely in the normal direction. In this case the normal component of the velocity, $v=\bm{v}\cdot\bm{n}$, satisfies $v\bm{n}=[\partial\bm{x}/\partial t]\rvert_{\theta^{\alpha}}$—and a point of constant $\theta^{\alpha}$ is unaffected by in-plane flow, as shown in Fig. 1(b). We accordingly refer to $\theta^{\alpha}$ as a surface-fixed parametrization. The in-plane contravariant velocity components $v^{\alpha}$ are defined as $v^{\alpha}:=[\partial\theta^{\alpha}/\partial t]\rvert_{\xi^{\beta}}$, such that the velocity $\bm{v}$ can be written as $\bm{v}=v^{\alpha}\bm{a}_{\alpha}+v\bm{n}$—indicating our definitions of the normal velocity $v$ and contravariant velocity components $v^{\alpha}$ are consistent with the geometric description of the surface.

The two different parametrizations introduced thus far offer perspectives analogous to the familiar Lagrangian and Eulerian formulations from standard continuum mechanics. A point of constant $\xi^{\alpha}$ is a material point, so the convected coordinates provide a Lagrangian perspective. A point of constant $\theta^{\alpha}$, on the other hand, is independent of the in-plane surface flow and so the surface-fixed coordinates provide an in-plane Eulerian perspective. The material time derivative is calculated as $\mathrm{d}(\,~{}\,)/\mathrm{d}t:=[\partial(\,~{}\,)/\partial t]\rvert_{\xi^{\beta}}$ in the Lagrangian perspective, and as $\mathrm{d}/\mathrm{d}t(\,~{}\,)=v^{\alpha}(\,~{}\,)_{,\alpha}+(\,~{}\,)_{,t}$ for scalar quantities in the in-plane Eulerian perspective, where the partial time derivative is defined as $(\,~{}\,)_{,t}:=[\partial(\,~{}\,)/\partial t]\rvert_{\theta^{\alpha}}$. We often denote the material time derivative of a quantity with a dot over that quantity, as in $\bm{v}=\dot{\bm{x}}$ or $v^{\alpha}=\dot{\theta}^{\alpha}$ . By applying the material time derivative to the basis vectors $\bm{a}_{\alpha}$ and $\bm{n}$, we find $\dot{\bm{a}}_{\alpha}=\bm{v}_{,\alpha}$ and $\dot{\bm{n}}=-(\bm{a}^{\alpha}\otimes\bm{n})\,\bm{v}_{,\alpha}$, where $\otimes$ denotes the dyadic or outer product. We calculate the material time derivatives of the metric and curvature components as $\dot{a}_{\alpha\beta}=\bm{v}_{,\alpha}\cdot\bm{a}_{\beta}+\bm{v}_{,\beta}\cdot\bm{a}_{\alpha}$ and $\dot{b}_{\alpha\beta}=\bm{v}_{;\alpha\beta}\cdot\bm{n}$. A more detailed description of the surface geometry and kinematics is provided in our past theoretical work [46].

2.3 ALE Description and Geometry

Thus far, the convected coordinates $\xi^{\alpha}$ and surface-fixed coordinates $\theta^{\alpha}$ were introduced to describe a surface with in-plane fluidity. Lagrangian numerical methods based on convected coordinates are conceptually simpler to develop, yet they cannot capture in-plane flows—which highly distort mesh elements, as discussed previously. Surface-fixed coordinates are most natural for a theoretical description of arbitrarily curved and deforming surfaces with in-plane fluidity, and their numerical implementation requires nodes to only move orthogonally to the surface (see Fig. 1(b)). A numerical method based on the $\theta^{\alpha}$ parametrization can model arbitrarily large in-plane flows, yet in rare instances the discretized nodes move towards one another as the surface deforms, also shown in Fig. 1(b). In the spirit of formulating a completely general numerical method to describe evolving surfaces with in-plane fluidity, we begin by developing an ALE surface parametrization, denoted by the ALE coordinate $\zeta^{\alpha}$, which endows the surface with a mesh. The main idea is to define $\zeta^{\alpha}$ such that the corresponding mesh deforms out-of-plane with the material, while its in-plane motion can be specified arbitrarily (see Fig. 1(c)). We now describe how the geometric description introduced in Sec. 2.1 is modified when expressing quantities in terms of the ALE parametrization $\zeta^{\alpha}$.

In general, the material at any ALE coordinate $\zeta^{\alpha}$ will change over time, as there is in-plane flow relative to the mesh, and the mapping from convected to ALE coordinates is expressed as $\zeta^{\alpha}=\zeta^{\alpha}(\xi^{\beta},\,t)$. The surface position can be described equivalently with convected coordinates, surface-fixed coordinates, or the newly introduced ALE coordinates. To this end, any surface position $\bm{x}$ can be written as

[TABLE]

where in the last equality and from now on a ‘check’ accent is used to denote the position expressed in terms of the $\zeta^{\alpha}$ parametrization (see Fig. 2).

Under a change of surface parametrization from $\theta^{\alpha}$ to $\zeta^{\alpha}$, the latter of which can be specified arbitrarily, our geometric description of the surface is modified—as is any quantity with a Greek index. However, quantities transform in such a way that any variable without an index is invariant to the change in parametrization. Thus, for example, $\bm{v}$ and $v$ are invariant quantities while $v^{\alpha}$ and $a_{\alpha\beta}$ are not. As we will see, the continuity equation and vector form of the equations of motion contain no free indices, and can be expressed in terms of any surface parametrization—which shows the utility of our differential geometric developments and notation.

In this work, quantities are expressed in terms of the $\zeta^{\alpha}$ parametrization by placing a ‘check’ accent over every Greek index, where checked and unchecked indices take the same value. For example, $\bm{a}_{\check{\alpha}}:=\partial\check{\bm{x}}/\partial\zeta^{\alpha}=\check{\bm{x}}_{,{\check{\alpha}}}$ are the new in-plane basis vectors, $a_{{\check{\alpha}}{\check{\beta}}}:=\bm{a}_{\check{\alpha}}\cdot\bm{a}_{\check{\beta}}$ are the new metric components, and $b_{{\check{\alpha}}{\check{\beta}}}:=\bm{n}\cdot\check{\bm{x}}_{,{\check{\alpha}}{\check{\beta}}}$ are the new curvature components. In this manner, the ALE parametrization is used throughout the rest of this manuscript. We used a similar technique in our previous Lagrangian surface description [75], where we transitioned from the $\theta^{\alpha}$ to the $\xi^{\alpha}$ parametrization. As all quantities in Ref. [75] were written in the $\xi^{\alpha}$ representation, all Greek indices should be interpreted as having a ‘hat’ accent to be consistent with our notation.

2.4 ALE Kinematics

While the material velocity $\bm{v}$ is an invariant quantity, it is expressed differently for different surface parametrizations. For example, $\bm{v}:=[\partial\hat{\bm{x}}/\partial t]\rvert_{\xi^{\beta}}$ for convected coordinates and $\bm{v}=v^{\alpha}\bm{a}_{\alpha}+v\bm{n}$ for surface-fixed coordinates (see Sec. 2.2). In this section, we characterize the kinematics of an arbitrarily curved and deforming surface when the surface is parametrized by the ALE coordinates $\zeta^{\alpha}$. Our developments mirror those of Sec. 2.2.

Using the mapping $\zeta^{\alpha}=\zeta^{\alpha}(\xi^{\beta},\,t)$, any surface position can be written as

[TABLE]

which is analogous to Eq. (1). The velocity of a point, $\bm{v}=[\partial\hat{\bm{x}}/\partial t]\rvert_{\xi^{\beta}}$, can be expressed in the $\zeta^{\alpha}$ representation as

[TABLE]

where in the second equality we substituted $\bm{a}_{\check{\alpha}}=\bm{x}_{,{\check{\alpha}}}$. The last term in Eq. (5) describes how the position of a mesh point changes in time, which we denote the mesh velocity $\bm{v}^{\mathrm{m}}$:

[TABLE]

In Eq. (6), the notation $(\,~{}\,)^{\prime}:=[\partial(\,~{}\,)/\partial t]\rvert_{\zeta^{\alpha}}$ indicates how a quantity at a mesh point evolves in time. The partial time derivative $(\,~{}\,)^{\prime}$ in the ALE parametrization is analogous to the partial time derivative $(\,~{}\,)_{,t}$ in the surface-fixed parametrization, as both describe how a quantity changes at a fixed value of the appropriate coordinates. Note the mesh velocity $\bm{v}^{\mathrm{m}}$ need not be orthogonal to the surface, in contrast to its surface-fixed analog $v\bm{n}=[\partial\bm{x}/\partial t]\rvert_{\theta^{\alpha}}$ (see Sec. 2.2). Defining

[TABLE]

and using Eq. (6), we express Eq. (5) as

[TABLE]

Equation (8) indicates that $\bm{c}$ is the relative velocity between the material and the mesh. Furthermore, since $\bm{c}\cdot\bm{n}=c^{\check{\alpha}}\bm{a}_{\check{\alpha}}\cdot\bm{n}$ and $\bm{a}_{\check{\alpha}}\cdot\bm{n}=0$, Eq. (7)2 shows that the relative velocity lies entirely in the tangent plane to the surface, as shown schematically in Fig. 1(c).

The material time derivative of any invariant quantity can be expressed in the $\zeta^{\alpha}$ representation as

[TABLE]

where the relation $c^{\check{\alpha}}=(\bm{v}-\bm{v}^{\mathrm{m}})\cdot\bm{a}^{\check{\alpha}}$, obtained from Eqs. (7)2 and (8), is used in the second equality. The acceleration of a point $\bm{x}$ is calculated using Eq. (9) as

[TABLE]

Finally, in our simulations, the mesh velocity $\bm{v}^{\mathrm{m}}$ is treated as a fundamental unknown. The mesh position is calculated by integrating the mesh velocity over time, formally written as

[TABLE]

where $\check{\bm{x}}(\zeta^{\alpha},\,t_{0})$ is the initial mesh position at time $t_{0}$.

2.5 Mesh Velocity Equations

With the introduction of three new unknowns, namely the three components of the mesh velocity $\bm{v}^{\mathrm{m}}$, three additional governing equations are required for the problem to be mathematically well-posed. One equation is found by taking the dot product of Eq. (8) with the normal vector $\bm{n}$ and recognizing $\bm{c}\cdot\bm{n}=0$, yielding

[TABLE]

which ensures the mesh and the surface always overlap with one another as the surface deforms. In this sense, all schemes considered are out-of-plane Lagrangian.

The remaining two equations required to close the problem come from specifying the relative velocity $\bm{c}$, or equivalently specifying the relationship between $\bm{v}^{\mathrm{m}}\cdot\bm{a}_{\check{\alpha}}$ and $\bm{v}\cdot\bm{a}_{\check{\alpha}}$. There are no restrictions on how $\bm{v}^{\mathrm{m}}\cdot\bm{a}_{\check{\alpha}}$ and $\bm{v}\cdot\bm{a}_{\check{\alpha}}$ are related, so one can specify their relationship arbitrarily. If we were to choose $\bm{v}^{\mathrm{m}}\cdot\bm{a}_{\check{\alpha}}=\bm{v}\cdot\bm{a}_{\check{\alpha}}$, for example, we implicitly set $\zeta^{\alpha}=\xi^{\alpha}$ and therefore recover a Lagrangian scheme in which the mesh velocity and material velocity coincide. If, on the other hand, we choose $\bm{v}^{\mathrm{m}}\cdot\bm{a}_{\check{\alpha}}=0$, we implicitly set $\zeta^{\alpha}=\theta^{\alpha}$ and recover an in-plane Eulerian scheme in which the mesh moves only in the direction normal to the surface. The theoretical developments of this section allow us to specify $\bm{v}^{\mathrm{m}}\cdot\bm{a}^{\check{\alpha}}$, or equivalently $\zeta^{\alpha}$, arbitrarily as is best-suited to solve the problem at hand (see Fig. 1(c)). This flexibility is analogous to that of a Cartesian ALE formulation [51], and for this reason we name our scheme ‘arbitrary Lagrangian–Eulerian.’

For the majority of this manuscript we consider an out-of-plane Lagrangian, in-plane Eulerian scheme, for which the mesh velocity satisfies

[TABLE]

which from now on is called a Lagrangian–Eulerian (LE) scheme. As shown in Sec. 6, the LE implementation can be easily modified to produce a pure Lagrangian scheme. The analysis of more general mesh velocity descriptions is left to a future study, as care must be taken to avoid well-known ALE issues arising in the discretized equations—such as the violation of the geometric conservation law [78, 79]. In the LE case, such issues do not arise, and one can condense the three constraints on the mesh velocity into a single vector equation, given by

[TABLE]

Equation (14) provides the three equations necessary to resolve the LE mesh motion, and concludes our theoretical ALE surface description.

3 Fluid Film Equations

The equations governing an arbitrarily curved and deforming fluid film can be obtained in the form presented below by starting with the lipid membrane equations of Ref. [46], obtained within the framework of irreversible thermodynamics, and setting the bending moduli to zero. Here and from now on, all equations are written in terms of the ALE parametrization by placing ‘check’ accents over all Greek indices, as described in the previous section.

3.1 Strong Formulation

The continuity equation for an arbitrarily curved, incompressible two-dimensional material is given by

[TABLE]

Equation (15) is also called the incompressibility constraint and is enforced by the Lagrange multiplier $\lambda=\lambda(\zeta^{\alpha},t)$—which is physically the surface tension of the material. The local form of the linear momentum balance is given by

[TABLE]

where $\rho$ is the density, $\dot{\bm{v}}$ is the acceleration, $\bm{b}$ are the body forces, and $\bm{T}^{\check{\alpha}}$ are the stress vectors, namely, the boundary tractions across curves of constant $\zeta^{\alpha}$. For fluid films, we find $\bm{T}^{\check{\alpha}}=\sigma^{{\check{\alpha}}{\check{\beta}}}\bm{a}_{\check{\beta}}$, where $\sigma^{{\check{\alpha}}{\check{\beta}}}$ are the in-plane stress components given by (see Appendix A.1 and Ref. [46] for more details)

[TABLE]

In Eq. (17), $\lambda$ is the surface tension enforcing areal incompressibility (15) and $\pi^{{\check{\alpha}}{\check{\beta}}}$ are the in-plane viscous stresses. For an isotropic and incompressible fluid film, the viscous stresses are found to be

[TABLE]

where $\zeta$ is the two-dimensional shear viscosity with units of force$\cdot$time/length or equivalently mass/time.

To write the equations of motion in component form, we decompose the body force $\bm{b}$ in the $\{\bm{a}^{\check{\alpha}},\bm{n}\}$ basis as $\bm{b}=b_{\check{\alpha}}\,\bm{a}^{\check{\alpha}}+p\,\bm{n}$, where $b_{\check{\alpha}}$ are the in-plane covariant components and $p$ is the pressure drop across the surface. The in-plane equations of motion are given by

[TABLE]

as shown in Appendix A.2. The left-hand side of Eq. (19) contains the inertial terms, while the right-hand side consists of the body forces, surface tension gradient, and divergence of the viscous stresses, respectively. The viscous forces clearly show the coupling between surface curvature and the in-plane and out-of-plane velocity components.

The out-of-plane equation of motion, also called the shape equation, is found to be (see Appendix A.2)

[TABLE]

Equation (20) is an extension of the Young–Laplace equation to fluid films with nonzero velocity. Indeed, by setting $\bm{v}=\bm{0}$, Eq. (20) simplifies to the familiar Young–Laplace equation, $p+2\lambda H=0$. The presence of the in-plane fluid viscosity $\zeta$ and in-plane velocity components $v_{\check{\alpha}}$ in Eq. (20) leads to nontrivial coupling between the in-plane and out-of-plane equations when the surface is curved, i.e. when $b_{{\check{\alpha}}{\check{\beta}}}\neq 0$.

In our ALE formulation, the surface shape is evolved with the mesh velocity, rather than the material velocity, according to Eq. (11). As such, there are seven unknowns to solve for: three components of the material velocity $\bm{v}$, three components of the mesh velocity $\bm{v}^{\mathrm{m}}$, and the surface tension $\lambda$. The corresponding equations are the three components of the equation of motion (19) and (20), the three components of the mesh equation (14), and the incompressibility constraint (15). These seven equations constitute the strong formulation of the problem.

3.2 Boundary Conditions

The decomposition of the equations of motion into in-plane and out-of-plane components allows us to determine possible boundary conditions. The first term in parenthesis on the right-hand side of Eq. (19), $a^{{\check{\beta}}{\check{\gamma}}}\,v_{{\check{\alpha}};{\check{\beta}}{\check{\gamma}}}$, contains two derivatives of the in-plane velocity components $v_{\check{\alpha}}$. As no higher derivative of $v_{\check{\alpha}}$ appears, we accordingly specify either Dirichlet or Neumann boundary conditions at every point on the patch boundary. In analogy to the boundary conditions of a bulk fluid in three dimensions, Dirichlet boundary conditions specify the in-plane velocity component $v_{\check{\alpha}}$, while Neumann boundary conditions specify the in-plane boundary tractions $\bm{T}=\sigma^{{\check{\alpha}}{\check{\beta}}}\,\nu_{\check{\alpha}}\,\bm{a}_{\check{\beta}}$, where $\nu_{\check{\alpha}}=\bm{\nu}\cdot\bm{a}_{\check{\alpha}}$ and $\bm{\nu}$ is the in-plane unit normal to the surface at its boundary [46]. Accordingly, the patch boundary ${\partial\!\mathscr{P}}$ is separated into a Dirichlet portion $\partial^{\text{v}}\!\mathscr{P}$ and a Neumann portion $\partial^{\text{t}}\!\mathscr{P}$, such that $\partial^{\text{v}}\!\mathscr{P}\cap\partial^{\text{t}}\!\mathscr{P}=\varnothing$ and $\overline{\partial^{\text{v}}\!\mathscr{P}\cup\partial^{\text{t}}\!\mathscr{P}}={\partial\!\mathscr{P}}$. In this manuscript, only traction-free boundary conditions are considered, for which $\bm{T}=\bm{0}$ on $\partial^{\text{t}}\!\mathscr{P}$. General traction boundary conditions will be considered in Part II.

We next consider the shape equation (20), which describes the out-of-plane behavior of the fluid film and also provides an evolution equation for the position $\bm{x}$ through the normal velocity $v$. The shape equation contains two spatial derivatives of the position through the curvature components $b^{{\check{\alpha}}{\check{\beta}}}$. Consequently, at each point on the boundary we specify either the normal velocity or its in-plane gradient along the $\bm{\nu}$ direction, perpendicular to the surface boundary. These boundary conditions are independent of the in-plane boundary conditions, and in this manuscript we always specify the normal velocity of the patch boundary.

Our boundary conditions on the equations of motion are succinctly written as

[TABLE]

where $\bar{\bm{v}}$ is a known velocity on the boundary and $\bar{v}$ is a known normal component of the velocity on the boundary. As the equations governing the mesh velocity (14) are algebraic equations which do not contain any derivatives, we do not specify any boundary conditions for the mesh velocity $\bm{v}^{\mathrm{m}}$.

4 Finite Element Formulation for Fluid Films

In this section, we determine the weak formulation of the governing equations presented in Sec. 3.1, subject to the boundary conditions in Sec. 3.2. The weak form is modified to remove numerical inf–sup instabilities arising from the incompressibility of the fluid film, with a method inspired by Dohrmann and Bochev [76]. We provide the function spaces in which the solution to the weak formulation resides, and then use the standard tools of finite element analysis to find an approximate numerical solution to the weak form of the governing equations.

While inertial terms are included in the strong and weak formulations for completeness, they are not included in our simulations of arbitrarily curved and deforming fluid films—as they are negligible in many physical problems of interest. Despite the absence of inertia, the equations of motion are nonlinear due to the many terms involving the surface geometry. Furthermore, time derivatives still remain in the problem because the rate of change of the surface position is contained in the mesh velocity. The only simulation including inertial terms is the lid-driven cavity benchmark problem, in which the mesh is constrained to be fixed and inertia is included purely to demonstrate the validity of our numerical method.

4.1 Weak Formulation

Here we derive the weak formulation of the strong form equations provided in Sec. 3.1. Let $\mathcal{V}$ be the space of functions for the material velocity $\bm{v}$ and mesh velocity $\bm{v}^{\mathrm{m}}$, and let $\Lambda$ be the space of functions for the surface tension $\lambda$. We consider the arbitrary variations $\delta\bm{v}\in\mathcal{V}_{0}$, $\delta\bm{v}^{\mathrm{m}}\in\mathcal{V}$, and $\delta\lambda\in\Lambda$, where all elements of $\mathcal{V}_{0}\subset\mathcal{V}$ vanish on the Dirichlet portion of the boundary. The variations are contracted with the appropriate strong form equations and integrated over the fluid surface to yield the weak formulation of the problem.

In our previous work [46], we presented surface integrals as being of the general form $\int_{\mathscr{P}}(\ldots)~{}\mathrm{d}a$, where $\mathrm{d}a$ is a differential areal element of the patch $\mathscr{P}$. While such a description is theoretically sound, it is not amenable to numerical integration. We define $\Omega:=\{(\zeta^{1},\,\zeta^{2})\}$ to be the space of all ALE coordinates $\zeta^{\alpha}$, shown in blue in Fig. 2, and map areal integrals to $\Omega$ according to

[TABLE]

In Eq. (22), $\mathrm{d}\Omega:=\mathrm{d}\zeta^{1}\,\mathrm{d}\zeta^{2}$ is a differential element of $\Omega$ and $J^{\mathrm{m}}:=\sqrt{\rule{0.0pt}{5.03749pt}}\overline{\det\,a_{{\check{\alpha}}{\check{\beta}}}}$ is the Jacobian determinant of the mapping $\check{\bm{x}}(\zeta^{\alpha}):\Omega\rightarrow\mathscr{P}$. In a similar way, integrals over the patch boundary ${\partial\!\mathscr{P}}$ are mapped to integrals over the parametric domain boundary $\Gamma:=\partial\Omega$ according to

[TABLE]

where $\mathrm{d}s$ is a differential line element of the patch boundary ${\partial\!\mathscr{P}}$, $\mathrm{d}\Gamma$ is a differential line element on the parametric domain boundary, and $J^{\mathrm{m}}_{\Gamma}$ is the Jacobian determinant of the mapping $\check{\bm{x}}_{\mathrm{b}}(\zeta^{\alpha}):\Gamma\rightarrow{\partial\!\mathscr{P}}$—in which $\check{\bm{x}}_{\mathrm{b}}$ refers to the position of the patch boundary. The Dirichlet and Neumann portions of $\Gamma$, denoted $\Gamma^{\text{v}}$ and $\Gamma^{\text{t}}$, map to the patch boundaries $\partial^{\text{v}}\!\mathscr{P}$ and $\partial^{\text{t}}\!\mathscr{P}$, respectively.

To obtain the weak formulation, we begin by contracting the equations of motion (16) with an arbitrary velocity variation $\delta\bm{v}\in\mathcal{V}_{0}$ and integrating over the patch $\mathscr{P}$ to obtain

[TABLE]

where all integrals are mapped to $\Omega$ through Eq. (22). Note that in our previous Lagrangian formulation [74, 75], we solved for the material position $\bm{x}$ and the weak form contained an arbitrary position variation $\delta\bm{x}$. In this case, however, the in-plane fluidity necessitates the velocity $\bm{v}$ to be the fundamental unknown, such that the weak form is calculated with an arbitrary velocity variation $\delta\bm{v}$. Applying the surface divergence theorem to the third term on the left-hand side in Eq. (24) yields

[TABLE]

and the boundary integral is simplified by recognizing the velocity variation $\delta\bm{v}=\bm{0}$ on $\Gamma^{\text{v}}$ by construction and $\bm{T}=\bm{0}$ on $\Gamma^{\text{t}}$ (21)2, such that

[TABLE]

Substituting Eq. (26) into Eq. (24) and recognizing $\bm{T}^{\check{\alpha}}=\sigma^{{\check{\alpha}}{\check{\beta}}}\bm{a}_{\check{\beta}}\,$, with $\sigma^{{\check{\alpha}}{\check{\beta}}}$ given by Eq. (17), we obtain

[TABLE]

The weak form of the mesh equation is found by contracting Eq. (14) with an arbitrary mesh velocity variation $\delta\bm{v}^{\mathrm{m}}$, multiplying by a constant $\alpha^{\mathrm{m}}$, and integrating over the surface, which yields

[TABLE]

The factor $\alpha^{\mathrm{m}}$ is introduced such that $\mathcal{G}_{\mathrm{m}}$ has units of power, in agreement with $\mathcal{G}_{\mathrm{v}}$ (27). In our simulations, $\alpha^{\mathrm{m}}$ is set to unity. While Eq. (28) corresponds to the LE mesh equation (14), a different mesh equation can be used in our ALE formulation by setting $\mathcal{G}_{\mathrm{m}}$ to be the corresponding weak form expression.

Finally, as $\lambda$ is the Lagrange multiplier corresponding to the incompressibility constraint (15), we multiply Eq. (15) by an arbitrary variation $\delta\lambda\in\Lambda$ and integrate over the patch to obtain

[TABLE]

Eqs. (27)–(29) are the weak forms corresponding to the strong forms given respectively by Eqs. (16), (14), and (15). By summing them together and introducing a shorthand for the vector of unknowns, $\bm{u}$, as $\bm{u}^{\mathrm{T}}:=(\bm{v}^{\mathrm{T}},\,(\bm{v}^{\mathrm{m}})^{\mathrm{T}},\,\lambda)^{\mathrm{T}}$, its variation $\delta\bm{u}$ as $\delta\bm{u}^{\mathrm{T}}:=(\delta\bm{v}^{\mathrm{T}},\,(\delta\bm{v}^{\mathrm{m}})^{\mathrm{T}},\,\delta\lambda)^{\mathrm{T}}$, and the space of arbitrary variations $\mathcal{U}_{0}$ as $\mathcal{U}_{0}:=\mathcal{V}_{0}\times\mathcal{V}\times\Lambda$, we obtain the overall weak formulation, given by

[TABLE]

Note the weak form is nonlinear due to the out-of-plane deformations of the fluid film, as well as the inertial terms. Introducing $\mathcal{G}$ as the direct Galerkin expression [80] corresponding to the left-hand side of Eq. (30), the weak form can be compactly written as

[TABLE]

4.2 Solution Spaces

With the weak formulation (30), in what follows, we define the infinite-dimensional spaces in which the surface tension $\lambda$, velocity $\bm{v}$, and mesh velocity $\bm{v}^{\mathrm{m}}$ reside.

Surface Tension Solution Space

The surface tension $\lambda$ enters the weak form (30) without any gradients, and so we require $\lambda$ only to be square-integrable. We define the space of all possible surface tension fields, $\Lambda$, as the space of square-integrable functions on the parametric domain $\Omega$, denoted $L^{2}(\Omega)$. Thus, $\Lambda$ is given by

[TABLE]

Velocity Solution Space

The weak formulation (30) contains a gradient of both the velocity variation $\delta\bm{v}$ and the material velocity $\bm{v}$, the latter of which is contained in the viscous stresses $\pi^{{\check{\alpha}}{\check{\beta}}}$ (18). The velocity and velocity variation are both elements of the space of functions $\mathcal{V}$, and so elements of $\mathcal{V}$ are required to be square-integrable and have square-integrable gradients in order for the weak formulation to remain bounded. Furthermore, as the lipid membrane weak form requires the second derivatives of elements of $\mathcal{V}$ to be square-integrable, we define the space of velocities $\mathcal{V}$ as

[TABLE]

Each of the three Cartesian components of the velocity lies in $H^{2}(\Omega)$: the Sobolev space of order two on $\Omega$, defined as

[TABLE]

We also define $\mathcal{V}_{0}$ as the space of functions in $\mathcal{V}$ which vanish on $\Gamma^{\text{v}}$, written as

[TABLE]

such that $\delta\bm{v}\in\mathcal{V}_{0}$.

Mesh Velocity Solution Space

The weak formulation (30) contains terms involving the gradient of both the velocity variation $\delta\bm{v}$ and the mesh velocity $\bm{v}^{\mathrm{m}}$. While the mesh velocity gradient is not easily recognized in Eq. (30), it is found once the weak form is linearized and discretized (see Eq. (C.25)). As a result, we require $\bm{v}^{\mathrm{m}}\in\mathcal{V}$ in order for the weak form to remain bounded.

4.3 Finite-Dimensional Subspaces

We now choose the finite-dimensional subspaces in which we seek $\bm{v}_{h}$, $\bm{v}^{\mathrm{m}}_{h}$, and $\lambda_{h}$, which are approximations of the true solutions $\bm{v}$, $\bm{v}^{\mathrm{m}}$, and $\lambda$, respectively. The approximate surface position $\bm{x}_{h}$ is chosen to lie in the same subspace as $\bm{v}_{h}$, as is standard in isoparametric finite element methods [81]. To this end, we discretize $\Omega$ into ne (number of elements) non-overlapping elements $\{\Omega^{1},\,\Omega^{2},\,\ldots,\,\Omega^{\texttt{ne}}\}$, such that $\overline{\Omega}=\overline{\cup_{e=1}^{\,\,\texttt{ne}}\,{\Omega^{e}}}$ and $\Omega^{j}\cap\Omega^{k}=\varnothing$ for $j\neq k$. In all cases considered the parametric domain is discretized with a rectangular grid, as required by our choice of basis functions (see Sec. 4.5), such that all elements have the same dimensions—which is denoted $h$. The partitioning of the parametric domain naturally leads to finite-dimensional subspaces in which functions are polynomials over single elements and have certain continuity requirements across element boundaries.

The finite-dimensional subspace of velocities, $\mathcal{V}_{h}\subset\mathcal{V}$, is defined as

[TABLE]

where $C^{m}(\Omega)$ denotes the space of scalar functions on $\Omega$ with at least $m$ continuous derivatives, and $\mathbb{Q}_{n}({\Omega^{e}})$ is the space of bi-polynomial functions of order $n$ on the parametric element ${\Omega^{e}}$. Accordingly, $\mathcal{V}_{h}$ is the space of piecewise bi-quadratic functions with continuous first derivatives over the entire domain $\Omega$. While in the present formulation for fluid films, first derivatives need not be continuous, they are required to be continuous when modeling more complex systems which resist bending, such as lipid membranes [75] and viscoelastic Kirchhoff–Love shells [82]. We define the subspace $\mathcal{V}_{0,h}\subset\mathcal{V}_{0}$ to be the space of functions in $\mathcal{V}_{h}$ which are also zero on $\Gamma^{\text{v}}$, formally written as

[TABLE]

In the finite element analysis of bulk fluids, it is well-known that choosing the Lagrange multiplier space to be of the same polynomial order as the velocity leads to an unstable matrix equation, an issue resulting from the inf–sup condition, also called the Ladyzhenskaya–Babuška–Brezzi (LBB) condition [83, 84, 85]. We refer the reader to Ref. [86] for an excellent analysis of this numerical instability, which is often avoided in practice by choosing the Lagrange multiplier basis functions to be one polynomial order lower than the velocity basis functions. As we chose our velocities to be piecewise bi-quadratic functions on $\Omega$ (36), the Lagrange multiplier subspace is accordingly chosen to be continuous, piecewise bi-linear functions on $\Omega$, written as

[TABLE]

Even with this choice of basis functions, our scheme is LBB-unstable. As a result, we invoke a projection method devised by Dohrmann and Bochev [76], as described below, to further stabilize our numerical method.

4.4 Inf–Sup Stabilization

In modeling fluids with finite element methods, it is well-known that LBB errors arise when the velocity and surface tension solution spaces are identical. We initially used bi-quadratic velocities and bi-linear surface tensions (32)–(34), which were successfully used in our previous work [74, 75]. However, in the present study our numerical scheme exhibited LBB instabilities (see Appendix B.6.1). Further inspection showed our past work may have unknowingly avoided such instabilities by prescribing the surface tension along the entire boundary. In the spirit of developing a completely general finite element formulation, we seek to modify the numerical method presented thus far to remove the LBB instability.

We note there are many techniques to overcome the LBB instability, for example (i) using lower-order shape functions for the Lagrange multipliers [87], (ii) reduced and selective integration of the Lagrange multiplier equations [88, 89], (iii) stabilization methods [90], and (iv) macroelement approaches [91, 92]. All of these methods are valid in different cases. In this section, we describe how our weak form is modified by a technique developed by Dohrmann and Bochev [76] to remove LBB instabilities given general polynomial spaces for the velocity and surface tension; the analysis of other methods is left to a future study. The main idea underlying the Dohrmann–Bochev method is to locally project the surface tension onto a space of discontinuous, piecewise linear functions, and to energetically penalize the difference between the projected and unprojected surface tensions. The Dohrmann–Bochev method is thus based on two equations: one which projects the surface tension, and another which penalizes surface tension deviations from the projected space in a manner suitable for finite element analysis. A thorough description of the Dohrmann–Bochev method is provided in Ref. [81], and we follow their notation in this manuscript. The details of our numerical implementation can be found in Appendix B.6.

4.4.1 Theory

We begin by specifying the space of piecewise linear, discontinuous basis functions $\breve{\Lambda}$ onto which the surface tension is projected, given by

[TABLE]

where $\mathbb{P}_{n}({\Omega^{e}})$ is the space of polynomial functions of order $n$ on the parametric element ${\Omega^{e}}$. Note that while piecewise bi-linear functions can be continuous on quadrilateral elements (38), piecewise linear functions cannot be. Accordingly, the space $\breve{\Lambda}$ is discontinuous and over a single element ${\Omega^{e}}$: \breve{\Lambda}\big{\rvert}_{\Omega^{e}}=\mathbb{P}_{1}({\Omega^{e}}).

We next introduce the projection of the surface tension and its arbitrary variation, denoted $\breve{\lambda}$ and $\delta\breve{\lambda}$, respectively, such that $\breve{\lambda},\,\delta\breve{\lambda}\in\breve{\Lambda}$. The $L^{2}$-projection of the surface tension, $\breve{\lambda}$, is defined by

[TABLE]

As $\delta\breve{\lambda}$ belongs to a space of linear functions which are discontinuous across elements, Eq. (40) can be considered separately for individual elements ${\Omega^{e}}$, and is equivalently expressed as

[TABLE]

Deviations between $\breve{\lambda}$ and $\lambda$ are penalized in the weak form by subtracting the term

[TABLE]

from the left-hand side of Eq. (30). In Eq. (42), $\zeta$ is the two-dimensional fluid shear viscosity and $\alpha^{\mathrm{DB}}$ is a computational parameter having units of $J^{\mathrm{m}}$ which, as in Ref. [76], is chosen to be unity. The units of $\alpha^{\mathrm{DB}}$ ensure that Eq. (42) is dimensionally consistent with the other terms in the weak form (30). The weak formulation is now given by

[TABLE]

and the direct Galerkin expression $\mathcal{G}$ in Eq. (31) is redefined such that

[TABLE]

We provide the details of our Dohrmann–Bochev implementation in Appendix B.6, within which Appendix B.6.1 demonstrates the success of this method.

4.5 Summary of Numerical Solution Method

At this point, we seek an approximate solution to the weak formulation, as presented in Eqs. (43) and (44). To this end, we discretize the fundamental unknowns and their arbitrary variations. We then obtain the residual equations, discretize them temporally, and solve the resulting equations via Newton–Raphson iteration. Our techniques are briefly summarized here, however, extensive details of our numerical implementation are provided in Appendix B.

We first introduce the basis functions $\{N_{I}(\zeta^{\alpha})\}$ and Lagrange multiplier basis functions $\{\bar{N}_{J}(\zeta^{\alpha})\}$ such that $\mathcal{V}_{h}=(\textrm{span}\,\{N_{I}(\zeta^{\alpha})\})^{3}$, and $\Lambda_{h}=\textrm{span}\,\{\bar{N}_{J}(\zeta^{\alpha})\}$. The fundamental unknowns are discretized as

[TABLE]

where $[\mathbf{N}]$ and $[\bar{\mathbf{N}}]$ are shape function matrices and $[\mathbf{v}(t)]$, $[\mathbf{v}^{\mathrm{m}}(t)]$, and $[\bar{\bm{\lambda}}(t)]$ are the velocity, mesh velocity, and surface tension degree of freedom vectors, respectively. We introduce the shorthand $[\mathbf{u}(t)]=([\mathbf{v}(t)]^{\mathrm{T}},[\mathbf{v}^{\mathrm{m}}(t)]^{\mathrm{T}},[\bar{\bm{\lambda}}(t)]^{\mathrm{T}})^{\mathrm{T}}$ as the collection of all degrees of freedom in the discretized system. The arbitrary variations $\delta\bm{v}$, $\delta\bm{v}^{\mathrm{m}}$, and $\delta\lambda$ are discretized with the same basis functions as

[TABLE]

according to a Bubnov–Galerkin approximation, with $[\delta\mathbf{v}]$, $[\delta\mathbf{v}^{\mathrm{m}}]$, and $[\delta\bar{\bm{\lambda}}]$ collectively gathered into $[\delta\mathbf{u}]$. The weak formulation (44) can then be written as $[\delta\mathbf{u}]^{\mathrm{T}}[\mathbf{r}(t)]\,=\,0$ for all $[\delta\mathbf{u}]$, which is satisfied when the residual vector

[TABLE]

To solve Eq. (47) at a set of $N$ discrete times $\{t_{1},\,t_{2},\ldots,t_{N}\}$, we assume $[\mathbf{u}(t_{n})]$ is known and solve for $[\mathbf{u}(t_{n+1})]:=[\mathbf{u}(t_{n})]+[\Delta\mathbf{u}(t_{n+1})]$ according to the Newton–Raphson method. Again, a detailed description of our numerical procedure is provided in Appendix B.

In our implementation, the spaces $\mathcal{V}_{h}$ (36) and $\mathcal{V}_{0,h}$ (37) involve $C^{1}$-continuous basis functions. We maintain basis function continuity across elements by using uniform B-spline basis functions, which have the advantage of naturally enforcing arbitrary continuity requirements yet the complication of non-interpolatory basis function coefficients, and the requirement of a rectangular parametric mesh, as well as basis functions spreading over multiple elements. The method of using B-spline basis functions within an isoparametric finite element framework is in the spirit of Ref. [93], and detailed in Ref. [94]. We use the algorithms described in Ref. [95] to efficiently calculate the basis functions and their derivatives.

We have now concluded our discussion of the LE finite element formulation, and a high-level overview of our code structure can be found in Algorithm 1 of Appendix C.6.

5 Numerical Simulations

We now present several results from our LE finite element formulation to validate the robustness of the method and demonstrate its capabilities. The numerical implementation of the method is tested with problems of increasing complexity, and results are compared to known analytical solutions whenever possible. The first test cases involve fluid flows on flat planes, and once several cases are validated we move on to study fluid flows on fixed, curved surfaces. In the scenarios mentioned thus far, the mesh is constrained to remain stationary and the mesh velocity is not solved for. In our last example, the entire LE implementation is tested by modeling an initially cylindrical fluid film, which is allowed to deform over time. The simulations show that fluid films are unstable with respect to long wavelength perturbations, which may explain why bubbles are often observed and long, cylindrical fluid films are not. Namely, in any experimental system, we expect the latter to break up and form bubbles. A linear stability analysis is performed to calculate the critical length above which the cylinder becomes unstable, as well as the time scale of the instability. We show that the analytical predictions for the critical length and the time scale of the fastest growing unstable mode agree quantitatively with our simulation results. Moreover, we find both theoretically and computationally that cylindrical fluid films are stable to non-axisymmetric perturbations.

5.1 Fluid Flow on a Flat Plane

We first consider the simplest test case, fluid flowing on a flat plane, for which the surface tension can be equivalently thought of as the negative surface pressure and the governing equations simplify to the two-dimensional incompressible Navier–Stokes equations (see Appendix A.2.1). While the problems considered in this section are easily solved using standard Cartesian finite element methods, we solve them using our nonlinear isoparametric finite element framework, in which differential geometry is used to express the surface position and fundamental unknowns in terms of curvilinear coordinates. Thus, even these simple problems serve as important benchmarks for our LE finite element implementation.

In our simulation of fluid flow on a flat plane, the mesh is constrained to be fixed, such that $\bm{v}^{\mathrm{m}}=\bm{0}$, and we do not solve for the mesh velocity. Furthermore, as the mesh is flat, there is no motion in the out-of-plane direction and we solve only for the $x$- and $y$-components of the fluid velocity as well as the surface tension. The three corresponding strong form equations are the continuity equation (15) and the two in-plane Navier–Stokes equations (A.22)3. We simulate five scenarios with increasingly complex solutions: (i) a hydrostatic fluid, with zero velocity and linear surface tension; (ii) Couette flow between parallel plates, with linear velocity and zero surface tension; (iii) Couette flow between parallel plates under the influence of a quadratic body force, with linear velocity and cubic surface tension; (iv) Hagen–Poiseuille flow in a channel, with quadratic velocity and linear surface tension; and finally (v) the lid-driven cavity problem, for which no analytical solution is known. Only the lid-driven cavity result is discussed in the main text, and the validation of the first four cases against analytical solutions is left to Appendix D.

For the numerical results presented in this manuscript, we neglect inertial terms in all cases except the lid-driven cavity, for which inertial terms are evaluated using a backwards Euler temporal discretization. The contribution of inertial terms to the tangent matrix and residual vector, for the limited case of a fixed surface, is provided in Appendix C.3.1. Inertial terms are included in the lid-driven cavity problem only to further validate our numerical implementation with flows at moderate Reynolds numbers.

5.1.1 Lid-Driven Cavity

A schematic of the problem is shown in Fig. 3(a). Fluid in a square cavity with stationary walls is driven by a top lid which moves to the right at constant speed $V$. We solve for the flow field and surface tension in the cavity, which is taken to be a unit square in the $x$–$y$ plane. The boundary condition $\bm{v}=\bm{0}$ is imposed on the sides and bottom of the square domain, and $\bm{v}=V\bm{e}_{x}$ is set on the top edge. There is a choice in what velocity to specify at the top two corners of the domain where the stationary edges meet the moving top lid. At these locations we set $\bm{v}=\bm{0}$, as done in Ref. [96]. Furthermore, as only gradients of the surface tension enter the equations governing a flat, two-dimensional fluid, the surface tension $\lambda$ is indeterminate up to a constant. Consequently, $\lambda$ is specified to be zero at the center of the domain, located at $(x,\,y)=(0.5,\,0.5)$.

We first set the inertial terms to zero and solve for the Stokes flow result at Reynolds number $\mbox{{Re}}:=\rho\,VL/\zeta=0$, for which the $x$-velocity is plotted in Fig. 3(b). The flow field is moving towards the right at the top of the domain, and the presence of the right wall requires the fluid to be recirculated towards the left side in the bulk of the domain. We then include inertial terms and set $\mbox{{Re}}=100$ by setting $\rho=1$ and $\zeta=0.01$ (we already have $V=1$ and $L=1$). The corresponding $x$-velocity is plotted in Fig. 3(c). Relative to Fig. 3(b), Fig. 3(c) shows that at the top of the domain, inertia pulls the fluid to the right in the direction of the moving top lid. In Figs. 3(d) and 3(e), the $x$-velocity and surface tension are plotted at different Re across a vertical cross-section through the domain. Our results for the lid-driven cavity problem agree with those of other numerical studies [97, 98], as plotted in Fig. 3(d).

Our final task for the lid-driven cavity problem is to analyze how our simulations converge to the true solution as we increase the number of elements. We consider the case with no inertia, i.e. $\mbox{{Re}}=0$, and calculate the error as a function of the length (and width) of a single element, denoted $h$. As the exact solution is unknown, we treat the solution on the finest, $256\times 256$ mesh as the true solution. We calculate the $L^{2}$-error for the velocity, denoted $||\bm{v}-\bm{v}_{h}||_{0}$ and defined as

[TABLE]

where $\bm{v}$ is the true velocity and $\bm{v}_{h}$ is the approximate velocity found on coarser meshes with elements of size $h\times h$. The $L^{2}$-error is denoted with a subscript ‘0’ because $L^{2}(\Omega)=H^{0}(\Omega)$, i.e. the Sobolev space of order zero. A plot of $||\bm{v}-\bm{v}_{h}||_{0}$ as a function of $1/h$ is shown in Fig. 3(f), in which we see a linear scaling of the $L^{2}$-error with the mesh size $h$. The numerical treatment of the corner nodes may explain why our simulations do not converge faster: at each mesh size, we set $\bm{v}=\bm{0}$ at the corner nodes and $\bm{v}=V\,\bm{e}_{z}$ at the adjacent nodes on the top edge, however, as the number of elements increases, these two nodes move closer together such that we solve a slightly different problem at each mesh size.

Due to the moving top lid being in contact with the stationary walls at the top two corners of the domain, the lid-driven cavity solution is known to have large surface tension spikes at the corners [96]. For this reason, numerical studies generally use meshes which are very finely discretized at the corners. We find the surface tension spikes vary significantly with the mesh size, for example, on a $16\times 16$ mesh $\lambda$ is $\mathcal{O}(10^{2})$ at the corners while on a $128\times 128$ mesh it is $\mathcal{O}(10^{3})$. These spikes dominate the error calculation, which even on the $128\times 128$ mesh is $\mathcal{O}(1)$ relative to the solution on a $256\times 256$ mesh. The surface tension spikes are plotted in Appendix B.6.1, in which the success of the Dohrmann–Bochev method in removing inf–sup instabilities is also shown. As the surface tension error does not converge in this situation, in Appendix D.3 we show how it converges in the benchmark problem of Couette flow with a body force.

With the numerical results provided thus far, as well as those of Appendix D, we conclude our analysis of fluid flow on a fixed, flat plane. Our general ALE finite element framework, based on a curvilinear coordinate description via the machinery of differential geometry, has successfully reproduced the classical benchmarks of hydrostatic flow, Couette flow, Hagen–Poiseuille flow, and lid-driven cavity flow. We are therefore confident our finite element method can model arbitrary flows on flat surfaces.

5.2 Fluid Flow on a Stationary, Curved Surface

After testing our code on the simplest case of fluid flow on a flat plane, we turn to study fluid flows on stationary, curved surfaces. In such problems the mesh is fixed, and thus once again the mesh velocities are not solved for. However, a complication arises because for a fixed surface, the normal component of the material velocity $v=\bm{v}\cdot\bm{n}=0$. In practice, there are two ways we could handle this problem. First, we could posit that our arbitrary velocity variation $\delta\bm{v}=\delta v^{\check{\alpha}}\bm{a}_{\check{\alpha}}$, where $\bm{a}_{\check{\alpha}}$ is known for the fixed surface. However, such a method cannot be easily extended to model a deforming fluid film: when the surface is deforming, the velocity $\bm{v}$ is represented as $\bm{v}=v^{\check{\alpha}}\bm{a}_{\check{\alpha}}+v\bm{n}$, and in a fully implicit numerical scheme $v^{\check{\alpha}}$, $\bm{a}_{\check{\alpha}}$, $v$, and $\bm{n}$ are all unknown.

We employ a different approach by representing the velocity and its arbitrary variation in a Cartesian basis, including both the in-plane and shape equations in our weak formulation, and enforcing $v=\bm{v}\cdot\bm{n}=0$ with a Lagrange multiplier—which is interpreted physically as the pressure drop required to constrain the fluid to the surface. In this section, we describe how both the strong and weak formulations are modified by the normal pressure $p$ being an unknown Lagrange multiplier field. Our method is then tested by modeling fluid flow on a fixed, bulged cylinder. Numerical results are compared with analytical solutions, and we calculate how our numerical error decreases on mesh refinement.

5.2.1 Strong Form Modification

In satisfying the constraint $v=0$ with a Lagrange multiplier, we include the shape equation (20) in our description of the fluid film. In the shape equation (20), the pressure drop $p$ is an unknown Lagrange multiplier field, which at every point on the surface takes the requisite value to satisfy $v=0$. There are thus five unknowns: the three components of the velocity, the surface tension, and the normal pressure drop, and five corresponding equations: the continuity equation (15), the two in-plane equations (19), the shape equation (20), and the constraint $v=0$.

5.2.2 Weak Form Modification

With the pressure drop $p$ being a fundamental unknown, the arbitrary pressure variation $\delta p$ is expected to enter the weak formulation. To understand how $\delta p$ will appear, we take the variation of the virtual work associated with moving the fluid film in the normal direction and obtain

[TABLE]

The first term on the right-hand side of Eq. (49) is already contained in Eq. (43) through the body force term. Assuming no in-plane body forces ($b_{\check{\alpha}}=0$), treating the pressure $p$ as a fundamental unknown, again recognizing inertia is negligible, and removing mesh velocity degrees of freedom yields a modified weak form (c.f. Eq. (43))

[TABLE]

where $\mathcal{P}$ is the space of pressure solutions and arbitrary pressure variations. The weak form (50) contains no gradients of pressure and thus it is theoretically sound for us to choose $\mathcal{P}$ as the space of square-integrable functions on $\Omega$, i.e. $L^{2}(\Omega)$. However, we simplify our finite element analysis by using the same basis functions for the velocity and pressure. In accordance with Eqs. (33) and (34), we define

[TABLE]

The structure of Eq. (50) indicates the pressure variation $\delta p$ enforces the normal constraint $\bm{v}\cdot\bm{n}=0$, in the same way the surface tension variation $\delta\lambda$ enforces the incompressibility constraint $\bm{a}^{\check{\alpha}}\cdot\bm{v}_{,{\check{\alpha}}}=0$. With the weak formulation (50), an identical procedure to that of Sec. 4 is followed to linearize and discretize the equation, calculate the tangent matrix and residual vector, and then iteratively solve for the unknowns via Newton–Raphson iteration. The approximate pressure $p_{h}$ is an element of the finite-dimensional subspace $\mathcal{P}_{h}\subset\mathcal{P}$, which is chosen to be

[TABLE]

in accordance with the finite-dimensional space of velocities $\mathcal{V}_{h}$ (36). The pressure can then be expressed in terms of the same set of basis functions, $\{N_{I}(\zeta^{\alpha})\}$, used for the velocities. As mentioned previously, our choice of $\mathcal{P}$ (51) and $\mathcal{P}_{h}$ (52) is purely for convenience in our numerical implementation. The details of the modifications to our finite element implementation are provided in Appendix C.3.

5.2.3 Flow on a Fixed Cylinder with a Bulge

In our numerical implementation, we consider fluid flowing on a fixed, bulged cylinder, as shown in Figs. 4(a)–4(c). Our boundary conditions are shown schematically in Fig. 4(a), where constant inflow and outflow of the fluid is prescribed at the entrance and exit of the cylinder, respectively. The surface tension is specified at a single point on the boundary, as only gradients of $\lambda$ enter the in-plane equations. The bulge in the center leads to nontrivial velocity and surface tension profiles due to the coupling between curvature and fluid flow, and the symmetry of the surface shape allows us to determine the analytical solution (see Appendix A.2.3). The bulged cylinder thus serves as a useful benchmark problem for our numerical method, in the study of flows on curved yet fixed surfaces.

The position of the bulged, cylindrical surface is given by

[TABLE]

where $\theta$ and $z$ are the polar angle and axial distance, respectively, of a standard cylindrical coordinate system. The radius $r(z)$ is independent of $\theta$ because the surface is axisymmetric. Our choice of radial profile is given in Appendix A.2.3 and shown in Fig. 4(b). We define $h_{\theta}$ and $h_{z}$ as the width of an element ${\Omega^{e}}$ in the $\theta$ and $z$ directions, respectively, and denote a mesh with 16 elements in the $\theta$-direction and 32 elements in the $z$-direction, for example, as a $16\times 32$ mesh. A coarse $10\times 10$ mesh of the bulged cylinder is shown in Fig. 4(c). The surface tension, $z$-velocity, and normal pressure are calculated numerically on a $128\times 128$ mesh, and compared to their analytical counterparts in Figs. 4(d)–4(f), respectively. In all cases, our numerical results show excellent agreement with the analytical calculations.

We conclude our analysis of the fixed, bulged cylinder by calculating the $L^{2}$-error of the velocity, surface tension, and normal pressure. We also calculate the error between our numerical surface position $\bm{x}_{h}$ and the true position $\bm{x}$ given in Eq. (53), which is a function of how closely our basis functions can represent the analytical geometry. Figures 4(g)–4(i) show three different plots of the error: changing $h_{\theta}$ and $h_{z}$ together (4(g)), changing $h_{z}$ with fixed $h_{\theta}$ (4(h)), and changing $h_{\theta}$ with fixed $h_{z}$ (4(i)). Figures 4(g) and 4(h) are nearly identical, indicating $\theta$-refinement has little effect on the errors, an observation which is confirmed with Fig. 4(i). As the analytical solution is only a function of $z$, it is unsurprising that refining in $\theta$ has no significant effect on the errors. In particular, $\theta$-refinement only makes observable changes to our numerical surface representation when there are few elements in the $\theta$-direction. In Figs. 4(g) and 4(h), the errors in the position, velocity, and normal pressure converge quadratically on mesh refinement. However, the error in $\lambda$ is approximately constant. As discussed in Appendix A.2.3, a cylinder length $L=4,\!000$ is chosen, with radius changes over extended length scales, to avoid discontinuities in the normal pressure. For such cylinders, $r^{\prime}(z)\ll 1$ and $r^{\prime\prime}(z)\ll 1$, and therefore $\lambda(z)\approx\lambda_{0}=1$ everywhere (see Eq. (A.37))—irrespective of the shape of the bulge. Analytically, the difference between $\lambda$ and $\lambda_{0}$ is $\mathcal{O}(10^{-4})$ and occurs when the radius changes due to nonzero $r^{\prime}(z)$, $r^{\prime\prime}(z)$, and $z$-velocity (A.37). Numerically, errors in $\lambda$ are also $\mathcal{O}(10^{-4})$ and occur when the radius changes due to the difference between the numerical and analytical position $\bm{x}$; these errors persist along the length of the cylinder (Fig. 4(d)). As can be seen from Fig. 4(h), the errors in position $\bm{x}$ are $\gtrsim 10^{-4}$ for all simulated meshes, and therefore the error in $\lambda$ remains approximately constant at $10^{-4}$. We expect the $\lambda$ error to decrease as one further refines the mesh, thus reducing errors in the position. However, these simulations are prohibitively expensive in our current implementation and require parallelization of our code. Our error calculations, captured in Figs. 4(g)–4(i), conclude our analysis of fluid flows on fixed, curved surfaces.

5.3 Curved and Deforming Fluid Films

In this section, we demonstrate the capabilities of our full LE finite element implementation by studying the stability of cylindrical fluid films which have a constant pressure drop across their surface. It is already known that given a surface tension and pressure drop, spherical fluid films of radius $R$ satisfy the Young–Laplace equation $p=2\lambda/R$ and are stable [99]. The cylindrical fluid films under consideration, however, are found to be unstable to axisymmetric perturbations; the large, nontrivial shape changes resulting from the instability are studied here. Our numerical results are compared with analytical results from a linear stability analysis, and the two are found to be in excellent agreement. Moreover, we find the fluid film instability is mediated by the in-plane flow resulting from the initial shape perturbation, thus demonstrating the importance of our ALE framework in studying two-dimensional surfaces with in-plane flow. Our general LE implementation is also employed to show cylinders of any length are stable to non-axisymmetric perturbations, as confirmed analytically. We end by showing our results are independent of our choice of time step and mesh size, thus demonstrating the necessary convergence of numerical solutions.

5.3.1 Stability of Axisymmetric Perturbations

To begin, the position of an unperturbed cylinder of radius $R$ and length $L$ is given by

[TABLE]

with the polar angle $\theta\in[0,2\pi)$ and the axial length $z\in[0,L]$. For the boundary conditions $\bm{v}=\bm{0}$ on both edges of the cylinder, a valid base state solution is given by $\bm{v}=\bm{0}$ and $\lambda=\lambda_{0}$ everywhere, as shown in Appendix A.3. In this case the shape equation (20) simplifies to the Young–Laplace equation $\lambda_{0}=p\,R$, with $p$ being the constant pressure drop imposed across the fluid surface.

At time $t=0$, we perturb the radius of the stationary cylinder such that the initial position is given by

[TABLE]

where the small dimensionless parameter $\epsilon$ is set to $0.01$ in our simulations. A schematic of the initial surface shape is shown in Fig. 5(a). We evolve the fluid film from its initial state using our LE method and observe if it is stable or unstable with respect to the initial perturbation. Over the course of our simulation, we maintain a constant pressure drop across the fluid surface. This pressure drop enters as a body force $\rho\bm{b}=p\bm{n}$, and its contribution to the tangent matrix and residual vector is provided in Appendix C.4. We find the cylinder is stable to the initial perturbation when its length $L<2\pi R$ and unstable when $L>2\pi R$. In the latter case, the initial perturbation continues to grow and eventually reaches a configuration that has spherical bulbs on the two ends (Figs. 5(b)–5(f), supplemental movie E.1, video provided at youtu.be/FUx8fGXuzqY). These spherical shapes are believed to form because they are the stable surfaces compatible with our boundary conditions and incompressibility constraint. We proceed to use simple physical arguments to understand why the initial perturbation is unstable.

5.3.2 Physical Explanation of the Instability

The instability arises because our initial shape perturbation changes the mean curvature of the surface, which in turn changes the surface tension through the Young–Laplace equation. The resultant surface tension gradients then drive in-plane fluid flow, as can be seen from the in-plane equations (19). When $L>2\pi R$, fluid flows from the narrow region of the cylinder to the wide region (see Fig. 5(a)), resulting in an unstable film. However, when $L<2\pi R$, the surface flow is directed from the wide region to the narrow region, which causes the initial bulge to shrink over time such that the surface returns to its cylindrical configuration.

To understand this general idea in more detail, we begin with the Young–Laplace equation, written as $\lambda=-p/(2H)$. For an unperturbed cylinder, $H=-1/(2R)$. The initial perturbation (55) alters the mean curvature of the film by modifying both radii of curvature: it changes the radius of a circular cross-section of the cylinder, and the sinusoidal shape along $z$ introduces a nonzero radius of curvature in the axial direction. Analytically, the mean curvature $H$ of the initially perturbed shape (55) is calculated according to Sec. 2.1 as

[TABLE]

Consider the quantity in square brackets in Eq. (56), which consists of two terms. The first term comes from the change in the circular cross-section, while the second term, which contains a factor of $L^{-2}$, comes from the change in radius along the $z$-direction. According to Eq. (56), when $L>2\pi R$, the mean curvature $H$ becomes less negative where the cylinder bulges outwards and more negative where it bulges inwards. The Young–Laplace equation then requires $\lambda$ to become larger at the outward bulge and smaller at the inward bulge (see Fig. 5(b)), and the in-plane equations (19) indicate fluid flows from regions of low $\lambda$ to high $\lambda$. As a result, fluid flows along the surface tension gradient which in this case causes the instability to grow (see Fig. 5(f)). When $L<2\pi R$, on the other hand, the effect on $\lambda$ is reversed and the in-plane flow serves to stabilize the perturbed cylindrical shape.

5.3.3 Instability Time Scale

In addition to describing the instability with simple physical arguments, we performed a linear stability analysis and found the fluid film equations are indeed unstable to the initial perturbation (55) when $L>2\pi R$ (see Appendix A.3). Our analysis also revealed a theoretical time scale $\tau_{\mathrm{theo}}$ for the instability which, when $L>2\pi R$, is given by

[TABLE]

Dimensionally, the ratio $\zeta/\lambda_{0}$ is expected to set the time scale of the deforming fluid film, as found in a study of topological transitions in two-dimensional dry foams [100]. Equation (57) indicates $\tau_{\mathrm{theo}}\rightarrow 4\zeta/\lambda_{0}$ as $L\rightarrow\infty$ and $\tau_{\mathrm{theo}}\rightarrow\infty$ as $L\rightarrow 2\pi R^{+}$, as shown by the green curve in Fig. 6.

We also calculate the time scale $\tau_{\mathrm{sim}}$ from our full LE simulations as a function of length $L$ and compare it to the analytical prediction (57). Assuming the unstable perturbation (55) initially grows exponentially in time, at small times the surface shape is given by

[TABLE]

where $\tilde{r}:=\sin(2\pi z/L)$ is defined for notational convenience. We seek the time $t^{*}$ for which the initial instability has grown by the chosen multiplicative factor $\chi$, such that the surface position is given by

[TABLE]

If $\chi=2$, for example, we numerically measure the time $t^{*}$ when the initial perturbation doubles in size. By comparing Eqs. (58) and (59), the time scale $\tau_{\mathrm{sim}}$ can be numerically calculated as

[TABLE]

We emphasize that in Eq. (60), $\chi$ is a chosen marker for the growth of the instability and $t^{*}$ is numerically measured to calculate $\tau_{\mathrm{sim}}$. Figure 6 shows the excellent agreement between numerically calculated (60) and analytical (57) time scales, as a function of cylinder length, and demonstrates the simulations correctly predict the limiting time scale $\tau=4\zeta/\lambda_{0}$ as $L\rightarrow\infty$.

5.3.4 Stability of Non-axisymmetric Perturbations

The linear stability analysis (Appendix A.3) also considered non-axisymmetric modes, and found the resultant surface tension gradients always stabilized the cylindrical configuration. Thus, all non-axisymmetric modes are stable. Here we consider the non-axisymmetric perturbation

[TABLE]

where we redefine $\tilde{r}$ as $\tilde{r}:=\sin(3\theta)\,\sin(2\pi z/L)$ for notational convenience. The corresponding mean curvature is calculated as

[TABLE]

According to Eq. (62), if the cylinder bulges outwards ($\tilde{r}>0$) then $H$ is more negative and the surface tension decreases. On the other hand, if the cylinder bulges inwards ($\tilde{r}<0$) then $H$ is less negative and the surface tension increases. The resulting in-plane flow goes from the outward bulges to the inward bulges, and always stabilizes the film shape. We confirm the theoretical result with our full non-axisymmetric simulations, as shown in Fig. 7, and note all non-axisymmetric perturbations are stable (see Appendix A.3).

5.3.5 Time Step and Mesh Refinement

Finally, we study the convergence of our axisymmetric numerical results in three cases: refining $\Delta t$, refining $h_{z}$, and refining $h_{\theta}$. In each convergence study, all simulations are run until time $t=35$, at which point the initially cylindrical film has undergone significant deformation. Figure 5(f) shows the shape of the fluid film at $t=35$; Fig. 8 shows axial profiles of the radius and fluid $z$-velocity at different snapshots in time. In the latter, the fluid film shape and velocity are visually indistinguishable for different choices of the time step, as in each case $\Delta t$ is considerably less than the time scale $\tau\approx 6.6$ seconds (Eq. (57) with $L=10$).

The $L^{2}$-error is calculated for the fluid film position $\bm{x}$ and material velocity $\bm{v}$ at time $t=35$, with the true solution being approximated as the finest simulation run (further details are provided in Fig. 9). Refining the time step $\Delta t$ shows linear scaling in the position and material velocity, as shown in Fig. 9(a). Both the position and velocity are expected to scale linearly because we used a backward Euler temporal discretization, in which the position is not a fundamental unknown but rather calculated from the mesh velocity according to Eq. (11). The $L^{2}$-error scales quadratically for both $z$-refinement (Fig. 9(b)) and $\theta$-refinement (Fig. 9(c)). Given our implementation of $C^{1}$-continuous bi-quadratic velocities (33) and $C^{0}$-continuous bi-linear surface tensions (32) with uniform B-splines, the quadratic scaling on mesh refinement is also expected [94]. Thus, our LE simulation results demonstrate the anticipated convergence behavior upon both time step and mesh refinement.

6 Lagrangian Implementation

As previously mentioned, our ALE formulation accommodates different mesh velocity equations. In this section, a Lagrangian scheme is implemented to demonstrate the generality of our formulation. First, the lid-driven cavity problem is reexamined, which demonstrates how Lagrangian methods are unsuitable for problems involving shear flows. The Lagrangian implementation is then used to model the unstable cylindrical fluid film described in Sec. 5.3, and is shown to correctly capture the time scale associated with the instability. However, in the Lagrangian simulations, mesh nodes move due to in-plane fluid flow. As a consequence, at later times certain regions of the surface contain only a few nodes and are spatially unresolved. Thus, Lagrangian simulations have errors larger than those of their LE counterparts.

Our finite element framework easily accommodates a Lagrangian formulation, as a Lagrangian scheme is recovered when

[TABLE]

which replaces Eq. (14) as the strong form of the mesh equation. Rather than satisfying Eq. (63) weakly, we simply enforce the velocity and mesh velocity degrees of freedom to be equal in our discretized system of equations at every Newton–Raphson step. Our numerical implementation is described in Appendix B.7.

A Lagrangian formulation is unsuitable for even simple situations with in-plane velocity gradients. For example, a Lagrangian scheme fails to reproduce the lid-driven cavity result of Sec. 5.1.1 (supplemental movie E.2, video at youtu.be/FCoShaa_FhM). We thus analyze our Lagrangian implementation by simulating a perturbed cylindrical fluid film, and comparing the results to those of the LE simulations, found in Sec. 5.3. We first repeat the calculation of the deformation time scale (see Sec. 5.3.3), and find excellent agreement between theory and Lagrangian simulations (Fig. 10(a)). Next, the distribution of nodes at time $t=35$ are compared for both the Lagrangian and LE implementations. While the LE nodes maintain roughly equal spacing, the Lagrangian nodes move with the in-plane flow and become more concentrated in the spherical bulbs (Fig. 10(b)). Relative to the LE results, in Lagrangian simulations the central tubular region is spatially underresolved at later times (Figs. 10(d)–10(f), supplemental movie E.3, video at youtu.be/wnXuK6d3WTQ). Accordingly, Lagrangian simulations have comparatively larger errors upon $z$-refinement (Fig. 10(c), c.f. Fig. 9(b)).

7 Conclusions

In this paper, we developed an ALE theory and formulation for two-dimensional materials which are arbitrarily curved, may be deforming over time, and exhibit in-plane fluidity. Within the setting of differential geometry, we introduced a new surface parametrization to endow the surface with a mesh whose in-plane motion can be specified arbitrarily. In particular, the mesh need not convect with the in-plane material velocity, and thus our framework does not suffer from the limitations of a Lagrangian formulation. With the ALE theory, and the results of a previous irreversible thermodynamic study [46], we used the standard tools of finite element analysis to devise an isoparametric, fully implicit finite element method for two-dimensional curved and deforming fluid films. In particular, we (i) developed a weak formulation of the problem, (ii) spatially discretized it with a standard Bubnov–Galerkin approximation, (iii) temporally discretized the residual equations with the backward Euler method, and (iv) linearized the resulting system of equations using the Newton–Raphson method. We also highlighted our use and implementation of the Dohrmann–Bochev method, which overcame the LBB instability. Finally, we provided several numerical examples to showcase the merits of our method. We showed how our formulation is easily adapted to solve for fixed-surface flows, with the added benefit of determining the pressure drop required to constrain the fluid to the surface. We also used both LE and pure Lagrangian numerical implementations to model cylindrical films which can deform over time. We found such films are unstable above a critical length, and demonstrated that our numerical results are in agreement with the predictions of a linear stability analysis.

Our main motivation behind this work is to develop a robust, theoretically sound, and sufficiently general method for modeling lipid membranes. To this end, we intend to implement a mesh velocity scheme which is not simply in-plane Eulerian (14) or Lagrangian (63) in our future work—which, in principle, should overcome the squeezing of the surface-fixed coordinate nodes depicted in Fig. 1(b). Furthermore, as fluid films are a computational precursor to lipid membranes, we plan to extend our implementation to model single-component lipid membranes by adding bending contributions to the weak formulation. This analysis will be presented as Part II of the current work [1]. With a numerical method to simulate single-component lipid membranes, our framework can be extended to describe various phenomena governing biological membranes—such as the coupling between the bulk fluid surrounding the membrane and intramembrane lipid flows [63] and in-plane species diffusion [62]. Moreover, one could analyze multi-component systems, which consist of different lipid and transmembrane proteins. Recent microscopic studies demonstrate how protein–lipid and protein–protein interactions play an important role in membrane bending [101, 102, 103, 104] and lipid domain formation [105]; experiments also demonstrate the coupling between lipid phase separation and membrane bending [106]. Modeling such phenomena [33] would be a natural extension of our work.

Acknowledgements

We are indebted to Prof. Panos Papadopoulos for discussions on LBB conditions, and for bringing the Dohrmann–Bochev method to our attention. We thank Michelle Liu for assistance in visualizing our numerical results, Joël Tchoufag for many useful discussions, and Prof. Howard Stone for helpful comments and suggestions. A.S. would also like to thank Prof. Phil Colella for his C++ makefile and suggested coding best practices.

A.S. is supported by the Computational Science Graduate Fellowship from the U.S. Department of Energy, as well as U.C. Berkeley. Y.A.D.O. would like to thank the support of U.C. Berkeley. R.A.S. acknowledges the support of the German Research Foundation (DFG) through Grant No. GSC 111, as well as the Aachen–California Network of Academic Exchange. K.K.M. acknowledges the support of U.C. Berkeley and U.S. Department of Energy Contract No. DE-AC02-05CH11231, FWP no. CHPHYS02.

Appendix A Analytical Fluid Film Calculations

In this appendix, we provide several theoretical calculations for two-dimensional curved and deforming fluid films. First, we determine the material stresses using an irreversible thermodynamic framework [46]. The Helmholtz free energy density for a two-dimensional fluid film is calculated by using a Lagrange multiplier $\lambda$ to enforce areal incompressibility. The in-plane and out-of-plane equations of motion of such a film are then determined, and simplified in three cases: a completely flat film, a cylindrical film, and a cylindrical film with a bulge in the center. Finally, a linear stability analysis is performed on an initially cylindrical fluid film that is allowed to deform, which is found to be unstable when its length exceeds its circumference. We end by calculating the time scale associated with the instability.

A.1 Irreversible Thermodynamics and Material Stresses

We begin by briefly overviewing how the stresses of a general two-dimensional material are obtained via the balance law formulation and irreversible thermodynamic framework of Ref. [46]. As discussed in Sec. 3.1, the stress vectors $\bm{T}^{\check{\alpha}}$ are the boundary tractions across curves of constant $\zeta^{\alpha}$. Without loss of generality, the stress vectors can be decomposed in the $\{\bm{a}_{\check{\alpha}},\bm{n}\}$ basis as

[TABLE]

where $N^{{\check{\alpha}}{\check{\beta}}}$ and $S^{\check{\alpha}}$ are the in-plane and out-of-plane (shear) traction components, respectively. An angular momentum balance reveals there are two conditions to be satisfied in relating $N^{{\check{\alpha}}{\check{\beta}}}$ and $S^{\check{\alpha}}$ to the couple-stress components $M^{{\check{\alpha}}{\check{\beta}}}$—the latter of which are related to the couples acting on the edge of a surface patch (see Ref. [46] for more details). The first condition requires

[TABLE]

to be symmetric, an analogous condition to the symmetry of the stress tensor for bulk systems. Physically, $\sigma^{{\check{\alpha}}{\check{\beta}}}$ contains the couple-free components of the in-plane stresses. The second condition requires

[TABLE]

which indicates gradients of the moments lead to shear forces, in analogy to well-known beam bending examples in which the shear force is proportional to the spatial derivative of the moment. For two-dimensional surfaces which cannot sustain moments, such as fluid films, the couple-stress components $M^{{\check{\alpha}}{\check{\beta}}}=0$, such that $N^{{\check{\alpha}}{\check{\beta}}}$ is symmetric and the shear forces $S^{\check{\alpha}}=0$.

To determine the stresses in a general two-dimensional material, an irreversible thermodynamic framework is employed. We introduce a Helmholtz free energy per unit mass, $\psi$, for the material, which as a thermodynamic state function captures only the reversible behavior of the material. For the two-dimensional materials of interest, we assume $\psi$ depends on only the covariant metric $a_{{\check{\alpha}}{\check{\beta}}}$, the curvature components $b_{{\check{\alpha}}{\check{\beta}}}$, and the temperature $T=T(\zeta^{\alpha},\,t)$ with the functional dependence

[TABLE]

With a general form of the Helmholtz free energy density satisfying Eq. (A.4), the irreversible thermodynamic framework developed in Ref. [46] shows that in the linear irreversible regime, the couple-free stress components $\sigma^{{\check{\alpha}}{\check{\beta}}}$ and the couple-stress components $M^{{\check{\alpha}}{\check{\beta}}}$ are given by

[TABLE]

In Eq. (A.5), the first terms on the right-hand side provide the elastic contribution to the stresses and couple-stresses from the free energy $\psi$, while $\pi^{{\check{\alpha}}{\check{\beta}}}$ and $\omega^{{\check{\alpha}}{\check{\beta}}}$ describe the stresses due to in-plane and out-of-plane dissipative phenomena, respectively. In all cases considered, we assume there is no dissipation from rates of change of curvature, such that $\omega^{{\check{\alpha}}{\check{\beta}}}=0$. Furthermore, for isotropic and incompressible materials with in-plane fluidity, the viscous stresses are given by Eq. (18). Note that studies based on the Euler–Lagrange equations for lipid membranes do not contain forces resulting from the viscous contribution $\pi^{{\check{\alpha}}{\check{\beta}}}$ to the stress components $\sigma^{{\check{\alpha}}{\check{\beta}}}$.

At this point, we have a general method for determining the governing equations for an incompressible, two-dimensional material with in-plane fluidity. By choosing the form of the Helmholtz free energy density $\psi$, we calculate the stress and couple-stress components according to Eq. (A.5), determine the components of the tractions (A.1) through Eqs. (A.2) and (A.3), and substitute the tractions into the linear momentum balance (16). In what follows, we determine the Helmholtz free energy density for a fluid film, and then carry out such a procedure.

A.1.1 Helmholtz Free Energy Density

In this section, we calculate the Helmholtz free energy density of an incompressible fluid film, which is a constraint enforcing areal incompressibility through the Lagrange multiplier $\lambda$. To this end, we introduce a reference configuration of the surface, write the total surface energy such that there is no area change between the current and reference configurations, and then determine the form of the Helmholtz free energy density $\psi$.

Consider an arbitrarily curved patch of surface $\mathscr{P}$ which is deforming over time, and whose position is given by $\check{\bm{x}}(\zeta^{\alpha},t)$. At some time $t_{0}$, the reference patch ${\mathscr{P}_{\!0}}$ is specified by its position $\check{\bm{x}}_{0}$, defined to be

[TABLE]

The reference patch ${\mathscr{P}_{\!0}}$ is static and its properties by definition do not change over time. The areal mass density of the reference patch is denoted $\rho_{0}$, and the Jacobian determinant $J$ between the reference and current patch is given by $J:=\sqrt{\rule{0.0pt}{5.03749pt}}\overline{\det\,a_{{\check{\alpha}}{\check{\beta}}}}$ . A differential areal element $\mathrm{d}A$ on the reference patch is related to the corresponding areal element $\mathrm{d}a$ on the current patch according to $\mathrm{d}a=J\,\mathrm{d}A$. The global form of the conservation of mass can be written as

[TABLE]

As the reference patch ${\mathscr{P}_{\!0}}$ is stationary, Eq. (A.7) shows

[TABLE]

Areal incompressibility is enforced through the equation $J=1$, which indicates areas do not change (A.8). We introduce the Helmholtz free energy of the patch, $W$, as

[TABLE]

where $\lambda$ is a Lagrange multiplier enforcing the incompressibility constraint $J=1$. As $\psi$ is the Helmholtz free energy per unit mass, the total Helmholtz free energy of the patch is also given by

[TABLE]

In Eq. (A.10), we used the relation $\mathrm{d}a=J\,\mathrm{d}A$ to map the integral to the reference patch, and then substituted Eq. (A.8). Due to the arbitrariness of ${\mathscr{P}_{\!0}}$, the integrands of Eq. (A.9) and the right-hand side of Eq. (A.10) are equal, and thus

[TABLE]

With the result of Eq. (A.11), as well as the identities $\partial\rho/\partial a_{{\check{\alpha}}{\check{\beta}}}=-\tfrac{1}{2}\,\rho\,a^{{\check{\alpha}}{\check{\beta}}}$ and $\partial\rho/\partial b_{{\check{\alpha}}{\check{\beta}}}=0$, we find $\sigma^{{\check{\alpha}}{\check{\beta}}}$ to be given by Eq. (17) and $M^{{\check{\alpha}}{\check{\beta}}}=0$, such that $\bm{T}^{\check{\alpha}}=\sigma^{{\check{\alpha}}{\check{\beta}}}\bm{a}_{\check{\beta}}$.

A.2 In-Plane and Out-of-Plane Equations of Motion

As described in Appendix A.1, given the in-plane stresses $\sigma^{{\check{\alpha}}{\check{\beta}}}$ (17), the stress vectors of a fluid film are given by $\bm{T}^{\check{\alpha}}=\sigma^{{\check{\alpha}}{\check{\beta}}}\bm{a}_{\check{\beta}}$. By substituting the stress vectors into Eq. (16), applying the product rule, and using the identity $\bm{a}_{{\check{\beta}};{\check{\alpha}}}=b_{{\check{\alpha}}{\check{\beta}}}\bm{n}$, we obtain the equations of motion as

[TABLE]

We proceed to show how Eqs. (19) and (20) can be found by taking the dot product of Eq. (A.12) with the unit normal $\bm{n}$ and the in-plane basis vectors $\bm{a}_{\check{\alpha}}$, respectively.

Using Eq. (17) for $\sigma^{{\check{\alpha}}{\check{\beta}}}$, realizing $H=\tfrac{1}{2}a^{{\check{\alpha}}{\check{\beta}}}b_{{\check{\alpha}}{\check{\beta}}}$, and taking the dot product of Eq. (A.12) with the unit normal $\bm{n}$ yields

[TABLE]

where the pressure drop $p$ across the surface is given by $p=\rho\bm{b}\cdot\bm{n}$. The component form of $\pi^{{\check{\alpha}}{\check{\beta}}}$ can be more conveniently written as [46]

[TABLE]

Substituting Eq. (A.14) into Eq. (A.13) and using the identity $b^{{\check{\alpha}}{\check{\beta}}}b_{{\check{\alpha}}{\check{\beta}}}=4H^{2}-2K$ leads to the out-of-plane shape equation

[TABLE]

The result of Eq. (A.15) is presented in the main text as Eq. (20). When there are no flows, Eq. (A.15) simplifies to the well-known Young–Laplace equation for fluid films, $p+2\lambda H=0$.

To find the in-plane equations of motion, Eq. (A.12) is contracted with the in-plane basis vector $\bm{a}_{\check{\alpha}}$ and Eq. (17) is substituted into the result to yield

[TABLE]

The last term in Eq. (A.16) is simplified by first taking the covariant derivative of the viscous stresses provided in Eq. (A.14) to find

[TABLE]

Using the identities $b^{{\check{\gamma}}{\check{\beta}}}_{;{\check{\gamma}}}=2\,H_{,{\check{\gamma}}}\,a^{{\check{\gamma}}{\check{\beta}}}$ and $v^{\check{\gamma}}_{;{\check{\mu}}{\check{\gamma}}}=v^{\check{\gamma}}_{;{\check{\gamma}}{\check{\mu}}}+K\,v^{\check{\gamma}}\,a_{{\check{\gamma}}{\check{\mu}}}$, along with Eq. (15), Eq. (A.17) can be written as

[TABLE]

Substituting Eq. (A.18) into Eq. (A.16), one obtains the in-plane equations

[TABLE]

which, upon exchange of dummy indices, is identical to Eq. (19). Equation (A.19) is the extension of the Navier–Stokes equations to curved and deforming fluid films. We note that Eqs. (A.15) and (A.19) were first presented in Refs. [16, 17].

A.2.1 Flat Plane Equations

A flat surface is parametrized by $\zeta^{1}:=x$ and $\zeta^{2}:=y$, such that the position $\bm{x}$ is given by

[TABLE]

Employing the results of Sec. 2.1 yields $\bm{a}_{\alpha}=\bm{e}_{\alpha}$, $\bm{n}=\bm{e}_{z}$, $a_{\alpha\beta}=\delta_{\alpha\beta}$, $b_{\alpha\beta}=0$, $H=0$, $K=0$, and $\Gamma^{\alpha}_{\lambda\mu}=0$, where $\delta_{\alpha\beta}$ is the standard Kronecker delta. The velocity $\bm{v}$ is written as

[TABLE]

There is no normal component of the velocity in the normal direction, $v=\bm{v}\cdot\bm{n}=0$, because the fluid is constrained to lie in the plane. We substitute the geometric quantities and velocity equations into the continuity equation (15), shape equation (A.15), and in-plane equations (A.19) to obtain

[TABLE]

The first equation in (A.22) is the familiar continuity equation, the second equation says there is no pressure drop across a fluid flowing on a flat plane, and the last equations are the two-dimensional Navier–Stokes equations in which the pressure has been replaced by the negative surface tension.

A.2.2 Fixed Cylinder Equations

A cylindrical surface is parametrized with $\zeta^{1}:=\theta$ and $\zeta^{2}:=z$, which are the standard polar angle and axial position of a cylindrical coordinate system. The surface position $\bm{x}$ on the cylinder is given by

[TABLE]

which yields $\bm{a}_{1}=R\bm{e}_{\theta}$, $\bm{a}_{2}=\bm{e}_{z}$, $\bm{n}=\bm{e}_{r}$, $a_{\alpha\beta}=\mathrm{diag}\,(R^{2},\,1)$, $a^{\alpha\beta}=\mathrm{diag}\,(R^{-2},\,1)$, $b_{\alpha\beta}=\mathrm{diag}\,(-R,\,0)$, $H=-1/(2R)$, $K=0$, and $\Gamma^{\alpha}_{\lambda\mu}=0$. The velocity $\bm{v}$ is written as

[TABLE]

Due to the different units of $\bm{a}_{1}$ and $\bm{a}_{2}$, $v^{\theta}$ has units of angular velocity and $v^{z}$ has units of velocity. As is standard for a fixed surface, the normal velocity $v=\bm{v}\cdot\bm{n}=0$. Substituting the geometric equations into the continuity equation (15), shape equation (A.15), and in-plane equations (A.19), and ignoring inertial terms, one obtains, respectively,

[TABLE]

In the absence of flows, the shape equation (A.26) again reduces to the familiar Young–Laplace equation for cylinders: $p=\lambda/R$.

A.2.3 Equations for a Cylinder with a Bulge

In this section, we analytically calculate the velocity, surface tension, and normal pressure solution of the bulged cylinder shown in Fig. 4(a). We assume the solution is axisymmetric, such that all unknowns depend only on the axial position $z$, and also neglect inertial terms. Our analytical solution is used to validate the numerical solution, as described in the main text.

The bulged cylinder is parametrized with the axial distance along the cylinder $z$ and polar angle $\theta$. The cylinder radius, which is now a function of axial position only, is denoted $r(z)$ and given by

[TABLE]

where $z^{*}:=z/L$. Defining the quantities $\zeta^{1}:=\theta$ and $\zeta^{2}:=z$, the surface position is given by

[TABLE]

From the surface position we calculate $\bm{a}_{1}=r\bm{e}_{\theta}$, $\bm{a}_{2}=r^{\prime}\bm{e}_{r}+\bm{e}_{z}$, $\bm{n}=(\bm{e}_{r}-r^{\prime}\bm{e}_{z})/(1+(r^{\prime})^{2})$, $a_{\alpha\beta}=\mathrm{diag}\,(r^{2},\,1+(r^{\prime})^{2})$, $a^{\alpha\beta}=\mathrm{diag}\,(r^{-2},\,1/(1+(r^{\prime})^{2}))$, and $b_{\alpha\beta}=1/(1+(r^{\prime})^{2})\cdot\mathrm{diag}\,(-r,\,r^{\prime\prime})$, where $(\,\cdot\,)^{\prime}$ denotes differentiation with respect to $z$. The mean and Gaussian curvatures are

[TABLE]

respectively. The Christoffel symbols are also calculated, and those which are nonzero are given by

[TABLE]

With the known geometric quantities, we solve for the velocity and surface tension on the fixed surface. The velocity $\bm{v}$ is first decomposed in the $\{\bm{a}_{1},\,\bm{a}_{2}\}$ basis as

[TABLE]

where again the normal velocity $\bm{v}\cdot\bm{n}=0$ because the fluid is constrained to flow on the fixed surface. Substituting the velocity decomposition (A.33) and Christoffel symbols (A.32) into the continuity equation (15) and setting $v^{\theta}=0$, due to our axisymmetric assumption, yields

[TABLE]

Eq. (A.34) is separable as the cylinder shape $r(z)$ is known, and by specifying the inlet velocity $v^{z}(z=0)$ to be a given value $V$, we find the $z$-velocity solution to be

[TABLE]

Substituting the known $z$-velocity (A.35) and geometric relations into the in-plane $v^{z}$ equation (A.19) and neglecting inertia, one obtains

[TABLE]

Since $v^{z}$ is known from the continuity equation (A.34), Eq. (A.36) determines the surface tension $\lambda=\lambda(z)$. By setting $\lambda(z=0)$ to be a constant, $\lambda_{0}$, and integrating Eq. (A.36), we find the surface tension is given by

[TABLE]

Finally, by substituting the $z$-velocity, surface tension, and geometric quantities into the shape equation (A.15), the pressure $p=p(z)$ is found to be

[TABLE]

Eqs. (A.33), (A.37), (A.38), and our assumption $v^{\theta}=0$ constitute a solution to the problem, which is used to analyze the numerical solutions in Sec. 5.2.3 in the main text.

Choice of Cylinder Length

In our numerical solution of the bulged cylinder, $C^{1}$-continuous piecewise quadratic basis functions are used to represent the surface (36). In this basis, $r^{\prime\prime}$ is discontinuous, and as a result, the mean curvature $H$ and the pressure drop $p$ are discontinuous as well. Consequently, they are poorly represented in our finite-dimensional function space. To avoid such discontinuities in Sec. 5.2.3, a large cylinder length $L=4,\!000$ is chosen such that $r^{\prime}$ and $r^{\prime\prime}$ are approximately zero. In this limit, $H$ is given approximately by $H(z)\approx-1/(2\,r(z))$, which is continuous and, as it is a function of $z$, leads to a nontrivial solution.

A.3 Stability Analysis of a Deforming Cylinder

In this section, we analyze the linear stability of an initially cylindrical fluid film, which is able to deform over time. The film is acted on by a constant pressure drop $p$ across its surface, and is prescribed to have zero velocity everywhere on its boundary. A base state solution to the cylindrical fluid equations of Sec. A.2.2 is provided, and then the first-order perturbed equations are presented. The perturbed quantities are expanded in a Fourier basis, and the dispersion relation is solved for. It is found that fluid films are unstable when their length $L$ is greater than their circumference, $2\pi R$. The instability analysis also reveals the time scale of the growing, unstable modes. We end by showing how an axisymmetric assumption, and the manipulation of the perturbed equations, leads to the same theoretical predictions.

A solution to the unperturbed cylindrical equations (A.25)–(A.28) is given by

[TABLE]

In Eq. (A.39), $V$ is a constant speed and $\lambda_{0}$ is related to the constant pressure drop $p$ according to $\lambda_{0}=R\,p$. A perturbation of the film surface position is introduced as

[TABLE]

where $\epsilon$ is a small parameter and $\tilde{r}(\theta,z,t)$ is an $\mathcal{O}(R)$ radial perturbation. Importantly, the perturbation introduced in Eq. (A.40) both changes the geometry of the surface, and also allows for a normal velocity $v=\bm{v}\cdot\bm{n}$ given by $v=\epsilon\,\tilde{r}_{,t}$. In this section, a ‘tilde’ $(\tilde{\,\cdot\,})$ is used over all perturbed variables, such that the $\theta$-velocity, $z$-velocity, and surface tension are written as

[TABLE]

Setting the perturbed quantities to zero in Eq. (A.41) recovers the base state solution (A.39).

By calculating the perturbed geometric quantities describing the surface, substituting them along with Eq. (A.41) into the governing equations, and keeping only terms of order $\epsilon$, we obtain the perturbed equations for an initially cylindrical fluid film. The details of this calculation are omitted, however, the results are presented as

[TABLE]

Eqs. (A.42)–(A.45) are the perturbed continuity, in-plane $\theta$, in-plane $z$, and shape equations, respectively.

The four equations (A.42)–(A.45) contain four unknowns: $\tilde{r}$, $\tilde{v}^{\theta}$, $\tilde{v}^{z}$, and $\tilde{\lambda}$. To perform a linear stability analysis, we follow the canonical treatment of Ref. [107], and expand all unknowns in a Fourier basis as

[TABLE]

where $m\in\mathbb{Z}$ such that all quantities are periodic in $\theta$ and $q$ are the allowed wavenumbers given by $q=n\pi/L$ for $n\in\mathbb{Z}$. The coefficients $\hat{r}$, $\hat{v}^{\theta}$, $\hat{v}^{z}$, and $\hat{\lambda}$ are the Fourier coefficients. Substituting Eq. (A.46) into the perturbed equations (A.42)–(A.45) and recognizing modes of different $m$ and $q$ are independent leads to the equations

[TABLE]

Eqs. (A.47)–(A.50) are linear in the unknowns and can be written as the matrix equation $\mathbf{A}\,\mathbf{y}=\bm{0}$, for $\mathbf{y}=(\hat{r},\,\hat{v}^{\theta},\,\hat{v}^{z},\,\hat{\lambda})^{\mathrm{T}}$ and coefficient matrix $\mathbf{A}$. For a nontrivial solution to exist, the matrix $\mathbf{A}$ cannot be invertible, which is guaranteed by setting $(\det\mathbf{A})$ equal to zero. Doing so yields the dispersion relation

[TABLE]

Given the form of our ansatz is $\mathrm{e}^{i(m\theta+qz-\omega t)}$ for all terms (A.46), our solution is unstable when $\mathrm{Im}\,\{\omega\}>0$, which only occurs when $m^{2}+R^{2}q^{2}<1$. Accordingly, unstable modes are axisymmetric $(m=0)$ and only occur when $q<1/R$, or equivalently $L>n\pi R$.

Note that our linear stability analysis predicts the $n=1$ mode is unstable when $L>\pi R$. However, for such an instability to grow, material is required to be drawn in to increase the area of the fluid film, a phenomena which is incompatible with our zero-velocity boundary conditions. Accordingly, given our choice of boundary conditions, instabilities only occur when $L>2\pi R$. For the $n=2$ mode, the time scale $\tau$ is given by $1/\mathrm{Im}\{\omega\}$, found to be

[TABLE]

The theoretical time scale $\tau$ in Eq. (A.52) is presented in Eq. (57) of the main text.

With the understanding that axisymmetric, sinusoidal perturbations lead to unstable solutions, a simpler stability analysis can be performed. Assuming all perturbed quantities are independent of $\theta$, no base state $z$-velocity $(V=0)$, and no perturbed $\theta$-velocity $(\tilde{v}^{\theta}=0)$, Eqs. (A.42)–(A.45) simplify to

[TABLE]

We consider an initial perturbation of the form

[TABLE]

where as before $q=n\pi/L$, and seek to determine how this perturbation evolves in time.

To understand the time evolution of the perturbation, captured in $\tilde{r}_{,t}$, a series of algebraic manipulations are carried out. The initial perturbation (A.56) is substituted into the perturbed shape equation (A.55), leading to

[TABLE]

Taking the partial derivative of Eq. (A.53) with respect to $z$ and subtracting Eq. (A.54) yields

[TABLE]

To Eq. (A.58) we add the partial derivative of Eq. (A.57) with respect to $z$ and obtain

[TABLE]

Integrating Eq. (A.59) with respect to $z$, recognizing the integration constant is zero, and substituting the result into Eq. (A.57) leads to

[TABLE]

Eq. (A.60) describes the time evolution of $\tilde{r}$ only in terms of $\tilde{r}$ itself, and thus describes when the perturbation will grow or decay in time. When $q^{2}R^{2}<1$, or equivalently when $L>n\pi R$, the initial perturbation grows in time and the cylinder is unstable. As before, the boundary conditions exclude the $n=1$ mode, and cylindrical fluid films are expected to be unstable when $L>2\pi R$.

Eq. (A.60) reveals the time scale $\tau$ for the evolution of the initial perturbation. For the $n=2$ mode, the time scale is given by

[TABLE]

in agreement with the full perturbation analysis time scale (A.52).

Appendix B Numerical Solution Method

In this section, we discretize the fluid surface and define a finite number of basis functions on the discretized domain. The fundamental unknowns, as well as their arbitrary variations, are expressed in terms of these basis functions according to the Bubnov–Galerkin approximation. The resulting residual equations are temporally discretized using the backward Euler method, and then solved using Newton–Raphson iteration. We end by presenting the details of our numerical implementation of the Dohrmann–Bochev method [76], which removes numerical inf–sup instabilities arising from the incompressibility of the fluid film.

B.1 Discretization and Bubnov–Galerkin Approximation

For the space $\mathcal{V}_{h}$ (36), we introduce the nn (number of nodes) basis functions $\{N_{I}(\zeta^{\alpha})\}$ such that $\mathcal{V}_{h}=(\textrm{span}\,\{N_{I}(\zeta^{\alpha})\})^{3}$, and an arbitrary velocity $\bm{v}_{h}\in\mathcal{V}_{h}$ can be expressed as

[TABLE]

In Eq. (B.1) and from now on, the subscript ‘$h$’ is dropped for notational convenience. All variables refer to the approximate solution unless otherwise noted. We introduced the $3\times(3\cdot\texttt{nn})$ matrix of shape function values, $[\mathbf{N}]$, and the $(3\cdot\texttt{nn})\times 1$ vector of velocity degrees of freedom at time $t$, $[\mathbf{v}(t)]$, given respectively by

[TABLE]

where $[\bm{1}]$ is the $3\times 3$ identity matrix. In practice, the entire matrix $[\mathbf{N}]$ is never computed. Rather, as discussed in Sec. B.5, the shape functions over a single element are calculated and stored in a local shape function matrix.

As the mesh velocity $\bm{v}^{\mathrm{m}}\in\mathcal{V}_{h}$, we similarly expand the mesh velocity as

[TABLE]

where the mesh velocity degree of freedom vector $[\mathbf{v}^{\mathrm{m}}(t)]$ is defined analogously to $[\mathbf{v}(t)]$ in Eq. (B.2). The identity matrix $[\bm{1}]$, which appears in Eqs. (B.1) and (B.3) via the shape function matrix $[\mathbf{N}]$ (B.2)1, is chosen to be the Cartesian identity matrix. Accordingly, the velocity and mesh velocity degrees of freedom contained in $[\mathbf{v}]$ and $[\mathbf{v}^{\mathrm{m}}]$ are Cartesian components as well. Thus, when $\bm{v}$ and $\bm{v}^{\mathrm{m}}$ are calculated through Eqs. (B.1) and (B.3), we find their representation in the canonical Cartesian basis.

The space of surface tensions $\Lambda_{h}$ is spanned by the nln (number of Lagrange multiplier nodes) Lagrange multiplier basis functions $\{\bar{N}_{J}(\zeta^{\alpha})\}$, written as $\Lambda_{h}=\textrm{span}\,\{\bar{N}_{J}(\zeta^{\alpha})\}$. Here and from now on, discretized quantities with a ‘bar’ accent correspond to the Lagrange multiplier $\lambda$. We introduce $[\bar{\mathbf{N}}]$ as the $1\times\texttt{nln}$ vector of Lagrange multiplier basis functions and $[\bar{\bm{\lambda}}(t)]$ as the corresponding $\texttt{nln}\times 1$ vector of Lagrange multiplier degrees of freedom, given by

[TABLE]

such that the surface tension can be expanded as

[TABLE]

Note the basis functions contained in $[\mathbf{N}]$ and $[\bar{\mathbf{N}}]$ are determined by the discretization of the domain $\Omega$, so the only remaining unknowns are contained in $[\mathbf{v}(t)]$, $[\mathbf{v}^{\mathrm{m}}(t)]$, and $[\bar{\bm{\lambda}}(t)]$, which are together called the degrees of freedom of the system.

The arbitrary variations $\delta\bm{v}(\zeta^{\alpha})$, $\delta\bm{v}^{\mathrm{m}}(\zeta^{\alpha})$, and $\delta\lambda(\zeta^{\alpha})$ are discretized with the same basis functions as their unknown counterparts, according to the Bubnov–Galerkin approximation, and are expanded as

[TABLE]

The variation vectors $[\delta\mathbf{v}]$, $[\delta\mathbf{v}^{\mathrm{m}}]$, and $[\delta\bar{\bm{\lambda}}]$ are given by prefixing a ‘$\delta$’ to every entry of their respective counterparts, $[\mathbf{v}]$, $[\mathbf{v}^{\mathrm{m}}]$, and $[\bar{\bm{\lambda}}]$.

B.2 Residual Vector Equations

As shown in Appendix C.1, by substituting the discretized unknown variations (B.6) into Eq. (31) and introducing the shorthand $\tilde{\mathcal{G}}(t):=\mathcal{G}(\bm{u}(\zeta^{\alpha},\,t),\,\delta\bm{u}(\zeta^{\alpha}))$, where $\bm{u}$ represents all unknowns and $\delta\bm{u}$ is its arbitrary variation, the weak form can be written as

[TABLE]

for any time $t$. Since the weak form is linear in the unknown variations, the global velocity, mesh velocity, and surface tension residual vectors are respectively given by

[TABLE]

The calculation of the residual vectors, according to Eqs. (B.7) and (B.8), is provided in Appendix C.1. Since the discretized unknown variations are arbitrary, Eq. (B.7) is equivalent to simultaneously requiring the global residual vectors to be zero, written as

[TABLE]

B.3 Time Integration

The relations provided in Eq. (B.9) are true at any time $t$, however, they are only solved numerically at a set of $N$ discrete times $\{t_{1},\,t_{2},\ldots,t_{N}\}$. Thus, assuming a known solution $\bm{u}(\zeta^{\alpha},\,t_{n})$ satisfying Eq. (B.9) at time ${t_{n}}$, we seek the unknown solution $\bm{u}(\zeta^{\alpha},\,{t_{n+1}})$ satisfying Eq. (B.9) at time ${t_{n+1}}$. To this end, the fundamental unknowns are expressed as

[TABLE]

where $\Delta t:={t_{n+1}}-{t_{n}}$. According to Eq. (B.13), changes in the mesh velocity affect the surface position, which in turn affects the various geometric terms found in the residual vector equations (B.9).

Discretizing the fundamental unknowns (B.1, B.3, B.5), and introducing the notation $[\mathbf{v}_{n}]=[\mathbf{v}({t_{n}})]$, for example, Eqs. (B.10)–(B.13) are equivalently expressed as

[TABLE]

Note that while the position $\check{\bm{x}}$ is not a fundamental unknown, it is discretized as in Eq. (B.17) such that the change in mesh velocity, $[\Delta\mathbf{v}^{\mathrm{m}}_{n+1}]$, is properly accounted for when solving for $\bm{u}(\zeta^{\alpha},\,{t_{n+1}})$. Substituting Eqs. (B.14)–(B.17) into the residual equations (B.9) at time ${t_{n+1}}$ leads to a set of nonlinear algebraic equations in the unknowns $[\Delta\mathbf{v}_{n+1}]$, $[\Delta\mathbf{v}^{\mathrm{m}}_{n+1}]$, and $[\Delta\bar{\bm{\lambda}}_{n+1}]$, which are collectively gathered in the vector $[\Delta\mathbf{u}_{n+1}]$. These equations are solved via the Newton–Raphson method.

B.4 Newton–Raphson Iteration

We begin by initially guessing the solution at time ${t_{n+1}}$ is equal to the known solution at time ${t_{n}}$, written as $[\mathbf{u}_{n+1}]_{0}=[\mathbf{u}_{n}]$, for which $[\Delta\mathbf{u}_{n+1}]_{0}=[\bm{0}]$. We then generate a sequence of iterative solutions $[\mathbf{u}_{n+1}]_{i+1}$, where

[TABLE]

Substituting Eqs. (B.14)–(B.18) into Eq. (B.9), evaluated at time ${t_{n+1}}$, and keeping terms to first order in $[\Delta\mathbf{u}_{n+1}]_{i+1}$ yields

[TABLE]

In Eq. (B.19), the nine components of the tangent matrix $[\mathbf{K}_{n+1}]_{i}$ are the partial derivatives of the residual vector equations (B.8) with respect to the unknown vectors. The stiffness matrix components for an LE implementation are provided in Appendix C.2. For a different mesh velocity equation, only $[\mathbf{K}^{\mathrm{mv}}]$, $[\mathbf{K}^{\mathrm{mm}}]$, $[\mathbf{K}^{\mathrm{m}\bar{\lambda}}]$, and $[\mathbf{r}^{\mathrm{m}}]$ are modified, however, other changes may be required to avoid well-known ALE issues associated with arbitrary mesh motions—such as the violation of the geometric conservation law [78, 79].

Equation (B.19) may be compactly written as

[TABLE]

where the tangent matrix $[\mathbf{K}_{n+1}]_{i}$ and residual vector $[\mathbf{r}_{n+1}]_{i}$ are both calculated from $[\mathbf{u}_{n+1}]_{i}$, which is known. The degrees of freedom are updated according to Eq. (B.18), and the process is repeated until the 2-norm of $[\Delta\mathbf{u}_{n+1}]_{i+1}$ falls below a specified iteration threshold $\epsilon_{\texttt{iter}}$. At this point, the solution at time ${t_{n+1}}$, $\bm{u}(\zeta^{\alpha},\,{t_{n+1}})$, is assumed to be specified by $[\mathbf{u}_{n+1}]_{i+1}$ through Eqs. (B.1), (B.3), and (B.5).

B.5 Local Element Calculations

The residual vector and tangent matrix in Eq. (B.20) consists of integrals over the parametric domain $\Omega$. However, the basis functions $N_{I}$ and $\bar{N}_{I}$ have compact support and are only nonzero over a fixed number of elements. In practice, integrals are calculated locally over a single element ${\Omega^{e}}$ to generate the elemental (or local) residual vector $[\mathbf{r}^{e}]$ and tangent matrix $[\mathbf{K}^{e}]$, which are then assembled to form their global counterparts $[\mathbf{r}]$ and $[\mathbf{K}]$, respectively [81].

Over any element, the fixed number of nonzero basis functions is denoted nen (number of element nodes) and the fixed number of nonzero Lagrange multiplier basis functions is denoted neln (number of element Lagrange multiplier nodes). The corresponding elemental basis function matrices are given by

[TABLE]

Over a single element, at time ${t_{n}}$, the fundamental unknowns are expanded as

[TABLE]

where $[\mathbf{v}^{{e}}_{n}]$, $[\mathbf{v}^{\mathrm{m},{e}}_{n}]$, and $[\bar{\bm{\lambda}}^{e}_{n}]$ are the local degrees of freedom at time ${t_{n}}$. Eq. (B.23) is the local analog of Eqs. (B.1), (B.3), and (B.5). Furthermore, as the shape function matrices in Eqs. (B.21) and (B.22) correspond to all nonzero degrees of freedom over a single element, integrals over $\Omega$ are calculated as

[TABLE]

where the matrices $[\mathbf{N}]$ and $[\bar{\mathbf{N}}]$ in the integrand on the left-hand side are replaced by $[\mathbf{N}^{{e}}]$ and $[\bar{\mathbf{N}}^{{e}}]$, respectively, on the right-hand side. Standard finite element techniques are used to assumble $[\mathbf{r}]$ and $[\mathbf{K}]$ from their local counterparts, $[\mathbf{r}^{e}]$ and $[\mathbf{K}^{e}]$, according to Eq. (B.24) [81]. From now on, only local residual vector and tangent matrix calculations are provided. The details of the numerical integration are left to Appendix C.5.

B.6 Dohrmann–Bochev Implementation

In this section, we demonstrate how the Dohrmann–Bochev contribution to the weak form (42) is incorporated into our finite element framework. Since the global projection criterion (40) can be simplified to a condition over individual elements (41), the Dohrmann–Bochev method naturally takes advantage of local element calculations. In particular, Eq. (41) allows one to express $\breve{\lambda}$ in terms of $\lambda$ over any element. With this result, the Dohrmann–Bochev direct Galerkin expression (42) can be rewritten such that its contributions to the residual vector and tangent matrix are easily calculated.

We begin by considering Eq. (41) over a single element ${\Omega^{e}}$. Local basis functions $\breve{N}_{i}^{e}$ are chosen for the space $\mathbb{P}_{1}({\Omega^{e}})$ (39). Note that in this manuscript, local basis functions are indexed by lower-case Latin letters, while global basis functions are indexed by upper-case Latin letters. During numerical integration, the element ${\Omega^{e}}$ is mapped onto the unit square $\Omega_{\square}$, which is parametrized by $(\xi,\eta)\in[-1,1]\times[-1,1]$ (see Appendix C.5). On the unit square, the basis functions of $\mathbb{P}_{1}({\Omega^{e}})$, $\breve{N}^{e}_{i}$, are defined to span the space of all planes over the domain. As any plane can be expressed as $a+b\xi+c\eta=0$ for constants $a,\,b,\,c\in\mathbb{R}$, the basis functions of $\mathbb{P}_{1}({\Omega^{e}})$ are given by $\breve{N}^{e}_{1}=1$, $\breve{N}^{e}_{2}=\xi$, and $\breve{N}^{e}_{3}=\eta$. We define the row vector of basis functions, the column vector of projected unknowns at time ${t_{n}}$, and the column vector of projected unknown variations, respectively, as

[TABLE]

On a single element ${\Omega^{e}}$ at time ${t_{n+1}}$, the projected quantities $\breve{\lambda}$ and $\delta\breve{\lambda}$ can be written as

[TABLE]

With Eq. (B.26), and the separation of the projection criterion over individual elements (41), $\breve{\lambda}$ can be expressed in terms of $\lambda$ alone. Substituting Eqs. (B.5), (B.6)3, and (B.26) into Eq. (41) and rearranging terms yields

[TABLE]

Defining the matrices $[\mathbf{G}^{\breve{\lambda}\bar{\lambda},{e}}]$ and $[\mathbf{H}^{\breve{\lambda}\breve{\lambda},{e}}]$ according to

[TABLE]

and owing to the arbitrariness of $[\delta\breve{\bm{\lambda}}^{e}]$, Eq. (B.27) simplifies to

[TABLE]

The mass matrix $[\mathbf{H}^{\breve{\lambda}\breve{\lambda},{e}}]$ in Eq. (B.29) is invertible, such that the projected surface tension coefficient matrix $[\breve{\bm{\lambda}}^{e}_{n+1}]$ is given in terms of the original surface tension coefficient matrix $[\bar{\bm{\lambda}}^{e}_{n+1}]$ as

[TABLE]

Eqs. (B.26)1 and (B.30) provide a way to calculate $\breve{\lambda}$ given $\lambda$, on any element ${\Omega^{e}}$ and at any time ${t_{n+1}}$, in accordance with the projection definition (41).

Our final step is to substitute our results into Eq. (42) and calculate the Dohrmann–Bochev contribution to the tangent matrix and residual vector. To this end, the Dohrmann–Bochev contribution to the weak formulation (42) is rewritten in a more convenient form. The integrand of Eq. (42) is expressed as $(\delta\lambda-\delta\breve{\lambda})\,(\lambda-\breve{\lambda})=(\lambda\,\delta\lambda-\breve{\lambda}\,\delta\breve{\lambda})-\breve{\lambda}(\delta\lambda-\delta\breve{\lambda})-\delta\breve{\lambda}(\lambda-\breve{\lambda})$, such that with the projection definition (41), the Dohrmann–Bochev contribution to the weak form (42) is given by

[TABLE]

where \mbox{\mathcal{G}}_{\raisebox{-1.0pt}{\text{\tiny{DB}}}} is expressed as a sum over elements to take advantage of Eq. (B.30). Substituting Eqs. (B.5), (B.6)3, (B.26), and (B.30) into Eq. (B.31) yields

[TABLE]

where the second line is simplified by recognizing the integral in parenthesis is $[\mathbf{H}^{\breve{\lambda}\breve{\lambda},{e}}]$ (B.28). As the result of Eq. (B.32) is linear in $[\delta\bar{\bm{\lambda}}^{e}]$ and $[\bar{\bm{\lambda}}^{e}_{n+1}]$, its contributions to the elemental tangent matrix $[\mathbf{K}^{\bar{\lambda}\bar{\lambda},{e}}_{n+1}]$ and residual vector $[\mathbf{r}^{\bar{\lambda},{e}}_{n+1}]$ are easily calculated, as described in Appendix C. The local tangent matrices and residual vectors are then assembled to form their global counterparts using standard procedures.

B.6.1 Dohrmann–Bochev Example

The lid-driven cavity problem demonstrates the need for inf–sup stabilization. When the Dohrmann–Bochev terms are turned off and the Stokes flow problem is solved on a coarse $16\times 16$ mesh, checkerboard patterns in the $x$-velocity (Fig. 11(a)) and surface tension (Fig. 11(b)) are observed relative to the corresponding solution when the stabilization is included (Figs. 11(c) and 11(d), respectively). Moreover, a single square in the checkerboard pattern has length and width ${}^{1\!}/_{16}$, and as the mesh size is changed, the checkerboard pattern changes in tandem. Such oscillations indicate the violation of the inf–sup condition [87]. Fig. 11(d) is also useful in highlighting the sharp surface tension spikes at the top two corners of the domain, which prevent us from having a meaningful error analysis of the surface tension in this problem—as discussed in Sec. 5.1.1.

B.7 Lagrangian Finite Element Formulation

As described in Sec. 2.5, a Lagrangian scheme is recovered when the mesh velocity and material velocity are equal, i.e.

[TABLE]

The velocities $\bm{v}^{\mathrm{m}}$ and $\bm{v}$ belong to the same space of functions (see Sec. 4.2), as do their approximate counterparts $\bm{v}^{\mathrm{m}}_{h}$ and $\bm{v}_{h}$ (36). We thus attain a Lagrangian scheme by setting all the nodal mesh and material velocities to be equal at all times $t$, written as $[\mathbf{v}^{\mathrm{m}}(t)]-[\mathbf{v}(t)]=[\bm{0}]$. Assuming $[\mathbf{v}^{\mathrm{m}}(t_{0})]-[\mathbf{v}(t_{0})]=[\bm{0}]$ at the initial time $t_{0}$, the discretized Lagrangian mesh equation is given by

[TABLE]

for all times $t_{n+1}$ and iterations $i+1$ (see Secs. B.3 and B.4).

To implement a Lagrangian scheme, Eq. (B.34) replaces the linearized and discretized mesh residual equation (B.19)2 of the ALE formulation. Accordingly, the Lagrangian analog to Eqs. (B.19) and (B.20) is given by

[TABLE]

which is equivalently expressed as

[TABLE]

With Eq. (B.36), we implement a Lagrangian scheme within our ALE framework. The tangent matrix and residual vector components in Eq. (B.36) are calculated as before (see Appendix C), with the velocity $\bm{v}$ replacing the mesh velocity $\bm{v}^{\mathrm{m}}$. Note that in a Lagrangian implementation, nonzero velocity boundary conditions lead to a moving boundary, according to Eqs. (11) and (B.33).

Appendix C Tangent Matrix & Residual Vector Calculations

In this section, the residual vectors and tangent matrices required for our LE finite element implementation are calculated. To do so, several geometric quantities are linearized. We then detail how the residual vectors and tangent matrices are modified when the surface is constrained to not deform (relevant to Secs. 5.1 and 5.2) and when the body force is a pressure drop acting everywhere normal to the surface (used in Sec. 5.3). In the former case, inertial terms are included in the calculation as well. However, these terms are not included when the fluid film can deform. We end with a brief description of how integrals are numerically calculated, and then present an overview of our code, which shows the structure of the residual vector and tangent matrix calculations.

C.1 Residual Vector

In this section, the residual vectors $[\mathbf{r}^{\mathrm{v}}(t)]$, $[\mathbf{r}^{\mathrm{m}}(t)]$, and $[\mathbf{r}^{\bar{\lambda}}(t)]$ are calculated such that they satisfy Eqs. (B.7) and (B.8). We adopt a convention where if $[\delta\mathbf{v}^{e}]$ is a column vector, then $\partial\mathcal{G}/\partial[\delta\mathbf{v}^{e}]$ is a column vector as well. Also, the matrix $[\mathbf{N}]_{,{\check{\alpha}}}$ contains the partial derivatives of the basis functions with respect to $\zeta^{\alpha}$, and will be frequently used in our implementation.

We substitute the discretized arbitrary variations (B.6) and simplified form of the Dohrmann–Bochev weak form (B.32) into the direct Galerkin expression (43, 44), remove inertial terms as they are assumed to be negligible, and rearrange the remaining terms to obtain

[TABLE]

where the first negative sign in the third line arises because \mbox{\mathcal{G}}_{\raisebox{-1.0pt}{\text{\tiny{DB}}}} is subtracted in Eq. (44). In Eq. (C.1), we used the relation (for example) $\delta\bm{v}\cdot\rho\bm{b}=[\delta\mathbf{v}]^{{\mathrm{T}}}[\mathbf{N}]^{{\mathrm{T}}}\,\rho\bm{b}$, where $[\mathbf{N}]$ is defined in Eq. (B.2)1. The variation coefficient vectors in Eq. (C.1) are independent of $\zeta^{\alpha}$, and are moved outside the integrals to yield

[TABLE]

Defining the velocity, mesh velocity, and surface tension residual vectors as

[TABLE]

Eq. (C.2) can be written as in Eq. (B.7).

C.2 Tangent Matrix

The tangent matrix components resulting from the linearization given in Eq. (B.19) are defined as

[TABLE]

with all components calculated from the known degree of freedom vector at iteration $i$, $[\mathbf{u}_{n+1}]_{i}$, alone. Before calculating the tangent matrix components, several geometric quantities are linearized for their subsequent use in the tangent matrix calculation.

C.2.1 Unknown Linearization

To linearize the residual vectors $[\mathbf{r}^{\mathrm{v}}_{n+1}]_{i+1}$, $[\mathbf{r}^{\mathrm{m}}_{n+1}]_{i+1}$, and $[\bar{\bm{\lambda}}_{n+1}]_{i+1}$ about the state $\bm{u}({t_{n+1}})_{i}$, we rewrite Eq. (B.18) as

[TABLE]

where the time dependence is not written for notational simplicity. From now on, we also omit the subscript $(i+1)$ from the change in fundamental unknowns. In this section, we calculate how various geometric quantities, namely the position $\bm{x}$, basis vectors $\bm{a}_{\check{\alpha}}$ and $\bm{a}^{\check{\alpha}}$, metric components $a_{{\check{\alpha}}{\check{\beta}}}$ and $a^{{\check{\alpha}}{\check{\beta}}}$, normal $\bm{n}$, and Jacobian $J^{\mathrm{m}}$, differ between iterations $i$ and $i+1$ given the relations in Eq. (C.7).

With the backward Euler time discretization (B.13), the surface position is iterated according to

[TABLE]

from which we define $\Delta\check{\bm{x}}:=\Delta t\,\Delta\bm{v}^{\mathrm{m}}$. With an equation for the position, and the definitions $\bm{a}_{\check{\alpha}}=\check{\bm{x}}_{,{\check{\alpha}}}$ and $a_{{\check{\alpha}}{\check{\beta}}}=\bm{a}_{\check{\alpha}}\cdot\bm{a}_{\check{\beta}}$, it is straightforward to calculate

[TABLE]

The shorthand $\Delta\bm{a}_{\check{\alpha}}:=\Delta t\,\Delta\bm{v}^{\mathrm{m}}_{,{\check{\alpha}}}$ and $\Delta a_{{\check{\alpha}}{\check{\beta}}}:=\Delta t\,\Delta\bm{v}^{\mathrm{m}}_{,{\check{\alpha}}}\cdot\bm{a}_{\check{\beta}}^{i}+\Delta t\,\Delta\bm{v}^{\mathrm{m}}_{,{\check{\beta}}}\cdot\bm{a}_{\check{\alpha}}^{i}$ will often be used in our linearization calculations.

To determine $a^{{\check{\alpha}}{\check{\beta}}}_{i+1}$, the identity $\partial a^{{\check{\alpha}}{\check{\beta}}}/\partial a_{{\check{\mu}}{\check{\nu}}}=-\tfrac{1}{2}(a^{{\check{\alpha}}{\check{\mu}}}\,a^{{\check{\beta}}{\check{\nu}}}+a^{{\check{\alpha}}{\check{\nu}}}\,a^{{\check{\beta}}{\check{\mu}}})$ is required, such that $\Delta a^{{\check{\alpha}}{\check{\beta}}}$ is calculated according to $\Delta a^{{\check{\alpha}}{\check{\beta}}}=(\partial a^{{\check{\alpha}}{\check{\beta}}}/\partial a_{{\check{\mu}}{\check{\nu}}})\,\Delta a_{{\check{\mu}}{\check{\nu}}}$. After rearranging terms and using the symmetry of the metric components, one obtains

[TABLE]

The calculation of $\bm{a}^{\check{\alpha}}_{i+1}$ is straightforward, as $\bm{a}^{\check{\alpha}}=a^{{\check{\alpha}}{\check{\beta}}}\,\bm{a}_{\check{\beta}}$; however terms are rearranged to simplify our numerical implementation. Introducing the identity tensor $\bm{1}$ in curvilinear coordinates as $\bm{1}=\bm{a}_{\check{\lambda}}\otimes\bm{a}^{\check{\lambda}}+\bm{n}\otimes\bm{n}$, we find

[TABLE]

In the third line of Eq. (C.12) the metric components were used to raise indices, in the fourth line the identity tensor was introduced and the last two terms were rewritten using dyadic products, and in the last line dummy indices were rearranged. With the result of Eq. (C.12), one obtains

[TABLE]

To calculate $\Delta\bm{n}$, note $\bm{n}\cdot\bm{n}=1$, which implies $\Delta\bm{n}\cdot\bm{n}=0$, so $\Delta\bm{n}$ is completely specified by its in-plane components: $\Delta\bm{n}=\bm{a}^{\check{\alpha}}(\Delta\bm{n}\cdot\bm{a}_{\check{\alpha}})$. Next, observe $\Delta\bm{n}\cdot\bm{a}_{\check{\alpha}}=\Delta(\bm{n}\cdot\bm{a}_{\check{\alpha}})-\bm{n}\cdot\Delta\bm{a}_{\check{\alpha}}=-\bm{n}\cdot\Delta\bm{a}_{\check{\alpha}}$, as $\bm{n}$ is orthogonal to $\bm{a}_{\check{\alpha}}$. Writing $\Delta\bm{n}$ as $\Delta\bm{n}=-(\bm{a}^{\check{\alpha}}\otimes\bm{n})\Delta\bm{a}_{\check{\alpha}}$, we find $\bm{n}_{i+1}$ to be given by

[TABLE]

Finally, to calculate $\Delta J^{\mathrm{m}}$, we begin with the definition of the normal vector $\bm{n}$ as $\bm{n}=(\bm{a}_{\check{1}}\times\bm{a}_{\check{2}})/J^{\mathrm{m}}$, such that $\Delta\bm{n}=(\Delta\bm{a}_{\check{1}}\times\bm{a}_{\check{2}})/J^{\mathrm{m}}+(\bm{a}_{\check{1}}\times\Delta\bm{a}_{\check{2}})/J^{\mathrm{m}}+(\bm{a}_{\check{1}}\times\bm{a}_{\check{2}})\,\Delta(1/J^{\mathrm{m}})$. Substituting $(\bm{a}_{\check{1}}\times\bm{a}_{\check{2}})\,\Delta(1/J^{\mathrm{m}})=-\bm{n}\,(\Delta J^{\mathrm{m}}/J^{\mathrm{m}})$ into the expression for $\Delta\bm{n}$, contracting both sides with $\bm{n}$, and remembering $\Delta\bm{n}\cdot\bm{n}=0$ leads to

[TABLE]

Substituting $\bm{n}=(\bm{a}_{\check{1}}\times\bm{a}_{\check{2}})/J^{\mathrm{m}}$ into the above equation, and using the vector triple product rule, yields

[TABLE]

where in the second line we substituted $a_{{\check{1}}{\check{1}}}=(J^{\mathrm{m}})^{2}\,a^{{\check{2}}{\check{2}}}$, $a_{{\check{1}}{\check{2}}}=a_{{\check{2}}{\check{1}}}=-(J^{\mathrm{m}})^{2}\,a^{{\check{2}}{\check{1}}}$, and $a_{{\check{2}}{\check{2}}}=(J^{\mathrm{m}})^{2}\,a^{{\check{1}}{\check{1}}}$ because $a^{{\check{\alpha}}{\check{\beta}}}$ is the matrix inverse of $a_{{\check{\alpha}}{\check{\beta}}}$, and in the third line simplified with $\bm{a}^{\check{\alpha}}=a^{{\check{\alpha}}{\check{\beta}}}\,\bm{a}_{\check{\beta}}$. Accordingly, $J^{\mathrm{m}}_{i+1}$ is found to be

[TABLE]

C.2.2 Local Component Calculations

The tangent matrix components are calculated by substituting Eq. (C.7), along with the corresponding changes in geometric quantities described above, into the residual equations (C.3)–(C.5) and rearranging terms to match the structure shown in Eq. (B.19). As described in Sec. B.5, in practice only local residual vectors and tangent matrices are numerically evaluated. We accordingly provide local results here, where elemental vector and matrix quantities are denoted with a superscript $e$.

We present the calculation of $[\mathbf{K}^{\bar{\lambda}\mathrm{v},{e}}]$, $[\mathbf{K}^{\bar{\lambda}\mathrm{m},{e}}]$, and $[\mathbf{K}^{\bar{\lambda}\bar{\lambda},{e}}]$ to demonstrate our procedure, and the remaining components of the tangent matrix will be provided without calculation. Substituting Eq. (C.7) into Eq. (C.5) and keeping only terms up to first order in $\Delta\bm{u}$ over the element ${\Omega^{e}}$, one obtains

[TABLE]

Substituting $\Delta\bm{a}^{\check{\alpha}}$ from Eq. (C.13) and $\Delta J$ from Eq. (C.17) into Eq. (C.18), recognizing the first term on the right-hand side and the Dohrmann–Bochev term on the last line involving $[\bar{\bm{\lambda}}^{e}_{n}]$ combine to be the residual at iteration $i$, and rearranging terms, we find

[TABLE]

Equation (C.19) is now expressed entirely in terms of known quantities at iteration $i$ and the changes to the fundamental unknowns. The changes in the fundamental unknowns are discretized as in Eqs. (B.14)–(B.16), substituted into Eq. (C.19), and the unknown vectors are moved outside the resultant integrals to give

[TABLE]

In comparing Eq. (C.20) with Eq. (B.19)3, the tangent matrices for this equation are given by

[TABLE]

The presence of the matrices $[\mathbf{G}^{\breve{\lambda}\bar{\lambda},{e}}]$ and $[\mathbf{H}^{\breve{\lambda}\breve{\lambda},{e}}]$ in Eqs. (C.20) and (C.23) require a minor modification to the structure of standard finite element codes. For each element, $[\mathbf{G}^{\breve{\lambda}\bar{\lambda},{e}}]$ and $[\mathbf{H}^{\breve{\lambda}\breve{\lambda},{e}}]$ are calculated as integrals over the element (B.28); they are themselves assembled during the iteration loop over Gauss quadrature points. Once the loop is complete, the contributions of the last terms in Eq. (C.5), as well as the terms in Eq. (C.23), are calculated and assembled. The organization of this calculation is highlighted in Appendix C.6.

An identical procedure is followed to determine the remaining local tangent matrix components. For Eq. (B.19)1 such an analysis yields

[TABLE]

where body forces $\rho\bm{b}$ are assumed to be independent of the surface geometry. For Eq. (B.19)2 we find

[TABLE]

C.3 Modification for a Fixed Surface

In this section, we continue to work only with the local residual vector and tangent matrix. Constraining the fluid surface to be fixed introduces a new unknown variable $p$, and corresponding variation $\delta p$, which are locally expanded in the $\{N_{i}(\zeta^{\alpha})\}$ basis as

[TABLE]

where $[\mathbf{p}^{e}_{n+1}]$ are the local pressure degrees of freedom at time ${t_{n+1}}$. As $p$ is a scalar quantity, $[\mathbf{p}^{e}_{n+1}]$ and $[\delta\mathbf{p}^{e}]$ are $\texttt{nen}\times 1$ vectors and are multiplied by the $1\times\texttt{nen}$ vector of basis functions, which is denoted by $[\mathbf{N}^{{e}}_{(1)}]$ in contrast to the $3\times(3\cdot\texttt{nen})$ matrix $[\mathbf{N}^{{e}}]$ introduced in Eq. (B.2).

As the mesh is constrained to remain stationary, $\bm{v}^{\mathrm{m}}=\bm{0}$, all geometric quantities are fixed, and there are no mesh velocity degrees of freedom. The weak form provided in Eq. (50) is rewritten at time ${t_{n+1}}$ as in Eq. (B.7), yielding

[TABLE]

for all variation vectors $[\delta\mathbf{v}^{e}]$, $[\delta\bar{\bm{\lambda}}^{e}]$, and $[\delta\mathbf{p}^{e}]$. The residual vector $[\mathbf{r}^{\mathrm{p},{e}}_{n+1}]$ is defined by

[TABLE]

for the direct Galerkin expression provided in Eqs. (50) and (C.31). Following the same procedure as in Appendix C.1 leads to

[TABLE]

Furthermore, the first term in Eq. (C.3) is modified by replacing $\rho\bm{b}$ with $p\bm{n}$, which gives

[TABLE]

The form of $[\mathbf{r}^{\bar{\lambda},{e}}_{n+1}]$ is unchanged from Eq. (C.5).

In removing the mesh velocity degrees of freedom and introducing the pressure degrees of freedom, there are once more nine tangent matrix components. The matrices $[\mathbf{K}^{\mathrm{vv},{e}}_{n+1}]$, $[\mathbf{K}^{\mathrm{v}\bar{\lambda},{e}}_{n+1}]$, $[\mathbf{K}^{\bar{\lambda}\mathrm{v},{e}}_{n+1}]$, and $[\mathbf{K}^{\bar{\lambda}\bar{\lambda},{e}}_{n+1}]$ are unchanged, and the weak form (50) will modify only $[\mathbf{K}^{\mathrm{pv},{e}}_{n+1}]$ and $[\mathbf{K}^{\mathrm{vp},{e}}_{n+1}]$. Accordingly, $[\mathbf{K}^{\mathrm{p}\bar{\lambda},{e}}_{n+1}]$, $[\mathbf{K}^{\bar{\lambda}\mathrm{p},{e}}_{n+1}]$, and $[\mathbf{K}^{\mathrm{pp},{e}}_{n+1}]$ are all zero. As Eq. (C.33) is linear in $\bm{v}$, the tangent matrix component $[\mathbf{K}^{\mathrm{pv},{e}}_{n+1}]$ is given by

[TABLE]

With the residual vectors (C.33, C.34) and tangent matrix components (C.35, C.36), our general ALE finite element framework has been adapted to solve for fixed-surface fluid flows.

C.3.1 Inertial Contributions

For simplicity, the subscript $\rho$ is used to denote a quantity which arises only due to inertia. In this section, the inertial contributions to the local residual vector and tangent matrix are calculated. The inertial contribution to the weak form component $\mathcal{G}_{\mathrm{v}}$ (27)—denoted $\mathcal{G}_{\mathrm{v},\rho}$—can be written as

[TABLE]

Substituting Eq. (10), with $\bm{v}^{\mathrm{m}}=\bm{0}$ by construction, into Eq. (C.37) leads to the inertial contribution to the local residual vector at time ${t_{n+1}}$, which is expressed as

[TABLE]

The time derivative $\bm{v}^{\prime}$ is approximated with an implicit backwards-Euler discretization, given by

[TABLE]

Substituting Eq. (C.39) into the residual vector in Eq. (C.38), evaluated at time ${t_{n+1}}$, approximating $\bm{v}(\zeta^{\alpha},{t_{n+1}})$ with $\bm{v}_{i+1}=\bm{v}_{i}+\Delta\bm{v}$, and keeping terms to first order in $\Delta\bm{v}$, we obtain

[TABLE]

where $[\mathbf{r}^{\mathrm{v},{e}}_{\rho,n+1}]_{i}$ is given by Eq. (C.38) evaluated at $\bm{v}(t+\Delta t)_{i}$, written as

[TABLE]

From Eq. (C.40), the tangent matrix contribution from the inertial terms is found to be given by

[TABLE]

C.4 Modification for a Normal Body Force

In this section, we discuss how the finite element formulation for a curved and deforming fluid film is modified when the constant body force $\rho\bm{b}$ is replaced by $p\bm{n}$, a pressure normal to the surface. The body force term in Eq. (C.3) is modified, and the local residual vector at time ${t_{n+1}}$ is rewritten as

[TABLE]

The linearization of the first term yields an additional contribution to the tangent matrix $[\mathbf{K}^{\mathrm{vm},{e}}_{n+1}]$, given by

[TABLE]

Denoting the additional tangent matrix contribution as $[\mathbf{K}^{\mathrm{vm},{e}}_{\mathrm{p},n+1}]$, and calculate it by substituting $\Delta\bm{n}$ (C.14) into Eq. (C.44) and moving the mesh velocity degrees of freedom outside the integral, leads to

[TABLE]

By modifying our residual vector and tangent matrix as in Eqs. (C.44) and (C.45), we simulate body forces normal to the surface as it deforms over time.

C.5 Numerical Integration

To numerically evaluate integrals over a single element, we first consider a one-dimensional integral over the non-empty portion of the element boundary $\Gamma^{\text{t}}_{e}$ for which $\zeta^{2}$ is fixed and $\zeta^{1}$ ranges from $\zeta^{1}_{\,a}$ to $\zeta^{1}_{\,b}$. We map the parametric domain $\zeta^{1}\in[\zeta^{1}_{\,a},\,\zeta^{1}_{\,b}]$ to the interval $[-1,\,1]$, which is parametrized by the scaled variable $\xi$, and on this rescaled domain we sum the values of the integrand at $\xi=\{-\sqrt{{}^{3}\!/_{\!5}},\,0,\,\sqrt{{}^{3}\!/_{\!5}}\}$ with corresponding weights $\{^{5}\!/_{\!9},\,^{8}\!\!/_{\!9},\,^{5}\!\!/_{\!9}\}$, as is standard for numerical integration using Gaussian quadrature [108].

To calculate two-dimensional areal integrals over an element of the parametric domain ${\Omega^{e}}$, we perform an analogous procedure. The element ${\Omega^{e}}$ is parametrized by $(\zeta^{1},\,\zeta^{2})=[\zeta^{1}_{\,a},\,\zeta^{1}_{\,b}]\times[\zeta^{2}_{\,a},\,\zeta^{2}_{\,b}]$, which we map to the unit square $\Omega_{\square}$, parametrized by $(\xi,\eta)\in[-1,\,1]\times[-1,\,1]$. The integrand is then calculated at the nine quadrature points with their corresponding weights.

C.6 Code Structure

We have now presented the computational details relevant to our ALE finite element formulation. A high-level overview of how our code is structured is shown in Algorithm 1.

Appendix D Numerical Benchmarks

In this section, we supplement the validation of the LE finite element method in Sec. 5 by considering four simple benchmark problems for which analytical solutions are known. In each case, the numerical calculation of the $x$-velocity and surface tension are plotted against their analytical results. The $L^{2}$-errors of the velocity and surface tension upon mesh refinement are also provided. All simulations are run on a stationary, flat mesh corresponding to the spatial domain $(x,\,y)\in[0,1]\times[0,1]$, and in three of the four cases the analytical solution lies in our finite-dimensional solution space $\mathcal{U}_{h}$. The remaining case demonstrates the convergence of the surface tension as the mesh is refined. Dohrmann–Bochev stabilization is used in all cases.

D.1 Hydrostatic Fluid

We first consider a fluid at rest which is acted upon by gravity, such that the body force is given by $\rho\bm{b}=-\rho g\,\bm{e}_{y}$. Units are chosen such that $\rho g=1$, and the surface tension is specified to be zero at the bottom of the domain. The analytical solution is given by $\bm{v}=\bm{0}$ and $\lambda=y$. We solve the problem numerically by specifying $\bm{v}=\bm{0}$ on all four boundaries and $\lambda=0$ on the bottom edge. Results from our numerical solution are shown in Fig. 12; note the scale of the $x$-axis in part (a).

D.2 Couette Flow

We next simulate Couette flow, in which fluid fills the space between two parallel plates; the bottom plate is stationary while the top plate moves tangentially at constant speed $V$. The bottom plate is located at $y=0$, the top plate at $y=1$, and furthermore units are chosen such that $V=1$. In this case, the analytical solution is given by $\bm{v}=y\,\bm{e}_{x}$ and $\lambda=0$.

In our numerical solution, we specify Dirichlet boundary conditions rather than formulating the degrees of freedom to be periodic on the left and right edges of the domain. Thus $\bm{v}=\bm{e}_{x}$ on the top edge, $\bm{v}=\bm{0}$ on the bottom edge, and $\bm{v}=y\,\bm{e}_{x}$ on the left and right edges. Furthermore, $\lambda$ is specified to be zero on the left edge. The results of our simulation are provided in Fig. 13; note the $y$-axis in (b) is multiplied by $10^{16}$.

D.3 Couette Flow with a Body Force

To demonstrate the convergence of the surface tension in our numerical method, we again consider the case of Couette flow, yet now with a body force $\rho\bm{b}=-3y^{2}\,\bm{e}_{y}$. In this case, the analytical solution is given by $\bm{v}=y\,\bm{e}_{x}$ and $\lambda=y^{3}$. The numerical solutions and convergence of the surface tension are shown in Fig. 14.

D.4 Hagen–Poiseuille Flow

The final benchmark problem considered is Hagen–Poiseuille flow, for which a surface tension change across the length of a stationary channel drives flow. No-slip boundary conditions on the top and bottom walls of the channel, located at $y=1$ and $y=0$, respectively, are assumed. Units are chosen such that the surface tension change per length $\Delta\lambda/L=8$. Furthermore setting the viscosity to be unity and arbitrarily choosing $\lambda=0$ at $x=0$, the analytical solution is given by $\bm{v}=4\,y\,(1-y)\,\bm{e}_{x}$ and $\lambda=8x$.

We specify Dirichlet boundary conditions on all four edges of the square domain, rather than incorporating the surface tension change across the domain into the weak form through the boundary traction. We set $\bm{v}=0$ on the top and bottom edges, $\bm{v}=4\,y\,(1-y)\,\bm{e}_{x}$ on the left and right edges, and specify $\lambda=0$ on the left edge. Our simulation finds a parabolic velocity profile throughout the domain and correctly calculates the surface tension change across the domain. Numerical results are shown in Fig. 15. Once again, there is excellent agreement with the analytical solution.

Appendix E Movies

1. Curved and deforming fluid film

The movie provided at youtu.be/FUx8fGXuzqY shows an initially cylindrical fluid film which is perturbed and then allowed to dynamically evolve over time. The color and color bar in the movie show the surface tension, whose variations mediate an instability via in-plane fluid flow. The LE simulation is run on a $10\times 40$ mesh, from time $t=0$ to $t=35$. Snapshots of the simulation are found in Fig. 5 in the main text.

2. Lagrangian scheme fails lid-driven cavity benchmark

This movie, which can be found at youtu.be/FCoShaa_FhM with snapshots presented in Fig. 16, shows a Lagrangian finite element scheme modeling the well-known lid-driven cavity problem. Due to in-plane velocity gradients, elements deform over time. The color and color bar in the movie indicate the $x$-component of the velocity, which is incorrect at late times. The cavity is represented by a unit square, the mesh is $32\times 32$ elements, and the time step $\Delta t=1/640$. After 20 time steps, nodes on the top edge move by a mesh spacing due to the velocity boundary condition at the top edge. Note in particular the top right corner, where after 20 time steps two nodes coincide, causing the simulation to fail.

3. Comparison of Lagrangian and LE schemes in modeling a fluid film

This movie, provided at youtu.be/wnXuK6d3WTQ, compares an LE scheme (top) and Lagrangian scheme (bottom) in modeling an unstable cylindrical fluid film. The movie runs from time $t=0$ to $t=35$, the color and color bars indicate the surface tension, and both simulations use a $30\times 60$ mesh. While both schemes capture the same solution, the Lagrangian nodes move in-plane due to the fluid flow. As a result, at later times there are comparatively more nodes in the left spherical bulb, and fewer nodes in the central tubular region. The underresolution of the tube in the Lagrangian scheme leads to larger errors. Snapshots from this movie are shown in Figs. 10(d)–10(f) in the main text.

Appendix F List of Important Symbols

[TABLE]

[TABLE]

[TABLE]

[TABLE]