2.1. Different surface parameterisations

The surface under consideration is denoted $\mathcal{S}$ and its surface points are described by the parametrisation

(2.1) \begin{equation} \begin{array}{l} {\boldsymbol{x}} = {\boldsymbol{x}}(\zeta ^\alpha,t). \end{array} \end{equation}

Here, $\zeta ^\alpha$ denotes a set of arbitrary surface coordinates ( $\alpha = 1,\,2$ ) that is associated with the ALE frame of reference. In a computational description, this frame of reference can be taken as the computational grid or mesh. It can also be associated with an observer motion (Nitschke & Voigt Reference Nitschke and Voigt2022).

Coordinate $\zeta ^\alpha$ is generally different from the convective coordinate $\xi ^\alpha$ , which is used to track material points and which defines the material time derivative

(2.2) \begin{equation} \begin{array}{l} \dot {(\ldots)} = \displaystyle \left.\frac {\partial {\ldots }}{\partial {t}}\right |_{\xi ^\alpha }. \end{array} \end{equation}

It is also different from the surface-fixed coordinate $\theta ^\alpha$ , whose material time derivative defines the tangential fluid velocity, i.e. $v^\alpha = \dot \theta ^\alpha$ . In other words

(2.3) \begin{equation} \begin{array}{l} {{\boldsymbol{x}}}=\hat {{\boldsymbol{x}}}(\xi ^\alpha,t), \end{array} \end{equation}

is a Lagrangian surface description (that follows the material motion), while

(2.4) \begin{equation} \begin{array}{l} {\boldsymbol{x}} = \tilde {\boldsymbol{x}}(\theta ^\alpha,t), \end{array} \end{equation}

is an in-plane Eulerian surface description. (This notation slightly differs from that of Sahu et al. (Reference Sahu, Omar, Sauer and Mandadapu2020), which places no tilde on $\boldsymbol{x}$ in (2.4), but places a check on $\boldsymbol{x}$ in (2.1). This then implies a check will appear on indices that can be written without accent here). Equation (2.1) on the other hand is an arbitrary Lagrangian–Eulerian surface description that contains the two special cases $\zeta ^\alpha =\xi ^\alpha$ and $\zeta ^\alpha =\theta ^\alpha$ . The mappings (2.1) and (2.3) are illustrated in figure 1. They induce a functional relationship between $\zeta ^\alpha$ and $\xi ^\alpha$ , i.e. $\zeta ^\alpha = \zeta ^\alpha (\xi ^\beta,t)$ and $\xi ^\alpha = \xi ^\alpha (\zeta ^\beta,t)$ .

Figure 1.

The ALE surface parameterisation: the surface $\mathcal{S}$ can be described by the ALE mapping ${\boldsymbol{x}}={\boldsymbol{x}}(\zeta ^\alpha,t)$ from the ALE parameter domain $\mathcal{P}$ , and the material mapping ${\boldsymbol{x}}=\hat {\boldsymbol{x}}(\xi ^\alpha,t)$ from the material parameter domain $\hat {\mathcal{P}}$ . From these mappings follow the tangent vectors ${\boldsymbol{a}}_\alpha:= {\boldsymbol{x}}_{\!,\alpha }$ and ${\boldsymbol{a}}_{\hat \alpha }:= \hat {\boldsymbol{x}}_{\!,\hat \alpha }$ .

Parameterisations (2.1), (2.3) and (2.4) define the three parametric derivatives

(2.5) \begin{equation} \begin{array}{l} \ldots _{,\alpha }:= \displaystyle \frac {\partial {\ldots }}{\partial {\zeta ^\alpha }},\quad \ldots _{,\hat \alpha }:= \displaystyle \frac {\partial {\ldots }}{\partial {\xi ^\alpha }},\quad \ldots _{,\tilde \alpha }:= \displaystyle \frac {\partial {\ldots }}{\partial {\theta ^\alpha }}. \end{array} \end{equation}

Applied to $\boldsymbol{x}$ , this defines the tangent vectors

(2.6a–c) \begin{equation} \begin{array}{l} {\boldsymbol{a}}_\alpha:= {\boldsymbol{x}}_{\!,\alpha },\quad {\boldsymbol{a}}_{\hat \alpha }:= \hat {\boldsymbol{x}}_{\!,\hat \alpha },\quad {\boldsymbol{a}}_{\tilde \alpha }:= \tilde {\boldsymbol{x}}_{\!,\tilde \alpha }. \end{array} \end{equation}

Together with the normal vector

(2.7) \begin{equation} \begin{array}{l} {\boldsymbol{n}}:= \displaystyle \frac {{\boldsymbol{a}}_1\times {\boldsymbol{a}}_2}{\|{\boldsymbol{a}}_1\times {\boldsymbol{a}}_2\|} = \displaystyle \frac {{\boldsymbol{a}}_{\hat 1}\times {\boldsymbol{a}}_{\hat 2}}{\|{\boldsymbol{a}}_{\hat 1}\times {\boldsymbol{a}}_{\hat 2}\|} = \displaystyle \frac {{\boldsymbol{a}}_{\tilde 1}\times {\boldsymbol{a}}_{\tilde 2}}{\|{\boldsymbol{a}}_{\tilde 1}\times {\boldsymbol{a}}_{\tilde 2}\|}, \end{array} \end{equation}

they can be used as a basis to decompose general vectors into tangential and normal components, i.e.

(2.8) \begin{equation} \begin{array}{l} {\boldsymbol{v}} = v^\alpha \,{\boldsymbol{a}}_{\alpha } + v\,{\boldsymbol{n}} = v^{\hat \alpha }\,{\boldsymbol{a}}_{\hat \alpha } + v\,{\boldsymbol{n}} = v^{\tilde \alpha }\,{\boldsymbol{a}}_{\tilde \alpha } + v\,{\boldsymbol{n}}. \end{array} \end{equation}

Tangent vectors ${\boldsymbol{a}}_{\alpha }$ , ${\boldsymbol{a}}_{\hat \alpha }$ and ${\boldsymbol{a}}_{\tilde \alpha }$ define the surface metrics

(2.9a–c) \begin{equation} \begin{array}{l} a_{\alpha \gamma }:= {\boldsymbol{a}}_{\alpha }\boldsymbol{\cdot} {\boldsymbol{a}}_{\gamma },\quad a_{\hat \alpha \hat \gamma }:= {\boldsymbol{a}}_{\hat \alpha }\boldsymbol{\cdot} {\boldsymbol{a}}_{\hat \gamma },\quad a_{\tilde \alpha \tilde \gamma }:= {\boldsymbol{a}}_{\tilde \alpha }\boldsymbol{\cdot} {\boldsymbol{a}}_{\tilde \gamma }, \end{array} \end{equation}

and the curvature tensor components

(2.10a–c) \begin{equation} \begin{array}{l} b_{\alpha \gamma }:= {\boldsymbol{a}}_{\alpha,\gamma }\boldsymbol{\cdot} {\boldsymbol{n}},\quad b_{\hat \alpha \hat \gamma }:= {\boldsymbol{a}}_{\hat \alpha,\hat \gamma }\boldsymbol{\cdot} {\boldsymbol{n}},\quad b_{\tilde \alpha \tilde \gamma }:= {\boldsymbol{a}}_{\tilde \alpha,\tilde \gamma }\boldsymbol{\cdot} {\boldsymbol{n}}. \end{array} \end{equation}

Through the inverse metrics

(2.11a–c) \begin{equation} \begin{array}{l} [a^{\alpha \gamma }]:= [a_{\alpha \gamma }]^{-1} ,\quad [a^{\hat \alpha \hat \gamma }]:= [a_{\hat \alpha \hat \gamma }]^{-1} ,\quad [a^{\tilde \alpha \tilde \gamma }]:= [a_{\tilde \alpha \tilde \gamma }]^{-1}, \end{array} \end{equation}

the dual tangent vectors

(2.12a–c) \begin{equation} \begin{array}{l} {\boldsymbol{a}}^\alpha:= a^{\alpha \gamma }{\boldsymbol{a}}_\gamma ,\quad {\boldsymbol{a}}^{\hat \alpha }:= a^{\hat \alpha \hat \gamma }{\boldsymbol{a}}_{\hat \gamma },\quad {\boldsymbol{a}}^{\tilde \alpha }:= a^{\tilde \alpha \tilde \gamma }{\boldsymbol{a}}_{\tilde \gamma }, \end{array} \end{equation}

are defined. (Following index notation, summation is implied on repeated indices within terms). They in turn define the Christoffel symbols

(2.13a–c) \begin{equation} \begin{array}{l} \Gamma ^\mu _{\alpha \gamma }:= {\boldsymbol{a}}^\mu \boldsymbol{\cdot} {\boldsymbol{a}}_{\alpha,\gamma },\quad \Gamma ^{\hat \mu }_{\hat \alpha \hat \gamma }:= {\boldsymbol{a}}^{\hat \mu }\boldsymbol{\cdot} {\boldsymbol{a}}_{\hat \alpha,\hat \gamma },\quad \Gamma ^{\tilde \mu }_{\tilde \alpha \tilde \gamma }:= {\boldsymbol{a}}^{\tilde \mu }\boldsymbol{\cdot} {\boldsymbol{a}}_{\tilde \alpha,\tilde \gamma }. \end{array} \end{equation}

Based on the preceding quantities, further quantities in surface differential geometry, e.g. see Sauer (Reference Sauer2018), can be formulated in each basis. In particular, the surface gradient and surface divergence of a general scalar $\phi$ and vector $\boldsymbol{v}$ become

(2.14) \begin{align} \begin{array}{l} \boldsymbol{\nabla}_{\!s}\phi = \phi _{,\alpha }\,{\boldsymbol{a}}^\alpha = \phi _{,\hat \alpha }\,{\boldsymbol{a}}^{\hat \alpha } = \phi _{,\tilde \alpha }\,{\boldsymbol{a}}^{\tilde \alpha }, \end{array}\end{align}

(2.15) \begin{align} \begin{array}{l} \boldsymbol{\nabla}_{\!s}{\boldsymbol{v}} = {\boldsymbol{v}}_{\!,\alpha }\otimes {\boldsymbol{a}}^\alpha = {\boldsymbol{v}}_{\!,\hat \alpha }\otimes {\boldsymbol{a}}^{\hat \alpha } = {\boldsymbol{v}}_{\!,\tilde \alpha }\otimes {\boldsymbol{a}}^{\tilde \alpha } \end{array}\end{align}

and

(2.16) \begin{equation} \begin{array}{l} \textrm{div}_{s} \,{\boldsymbol{v}} = {\boldsymbol{v}}_{\!,\alpha }\boldsymbol{\cdot} {\boldsymbol{a}}^\alpha = {\boldsymbol{v}}_{\!,\hat \alpha }\boldsymbol{\cdot} {\boldsymbol{a}}^{\hat \alpha } = {\boldsymbol{v}}_{\!,\tilde \alpha }\boldsymbol{\cdot} {\boldsymbol{a}}^{\tilde \alpha }. \end{array} \end{equation}

The comma here denotes the parametric derivative according to (2.5). In case of $\boldsymbol{v}$ , it applies to the full vector and not just its components.

The description following in § 3 exclusively uses basis $\{{\boldsymbol{a}}_1,\,{\boldsymbol{a}}_2,\,{\boldsymbol{n}}\}$ induced by (2.1). Basis $\{\hat {\boldsymbol{a}}_1,\,\hat {\boldsymbol{a}}_2,\,{\boldsymbol{n}}\}$ induced by (2.3) and basis $\{\tilde {\boldsymbol{a}}_1,\,\tilde {\boldsymbol{a}}_2,\,{\boldsymbol{n}}\}$ induced by (2.4) are not used further, apart from some derivations. But it is important to note that they exist and provide alternative descriptions. Description (2.1) is the most general one, that contains the other two as special cases. Appendix A contains the transformation rules that can be used to adapt an expression from one parameterisation to another. They can be used to show that tensor invariants, such as the curvature tensor invariants $H:= a^{\alpha \beta }b_{\alpha \beta }/2$ and $\kappa:= \det [b_{\alpha \beta }]/\det [a_{\alpha \beta }]$ are invariant with respect to the choice of basis. Essentially, all quantities without free index are invariant with respect to the parametrisation, which includes its change – a property often denoted reparametrisation invariance, e.g. see Guven & Vázquez-Montejo (Reference Guven and Vázquez-Montejo2018).

An important aspect to note, is that the material time derivative (2.2) only commutes with the parametric derivative with respect to $\xi ^\alpha$ , but not with the other two parametric derivatives in (2.5). This is needed in the following.

2.2. The ALE expressions for various material time derivatives

The material velocity is defined as the material time derivative of $\boldsymbol{x}$ , i.e.

(2.17) \begin{equation} \begin{array}{l} {\boldsymbol{v}} = \dot {\boldsymbol{x}} = \displaystyle \left.\frac {\partial {\hat {\boldsymbol{x}}}}{\partial {t}}\right |_{\xi ^\alpha }. \end{array} \end{equation}

Changing to $\zeta ^\alpha$ , this expands into

(2.18) \begin{equation} \begin{array}{l} {\boldsymbol{v}} = \displaystyle \left.\frac {\partial {{\boldsymbol{x}}}}{\partial {t}}\right |_{\zeta ^\alpha } + \frac {\partial {{\boldsymbol{x}}}}{\partial {\zeta ^\alpha }} \left.\frac {\partial {\zeta ^\alpha }}{\partial {t}}\right |_{\xi ^\beta }. \end{array} \end{equation}

Introducing the velocity of the ALE frame,

(2.19) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{{m}}:= \displaystyle \left.\frac {\partial {{\boldsymbol{x}}}}{\partial {t}}\right |_{\zeta ^\alpha }, \end{array} \end{equation}

which in computational methods is also referred to as the mesh velocity, and using (2.6a ), (2.18) can be rewritten into

(2.20) \begin{equation} \begin{array}{l} {\boldsymbol{v}} = {\boldsymbol{v}}_{{m}} + \dot \zeta ^\alpha \, {\boldsymbol{a}}_\alpha . \end{array} \end{equation}

This expression admits the two special cases:

1. (i) In-plane Lagrangian description: $\zeta ^\alpha = \xi ^\alpha$ , for which $\dot \zeta ^\alpha =0$ and thus ${\boldsymbol{v}}_{{m}}={\boldsymbol{v}}$ ; and

2. (ii) In-plane Eulerian description: $\zeta ^\alpha = \theta ^\alpha$ , for which $\dot \zeta ^\alpha =v^\alpha$ and thus ${\boldsymbol{v}}_{{m}}=v{\boldsymbol{n}}$ due to (2.8).

Due to these properties, $\xi ^\alpha$ is also referred to as the convected coordinate, while $\theta ^\alpha$ is referred to as the surface-fixed coordinate.

In general, $\dot \zeta ^\alpha$ characterises the tangential velocity difference

(2.21) \begin{equation} \begin{array}{l} \dot \zeta ^\alpha = {\boldsymbol{a}}^\alpha \boldsymbol{\cdot} ({\boldsymbol{v}}-{\boldsymbol{v}}_{{m}}), \end{array} \end{equation}

according to (2.20). This leads to the following interpretation of (2.20): the (absolute) material velocity $\boldsymbol{v}$ is composed of the (absolute) mesh velocity ${\boldsymbol{v}}_{{m}}$ plus the relative material velocity with respect to the mesh, $\dot \zeta ^\alpha \, {\boldsymbol{a}}_\alpha$ . A graphical representation of this is shown in the example of § 6.1, see figure 2.

Likewise to (2.17), the material acceleration is defined as

(2.22) \begin{equation} \begin{array}{l} \dot {{\boldsymbol{v}}} = \displaystyle \left.\frac {\partial {{\boldsymbol{v}}}}{\partial {t}}\right |_{\xi ^\alpha }. \end{array} \end{equation}

Changing to $\zeta ^\alpha$ , this expands into

(2.23) \begin{equation} \begin{array}{l} \dot {{\boldsymbol{v}}} = \displaystyle \left.\frac {\partial {{\boldsymbol{v}}}}{\partial {t}}\right |_{\zeta ^\alpha } + \frac {\partial {{\boldsymbol{v}}}}{\partial {\zeta ^\alpha }} \left.\frac {\partial {\zeta ^\alpha }}{\partial {t}}\right |_{\xi ^\beta }. \end{array} \end{equation}

Introducing

(2.24) \begin{equation} \begin{array}{l} (\ldots)^{\prime}:= \displaystyle \left.\frac {\partial {\ldots }}{\partial {t}}\right |_{\zeta ^\alpha }, \end{array} \end{equation}

and using (2.21) and (2.15), (2.23) can be rewritten into

(2.25) \begin{equation} \begin{array}{l} \dot {{\boldsymbol{v}}} = {\boldsymbol{v}}^{\prime} + \boldsymbol{\nabla}_{\!s}{\boldsymbol{v}}\,\left ({\boldsymbol{v}}-{\boldsymbol{v}}_{{m}}\right)\kern-2.5pt, \end{array} \end{equation}

which is the analogous surface version of the well-known three-dimensional ALE equation (Donea & Huerta Reference Donea and Huerta2003). Using (2.15), (2.25) can also be written as

(2.26) \begin{equation} \begin{array}{l} \dot {{\boldsymbol{v}}} = {\boldsymbol{v}}^{\prime} + {\boldsymbol{v}}_{\!,\alpha }\,\left (v^\alpha -v^\alpha _{m}\right)\kern-2.5pt, \end{array} \end{equation}

where $v^\alpha:= {\boldsymbol{a}}^\alpha \boldsymbol{\cdot} {\boldsymbol{v}}$ and $v^\alpha_{m}:= {\boldsymbol{a}}^\alpha \boldsymbol{\cdot} {\boldsymbol{v}}_{{m}}$ . It is emphasised that the temporal and spatial differentiation in ${\boldsymbol{v}}_{\!,\alpha }$ is generally not exchangeable, i.e. ${\boldsymbol{v}}_{\!,\alpha } \neq \dot {\boldsymbol{a}}_\alpha$ . (The identity ${\boldsymbol{v}}_{,\alpha } = \dot {\boldsymbol{a}}_\alpha$ , used in Rangamani et al. (Reference Rangamani, Agrawal, Mandadapu, Oster and Steigmann2013) and Sahu et al. (Reference Sahu, Sauer and Mandadapu2017, Reference Sahu, Omar, Sauer and Mandadapu2020) for the Eulerian frame, is not general as it relies on the special case $\dot \zeta ^\gamma _{,\alpha }=0$ . This is only the case for particular ALE descriptions such as the Lagrangian description $\zeta ^\gamma = \xi ^\gamma$ , or when $\dot \zeta ^\gamma _{,\alpha }$ is negligibly small). Instead

(2.27) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{\!,\alpha } = \dot {\boldsymbol{a}}_\alpha + \dot \zeta ^\gamma _{,\alpha }\,{\boldsymbol{a}}_\gamma , \end{array} \end{equation}

see Appendix B. Here,

(2.28) \begin{equation} \begin{array}{l} \dot \zeta ^\gamma _{,\alpha }:= \displaystyle \frac {\partial {\dot \zeta ^\gamma }}{\partial {\zeta ^\alpha }} = \displaystyle \frac {\partial {}}{\partial {\zeta ^\alpha }}\left (\frac {\partial {\zeta ^\gamma }}{\partial {t}}\right)_{\!\!\xi ^\beta }. \end{array} \end{equation}

Also here the order of differentiation cannot be exchanged (i.e. $\dot \zeta ^\gamma _{,\alpha }$ is generally not equal to $(\zeta ^\gamma _{,\alpha }\dot)=0$ ). Equation (2.27) admits the two special cases: 1. $\zeta ^\alpha =\xi ^\alpha$ , for which $\dot \zeta ^\gamma _{,\alpha }=0$ , and 2. $\zeta ^\alpha =\theta ^\alpha$ , for which $\dot \zeta ^\gamma _{,\alpha } = v^\gamma _{,\alpha }$ . The former case implies that

(2.29) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{\!,\hat \alpha } = \dot {\boldsymbol{a}}_{\hat \alpha } \end{array} \end{equation}

The expansion used in (2.18) and (2.23) is a fundamental relation. It generalises to

(2.30) \begin{equation} \begin{array}{l} \boxed{\dot {(\ldots)} = (\ldots)^{\prime} + \ldots _{,\alpha}{\dot \zeta}^{\alpha}}, \end{array} \end{equation}

and is denoted the fundamental surface ALE equation in the following. As an example, the material time derivative of a scalar field $\phi$ , such as the density $\rho$ , then follows as

(2.31) \begin{equation} \begin{array}{l} \dot \phi:= \displaystyle \left.\frac {\partial {\phi }}{\partial {t}}\right |_{\xi ^\alpha } = \phi ^{\prime} + \phi _{,\alpha }\,\dot \zeta ^\alpha , \end{array} \end{equation}

which, in view of (2.21), leads to

(2.32) \begin{equation} \begin{array}{l} \dot \phi = \phi ^{\prime} + \phi _{,\alpha }\,\left (v^\alpha -v_{m}^\alpha \right)\kern-2.5pt, \end{array} \end{equation}

or equivalently

(2.33) \begin{equation} \begin{array}{l} \dot \phi = \phi ^{\prime} + \boldsymbol{\nabla}_{\!s}\phi \boldsymbol{\cdot} \left ({\boldsymbol{v}}-{\boldsymbol{v}}_{{m}}\right)\kern-2.5pt. \end{array} \end{equation}

2.3. Velocity gradient

An important object for characterising surface flows is the symmetric surface velocity gradient (also known as the rate-of-deformation tensor) $\boldsymbol{d}:= (\boldsymbol{\nabla}_{\!s}{\boldsymbol{v}} + \boldsymbol{\nabla}_{\!s}{\boldsymbol{v}}^{\textrm{T}})/2$ , which in view of (2.15) and (2.27) becomes

(2.34) \begin{equation} \begin{array}{l} 2\boldsymbol{d} = {\boldsymbol{v}}_{\!,\alpha } \otimes {\boldsymbol{a}}^\alpha + {\boldsymbol{a}}^\alpha \otimes {\boldsymbol{v}}_{\!,\alpha }\end{array} \end{equation}

and

(2.35) \begin{equation} \begin{array}{l} 2\boldsymbol{d} = \dot {\boldsymbol{a}}_\alpha \otimes {\boldsymbol{a}}^\alpha + {\boldsymbol{a}}^\alpha \otimes \dot {\boldsymbol{a}}_\alpha + \dot \zeta ^\gamma _{,\alpha }\left ( {\boldsymbol{a}}_\gamma \otimes {\boldsymbol{a}}^\alpha + {\boldsymbol{a}}^\alpha \otimes {\boldsymbol{a}}_\gamma \right)\kern-2.5pt. \end{array} \end{equation}

From this one finds the in-plane components

(2.36) \begin{equation} \begin{array}{l} 2d_{\alpha \beta } = \dot a_{\alpha \beta } + \dot \zeta ^\gamma _{,\alpha }\,a_{\gamma \beta } + a_{\alpha \gamma }\,\dot \zeta ^\gamma _{,\beta }, \end{array} \end{equation}

and

(2.37) \begin{equation} \begin{array}{l} 2d^{\alpha \beta } = -\dot a^{\alpha \beta } + \dot \zeta ^\alpha _{,\gamma }\,a^{\gamma \beta } + a^{\alpha \gamma }\,\dot \zeta ^\beta _{,\gamma }, \end{array} \end{equation}

that are connected through $\dot a^{\alpha \beta } = -a^{\alpha \gamma } \dot a_{\gamma \delta }\, a^{\delta \beta }$ . In general, $\boldsymbol{d}$ also has out-of-plane components, but those are not relevant to the constitutive models considered in § 4. The in-plane components can be arranged in the in-plane tensor

(2.38) \begin{equation} \begin{array}{l} \boldsymbol{d}_{{s}}:= d_{\alpha \beta }\,{\boldsymbol{a}}^\alpha \otimes {\boldsymbol{a}}^\beta = d^{\alpha \beta }{\boldsymbol{a}}_\alpha \otimes {\boldsymbol{a}}_\beta . \end{array} \end{equation}

2.4. Vorticity

Another important object for characterising flows is the vorticity. It derives from the curl of the velocity. To characterise surface flows we thus introduce the surface curl of the velocity

(2.39) \begin{equation} \begin{array}{l} \textrm{curl}_{{s}}{\boldsymbol{v}}:= {\boldsymbol{a}}^\alpha \times {\boldsymbol{v}}_{\!,\alpha }, \end{array} \end{equation}

analogous to the definitions of surface gradient and surface divergence in (2.15) and (2.16). The surface vorticity $\omega$ is then defined as the normal component of $\textrm{curl}_{{s}}{\boldsymbol{v}}$ , i.e.

(2.40) \begin{equation} \begin{array}{l} \omega:= \textrm{curl}_{{n}}{\boldsymbol{v}}:= ({\boldsymbol{a}}^\alpha \times {\boldsymbol{v}}_{\!,\alpha })\boldsymbol{\cdot} {\boldsymbol{n}}. \end{array} \end{equation}

It can be shown that

(2.41) \begin{equation} \begin{array}{l} \textrm{curl}_{{n}}{\boldsymbol{v}} = \textrm{div}_{s} \,({\boldsymbol{v}}\times {\boldsymbol{n}}), \end{array} \end{equation}

and

(2.42) \begin{equation} \begin{array}{l} \textrm{curl}_{{n}}{\boldsymbol{v}} = ({\boldsymbol{n}}\times {\boldsymbol{a}}^\alpha)\boldsymbol{\cdot} {\boldsymbol{v}}_{\!,\alpha }. \end{array} \end{equation}

The latter identity motivates the definition of the surface curl of a scalar $\phi$ by

(2.43) \begin{equation} \begin{array}{l} \textrm{curl}_{{s}}\phi:= ({\boldsymbol{n}}\times {\boldsymbol{a}}^\alpha)\,\phi _{,\alpha } = {\boldsymbol{n}}\times \boldsymbol{\nabla}_{\!s}\phi , \end{array} \end{equation}

which is a vector like $\textrm{curl}_{{s}}{\boldsymbol{v}}$ . The surface curl of $\phi$ can also be written as

(2.44) \begin{equation} \begin{array}{l} \textrm{curl}_{{s}}\phi = \phi _{,\alpha }\,\epsilon ^{\alpha \beta }\,{\boldsymbol{a}}_\beta , \end{array} \end{equation}

where $[\epsilon ^{\alpha \beta }] = [0\,1;-1\,0]/\sqrt {\det [a_{\gamma \delta }]}$ is the scaled permutation (or Levi-Civita) symbol.

2.5. Surface stretch

The surface stretch $J$ measures the local area change with respect to the initial configuration of the surface, denoted $\mathcal{S}_0$ . The current surface point ${\boldsymbol{x}}\in \mathcal{S}$ has the initial location $\boldsymbol{X}={\boldsymbol{x}}|_{t=0}\in \mathcal{S}_0$ . Based on the three parameterisations (2.1), (2.3) and (2.4), this leads to the basis vectors $\boldsymbol{A}_\alpha ={\boldsymbol{a}}_\alpha |_{t=0}$ , $\boldsymbol{A}_{\hat \alpha }={\boldsymbol{a}}_{\hat \alpha }|_{t=0}$ and $\boldsymbol{A}_{\tilde \alpha }={\boldsymbol{a}}_{\tilde \alpha }|_{t=0}$ , and the corresponding surface metrics $A_{\alpha \gamma }=a_{\alpha \gamma }|_{t=0}$ , $A_{\hat \alpha \hat \gamma }=a_{\hat \alpha \hat \gamma }|_{t=0}$ and $A_{\tilde \alpha \tilde \gamma }= a_{\tilde \alpha \tilde \gamma }|_{t=0}$ . The surface stretch is given by

(2.45) \begin{equation} \begin{array}{l} J = \sqrt {\det [a_{\hat \alpha \hat \gamma }]}\left /\sqrt {\det [A_{\hat \alpha \hat \gamma }]}\right.. \end{array} \end{equation}

It is obtained in the Lagrangian frame, since only this frame is tracking the physical material motion. The quantity

(2.46) \begin{equation} \begin{array}{l} J_{{m}}:= \sqrt {\det [a_{\alpha \gamma }]}\left /\sqrt {\det [A_{\alpha \gamma }]}\right., \end{array} \end{equation}

on the other hand, tracks (non-physical) area-changes of the ALE frame. Generally $J_{{m}}\,{\neq}\,J$ .

As was already mentioned, the following presentation exclusively uses basis $\{{\boldsymbol{a}}_1,\,{\boldsymbol{a}}_2,\,{\boldsymbol{n}}\}$ . Expression (2.45) therefore becomes impractical. Instead, $J$ can be determined from the surface velocity $\boldsymbol{v}$ through the evolution law

(2.47) \begin{equation} \begin{array}{l} \displaystyle \frac {\dot J}{J} = \textrm{div}_{s} \,{\boldsymbol{v}}, \end{array} \end{equation}

following from (2.16) and $\dot {\boldsymbol{a}}_{\hat \alpha }\boldsymbol{\cdot} {\boldsymbol{a}}^{\hat \alpha } = \dot J/J$ , e.g. see Sauer (Reference Sauer2018). Since $ \textrm{div}_{s} \,{\boldsymbol{v}} = \textrm{tr}\,\boldsymbol{d}$ , one can also write

(2.48) \begin{equation} \begin{array}{l} \displaystyle \frac {\dot J}{J} = d^{\alpha \beta } a_{\alpha \beta }. \end{array} \end{equation}

Using (2.16), (2.27) and $\dot {\boldsymbol{a}}_\alpha \boldsymbol{\cdot} {\boldsymbol{a}}^\alpha = \dot J_{{m}}/J_{{m}}$ , one can further write

(2.49) \begin{equation} \begin{array}{l} \textrm{div}_{s} \,{\boldsymbol{v}} = \displaystyle \frac {\dot J_{{m}}}{J_{{m}}} + \dot \zeta ^\alpha _{,\alpha }. \end{array} \end{equation}

For area-incompressible flows, where $\dot J=0$ , (2.47) leads to the velocity constraint

(2.50) \begin{equation} \begin{array}{l} \textrm{div}_{s} \,{\boldsymbol{v}} = 0. \end{array} \end{equation}

Associated with this constraint is an unknown surface tension, the Lagrange multiplier $q$ . It can be different to the physical surface tension $\gamma$ , as is seen in § 4.2.

3.1. Mass balance

The surface fluid density $\rho =\rho (\zeta ^\alpha,t)$ is governed by the the continuity equation

(3.1) \begin{equation} \begin{array}{l} \dot \rho + \rho \,\textrm{div}_{s} \,{\boldsymbol{v}} = 0, \end{array} \end{equation}

which follows from surface mass balance (Marsden & Hughes Reference Marsden and Hughes1994; Sahu, Sauer & Mandadapu Reference Sahu, Sauer and Mandadapu2017). (Both references write the surface divergence of ${\boldsymbol{v}} = v^\alpha {\boldsymbol{a}}_\alpha + v\,{\boldsymbol{n}}$ in the form $\textrm{div}_{s} \,{\boldsymbol{v}} = v^\alpha _{;\alpha } - 2Hv$ ). In the Lagrangian frame this is a first-order ordinary differential equation (ODE) that only requires the initial condition $\rho (\zeta ^\alpha,0) = \rho _0(\zeta ^\alpha)$ , where $\rho _0$ is the given initial density distribution. In the ALE frame, where $\dot \rho$ can be expanded according to (2.32), this is a first-order PDE that additionally requires a boundary condition to capture the mass influx

(3.2) \begin{equation} \begin{array}{l} j^\alpha:= \rho \,\left ({\boldsymbol{v}}_{{m}}-{\boldsymbol{v}}\right)\boldsymbol{\cdot} {\boldsymbol{a}}^\alpha, \end{array} \end{equation}

on all boundaries where ${\boldsymbol{v}}_{{m}}\neq {\boldsymbol{v}}$ .

Remark 3.1. Inserting (2.47) into (3.1) shows that ODE (3.1) is solved by

(3.3) \begin{equation} \begin{array}{l} \rho = \rho _0/J. \end{array} \end{equation}

This is convenient in a Lagrangian description, where $J$ is easily obtained from (2.45). In an ALE description, (2.45) requires reconstructing basis ${\boldsymbol{a}}_{\hat \alpha }$ . To avoid this, one can either solve (2.47) for $J$ and then get $\rho$ from (3.3), or – equivalently – solve (3.1) for $\rho$ and then get $J$ from (3.3).

3.2. Momentum balance

The surface fluid velocity ${\boldsymbol{v}}={\boldsymbol{v}}(\zeta ^\alpha,t)$ is governed by the equation of motion

(3.4) \begin{equation} \begin{array}{l} \rho \,\dot {{\boldsymbol{v}}} = \boldsymbol{T}^\alpha _{;\alpha } + \boldsymbol{f}, \end{array} \end{equation}

which follows from linear surface momentum balance, see e.g. Sauer & Duong (Reference Sauer and Duong2017). Equation (3.4) is a very general expression, that includes surface flow on fixed and evolving surfaces, as well as flows without and with bending resistance (e.g. Wilmore flows). Here, $\boldsymbol{f} = p\,{\boldsymbol{n}} + f_\alpha \,{\boldsymbol{a}}^\alpha$ is a surface load, while

(3.5) \begin{equation} \begin{array}{l} \boldsymbol{T}^\alpha = {\boldsymbol \sigma}^{\textrm{T}} {\boldsymbol{a}}^\alpha, \end{array} \end{equation}

is the traction vector on a cut through $\mathcal{S}$ that is perpendicular to ${\boldsymbol{a}}^\alpha$ . It depends on the Cauchy stress

(3.6) \begin{align} \begin{array}{l} {\boldsymbol{\sigma }} = N^{\alpha \beta } {\boldsymbol{a}}_\alpha \otimes {\boldsymbol{a}}_\beta + S^\alpha \,{\boldsymbol{a}}_\alpha \otimes {\boldsymbol{n}}, \end{array} \end{align}

that has the components

(3.7a) \begin{align} N^{\alpha \beta} & = \sigma ^{\alpha \beta } + b^\beta _\gamma \,M^{\gamma \alpha},\end{align}

(3.7b) \begin{align} S^\alpha &=- M^{\beta \alpha }_{\,;\beta},\end{align}

for Kirchhoff–Love shells (Naghdi Reference Naghdi1972; Steigmann Reference Steigmann1999). Here, the in-plane membrane stresses $\sigma ^{\alpha \beta }$ and the bending stress couples $M^{\alpha \beta }$ are defined by constitution, see § 4. The out-of-plane shear stress $S^\alpha$ on the other hand, follows from $M^{\alpha \beta }$ , which is a consequence of assuming Kirchhoff–Love kinematics (i.e. neglecting out-of-plane shear deformations). The stress $\sigma ^{\alpha \beta }$ is also referred to as the effective stress (Simo & Fox Reference Simo and Fox1989). On the other hand, $N^{\alpha \beta }$ corresponds to the physical stress according to Cauchy’s theorem (Sahu et al. Reference Sahu, Sauer and Mandadapu2017).

Using (2.25) or (2.26), PDE (3.4) can be expressed in the ALE frame. In order to solve PDE (3.4), an initial condition and boundary conditions on ${\boldsymbol{v}}={\boldsymbol{v}}(\zeta ^\alpha,t)$ are needed. In PDE (3.4), one can also replace $\boldsymbol{T}^\alpha _{;\alpha }$ by

(3.8) \begin{equation} \begin{array}{l} \boldsymbol{T}^\alpha _{;\alpha } = {\boldsymbol \sigma }_{\!;\alpha }^{\textrm{T}}\,{\boldsymbol{a}}^\alpha =: \textrm{div}_{s} \,{\boldsymbol \sigma }^{\textrm{T}} , \end{array} \end{equation}

which follows from (3.5) and ${\boldsymbol \sigma }^{\textrm{T}}{\boldsymbol{a}}^\alpha _{;\alpha }=2H{\boldsymbol \sigma }^{\textrm{T}}{\boldsymbol{n}}=\textbf{0}$ , which in turn follows from ${\boldsymbol{a}}^\alpha _{;\beta } = b^\alpha _\beta \,{\boldsymbol{n}}$ , $b^\alpha _\alpha =2H$ and (3.6).

In the expressions above, $(\ldots)_{;\alpha }$ denotes the so-called covariant derivative. It is equal to the parametric derivative $(\ldots)_{,\alpha }$ for objects without free index, such as $\boldsymbol{v}$ and $\boldsymbol \sigma$ , e.g. see Sauer (Reference Sauer2018).

3.3. Mesh ‘balance’

The mesh velocity ${\boldsymbol{v}}_{{m}}={\boldsymbol{v}}_{{m}}(\zeta ^\alpha,t)$ is governed by the following set of equations. Firstly, the condition

(3.9) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{{m}}\boldsymbol{\cdot} {\boldsymbol{n}} = {\boldsymbol{v}}\boldsymbol{\cdot} {\boldsymbol{n}}, \end{array} \end{equation}

which states that ${\boldsymbol{v}}_{{m}}$ must have the same normal component as $\boldsymbol{v}$ , and which essentially corresponds to an Lagrangian out-of-plane fluid description. Secondly, an equation for $v_{{m}}^\alpha (\zeta ^\alpha,t)$ , the in-plane component of ${\boldsymbol{v}}_{{m}}$ , is needed. The following options will be considered here:

1. (i) Prescribed mesh motion. In this case $v_{m}^\alpha$ is prescribed either directly, or indirectly by choosing $\zeta ^\alpha =\zeta ^\alpha (\xi ^\beta,t))$ and then calculating $v_{m}^\alpha$ from (2.21). Examples for the latter choice are given in § 6. A special case is prescribing the following suboption:

2. (ia) Zero in-plane mesh velocity. In this case

   (3.10) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{{m}}\boldsymbol{\cdot} {\boldsymbol{a}}_\alpha = 0. \end{array} \end{equation}
   This corresponds to an Eulerian in-plane fluid description. If desired, (3.9) and (3.10) can be combined into
   (3.11) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{{m}} = ({\boldsymbol{n}}\otimes {\boldsymbol{n}})\,{\boldsymbol{v}}. \end{array} \end{equation}
   It is noted that for evolving surfaces, an Eulerian description, just like a Lagrangian description, can lead to large mesh distortion. An example is shown in § 6.2, see figure 4.

3. (ii) Mesh motion defined by membrane elasticity. In this case $v_{{m}}^\alpha$ is characterised by the membrane PDE

   (3.12) \begin{equation} \begin{array}{l} \sigma ^{\alpha \beta }_{{m}\,;\beta } = 0, \end{array} \end{equation}
   that can be derived from (3.4) for the special choice $\rho =0$ , $\sigma ^{\alpha \beta }=\sigma ^{\alpha \beta }_{m}$ , $M^{\alpha \beta }=0$ and $\boldsymbol{f}=\textbf{0}$ . (Contracting (3.4) by ${\boldsymbol{a}}^\alpha$ and using (3.6)–(3.8) immediately leads to (3.12)). A possible elasticity model for the mesh is
   (3.13) \begin{equation} \begin{array}{l} \sigma ^{\alpha \beta }_{m} = \displaystyle \frac {\mu _{{m}}}{J_{{m}}}\left (A^{\alpha \beta }-a^{\alpha \beta }\right) \end{array}; \end{equation}
   (Sauer & Duong Reference Sauer and Duong2017). The parameter $\mu _{{m}}$ is physically irrelevant, but it can be used to improve the numerical behaviour. Equation (3.12) requires boundary conditions for the mesh position. An example is taking $v^\alpha _{m} = 0$ on the boundary. Equations (3.12)–(3.13) can then be solved at every time step for the current in-plane mesh position, see § 5.

The choice of option depends on the application at hand. If the surface does not deform much, or is even fixed, the in-plane Eulerian description (option ia) is best, as it is simplest. Examples are shown in §§ 6.3 and 6.5. On the other hand, if large surface deformations occur, the Eulerian description can become very inaccurate due to large local mesh distortions, and so the elastic mesh description (option ii) is best. This is demonstrated in the examples of §§ 6.2 and 6.6.

4.1. Pure in-plane flow

The first case considers in-plane flow without (out-of-plane) bending resistance ( $M^{\alpha \beta}\,{=}\,0$ ). From $\Psi _0 = q\,g$ , where $g:= J-1$ is the area-incompressibility constraint and $q$ its Lagrange multiplier, and the Newtonian surface fluid (sometimes also denoted a Boussinesq–Scriven surface fluid) model

(4.5) \begin{equation} \begin{array}{l} \sigma _{{inel}}^{\alpha \beta } = 2\eta \,d^{\alpha \beta }, \end{array} \end{equation}

follows

(4.6) \begin{equation} \begin{array}{l} \sigma ^{\alpha \beta } = q\,a^{\alpha \beta } + 2\eta \,d^{\alpha \beta }, \end{array} \end{equation}

see Sauer (Reference Sauer2016). Since $M^{\alpha \beta }=0$ , we find $N^{\alpha \beta }=\sigma ^{\alpha \beta }$ and $S^\alpha =0$ from (3.7). Hence the surface tension

(4.7) \begin{equation} \begin{array}{l} \gamma:= \frac {1}{2}\,N^{\alpha \beta } a_{\alpha \beta } , \end{array} \end{equation}

simply is $\gamma = q$ due to (2.48) and $\dot J = 0$ . With ${\boldsymbol \sigma } = \sigma ^{\alpha \beta }\,{\boldsymbol{a}}_\alpha \otimes {\boldsymbol{a}}_\beta$ and (2.38), the tensorial form of constitutive model (4.6) becomes

(4.8) \begin{equation} \begin{array}{l} {\boldsymbol \sigma } = \gamma \,\boldsymbol{i} + 2\eta \,\boldsymbol{d}_{{s}}. \end{array} \end{equation}

Here $\boldsymbol{i}:={\boldsymbol{a}}_\alpha \otimes {\boldsymbol{a}}^\alpha$ is the surface identity on $\mathcal{S}$ . Since $\boldsymbol{i}_{,\alpha }\,{\boldsymbol{a}}^\alpha = 2H{\boldsymbol{n}}$ , (3.8) then yields

(4.9) \begin{equation} \begin{array}{l} \boldsymbol{T}^\alpha _{;\alpha } = \boldsymbol{\nabla} _{\!s}\gamma + 2H\gamma \,{\boldsymbol{n}} + 2\eta \,\textrm{div}_{s} \,\boldsymbol{d}_{{s}}, \end{array} \end{equation}

with $\textrm{div}_{s} \,\boldsymbol{d}_{{s}} = \boldsymbol{d}_{{s},\alpha }\,{\boldsymbol{a}}^\alpha$ . This expression is used in some of the analytical examples of § 6.

4.2. In-plane flow with bending elasticity

The second case is like the first but with an additional out-of-plane bending energy based on the model of Helfrich (Reference Helfrich1973). In this case

(4.10a,b) \begin{equation} \begin{array}{l} \Psi _0 = q\,g + J\left (k (H-H_0)^2 + k_{g}\,\kappa \right)\kern-2.5pt,\quad g:= J-1. \end{array} \end{equation}

Here, $k$ and $k_{g}$ are bending moduli, and $H_0$ is the so-called spontaneous curvature – the value of the mean curvature $H$ that is energetically favourable. From (4.1)–(4.3) and (4.5) then follow

(4.11a) \begin{align}\sigma ^{\alpha \beta}&=\big(q + k\,\Delta H^2 - k_{g}\,\kappa \big)a^{\alpha \beta } - 2k\,\Delta H\,b^{\alpha \beta } + 2\eta \,d^{\alpha \beta }, \end{align}

(4.11b) \begin{align}M^{\alpha \beta}&=\left (k\,\Delta H+2k_{g} H\right)a^{\alpha \beta } - k_{g}\,b^{\alpha \beta }, \end{align}

(Sauer et al. Reference Sauer, Duong, Mandadapu and Steigmann2017) where $\Delta H:=H-H_0$ . Inserting this into (3.7a ) then gives

(4.12) \begin{equation} \begin{array}{l} N^{\alpha \beta } = \left (q + k\,\Delta H^2 \right)a^{\alpha \beta } - k\,\Delta H\,b^{\alpha \beta } + 2\eta \,d^{\alpha \beta } . \end{array} \end{equation}

In this case the physical surface tension defined in (4.7) is

(4.13) \begin{equation} \begin{array}{l} \gamma = q - k\,H_0\,\Delta H. \end{array} \end{equation}

The stress tensor can then be determined from (3.6).

It can be shown that for a sphere, where $b^{\alpha \beta }=Ha^{\alpha \beta }$ , (4.8) and (4.9) are still valid (but now need to be used with the new $\gamma$ ). From (4.11b ) follows for a sphere

(4.14) \begin{equation} \begin{array}{l} M^{\alpha \beta } = \left (k\,\Delta H + k_{g}\,H\right)\,a^{\alpha \beta }. \end{array} \end{equation}

5.1. Weak form for the fluid velocity

Multiplying PDE (3.4) with the admissible variation $\delta {\boldsymbol{x}}\in \mathcal{V}$ , integrating over the current surface and using the surface divergence theorem, leads to the weak form

(5.1) \begin{equation} \begin{array}{l} G:= G_{{in}} + G_{{int}} - G_{{ext}} = 0 \quad \forall \,\delta {\boldsymbol{x}}\in \mathcal{V}, \end{array} \end{equation}

with

(5.2a) \begin{align} G_{{in}} &= \displaystyle \int _{\mathcal{S}} \delta {\boldsymbol{x}}\boldsymbol{\cdot} \rho \,\dot {{\boldsymbol{v}}}\,\textrm{d} a, \end{align}

(5.2b) \begin{align} G_{{int}} &= \displaystyle \int _{\mathcal{S}} \sigma ^{\alpha \beta }\,\delta {\boldsymbol{x}}_{,\alpha }\boldsymbol{\cdot} {\boldsymbol{a}}_\beta \, \textrm{d} a + \int _{\mathcal{S}} M^{\alpha \beta }\,\delta {\boldsymbol{x}}_{;\alpha \beta }\boldsymbol{\cdot} {\boldsymbol{n}} \, \textrm{d} a,\end{align}

(5.2c) \begin{align} G_{{ext}} &= \displaystyle \int _{\mathcal{S}}\delta {\boldsymbol{x}}\boldsymbol{\cdot} \boldsymbol{f}\,\textrm{d} a + \int _{\partial \mathcal{S}}\delta {\boldsymbol{x}}\boldsymbol{\cdot} \boldsymbol{T}\,\textrm{d} s + \int _{\partial \mathcal{S}}\delta {\boldsymbol{n}}\boldsymbol{\cdot} {\boldsymbol{M}}\,\textrm{d} s,\end{align}

and

(5.3) \begin{equation} \begin{array}{l} \delta {\boldsymbol{x}}_{;\alpha \beta } = \delta {\boldsymbol{x}}_{,\alpha \beta } - \Gamma ^\gamma _{\alpha \beta }\,\delta {\boldsymbol{x}}_{,\gamma }, \end{array} \end{equation}

see Appendix D. Here, the integration is considered over the ALE parameterised surface ${\boldsymbol{x}}={\boldsymbol{x}}(\zeta ^\alpha,t)$ as this is required for a computational formulation in the ALE frame. Accordingly the surface variation $\delta {\boldsymbol{x}}$ is considered at fixed $\zeta ^\alpha$ , such that $\delta {\boldsymbol{x}}_{,\alpha } = \delta {\boldsymbol{a}}_\alpha$ , i.e. variation and parametric differentiation with respect to $\zeta ^\alpha$ commute. Thus, variation $\delta (\ldots)$ is mathematically equivalent to the time derivative $(\ldots)^{\prime}$ . An alternative would be to integrate over the Lagrangian surface parameterisation ${\boldsymbol{x}}=\hat {\boldsymbol{x}}(\xi ^\alpha,t)$ using a surface variation at fixed $\xi ^\alpha$ . Such a variation would be mathematically equivalent to the material time derivative $\dot {(\ldots)}$ , see Appendix D.

The second term in $G_{{int}}$ contains second derivatives (in both $M^{\alpha \beta }$ and $\delta {\boldsymbol{x}}_{;\alpha \beta }$ ) that can only be properly captured by at least $C^1$ -continuous ansatz functions. Examples in FEMs for the Helfrich model are subdivision surfaces (Feng & Klug Reference Feng and Klug2006) and NURBS – non-uniform rational B-splines (Sauer et al. Reference Sauer, Duong, Mandadapu and Steigmann2017). NURBS are used in the example of § 6.6.

Inserting (2.26), $G_{{in}}$ can be decomposed into the transient and convective parts

(5.4) \begin{equation} \begin{array}{r@{\,}l@{\,}l} G_{{trans}}&=&\displaystyle \int _{\mathcal{S}} \delta {\boldsymbol{x}}\boldsymbol{\cdot} \rho \,{\boldsymbol{v}}^{\prime}\,\textrm{d} a, \\[4mm] G_{{conv}} &=& \displaystyle \int _{\mathcal{S}} \delta {\boldsymbol{x}}\boldsymbol{\cdot} \rho \,{\boldsymbol{v}}_{,\alpha }\,\left (v^\alpha -v^\alpha _{m}\right)\,\textrm{d} a, \end{array} \end{equation}

such that $G_{{in}} = G_{{trans}} + G_{{conv}}$ .

If desired the integrals can be mapped to the reference configuration using $\textrm{d} a = J_{{m}}\,\textrm{d} A$ . If the fluid film does not resist bending, the internal and external bending moments $M^{\alpha \beta }$ and $\boldsymbol{M}$ are excluded from (5.2b ,5.2c ).

5.2. Weak form for the area-incompressibility constraint

Multiplying the area-incompressibility constraint (2.50) with the admissible variation $\delta q\in \mathcal{Q}$ and integrating over the current surface leads to the weak form

(5.5) \begin{equation} \begin{array}{l} \bar G:= \displaystyle \int _{\mathcal{S}_0}\delta q\,\textrm{div}_{s} \,{\boldsymbol{v}}\,\textrm{d} a = 0,\quad \forall \,\delta q \in \mathcal{Q}. \end{array} \end{equation}

5.3. Weak form for the mesh motion

The two components of the mesh velocity generally need to be treated differently:

1. (i) Out-of-plane mesh velocity: this component is defined by the strong form equation (3.9), which can either be enforced directly or converted into the weak form

   (5.6) \begin{equation} \begin{array}{l} \tilde G_{{o}}:= \displaystyle \int _{\mathcal{S}_0}w\,{\boldsymbol{n}}\boldsymbol{\cdot} \left ({\boldsymbol{v}}_{{m}}-{\boldsymbol{v}}\right)\,\textrm{d} A = 0 ,\quad \forall \,w \in \mathcal{W}_{{n}}. \end{array} \end{equation}
   Here, the integration is written over the reference configuration as this simplifies computations. Alternatively one can also chose to integrate over the current configuration.

2. (ii) In-plane mesh velocity: this component can either be prescribed, e.g. as zero, see (3.10) or it can be obtained from membrane equilibrium (3.12). If (3.12) is used, it is advantageous to derive the weak form by integration over the current ALE surface using $\textrm{d} a = J_{{m}}\,\textrm{d} A$ . This leads to the weak form

   (5.7) \begin{equation} \begin{array}{l} \tilde G_{{iel}}:= \displaystyle \int _{\mathcal{S}_0}w_{\alpha;\beta }\,\sigma _{{m}}^{\alpha \beta }\,J_{{m}}\,\textrm{d} A = 0 ,\quad \forall \,w_\alpha \in \mathcal{W}_\alpha . \end{array} \end{equation}
   This follows from (5.2b ) when $\delta {\boldsymbol{x}}$ is replaced by the test function $\boldsymbol{w}=w_\alpha \,{\boldsymbol{a}}^\alpha$ , see also (187.3) in Sauer & Duong (Reference Sauer and Duong2017). For the elastic membrane stress in (3.13), $J_{{m}}$ can then be conveniently cancelled. If (3.10) is used, it can be either prescribed in strong form or by the weak form
   (5.8) \begin{equation} \begin{array}{l} \tilde G_{{i}0}:= \displaystyle \int _{\mathcal{S}_0}w_\alpha \,{\boldsymbol{a}}^\alpha \boldsymbol{\cdot} {\boldsymbol{v}}_{{m}}\,\textrm{d} A = 0 ,\quad \forall \,w_\alpha \in \mathcal{W}_\alpha . \end{array} \end{equation}

The last statement combines with (5.6) to

(5.9) \begin{equation} \begin{array}{l} \tilde G_0:= \displaystyle \int _{\mathcal{S}_0}\boldsymbol{w}\boldsymbol{\cdot} \left ({\boldsymbol{v}}_{{m}}-({\boldsymbol{n}}\otimes {\boldsymbol{n}})\,{\boldsymbol{v}}\right)\,\textrm{d} A = 0 ,\quad \forall \,\boldsymbol{w}\in \mathcal{W}, \end{array} \end{equation}

which is the weak form of (3.11). Combining (5.6) and (5.7), on the other hand, yields

(5.10) \begin{equation} \begin{array}{l} \tilde G_{{el}}:= \alpha _{{m}} \displaystyle \int _{\mathcal{S}_0}w\,{\boldsymbol{n}}\boldsymbol{\cdot} \left ({\boldsymbol{v}}_{{m}}-{\boldsymbol{v}}\right)\,\textrm{d} A + \int _{\mathcal{S}_0}w_{\alpha;\beta }\,\tau _{{m}}^{\alpha \beta }\,\textrm{d} A = 0 ,\quad \forall \,\boldsymbol{w}\in \mathcal{W}, \end{array} \end{equation}

where $\tau _{m}^{\alpha \beta }:=J_{{m}}\,\sigma _{m}^{\alpha \beta }$ . Here, the factor $\alpha _{{m}}$ is included to ensure dimensional consistency between the two terms, which is needed in case $\tau ^{\alpha \beta }_{m}$ (via $\mu _{{m}}$ ) is chosen to have units of membrane stress (force/length). In (5.6)–(5.10), $\boldsymbol{w}=w_\alpha \,{\boldsymbol{a}}^\alpha + w\,{\boldsymbol{n}}$ is the test function corresponding to mesh velocity ${\boldsymbol{v}}_{{m}}$ .

5.4. Summary

The above weak form statements can be solved for the three fields ${\boldsymbol{v}}(\zeta ^\alpha,t)$ , $q(\zeta ^\alpha,t)$ and ${\boldsymbol{v}}_{{m}}(\zeta ^\beta,t)$ . In Appendix F a stabilised FEM is outlined that uses quadratic shape functions for the spatial discretisation and the implicit trapezoidal rule for the temporal discretisation. This leads to an algebraic system of three coupled equations that can be solved monolithically using the Newton–Raphson method. If the mesh velocity is prescribed or eliminated via (3.9) and (3.10), the algebraic system only consists of two coupled equations.

The fields $\boldsymbol{v}$ , $q$ and ${\boldsymbol{v}}_{{m}}$ , along with their variations $\delta {\boldsymbol{x}}$ , $\delta q$ and $\boldsymbol{w}$ , have to be picked such that all integrals above are well defined. This is the case for piecewise polynomial functions, as are classically considered in FEMs.

It is noted that in the above formulation no unknown mesh pressure appears. Such a variable can be used to enforce constraint (3.9) in the framework of Lagrange multiplier methods (Sahu Reference Sahu2024). It adds an extra field to the system that needs to be accounted for in the out-of-plane mesh equation. Here instead, constraint (3.9) is enforced through weak form (5.6) without introducing such an extra variable. Consequently, the out-of-plane mesh equation is not needed either.

6.1. One-dimensional fluid flow

The first example considers the one-dimensional fluid motion shown in figure 2. It consists of a flat fluid sheet extruded with constant velocity from the left boundary. Even though the example is very simple, it illustrates the proposed ALE formulation in the curvilinear setting. It is also a first step towards the following bubble example. The motion in figure 2 can be described by

(6.1) \begin{equation} \begin{array}{l} {\boldsymbol{x}} = x\,{\boldsymbol{e}}_1 + y\,{\boldsymbol{e}}_2, \end{array} \end{equation}

where position $x = \hat x(\xi ^1,t) = x(\zeta ^1,t) = \tilde x(\theta ^1,t)$ depends on the parameterisation, while $y$ is equal for all, i.e. $y=\xi ^2 = \zeta ^2 = \theta ^2$ . To shorten notation we let $\xi:=\xi ^1$ , $\zeta:=\zeta ^1$ and $\theta = \theta ^1$ . The Lagrangian motion is given by

(6.2) \begin{equation} \begin{array}{l} x = \hat x(\xi,t) = x_0(t) + (1-\xi)\,H_0, \end{array} \end{equation}

with

(6.3) \begin{equation} \begin{array}{l} x_0(t) = v_{{in}}\,t, \end{array} \end{equation}

such that

(6.4) \begin{equation} \begin{array}{l} X(\xi) = \hat x(\xi,0) = (1-\xi)\,H_0 \end{array} \end{equation}

is the initial position. Here, length $H_0$ , inflow velocity $v_{{in}}$ and directions ${\boldsymbol{e}}_\alpha$ are constants. Choosing the ALE coordinate

Figure 2.

One-dimensional fluid flow: (a) initial and current configuration, and illustration of (2.20). The Lagrangian and ALE surface parameterisations are shown in blue and green (coordinates $\xi$ and $\zeta$ ), respectively. (b) Position $x = \hat x(\xi)=x(\zeta)=\tilde x(\theta)$ for $x_0 = H_0$ (at $t= H_0/v_{{in}}$ ). Note that coordinates $\zeta$ , $\xi$ and $\theta$ run from right to left here.

(6.5) \begin{equation} \begin{array}{l} \zeta = \displaystyle \frac {H_0}{h_0}\,\xi , \end{array} \end{equation}

leads to the ALE parametrisation

(6.6) \begin{equation} \begin{array}{l} x(\zeta,t) = (1-\zeta)\,h_0(t), \end{array} \end{equation}

with

(6.7) \begin{equation} \begin{array}{l} h_0(t):= H_0 + x_0(t), \end{array} \end{equation}

such that

(6.8) \begin{equation} \begin{array}{l} X(\zeta) = x(\zeta,0) = (1-\zeta)\,H_0. \end{array} \end{equation}

Since $\dot h_0 = h_0^{\prime} = x_0^{\prime} = \dot x_0 = v_{{in}}$ , one finds the velocities

(6.9) \begin{align} \begin{array}{l} {\boldsymbol{v}} = \dot {\boldsymbol{x}} = v_{{in}}\,{\boldsymbol{e}}_1, \end{array} \end{align}

(6.10) \begin{align} \begin{array}{l} {\boldsymbol{v}}_{{m}} = {\boldsymbol{x}}^{\prime} = (1-\zeta)\,v_{{in}}\,{\boldsymbol{e}}_1 \end{array} \end{align}

and

(6.11) \begin{equation} \begin{array}{l} \dot \zeta ^1 = \dot \zeta = -\zeta \displaystyle \frac {v_{{in}}}{h_0},\quad \dot \zeta ^2 = 0, \end{array} \end{equation}

from (6.1), (6.2), (6.6) and (6.5), respectively. The tangent vectors

(6.12) \begin{equation} \begin{array}{l} {\boldsymbol{a}}_1 = \displaystyle \frac {\partial {{\boldsymbol{x}}}}{\partial {\zeta ^1}} = -h_0\,{\boldsymbol{e}}_1,\quad {\boldsymbol{a}}_2 = \displaystyle \frac {\partial {{\boldsymbol{x}}}}{\partial {\zeta ^2}} = {\boldsymbol{e}}_2, \end{array} \end{equation}

then give

(6.13) \begin{equation} \begin{array}{l} \dot \zeta ^\alpha {\boldsymbol{a}}_\alpha = \zeta \,v_{{in}}\,{\boldsymbol{e}}_1. \end{array} \end{equation}

The velocities $\boldsymbol{v}$ , ${\boldsymbol{v}}_{{m}}$ and $\dot \zeta ^\alpha {\boldsymbol{a}}_\alpha$ are visualised in figure 2 (by showing the displacements ${\boldsymbol{v}}\,t$ , ${\boldsymbol{v}}_{{m}}\,t$ and $\dot \zeta ^\alpha {\boldsymbol{a}}_\alpha \,t$ ). It is easy to confirm that they satisfy (2.20), which states that the material velocity $\boldsymbol{v}$ is composed of the frame velocity ${\boldsymbol{v}}_{{m}}$ and the frame-relative material velocity $\dot \zeta ^\alpha {\boldsymbol{a}}_\alpha$ . From ${\boldsymbol{v}}_{\!,\alpha }=\textbf{0}$ and $\dot {\boldsymbol{a}}_1 = -v_{{in}}\,{\boldsymbol{e}}_1$ and $\dot {\boldsymbol{a}}_2=\textbf{0}$ and

(6.14) \begin{equation} \begin{array}{l} \dot \zeta ^1_{,1} = -\displaystyle \frac {v_{{in}}}{h_0}, \end{array} \end{equation}

one can also confirm that (2.27) is satisfied. This underlines that ${\boldsymbol{v}}_\alpha \neq \dot {\boldsymbol{a}}_\alpha$ here. The equality of ${\boldsymbol{v}}_\alpha$ and $\dot {\boldsymbol{a}}_\alpha$ , used in the curvilinear ALE formulation of Sahu et al. (Reference Sahu, Omar, Sauer and Mandadapu2020), is therefore not generally satisfied. Instead they differ here proportional to the inflow velocity $v_{{in}}$ .

The Eulerian motion is given by

(6.15) \begin{equation} \begin{array}{l} x = \tilde x(\theta) = (1-\theta)\,H_0, \end{array} \end{equation}

since it leads to ${\boldsymbol{v}}_{{m}}=\textbf{0}$ . Figure 2(b) illustrates the Lagrangian, ALE and Eulerian description of motion.

6.2. Inflation of a two-dimensional soap bubble

Next, the inflation of a two-dimensional soap bubble as shown in figure 3 is considered. Apart from the kinematics of inflation, the example provides insight to how Lagrangian and Eulerian frame motion leads to severely distorted surface grids, even in simple flow examples, while these distortions can be avoided with the proposed ALE formulation, which then leads to superior numerical methods, as is also shown.

Figure 3.

Inflation of two-dimensional soap bubble: (a) initial and (b) current configurations. Fluid inflow is considered at the left boundary such that $J\equiv 1$ (i.e. $\textrm{div}_{s} \,{\boldsymbol{v}}\equiv 0$ ) over time. The inflation can be driven by prescribing $v_{{in}}$ or $V(t)$ at the opening. The Lagrangian and ALE surface parameterisations are shown in blue and green (coordinates $\xi$ and $\zeta$ ), respectively.

6.2.1. Surface motion

As in the example before, no flow in ${\boldsymbol{e}}_3$ direction occurs and we again use $\xi =\xi ^1$ , $\zeta =\zeta ^1$ and $\theta =\theta ^1$ as shorthand for the main parameter, while $z=\xi ^2=\zeta ^2=\theta ^2$ . The motion of the bubble is described by

(6.16) \begin{equation} \begin{array}{l} {\boldsymbol{x}} = x_0\,{\boldsymbol{e}}_1 + r\,{\boldsymbol{e}}_r + z\,{\boldsymbol{e}}_3, \end{array} \end{equation}

where the centre

(6.17) \begin{equation} \begin{array}{l} x_0 = -r\,\cos \psi _{\textrm{o}}, \end{array} \end{equation}

and radius

(6.18) \begin{equation} \begin{array}{l} r = \displaystyle \frac {L}{\sin \psi _{\textrm{o}}}, \end{array} \end{equation}

depend on the fixed opening height $L$ and the varying opening angle $\psi _{\textrm{o}}=\psi _{\textrm{o}}(t)$ , and therefore are only functions of time, i.e. $x_0=x_0(t)$ and $r=r(t)$ . The only spatially varying objects in (6.16) are $z$ and

(6.19a,b) \begin{equation} \begin{array}{l} {\boldsymbol{e}}_r = \cos \psi \,{\boldsymbol{e}}_1 + \sin \psi \,{\boldsymbol{e}}_2,\quad 0\leq \psi \leq \psi _{\textrm{o}}\lt \pi . \end{array} \end{equation}

The flow is fully described by the Lagrangian parameterisation of angle $\psi$

(6.20) \begin{equation} \begin{array}{l} \psi = \hat \psi (\xi,t) = \displaystyle \frac {R}{r}\Psi _{\textrm{o}}\,\xi , \end{array} \end{equation}

where $R = r(0)$ and $\Psi _{\textrm{o}}=\psi _{\textrm{o}}(0)$ , see figure 3. A possible ALE parameterisation that satisfies $\psi (0,t)=0$ and $\psi (1,t)=\psi _{\textrm{o}}$ is

(6.21) \begin{equation} \begin{array}{l} \psi = \psi (\zeta,t) = \psi _{\textrm{o}}\,\zeta . \end{array} \end{equation}

From this follows

(6.22) \begin{equation}\zeta (\xi,t) = \displaystyle \frac {R\Psi _{\textrm{o}}}{r\psi _{\textrm{o}}}\xi, \end{equation}

and

(6.23) \begin{equation} \begin{array}{l} {\boldsymbol{a}}_1 = r\,\psi _{\textrm{o}}\,{\boldsymbol{e}}_\psi , \end{array} \end{equation}

where

(6.24) \begin{equation} \begin{array}{l} {\boldsymbol{e}}_\psi = -\sin \psi \,{\boldsymbol{e}}_1 + \cos \psi \,{\boldsymbol{e}}_2, \end{array} \end{equation}

and (E1) from Appendix E have been used. Basis ${\boldsymbol{a}}_2$ on the other hand is constant and equal to ${\boldsymbol{e}}_3$ . The function $\zeta (\xi,t)$ defines the in-plane material motion relative to the ALE frame, just like in the example before.

Remark 6.1. Coordinates $\xi$ and $\zeta$ parameterise angle $\psi$ according to (6.20) and (6.21). Thus, they also parameterise the arc length $r\psi$ , but are not necessarily equal to it, i.e. generally $R\Psi _{\operatorname{o}}\neq 1$ and $r\Psi _{\operatorname{o}}\neq 1$ .

6.2.2. Surface velocity

Since $\psi _{\textrm{o}}$ , $r$ and $x_0$ are only functions of time, there is no difference between their time derivatives $\dot {(\ldots)}$ and $(\ldots)^{\prime}$ . Thus, not all time derivatives differ between the Lagrangian and ALE frame. Given $\dot \psi _{\textrm{o}}$ (that can be prescribed to control the inflation, see below), one finds

(6.25) \begin{equation} \begin{array}{l} \dot x_0 = \displaystyle \frac {1}{\sin \psi _{\textrm{o}}}\,r\,\dot \psi _{\textrm{o}} , \end{array} \end{equation}

and

(6.26) \begin{equation} \begin{array}{l} \dot r = -\displaystyle \frac {\cos \psi _{\textrm{o}}}{\sin \psi _{\textrm{o}}}\,r\,\dot \psi _{\textrm{o}}, \end{array} \end{equation}

from (6.17) and (6.18). Further, (6.20) and (6.21) lead to

(6.27) \begin{align} \begin{array}{l} \dot \psi = -\displaystyle \frac {\dot r}{r}\psi , \end{array} \end{align}

(6.28) \begin{align} \begin{array}{l} \psi ^{\prime} = \dot \psi _{\textrm{o}}\,\zeta \end{array} \end{align}

and

(6.29) \begin{equation} \begin{array}{l} \dot \zeta = \displaystyle \frac {\dot \psi }{\psi _{\textrm{o}}}-\zeta \,\frac {\dot \psi _{\textrm{o}}}{\psi _{\textrm{o}}}. \end{array} \end{equation}

It can be confirmed that (6.27)–(6.29) satisfy (2.30) as $\partial \psi /\partial \zeta =\psi _{\textrm{o}}$ . Inserting (6.27) and (6.21), derivative $\dot \zeta$ , which describes the in-plane material velocity relative to the ALE frame, can also be written as

(6.30) \begin{equation} \begin{array}{l} \dot \zeta = -\displaystyle \frac {\dot {(r\psi _{\textrm{o}})}}{r\psi _{\textrm{o}}}\,\zeta . \end{array} \end{equation}

Using (E2), and writing

(6.31) \begin{equation} \begin{array}{l} {\boldsymbol{e}}_1 = \cos \psi \,{\boldsymbol{e}}_r - \sin \psi \,{\boldsymbol{e}}_\psi , \end{array} \end{equation}

the material velocity field follows as

(6.32) \begin{equation} \begin{array}{l} {\boldsymbol{v}} = \dot {\boldsymbol{x}} = (\dot x_0\,\cos \psi + \dot r)\,{\boldsymbol{e}}_r + (r\,\dot \psi -\dot x_0\,\sin \psi)\,{\boldsymbol{e}}_\psi , \end{array} \end{equation}

while the ALE frame velocity becomes

(6.33) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{{m}} = {\boldsymbol{x}}^{\prime} = (\dot x_0\,\cos \psi + \dot r)\,{\boldsymbol{e}}_r + (r\,\psi ^{\prime}-\dot x_0\,\sin \psi)\,{\boldsymbol{e}}_\psi . \end{array} \end{equation}

Together with (6.23), (6.28), (6.29) and $\dot \zeta ^2=0$ , they correctly satisfy (2.20). As seen, $\boldsymbol{v}$ and ${\boldsymbol{v}}_{{m}}$ share the same normal velocity component

(6.34) \begin{equation} \begin{array}{l} v_r = (\cos \psi - \cos \psi _{\textrm{o}})\,\displaystyle \frac {r\,\dot \psi _{\textrm{o}}}{\sin \psi _{\textrm{o}}}, \end{array} \end{equation}

along ${\boldsymbol{e}}_r$ (as they should), while they differ in their tangential velocity components

(6.35) \begin{equation} \begin{array}{l} v_\psi = r\,\psi _{\textrm{o}}\,\dot \zeta + \left (\zeta -\displaystyle \frac {\sin \psi }{\sin \psi _{\textrm{o}}}\right)\,r\,\dot \psi _{\textrm{o}}, \end{array} \end{equation}

and

(6.36) \begin{equation} \begin{array}{l} v_{{m}\psi } = \left (\zeta -\displaystyle \frac {\sin \psi }{\sin \psi _{\textrm{o}}}\right)\,r\,\dot \psi _{\textrm{o}}, \end{array} \end{equation}

along ${\boldsymbol{e}}_\psi$ . These expressions follow from inserting (6.25)–(6.26) and (6.28)–(6.29) into (6.32)–(6.33). At the inflow, where $\zeta =1$ and $\psi =\psi _{\textrm{o}}$ , the velocity components $v_{{m}\phi }$ and $v_r$ are zero (the latter since $L$ is fixed), while the tangential velocity $-v_\psi$ is equal to the inflow velocity $v_{{in}}$ (see figure 3), giving

(6.37) \begin{equation} \begin{array}{l} v_{{in}} = \left.-r\,\psi _{\textrm{o}}\,\dot \zeta \right |_{\zeta =1}, \end{array} \end{equation}

which in light of (6.30) becomes

(6.38) \begin{equation} \begin{array}{l} v_{{in}} = \dot {(r\psi _{\textrm{o}})}. \end{array} \end{equation}

Figure 4.

Inflation of a two-dimensional soap bubble: (a) Lagrangian, (b) ALE and (c) Eulerian surface motion and corresponding surface velocity fields $\boldsymbol{v}$ , ${\boldsymbol{v}}_{{m}}$ and ${\boldsymbol{v}}_{{e}}$ . Shown is the motion of 17 initially equidistant surface points. They remain equidistant in the Lagrangian and ALE motion, but only the former maintains their original distance. In case of the Lagrangian motion, new material points are drawn in at the boundary (shown by path lines three times more spaced). The Eulerian motion, which is always normal to the surface, leads to unequal distances between points. The initial opening angle is chosen as $\Psi_{\textrm{o}} = 15^\circ$ . The black lines show the bubble at $\psi_{\textrm{o}} = [\Psi_{\textrm{o}},\,25^\circ,\,40^\circ,\,57.5^\circ,\,80^\circ,\,105^\circ,\,125^\circ,\,140^\circ,\,150^\circ]$ . Here, $L$ is used for normalising the geometry.

Figure 5.

Inflation of two-dimensional soap bubble: (a) radial and tangential velocity components $v_r(\zeta)$ , $v_\psi (\zeta)$ and $v_{{m}\psi }(\zeta)$ , and (b) angle $\psi$ in Lagrangian, $\hat \psi (\xi)$ , ALE, $\psi (\zeta)$ and Eulerian, $\tilde \psi (\theta)$ , description. All for $\psi_{\textrm{o}} = 150^\circ$ . The coordinates are running from right to left as in figure 3. The Lagrangian coordinate continues until $\xi = \psi_{\textrm{o}}\sin \Psi_{\textrm{o}}/(\Psi_{\textrm{o}}\sin \psi_{\textrm{o}}) \approx 5.176$ according to (6.22) and (6.18). The left figure shows how the velocity components vary across the surface: $v_r$ is maximum at the tip, and $-v_\psi$ at the inflow. Note $v_\psi = v_{{m}\psi } + \dot \zeta \,r\,\psi_{\textrm{o}}$ . The circles in (b) mark the positions shown in figure 4.

The Lagrangian and the ALE motion are visualised in figure 4, while figure 5 shows their velocity components. Figure 4 also shows the Eulerian surface motion, which is determined from the condition $v_{{m}\psi }=0$ for all $t$ . This leads to a nonlinear ODE for $\psi =\tilde \psi (\theta,t)$ that together with (6.16) and (6.19a , b ) then defines the Eulerian surface motion: setting the tangential component of (6.33) to zero and inserting (6.25) yields

(6.39) \begin{equation} \begin{array}{l} \displaystyle \frac {\psi ^{\prime}}{\sin \psi } = \displaystyle \frac {\psi _{\textrm{o}}^{\prime}}{\sin \psi _{\textrm{o}}}, \end{array} \end{equation}

since $\dot \psi _{\textrm{o}} = \psi _{\textrm{o}}^{\prime}$ . Multiplying by $\textrm{d} t$ and considering $\zeta =\theta$ fixed then gives

(6.40) \begin{equation} \begin{array}{l} \displaystyle \frac {\textrm{d}\psi }{\sin \psi } = \displaystyle \frac {\textrm{d}\psi _{\textrm{o}}}{\sin \psi _{\textrm{o}}}, \end{array} \end{equation}

which can be integrated on both sides to yield the ODE

(6.41) \begin{equation} \begin{array}{l} \displaystyle \frac {1+\cos \psi }{\sin \psi } = c(\theta)\,\displaystyle \frac {1+\cos \psi _{\textrm{o}}}{\sin \psi _{\textrm{o}}}. \end{array} \end{equation}

The integration constant $c(\theta)$ follows from evaluating (6.41) at $t=0$ , where $\psi _{\textrm{o}} = \Psi _{\textrm{o}}$ and $\psi = \Psi _{\textrm{o}}\,\theta:=\Psi (\theta)$ according to (6.21). This gives

(6.42) \begin{equation} \begin{array}{l} c(\theta) = \displaystyle \frac {1+\cos \Psi (\theta)}{1+\cos \Psi _{\textrm{o}}}\frac {\sin \Psi _{\textrm{o}}}{\sin \Psi (\theta)}. \end{array} \end{equation}

Picking a $\theta$ , one can then evaluate $c(\theta)$ and solve nonlinear (6.41) for $\psi = \tilde \psi (\theta,t)$ at every $t$ . The resulting angle $\tilde \psi (\theta,t)$ and motion $\tilde {\boldsymbol{x}}(\theta,t)$ are shown in figure 5(b) and 4(c). As the latter figure shows, the Eulerian motion is normal to the surface, and hence it has no tangential component, i.e. the Eulerian frame velocity is ${\boldsymbol{v}}_{{e}} = v_r\,{\boldsymbol{e}}_r$ . Due to this, grid points become widely separated at the tip of the bubble, as figure 4(c) shows. This causes problems in numerical methods as is seen in § 6.2.4 below. Also Lagrangian motion is problematic in numerical methods, as new grid points need to be drawn in at the inflow boundary, as figure 4(a) shows. On the other hand, ALE motion achieves equidistant grid points during inflation, which allows for superior accuracy in numerical methods (see § 6.2.4), which in turn allows for more accurate studies of evolving surface flows.

Next the velocity gradient and its consequences are examined. From (6.32), (6.27), (E1) and $\psi _{,1} = \psi _{\textrm{o}}$ follows

(6.43) \begin{equation} \begin{array}{l} {\boldsymbol{v}}_{\!,1} = \dot r\,\psi \,\psi _{\textrm{o}}\,{\boldsymbol{e}}_r. \end{array} \end{equation}

Further,

(6.44) \begin{equation} \begin{array}{l} \dot {\boldsymbol{a}}_1 = \dot r\, \psi \,\psi _{\textrm{o}}\,{\boldsymbol{e}}_r + v_{{in}}\,{\boldsymbol{e}}_\psi, \end{array}\end{equation}

and

(6.45) \begin{equation} \begin{array}{l} \dot \zeta _{,1} = -\displaystyle \frac {v_{{in}}}{r\psi _{\textrm{o}}}, \end{array}\end{equation}

follow from (6.23), (6.27), (6.38), (E2) and (6.30). On the other hand ${\boldsymbol{v}}_{\!,2}=\dot {\boldsymbol{a}}_2=\textbf{0}$ and $\dot \zeta _{,2}=0$ . Equation (2.27) is thus satisfied. With

(6.46) \begin{equation} \begin{array}{l} {\boldsymbol{a}}^{\prime}_1 = - r\,\psi \,\dot \psi _{\textrm{o}}\,{\boldsymbol{e}}_r + v_{{in}}\,{\boldsymbol{e}}_\psi , \end{array} \end{equation}

${\boldsymbol{a}}_{1,1}=-r\psi ^2_{\textrm{o}}\,{\boldsymbol{e}}_r$ and ${\boldsymbol{a}}_2^{\prime}={\boldsymbol{a}}_{1,2}={\boldsymbol{a}}_{2,2}=\textbf{0}$ , one can also confirm that (2.30) is satisfied for ${\boldsymbol{a}}_\alpha$ . This follows from (6.21), (6.23), (6.28), (6.30), (6.38), (E1) and (E2).

With (6.43), ${\boldsymbol{a}}^1 = {\boldsymbol{e}}_\psi /(r\psi _{\textrm{o}})$ and (2.16), one can then confirm that the flow satisfies (2.50), i.e. it is area incompressible, and one can find the rate-of-deformation tensor

(6.47) \begin{equation}2\boldsymbol{d} = -\dot \psi\big({\boldsymbol{e}}_r\otimes {\boldsymbol{e}}_\psi + {\boldsymbol{e}}_\psi \otimes {\boldsymbol{e}}_r \big),\end{equation}

according to (2.34) and (6.27), which implies $\boldsymbol{d}_{{s}} = \textbf{0}$ . (Tangential velocity fields that lead to $\boldsymbol{d}_{{s}} = \textbf{0}$ are sometimes referred to (Wilhelm) Killing vector fields, e.g. see Jankuhn et al. (Reference Jankuhn, Olshanskii and Reusken2018)). So the stress according to constitutive model (4.8) will only consist of the surface tension $\gamma$ that equilibrates the internal pressure $p$ through the well known relation $p=\gamma /r$ . Thus, if inertia is negligible and no additional surface forces act, the soap bubble remains perfectly circular during inflation.

6.2.3. Inflation control

The inflation can be controlled in various ways. One can prescribe inflow velocity $v_{{in}}$ , giving the opening angle change

(6.48) \begin{equation} \begin{array}{l} \dot \psi _{\textrm{o}} = \displaystyle \frac {\sin \psi _{\textrm{o}}}{r\,(\sin \psi _{\textrm{o}} - \psi _{\textrm{o}}\,\cos \psi _{\textrm{o}})}\,v_{{in}}, \end{array}\end{equation}

from (6.38) and (6.26). This then specifies the bubble motion according to the equations above. Prescribing $v_{{in}}$ corresponds to pushing material into the bubble, which may not be a stable way to inflate the bubble numerically. A stable alternative is to prescribe the bubble volume (per unit depth) $V$ since this pulls, instead of pushes, material into the bubble. V(t) is related to the opening angle $\psi _{\textrm{o}}$ via

(6.49) \begin{equation} \begin{array}{l} V = r^2\,(\psi _{\textrm{o}} - \sin \psi _{\textrm{o}}\,\cos \psi _{\textrm{o}}), \end{array} \end{equation}

according to figure 3(b). Taking a time derivative of this and using (6.26) then gives

(6.50) \begin{equation} \begin{array}{l} \dot \psi _{\textrm{o}} = \displaystyle \frac {\sin \psi _{\textrm{o}}}{2r^2\,(\sin \psi _{\textrm{o}}-\psi _{\textrm{o}}\,\cos \psi _{\textrm{o}})}\,\dot V. \end{array} \end{equation}

Comparing (6.48) and (6.50) shows that volume and inflow velocity are related by $\dot V=2r\,v_{{in}}$ . Figure 6 shows the change $\psi _{\textrm{o}}$ and $r$ due to given $V$ .

Figure 6.

Inflation of a two-dimensional soap bubble: (a) $\psi _{\textrm{o}}(V)$ and (b) $r(\psi _{\textrm{o}})$ vs. $V(\psi _{\textrm{o}})$ . The circles mark the configurations shown in figure 4.

Instead of the volume, the internal pressure $p(t)$ can also be prescribed. However, in contrast to $V(t)$ , $p(t)$ generally cannot be chosen as a monotonically increasing function. This is due to the minimum attained by $r$ during the inflation process (see figure 6), which leads to a maximum in $p$ for constant surface tension $\gamma$ , due to the relation $\gamma = \textit{rp}$ . Prescribing $p$ numerically therefore can also become unstable.

6.2.4. Numerical solution

It is finally shown that the equidistant grid spacing of the ALE description leads to much more accurate numerical results than the unequal grid spacing of the Eulerian description. Therefore the two-dimensional soap bubble is solved with the FE method using either Eulerian mesh motion (described by weak form (5.9)) or elastic mesh motion (described by weak form (5.10)). The latter leads to an equidistant FE mesh for this problem. Quadratic Lagrange interpolation is used for all fields as described in Appendix F. The physical soap bubble properties in the simulation are taken as $L = 10$ mm, $\gamma = 25\,\unicode{x03BC}$ N mm⁻¹, $\eta = 0.04\,\unicode{x03BC}$ N ms mm)⁻¹ and $\rho = 0.2\,\unicode{x03BC}$ Nms $^2$ mm $^{-3}$ , while the mesh stiffness is taken as $\mu _{{m}} = 0.5\,$ nN mm $^{-2}$ .

Figure 7(a) shows the FE results for velocity components $v_r$ and $v_\psi$ in comparison with the analytical result of (6.34) and (6.35). As seen, the elastic ALE result is much more accurate, even for a coarser mesh. This can also be seen in the error plot of figure 7(b): the ALE result is around 10 times more accurate for $v_r$ and more than 100 times more accurate for $v_\psi$ for coarse meshes. The difference decreases for finer meshes, but it remains positive. Interesting, the convergence rate for the tangential velocity $v_\psi$ decreases for finer meshes, which calls for further study. Also interesting is that the convergence behaviour of $v_\psi$ in the Eulerian case is reminiscent of membrane locking in thin shells (Sauer, Zou & Hughes Reference Sauer, Zou and Hughes2024).

Figure 7.

Inflation of a two-dimensional soap bubble: (a) FE result and (b) FE convergence of velocity components $v_r$ and $v_\psi$ for Eulerian and elastic ALE mesh motion. The latter is around 10 times more accurate in $v_r$ and over 100 times more accurate in $v_\psi$ for coarse meshes. The vertical grey lines in (a) show the FE boundaries: dashed for the Eulerian mesh with 32 FE, solid for the ALE mesh with 16 FE (which is exactly the mesh in figure 4 b). The initial FE convergence rate is at least $O(n_{{el}}^{-2}) = O(h^2)$ .

6.3. Shear flow on a sphere

The third example considers simple shear flow on a sphere. The example provides insight into the coupling of in-plane surface flow with out-of-plane surface forces and deformations. The latter remain small here and hence the ALE mesh motion can be set to zero here, i.e. the ALE frame coincides with an in-plane Eulerian one. It is further shown that this is a much better way to treat the example than using a Lagrangian frame.

6.3.1. Surface description

The spherical surface can be described by

(6.51) \begin{equation} \begin{array}{l} {\boldsymbol{x}} = r\,{\boldsymbol{e}}_r, \end{array} \end{equation}

where radius $r$ is now considered fixed and the direction vectors

(6.52a) \begin{align}{\boldsymbol{e}}_r &:= \cos \zeta ^2\,{\boldsymbol{e}}_{r_0} + \sin \zeta ^2\,{\boldsymbol{e}}_3,\quad -\frac {\pi }{2}\leq \zeta ^2\leq \frac {\pi }{2},\end{align}

(6.52b) \begin{align} {\boldsymbol{e}}_{r_0}&:= \cos \zeta ^1\,{\boldsymbol{e}}_1 + \sin \zeta ^1\,{\boldsymbol{e}}_2,\quad 0\leq \zeta ^1\leq 2\pi,\end{align}

generally depend on $\zeta ^\alpha = \zeta ^\alpha (\xi ^\beta,t)$ , see figure 16 in Appendix E. Here, $\{{\boldsymbol{e}}_1,{\boldsymbol{e}}_2,{\boldsymbol{e}}_3\}$ is the basis of a fixed Cartesian reference coordinate system. Before specifying $\zeta ^\alpha (\xi ^\beta,t)$ , which defines the material flow relative to the ALE frame (= Eulerian frame here), one can already analyse the geometry. Introducing the local orthonormal basis $\{{\boldsymbol{e}}_r,{\boldsymbol{e}}_\phi,{\boldsymbol{e}}_\theta \}$ with

(6.53a) \begin{align}{\boldsymbol{e}}_\phi &:=-\sin \phi \,{\boldsymbol{e}}_1 + \cos \phi \,{\boldsymbol{e}}_2, \quad \phi := \zeta ^1,\end{align}

(6.53b) \begin{align} {\boldsymbol{e}}_\theta &:=-\sin \theta \,{\boldsymbol{e}}_{r_0} + \cos \theta \,{\boldsymbol{e}}_3, \quad \theta := \zeta ^2, \end{align}

and using (E4) from Appendix E, one finds the tangent vectors

(6.54) \begin{equation} \begin{array}{l} {\boldsymbol{a}}_1 = r\,\cos \theta \,{\boldsymbol{e}}_\phi ,\quad {\boldsymbol{a}}_2 = r\,{\boldsymbol{e}}_\theta , \end{array} \end{equation}

surface metric

(6.55) \begin{equation} \begin{array}{l} \left [a_{\alpha \gamma }\right] = r^2\left [\begin{array}{cc} \cos ^2\theta & 0 \\[4pt] 0 & 1 \end{array}\right]\kern-2.5pt, \end{array} \end{equation}

dual tangent vectors

(6.56) \begin{equation} \begin{array}{l} {\boldsymbol{a}}^1 = \displaystyle \frac {1}{r\,\cos \theta }\,{\boldsymbol{e}}_\phi ,\quad {\boldsymbol{a}}^2 = \frac {1}{r}\,{\boldsymbol{e}}_\theta , \end{array} \end{equation}

surface normal

(6.57) \begin{equation} \begin{array}{l} {\boldsymbol{n}} = {\boldsymbol{e}}_r, \end{array} \end{equation}

and curvature components

(6.58) \begin{equation} \begin{array}{l} \left [b_{\alpha \gamma }\right] = -r\left [\begin{array}{cc} \cos ^2\theta & 0 \\[4pt] 0 & 1 \end{array}\right]\kern-2.5pt, \end{array}\end{equation}

in the ALE frame. Consequently, $H = -1/r$ and $\kappa = 1/r^2$ . Further, (6.51) implies ${\boldsymbol{v}}_{{m}} = \textbf{0}$ .

6.3.2. Surface flow

Now consider the shear flow defined by

(6.59a) \begin{align} \zeta ^1&= \theta ^1 = \xi ^1 + \omega _0 t \sin \xi ^2 ,\end{align}

(6.59b) \begin{align}\zeta ^2 & =\theta ^2 = \xi ^2 =: \theta,\end{align}

which leads to

(6.60a) \begin{align}\dot \zeta ^1& = \omega _0 \sin \theta,\end{align}

(6.60b) \begin{align}\dot \zeta ^2&=0,\end{align}

(6.61) \begin{align} \left [\displaystyle \frac {\partial {\zeta ^\alpha }}{\partial {\xi ^\gamma }}\right] & = \left [\begin{array}{cc} 1 & \omega t\cos \theta \\[4pt] 0 & 1 \end{array}\right]\kern-2.5pt,\end{align}

and

(6.62) \begin{equation} \begin{array}{l} \left [\dot \zeta ^\alpha _{,\gamma }\right] = \left [\begin{array}{cc} 0 & \omega _0\cos \theta \\ [6pt]0 & 0 \end{array}\right]\kern-2.5pt, \end{array} \end{equation}

for constant angular velocity $\omega _0$ . Using (E5) from Appendix E, then leads to the frame velocities

(6.63) \begin{equation} \begin{array}{l} \dot {\boldsymbol{a}}_1 = -r\omega _0\,\sin \theta \cos \theta \,{\boldsymbol{e}}_{r_0},\quad \dot {\boldsymbol{a}}_2 = -r\omega _0\sin ^2\theta \,{\boldsymbol{e}}_\phi , \end{array} \end{equation}

the material velocity

(6.64) \begin{equation} \begin{array}{l} {\boldsymbol{v}} = r\omega _0\sin \theta \cos \theta \,{\boldsymbol{e}}_\phi \end{array} \end{equation}

and the material acceleration

(6.65) \begin{equation} \begin{array}{l} \dot {{\boldsymbol{v}}} = -r\omega _0^2\sin ^2\theta \cos \theta \,{\boldsymbol{e}}_{r_0}. \end{array} \end{equation}

From $\boldsymbol{v}$ and (E4) follow

(6.66a) \begin{align} {\boldsymbol{v}}_{\!,1}&=-r\omega _0\sin \theta \cos \theta \,{\boldsymbol{e}}_{r_0},\end{align}

(6.66b) \begin{align} {\boldsymbol{v}}_{\!,2}&= r\omega _0\,(2\cos ^2\theta -1)\,{\boldsymbol{e}}_\phi,\end{align}

which can also be obtained equivalently from (2.27). From the preceding equations follows, that the ALE equation (2.26) decomposes into the separate equations

(6.67) \begin{equation} \begin{array}{l} \dot {{\boldsymbol{v}}} = {\boldsymbol{v}}_{\!,\alpha }\,v^\alpha ,\quad {\boldsymbol{v}}^{\prime} = {\boldsymbol{v}}_{\!,\alpha }\,v^\alpha _{m}, \end{array}\end{equation}

for this example. The velocity field and its non-zero component $v_\phi:= {\boldsymbol{v}}\boldsymbol{\cdot} {\boldsymbol{e}}_\phi$ are visualised in figure 8 and compared with the numerical solution following in § 6.3.4. The velocity field leads to the surface vorticity

(6.68) \begin{equation} \begin{array}{l} \omega = \omega _0\,(2\sin ^2\theta -\cos ^2\theta), \end{array}\end{equation}

according to definition (2.40) and relation (E3).

Figure 8.

Shear flow on a sphere: (a) analytical and (b) numerical flow field $\boldsymbol{v}$ of a fixed surface, and (c) numerical flow field of a freely evolving surface; (d) analytical velocity profile $v_\phi (\theta)$ ; (e) FE error for fluid velocity $\boldsymbol{v}$ , surface tension $q$ , mesh velocity ${\boldsymbol{v}}_{{m}}$ and surface position $\boldsymbol{x}$ , converging at optimal rates: $O(n^{-1.5}_{{el}})=O(h^3)$ for fixed surface $\boldsymbol{v}$ and free surface $\boldsymbol{x}$ , and $O(n^{-1}_{{el}})=O(h^2)$ otherwise.

6.3.3. Surface forces

From (2.34), (6.56), (6.66) and (E3) now follows the rate-of-deformation tensor

(6.69) \begin{equation} \begin{array}{l} 2\boldsymbol{d} = \omega _0\cos ^2\theta \,\left ({\boldsymbol{e}}_\phi \otimes {\boldsymbol{e}}_\theta + {\boldsymbol{e}}_\theta \otimes {\boldsymbol{e}}_\phi \right) -\omega _0\sin \theta \cos \theta \,\left ({\boldsymbol{e}}_\phi \otimes {\boldsymbol{e}}_r + {\boldsymbol{e}}_r\otimes {\boldsymbol{e}}_\phi \right)\kern-2.5pt, \end{array} \end{equation}

which has the in-plane part

(6.70) \begin{equation} \begin{array}{l} 2\boldsymbol{d}_{{s}} = \omega _0\cos ^2\theta \,\left ({\boldsymbol{e}}_\phi \otimes {\boldsymbol{e}}_\theta + {\boldsymbol{e}}_\theta \otimes {\boldsymbol{e}}_\phi \right)\kern-2.5pt, \end{array}\end{equation}

and results in the stress

(6.71) \begin{equation} \begin{array}{l} {\boldsymbol \sigma } = \gamma \,\boldsymbol{i} + \eta \,\omega _0\cos ^2\theta \,\left ({\boldsymbol{e}}_\phi \otimes {\boldsymbol{e}}_\theta + {\boldsymbol{e}}_\theta \otimes {\boldsymbol{e}}_\phi \right)\kern-2.5pt, \end{array} \end{equation}

according to constitutive model (4.8) (that is valid for both the material models in § 4). From (4.9) then follows

(6.72) \begin{equation} \begin{array}{l} \boldsymbol{T}^\alpha _{;\alpha } = \left (\gamma _{,1} - 4\eta \,\omega _0\sin \theta \cos ^2\theta \right)\,{\boldsymbol{a}}^1 + \gamma _{,2}\,{\boldsymbol{a}}^2 - \displaystyle \frac {2\gamma }{r}\,{\boldsymbol{n}}, \end{array} \end{equation}

since $\textrm{div}_{s} \,\boldsymbol{d}_{{s}} = \boldsymbol{d}_{{s},\alpha }\,{\boldsymbol{a}}^\alpha$ with

(6.73) \begin{equation} \begin{array}{l} \boldsymbol{d}_{{s},1}\,{\boldsymbol{a}}^1 = \boldsymbol{d}_{{s},2}\,{\boldsymbol{a}}^2 = -\omega _0\,\sin \theta \cos ^2\theta \,{\boldsymbol{a}}^1, \end{array}\end{equation}

in this example. Due to rotational symmetry $\gamma _{,1} = 0$ . The equation of motion (3.4) is then satisfied for the body force components

(6.74a) \begin{align}f_1& = 4\eta \,\omega _0\sin \theta \cos ^2\theta,\end{align}

(6.74b) \begin{align}f_2 & = \rho \,r^2\,\omega _0^2\,\sin ^3\theta \cos \theta - \gamma _{,2},\end{align}

(6.74c) \begin{align} p & = \displaystyle \frac {2\gamma }{r} - \rho \,r\,\omega _0^2\,\sin ^2\theta \cos ^2\theta,\end{align}

that admit the two special cases:

Case 1: $f_2 = 0$ . Then (6.74( b )) can be integrated to yield

(6.75) \begin{equation} \begin{array}{l} \gamma = \displaystyle \frac {\rho \,r^2\,\omega _0^2}{4}\sin ^4\theta + \gamma _0. \end{array} \end{equation}

In this case, we can also write $\boldsymbol{f} = \eta _{{s}}\,{\boldsymbol{v}} + p\,{\boldsymbol{n}}$ , where $\eta _{{s}}:= 4\eta /r^2$ can be interpreted as a viscous surface friction coefficient.

Case 2: $\gamma = \gamma _0 =$ const. Thus $\gamma _{,2} = 0$ and we can also write $\boldsymbol{f} = \eta _{{s}}\,{\boldsymbol{v}} + \rho \,\dot {{\boldsymbol{v}}} + p_0\,{\boldsymbol{n}}$ , where $p_0 = 2\gamma _0/r$ .

Figure 9 shows the stress components $\sigma _{\phi \theta } = {\boldsymbol{e}}_\phi \boldsymbol{\cdot} {\boldsymbol \sigma }{\boldsymbol{e}}_\theta$ and $\gamma$ , as well as the surface pressure $p$ for case 1 considering $\gamma _0 = \rho r^2\omega _0^2/4$ . Equation (6.74) and its two cases provide two analytical solutions for the steady flow in (6.64). It is emphasised that this solution contains out-of-plane body forces. Only the precise pressure distribution of (6.74c) keeps the surface in spherical shape. Other distributions will lead to shape changes, as the following subsection shows.

Figure 9.

Shear flow on a rigid sphere: (a) shear stress $\sigma _{\phi \theta }(\theta)$ and surface tension $\gamma (\theta)$ , and (b) surface pressure $p(\theta)$ for case 1 considering $\gamma _0 = \rho r^2\omega _0^2/4$ .

6.3.5. Lagrangian parametrisation

Alternatively (but not advantageously), the system can be described in the material frame. Using (E6), one finds

(6.76) \begin{equation} \begin{array}{l} {\boldsymbol{a}}_{\hat 1} = {\boldsymbol{a}}_1,\quad {\boldsymbol{a}}_{\hat 2} = \omega _0 t\cos \theta \,{\boldsymbol{a}}_1 + {\boldsymbol{a}}_2, \end{array}\end{equation}

and

(6.77) \begin{equation} \begin{array}{l} \left [a_{\hat \alpha \hat \gamma }\right] = r^2\left [\begin{array}{cc} \cos ^2\theta & \omega _0 t\cos ^3\theta \\[4pt] \omega _0 t\cos ^3\theta & 1 + \omega _0^2 t^2 \cos ^4\theta \end{array}\right]\kern-2.5pt. \end{array} \end{equation}

From this follows

(6.78) \begin{equation} \begin{array}{l} \left [\dot a_{\hat \alpha \hat \gamma }\right] = r^2\left [\begin{array}{cc} 0 & \omega _0\cos ^3\theta \\[4pt] \omega _0\cos ^3\theta & 2\omega _0^2 t \cos ^4\theta \end{array}\right]\kern-2.5pt, \end{array}\end{equation}

and

(6.79) \begin{equation} \begin{array}{l} \left [a^{\hat \alpha \hat \gamma }\right] = \displaystyle \frac {1}{r^2}\left [\begin{array}{cc} \cos ^{-2}\theta + \omega _0^2 t^2 \cos ^2\theta & -\omega _0 t\cos \theta \\[4pt] -\omega _0 t\cos \theta & 1 \end{array}\right]\kern-2.5pt, \end{array} \end{equation}

so that

(6.80) \begin{equation} \begin{array}{l} {\boldsymbol{a}}^{\hat 1} = \displaystyle \frac {1}{r\cos \theta }{\boldsymbol{e}}_\phi - \frac {\omega _0 t \cos \theta }{r}\,{\boldsymbol{e}}_\theta ,\quad {\boldsymbol{a}}^{\hat 2} = {\boldsymbol{a}}^2. \end{array} \end{equation}

From (6.76) and (6.63) then follow

(6.81a) \begin{align} \dot {\boldsymbol{a}}_{\hat 1} & = -r\omega _0\sin \theta \cos \theta \,{\boldsymbol{e}}_{r_0},\end{align}

(6.81b) \begin{align} \dot {\boldsymbol{a}}_{\hat 2} &= r\omega _0 \left (\cos ^2\theta -\sin ^2\theta \right)\,{\boldsymbol{e}}_\phi -r\omega _0^2t\sin \theta \cos ^2\theta \,{\boldsymbol{e}}_{r_0}.\end{align}

Together with (2.15), (2.29) and (6.80) this leads to the same $\boldsymbol{d}$ as in (6.69), as it is supposed to. This shows that for shear flows the Lagrangian parameterisation, although physically equivalent, is mathematically much more complicating than the ALE one.

6.4. Spinning sphere

The previous example can be easily modified to describe the rigid body rotation of a spinning sphere. In that case, (6.59a ) needs to be simply replaced by

(6.82) \begin{equation} \begin{array}{l} \zeta ^1 = \theta ^1 = \xi ^1 + \omega _0 t, \end{array} \end{equation}

leading to $\dot \zeta ^1 = \omega _0$ , $\dot \zeta ^2 = 0$ , $\zeta ^\alpha _{,\hat \gamma } = \delta ^\alpha _{\hat \gamma }$ and $\dot \zeta ^\alpha _{,\gamma }=0$ . The ALE description in (6.51)–(6.58) remains unchanged. But now velocity and acceleration become

(6.83) \begin{equation} \begin{array}{lll} {\boldsymbol{v}} \! \! \! &=& \! \! \! r\omega _0\cos \theta \,{\boldsymbol{e}}_\phi \\[1mm] \dot {{\boldsymbol{v}}} \! \! \! &=& \! \! \! -r\omega _0^2\cos \theta \,{\boldsymbol{e}}_{r_0}, \end{array} \end{equation}

and lead to the gradients

(6.84a) \begin{align} {\boldsymbol{v}}_{\!,1}& = \dot {\boldsymbol{a}}_1 = -r\omega \cos \theta \,{\boldsymbol{e}}_{r_0}, \end{align}

(6.84b) \begin{align} {\boldsymbol{v}}_{\!,2} & = \dot {\boldsymbol{a}}_2 = -r\omega \sin \theta \,{\boldsymbol{e}}_\phi.\end{align}

As a consequence, the vorticity is

(6.85) \begin{equation} \begin{array}{l} \omega = 2\omega _0\sin \theta , \end{array} \end{equation}

which is maximum at the poles, while the rate-of-deformation tensor becomes

(6.86) \begin{equation} \begin{array}{l} 2\boldsymbol{d} = -\omega _0\cos \theta \,\left ({\boldsymbol{e}}_\phi \otimes {\boldsymbol{e}}_r + {\boldsymbol{e}}_r\otimes {\boldsymbol{e}}_\phi \right)\kern-2.5pt, \end{array} \end{equation}

which has zero in-plane part ( $\boldsymbol{d}_{{s}}=\textbf{0}$ ). Therefore the stress only consists of surface tension, i.e.

(6.87) \begin{equation} \begin{array}{l} {\boldsymbol \sigma } = \gamma \,\boldsymbol{i}. \end{array} \end{equation}

The body force required to equilibrate this stress under dynamic conditions now becomes

(6.88a) \begin{align}f_1 & = 0,\end{align}

(6.88b) \begin{align}f_2 & = \rho \,r^2\,\omega ^2_0\,\sin \theta \cos \theta - \gamma _{,2},\end{align}

(6.88c) \begin{align} p& = \displaystyle \frac {2\gamma }{r} - \rho \,r\,\omega ^2_0\,\cos ^2\theta,\end{align}

since still $\gamma _{,1} = 0$ due to rotational symmetry. The special case $f_2 = 0$ can then be maintained for the surface tension

(6.89) \begin{equation} \begin{array}{l} \gamma = \displaystyle \frac {\rho \,r^2\,\omega _0^2}{2}\sin ^2\theta + \gamma _0. \end{array} \end{equation}

The flow field, surface tension and surface pressure are visualised in figure 10 for the choice $\gamma _0 = \rho r^2\omega _0^2/4$ . This shows that even a simple spinning sphere requires varying body forces and surface tension to stay in spherical shape.

Figure 10.

Spinning sphere: (a) flow field $\boldsymbol{v}$ and (b) profiles of surface velocity $v_\phi (\theta)$ , surface tension $\gamma (\theta)$ and surface pressure $p(\theta)$ considering $\gamma _0 = \rho r^2\omega _0^2/4$ .

6.5. Octahedral vortex flow on a sphere

The shear flow solution of § 6.3 has no $\phi$ -dependency and is steady for fixed surface and mesh. The following example therefore presents a transient example with $\phi$ -dependency. It considers eight counter-rotating vortices on a sphere arranged in the pattern of an octahedron. The example is motivated by the surface flow induced by large-scale surface budding on spherical vesicles (Sauer et al. Reference Sauer, Duong, Mandadapu and Steigmann2017). It is adapted from a comparable example of Nitschke et al. (Reference Nitschke, Voigt and Wensch2012). The spherical surface is described by the same parameterisation from before, see (6.51)–(6.58). As before, $\phi = \zeta ^1$ and $\theta = \zeta ^2$ are used. The flow in this example is derived from a streamfunction and thus automatically satisfies the area-incompressibility constraint (2.50). The chosen streamfunction is

(6.90) \begin{equation} \begin{array}{l} \psi = v_0\,r\,\sin 2\phi \,\sin \theta \,\cos ^2\!\theta . \end{array} \end{equation}

The (tangential) surface velocity then follows from the surface curl

(6.91) \begin{equation} \begin{array}{l} {\boldsymbol{v}} = \textrm{curl}_{{s}}\psi, \end{array}\end{equation}

defined in (2.43). This gives the components

(6.92a) \begin{align} v^1 & = \dot \zeta ^1 = \displaystyle \frac {v_0}{r}\sin 2\phi \, (2\sin ^2\!\theta - \cos ^2\!\theta ),\end{align}

(6.92b) \begin{align} v^2 & = \dot \zeta ^2 = \displaystyle \frac {v_0}{r}\cos 2\phi \,\sin 2\theta,\end{align}

in the $\{{\boldsymbol{a}}_\alpha \}$ basis, i.e. ${\boldsymbol{v}} = v^\alpha \,{\boldsymbol{a}}_\alpha$ . They can be easily transformed into the components

(6.93a) \begin{align} v_\phi&= v_0\,\sin 2\phi (2\sin ^2\!\theta - \cos ^2\!\theta)\cos \theta,\end{align}

(6.93b) \begin{align} v_\theta &= v_0\,\cos 2\phi \,\sin 2\theta,\end{align}

of the spherical basis $\{{\boldsymbol{e}}_\phi,\,{\boldsymbol{e}}_\theta \}$ using (6.54), i.e. ${\boldsymbol{v}} = v_\phi \,{\boldsymbol{e}}_\phi + v_\theta \,{\boldsymbol{e}}_\theta$ . Figure 11 shows the magnitude and components of the flow field as well as the streamfunction $\psi$ and surface vorticity $\omega$ . The latter follows from (2.40), which in view of (6.91) is equal to the surface Laplacian of $\psi$ ,

(6.94) \begin{equation} \begin{array}{l} \omega = \varDelta _{{s}}\psi = a^{\alpha \beta }\,\psi _{;\alpha \beta }, \end{array} \end{equation}

and turns out to be equal to $\omega = -12\psi /r^2$ in the present example. It is noted that (6.92) is a system of ODEs from which one can reconstruct $\zeta ^\alpha = \zeta ^\alpha (\xi ^\beta,t)$ if desired.

Figure 11.

Octahedral vortex flow on a sphere: (a) analytical flow field $\boldsymbol{v}$ and (b) vorticity $\omega$ and streamfunction $\psi$ normalised by their maximum values $v_0$ , $\psi _{{max}}$ and $\omega _{{max}}$ . The latter two occur at $\phi =\pi /4$ and $\theta =\arctan (\sqrt {2}/2)$ and are equal to $\psi _{{max}} = 2v_0 r/(3\sqrt {3})$ and $\omega _{{max}} = 8v_0/(\sqrt {3}r)$ . The analytical flow field in (a) is compared with the numerical flow field for a fixed (c) and free surface (d). The radial surface displacement in (d) ranges between −1.18 % and 1.13 %. It is scaled up by a factor of 10 to increase its visibility.

From

(6.95a) \begin{align} {\boldsymbol{v}}_{\!,1}& = -c_\theta \big[2v_0\,c_{2\phi }\,c_{2\theta }\,{\boldsymbol{e}}_\phi + v_0\,s_{2\phi }\big(5c^2_\theta +2s^2_\theta \big)s_\theta \,{\boldsymbol{e}}_\theta + v_\phi \,{\boldsymbol{e}}_r \big],\end{align}

(6.95b) \begin{align}{\boldsymbol{v}}_{\!,2} & = v_0\,s_{2\phi }\big(7c^2_\theta -2s^2_\theta \big)s_\theta \,{\boldsymbol{e}}_\phi + 2v_0\,c_{2\phi }\,c_{2\theta }\,{\boldsymbol{e}}_\theta - v_\theta \,{\boldsymbol{e}}_r,\end{align}

one can find the in-plane rate-of-deformation tensor

(6.96) \begin{equation}\begin{array}{l} 2\boldsymbol{d}_{{s}} = \displaystyle \frac {v_0}{r}\left [4c_{2\phi }\,c_{2\theta }\,\left ({\boldsymbol{e}}_\theta \otimes {\boldsymbol{e}}_\theta - {\boldsymbol{e}}_\phi \otimes {\boldsymbol{e}}_\phi \right) + s_{2\phi }\,(3c_{2\theta }-1)\,s_\theta \,\left ({\boldsymbol{e}}_\phi \otimes {\boldsymbol{e}}_\theta +{\boldsymbol{e}}_\theta \otimes {\boldsymbol{e}}_\phi \right) \right]\kern-2.5pt. \end{array} \end{equation}

Here, the abbreviations $s_{\ldots }:= \sin \ldots$ and $c_{\ldots }:= \cos {\ldots }$ have been used. This leads to the stress components

(6.97a) \begin{align} \sigma _{\phi \phi }&= \gamma - 4\displaystyle \frac {v_0\eta }{r}\,\cos 2\phi \,\cos 2\theta ,\end{align}

(6.97b) \begin{align} \sigma _{\theta \theta } &=\gamma + 4\displaystyle \frac {v_0\eta }{r}\,\cos 2\phi \,\cos 2\theta , \end{align}

(6.97c) \begin{align} \sigma _{\phi \theta } &=\displaystyle \frac {v_0\eta }{r}\,\sin 2\phi \,(3\cos 2\theta -1)\sin \theta , \end{align}

according to (4.8). A meaningful quantity to examine is the surface shear stress (Zimmermann et al. Reference Zimmermann, Toshniwal, Landis, Hughes, Mandadapu and Sauer2019)

(6.98) \begin{equation} \begin{array}{l} s = \sqrt {\sigma ^{\alpha \beta }_{{\textit{de}{\kern-0.5pt}v}}\,\sigma _{\alpha \beta }^{\textit{de}{\kern-0.5pt}v}},\quad \sigma ^{\alpha \beta }_{\textit{de}{\kern-0.5pt}v}:= \sigma ^{\alpha \beta } - \gamma \,a^{\alpha \beta }, \end{array} \end{equation}

which becomes

(6.99) \begin{equation} \begin{array}{l} s = \sqrt {\sigma _{\phi \phi }^{\textit{de}{\kern-0.5pt}v}\sigma _{\phi \phi }^{\textit{de}{\kern-0.5pt}v} + \sigma _{\theta \theta }^{\textit{de}{\kern-0.5pt}v}\sigma _{\theta \theta }^{\textit{de}{\kern-0.5pt}v} + 2\sigma _{\phi \theta }^{\textit{de}{\kern-0.5pt}v}\sigma _{\phi \theta }^{\textit{de}{\kern-0.5pt}v}}, \end{array}\end{equation}

in spherical coordinates and is shown in figure 12. It confirms the octahedral symmetry of the present flow example.

Figure 12.

Octahedral vortex flow on a sphere: (a) shear stress field $s$ and (b) stress components $\sigma _{\phi \phi }^{\textit{de}{\kern-0.5pt}v}$ , $\sigma _{\theta \theta }^{\textit{de}{\kern-0.5pt}v}$ , $\sigma _{\phi \theta }$ and $s$ vs. $\phi$ for $\theta =30^\circ$ . The maximum of $s$ is $s_{{max}} = \sqrt {32}v_0\eta /r$ .

From a tedious but straightforward calculation one can find the simple relation

(6.100) \begin{equation} \begin{array}{l} \boldsymbol{d}_{{s},\alpha }\,{\boldsymbol{a}}^\alpha = - \displaystyle \frac {5}{v_0}{\boldsymbol{v}}, \end{array} \end{equation}

such that

(6.101) \begin{equation} \begin{array}{l} \boldsymbol{T}^\alpha _{;\alpha } = \gamma _{,\alpha }\,{\boldsymbol{a}}^\alpha + 2H\gamma \,{\boldsymbol{n}} - \displaystyle \frac {10\eta }{v_0}{\boldsymbol{v}}, \end{array} \end{equation}

according to (4.9). Choosing $\gamma =$ const., this satisfies PDE (3.4) for the manufactured body force

(6.102) \begin{equation} \begin{array}{l} \boldsymbol{f} = \rho \,\dot {{\boldsymbol{v}}} + \eta _{{s}}\,{\boldsymbol{v}} + p_0\,{\boldsymbol{n}} , \end{array} \end{equation}

with the viscous friction coefficient $\eta _{{s}}:=10\eta /r^2$ and static surface pressure $p_0:= 2\gamma /r$ . Taking the material time derivative of ${\boldsymbol{v}} = v_\phi \,{\boldsymbol{e}}_\phi + v_\theta \,{\boldsymbol{e}}_\theta$ , and providing $\dot v_\phi$ , $\dot v_\theta$ , $\dot {\boldsymbol{e}}_\phi$ and $\dot {\boldsymbol{e}}_\theta$ through (6.93) and (E5), leads to the material acceleration

(6.103) \begin{equation} \begin{array}{l} \dot {{\boldsymbol{v}}} = {\boldsymbol{a}}_{{s}} + a_{{n}}\,{\boldsymbol{n}}, \end{array} \end{equation}

with the components

(6.104a) \begin{align} {\boldsymbol{a}}_{{s}}&= a_\phi \,{\boldsymbol{e}}_\phi + a_\theta \,{\boldsymbol{e}}_\phi ,\quad a_{{n}} = -\displaystyle \frac {\|{\boldsymbol{v}}\|^2}{r},\end{align}

(6.104b) \begin{align} a_\phi&=\displaystyle \frac {\dot v_0}{v_0}v_\phi - 2v_0\,c_{2\phi }\,c_{2\theta }\,c_\theta \,\dot \phi + v_0\,s_{2\phi }\,\left (7c_\theta ^2-2s_\theta ^2\right)\,s_\theta \,\dot \theta ,\end{align}

(6.104c) \begin{align} a_\theta&= \displaystyle \frac {\dot v_0}{v_0}v_\theta + 2v_0\,c_{2\phi }\,c_{2\theta }\,\dot \theta - v_0\,s_{2\phi }\,\left (5c_\theta ^2+2s_\theta ^2\right)\,s_\theta \,c_\theta \,\dot \phi . \end{align}

Comparing (6.103) with (2.26) and using (6.92) and (6.95), reveals the transient part

(6.105) \begin{equation} \begin{array}{l} {\boldsymbol{v}}^{\prime} = \displaystyle \frac {\dot v_0}{v_0}{\boldsymbol{v}} + {\boldsymbol{v}}_{,\alpha }\,v^\alpha _{m} . \end{array} \end{equation}

The acceleration $\dot {{\boldsymbol{v}}}$ thus contains a transient part that is proportional to $\dot v_0$ and is only in-plane, while all the remaining terms are steady-state terms and are proportional to  $v_0^2$ . The acceleration becomes negligible for small $\dot v_0$ and small $v_0$ . Figure 13 shows the steady-state part of the in-plane and out-of-plane acceleration components. Note that $a_{{n}}$ is negative, i.e. the normal acceleration $a_{{n}}\,{\boldsymbol{n}}$ points inward. It has to be equilibrated by the varying surface pressure $p = p_0 + \rho \,a_{{n}}$ for the surface to remain spherical. If the pressure is taken constant, the shape changes as the numerical solution in figure 11(d) shows for $p = 3p_0$ . As seen, the surface bulges out, where the velocity is high, and in, where the velocity is small. Further, the outward bulges are shifted in flow direction, which agrees with intuition. This solution has been obtained by a FE computation using the discretisation described in Appendix F. It has been confirmed that the FE convergence rates of figure 8(e) are maintained here, both for the fixed surface solution in figure 11(c) and the free surface solution in figure 11(d).

Figure 13.

Octahedral vortex flow on a sphere: (a) normal component $a_{{n}}$ and (b) magnitude of the tangential component $a_{{s}}:=\|{\boldsymbol{a}}_{{s}}\|$ of the steady-state part of the material acceleration.

It is finally noted that the FE code used here has been debugged and verified using the analytical examples proposed above. They are characterised by increasing complexity and thus suitable for code development. They test all flow components, in-plane as well as out-of-plane. The last three examples are derived with zero mesh motion, but a mesh motion can be easily superimposed for further testing. This will be discussed in future work.

6.6. Vesicle budding

As a final numerical example, the budding of a spherical vesicle is shown. The example uses the Helfrich bending model of § 4.2 and the hemispherical FE set-up from Sauer et al. (Reference Sauer, Duong, Mandadapu and Steigmann2017). The initial vesicle radius is taken as $R=1\,\unicode{x03BC}$ m, its bending stiffness and viscosity are taken from Omar et al. (Reference Omar, Sahu, Sauer and Mandadapu2020), i.e. $k = 0.1236$ nNnm, $k_{g} = -0.83\,k$ and $\eta = 10^{-8}\,$ Ns m⁻¹, and the mesh stiffness is chosen as $\mu _{{m}} = 10^{-4}\,kR^2$ . As inertia is expected to play no role here, steady-state Stokes flow is considered. Budding is induced by prescribing the spontaneous curvature $H_0$ at a constant rate $\dot H_0=-0.05\,$ nm s⁻¹ within a spherical region centred around the vesicle tip. (The region has radius $r=0.4R$ ; full $H_0$ is applied within $2r/3$ and ramped down to zero between $2r/3$ and $r$ ). Contrary to Sauer et al. (Reference Sauer, Duong, Mandadapu and Steigmann2017), who consider this region Lagrangian (i.e. fixed in the reference configuration), it is treated Eulerian here (i.e. fixed in the current configuration). Figure 14 shows the bud formation over time together with the fluid velocity and surface tension $q$ .

Figure 14.

Vesicle budding: accurate reference solution for the fluid velocity $\|{\boldsymbol{v}}\|$ (a) and surface tension $q$ (b) at time $t=[40,\,80,\,120,\,160,\,200]$ ms (left to right), which corresponds to prescribed curvature $H_0 = -[2,\,4,\,6,\,8,\,10]/R$ . The result is obtained with 49 152 NURBS FEs.

Figure 15.

Vesicle budding: Eulerian (a), ALEu (b) and ALEd solution (c) for $n_{{el}} = 3072$ in comparison with the accurate reference solution (d), all for $H_0 = -6/R$ . Evolution of the tip velocity $\|{\boldsymbol{v}}\|$ (e) and tip surface tension $q$ (f) during budding for the different formulations and FE meshes. The dots mark the location where the relative error exceeds 2 %. As seen the two ALE formulations are much more accurate than the Eulerian one.

The local kinematics of bud formation resemble that of bubble inflation. Therefore, the Eulerian surface formulation can be expected to yield much less accurate results than an ALE surface formulation. This is confirmed by figure 15: the Eulerian formulation becomes inaccurate at the bud due to large mesh distortion. These distortions are much smaller for the proposed ALE formulation, resulting in much higher accuracy. Two cases are shown: uniform mesh elasticity (ALEu) and distributed mesh elasticity (ALEd – where $\mu _{{m}}$ is lowered to 10 % for points initially beyond 0.2R from the tip). The latter allows to reduce mesh distortions at the bud even further and hence allows to obtain highly accurate results as figure 15 shows.

The example demonstrates that the proposed ALE formulation also works well in the presence of bending elasticity and large local surface deformations, and it can be tailored to them. It is interesting to note that here the chosen Eulerian prescription of $H_0$ ensures that the bud remains axisymmetric, which is not the case for a Lagrangian prescription as is seen in Sauer et al. (Reference Sauer, Duong, Mandadapu and Steigmann2017). A more detailed study and analysis of vesicle budding, also for non-axisymmetric flows, is planned for future work.