An ALE method for simulations of axisymmetric elastic surfaces in flow

The dynamics of membranes, shells and capsules in fluid flow has become an active research area in computational physics and computational biology. The small thickness of these elastic materials enables their efficient approximation as a hypersurface which exhibits an elastic response to in-plane stretching and out-of-plane bending, possibly accompanied by a surface tension force. In this work, we present a novel ALE method to simulate such elastic surfaces immersed in Navier-Stokes fluids. The method combines high accuracy with computational efficiency, since the grid is matched to the elastic surface and can therefore be resolved with relatively few grid points, and the axisymmetric setting reduces the system effectively to a two-dimensional problem. We show in several numerical test cases that accurate results can be achieved at computational times on the order of minutes on a single core CPU. We further illustrate that the method can be easily extended to account for very different physics in both domains surrounding the surface. As an example, we present first simulations of the observed shape oscillations of novel microswimming shells and the uniaxial compression of biological cells filled with cytoplasm during an AFM experiment.


Introduction
The dynamics of elastic membranes, shells and capsules in fluid flow has become an active research area in computational physics and computational biology. Example systems include vesicle membranes immersed in fluids [1,2,3,4], red and white blood cells transported with the blood plasma [5,6,7,8], general biological cells including cytoplasmic flows [9,10], or even man-made elastic thin shells to deliver cargo though a fluid [11]. Therefore, there exists a broad interest in advanced modeling and simulation technologies that enable the understanding of such systems.
In all these examples, the elastic material is typically very thin such that its direct numerical resolution with continuum-based discretization techniques requires prohibitively fine mesh sizes [12]. This problem can be circumvented by replacing volumes of thin layers by dimensionally reduced surfaces, i.e., the elastic material is approximated as a hypersurface of zero off one of the fluid phases, which offers a further advantage w.r.t. Immersed Boundary Methods. Similarly, the method can be easily extended to account for completely different physics in both surrounding domains. We illustrate this by simulating the uniaxial compression of a biological cell filled with intracellular fluid and the propulsion of novel microswimmers having compressible/incompressible fluids inside/outside their elastic shell.
The rest of this article is structured as follows. In Sec. 2 the governing equations for the model are introduced along with the suitable boundary conditions to ensure force balance across the membrane. Details about the membrane forces are provided in Sec. 3. The numerical realization of the problem is presented in Sec. 4 including details about time and space discretization, the weak form of the resulting finite element problem in each time step, and the ALE mesh movement in each time step, respectively. Numerical test cases are presented in Sec. 5, where three different parameter configurations have been used in order to verify the correctness of the imposed surface tension, bending stiffness and in-plane stretching on the membrane, respectively. This is followed by a time step study, a mesh resolution study, some remarks on how to stabilize the in-plane stretching force numerically, and the presentation of the previously mentioned simulations of uniaxial compression of biological cells and the propulsion of novel microswimmers. Some further implementation details can be found in the appendix.

Governing equations
Consider a (closed) surface of revolution immersed in a cylindrical 3D domain, representing an elastic membrane immersed in a fluid ( Fig. 1(a)). Half of the domains cross section acts as the two-dimensional computational domain Ω ( Fig. 1(b)). Ω is composed of two separate parts, the exterior Ω 0 (t) and the interior Ω 1 (t), both describing incompressible viscous fluids: Ω = int(Ω 0 ∪ Ω 1 ). The interface Γ(t) = Ω 0 (t) ∩ Ω 1 (t), which separates both domain parts, corresponds physically to an elastic shell and is described by the parametric form X(s, t) = (X(s, t), R(s, t)) T , with the arc length parameter s. The x−axis is chosen to be the symmetry axis and the r−axis is the distance to the symmetry axis in the considered 2D meridian x − r domain. With the divergence operator for cylindrical coordinate systems of the form [38] 3 the Navier-Stokes equations for the external (i = 0) and internal (i = 1) fluid read with the velocity v = (v x , v r ) T in Ω ∪ Γ, pressures p i , viscosities η i , and densities ρ i in Ω i , respectively. The velocity is continuous across Γ and can hence be defined in the whole domain, while the pressure is discontinuous across Γ and has therefore to be defined separately in both phases. The following jump condition holds at the interface to ensure the force balance at the elastic membrane where n is the interface normal pointing to Ω 0 and f = f 0 − f 1 denotes the jump operator across the interface Γ. The surface force term ∂E ∂Γ is defined in the following. At the symmetry axis we specify the usual free slip condition.

Shell membrane forces
The membrane is assumed to be an isotropic thin shell with a shell thickness d. The force response of a membrane to elastic deformations can be described with in-plane and out-of-plane forces acting to minimize the corresponding elastic energies. The in-plane energy is also referred as stretching energy (E stretch ) and the out-of-plane energy as bending energy (E bend ). Additionally, a surface tension energy (E tension ) can be present. The membrane forces are then given by the first variation of these energies with respect to changes in Γ The surface tension energy, tending to minimize the membrane area, reads [9] E tension = Γ γ dA , with the material specific surface tension γ [N/m], which is assumed to be constant on the whole membrane.
Bending stiffness is the resistance of the membrane against bending (out-of-plane) deformations. The bending energy tends to minimize the deviation of the local curvature from the material's reference curvature (also termed spontaneous curvature). The bending energy reads with the material specific bending stiffness c b [Nm] and the total curvature κ = ∇ · n, which is twice the mean curvature of the membrane. The spontaneous curvature κ ref is in many practical cases either zero or the total curvature in the initial state, depending on the physical context of the problem. The stretching energy E stretch minimizes in-plane stretching and compression of the membrane compared to the reference state. In the axisymmetric setting, it is useful to describe the stretching energy in terms of the two principal stretches λ 1 and λ 2 , which provide information about relative changes of surface lengths in lateral and circumferential direction, respectively. An illustration can be found in Fig. 2(a). The principal stretches read with the subscript ref corresponding to the quantities at the same material point in the reference state.
In 3D elasticity theory, the response of an isotropic elastic body to elastic deformations can be described by two material specific parameters: Young's modulus E and the Poisson ratio ν. For a thin elastic material of thickness d, these parameters are typically reformulated into surface parameters, for example the the area dilation modulus K A and area shear modulus K S [42]. Considering a rectangular surface element, K A describes the response of the membrane to in-plane area changes with constant aspect ratio of the surface element ( Fig. 2(b)). K S provides information about the response to in-plane shear deformations with constant area of the surface element ( Fig. 2(c)). The two parameters together with c b can be calculated directly from Young's modulus, Poisson ratio and shell thickness: In terms of K A and K S , the stretching energy can be written as [9] The corresponding forces to the respective energies can be obtained by calculating the first variations of the three surface energies: where K is the Gaussian curvature, ∆ Γ the surface Laplacian, and ∇ Γ the surface gradient, respectively. Derivations can be found in [9,35].

Time discretization
The problem is discretized in time with equidistant time steps of size τ. The complete system is split into several subproblems which are solved subsequently in each time step. At first the membrane forces are computed from the current membrane shape. Afterwards, the Navier-Stokes equations are solved by an implicit Euler method using the previously computed membrane forces. The resulting velocity is then used to advect the surface and the domains.
Using the ALE approach, the material derivative of the velocity ∂ • t v = d t v+v·∇v is discretized in the n-th time step as follows where v n−1 grid is the velocity of the grid movement from the previous time step calculated with one of the mesh smoothing algorithms presented in Sec. 4.4. This term has to be subtracted from the convection term due to the mesh update. The velocity v n−1 moved is the velocity of the last time step, but after the mesh update, i.e. the grid point coordinates have been moved without changing the velocity values in each degree of freedom (DOF).

Space Discretization
We use a Finite-Element method where the grids to represent the domains are matched at the immersed interface, i.e. they share the same grid points. Accordingly, let T h,i , i = 0, 1 be the triangulations of Ω i such that T h = T h,0 ∪ T h,1 is a conforming triangulation of Ω. The triangulation of the membrane is given by Γ h = T h,0 ∩ T h,1 . An example for the corresponding numerical mesh is shown in Sec. 4.2. The triangulations T h,0 and T h,1 are separated, i.e. Γ h acts as a boundary for both. This definition of the mesh ensures the possibility of continuous velocity and non-zero pressure jump across Γ. See the appendix for further details on the technical realization.
In order to obtain the membrane force, the tangential t and the values for n, κ, K, ∆ Γ κ, λ 1 , and λ 2 have to be calculated. Consider the parametrization X(s) with the arc length parameter s (hence, X s = 1) we obtain from [24] In the discrete case, X i = X i , R i for i = 0, ..., N − 1 is the sequence of membrane grid points ordered counterclockwise, approximating Γ by piece wise linear line segments X i , X i+1 . Then, the derivatives of X i , R i and κ i with respect to s can be calculated with finite differences, here shown for X i , where i = 1, ..., N − 2 The approximations for the principal stretches read The usage of λ i 1 to compute the derivative ∂ s λ 1 turned out to be numerically unstable. Hence, instead of the vertex-based values λ i 1 we use the following values which represent the stretching at the midpoints of the two neighboring edges of vertex i With these quantities the derivatives of the principal stretches are computed, To calculate the surface quantities on the symmetry axis (e.g. i = 0 and i = N − 1), one takes advantage of the fact that R = R ss = R ssss = 0 and X s = X sss = 0 and applies l'Hospital's rule, see [24, Sec. 3.1] for details.

Weak form
The fully discrete system in weak form is presented in the following. The finite element (FE) spaces read where V h is the FE space for the velocity. The resulting velocity will be continuous across Γ h . The FE spaces M h,i for the pressures p i are defined separately in T h,i to allow discontinuous pressure across Γ h . Assuming constant viscosities η i in T h,i , the weak form of the system given in Sec. 2 reads: with surface forces ∂E · ∂Γ defined by Eqs. (12)- (20).

Mesh movement
The basic idea of the ALE approach is to move the grid of the elastic structure with the material velocity, while the fluid grid is moved with an arbitrary velocity keeping the mesh in a proper shape. Accordingly, in the present work the membrane grid points are displaced with the velocity v in every time step. As the membrane is moving, it is necessary to extend this movement to every grid point in T h to keep the mesh well in shape. Two different strategies for rearranging grid points on a mesh with given boundary movement have been used for the simulations in the present work.

Solving Laplace problem
The first strategy involves the harmonic extension of interface movement by solving the Laplace problem where v grid is the velocity of the grid points in T h . The mesh is then moved with v grid . In some cases, it can be helpful to use the Bilaplace problem instead: Both approaches can be used in most cases. However, for strong and/or periodic deformations, some problems can occur using the Laplace Problem for mesh smoothing, e.g. elements near the interface can degenerate slowly and cause a crash of the simulation. In this case, the mesh smoothing approach shown in the following may be a better choice.

Element area and length conservation
The second strategy is motivated by the fact that solving the Laplace problem Eq. (24) or the Bilaplace problem Eq. (25) can lead to degenerated elements. To avoid this, it is helpful to penalize changes in area and edge length of the elements in order to prevent large deformations in both, the area and the aspect ratio of the elements. How to preserve element surface areas or grid point distances is described in [43].
Assume that the boundary has been moved already. Iterate over all grid points x i . Check whether the length of each edge and/or the area of each element, which has x i as a vertex, has changed. If, for example, the area of an element, which has the vertex x i , decreased after the movement of Γ, x i will be moved in order to re-increase the area of this element. Since this has to be done for all elements that share x i as a vertex, a mean value of the calculated displacements is computed. The scheme for calculating the grid point movement in order to get the new grid reads: 1. Calculate areas A and edge lengths l a , l b , l c of all elements using the point coordinates.

Compute the velocity of every grid point
where shows an example for both, the Laplace smoothing approach and the area and length conservation approach. In both cases, a strong deformation has been imposed within a total amount of 500 time steps. The membrane shape at the end is circular with perfect agreement in both cases. The Laplace smoothing result ( Fig. 5(b)) seems a little more even in the internal fluid. However, the area and length conservation approach (Fig. 5(c)) produces a better result: In the external fluid, the triangles near the membrane are less deformed. Around the symmetry axis, elements have been less compressed in x direction where at the poles, they have been less stretched in r direction (the latter is also visible in the internal fluid). These rather small improvements of the second approach can be quite important, e.g. when due to a periodic deformation small errors that occur in the Laplace smoothing accumulate over time. Nevertheless, in this paper, if not mentioned otherwise, the Laplace smoothing approach is used as it is independent of additional problem-specific parameters (like c a , c l ).

Numerical Tests
Test case simulations were performed to show the quality and capability of the presented model, to verify the correctness of the interfacial forces and to analyze the mesh and time step stability.

Verification of the interfacial forces
In the following, we prescribe the elastic membrane as an initially oblate-shaped object, i.e. its cross-section as the combination of two parallel lines and a semicircle ( Fig. 7(a)). The radius of the semicircle amounts to 0. In case 1, the membrane is expected to evolve into a sphere, since surface tension tends to minimize the surface area. This behavior is verified in the simulations, see Fig. 7(b). In case 2 the stationary state is expected to yield a red blood cell like shape, since bending stiffness tends to minimize the curvature locally while the additional influence of area dilation prevents the membrane to deform (or stretch) strong enough to get spherical. In this sense, case 2 is similar to a lipid vesicle, as the finite K A leads to approximate conservation of surface area. The stationary shapes of these simulations fit qualitatively well with theory Fig. 7(c) [24]. 11 In case 3, in-plane elasticity penalizes stretching in tangential direction to the surface. Hence, the initial condition of λ 1 and λ 2 being larger than 1.0 causes a deformation of the membrane such that the radius of the oblate should decrease, where the thickness should somehow increase a bit. Note, that in the absence of volume conservation, the membrane would contract in both directions such that the principal stretches would approach 1.0 everywhere on the membrane. However, the conservation of enclosed volume prevents the principal stretches from reaching 1.0, in general. The stationary state of the stretching dominant case is shown in (Fig. 7(c)). The principal stretches in the stationary state are illustrated in Fig. 8(a). The equilibrium configuration of the elastic surface shows a significant stretch (λ 1 > 1) in lateral direction, which is necessary to accommodate the excess volume. This is accompanied by a compression in circumferential direction (λ 2 < 1) to minimize surface dilation. In the numerical model, the Navier-Stokes equations are solved with separately defined pressures in Ω 0 and Ω 1 , which leads to a discontinuous pressure field along Γ. The pressure field together with velocity glyphs is shown in Fig. 8(b) for the bending stiffness dominant case.

Mesh resolution study
Three different mesh resolutions have been chosen to investigate the dependence of the simulation results on the mesh (Fig. 5.2). In the following we denote the mesh size by h i and number of surface grid points by N i for i ∈ {1, 2, 3}. The coarsest mesh has a mesh size of h 1 = 0.055 at the interface, hence the membrane is resolved by N 1 = 23 grid points. The complete mesh (a) base mesh, h 1 (b) one refinement step, h 2 (c) two refinement steps, h 3 Figure 9: The three (initial) meshes used for the mesh resolution study. Close up views around the interface. Blue refers to Ω 0 , green refers to Ω 1 .
has 228 grid points. Refining this mesh by two triangle bisections leads to an intermediate mesh (h 2 = 0.0275, N 2 = 45, 820 total grid points). Two further bisections lead to the finest mesh (h 3 = 0.01375, N 3 = 89, 3102 total grid points). Fig. 5.2 shows the membrane points of the respective case for all chosen mesh resolutions. As can be seen, the method produces accurate results even with relatively coarse meshes. There is no visible difference between the shapes of all three meshes. To quantify that, we introduce two different error measures in the following.
The mean distance error between membrane points on the h 1 -and h 2 -mesh and the corresponding points on the next finer mesh is defined as For the h 2 (or h 3 ) mesh, every other (or 4th) membrane point is used, s.t. only corresponding grid points (existing in the coarsest mesh) are compared. As a second error indicator we use the membrane cross section perimeter. Given that membrane points are sorted, we define the perimeter where X h i N i := X h i 0 . Then, the error is calculated using With these values, the experimental order of convergence (EOC E and EOC P ) can be obtained The obtained values are shown in Eq. 5.3. The point coordinates converge with order 1 and the areas/perimeters converge with order 2.

Time step study
According to the previous section, it seems sufficient to do further studies using the h 1 -mesh. The three different cases are now being analyzed using different time step sizes. Fig. 5.3(ac) illustrates the shape evolution of the surface tension case (tested with timesteps τ = 2.5 · 10 −4 , 2.5 · 10 −5 , 2.5 · 10 −6 ). During the evolution, slight differences in the shapes can be observed. The stationary state shows very good agreement with the expected spherical shape (Fig. 11(c)), with only slight tangential differences (surface tension works in normal direction only), for larger time steps. The simulation with the largest time step required only 80 time steps to reach the stationary state. The necessary compute time of less than 1 minute on a single core CPU (Intel Haswell, 2.50 GHz) illustrates the efficiency of the proposed method.
In the bending stiffness dominant case (case 2) the shape differences are also comparatively small. Time step sizes 10 −5 , 10 −6 , 10 −7 have been tested. Fig. 5.3(d-f) shows the time evolution of the cross section shape. The agreement of point coordinates is good enough to have no visible difference between the different time step sizes. This is also observable when analyzing the membrane energy (Fig. 5.3(j,k)) Fig. 11(j) shows E membrane for the bending stiffness dominant case for all tested time step sizes. A close up view of the increase at t = 0.001 is shown in Fig. 11(k), for the largest and smallest time step. Additionally, the time step τ = 1.2 · 10 −5 is included. This value was the largest time step size, where stable simulations were possible, as the explicit coupling between flow and membrane elasticity introduces a numerical stiffness. The membrane energy in the case of τ = 1.2·10 −5 oscillates for t < 0.0015. These oscillations smooth out when the energy dissipates, leading to the exact same behaviour as for the smallest time step. Consequently, even with large time steps, the membrane energy dissipates in the same manner as with small time steps and the resulting shapes show nearly no differences. However, the coupling of bending stiffness and surface elasticity makes the membrane movement more subtle, such that it takes relatively long to reach the stationary state. With the largest time step, 2000 time steps were required. The total compute time in this case was ≈ 30 minutes. The stretching dominant case has been tested with time steps τ = 4 · 10 −5 , 4 · 10 −6 , 4 · 10 −7 . Images are shown in Fig. 5.3(g-i). As for the bending stiffness dominant case, the membrane energy dissipates in the same manner for small and large time steps and there are only very small shape differences between the shapes for the respective time step sizes. A convergence study is done also for the time step analysis. The results can be seen in Eq. 5.3 and show the expected first order convergence.

Timestep stabilization of the stretching force
In some applications, the time step restrictions due to the surface forces can be rather strict. In general, such restrictions appear when the elastic force is stiff, and the forces are explicitly calculated from the previous membrane position. To overcome this problem, implicit and semiimplicit schemes have been presented to compute the surface tension and bending forces [44]. The idea is to include the interface movement monolithically in the flow solver to evaluate the surface forces implicitly at the newly computed time step. This approach has been shown to work extremely well in the immersed boundary method, and can in principle be combined with our present ALE method.
However, while this methodology has been proposed for surface tension and bending forces, we are not aware of a similar method for the in-plane elastic stretching force. Accordingly, we present a novel approach to monolithically couple the in-plane stretching force to the flow in this section.
The crucial issue in the stretching force (Eq. (14)) is the implicit treatment of the principal stretches λ 1 and λ 2 . A prediction for the values of the principal stretches in the upcoming time step can be obtained from the following evolution equations [9,Appendix]: where ∇ Γ is the (non-axisymmetric) surface divergence in the two-dimensional domain. Accordingly, we can approximate the new values of the principle stretches by where the old solution λ n−1 i is calculated, as before, from the current point coordinates (see Eq. (18)) on Γ.
These equations can now be solved together with the Navier-Stokes equation in one system. To complete the monolithic coupling, the term which includes the implicitly calculated values of the principal stretches. The capability of this implicit strategy has been tested with the parameters of the stretching dominant case. Fig. 5.4 shows the grid points of the membrane for the three different meshes with sizes h 1 , h 2 , and h 3 and different time step sizes. The time step sizes have been chosen to be maximal in each case, i.e. such that the usage of a larger time step would lead to crashing simulations due to oscillations on the membrane points.
As seen in Fig. 5.4 the implicit approach permits an enlargement of the time step size by a factor of 8 − 10. Also notable is that the maximum time step is roughly proportional to the grid size, for both, implicit and explicit approach. Only for the coarsest mesh ( Fig. 12(a)), differences between the explicit and implicit approach are visible, which can be attributed to the large time step sizes. Hence, in particular for higher mesh resolution at the membrane, the implicit approach can lead to high quality results with up to 10 times faster simulations.

Uniaxial compression of a biological cell
As mentioned in the Introduction Sec. 1, the ALE method presented in this work offers the opportunity to restrict the simulation to only one of the fluid phases. This capability can be very useful if the viscosity ratio between the two fluid phases is large such that the phase with the smaller viscosity has little influence on the results and can be neglected. This does not only speed up the simulations but can also greatly simplify the coupling of the elastic surface to exterior forces, as for example in a contact problem.
We illustrate this in the following, as we consider a fluid-filled elastic shell under external compression. The application behind this test case is the uniaxial compression of a biological cell during an atomic force microscopy (AFM) experiment, which is used to measure cell mechanical properties. The viscosity of the intracellular fluid is typically several orders of magnitude larger than the viscosity of the surrounding medium. A detailed description of the problem can be found in [42].
Consider a rounded biological cell confined between two parallel plates. During the experiment, the upper plate moves with a prescribed velocity v compress until the initial plate distance h(t) = h 0 is decreased over time to h(t) = h 0 − ∆h, see Fig. 5.5(a). The necessary force F is constantly measured during compression. Matching experimental data with simulations can be used to extract the surface properties (i.e. K A and K S ) of the cell's elastic shell (the actin cortex). The experimental setup along with the simulation domain is illustrated in Fig. 5.5(a,b). In the simulations, the interior of the cell is denoted by Ω 1 which is bounded by the elastic cell cortex Γ and the symmetry axis. Γ is subdivided into the area touching the plates Γ p and the free surface area Γ f . During compression a part of the free surface will touch the plate, accordingly Γ p and Γ f are time-dependent: A simple contact algorithm is implemented in the numerical simulation: Surface grid points are initially marked to belong to either Γ p or Γ f , as soon as a grid point of Γ f touches the upper plate (x ≥ h(t)), it is shifted to Γ p . The interface curve of Γ for the initial meshes with plate distance h = h 0 is given by a minimal surface calculated according to equations described in [13]. The surrounding medium is neglected to avoid the complicated numerical handling of it being squeezed out of the contact region. Accordingly, the Navier-Stokes equations (3)-(4) are only solved in Ω 1 . The corresponding boundary conditions Eq. (5) have to be adapted where t is the tangential vector to Γ p . Numerical studies have been conducted in the realistic physical parameter regime h 0 ∈  Fig. 13(c) shows exemplary streamlines during compression, as fluid is driven towards the free boundary Γ f , extending the cell's radius. The simulations reproduce the characteristic force response of biological cells, matching simulations with experiments can provide new insights in cell mechanical properties [42].

Shape oscillations of novel microswimming shells
Finally, we illustrate the versatility of our method by simulating hollow elastic shells immersed in a fluid, i.e. a combination of a compressible fluid inside and an incompressible fluid outside an elastic shell. Such shells have been recently proposed in [11] as a very powerful new mechanism for microscopic swimming, since they are extremely fast, simple in shape and controllable from the outside (e.g. by ultrasound driven pressure oscillations).
The idea of the propulsion mechanism is as follows. Consider a microscopic spherical elastic shell, filled with air and immersed in a liquid in the low Reynolds number regime. The shell has a weak spot, i.e. a small area, where its (otherwise uniform) thickness is slightly reduced. The outer fluid pressure is assumed to be controllable, e.g. via ultrasound. If the outer pressure is increased, the air inside the shell will be compressed such that the shell shrinks but remains spherical. When the pressure difference exceeds a shell-specific threshold, the shell will buckle at the weak spot [11]. Decreasing the outer pressure leads to inflation of the shell and debuckling. The motion of the shell during such a buckling cycle induces a fluid flow moving the shell in direction of the weak spot. The experimental setting is illustrated in Fig. 14(a). Images of different stages of a buckling cycle are shown in Fig. 14(c).
Here, we present the first numerical simulations of this process, which can help to gain a better understanding of the buckling dynamics and the influence of parameters on the swimming efficiency. The sudden release of stored elastic energy during the buckling of the shell leads to very high flow rates such that the hydrodynamics are significantly influenced by inertial forces. Hence, even at the microscale, the full Navier-Stokes equations are necessary to describe the process [11]. The simulation of the inner phase Ω 1 is not necessary due to the low density and viscosity of air. Instead, we assume a homogeneous air pressure p 1 inside and use adiabatic gas theory, to relate this pressure to the inner shell volume V by p 1 = p 1,0 (V 0 /V) 1.4 , where p 1,0 and V 0 denote the respective initial values. Accordingly, the stress exerted by the air is reduced to p 1 n whereupon the stress boundary condition (5) becomes The weak spot is introduced by slightly decreasing the elastic surface moduli locally around the membrane point touching the symmetry axis (R = 0) on the right. Note, that the numerical results are independent of the exact amount of this reduction, as the weak spot only serves as a locator for the occurring buckling instability. Numerical parameters are as in [11] where a macroscopic shell diameter of ≈ 5 cm and increased viscosity (glycerol η = 1 Pa or oil η = 37.5 Pa) was used to mimic flow conditions at the microscale. The initial air pressure is p 1,0 = 1 bar.
During the simulation, a sudden jump of the external pressure is imposed at t = 0 by setting p 0 = 1.4 bar at the outer boundaries of the computational domain. This pressure jump leads to sudden shrinkage and buckling of the shell. The concomitant complex flow patterns induce a thrust of the shell, propelling it to the right. Debuckling is imposed at t = 100ms leading to inflation and debuckling and again a propulsion to the right, see Fig. 14 could be discovered in the experiments [11]. Our simulations shall be used in the future in collaboration with the authors of [11] to gain understanding of the complex dynamical coupling of surface shape deformations, pressure oscillations and shell propulsion to create efficient novel microswimmers.

Conclusion
In this paper we have presented a novel ALE method to simulate axisymmetric elastic surfaces immersed in Navier-Stokes fluids. As inherent to ALE methods, the grid is matched to the elastic material which can therefore be resolved with relatively few grid points. The axisymmetric setting reduces the system effectively to a two-dimensional problem. Accordingly, the method combines high accuracy with computational efficiency. This is confirmed in several numerical test cases dominated either by surface tension, bending stiffness or in-plane elasticity. In all these cases we find that numerical errors are relatively small even for relatively coarse grids and converge with order 1-2 with respect to grid size and time step size. The computational times are on the order of minutes on a single core CPU.
Further, it is straightforward to account for completely different physics in both surrounding domains. As an example, we present first simulations of the observed shape oscillations of novel microswimming shells [11] and the uniaxial compression of biological cells filled with cytoplasm [42]. In collaboration with (bio)physicists, the method is currently used to get more insight into the dynamics of microswimming shells and the cellular cortex. As a next step we plan to extend the method to full three dimensions, which requires a new definition of the surface forces. 20 Figure 15: Mesh as it is imported by the simulation.

Implementation details I -mesh and boundary conditions
The presented method is implemented in the finite element toolbox AMDiS [46,47]. While AMDiS does support Taylor-Hood (P2/P1) elements, it does not support different solution variables being defined on different meshes, as necessary for the pressure p i here. Accordingly, we use some technical workarounds to implement the proposed method, described in the following.
A single mesh that contains both domains Ω 0 , Ω 1 as seperate disconnected parts is constructed in gmsh [48]. To easily distinguish between internal and external fluid, the internal part is translated below the symmetry axis (r < 0). Sec. 7.1 shows the mesh in the initial state. At the beginning of the simulation an indicator function is created discriminating between the two mesh parts based on the r−value (grid points with negative r belong to the internal fluid). Immediately afterwards, the internal fluid mesh is shifted upwards to obtain a spatially matched grid, which is yet unconnected at the interface.
The Navier-Stokes equations are assembled on both mesh partitions separately using Lagrangian P2/P1 elements. The implementation of interfacial stress balance Eq. (5) and velocity continuity are technically realized as follows. The discrete membrane Γ h = Γ h,0 ∪ Γ h,1 is composed of the membrane boundaries of the mesh for the external and internal fluid, respectively. Due to the disconnected grid, we can only enforce one-sided stress conditions, which we prescribe as follows −p 0 I + η 0 ∇v + (∇v) T · n = − ∂E ∂Γ on Γ h,0 − −p 1 I + η 1 ∇v + (∇v) T · n = 0 on Γ h,1 Define n as the number of (P2−)DOFs on the membrane. Let j = 0, . . . , n − 1, be an ordered (e.g. counterclockwise) numbering of the DOFs on Γ h,0 and Γ h,1 , such that the numbers of both boundaries are equal. At the coordinate of each membrane DOF j there exist exactly two P2-test-functions ψ j 0 , ψ j 1 which are non-zero, with supp(ψ j 0 ) ⊂ Ω 0 and supp(ψ j 1 ) ⊂ Ω 1 . These test-functions are each used twice, once for testing with the v x -equation and once for testing with the v r -equation. Each of these test cases corresponds to a single row in the assembled Finite-Element matrix. The correct stress jump condition can be realized by adding the two rows related to ψ j 1 to the two respective rows related to ψ j 0 , for every j = 0, . . . , n − 1. This adds both one-sided test-functions and effectively mimics that the momentum equations were tested with a single test-function of an interfacially connected grid. Also the boundary conditions Eq. added up to recover the correct stress jump condition Eq. (5). As an illustration, Eq. 7.1 shows the test-functions in the case of an 1D mesh.
The obsolete rows related to ψ j 1 can then be overwritten to enforce the continuity of the velocities at the interface v j 1 − v j 0 = 0, on Γ h,1 .
As v = (v x , v r ) the above describes two continuity conditions which can be assembled in the v xand v r -rows tested with ψ j 1 . The continuity of velocity ensures in particular, that corresponding DOFs on Γ h,0 and Γ h,1 share the same coordinate points for all times.

Implementation details II -extension of surface equations to Ω
In AMDiS, it is not possible to solve 2D bulk equations and 1D surface equations all in one system. To compute the principal stretches implicitly, it is therefore necessary to extend Eq. (32) to Ω. The equations for the principal stretches then read where v Γ,ext is the extension of the interfacial velocity to Ω constant in normal direction. R ext is the distance to the symmetry axis for any point in Ω. The equations for the principal stretches are solved using P1 elements. The calculation of the normal extension is based on a Hopf-Lax algorithm described in [49]. The surface divergence then reads ∇ Γ · v Γ,ext = P : ∇v Γ,ext , with the projection onto the tangent space P = I − nn T .