2.2.1. The Dohrmann–Bochev method: Numerical stabilisation

The Dohrmann–Bochev method prevents spurious surface tension oscillations by projecting the surface tension onto a space of discontinuous, piecewise linear functions. The function space, denoted $ \mkern 1.5mu\breve {\mkern -1.5mu L}$ , is defined in (3.7) following a discussion of the surface discretisation. The $ L^2$ projection of a given surface tension field $ \lambda \in L^2 (\varOmega )$ onto $ \mkern 1.5mu\breve {\mkern -1.5mu L}$ is denoted $ \breve {\lambda }$ and is defined through the relation (Dohrmann & Bochev Reference Dohrmann and Bochev2004)

(2.9) \begin{equation} \int _\varOmega \delta \breve {\lambda } \, \big ( \lambda - \breve {\lambda } \big ) \,{\textrm{d}} \varOmega = 0 \quad \text{for all}\quad \delta \breve {\lambda } \in \mkern 1.5mu\breve {\mkern -1.5mu L} . \end{equation}

In practice, the method of Dohrmann & Bochev (Reference Dohrmann and Bochev2004) is implemented by quadratically penalising surface tension deviations from the space $ \mkern 1.5mu\breve {\mkern -1.5mu L}$ , for which the quantity

(2.10) \begin{equation} \dfrac {\alpha ^{\textit{DB}}}{\zeta } \int _\varOmega \big ( \delta \lambda - \delta \breve {\lambda } \big ) \big ( \lambda - \breve {\lambda } \big ) \, {\textrm{d}} \varOmega \end{equation}

is subtracted from (2.8). Here, $ \zeta$ is the 2-D intramembrane viscosity: a material property with dimensions of energy × time/area. In addition, $ \alpha ^{\textit{DB}}$ is a user-chosen parameter with dimensions of area, such that (2.8) and (2.10) have the same dimensions – though the value of $ \alpha ^{\textit{DB}}$ is not observed to affect simulation results. The numerically stabilised weak formulation of the incompressibility constraint is written as

(2.11) \begin{equation} \mathcal{G}_\lambda = 0 \quad \text{for all}\quad \delta \lambda \in L^2 (\varOmega ) , \end{equation}

where $ \mathcal{G}_\lambda$ is the surface tension contribution to the direct Galerkin expression (Zienkiewicz & Taylor Reference Zienkiewicz and Taylor2014) given by

(2.12) \begin{equation} \mathcal{G}_\lambda \ := \mkern 1mu \int _\varOmega \delta \lambda \, \big ( \boldsymbol{a}^{\check {\alpha }} \boldsymbol{\cdot } \boldsymbol{v}_{, {\check {\alpha }}} \big ) \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega \ - \ \dfrac {\alpha ^{\textit{DB}}}{\zeta } \int _\varOmega \big ( \delta \lambda - \delta \breve {\lambda } \big ) \big ( \lambda - \breve {\lambda } \big ) \, {\textrm{d}} \varOmega . \end{equation}

In (2.12), $ \breve {\lambda }$ and $ \delta \breve {\lambda }$ are understood to be the $ L^2$ projections of their respective counterparts $ \lambda$ and $ \delta \lambda$ onto $ \mkern 1.5mu\breve {\mkern -1.5mu L}$ , according to (2.9).

2.3.1. The boundary conditions

One cannot determine a well-posed set of boundary conditions to (2.18) and (2.19) by inspection. A series of systematic developments for elastic shells (Green & Naghdi Reference Green and Naghdi1968; Steigmann Reference Steigmann1998, Reference Steigmann1999) underlie the formulation of the general lipid membrane boundary conditions (Capovilla, Guven & Santiago Reference Capovilla, Guven and Santiago2002; Arroyo & DeSimone Reference Arroyo and DeSimone2009; Rangamani et al. Reference Rangamani, Agrawal, Mandadapu, Oster and Steigmann2012; Sahu et al. Reference Sahu, Sauer and Mandadapu2017; Sauer & Duong Reference Sauer and Duong2017). In what follows, we highlight possible boundary conditions and provide their physical justification. Details of our own derivations, which reproduce earlier results, are provided in Chapter V, § 5(d) of Sahu (Reference Sahu2022).

The in-plane equations governing lipid flows (2.18) are identical to those governing a 2-D Newtonian fluid (Scriven Reference Scriven1960). We thus expect the boundary conditions to be similar to those of a fluid, in which one specifies either the velocity $ \boldsymbol{v}$ or (for a surface) the force per length $ \boldsymbol{F}$ on the boundary. For general 2-D materials, the force per length on the patch boundary is given by

(2.20) \begin{equation} \boldsymbol{F} = \boldsymbol{T}^{\check {\alpha }} \nu ^{}_{\check {\alpha }} - \left ( M^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu \nu ^{}_{\check {\alpha }} \mkern 1mu \tau ^{}_{\check {\beta }} \, \boldsymbol{n} \right )_{\! , {\check {\mu }}} \, \tau ^{\check {\mu }} \, , \end{equation}

which – for the case of lipid membranes – has bending, tensile and viscous contributions (recall $ \boldsymbol{\nu }$ and $ \boldsymbol{\tau }$ are boundary basis vectors; see figure 1). In our numerical implementation, on each edge of the membrane patch we specify a component of either the velocity $ \boldsymbol{v}$ or the force per length $ \boldsymbol{F}$ in each of the three Cartesian directions. We denote $ \varGamma _{\!\! F}^{\mkern 1mu j}$ and $ \varGamma _{\! v}^{\mkern 1mu j}$ as the respective portions of the boundary where $ F_{\! j} := \boldsymbol{F} \boldsymbol{\cdot } \boldsymbol{e}_{\!j}$ and $ v^{}_{\mkern -1mu j} := \boldsymbol{v} \boldsymbol{\cdot } \boldsymbol{e}_{\!j}$ are specified. For $ j \in \{ 1, 2, 3 \}$ , the intersection $ \varGamma _{\! v}^{\mkern 1mu j} \cap \varGamma _{\!\! F}^{\mkern 1mu j} = \varnothing$ and the closure of the union $ \overline \varGamma _{\! v}^{\mkern 1mu j} \cup \varGamma _{\!\! F}^{\mkern 1mu j} = \varGamma$ . The boundary conditions are expressed as

(2.21) \begin{equation} v^{}_{\mkern -1mu j} = \bar {v}^{}_{\mkern -1mu j} \quad \text{on } \varGamma _{\! v}^{\mkern 1mu j} \qquad \text{and} \qquad F_{\! j} = \mkern 4.5mu\overline {\mkern -4.5mu F \mkern -0.5mu}\mkern 0.5mu_{\! j} \quad \text{on } \varGamma _{\!\! F}^{\mkern 1mu j} , \end{equation}

where $ \bar {v}^{}_{\mkern -1mu j}$ and $ \mkern 4.5mu\overline {\mkern -4.5mu F \mkern -0.5mu}\mkern 0.5mu_{\! j}$ are known quantities that we prescribe.

The boundary conditions in (2.21) are necessary but not sufficient for a mathematically well-posed scenario. To see why, note that the membrane bending energy gives rise to the $ \varDelta _{{s}} H = a^{{\check {\mu }} {\check {\nu }}} \mkern 1mu H_{; {\check {\mu }} {\check {\nu }}}$ term in the shape equation (2.19). Since the mean curvature contains two spatial derivatives of the surface position through the curvature components $ b_{{\check {\alpha }} {\check {\beta }}} = \boldsymbol{n} \boldsymbol{\cdot } \boldsymbol{x}_{, {\check {\alpha }} {\check {\beta }}}$ , the $ \varDelta _{{s}} H$ term contains four derivatives of the surface position. Following canonical developments in the theory of beam bending (Timoshenko Reference Timoshenko1921, Reference Timoshenko1922), we expect to specify two conditions along the entire boundary: either the out-of-plane velocity or force per length, as well as either the slope of the surface or the boundary moment

(2.22) \begin{equation} M \, := \, M^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu \nu ^{}_{\check {\alpha }} \mkern 1mu \nu ^{}_{\check {\beta }} . \end{equation}

Since the velocity and force boundary conditions are already contained in (2.21), we need only specify one of the latter pair. More precisely, we partition the boundary into the disjoint sets $ \varGamma _{\! n}$ and $ \varGamma _{\! M}$ , with $ \varGamma _{\! n} \cap \varGamma _{\! M} = \varnothing$ and $ \overline { \mkern 1mu \varGamma _{\! n} \cup \varGamma _{\! M} } = \varGamma$ , and prescribe

(2.23) \begin{equation} \boldsymbol{n} \boldsymbol{\cdot } \boldsymbol{v}_{, {\check {\alpha }}} \, \nu ^{{\check {\alpha }}} = \bar {v}_\nu \quad \text{on } \varGamma _{\! n} \qquad \text{and} \qquad M = \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu \quad \text{on } \varGamma _{\! M} , \end{equation}

where $ \bar {v}_\nu$ and $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu$ are prescribed quantities. With (2.20)–(2.23), we have a mathematically well-posed set of boundary conditions for the governing equations.

2.3.2. The weak formulation

The weak formulation of the balance of linear momentum is obtained by first contracting (2.13) with an arbitrary velocity variation $ \delta \boldsymbol{v}$ . At this point, we recognise the four spatial derivatives contained in the shape equation (2.19) through the $ \varDelta _{{s}} H$ term will, in the subsequent weak form, yield two spatial derivatives of both the velocity variation and the surface position. Since $ \boldsymbol{v}$ and $ \boldsymbol{x}$ are assumed to lie in the same space of functions, both are elements of $ \boldsymbol{H}^2 (\varOmega )$ – for which higher-order basis functions, with continuous first derivatives, are required in the numerical implementation. In accordance with the boundary conditions in (2.21) and (2.23), the space of admissible material velocity variations $ \boldsymbol{\mathcal{V}}^{}_{\! 0}$ is expressed as

(2.24) \begin{equation} \boldsymbol{\mathcal{V}}^{}_{\! 0} \mkern 1mu := \mkern 1mu \Big \{\, \boldsymbol{u} \big(\zeta ^{\check {\alpha }} \big) \colon \varOmega \rightarrow \mathbb{R}^3 \text{ such that } \boldsymbol{u} \in \boldsymbol{H}^2 (\varOmega ), \, u^{}_{\!j} \big \rvert _{\varGamma _{\! v}^{\mkern 1mu j}} = 0, \, \big ( \boldsymbol{n} \boldsymbol{\cdot } \boldsymbol{u}_{, {\check {\alpha }}} \, \nu ^{\check {\alpha }} \big ) \big \rvert _{\varGamma _{\! n}} = 0 \,\Big \} \end{equation}

for $ j \in \{ 1, 2, 3 \}$ . We integrate the result of the contraction over the membrane area, and substitute the general form of the stress vectors (2.14) to obtain

(2.25) \begin{equation} \int _\varOmega \! \delta \boldsymbol{v} \boldsymbol{\cdot } \Big [\, \sigma ^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu \boldsymbol{a}^{}_{\check {\beta }} - \big ( M^{{\check {\beta }} {\check {\alpha }}} \mkern 1mu \boldsymbol{n} \big )_{; {\check {\beta }}} \,\Big ]_{; \alpha } J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega + \int _\varOmega \! \delta \boldsymbol{v} \boldsymbol{\cdot } \boldsymbol{f} \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega = 0 \quad \text{for all}\quad \delta \boldsymbol{v} \in \boldsymbol{\mathcal{V}}^{}_{\! 0} . \end{equation}

Starting with (2.25), a series of algebraic manipulations are required to determine the weak form of the linear momentum balance. The calculations rely on developments by Green & Naghdi (Reference Green and Naghdi1968) and Steigmann (Reference Steigmann1998), and can be found in Sauer et al. (Reference Sauer, Duong, Mandadapu and Steigmann2017), Sauer & Duong (Reference Sauer and Duong2017). Following the notation and development of our prior efforts (Sahu Reference Sahu2022, Chapter V, $\mkern 1mu$ § 5(d)), the weak formulation of the balance of linear momentum is expressed as

(2.26) \begin{equation} \mathcal{G}_{{v}} = 0 \quad \text{for all}\quad \delta \boldsymbol{v} \in \boldsymbol{\mathcal{V}}^{}_{\! 0} , \end{equation}

where

(2.27) \begin{align} \mathcal{G}_{{v}} &:= \mkern 1mu \int _\varOmega \dfrac {1}{2} \mkern 1mu \Big ( \delta \boldsymbol{v}_{, {\check {\alpha }}} \boldsymbol{\cdot } \boldsymbol{a}^{}_{\check {\beta }} \mkern 1mu + \mkern 1mu \delta \boldsymbol{v}_{, {\check {\beta }}} \boldsymbol{\cdot } \boldsymbol{a}^{}_{\check {\alpha }} \Big ) \, \sigma ^{{\check {\alpha }} {\check {\beta }}} J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega \, + \mkern 1mu \int _\varOmega \dfrac {1}{2} \mkern 1mu \Big ( \delta \boldsymbol{v}_{; {\check {\alpha }} {\check {\beta }}} \mkern 1mu + \mkern 1mu \delta \boldsymbol{v}_{; {\check {\beta }} {\check {\alpha }}} \Big ) \boldsymbol{\cdot } \boldsymbol{n} \, M^{{\check {\alpha }} {\check {\beta }}} J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega \nonumber\\[5pt] &\quad - \, \int _\varOmega \delta \boldsymbol{v} \boldsymbol{\cdot } \boldsymbol{f} \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega \ - \ \sum _{j = 1}^3 \mkern 1mu \int _{\varGamma _{\!\! F}^{\mkern 1mu j}} \delta v_j \, \mkern 4.5mu\overline {\mkern -4.5mu F \mkern -0.5mu}\mkern 0.5mu_{\! j} \, J^{}_{\varGamma } \, {\textrm{d}} \varGamma - \int _{\varGamma _{\! M}} \delta \boldsymbol{v}_{, {\check {\alpha }}} \, \nu ^{\check {\alpha }} \boldsymbol{\cdot } \boldsymbol{n} \, \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu \, J^{}_{\varGamma } \, {\textrm{d}} \varGamma . \end{align}

Note that jump forces can arise at corners of the patch and enter (2.27). They are assumed here to be zero. In (2.27), $ \delta v_j := \delta \boldsymbol{v} \boldsymbol{\cdot } \boldsymbol{e}_{\!j}$ is the $j$ th Cartesian component of the arbitrary velocity variation. It is useful to recognise that (2.27) is general to any material whose stress vectors can be expressed as in (2.14), with the corresponding boundary conditions in (2.21) and (2.23).

2.7.1. The weak formulation of the kinematic constraint

The mesh pressure is an unknown Lagrange multiplier field enforcing the ALE kinematic constraint in (2.6). The weak form of (2.6) is obtained by multiplying it with an arbitrary variation $ \delta p^{{{m}}} \in L^2 (\varOmega )$ , and integrating the result over the membrane surface to obtain

(2.37) \begin{equation} \int _\varOmega \delta p^{{{m}}} \, \boldsymbol{n} \boldsymbol{\cdot } \big ( \boldsymbol{v}^{{{m}}} - \, \boldsymbol{v} \big ) \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega = 0 \quad \text{for all}\quad \delta p^{{{m}}} \in L^2 (\varOmega ) . \end{equation}

Equation (2.37) bears some resemblance to (2.8), which was obtained from the incompressibility constraint. In both equations, the variation of a scalar Lagrange multiplier interacts with a vector velocity. Moreover, when (2.37) is discretised, our numerical simulations exhibit an instability reminiscent of the checkerboarding that arises when the LBB condition is violated. We thus hypothesise that both numerical instabilities are similar, and attempt to stabilise the mesh pressure with the Dohrmann–Bochev method. Following the results of § 2.2.1, we express the weak formulation of the kinematic constraint (2.6) as

(2.38) \begin{equation} \mathcal{G}_{{p}} = 0 \quad \text{for all}\quad \delta p^{{{m}}} \in L^2 (\varOmega ) , \end{equation}

where the numerically stabilised mesh pressure contribution to the direct Galerkin expression is given by (cf. (2.12), (2.37))

(2.39) \begin{equation} \mathcal{G}_{{p}} \ := \ - \int _\varOmega \delta p^{{{m}}} \, \boldsymbol{n} \boldsymbol{\cdot } \big ( \boldsymbol{v}^{{{m}}} - \, \boldsymbol{v} \big ) \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega \ - \ \dfrac {\alpha ^{DB} \ell ^2}{\zeta } \int _\varOmega \big ( \delta p^{{{m}}} - \delta \breve {p}^{\mkern 1mu {{m}}} \big ) \big ( p^{{{m}}} - \breve {p}^{\mkern 1mu {{m}}} \big ) \, {\textrm{d}} \varOmega . \end{equation}

In (2.39), $ \breve {p}^{\mkern 1mu {{m}}}$ and $ \delta \breve {p}^{\mkern 1mu {{m}}}$ are the $ L^2$ projections of $ p^{{{m}}}$ and $ \delta p^{{{m}}}$ onto the space $ \mkern 1.5mu\breve {\mkern -1.5mu L}$ , as defined through (2.9). The factor of $ \ell ^2$ , for a characteristic length $ \ell$ , is required for dimensional consistency. Our choice to employ the Dohrmann–Bochev method to stabilise $ p^{{{m}}}$ does not yet sit on a firm theoretical footing, and currently can only be justified with empirical success. We hope to investigate the mathematical nature of the stabilisation of $ p^{{{m}}}$ in a future paper.

2.7.2. The case where the mesh dynamics is purely viscous

For the case of purely viscous mesh dynamics, quantities are denoted with a superscript or subscript ‘ ${v}$ ’. The scheme is referred to as ‘ALE-viscous’, with the shorthand ‘Av’ or ‘ALE-v’. Since the mesh velocity field is area-compressible, the mesh motion resists both shearing and dilatation.

Strong formulation: we prescribe for the mesh analogue of the in-plane stresses and couple stresses to be (cf. (2.15)–(2.17))

(2.40) \begin{equation} \sigma ^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {v}} \mkern 1mu = \pi ^{{\check {\alpha }} {\check {\beta }}}_{{m}} \mkern 1mu := \zeta ^{{{m}}} \, \boldsymbol{v}^{{{m}}}_{\,\, , \mkern 1mu {\check {\mu }}} \boldsymbol{\cdot } \Big ( \boldsymbol{a}^{\check {\alpha }} \mkern 1mu a^{{\check {\mu }} {\check {\beta }}} + \boldsymbol{a}^{\check {\beta }} \mkern 1mu a^{{\check {\mu }} {\check {\alpha }}} \Big ) \quad \text{and} \quad \ M^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {v}} \mkern 1mu = 0 , \end{equation}

where $ \zeta ^{{{m}}}$ is a user-specified parameter that we refer to as the mesh viscosity. Importantly, since the mesh velocity is area-compressible, our choice in (2.40) resists both shear and dilatation of the mesh. By substituting (2.35) and (2.40) into (2.36), the dynamical equations governing the mesh are found to be (Sahu Reference Sahu2022, Chapter V, $\mkern 1mu$ § 5(c))

(2.41) \begin{align} \zeta ^{{{m}}} \Big (\! \varDelta _{{s}} v^{{{m}}}_{{\check {\alpha }}} + K \mkern 1mu v^{{{m}}}_{\check {\alpha }} + 2 \mkern 1mu (v^{{{m}}})_{, {\check {\alpha }}} \mkern 1mu H - 2 \mkern 1mu (v^{{{m}}})_{, {\check {\beta }}} \, b^{\check {\beta }}_{\check {\alpha }} - 2 \mkern 1mu v^{{{m}}} H_{, {\check {\alpha }}} + \big ( a^{{\mkern 3mu\check {\mkern -3mu\gamma }} {\check {\beta }}} \mkern 1mu (v^{{{m}}}_{\mkern 3mu\check {\mkern -3mu\gamma }})_{; {\check {\beta }}} - 2 \mkern 1mu v^{{{m}}} H \big )_{ ; {\check {\alpha }}} \Big ) = 0 \end{align}

and

(2.42) \begin{align} p^{{{m}}} + \zeta ^{{{m}}} \, \Big ( 2 \mkern 1mu b^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu (v^{{{m}}}_{{\check {\alpha }}})_{; {\check {\beta }}} - 8 \mkern 1mu v^{{{m}}} H^2 + 4 \mkern 1mu v^{{{m}}} K \Big ) = 0 , \end{align}

where $ \boldsymbol{v}^{{{m}}} = v^{{{m}}}_{{\check {\alpha }}} \mkern 1mu \boldsymbol{a}^{\check {\alpha }} + v^{{{m}}} \mkern 1mu \boldsymbol{n}$ . The solution to (2.41) is independent of $ \zeta ^{{{m}}}$ , and the mesh pressure in (2.42) will adjust so as to satisfy the kinematic constraint (2.6). The mesh dynamics accordingly do not depend on the choice of $ \zeta ^{{{m}}}$ . In all numerical results presented, we choose $ \zeta ^{{{m}}} = \zeta$ for simplicity.

Boundary conditions: since the mesh velocity satisfies differential equations (see (2.41) and (2.42)) rather than an algebraic equation (cf. (2.30), (2.32)), boundary conditions for the mesh velocity are required. To this end, we introduce the mesh analogue of the force per length on the boundary:

(2.43) \begin{equation} \boldsymbol{F}^{{m}}_{ {v}} = \boldsymbol{T}^{\check {\alpha }}_{{{m}}, \mkern 1mu {v}} \, \nu ^{}_{\check {\alpha }} = \sigma ^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {v}} \, \nu ^{}_{\check {\alpha }} \mkern 1mu \boldsymbol{a}^{}_{\check {\beta }} . \end{equation}

Recalling $ M^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {v}} = 0$ according to (2.40), the purely viscous mesh dynamics cannot sustain moments – as is the case for a fluid film (Sahu et al. Reference Sahu, Omar and Mandadapu2020b ). Consequently, neither slope nor moment boundary conditions are required on the boundary, since the mesh analogue of the moment $ M^{{m}} = 0$ identically throughout. With this understanding, the parametric boundary $ \varGamma$ is partitioned into the portions $ \varGamma _{\! v}^{{{m}} j}$ and $ \varGamma _{\!\! F}^{{{m}} j}$ , on which $ v^{{{m}}}_{\mkern -1mu j} := \boldsymbol{v}^{{{m}}} \boldsymbol{\cdot } \boldsymbol{e}_{\!j}$ and $ F^{{m}}_{\! j} := \boldsymbol{F}^{{m}} \boldsymbol{\cdot } \boldsymbol{e}_{\!j}$ are respectively specified. For simplicity in our numerical formulation, the mesh is assumed to be either static or force-free on its entire boundary, for which (cf. (2.21))

(2.44) \begin{equation} v^{{{m}}}_{\mkern -1mu j} = 0 \quad \text{on } \quad \varGamma _{\! v}^{{{m}} j} \quad \text{and} \quad F^{{m}}_{\! j} = 0 \quad \text{on} \quad \varGamma _{\!\! F}^{{{m}} j} . \end{equation}

Weak formulation: the weak formulation of the purely viscous ALE mesh behaviour is obtained by drawing analogy to the developments in § 2.3.2. First, the space of arbitrary mesh velocity variations is defined as

(2.45) \begin{equation} \boldsymbol{\mathcal{V}}^{{{m}}}_{\! 0, {v}} := \Big \{\, \boldsymbol{u} \big(\zeta ^{\check {\alpha }} \big) \colon \varOmega \rightarrow \mathbb{R}^3 \text{ such that } \boldsymbol{u} \in \boldsymbol{H}^2 (\varOmega ), u^{}_{\!j} \big \rvert _{\varGamma _{\! v}^{{{m}} j}} = 0 \,\Big \} . \end{equation}

By contracting (2.36) with an arbitrary mesh velocity variation $ \delta \boldsymbol{v}^{{{m}}} \in \boldsymbol{\mathcal{V}}^{{{m}}}_{\! 0, {v}}$ and repeating the steps of § 2.3.2, the weak formulation of the ALE-viscous mesh motion is found to be (cf. (2.26) and (2.27))

(2.46) \begin{equation} \mathcal{G}_{{m}}^{ {Av}} = 0 \quad \text{for all}\quad \delta \boldsymbol{v}^{{{m}}} \in \boldsymbol{\mathcal{V}}^{{{m}}}_{\! 0, {v}} , \end{equation}

where

(2.47) \begin{equation} \mathcal{G}_{{m}}^{{Av}} \ := \, \int _\varOmega \frac {1}{2} \, \Big( \delta \boldsymbol{v}^{{{m}}}_{, {\check {\alpha }}} \boldsymbol{\cdot } \boldsymbol{a}^{}_{\check {\beta }} + \delta \boldsymbol{v}^{{{m}}}_{, {\check {\beta }}} \boldsymbol{\cdot } \boldsymbol{a}^{}_{\check {\alpha }} \Big) \, \sigma ^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {v}} \, J^{}_{\! \varOmega }.\end{equation}

In (2.46) and (2.47), the superscript ‘Av’ refers to the choice of a purely viscous ALE mesh motion. The direct Galerkin expression for this motion is then written as

(2.48) \begin{align} & \mathcal{G}^{ {Av}} \, := \,\mkern 1mu \mathcal{G}_\lambda + \mathcal{G}_{{v}} + \mathcal{G}_{{m}}^{ {Av}} + \mathcal{G}_{{p}} = \mkern 1mu 0 \notag\\[5pt] &\quad \text{for all}\quad \delta \lambda \in L^2 (\varOmega ) , \, \delta \boldsymbol{v} \in \boldsymbol{\mathcal{V}}^{}_{\! 0} , \, \delta \boldsymbol{v}^{{{m}}} \in \boldsymbol{\mathcal{V}}^{{{m}}}_{\! 0, {v}} , \, \delta p^{{{m}}} \in L^2 (\varOmega ) . \end{align}

2.7.3. The case where the mesh dynamics is viscous and resist bending

We now consider the case where the mesh dynamics is viscous in response to in-plane shears and area dilatations, and additionally resist bending. The mesh motion is referred to as ‘ALE-viscous-bending’, with shorthand ‘Avb’ or ‘ALE-vb’. Corresponding quantities are denoted with a subscript or superscript ‘vb’.

Strong formulation: we choose for the bending resistance of the mesh to arise in the same manner as the membrane itself, with the viscous resistance identical to that of the ALE-viscous scenario. To this end, we choose (cf. (2.15)– (2.17) and (2.40))

(2.49) \begin{align} M^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {vb}} = k^{{m}}_{{b}} \mkern 1mu H \mkern 1mu a^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu + \, k^{{m}}_{{g}} \mkern 1mu \Big ( 2 \mkern 1mu H \mkern 1mu a^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu - \mkern 1mu b^{{\check {\alpha }} {\check {\beta }}} \Big ) \end{align}

and

(2.50) \begin{align} \sigma ^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {vb}} = k^{{m}}_{{b}} \mkern 1mu \Big ( H^2 \mkern 1mu a^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu - \mkern 1mu 2 \mkern 1mu H \, b^{{\check {\alpha }} {\check {\beta }}} \Big ) - k^{{m}}_{{g}} \mkern 1mu K \mkern 1mu a^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu + \, \pi ^{{\check {\alpha }} {\check {\beta }}}_{{m}} , \end{align}

where $ \pi ^{{\check {\alpha }} {\check {\beta }}}_{{m}}$ is defined in (2.40). In (2.49) and (2.50), $ k^{{m}}_{{b}}$ and $ k^{{m}}_{{g}}$ are user-specified parameters; we always choose $ k^{{m}}_{{b}} = k_{{b}}$ and $ k^{{m}}_{{g}} = k_{{g}}$ for simplicity. It is well known that bending terms from the Helfrich free energy do not enter the in-plane equations (Rangamani et al. Reference Rangamani, Agrawal, Mandadapu, Oster and Steigmann2012), and so the in-plane dynamical mesh equations in the ALE-vb case are given by (2.41). The out-of-plane dynamical equation is expressed as (cf. (2.19) and (2.42))

(2.51) \begin{equation} p^{{{m}}} + \zeta ^{{{m}}} \, \Big ( 2 \mkern 1mu b^{{\check {\alpha }} {\check {\beta }}} \mkern 1mu v_{{\check {\alpha }}; {\check {\beta }}} - 8 \mkern 1mu v \mkern 1mu H^2 + 4 \mkern 1mu v \mkern 1mu K \Big ) - k^{{m}}_{{b}} \Big ( 2 \mkern 1mu H^3 \mkern 1mu - \mkern 1mu 2 \mkern 1mu H \mkern 1mu K + \varDelta _{{s}} H \Big ) = 0 . \end{equation}

Contrary to the ALE-v scenario (2.42), here the relative magnitude of $ \zeta ^{{{m}}}$ and $ k^{{m}}_{{b}}$ do affect the resultant mesh pressure $ p^{{{m}}}$ . Investigating the relationship between $ p^{{{m}}}$ , $ \zeta ^{{{m}}}$ and $ k^{{m}}_{{b}}$ – as well as the ensuing mesh dynamics – is left to a later study.

Boundary conditions: we once again introduce the mesh analogue of the force per length on the boundary, which in this case is written as (cf. (2.20) and (2.43))

(2.52) \begin{equation} \boldsymbol{F}^{{m}}_{vb} = \boldsymbol{T}^{\check {\alpha }}_{{{m}}, \mkern 1mu {vb}} \, \nu ^{}_{\check {\alpha }} - \Big ( M^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {vb}} \, \nu ^{}_{\check {\alpha }} \mkern 1mu \tau ^{}_{\check {\beta }} \, \boldsymbol{n} \Big )_{\! , {\check {\mu }}} \, \tau ^{\check {\mu }} . \end{equation}

In principle, one could specify boundary forces for the mesh dynamics. However, as in the ALE-viscous case, we choose simple boundary conditions where the mesh is either static or force-free on the entire boundary (2.44). In addition, since the mesh couple stresses are non-zero, either slope or moment boundary conditions are required for $ \boldsymbol{v}^{{{m}}}$ – where the analogue of the boundary moment is calculated as $ M^{{{m}}} = M^{{\check {\alpha }} {\check {\mu }}}_{{m}} \mkern 1mu \nu ^{}_{\check {\alpha }} \mkern 1mu \nu ^{}_{\check {\mu }}$ (cf. (2.22)). At present, we take the simpler approach and specify zero-slope conditions on the entire boundary (cf. (2.23)):

(2.53) \begin{equation} \boldsymbol{n} \boldsymbol{\cdot } \boldsymbol{v}^{{{m}}}_{, {\check {\alpha }}} \, \nu ^{{\check {\alpha }}} = 0 \quad \text{on } \varGamma . \end{equation}

The ability to specify the slope of the mesh at a boundary will be relevant in the tether pulling scenario of § 4.2.

Weak formulation: at this point, it is straightforward to write the weak formulation of the mesh dynamics. The space of arbitrary mesh variations is given by

(2.54) \begin{equation} \boldsymbol{\mathcal{V}}^{{{m}}}_{\! 0, {vb}} := \Big \{\, \boldsymbol{u} \big(\zeta ^{\check {\alpha }} \big) \colon \varOmega \rightarrow \mathbb{R}^3 \text{ such that } \boldsymbol{u} \in \boldsymbol{H}^2 (\varOmega ), u^{}_{\!j} \big \rvert _{\varGamma _{\! v}^{{{m}} j}} = 0, \big ( \boldsymbol{n} \boldsymbol{\cdot } \boldsymbol{u}_{, {\check {\alpha }}} \, \nu ^{\check {\alpha }} \big ) \big \rvert _{\varGamma } = \, 0 \,\Big \} . \end{equation}

The weak formulation is then expressed as

(2.55) \begin{equation} \mathcal{G}_{{m}}^{\textit{A}vb} = 0 \quad \text{for all}\quad \delta \boldsymbol{v}^{{{m}}} \in \boldsymbol{\mathcal{V}}^{{{m}}}_{\! 0, {vb}} , \end{equation}

where we define (cf. (2.27))

(2.56) \begin{align} \mathcal{G}_{{m}}^{\textit{A}vb} &:= \ \int _\varOmega \dfrac {1}{2} \, \Big ( \delta \boldsymbol{v}^{{{m}}}_{, {\check {\alpha }}} \boldsymbol{\cdot } \boldsymbol{a}^{}_{\check {\beta }} + \delta \boldsymbol{v}^{{{m}}}_{, {\check {\beta }}} \boldsymbol{\cdot } \boldsymbol{a}^{}_{\check {\alpha }} \Big ) \, \sigma ^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {vb}} \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega \notag\\[5pt] & \quad + \, \int _\varOmega \dfrac {1}{2} \, \Big ( \delta \boldsymbol{v}^{{{m}}}_{; {\check {\alpha }} {\check {\beta }}} + \delta \boldsymbol{v}^{{{m}}}_{; {\check {\beta }} {\check {\alpha }}} \Big ) \boldsymbol{\cdot } \boldsymbol{n} \ M^{{\check {\alpha }} {\check {\beta }}}_{{{m}}, \mkern 1mu {vb}} \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega - \int _\varOmega \big ( \delta \boldsymbol{v}^{{{m}}} \boldsymbol{\cdot } p^{{{m}}} \mkern 1mu \boldsymbol{n} \big ) \, J^{}_{\! \varOmega } \, {\textrm{d}} \varOmega . \end{align}

The direct Galerkin expression for the scenario with the ALE-vb mesh motion is written as

(2.57) \begin{align} & \mathcal{G}^{\textit{A}vb} \, := \,\mkern 1mu \mathcal{G}_\lambda + \mathcal{G}_{{v}} + \mathcal{G}_{{m}}^{\textit{A}vb} + \mathcal{G}_{{p}} = \mkern 1mu 0 \notag\\[5pt] &\quad \text{for all}\quad \delta \lambda \in L^2 (\varOmega ) , \, \delta \boldsymbol{v} \in \boldsymbol{\mathcal{V}}^{}_{\! 0} , \, \delta \boldsymbol{v}^{{{m}}} \in \boldsymbol{\mathcal{V}}^{{{m}}}_{\! 0, {vb}} , \, \delta p^{{{m}}} \in L^2 (\varOmega ) . \end{align}

3.1.1. The membrane tension, mesh pressure and Dohrmann–Bochev projection

In general, the Lagrange multipliers $ \lambda$ and $ p^{{{m}}}$ are elements of $ L^2 (\varOmega )$ , and their finite-dimensional counterparts $ \lambda _h$ and $ p^{{{m}}}_h$ need not be restricted to lie in $ \mathcal{U}_h$ . However, for convenience in the numerical implementation, we do in fact choose

(3.6) \begin{equation} \lambda _h \in \, \mathcal{U}_h \quad \text{and} \quad p^{{{m}}}_h \in \, \mathcal{U}_h . \end{equation}

Since the same basis functions are used for the velocity and surface tension, the well-known inf–sup condition is violated and the LBB instability arises (Ladyzhenskaya Reference Ladyzhenskaya1969; Babuška Reference Babuška1973; Brezzi Reference Brezzi1974). We apply the method of Dohrmann & Bochev (Reference Dohrmann and Bochev2004) to stabilise unphysical oscillations in both the surface tension and mesh pressure. To do so, $ \lambda _h$ and $ p^{{{m}}}_h$ are projected onto the space

(3.7) \begin{equation} \mkern 1.5mu\breve {\mkern -1.5mu L} := \Big \{\, u \big(\zeta ^{\check {\alpha }} \big) : \varOmega \rightarrow \mathbb{R} \text{ such that } u \big \rvert _{\varOmega ^e} \in \mathbb{P}_1 (\varOmega ^e) \text{ for all }\varOmega ^e \in \mathcal{T}_h \,\Big \} ,\end{equation}

where $ \mathbb{P}_n (\varOmega ^e)$ denotes the space of polynomials of order $ n$ on $ \varOmega ^e$ . Since functions in $ \mkern 1.5mu\breve {\mkern -1.5mu L}$ form a plane over each rectangular element $ \varOmega ^e$ , they are discontinuous across finite elements.

4.1.1. The problem set-up

In our code, the membrane patch is initially a square of side length $ \ell$ in the $x$ – $y$ plane. The initial velocity and position are respectively given by

(4.2) \begin{equation} \boldsymbol{v} \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = \boldsymbol{0} \quad \text{and} \quad \boldsymbol{x} \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = x \mkern 1mu \boldsymbol{e}_x + y \mkern 1mu \boldsymbol{e}_y , \end{equation}

where $ \zeta ^{\check {\alpha }} \in \varOmega := [\mkern 1mu 0, 1] \times [\mkern 1mu 0, 1]$ and the real-space coordinates $x$ and $y$ are parametrised as

(4.3) \begin{equation} x \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = \ell \, \zeta ^{\check {1}} \quad \text{and} \quad y \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = \ell \, \zeta ^{\check {2}} . \end{equation}

The magnitude of the applied moment, denoted $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu (t)$ , is linearly ramped up over time until $ t = t_{{M}}$ – after which the boundary moment is held constant at the final value $\mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}$ (see figure 2 b):

(4.4) \begin{equation} \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu (t) = \begin{cases} 0, & \, t \le 0, \\ (t / t_{{M}}) \, \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}, & \, 0 \lt t \lt t_{{M}}, \\ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}, & \, t \ge t_{{M}} . \end{cases} \end{equation}

For times $ t \gt t_{{M}}$ , the applied boundary moment is constant in time and the membrane is bent into a portion of a cylinder, as shown in figure 2(c). The boundary forces and moments of this stationary solution are calculated in Appendix B and summarised in table 1. Here, $ r_{\mkern -1mu {{c}}}$ and $ \lambda _0$ are respectively the cylinder radius and surface tension, the latter of which is constant in the equilibrium configuration. We choose the Gaussian bending modulus $ k_{{g}} = - k_{{b}} / 2$ such that $ M = 0$ on the top and bottom edges, and no body force ( $ \boldsymbol{f} = \boldsymbol{0}$ ) such that $ \lambda _0 = k_{{b}} / (4 \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2})$ and $ \boldsymbol{F} = \boldsymbol{0}$ on the left and right edges (see table 1). In this case, the cylinder radius $ r_{\mkern -1mu {{c}}}$ is related to the final boundary moment $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}$ according to

(4.5) \begin{equation} r_{\mkern -1mu {{c}}} = \dfrac {k_{{b}}}{2 \mkern 1mu \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}} . \end{equation}

All of the boundary conditions prescribed in our numerical implementation for Lagrangian and Eulerian mesh motions are provided in table 2. For the ALE simulations, we choose the ALE-v mesh motion discussed in § 2.7.2. The first three rows of table 2 are then repeated for the mesh counterparts $ \bar {v}^{{{m}}}_{\mkern -1mu j}$ and $ \mkern 4.5mu\overline {\mkern -4.5mu F \mkern -0.5mu}\mkern 0.5mu^{{m}}_{\! j}$ . Recall that $ M^{{{m}}} = 0$ by construction, so there is no mesh analogue of the moment boundary condition. As a result, the membrane deforms only due to the physically applied moment $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu (t)$ – which is also the case for the Lagrangian and Eulerian simulations.

Table 1.

Moments and forces on the boundary of a portion of a cylindrical membrane with radius $ r_{\mkern -1mu {{c}}}$ , as shown in figure 2 and calculated in Appendix B. The constant surface tension $ \lambda _0$ is set by the net body force per area in the normal direction, $ f := \boldsymbol{f} \boldsymbol{\cdot } \boldsymbol{n}$ . In figure 2(a) the cylindrical basis vector $ \boldsymbol{e}^{}_\theta$ is tangent to the black curve and points to the left.

Table 2.

Boundary conditions prescribed in the numerical implementation. We choose to set $ k_{{g}} = - k_{{b}} / 2$ and $ \boldsymbol{f} = \boldsymbol{0}$ ; the latter yields $ \lambda _0 = k_{{b}} / (4 \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2})$ . The prescribed moment $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu(t)$ is given by (4.4). The first three rows are repeated for $ \bar {v}^{{{m}}}_{\mkern -1mu j}$ and $ \mkern 4.5mu\overline {\mkern -4.5mu F \mkern -0.5mu}\mkern 0.5mu^{{m}}_{\! j}$ in the case of a viscous ALE mesh motion.

4.1.2. Non-dimensionalisation

For the case of pure bending, there are five membrane parameters: the bending modulus $ k_{{b}}$ , final bending moment $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}$ , ramp-up time $ t_{{M}}$ , patch length $ \ell$ and intramembrane viscosity  $ \zeta$ . These parameters all dimensionally consist of mass, length and time. The pure bending scenario is thus completely described by two dimensionless quantities. The first is the Föppl–von Kármán number, which compares surface tension to bending forces (Sahu et al. Reference Sahu, Glisman, Tchoufag and Mandadapu2020a ). For a cylinder with $ \boldsymbol{f} = \boldsymbol{0}$ , the membrane surface tension is constant: $ \lambda = \lambda _0 := k_{{b}} / (4 \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2}) = \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}^2 / k_{{b}}$ , where (4.5) was substituted in the last equality (Evans & Yeung Reference Evans and Yeung1994). The Föppl–von Kármán number $ {\varGamma }$ is then given by

(4.6) \begin{equation} {\varGamma } = \dfrac {\lambda _0 \, \ell ^2}{k_{{b}}} = \bigg ( \dfrac {\mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}} \, \ell }{k_{{b}}} \bigg )^{\! 2} . \end{equation}

The second dimensionless quantity is the Scriven–Love number $ {{S\scriptstyle {L}}}$ , which captures dynamical effects and in this scenario is set by $ t_{{M}}$ (Sahu et al. Reference Sahu, Glisman, Tchoufag and Mandadapu2020a ). As the characteristic velocity scale of membrane deformations is given by $ \ell / t_{{M}}$ , the Scriven–Love number is expressed as

(4.7) \begin{equation} {{S\scriptstyle {L}}} = \dfrac {\ell ^2 \mkern 1mu \zeta }{k_{{b}} \, t_{{M}}} \, , \end{equation}

which can be understood as a ratio between the fundamental membrane time scale $ \ell ^2 \zeta / k_{{b}}$ and the ramp-up time $ t_{{M}}$ . For the results presented, we run simulations with $ \zeta = 1.0$ , $ k_{{b}} = 1.0$ and $ \ell = 1.0$ . The Föppl–von Kármán and Scriven–Love numbers are then set by choosing $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}}$ and $ t_{{M}}$ , respectively.

Figure 3.

Errors in the (a) surface tension and (b) mean curvature at the final time $ t_{{\kern-1.5pt f}}$ , according to (4.1), when the membrane patch is subjected to pure bending boundary conditions. The mesh consists of $ \textit {nel} = 1 / h^2$ parametric area elements, each of which is a square with side length $ h$ . Here, $ 1 / h$ ranges from $ 2^1$ to $ 2^7$ in powers of two. The labels ‘E’, ‘Av’ and ‘L’ refer to Eulerian, ALE-viscous and Lagrangian mesh motions. In all cases, the convergence of the error confirms our numerical implementation is working as expected. Relevant parameters are specified in § 4.1.3; we also choose $ \zeta = 1.0$ , $ k_{{b}} = 1.0$ and $ \ell = 1.0$ .

4.1.3. Results

We simulate pure bending scenarios where $ {\varGamma } = 0.25$ and $ {{S\scriptstyle {L}}} = 0.5$ , for which $ \mkern 4.5mu\overline {\mkern -4.5mu M \mkern -1.5mu}\mkern 1.5mu_{\! {{f}}} = k_{{b}} / (2 \mkern 1mu \ell )$ and $ t_{{M}} = 2 \mkern 1mu \ell ^2 \zeta / k_{{b}}$ (see (4.6) and (4.7)). Lagrangian, Eulerian and ALE-viscous simulations are carried out on meshes ranging from $ 2 \times 2$ to $ 128 \times 128$ . Other relevant parameters are $ \alpha ^{\textit{DB}} = \ell ^2$ , $ \alpha _{{{m}}}^{\textit {E}} = \zeta / \ell ^2$ , $ \Delta t = 0.1 \mkern 1mu \ell ^2 \zeta / k_{{b}}$ and $ t_{{\kern-1.5pt f}} = 4 \mkern 1mu t_{{M}}$ . The final configuration from one such simulation, at time $ t_{{\kern-1.5pt f}}$ , is shown in figure 2(c). The $L^2$ error in surface tension and mean curvature, relative to the analytical solution and calculated at time $ t_{{\kern-1.5pt f}}$ according to (4.1), is plotted as a function of mesh size in figure 3. Note that the analytical solution is valid only for a static patch; we thus deform the membrane slowly relative to the intrinsic time scale $ \ell ^2 \zeta / k_{{b}}$ and calculate errors after it has relaxed sufficiently. All three mesh motions converge towards the analytical solution upon mesh refinement, indicating our numerical implementation is working as expected. We believe the error from the Eulerian mesh motion is larger than that of the other two schemes because boundary conditions are not prescribed for the mesh velocity. Instead, material velocity boundary conditions are prescribed and weakly communicated to the mesh velocity through (2.33); edges are also the location where errors are the largest (see figure 2 c).

4.2.1. The problem set-up

In numerical simulations the membrane patch is initially a square of side length $ \ell$ in the $x$ – $y$ plane, centred at the origin. We set $ v_z = 0$ on the entire patch boundary, and also pin the centre node of each edge ( $ \boldsymbol{v} = \boldsymbol{0}$ ) to remove rigid body rotations and translations. In addition, zero-slope boundary conditions are enforced by constraining $ v_z = 0$ for nodes adjacent to the boundary – hereafter referred to as inner boundary nodes. Finally, an in-plane force/length $ \boldsymbol{\bar {F}}_\parallel = \lambda _0 \, \boldsymbol{\nu }$ is applied on the boundary, where $ \boldsymbol{\nu }$ is the in-plane unit normal at the patch boundary and $ \lambda _0$ is the static surface tension – which we choose. All boundary conditions are summarised in table 3, including those for the Eulerian and ALE schemes where $ \boldsymbol{v}^{{{m}}}$ is also a fundamental unknown.

Table 3.

Boundary conditions prescribed in our numerical implementation to pull a tether, with Lagrangian (LAG), Eulerian (EUL) and ALE-viscous-bending (ALE-vb) mesh motions. Here, ‘Bdry’ refers to nodes on the boundary, ‘Inner’ refers to inner nodes adjacent to the boundary, ‘Pull’ corresponds to nodes associated with the central element $ \varOmega ^e_{{p}}$ and ‘Pin’ refers to nodes at the centre of each edge. The latter conditions are required to prevent rigid body rotations and translations. The ‘Inner’ column enforces zero-slope boundary conditions.

The initial membrane state is given by

(4.8) \begin{equation} \boldsymbol{v} \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = \boldsymbol{0} \quad \text{and} \quad \lambda \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = \lambda _0 , \end{equation}

where $ \zeta ^{\check {\alpha }} \in \varOmega := [\mkern 1mu 0, 1] \times [\mkern 1mu 0, 1]$ . The corresponding membrane position is expressed as

(4.9) \begin{equation} \boldsymbol{x} \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = x \mkern 1mu \boldsymbol{e}_x + y \mkern 1mu \boldsymbol{e}_y \, , \end{equation}

with

(4.10) \begin{equation} x \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = \big ( \zeta ^{\check {1}} - 0.5 \big ) \, \ell \quad \text{and} \qquad y \big(\zeta ^{\check {\alpha }}, t \le 0 \big) = \big ( \zeta ^{\check {2}} - 0.5 \big ) \, \ell . \end{equation}

At time $ t = 0$ , we seek to vertically displace the centre of the membrane patch, where $ x = 0$ and $ y = 0$ . However, as shown in Appendix C, there is no unique way to specify the velocity at a single point given our use of B-spline basis functions, as they are not interpolatory. Instead, we vertically displace the portion of the membrane corresponding to the entire parametric element $ \varOmega ^e_{{p}}$ containing the point $ \zeta _{{p}} \hspace {-3.7pt} {}^{{\check {\alpha }}} := (0.5, 0.5)$ at the centre of the parametric domain:

(4.11) \begin{equation} \boldsymbol{v} \big(\zeta ^{\check {\alpha }}, t \ge 0 \big) = \boldsymbol{v}^{}_{{p}} (t) \quad \text{for all}\quad \zeta ^{\check {\alpha }} \in \varOmega ^e_{{p}} . \end{equation}

In (4.11), $ \boldsymbol{v}^{}_{{p}} (t)$ is a known function of time. It is given by $ v^{}_{\!{p}} \mkern 1mu \boldsymbol{e}_z$ , for constant $ v^{}_{\!{p}}$ , unless otherwise specified. Following the derivations in Appendix C, (4.11) is enforced by setting all nodal velocity degrees of freedom associated with $ \varOmega ^e_{{p}}$ to $ \boldsymbol{v}^{}_{{p}} (t)$ . The resultant pull force $ \boldsymbol{\mathcal{F}}_{\! {{p}}} (t)$ is calculated as the sum of the corresponding components of the residual vector.

4.2.2. Non-dimensionalisation

In the scenario under consideration, there are five membrane parameters: the bending modulus $ k_{{b}}$ , equilibrium surface tension $ \lambda _0$ , 2-D intramembrane viscosity $ \zeta$ , patch length $ \ell$ and speed of tube drawing $ v^{}_{\!{p}}$ . We hope to include the hydrodynamics of the surrounding fluid, including membrane permeability (Alkadri & Mandadapu Reference Alkadri and Mandadapu2024), and the effects of cytoskeletal contacts (Shi et al. Reference Shi, Graber, Baumgart, Stone and Cohen2018) in a future effort – and better compare our results with experiments (Cuvelier et al. Reference Cuvelier, Derényi, Bassereau and Nassoy2005; Brochard-Wyart et al. Reference Brochard-Wyart, Borghi, Cuvelier and Nassoy2006). Since the five quantities dimensionally involve only mass, length and time, two dimensionless numbers once again determine the evolution of the system. The Föppl–von Kármán number $ {\varGamma }$ is given by (Sahu et al. Reference Sahu, Glisman, Tchoufag and Mandadapu2020a )

(4.12) \begin{equation} {\varGamma } \, := \, \dfrac {\lambda _0 \, \ell ^2}{k_{{b}}} \, , \end{equation}

and can be interpreted in two ways given the morphological changes that occur when pulling a tether. First, when the membrane is nearly planar and height deflections are small, gradients in membrane shape occur over the length scale $ \ell$ ; $ {\varGamma }$ then quantifies the relative magnitude of tensile and bending forces. In our simulations, $ {\varGamma }^{-1}$ is small, and consequently a boundary layer of characteristic width $ \sqrt { k_{{b}} / \lambda _0 } \ll \ell$ develops at the point of application of the pull force (Powers et al. Reference Powers, Huber and Goldstein2002). As the membrane is pulled further, a tether grows from the boundary layer region. The tether is close in shape to a cylinder, and the radius of a cylindrical membrane is well known to be (Zhong-can & Helfrich Reference Zhong-can and Helfrich1989)

(4.13) \begin{equation} r_{\mkern -1mu {{c}}} \, := \, \sqrt { \dfrac {k_{{b}}}{4 \mkern 1mu \lambda _0} \, } . \end{equation}

Equation (4.13) can be obtained from energetic arguments alone, or equivalently, from a balance of bending and tensile forces. Note that in situations where there is a pressure jump across the membrane due to the surrounding fluid and $ f = \boldsymbol{f} \boldsymbol{\cdot } \boldsymbol{n} \ne 0$ , (4.13) and (4.17) are not valid. See the discussion in Chapter IX, § 1(a) of Sahu (Reference Sahu2022). With (4.13), the Föppl–von Kármán number (4.12) is understood as the ratio

(4.14) \begin{equation} {\varGamma } = \bigg ( \dfrac {\ell }{2 \mkern 1mu r_{\mkern -1mu {{c}}}} \bigg )^{\! 2} . \end{equation}

Equations (4.12)–(4.14) confirm that the Föppl–von Kármán number describes membrane energetics, as only lengths and the parameters $ k_{{b}}$ and $ \lambda _0$ – which respectively have dimensions of energy and energy per area – are involved.

The second dimensionless quantity is the Scriven–Love number $ {{S\scriptstyle {L}}}$ , which compares the magnitude of viscous and bending forces in shaping the membrane (Sahu et al. Reference Sahu, Glisman, Tchoufag and Mandadapu2020a ). In this scenario, we define the Scriven–Love number as

(4.15) \begin{equation} {{S\scriptstyle {L}}} \, := \, \dfrac {r_{\mkern -1mu {{c}}} \, \zeta \mkern 1mu v^{}_{\!{p}}}{k_{{b}}} . \end{equation}

We chose to include $ r_{\mkern -1mu {{c}}}$ – rather than $ \ell$ – in (4.15) because we are primarily concerned with the behaviour of the tether, rather than the entire patch. Upon substituting (4.13) into (4.15), rearranging terms and recognising $ \zeta / \lambda _0$ is the fundamental time scale of lipid flows (Sahu et al. Reference Sahu, Glisman, Tchoufag and Mandadapu2020a ), we find the Scriven–Love number can be equivalently expressed as

(4.16) \begin{equation} {{S\scriptstyle {L}}} \, := \, \dfrac {v^{}_{\!{p}}}{4 \mkern 1mu r_{\mkern -1mu {{c}}} \mkern 1mu ( \zeta / \lambda _0 )^{-1} } \, , \end{equation}

i.e. a ratio of the speed at which the tether is pulled to the natural velocity scale of lipid flows in the tube (Sahu et al. Reference Sahu, Glisman, Tchoufag and Mandadapu2020a ). Thus, as $ {{S\scriptstyle {L}}}$ tends to zero, we approach the quasi-static limit.

In our source code, when running tether pulling simulations, we choose $ \zeta = 1.0$ , $ k_{{b}} = 1.0$ and $ \lambda _0 = 0.25$ . Our choice yields a tube radius $ r_{\mkern -1mu {{c}}} = 1.0$ according to (4.13). The Föppl–von Kármán number $ {\varGamma }$ is altered by varying the patch size $ \ell$ , while the Scriven–Love number $ {{S\scriptstyle {L}}}$ is modified by changing the pull speed $ v^{}_{\!{p}}$ .

4.2.3. The comparison of different mesh motions

We begin by pulling a tether in a scenario with Föppl–von Kármán number $ {\varGamma } = 1024$ and Scriven–Love number $ {{S\scriptstyle {L}}} = 0.1$ . Given these dimensionless numbers, $ \ell = 64 \mkern 1mu r_{\mkern -1mu {{c}}}$ and $ v^{}_{\!{p}} = 0.1 \mkern 1mu k_{{b}} / (\zeta \mkern 1mu r_{\mkern -1mu {{c}}})$ according to (4.14) and (4.15). Snapshots from tether pulling simulations with each mesh motion are shown in figure 4, with corresponding videos in the MembraneAleFem.jl package repository (Sahu Reference Sahu2024). A zoomed-in view of the late-time snapshots, with and without the underlying mesh, are shown in figure 5 to emphasise the advantages of the ALE scheme. In addition, figure 6 shows how the $ z$ component of the pull force, $ \mathcal{F}_{\! {{p}}, \mkern 1mu z} := \boldsymbol{\mathcal{F}}_{\! {{p}}} \boldsymbol{\cdot } \boldsymbol{e}_z$ , varies as a function of the vertical displacement $ z^{}_{{p}}$ . We note that in the quasi-static limit ( $ {{S\scriptstyle {L}}} \rightarrow 0$ ), the steady-state pull force of a perfect cylinder is given by Evans & Yeung (Reference Evans and Yeung1994)

(4.17) \begin{equation} \mathcal{F}_{\! {eq}} \, := \, \dfrac {\pi \mkern 1mu k_{{b}}}{r_{\mkern -1mu {{c}}}} = 2 \mkern 1mu \pi \sqrt { k_{{b}} \mkern 1mu \lambda _0 } . \end{equation}

Note that some studies define $ \kappa := k_{{b}} / 2$ to be the bending modulus, in which case $ \mathcal{F}_{\! {eq}} = 2 \pi \sqrt { 2 \mkern 1mu \kappa \mkern 1mu \lambda _0 }$ . Since a tether pulled from a flat patch deviates slightly from a cylinder (Powers et al. Reference Powers, Huber and Goldstein2002), (4.17) serves as an excellent approximation for $ \mathcal{F}_{\! {{p}}, \mkern 1mu z}$ when $ {{S\scriptstyle {L}}} = 0$ .

In what follows, we comment on the efficacy of the three mesh motions: Lagrangian (LAG), Eulerian (EUL) and ALE-viscous-bending (ALE-vb). Note that additional variables are used to solve for membrane behaviour in the latter two schemes: the fundamental variables are $ \boldsymbol{v}$ and $ \lambda$ (LAG); $ \boldsymbol{v}$ , $ \boldsymbol{v}^{{{m}}}$ and $ \lambda$ (EUL); and $ \boldsymbol{v}$ , $ \boldsymbol{v}^{{{m}}}$ , $ \lambda$ and $ p^{{{m}}}$ (ALE-vb). We present the number of degrees of freedom for each motion, as well as the wall clock run time, in table 4. When appropriate, our results are compared to those of prior theoretical and numerical developments.

Table 4.

Number of degrees of freedom (DOFs) and wall clock run time for the three different mesh motions, corresponding to the results presented in figure 4. All computations were carried out on a single node of the Perlmutter system at the National Energy Research Scientific Computing Center, with area element calculations distributed across 32 cores.

$^{\textit{a}}$ Scaled to the total number of time steps, as the simulation failed.

Figure 4.

Snapshots from tether pulling simulations with a $ 65 \times 65$ mesh, in which three different mesh motions are employed. The colour bar indicates the surface tension, in units of $ k_{{b}} / r_{\mkern -1mu {{c}}}^2$ , for all snapshots. Times are measured in units of $ \zeta r_{\mkern -1mu {{c}}}^2 / k_{{b}}$ . (Left) Lagrangian simulations (LAG) successfully generate a tether. The morphological shape change from a tent to a tube occurs around time $ t = 120$ . Since the membrane is area incompressible, the patch boundary must be pulled inwards to accommodate the tether surface area. (Centre) The Eulerian mesh motion (EUL) is not able to form a tube. Rather, the tent morphology persists and bulges outward until the method fails around time $ t = 143$ . The material velocity degrees of freedom of the central element are constrained to move upwards; no such constraint is placed on the mesh velocity degrees of freedom. (Right) An ALE mesh motion that is viscous and resists bending (ALE-vb) successfully forms a tether. Both material and mesh velocities of the central element are constrained to move upwards. Since the mesh is area compressible, the patch boundary can be constrained to remain stationary as the tether develops. For all three mesh motions, $ \Delta t = 0.5$ , $ {{S\scriptstyle {L}}} = 0.1$ and $ {\varGamma } = 1024$ , for which $ \ell = 64 \, r_{\mkern -1mu {{c}}}$ and $ v^{}_{\!{p}} = 0.1 \, k_{{b}} / (\zeta r_{\mkern -1mu {{c}}} )$ according to (4.14) and (4.15). See also supplementary movies 1–3.

Figure 5.

Zoomed-in views of the $ t = 240 \, \zeta r_{\mkern -1mu {{c}}}^2 / k_{{b}}$ snapshots from figure 4, with the underlying mesh shown (top) and hidden (bottom). In the Lagrangian simulations (left), the mesh is drawn into the tube along with the lipids due to the areal incompressibility of the membrane. Mesh elements close to the diagonal of the square patch, where $ y = \pm \mkern 1mu x$ , become highly distorted – which leads to an artificially low surface tension in the transition region between the tether and surrounding membrane. Numerical artefacts are visible in the striation of the tension. The ALE-viscous-bending result (right), in contrast, shows a less distorted mesh because the choice of mesh constitution resists both shear and dilatation. No surface tension artefacts are visible, and a smooth tension gradient is observed from the flat patch into the tether.

Figure 6.

The $ z$ component of the pull force ( $ \mathcal{F}_{\! {{p}}, \mkern 1mu z}$ ) as a function of the $ z$ displacement. The dashed cyan line is the result from the linear theory, as presented in (4.20). The Eulerian simulation (blue circles, dotted line) is unable to form a tether, and is unphysical. The Lagrangian (red triangles, dashed line) and ALE (black squares, solid line) simulations capture the tent-to-tube transition, after which the Lagrangian steady-state pull force $ \mathcal{F}_{\! {\textit{ss}}}$ is slightly larger. Both overshoot the equilibrium pull force $ \mathcal{F}_{\! {eq}}$ (4.17) due to dynamical effects from tether pulling (see § 4.2.4). Numerical parameters are those specified in figure 4.

Lagrangian scheme: we first consider results from simulations with a Lagrangian mesh motion, as shown in the left column of figures 4 and 5 (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10553). A tent-like structure develops and grows until the vertical displacement $ z^{}_{{p}} := v^{}_{\!{p}} \mkern 1mu t \approx 8 \mkern 1mu r_{\mkern -1mu {{c}}}$ , for which $ t \approx 80 \, \zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2} / k_{{b}}$ . Around this point, the linear relationship between $ \mathcal{F}_{\! {{p}}, \mkern 1mu z}$ and $ z^{}_{{p}}$ breaks down, as shown in figure 6. The membrane undergoes a morphological transition between $ t \approx 80 \, \zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2} / k_{{b}}$ and $ t \approx 120 \, \zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2} / k_{{b}}$ , during which $ \mathcal{F}_{\! {{p}}, \mkern 1mu z}$ reaches its largest value. After the shape transition, a tether continues to be drawn at an approximately constant force, which we refer to as the steady-state force $ \mathcal{F}_{\! {\textit{ss}}}$ (see figure 6). However, a pronounced area of low surface tension develops in the region where the tether meets the flat patch (see figures 4 and 5). We comment on the unphysical nature of this development following a presentation of our ALE results.

The majority of our findings from Lagrangian tether pulling simulations, as well as the Lagrangian implementation itself, are not new. Powers et al. (Reference Powers, Huber and Goldstein2002) provided a detailed account of the axisymmetric, quasi-static membrane shape in response to vertical displacements, and along with Derényi et al. (Reference Derényi, Jülicher and Prost2002) numerically calculated the force versus displacement curve in the quasi-static limit. More recently, rate effects were included in axisymmetric simulations that focused primarily on the slip between two monolayer leaflets (Rahimi & Arroyo Reference Rahimi and Arroyo2012) – an effect not considered in the present work. General non-axisymmetric Lagrangian formulations were developed since then. The implementation in Rodrigues et al. (Reference Rodrigues, Ausas, Mut and Buscaglia2015) included rate effects when pulling a tether, but the force was not provided as a function of displacement and it is unclear if the membrane underwent the morphological transition from a tent to a tube in simulations. In contrast, Sauer et al. (Reference Sauer, Duong, Mandadapu and Steigmann2017) pulled a tether from the centre of an initially axisymmetric mesh and illustrated the tent-to-tube transition. While this implementation seems to be capable of capturing rate effects, the reported tether pulling results employed a numerical stabilisation scheme that did not capture the coupling between in-plane lipid flows and out-of-plane forces. Thus, to the best of our knowledge, we present the first non-axisymmetric Lagrangian simulation capturing the dynamics of a tether pulled from a membrane patch – though we do not consider the results to be novel, as Rodrigues et al. (Reference Rodrigues, Ausas, Mut and Buscaglia2015) and Sauer et al. (Reference Sauer, Duong, Mandadapu and Steigmann2017) may have been able to do the same.

Eulerian scheme: let us next examine results from an Eulerian mesh motion. As shown in supplementary movie 2 and the centre column of figure 4, the numerical implementation fails to form a tether. In addition, unphysical gradients develop in the surface tension, and the code fails at time $ t = 143 \, \zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2} / k_{{b}}$ . In what follows, we discuss why such a failure is to be expected, as a tether cannot form when the mesh motion is Eulerian.

We begin by introducing $ J$ and $ J^{{{m}}}$ as the relative area dilatations of the membrane and mesh, respectively. These dilatations are related to the corresponding flow fields according to (Sahu Reference Sahu2022, Chapter V, § 1(c))

(4.18) \begin{equation} \dfrac {\dot {J}}{J} = v^{\check {\alpha }}_{; {\check {\alpha }}} - 2 \mkern 1mu v H \quad \text{and} \quad \dfrac {\dot {J}^{{m}}}{J^{{{m}}}} = \big(v^{\check {\alpha }}_{{{m}}}\big)_{ ; \, {\check {\alpha }}} - 2 \mkern 1mu v^{{{m}}} \mkern -1mu H . \end{equation}

Since the membrane is incompressible, $ \dot {J} / J = 0$ and $ v^{\check {\alpha }}_{; {\check {\alpha }}} = 2 \mkern 1mu v H$ . To understand how the mesh dilates or compresses, we recognise that $ v^{{{m}}} = v$ and $ v^{\check {\alpha }}_{{{m}}} = 0$ according to (2.32). We thus find that

(4.19) \begin{equation} \dfrac {\dot {J}^{{m}}}{J^{{{m}}}} = - 2 \mkern 1mu v H \, \ne \, 0 . \end{equation}

At this point, we make three observations regarding the geometry and dynamics of the initial tent formation: (i) $ v \ge 0$ everywhere, (ii) $ H \lt 0$ in a central region where the surface is concave down, and (iii) $ H \ge 0$ elsewhere. Accordingly, $ \dot {J}^{{m}} / J^{{{m}}} \gt 0$ in the centre of the tent, and the mesh continuously dilates. In this manner, the mesh expands laterally and under-resolves the region where the morphological transition would occur.

Mesh dilatation at the patch centre, where $ \zeta ^{\check {1}} = 0.5$ and $ \zeta ^{\check {2}} = 0.5$ (see (4.10)), is quantified in figure 7. Here, the relative area change of the mesh $ J^{{{m}}}$ is plotted as a function of the membrane displacement $ z^{}_{{p}}$ . The dilatation by two orders of magnitude suggests an Eulerian simulation with a similar increase in the number of elements. Such a patch, even after the large dilatation, may be sufficiently resolved to undergo the morphological transition to a tube. However, even if a tether was formed, it likely could not be extended with an Eulerian mesh motion: lipid flows are primarily in-plane during tether elongation, so the Eulerian mesh would remain stationary and the tether would be poorly resolved. We thus believe the inability to pull a tether is a general failure of Eulerian methods, including those implemented previously (Reuther et al. Reference Reuther, Nitschke and Voigt2020; Sahu et al. Reference Sahu, Omar and Mandadapu2020b ). Moreover, since the ALE implementation of Torres-Sánchez et al. (Reference Torres-Sánchez, Millán and Arroyo2019) employed a mesh motion that was close to Eulerian, we are unsure if it would be able to successfully pull a tether.

Figure 7.

Mesh dilatation of the Eulerian simulation shown in figure 4. (a) The relative area change at the patch centre, $ J^{{{m}}}$ , is plotted as a function of the displacement $ z^{}_{{p}}$ . The dashed vertical line is the displacement at which a tent is expected to transition to a tether, according to Lagrangian and ALE data in figure 6. (Inset) At small displacements, simulation data agrees with an approximate solution to (4.19). Here $ v = \dot {z}_{{{p}}}$ ; if the mean curvature is linear in $ z^{}_{{p}}$ then $ H \sim - z^{}_{{p}} / r_{\mkern -1mu {{c}}}^{\, 2}$ and integrating (4.19) yields $ \ln J^{{{m}}} \sim (z^{}_{{p}} / r_{\mkern -1mu {{c}}})^2$ . The number 16 in the analytical expression is a fitting parameter. (b) The same data in a log–log plot suggests a power-law growth of the dilatation at large displacements, prior to the expected morphological transition (dotted vertical line) – though the data does not span even a single decade.

ALE scheme: we close by discussing the results of our ALE mesh motion. As shown in figures 4–6 and supplementary movie 3, a tether was successfully pulled with the ALE-vb scheme, and the surface tension in the region where the tether meets the flat patch was approximately constant. This latter point is quantified by considering the range of surface tension values at time $ t = 240 \, \zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2} / k_{{b}}$ . In the ALE simulation, the lowest (dimensionless) value of the surface tension is $ 0.249$ – approximately equal to the magnitude of the in-plane force per length $ \lvert \boldsymbol{\bar {F}}_\parallel \rvert = \lambda _0 = 0.25$ maintained on the boundary. Larger tension values, in this case up to $ 0.345$ , arise in the tether to draw in lipids during the dynamic process of pulling. In contrast, the minimum value of the tension in Lagrangian simulations is $ 0.191$ : well below $ \lvert \boldsymbol{\bar {F}}_\parallel \rvert$ and also less than the larger values (up to $ 0.446$ ) attained on the tether. Since lipids flow from regions of low to high tension, there does not seem to be a smooth flow of lipids from the flat patch into the tether. We thus find the ALE-viscous-bending mesh motion to be superior to its Lagrangian counterpart, and report only ALE-vb results for the remainder of our tether pulling analysis.

It is important to note that while the ALE-vb scheme successfully pulls a tether, the purely viscous ALE mesh motion – employed in § 4.1 – does not: as the centre of the patch is translated upwards, the tether tapers to a point. We believe an angular shape arises in the ALE-v simulation because of incompatible boundary conditions for the material and mesh. More specifically, since all nodes over a single element are translated upwards, a zero-slope condition is implicitly prescribed at the element boundary. Both the membrane and ALE-vb mesh can sustain such a requirement due to their inherent bending stiffness. In contrast, the ALE-v mesh has no bending rigidity and so cannot maintain zero slope on the element boundary. The tether resulting from the ALE-v scheme is consequently unphysical.

4.2.4. The geometry and dynamics of tether pulling

In comparing different mesh motions in § 4.2.3, all simulations were carried out at a single choice of $ {\varGamma }$ and $ {{S\scriptstyle {L}}}$ . We now investigate the effects of altering the patch size relative to the tether radius, as well as changing the speed of tether pulling – which respectively modify the Föppl–von Kármán and Scriven–Love numbers. We confirm that $ {\varGamma }$ dictates the initial slope of the force versus displacement curve, as previously observed (Derényi et al. Reference Derényi, Jülicher and Prost2002; Powers et al. Reference Powers, Huber and Goldstein2002; Sauer et al. Reference Sauer, Duong, Mandadapu and Steigmann2017), while $ {{S\scriptstyle {L}}}$ captures the overshoot of the tether pull force relative to the equilibrium (or quasi-static) result.

We first investigate the dependence of the pull force on the patch geometry. Figure 8(a) plots $ \mathcal{F}_{\! {{p}}, \mkern 1mu z}$ as a function of $ z^{}_{{p}}$ for different values of $ {\varGamma }$ , at fixed $ {{S\scriptstyle {L}}} = 0.1$ . We observe that $ {\varGamma }$ alters the initial slope of the force versus displacement curve, but does not affect the steady-state pull force after the tether is formed – the latter of which is expected to be independent of the patch size (see (4.17) and Powers et al. Reference Powers, Huber and Goldstein2002; Derényi et al. Reference Derényi, Jülicher and Prost2002). To calculate the initial dependence of $ \mathcal{F}_{\! {{p}}, \mkern 1mu z}$ on $ z^{}_{{p}}$ , we analyse the tent-like membrane shape when deformations are small. In this limit, the shape equation (2.19) is decoupled from in-plane lipid flows, and identical to its quasi-static counterpart (Sahu et al. Reference Sahu, Glisman, Tchoufag and Mandadapu2020a ). We thus take the small-deformation, quasi-static membrane tent result from Powers et al. (Reference Powers, Huber and Goldstein2002) and calculate the pull force to be given by (cf. (4.17))

(4.20) \begin{equation} \dfrac {\mathcal{F}_{\! {{p}}, \mkern 1mu z}}{k_{{b}} / r_{\mkern -1mu {{c}}}} = \dfrac {\pi }{\ln ({\varGamma } \! / 2)} \boldsymbol{\cdot }\dfrac {z^{}_{{p}}}{r_{\mkern -1mu {{c}}}} , \quad \text{or equivalently} \quad \quad \dfrac {\mathcal{F}_{\! {{p}}, \mkern 1mu z}}{\mathcal{F}_{\! {eq}}} = \dfrac {1}{\ln ({\varGamma } \! / 2)} \boldsymbol{\cdot }\dfrac {z^{}_{{p}}}{r_{\mkern -1mu {{c}}}} , \end{equation}

which is shown as the dashed cyan line in figure 6. Our calculation of (4.20) is provided in Appendix D, and the collapse of force versus displacement curves is shown in figure 8(b).

The dependence of the pull force on the speed of tether pulling is investigated next. Results from simulations with variable $ {{S\scriptstyle {L}}}$ and fixed $ {\varGamma } = 1024$ are plotted in figure 8(c). The observed increase in pull force with increasing Scriven–Love number can be justified as follows. During tether pulling, lipids in the surrounding region flow inwards towards the tether. A mass balance, with the assumption of axisymmetry, indicates that the flow speed is approximately $ v^{}_{\!{p}} \mkern 1mu r_{\mkern -1mu {{c}}} / r$ , where $ r$ is the distance from the $ z$ axis. Since the in-plane velocity grows as we approach the tether from its surroundings, a surface tension gradient is required to sustain the flow of lipids. As the tether is pulled more quickly, larger velocities and, thus, larger tension gradients ensue. A greater surface tension in the tether results, and leads to the larger pull force observed in simulations.

At present, we are unable to determine a general functional form for the long-time, steady-state pull force as a function of $ {{S\scriptstyle {L}}}$ due to the high degree of nonlinearity in the governing equations. The main difficulty is that the surface tension and membrane geometry enter both the in-plane and shape equations. Moreover, the viscous–curvature coupling forces – which arise due to the flow of lipids – lead to $ O({{S\scriptstyle {L}}})$ changes of the membrane shape in the tent-like region. Since we are unable to approximate how the steady-state pull force depends on the Scriven–Love number, we choose not to collapse the data contained in figure 8(c). Instead, we plot the steady-state pull force $ \mathcal{F}_{\! {\textit{ss}}}$ as a function of the Scriven–Love number in figure 8(d). We expect $ \mathcal{F}_{\! {\textit{ss}}} \approx \mathcal{F}_{\! {eq}}$ when $ {{S\scriptstyle {L}}} \rightarrow 0$ . In the limit where $ {{S\scriptstyle {L}}} \ll 1$ , we approximate dynamical effects by calculating the change in surface tension under the assumption of no membrane shape changes. We find that

(4.21) \begin{equation} \mathcal{F}_{\! {\textit{ss}}} ({{S\scriptstyle {L}}}) = ( 1 + 2 \mkern 1mu {{S\scriptstyle {L}}}) \, \mathcal{F}_{\! {eq}} \quad \text{for } {{S\scriptstyle {L}}} \, \ll \, 1 , \end{equation}

which is shown as the dotted grey line in figure 8(d). Evidently, (4.21) is a reasonable approximation when $ {{S\scriptstyle {L}}} \lt 0.1$ .

Figure 8.

The $ z$ component of the pull force ( $ \mathcal{F}_{\! {{p}}, \mkern 1mu z}$ ) as a function of $ z^{}_{{p}}$ for different values of $ {\varGamma }$ and $ {{S\scriptstyle {L}}}$ . (a) When $ {{S\scriptstyle {L}}} = 0.1$ and $ {\varGamma }$ is varied, the initial slope of the force vs displacement curve is altered according to the linear theory (see (4.20)). (b) By scaling the $ z$ displacement appropriately, the data collapses – with the steady-state pull force independent of $ {\varGamma }$ after tether formation. The cyan line has slope one. (c) When $ {\varGamma } = 1024$ and $ {{S\scriptstyle {L}}}$ is varied, the initial slope of the force vs displacement curve is unchanged. The tent-to-tube transition occurs at larger displacements, and the long-time pull force increases with $ {{S\scriptstyle {L}}}$ . (d) Long-time pull force, $ \mathcal{F}_{\! {\textit{ss}}}$ , as a function of $ {{S\scriptstyle {L}}}$ for $ {\varGamma } = 1024$ . The dotted grey line is the linear prediction from (4.21), which agrees with simulation results when $ {{S\scriptstyle {L}}} \ll 1$ – as highlighted by the zoomed-in inset. The nonlinear dependence of $ \mathcal{F}_{\! {\textit{ss}}}$ on $ {{S\scriptstyle {L}}}$ arises from the coupling between in-plane flows and out-of-plane shape deformations.

Figure 9.

Tether translation in the $ + \boldsymbol{e}_x$ (left column) and then $ - \boldsymbol{e}_x$ (right column) directions. Times are measured in units of $ \zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2} / k_{{b}}$ , the pull velocity is specified in (4.22), and all other parameters are identical to those in figure 4. The vertical red lines are a visual aid to highlight the lateral translation of the tether. For a video of the simulation results, see supplementary movie 4.

4.2.5. The lateral translation of a membrane tether

Once a tube is pulled from a patch of membrane, a lateral force applied at the end of the tether causes it to translate relative to the surrounding region. Lipids quickly flow and readjust to accommodate the translating tether due to the in-plane fluidity of the membrane. While tether translation via a lateral force is observed in in vitro studies (Datar et al. Reference Datar, Bornschlögl, Bassereau, Prost and Pullarkat2015; Gomis Perez et al. Reference Gomis Perez, Dudzinski, Rouches, Landajuela, Machta, Zenisek and Karatekin2022; Shi, Innes-Gold & Cohen Reference Shi, Innes-Gold and Cohen2022), the physics of tether translation remains poorly understood. A theoretical description of tether translation is difficult because the system is no longer axisymmetric, and the greatly simplified axisymmetric equations (Derényi et al. Reference Derényi, Jülicher and Prost2002; Powers et al. Reference Powers, Huber and Goldstein2002; Agrawal & Steigmann Reference Agrawal and Steigmann2009,Reference Agrawal and Steigmann2011; Omar et al. Reference Omar, Sahu, Sauer and Mandadapu2020) no longer apply. Tether translation also cannot be captured with general, non-axisymmetric Lagrangian numerical methods – as laterally translating the tether induces a rigid body translation of the entire patch (see Appendix E). Arbitrary Lagrangian–Eulerian finite element methods are thus required to simulate tether translation.

In figure 9 we present snapshots from a simulation of tether translation, which employed the ALE-vb mesh motion. Starting with the $ t = 240 \, \zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2} / k_{{b}}$ configuration shown in figures 4 and 5, we prescribe a lateral velocity in the $ x$ direction given by

(4.22) \begin{equation} \boldsymbol{v}^{}_{{p}} (t) = \begin{cases} \, + \mkern 1mu v^{}_{\!{p}} \mkern 1mu \boldsymbol{e}_x & \, \text{ for } \, 240 \lt t \boldsymbol{\cdot }k_{{b}} / \left(\zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2}\right) \lt 300, \\[2pt] \, - \mkern 1mu v^{}_{\!{p}} \mkern 1mu \boldsymbol{e}_x & \, \text{ for } \, 300 \lt t \boldsymbol{\cdot }k_{{b}} / \big(\zeta \mkern 1mu r_{\mkern -1mu {{c}}}^{\, 2}\big) \lt 420 , \end{cases} \end{equation}

where in both cases $ v^{}_{\!{p}} = 0.1 \mkern 1mu k_{{b}} / (\zeta \mkern 1mu r_{\mkern -1mu {{c}}})$ to be consistent with our choice $ {{S\scriptstyle {L}}} = 0.1$ . In this simulation, $ {\varGamma } = 1024$ is unchanged. We observe that the surface tension does not appreciably change, as the tether is translated slowly relative to the fundamental time scale of lipid rearrangements. Our results demonstrate that ALE methods can be used in scenarios where Lagrangian methods fail, and set the stage for future investigations of the forces, geometry and dynamics of tether translation. A video of the laterally pulled tether simulation is provided in the software repository (Sahu Reference Sahu2024).