The solution of boundary value problems on curved surfaces has many practical applications in mathematics, physics, and engineering. For example, there are transport processes on interfaces, e.g., in foams, biomembranes and bubble surfaces [16, 25, 45], or structure-related phenomena such as in membranes and shells [2, 7]. Herein, Stokes and incompressible Navier-Stokes flows on curved, two-dimensional manifolds are considered. The governing equations for flows on moving surfaces are discussed in [3, 33] based on fundamental surface continuum mechanics and conservation laws and in , an energetic approach is presented. Earlier works in a similar context may be traced back to [15, 26, 43, 47]. For an excellent overview, the reader is refered to . The references given above often focus on mathematical properties such as the existence and uniqueness of the solutions or stabilitity analyses. Applications are often two-phase flows where the fluid field in the bulk and on the moving interface are coupled. However, it is also worthwhile to consider the situation for fixed manifolds, e.g., related to meterology and oceanography where the flows take place on (part of) a sphere. Special geometries such as hyperbolic planes and spheres are discussed in [6, 34, 36].
Herein, the focus is on the approximation of stationary and instationary (Navier-)Stokes flows on fixed manifolds based on the surface finite element method as outlined in [11, 13, 14]. The governing equations resemble the three-dimensional (Navier-)Stokes equations where the classical gradient and divergence operators are replaced by their tangential counterparts derived from tangential differential calculus . The equations are formulated in the classical stress-divergence form, contrasted to the approach in . An additional constraint is required to enforce that the velocities remain in the tangent space of the manifold; it is labelled “tangential velocity constraint”. The models are first given in strong form and are then transformed to the weak form to enable a numerical solution based on the surface FEM. Finite element spaces of different orders are employed for the approximation of the geometry and of the involved physical fields, i.e., the velocities, pressure and the Lagrange multiplier field required to enforce the tangential velocity constraint. It is found that the balance of these element orders is critical for the accuracy and conditioning of the system of equations. In particular, the well-known Babuška-Brezzi condition applies [1, 4, 17] as both, the incompressibility constraint and the tangential velocity constraint are enforced using Lagrange multipliers. For the case of the instationary Navier-Stokes equations, the Crank-Nicolson time stepping scheme is employed for the semi-discrete sytem of equations resulting from using the surface FEM in space. Surface FEM based on linear elements is used in the recent work , where the penalty method is employed to enforce tangential velocities and a projection method rather than a monolithic approach is suggested to solve for the different physical fields. Alternatives for the surface FEM are the TraceFEM [8, 22, 40] and CutFEM [27, 28], where the basis functions are generated from a background mesh in the bulk surrounding the manifold of interest.
Using the FEM for the Navier-Stokes flows at large Reynolds numbers requires stabilization. Herein, the streamline-upwind Petrov-Galerkin (SUPG) approach is used [5, 49]. Alternatively, other variants such as the Galerkin least squares stabilization  and variational multiscale approaches [29, 23] may also be employed. Stabilization for advection-diffusion applications on manifolds are considered in .
The numerical results show that higher-order convergence rates are achieved provided that the finite element spaces are properly chosen. Also the conditioning of the system of equations depends on the element orders employed for the approximation of the individual physical fields. The presented results are based on well-known benchmark test cases in two dimensions such as driven cavity flows and cylinder flows with vertex shedding which, herein, are extended to curved surfaces. Due to the higher-order elements, the results are highly accurate and may serve as future benchmarks in the context of (Navier-)Stokes flows on manifolds. Most test cases are carried out on parametrized surfaces, however, also the situation of flows on zero-isosurfaces is covered herein.
To the best of our knowledge, this is the first time, where (i) general higher-order surface FEM is used for the (in)stationary (Navier-)Stokes equations on manifolds including stabilization, (ii) numerical convergence studies are presented confirming higher-order convergence rates, and (iii) benchmark test cases are proposed and solutions presented. Furthermore the notation employed is closely related to the typical engineering literature and aims to provide a bridge from the mathematical to the engineering community.
The paper is organized as follows: In Section 2, some requirements and properties of surfaces are described and tangential differential operators are defined based on [10, 14]. Section 3 covers the governing equations for (i) Stokes flow, (ii) stationary, and (iii) instationary Navier-Stokes flows on two-dimensional manifolds. They are given in strong form, weak form, and discretized weak form according to the surface FEM. Numerical results are presented in Section 4. Convergence studies are performed for a test case for which an analytic solution is available and it is shown that higher-order convergence rates can be achieved. For the other test cases where no analytic solutions are available it is confirmed that in the flat two-dimensional case, well-known reference solutions are reproduced. Various meshes with different orders and resolutions have been employed to obtain highly accurate results on curved surfaces. Finally, a summary and outlook are given in Section 5.
The task is to solve a boundary value problem (BVP) on an arbitrary surface in three dimensions. Let the surface be fixed in space over time, possibly curved, sufficiently smooth, orientable, connected (so that there is only one
surface), and feature a finite area. There is a unit normal vectoron . The surface may be compact, i.e., without a boundary, , see Figs. 1(a) and (b) for examples. Otherwise, it may be bounded by as shown in Figs. 1(c) and (d). Then, associated with , there is a tangential vector pointing in direction of and a co-normal vector pointing “outwards” and being normal to and tangent to . The surface may be given in parametrized form or implied, e.g., based on the level-set method; both situations are considered herein. For the equivalence of these two cases and more mathematical details, see, e.g., .
2.2 Surface operators
2.2.1 The tangential projector
On the manifold , the tangential projector is defined by the normal vector as
Some important properties are: (i) , (ii) , and (iii) .
2.2.2 Surface gradient of scalar quantities
The tangential gradient operator of a differentiable scalar function on the manifold is given by
where is the standard gradient operator, and is a smooth extension of in a neighborhood of the manifold . Of course, may also be some given function (rather than an arbitrary extension) in global coordinates, i.e., . For the case of parametrized surfaces defined by the map , and a given scalar function , the tangential gradient may be determined without explicitly computing an extension using
with being the ()-Jacobi matrix and
being the metric tensor (first fundamental form). Equation (2.2) shall be used later in the context of the FEM to determine tangential gradients of shape functions. It is noteworthy that is in the tangent space of and, thus, and . The components of the tangential gradient are denoted by
representing the first-order partial derivatives on . Second-order partial derivatives may be denoted by
where is the tangential Hessian matrix. In the context of manifolds, this matrix is not symmetric , i.e., for mixed second derivatives for .
2.2.3 Surface gradient of vector quantities
Next, operators for vector quantitites are considered. The “directional gradient” of is the tensor of tangential derivatives and defined as
In contrast, the covariant derivatives are
One has to carefully distinguish these two different gradient operators. It is noted that appears frequently in the modeling of physical phenomena on manifolds, i.e., in the governing equations. On the other hand, is often used in straightforward extensions of identities such as product rules and divergence theorems. For example, we have for a scalar function and vector functions ,
however, the relations are less straightforward for the covariant counterparts and , respectively. Later on, in the context of FEM implementations, it proves useful to transform covariant derivatives systematically to directional ones. This allows the computation of directional derivatives of FE shape functions with respect to independent of the integration of the weak form of the governing equations.
2.2.4 Divergence operators and divergence theorem
The divergence of a vector function is given as
For a tensor function , there holds
It may be shown that with being the mean curvature and being the second fundamental form.
3 Governing equations
In the following, we consider (i) stationary Stokes flow, (ii) stationary Navier-Stokes flow, and (iii) instationary Navier-Stokes flow on fixed manifolds. The governing equations are first given in strong and weak forms. The surface FEM is then applied for the discretization of the weak forms. As mentioned above, these models are also considered, e.g., in [3, 33, 35] among others.
3.1 Flow models in strong form
3.1.1 Stationary Stokes flow
Starting point is stationary Stokes flow on a manifold. Let be the three-dimensional velocity field on the surface , a pressure field, and a tangential body force, e.g., with unit . The governing field equations (in stress-divergence-form ) to be fulfilled are
Equation (3.1) expands to three momentum equations, equation (3.2) is the incompressibility constraint and equation (3.3) represents the tangential velocity constraint that restricts the velocities to the tangent space of . Two different strain tensors are introduced,
which are related to each other as . The stress tensor is then defined as
where is the (constant) dynamic viscosity. It is easily shown that
Suppose there exists a boundary of the manifold that consists of two non-overlapping parts, the Dirichlet boundary, , and the Neumann boundary, . The corresponding boundary conditions are given as
where the prescribed velocities and tractions are in the tangent space of , i.e.,
Note that, in general, there are no explicit boundary conditions needed for the pressure . In cases where no Neumann boundary is present, i.e., and , the pressure is defined up to a constant [12, 24]. This includes compact manifolds where . In such situations, the pressure may be prescribed at a given point on or it is imposed by a constraint in the form of .
Vorticity on manifolds.
The vorticity is a physical quantity frequently computed in flow problems. In the context of manifolds, we shall define
Note that is co-linear to the normal vector , hence, . Therefore, it is useful to determine the signed magnitude of , that is, the scalar function
This scalar quantity may also be obtained using directional derivatives, i.e., .
3.1.2 Stationary Navier-Stokes flow
For stationary Navier-Stokes flow, a non-linear advection term is added to equation (3.1) resulting into
where is the (constant) fluid density with unit and . It is quite common to express the body force in the form where may consider gravity as for instance. The remaining equations (3.2) and (3.3) and the boundary conditions (3.6) remain unchanged. The solution of the non-linear governing equations can be obtained iteratively based on the Newton-Raphson method or other fixed-point iterations such as Picard iterations. Because the advection operator is not self-adjoint, well-known stability issues may arise for large Reynolds numbers in a numerical context.
3.1.3 Instationary Navier-Stokes flow
For instationary Navier-Stokes flow, the momentum equation (3.1) changes to
The functions representing the physical fields live in space (on ) and time, i.e., in the time interval . Therefore, Eqs. (3.10), (3.2), and (3.3) have to be solved in the space-time domain . Herein, we restrict ourselves to spatially fixed manifolds .
The boundary conditions (3.6) also extend in time dimension, hence, there are prescribed velocities along and tractions along . Furthermore, an initial condition is needed,
3.2 Flow models in weak form
The following trial and test function spaces are introduced,
As mentioned previously, if no Neumann boundary exists, i.e., , the pressure is defined up to a constant and one may replace by
3.2.1 Stationary Stokes flow
The weak form of the Stokes problem becomes: Given viscosity , body force in , and traction on , find the velocity field , pressure field , and Lagrange multiplier field such that for all test functions , there holds in
The following relations are easily derived:
It is readily verified that solutions of the strong form also fulfill the weak form from above. This is obvious for Eqs. (3.17) and (3.18) due to Eqs. (3.2) and (3.3), respectively. For the momentum equations, it is noted that (3.16) is fufilled for . Restricting this to the tangential space by multiplication with the projector yields the strong form of the momentum equations (3.1) because . It is thus also seen that the Lagrange multiplier field may be physically interpreted as a force in normal direction.
3.2.2 Stationary Navier-Stokes flow
The weak form of the stationary Navier-Stokes equations is similar to the Stokes problem from above, however, Eq. (3.16) is replaced by
where the added advection term is readily identified.
3.2.3 Instationary Navier-Stokes flow
The weak form of the instationary Navier-Stokes problem is: Given density , viscosity , body force in , traction on , and initial condition on at according to (3.11), find the velocity field , pressure field , and Lagrange multiplier field such that for all test functions , there holds in
3.3 Surface FEM for flows on manifolds
3.3.1 Surface meshes
Assume that a suitable surface mesh composed by higher-order triangular or quadrilateral Lagrange elements of order may be generated with desired element sizes and all nodes on . Well-known, necessary requirements of meshes such as the shape regularity of the elements and bounds on inner angles, are fulfilled. The shape of each (physical) element in the mesh results from a map of the corresponding reference element with nodes,
are classical Lagrangean shape functions of order in reference coordinates and are the nodal coordinates. The resulting mesh is an approximation of the exact surface . Clearly, is defined parametrically through the map (3.23) even if the original was implicitly given, e.g., by the zero-isosurface of a level-set function. See [18, 19, 20] for the automatic generation of higher-order meshes on zero-isosurfaces. The discrete unit normal vector is
and is not smooth across element edges due to the -continuity of the surface mesh. The discrete tangent and co-normal vectors and are easily obtained along the element edges on the boundary of . The definitions of the surface operators from Section 2.2 readily extend to the case of a discrete manifold and are not repeated here.
3.3.2 Surface FEM
We use higher-order surface FEM as detailed, e.g., in [11, 14] for the discretization of the weak forms from above. Finite element spaces of different orders are involved. As mentioned before, suitable surface meshes of order may be generated defining approximations of Let there be a “geometry mesh” of order with the sole purpose to approximate the geometry of the manifold and define the element maps (3.23). In particular, this mesh is not used to imply a finite element space for the approximation of the weak forms.
Next, a finite element space of order is generated on for which it is assumed that there is a second mesh of order . The two meshes feature the same element types and number of elements with identical coordinates at the corners, however, the total number of nodes differs due to the individual orders. It is emphasized that the coordinates of the nodes in the -th order mesh are, in fact, never needed and it is only the connectivity which is required to set up the finite element space.
Associated to triangular or quadrilateral elements in the -th order mesh, there is a fixed set of local basis functions defined in a reference element with and being the number of nodes per element. Classical Lagrange basis functions with are used herein. Based on the map (3.23) which is completely determined by the geometry mesh, one may generate for all and tangential derivatives are determined based on Eq. (2.2). This is only an iso-parametric map when . Summing up the element contributions for nodes belonging to several elements, this generates a set of global, -continuous basis functions in with and being the number of nodes of the -th order surface mesh. Note that to generate the nodal basis , only the coordinates of the geometry mesh are needed, however, not from the -th order mesh. A general finite element space of order is now defined by
Based on this, the following discrete trial and test function spaces are defined,
Although shape functions for the pressure and the Lagrange multiplier for enforcing the tangential velocity constraint may be discontinuous, we restrict ourselves to classical -continuous approximations. Note that individual orders , , and are associated to the approximations of velocities , pressure , and Lagrange multiplier field , respectively. Analogous to the continuous case, one may impose that the functions in have to fulfill if no Neumann boundary is present.
3.3.3 Stationary Stokes flow
The discrete weak form of the Stokes problem reads: Given viscosity , body force in , and traction on , find the velocity field , pressure field , and Lagrange multiplier field such that for all test functions , there holds in
The usual element assembly yields a linear system of equations in the form
with being the sought nodal values of the velocity components, pressure, and Lagrange multiplier. For the implementation, it is interesting to compare the system (3.31) with the system obtained for a classical three-dimensional Stokes problem,
Assume a function which generates and based on three-dimensional FE shape functions (including classical partial derivatives with respect to ) evaluated at given integration points in 3D. The same function may be used for generating and provided that (i) the integration points are restricted to with proper weights, (ii) the classical partial derivatives in are replaced by the tangential derivatives as in , and (iii) the contribution to at the current integration point, , is projected as which is due to Eq. (3.19). The same shall later hold for the advection matrix in the Navier-Stokes equations.
As expected in the context of the Lagrange multiplier method, the matrix in Eq. (3.31) has a saddle-point structure and is typical for a mixed FEM. The well-known Babuška-Brezzi condition [1, 4, 17] must be fulfilled to obtain useful solutions for all involved fields. This may be achieved by adjusting the orders of the approximation spaces for the different fields and is further detailed in the numerical results. It is noted that stabilization may be employed to circumvent the Babuška-Brezzi condition rather than to fulfill it, see, e.g., [17, 30, 31] which is, however, beyond the scope of this work.
3.3.4 Stationary Navier-Stokes flow
The discrete weak form of the stationary Navier-Stokes problem reads: Given density , viscosity , body force in , and traction on , find the velocity field , pressure field , and Lagrange multiplier field such that for all test functions , there holds in
The equations related to the different field equations were added up for brevity. The last row adds a stabilization term which is needed to obtain stable solutions for flows at high Reynolds numbers [12, 24]. In particular, the streamline upwind Petrov-Galerkin (SUPG) method is used for the stabilization. Different definitions of the stabilization parameter are found [44, 49, 48] and
is used herein with element-averaged velocity , element length and for the stationary case. When stabilization is not necessary because no oscillations occur, Note that in the stabilization term, second-order derivatives appear (only in the element interiors). The definition of tangential second-order derivatives is given, e.g., in .
Element assembly results in a non-linear system of equations of the form
which is no longer symmetric (partly) due to the advection matrix . The distinguishing feature of and (compared to and of the Stokes problem) are the added SUPG-stabilization terms. The issues related to mixed FEMs and the Babuška-Brezzi condition remain relevant.
3.3.5 Instationary Navier-Stokes flow
The discrete weak form of the instationary Navier-Stokes problem is: Given density , viscosity , body force in , traction on , and initial condition on at according to (3.11), find the velocity field , pressure field , and Lagrange multiplier field such that for all test functions , there holds in
This yields a system of non-linear semi-discrete equations for
with initial condition . This system may be advanced in time by using finite difference schemes and the Crank-Nicolson method is employed herein.
4 Numerical results
The following error measures are computed in the convergence studies. When analytic (exact) velocity and pressure fields, and , are known, the velocity error is determined by
and the pressure error calculated as
When analytic solutions are not available, it is useful to evaluate the error of the FE approximations in the strong form of the momentum or continuity equations, integrated over the domain. For the example of stationary Stokes flow, the corresponding residual errors are defined as
This can be easily extended to the case of Navier-Stokes flows where the advection term is added to the integrand in (4.3). Also the error in the tangential velocity constraint from Eq. (3.3) may be computed in a similar manner. The evaluation of the error involves second-order derivatives and convergence can only be expected for higher-order elements and sufficiently smooth solutions.
4.1 Stokes flow on an axisymmetric surface
A test case is developed for which analytic solutions are available. An axisymmetric surface with height and radius
is generated as illustrated in Fig. 2(a). Let and . In parametrized form, one may also define based on the map ,