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.


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.

Mathematical modelling
We start from a slightly more general Navier-Stokes equation in the thin film Ω ξ (t) = S(t) × [−ξ/2, ξ/2] ⊂ R 3 with a regular evolving surface S(t), film thickness ξ and surface normal ν.In the Eulerian description, it reads where P is the pressure, π = I − ν ⊗ ν and φ is an additional variable.The choice of φ = P results in the usual pressure gradient term divΣ P = −∇P, which allows actions resulting from pressure differences in normal direction, whereas for φ = 0 actions in normal direction resulting from pressure differences are omitted, leading to divΣ P = −π ∇P − PHν with mean curvature H = trB and B = −∇ S ν, the Weingarten mapping with covariant derivative ∇ S .We consider a surface observer parametrization X and a thin film observer parametrization X ξ for which the surface observer velocity reads ∂ t X = ∂ t X ξ | S = Vν, with V the normal velocity of the surface.This means that the surface observer velocity and the surface velocity V | S differ by the surface tangential velocity 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 ξ → 0 the surface Navier-Stokes equations for tangential and normal components of the surface velocity u = v + Vν and surface pressure p.
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 u.For φ = 0, it reads 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 V = 0 (see e. duality of fluid deformable surfaces is considered by supplementing the evolution equations with the contribution from a Helfrich energy (1/Be) S (H − H 0 ) 2 dS to account for bending forces (Helfrich 1973), with Be the bending capillary number and H 0 the spontaneous curvature.We will here only consider the case 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 φ = 0 to with Δ S Laplace-Beltrami operator.

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)

Time discretization
be a partition of the time with time step width τ 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, . . .do (i) Move the geometry according to ∂ t X = Vν, which reads in the time-discrete setting with the parametrization of the initial geometry X 0 and corresponding initial normal vector ν 0 .(ii) Update the normal vector according to ∂ t ν = −∇ S V (see Huisken (1984), Kovacs, Li & Lubich (2019)).This reads in the time-discrete setting g , and We linearize all nonlinear terms in p m , V m and v around the solutions at t m−1 , e.g. (3.6)

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 S h = S h (t)| t=t m be an interpolation of the surface S = S(t)| t=t m at time t m such that S h := T∈T T, where T denotes a conforming triangulation.Furthermore, the finite element space is introduced as V(S h ) := {v ∈ C 0 (S h ) : v| T ∈ P 1 (T), ∀v ∈ T } with C k (S h ) the space of k-times continuously differentiable functions on S h and P l (T) polynomials of degree l on the triangle T ∈ T .We use the finite element space V(S h ) twice as trial and as test space and additionally introduce the L 2 inner product on S h , as (a, b) := S h a, b dS.Thus, the finite element approximations of (3.3)-(3.5)read: find (3.9)  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 ω 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 v = (v x e x , v y e y , v z e z ) and the extended operators (see Reuther & Voigt (2018b)).We further use div with ω t > 0, to penalize normal components of the extended tangential velocity.For convergence studies in ω t , we refer to Nestler et al. (2018).The resulting equations for the components v x , v y and v z are solved by surface finite elements.From these fields, 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).

Simulation results
All examples are chosen to demonstrate the tight coupling between v and 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 τ .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, v and V, and the geometric quantity H.In the following simulations, we use D = 62.5, ω t = 1 × 10 5 and ω a = 1 × 10 3 .We further use an analytic form for ν 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)).

Relaxation of perturbed sphere
Let X S (φ, θ ) be the standard parametrization of the unit sphere with standard parametrization angles φ, θ.We use the parametrization X (φ, θ ) = r(φ, θ )X S (φ, θ ) with a space-dependent radius r(φ, θ ) = 1 + r 0 cos(φ) sin(3θ).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 v and 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 × 10 −2 and τ = 4.9 × 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.

Killing vector field
The initial tangential velocity on the unit sphere is given by the Killing vector field v| t=0 = (−z, 0, x) T with X = (x, y, z) T ∈ 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 × 10 −2 and τ = 4.9 × 10 −3 .Convergence studies with respect to h and τ are considered, indicating almost second-order convergence in h and first-order in τ .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 × 10 −2 and τ = 4.9 × 10 −3 .The results are similar to figure 2.

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 τ (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 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) 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 V seeReuther & 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 ofScriven (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 900 R8-3 https://doi.org/10.1017/jfm.2020.

900
FIGURE 2. Deviation from a sphere σ 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 σ S / σ S (b,d) with σ S the time average of σ S .Pure Helfrich flow is shown for comparison as a dashed line and is indicated as 'no flow'.

900
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 t (c).(d,e) Experimental order of convergence (EOC) for different mesh sizes h and time step widths τ with a constant ratio h 2 /τ and the measure of the error e := div S v − VH 2 (d).Deviation from a sphere σ S against time t with σ S := S (H − H S ) 2 dS and H S the mean curvature of a sphere with equal surface area for different mesh sizes h (e).