A numerical approach for fluid deformable surfaces

Abstract

Fluid deformable surfaces show a solid–fluid duality which establishes a tight interplay between tangential flow and surface deformation. We derive the governing equations as a thin film limit and provide a general numerical approach for their solution. The simulation results demonstrate the rich dynamics resulting from this interplay, where, in the presence of curvature, any shape change is accompanied by a tangential flow and, vice versa, the surface deforms due to tangential flow. However, they also show that the only possible stable stationary state in the considered setting is a sphere with zero velocity.

1. Introduction

Fluid deformable surfaces are ubiquitous interfaces in biology, playing an essential role in processes from the subcellular to the tissue scale. Examples are lipid bilayers, the cellular cortex or epithelia monolayers. They all can be considered as fluidic thin sheets. From a mechanical point of view, they are soft materials exhibiting a solid–fluid duality: while they store elastic energy when stretched or bent, as solid shells, under in-plane shear, they flow as viscous two-dimensional fluids. This duality has several consequences: it establishes a tight interplay between tangential flow and surface deformation. In the presence of curvature, any shape change is accompanied by a tangential flow and, vice versa, the surface deforms due to tangential flow. The dynamics of this interplay strongly depends on the relation between fluid- and solid-like properties of the thin sheets. The growing interest in these phenomena in biology is in contrast with the available tools to numerically solve the governing equations. Even for surface fluids on stationary surfaces, where the governing equations have been known since the pioneering work of Scriven (1960), numerical tools have only been developed recently (see Nitschke, Voigt & Wensch (2012), Gross & Atzberger (2018) for simply-connected surfaces and Nitschke, Praetorius & Voigt (2017), Reuther & Voigt (2018b), Fries (2018), Lederer, Lehrenfeld & Schöberl (2019) for general surfaces). The governing equations for fluid deformable surfaces have been more recently derived in a different context (see Arroyo & DeSimone (2009), Salbreux & Jülicher (2017), Miura (2018)), but have never been solved in a general setting. Recent approaches (Mietke, Jülicher & Sbalzarini 2019; Torres-Sanchez, Millan & Arroyo 2019; Sahu et al. 2020) are restricted to the Stokes limit, simply-connected surfaces or axisymmetric settings. We overcome these limitations and provide a general numerical approach for fluid deformable surfaces.

The motivation to consider fluid deformable surfaces without the surrounding bulk phases results from the theoretical interest to explore them without any additional influence and the limit of a large Saffman–Delbrück number. This number describes the relation between the viscosities of the surface and the typically less viscous bulk fluid and allows the decoupling of surface and bulk flows (Saffman & Delbrück 1975).

The paper is structured as follows. In § 2, we sketch the derivation of the governing equations as a thin film limit and compare them with existing models for special cases. Section 3 describes the numerical approach, which is based on evolution of geometric quantities and a generic finite element formulation for tensor-valued surface partial differential equations. Numerical examples demonstrating the tight interplay between tangential flow and surface deformation as well as convergence tests for the numerical approach are provided in § 4. Conclusions are drawn in § 5.

2. Mathematical modelling

We start from a slightly more general Navier–Stokes equation in the thin film $\Omega _{\xi }(t) = \mathcal {S}(t) \times [-{\xi }/{2},\ {\xi }/{2}] \subset \mathbb {R}^3$ with a regular evolving surface $\mathcal {S}(t)$, film thickness $\xi$ and surface normal $\boldsymbol {\nu }$. In the Eulerian description, it reads

(2.1)\begin{gather} \partial_t\boldsymbol{V} + \boldsymbol{\nabla}_{\!\!\boldsymbol{V}}{\boldsymbol{V}} = \textrm{div}\boldsymbol{\Sigma}_{P} + \frac{1}{{Re}}\boldsymbol{\Delta} \boldsymbol{V}, \end{gather}

(2.2)\begin{gather}\textrm{div}\boldsymbol{V} = 0, \end{gather}

with $\boldsymbol {V}$ velocity, $\boldsymbol {\nabla }_{\!\!\boldsymbol {V}}{}$ directional derivative, ${Re}$ Reynolds number, $\boldsymbol {\Sigma }_{P} = -P\boldsymbol {\pi } - \phi \boldsymbol {\nu }\otimes \boldsymbol {\nu }$, where $P$ is the pressure, $\boldsymbol {\pi } = \boldsymbol {I} - \boldsymbol {\nu }\otimes \boldsymbol {\nu }$ and $\phi$ is an additional variable. The choice of $\phi = P$ results in the usual pressure gradient term $\textrm{div}\boldsymbol {\Sigma }_{P} = -\nabla P$, which allows actions resulting from pressure differences in normal direction, whereas for $\phi = 0$ actions in normal direction resulting from pressure differences are omitted, leading to $\textrm{div}\boldsymbol {\Sigma }_{P} = -\boldsymbol {\pi }\nabla P - P\mathcal {H}\boldsymbol {\nu }$ with mean curvature $\mathcal {H} = \textrm{tr}\mathcal {B}$ and $\mathcal {B} = -\nabla _{\!\mathcal {S}}\boldsymbol {\nu }$, the Weingarten mapping with covariant derivative $\nabla _{\!\mathcal {S}}$. We consider a surface observer parametrization $\boldsymbol {X}$ and a thin film observer parametrization $\boldsymbol {X}_{\xi }$ for which the surface observer velocity reads $\partial _t\boldsymbol {X} = \partial _t\boldsymbol {X}_{\xi }|_{\mathcal {S}} = \mathcal {V}\boldsymbol {\nu }$, with $\mathcal {V}$ the normal velocity of the surface. This means that the surface observer velocity and the surface velocity $\boldsymbol {V}|_{\mathcal {S}}$ differ by the surface tangential velocity $\boldsymbol {v}$. This corresponds to an Eulerian description in the tangential space and a Lagrangian description in the normal direction. With slight modifications of the analysis in Nitschke, Reuther & Voigt (2019a), we obtain as the thin film limit $\xi \to 0$ the surface Navier–Stokes equations for tangential and normal components of the surface velocity $\boldsymbol {u} = \boldsymbol {v} + \mathcal {V}\boldsymbol {\nu }$ and surface pressure $p$.

(2.3)\begin{gather} \boldsymbol{\pi}_\mathcal{S}\partial_t\boldsymbol{v}\,{+}\, \boldsymbol{\nabla}_{\!\!\boldsymbol{v}}{\boldsymbol{v}}\,{-}\, \mathcal{V}(\mathcal{B}\boldsymbol{v}\,{+}\, \nabla_{\!\mathcal{S}}\mathcal{V})\,{=}\,-\nabla_{\!\mathcal{S}}p\,{+}\,\frac{1}{{Re}} ( \boldsymbol{\Delta}^{B}\boldsymbol{v}\,{+}\, \mathcal{K}\boldsymbol{v}\,{+}\, \mathcal{H}\nabla_{\!\mathcal{S}}\mathcal{V}\,{-}\, \mathcal{V}\nabla_{\!\mathcal{S}}\mathcal{H}\,{-}\, 2\mathcal{B}\nabla_{\!\mathcal{S}}\mathcal{V} ), \end{gather}

(2.4)\begin{gather} \textrm{div}_{\!\mathcal{S}}\boldsymbol{v} = \mathcal{V}\mathcal{H} , \end{gather}

(2.5)\begin{gather} \partial_t \mathcal{V} + 2\boldsymbol{\nabla}_{\!\!\boldsymbol{v}}{\mathcal{V}} + \left\langle \mathcal{B}\boldsymbol{v} , \boldsymbol{v} \right\rangle_{\boldsymbol{g}} = - \partial_{\boldsymbol{\nu}}\phi|_{\mathcal{S}} - \mathcal{H}( p - \phi|_{\mathcal{S}}) + \frac{2}{{Re}} ( \left\langle \mathcal{B} , \nabla_{\!\mathcal{S}}\boldsymbol{v} \right\rangle_{\boldsymbol{g}} - \mathcal{V}\|\mathcal{B}\|^2 ), \end{gather}

with $\boldsymbol {\pi }_\mathcal {S}$ the surface projection, $\boldsymbol {\nabla }_{\!\!\boldsymbol {v}}{}$ the directional derivative, $\boldsymbol {\Delta }^{B}$ the Bochner Laplacian, $\mathcal {K} = \textrm{det}\mathcal {B}$ the Gaussian curvature and $\partial _{\boldsymbol {\nu }}\phi |_{\mathcal {S}}$ the normal derivative of $\phi$ in the thin film evaluated on the surface $\mathcal {S}$. The time derivative $\partial _t \boldsymbol {v}$ needs to be interpreted through an extension of $\boldsymbol {v}$ off the surface. Note that the left-hand sides of (2.3) and (2.5) are the tangential and normal components of the material derivatives. Equations (2.3) and (2.4) are independent of $\phi$. For given $\mathcal {V}$, these equations have also been previously derived by various approaches (see Arroyo & DeSimone (2009) (with corrected acceleration term Yavari, Ozakin & Sadik 2016) and Koba, Liu & Giga (2017), Jankuhn, Olshanskii & Reusken (2018), Miura (2018), Nitschke et al. (2019a)). For $\phi = 0$, also (2.5) is the same as the equation derived in Koba et al. (2017), Jankuhn et al. (2018). The surface Navier–Stokes equations (2.3)–(2.5) nicely show the tight coupling between $\boldsymbol {v}$ and $\mathcal {V}$ in the presence of curvature. Most prominently, (2.4) forces any shape change to be accompanied by a tangent flow and the rate-of-deformation tensor $\boldsymbol{\mathsf{d}} = (\nabla _{\!\mathcal {S}}\boldsymbol {v} + (\nabla _{\!\mathcal {S}}\boldsymbol {v})^T )/{2} - \mathcal {V}\mathcal {B}$, with

(2.6)\begin{align} 2\textrm{div}_{\!\mathcal{S}}\boldsymbol{\mathsf{d}} &= \boldsymbol{\Delta}^{B}\boldsymbol{v} + \mathcal{K}\boldsymbol{v} + \nabla_{\!\mathcal{S}}(\mathcal{V}\mathcal{H}) - 2\textrm{div}_{\!\mathcal{S}}(\mathcal{V} \mathcal{B}) \nonumber\\ &= \boldsymbol{\Delta}^{B}\boldsymbol{v} + \mathcal{K}\boldsymbol{v} + \mathcal{H}\nabla_{\!\mathcal{S}}\mathcal{V} - \mathcal{V}\nabla_{\!\mathcal{S}}\mathcal{H} - 2\mathcal{B}\nabla_{\!\mathcal{S}}\mathcal{V} , \end{align}

in (2.3) and (2.5) forces the surface to deform due to tangential flows. However, additional coupling terms are also present in the inertial terms.

Equations (2.3)–(2.5) assume fluid-like behaviour in tangential and normal directions and can also be written in a more compact formulation for $\boldsymbol {u}$. For $\phi = 0$, it reads

(2.7)\begin{gather} \partial_t\boldsymbol{u} + \boldsymbol{\nabla}_{\!\!\boldsymbol{u}}{\boldsymbol{u}} = -\nabla_{\!\mathcal{S}}p + \frac{2}{{Re}}\textrm{div}_{\!\mathcal{S}}\boldsymbol{\mathsf{d}} + p\mathcal{H}\boldsymbol{\nu} , \end{gather}

(2.8)\begin{gather}\textrm{div}_{\!\mathcal{S}}\boldsymbol{u} = 0 \end{gather}

(see Jankuhn et al. (2018)). However, this formulation hides the tight interplay between tangential and normal velocity components and is thus less suited to explore the resulting phenomena. Numerical approaches for (2.3)–(2.5) or (2.7) and (2.8) only exist for special cases. Most work, including also numerical analysis, is concerned with the Stokes limit on stationary surfaces $\mathcal {V} = 0$ (see e.g. Olshanskii et al. (2018), Reusken (2020)). For the surface Navier–Stokes equations in this situation see, for example, Nitschke et al. (2012), Reuther & Voigt (2018b), Fries (2018) and for their extension to evolving surfaces with prescribed $\mathcal {V}$ see Reuther & Voigt (2015), Reuther & Voigt (2018a), Nitschke et al. (2019a). The Stokes limit of (2.3)–(2.5) or (2.7) and (2.8) corresponds to the classical model of Scriven (1960) and resamples, if coupled with bulk flow, with the Boussinesq–Scriven boundary condition in multiphase flow problems (see e.g. Barrett, Garcke & Nürnberg (2015a,b)).

We are only concerned with surface phenomena but are interested in an extended model, which in addition accounts for solid-like properties in normal direction. Such solid–fluid duality of fluid deformable surfaces is considered by supplementing the evolution equations with the contribution from a Helfrich energy $({1}/{{Be}}) \int _{\mathcal {S}}(\mathcal {H} - \mathcal {H}_0)^2\,\text {d}\mathcal {S}$ to account for bending forces (Helfrich 1973), with ${Be}$ the bending capillary number and $\mathcal {H}_0$ the spontaneous curvature. We will here only consider the case $\mathcal {H}_0 = 0$. Within the Stokes limit, the resulting equations have been derived in Arroyo & DeSimone (2009), Salbreux & Jülicher (2017), Torres-Sanchez et al. (2019) and are numerically solved for simply-connected and axisymmetric surfaces in Torres-Sanchez et al. (2019) and Arroyo & DeSimone (2009), Mietke et al. (2019), respectively. Barrett et al. (2015a,b) consider these equations coupled with bulk flow. We will consider the full surface Navier–Stokes equations and provide a numerical approach for general surfaces (not necessarily simply-connected). Equations (2.3) and (2.4) are not affected by the considered extensions but (2.5) changes for $\phi = 0$ to

(2.9)\begin{align} \partial_t\mathcal{V} + 2\boldsymbol{\nabla}_{\!\!\boldsymbol{v}}{\mathcal{V}} + \left\langle \mathcal{B}\boldsymbol{v} , \boldsymbol{v} \right\rangle_{\boldsymbol{g}} &= -p\mathcal{H} + \frac{2}{{Re}}( \left\langle \mathcal{B} , \nabla_{\!\mathcal{S}}\boldsymbol{v} \right\rangle_{\boldsymbol{g}} - \mathcal{V}\|\mathcal{B}\|^2) \nonumber\\ &\quad + \frac{1}{{Be}} \left(-\Delta_{\mathcal{S}}\mathcal{H} - \frac{1}{2}\mathcal{H}^3 + 2\mathcal{H}\mathcal{K}\right), \end{align}

with $\Delta _{\mathcal {S}}$ Laplace-Beltrami operator.

3. Numerical approach

To numerically solve (2.3), (2.4) and (2.9), we consider a semi-implicit Euler time stepping scheme, a Chorin-like projection approach for (2.3) and (2.4), similar to Reuther & Voigt (2018b), Nitschke et al. (2019a), evolution of geometric quantities and the generic finite element approach proposed in Nestler, Nitschke & Voigt (2019). The latter is based on a reformulation of all operators and quantities in Cartesian coordinates and penalization of normal components. Other applications of this approach can be found in, for example, Nestler et al. (2018), Jankuhn et al. (2018), Olshanskii et al. (2018), Nitschke et al. (2018), Groß et al. (2018), Hansbo, Larson & Larsson (2020) for stationary and Nitschke, Reuther & Voigt (2019b) for evolving surfaces.

3.1. Time discretization

Let $0 = t^0 < t^1 < t^2 < \cdots$ be a partition of the time with time step width $\tau ^m:=t^{m}-t^{m-1}$. Each variable/quantity with a superscript index $m$ corresponds to the respective variable/quantity at time $t^m$. The overall algorithm for (2.3), (2.4) and (2.9) reads as follows: for $m=1,2,\dots$ do

1. (i) Move the geometry according to $\partial _t\boldsymbol {X} = \mathcal {V}\boldsymbol {\nu }$, which reads in the time-discrete setting

   (3.1)\begin{equation} \boldsymbol{X}^{m} = \boldsymbol{X}^{m-1} + \tau^m\mathcal{V}^{m-1}\boldsymbol{\nu}^{m-1}, \end{equation}
   with the parametrization of the initial geometry $\boldsymbol {X}^{0}$ and corresponding initial normal vector $\boldsymbol {\nu }^{0}$.

2. (ii) Update the normal vector according to $\partial _t\boldsymbol {\nu } = -\nabla _{\!\mathcal {S}}\mathcal {V}$ (see Huisken (1984), Kovacs, Li & Lubich (2019)). This reads in the time-discrete setting

   (3.2)\begin{equation} \boldsymbol{\nu}^{m} = \boldsymbol{\nu}^{m-1} - \tau^m\nabla_{\!\mathcal{S}}\mathcal{V}^{m-1}. \end{equation}

3. (iii) Update all other geometric quantities, e.g. the mean curvature $\mathcal {H}^{m}$, the Gaussian curvature $\mathcal {K}^{m}$, the projection $\boldsymbol {\pi }_\mathcal {S}^{m}$ and the shape operator $\mathcal {B}^{m}$, by computing derivatives of the normal vector $\boldsymbol {\nu }^{m}$. For convergence tests of this approach, we refer to Nitschke et al. (2019b).

4. (iv) Solve for normal velocity $\mathcal {V}^m$, intermediate tangential velocity $\boldsymbol {v}^\star$ and pressure $p^m$

   (3.3)\begin{equation} {d}_{\mathcal{V}}^{m} + 2\boldsymbol{\nabla}_{\!\!\boldsymbol{v}^{m}}{\mathcal{V}^m} = - \frac{2}{{Re}}\mathcal{V}^m\|\mathcal{B}^m\|^2 + g^m, \end{equation}
   (3.4)\begin{equation}\boldsymbol{d}_{\boldsymbol{v}}^* + \boldsymbol{\nabla}_{\!\!\boldsymbol{v}^{*}}{\boldsymbol{v}^\star} - \mathcal{V}^m\mathcal{B}^m\boldsymbol{v}^\star = \frac{1}{{Re}} ( \boldsymbol{\Delta}^{B}\boldsymbol{v}^\star + \mathcal{K}^m\boldsymbol{v}^\star) + \boldsymbol{f}^m , \end{equation}
   (3.5)\begin{equation}-\tau^m \Delta_{\mathcal{S}}p^m = -\textrm{div}_{\!\mathcal{S}}\boldsymbol{v}^\star + \mathcal{V}^m\mathcal{H}^m, \end{equation}
   with discrete time-derivatives ${d}_{\mathcal {V}}^{m} = ({1}/{\tau ^m}) (\mathcal {V}^m - \mathcal {V}^{m-1})$ and $\boldsymbol {d}_{\boldsymbol {v}}^* = ({1}/{\tau ^m}) \boldsymbol {\pi }_\mathcal {S}^m (\boldsymbol {v}^\star - \boldsymbol {v}^{m-1})$ and coupling terms $g^m = - p^{m}\mathcal {H}^m + ({2}/{{Re}})\left \langle \mathcal {B}^m , \nabla _{\!\mathcal {S}} \boldsymbol {v}^{m} \right \rangle _{\boldsymbol {g}} + ({1}/{{Be}}) (-\Delta _{\mathcal {S}} \mathcal {H}^m - (\mathcal {H}^m)^3/{2} + 2\mathcal {H}^m\mathcal {K}^m) - \left \langle \mathcal {B}^m\boldsymbol {v}^{m} , \boldsymbol {v}^{m} \right \rangle _{\boldsymbol {g}}$, and $\boldsymbol {f}^m = ({1}/{{Re}})( \mathcal {H}^m\nabla _{\!\mathcal {S}}\mathcal {V}^m - \mathcal {V}^m\nabla _{\!\mathcal {S}} \mathcal {H}^m - 2\mathcal {B}^m\nabla _{\!\mathcal {S}}\mathcal {V}^m ) + \mathcal {V}^m\nabla _{\!\mathcal {S}}\mathcal {V}^m$. We linearize all nonlinear terms in $p^m$, $\mathcal {V}^m$ and $\boldsymbol {v}^\star$ around the solutions at $t^{m-1}$, e.g. $\left \langle \mathcal {B}^m\boldsymbol {v}^{m} , \boldsymbol {v}^{m} \right \rangle _{\boldsymbol {g}} = \left \langle \mathcal {B}^m\boldsymbol {v}^{m} , \boldsymbol {v}^{m-1} \right \rangle _{\boldsymbol {g}} + \left \langle \mathcal {B}^m\boldsymbol {v}^{m-1} , \boldsymbol {v}^{m} \right \rangle _{\boldsymbol {g}} - \left \langle \mathcal {B}^m\boldsymbol {v}^{m-1} , \boldsymbol {v}^{m-1} \right \rangle _{\boldsymbol {g}}$, and the tangential velocity $\boldsymbol {v}^{m}$ follows from (3.6).

5. (v) Update tangential velocity $\boldsymbol {v}^m$

   (3.6)\begin{equation} \boldsymbol{v}^m = \boldsymbol{v}^\star - \tau^m\nabla_{\!\mathcal{S}}p^m. \end{equation}

3.2. Space discretization

The remaining step is to discretize (3.3)–(3.5) from the above algorithm in space by using either the generic surface finite element method for tensor-valued surface partial differential equations (PDEs) proposed in Nestler et al. (2019) or the surface finite element method for scalar-valued surface PDEs from Dziuk & Elliott (2013). Let $\mathcal {S}_h=\mathcal {S}_h(t)|_{t=t^m}$ be an interpolation of the surface $\mathcal {S}=\mathcal {S}(t)|_{t=t^m}$ at time $t^m$ such that $\mathcal {S}_h := \bigcup _{T\in \mathcal {T}}T$, where $\mathcal {T}$ denotes a conforming triangulation. Furthermore, the finite element space is introduced as $\mathbb {V}(\mathcal {S}_h) := \lbrace v\in \mathcal {C}^0(\mathcal {S}_h) : v|_T\in \mathcal {P}^1(T), \forall v\in \mathcal {T} \rbrace$ with $\mathcal {C}^k(\mathcal {S}_h)$ the space of $k$-times continuously differentiable functions on $\mathcal {S}_h$ and $\mathcal {P}^l(T)$ polynomials of degree $l$ on the triangle $T\in \mathcal {T}$. We use the finite element space $\mathbb {V}(\mathcal {S}_h)$ twice as trial and as test space and additionally introduce the $L_2$ inner product on $\mathcal {S}_h$, as $(a,b) := \int _{\mathcal {S}_h}\langle a, b\rangle dS$. Thus, the finite element approximations of (3.3)–(3.5) read: find $\mathcal {V}^{m}\in \mathbb {V}(\mathcal {S}_h)$, $\boldsymbol {v}^{m}\in \mathbb {V}(\mathcal {S}_h)^3$ and $p^{m} \in \mathbb {V}(\mathcal {S}_h)$ such that $\forall \alpha \in \mathbb {V}(\mathcal {S}_h)$, $\boldsymbol {\alpha } \in \mathbb {V}(\mathcal {S}_h)^3$ and $\beta \in \mathbb {V}(\mathcal {S}_h)$

(3.7)\begin{align} ({d}_{\mathcal{V}}^{m},\alpha) + 2 (\boldsymbol{\nabla}_{\!\!\boldsymbol{v}^{m}}{\mathcal{V}^m},\alpha) &= -\frac{2}{{Re}}\left(\|\mathcal{B}^m\|^2\mathcal{V}^m + g^m + \frac{\omega_{{a}}}{A_0^2}(A^{m}-A_0) \mathcal{H}^{m}, \alpha\right) \nonumber\\ &\quad + (D(\nabla_{\!\mathcal{S}} \mathcal{V}^{m-1} - \nabla_{\!\mathcal{S}} \mathcal{V}^{m}), \nabla_{\!\mathcal{S}} \alpha), \end{align}

(3.8)\begin{align} (\boldsymbol{d}_{\boldsymbol{v}}^{*}, \boldsymbol{\alpha}) + (\boldsymbol{\nabla}_{\!\!\boldsymbol{v}^{*}}{\boldsymbol{v}^\star},\boldsymbol{\alpha}) - (\mathcal{V}^m\mathcal{B}^m\boldsymbol{v}^\star,\boldsymbol{\alpha}) &= -\frac{1}{{Re}} (\nabla_{\!\mathcal{S}}\boldsymbol{v}^\star, \nabla_{\!\mathcal{S}}\boldsymbol{\alpha}) + \frac{1}{{Re}} (\mathcal{K}^m\boldsymbol{v}^\star, \boldsymbol{\alpha}), \nonumber\\ &\quad + (\,\boldsymbol{f}^m + \omega_{\boldsymbol{t}} (\boldsymbol{v}^\star \boldsymbol{\cdot} \boldsymbol{\nu}^m) \boldsymbol{\nu}^m, \boldsymbol{\alpha}), \end{align}

(3.9)\begin{equation} \tau^m (\nabla_{\!\mathcal{S}}p^m , \nabla_{\!\mathcal{S}}\beta) = (\boldsymbol{v}^\star, \nabla_{\!\mathcal{S}}\beta) + (\mathcal{V}^m\mathcal{H}^m, \beta). \end{equation}

Equation (3.7) is stabilized by artificial diffusion with coefficient $D$, following ideas of Smereka (2003) for surface diffusion. As in Nitschke et al. (2019b), an additional term for global surface area conservation is included with a penalty parameter $\omega _{{a}}>0$, initial surface area $A_0$ and actual surface area $A^m$. In (3.8), we have used the same symbols for the extended tangential velocity field $\boldsymbol {v} = (\boldsymbol {v}_{x} \boldsymbol {e}^x, \boldsymbol {v}_{y} \boldsymbol {e}^y, \boldsymbol {v}_{z} \boldsymbol {e}^z)$ and the extended operators (see Reuther & Voigt (2018b)). We further use $\textrm{div}_{\!\mathcal {S}} \boldsymbol {v} = \nabla \boldsymbol{\cdot} \boldsymbol {v} - \boldsymbol {\nu } \boldsymbol{\cdot} (\nabla \boldsymbol {v} \boldsymbol{\cdot} \boldsymbol {\nu })$ and introduce the additional term $(\omega _{\boldsymbol {t}} (\boldsymbol {v}^\star \boldsymbol{\cdot} \boldsymbol {\nu }^m) \boldsymbol {\nu }^m,\boldsymbol {\alpha })$, with $\omega _{\boldsymbol {t}} > 0$, to penalize normal components of the extended tangential velocity. For convergence studies in $\omega _{\boldsymbol {t}}$, we refer to Nestler et al. (2018). The resulting equations for the components $\boldsymbol {v}^\star _{x}, \boldsymbol {v}^\star _{y}$ and $\boldsymbol {v}^\star _{z}$ are solved by surface finite elements. From these fields, $\boldsymbol {v}^m$ can be computed. For more details, especially for evaluating the local inner products in the $L_2$ inner products for the extended tangential velocity, we refer to Nestler et al. (2019). All equations are solved using the adaptive finite element toolbox AMDiS (Vey & Voigt 2007; Witkowski et al. 2015). The software is open-source and the source-code used to produce the following simulation results is provided in Reuther (2020).

4. Simulation results

All examples are chosen to demonstrate the tight coupling between $\boldsymbol {v}$ and $\mathcal {V}$ in the presence of curvature. The first considers a perturbed sphere with zero velocity as the initial condition. The Helfrich term induces a normal velocity, which generates tangential flow. The final configuration is a sphere with zero velocity. The second example considers a rotating Killing vector field on a sphere as the initial condition. The tangential flow induces a normal velocity and with it dissipation. The final configuration is again a sphere with zero velocity. We compare the dynamics of both examples with respect to Reynolds number ${Re}$ and bending capillary number ${Be}$. We further consider convergence studies in mesh size $h$ and time step width $\tau$. Besides several coarse-grained measures, such as energy components and eccentricity, (2.4) is used for convergence studies. It provides a severe measure for the accuracy of the algorithm, as it is never used in the approach and contains the tangential and normal parts of the velocity, $\boldsymbol {v}$ and $\mathcal {V}$, and the geometric quantity $\mathcal {H}$. In the following simulations, we use $D=62.5$, $\omega _{\boldsymbol {t}}=1 \times 10^{5}$ and $\omega _{{a}}=1\times 10^{3}$. We further use an analytic form for $\boldsymbol {\nu }^{0}$ and have chosen the examples such that mesh distortions do not alter the simulations’ results. More extreme examples will require a redistribution of mesh points (see e.g. Mikula et al. (2014)).

4.1. Relaxation of perturbed sphere

Let $\boldsymbol {X}_{S}(\phi , \theta )$ be the standard parametrization of the unit sphere with standard parametrization angles $\phi ,\theta$. We use the parametrization $\boldsymbol {X}(\phi , \theta ) = r(\phi , \theta )\boldsymbol {X}_{S}(\phi , \theta )$ with a space-dependent radius $r(\phi , \theta ) = 1 + r_0\cos (\phi )\sin (3\theta )$. Figure 1 shows the evolution for $r_0 = 0.4$ and zero initial velocity. The dynamics of the induced tangential flow field and shape changes are clearly visible. The correspondence between $\boldsymbol {v}$ and $\mathcal {V}$ is further highlighted in kinetic energy plots, with a strong increase in normal kinetic energy and an induced but delayed response of the tangent kinetic energy at the beginning. The later relaxation towards a sphere corresponds to a more intermediate coupling. The results correspond to ${Re} = 1$ and ${Be} = 2$, and the simulations are performed with $h = 4.68 \times 10^{-2}$ and $\tau = 4.9\times 10^{-3}$. The dependency on ${Re}$ (${Be} = 2$) and ${Be}$ (${Re} = 1$) is considered in figure 2. The strongest oscillations are observed for large ${Re}$ and small ${Be}$. However, also for small ${Re}$, the dynamics significantly differs from pure Helfrich flow, which is shown for comparison.

Figure 1. (a) Relaxation of a perturbed sphere for $t = 0$, $0.25$, $0.5$, $0.75$, $1$, $1.25$, $1.75$, $2$, $2.5$, $7.5$ (left to right, top to bottom); the tangential flow field is visualized by line integral convolution (LIC). Tangent, normal and overall kinetic energy $E_{{T}} := (\int _\mathcal {S}\langle \boldsymbol {v},\boldsymbol {v}\rangle \,\text {d}\mathcal {S})/{2}$, $E_{{N}} := (\int _\mathcal {S}\mathcal {V}^2\,\text {d}\mathcal {S})/{2}$ and $E_{{T}}+E_{{N}}$, respectively, against time $t$ (b) and tangent kinetic energy $E_{{T}}$ against normal kinetic energy $E_{{N}}$ (c).

Figure 2. Deviation from a sphere $\sigma _{{S}}$ against time $t$ for different Reynolds numbers ${Re}$ (a,b) and Bending capillary numbers ${Be}$ (c,d) for the perturbed sphere simulation (a,c) and power spectrum of the normalized deviation from a sphere $\sigma _{{S}} / \langle \sigma _{{S}}\rangle$ (b,d) with $\langle \sigma _{{S}}\rangle$ the time average of $\sigma _{{S}}$. Pure Helfrich flow is shown for comparison as a dashed line and is indicated as ‘no flow’.

4.2. Killing vector field

The initial tangential velocity on the unit sphere is given by the Killing vector field $\boldsymbol {v}|_{t=0}=(-z,0,x)^T$ with $\boldsymbol {X}=(x,y,z)^T\in \mathcal {S}$. The tangential velocity induces deformations towards ellipsoidal-like shapes. Due to the induced normal velocity, energy dissipates. Theoretically, a force balance with the bending forces of the Helfrich energy can be established. Using the axisymmetric setting, an ordinary differential equation for these meta-stable states can be derived. However, these states can never be reached during evolution. The shape instead overshoots, oscillates and further dissipates energy, which decreases the driving force and leads to a relaxation back to a sphere with zero tangential velocity (see figure 3). The results correspond to ${Re} = 1$ and ${Be} = 2$, and the simulations are performed with $h = 4.68\times 10^{-2}$ and $\tau = 4.9\times 10^{-3}$. Convergence studies with respect to $h$ and $\tau$ are considered, indicating almost second-order convergence in $h$ and first-order in $\tau$. However, also number, time and strength of shape oscillations change with refinement (see figure 3e). The dependency of the dynamics on ${Re}$ (${Be} = 2$) and ${Be}$ (${Re} = 1$) is shown in figure 4, again with $h = 4.68\times 10^{-2}$ and $\tau = 4.9\times 10^{-3}$. The results are similar to figure 2.

Figure 3. (a) Relaxation of Killing vector field for $t = 0$, $1.5$, $3.5$, $15$, $20$ (left to right); the tangential flow field is visualized by LIC. (b,c) Contour plot of the sliced geometry in the time interval $[0,10]$ with ascending grey scale indicating increasing time (b) and plot of the $x$-/$y$-coordinate of the geometry against time $t$ (c). (d,e) Experimental order of convergence (EOC) for different mesh sizes $h$ and time step widths $\tau$ with a constant ratio $h^2/\tau$ and the measure of the error $e := \|\textrm{div}_{\!\mathcal {S}}\boldsymbol {v} - \mathcal {V}\mathcal {H}\|_2$ (d). Deviation from a sphere $\sigma _{{S}}$ against time $t$ with $\sigma _{{S}} := \int _{\mathcal {S}}(\mathcal {H} - \mathcal {H}_{{S}})^2\,\text {d}\mathcal {S}$ and $\mathcal {H}_{{S}}$ the mean curvature of a sphere with equal surface area for different mesh sizes $h$ (e).

Figure 4. Deviation from a sphere $\sigma _{{S}}$ against time $t$ for different Reynolds numbers ${Re}$ (a,b) and Bending capillary numbers ${Be}$ (c,d) for the Killing vector field relaxation (a,c) and power spectrum of the normalized deviation from a sphere $\sigma _{{S}} / \langle \sigma _{{S}}\rangle$ (b,d) with $\langle \sigma _{{S}}\rangle$ the time average of $\sigma _{{S}}$.

5. Conclusion

With the considered thin film limit, we have provided a new approach to derive the governing equations of fluid deformable surfaces. They consider fluid-like behaviour in the tangential and normal directions beyond the Stokes limit and are supplemented by a Helfrich energy to model solid-like (bending) behaviour in the normal direction. The splitting of the surface velocity into tangential and normal components shows their tight interplay with geometric quantities of the surface. This is known for the rate-of-deformation tensor. However, additional coupling terms are also present in the inertial terms. The considered numerical approach to solve these equations, which combines evolution of geometric quantities with surface finite elements and a general finite element method for tangential tensor-valued surface partial differential equations, is applicable to general surfaces (not restricted to simply-connected surfaces) and shows reasonable convergence properties with respect to mesh size $h$ (second-order) and time step width $\tau$ (first-order). The computational examples are chosen to demonstrate the coupling between tangential and normal velocities, where in the presence of curvature any shape change is accompanied by a tangential flow and, vice versa, the surface deforms due to tangential flow. The dynamics of the relaxation strongly depends on the fluid and solid properties. However, the simulations also show that Killing vector fields are only possible as meta-stable states, in situations where the viscous force is balanced by the bending force. The only possible stable stationary state in the considered setting is a sphere with zero velocity.

The computational examples can provide benchmark problems for other numerical approaches, which can be extended to the considered model, e.g. Nitschke et al. (2017), Olshanskii et al. (2018), Torres-Sanchez, Santos-Olivan & Arroyo (2020), Lederer et al. (2019). They also form the basis for more complex models, which include coupling with concentration fields for proteins and dependency of $\mathcal {H}_0$ on concentration in lipid bilayers, or coupling with liquid crystal theory as in Nitschke et al. (2019a) for Erickson-Leslie type models or with Landau-de Gennes theory on surfaces (Nitschke et al. 2019b) for Beris-Edwards type models, which also can be extended by active contributions to model, e.g. phenomena as considered in Keber et al. (2014). However, any quantitative comparison in these applications will require to also consider the surrounding bulk phases, as, for example, considered in Barrett et al. (2015a), Barrett et al. (2015b), or at least a constraint for the enclosed volume. Even if the approach is applicable for general surfaces, it cannot handle topological changes. This would require a reformulation of the equations using an implicit description, e.g. the diffuse interface approach (Rätz & Voigt 2006).