A Mimetic Spectral Element Method for Free Surface Flows

We derive a compatible high order discretization method that relies heavily on the underlying geometric structure and obeys the topological sequences for 2D free surface flows with fully submerged structures. The flow is assumed inviscid and irrotational. This allows us to write the continuity and Navier–Stokes equations using a velocity potential as a Laplace equation and the differential form of the unsteady Bernoulli’s equation. We use a mixed variational formulation and results confirm the theoretical results that we obtain global and local divergence free solutions.

Using the vector identity 1 2 ∇(u · u) = (u · ∇)u + u × (∇ × u) and using that the fluid is irrotational (∇ × u = 0), we can rewrite the momentum and continuity equations, where φ is a vector potential defined as u = ∇φ, φ = φ(x, z, t). We can now rewrite the momentum equation as, which we can integrate in space to obtain the time dependent Bernoulli's equation.
where C(t) is an arbitrary function of integration. We assign C(t) = 0 by recalling that φ and φ + C(t)dt yield exactly the same flow. Redefining φ and retaining the symbol φ := φ + C(t)dt we obtain the time dependent Bernoulli's equation for the problem as, The governing equations for inviscid and irrotational flows for an incompressible fluid are stated through (2) and (3), where the unknowns are the velocity potential, φ, and the pressure, p. Equations (2) and (3) together with proper boundary conditions constitute a well-posed problem. The velocity potential, φ, can be solved from the Laplace equation and then substituted into the Bernoulli's equation to obtain the pressure field.

Boundary Conditions
The physical domain is shown in Fig. 1, where the notations are also illustrated.
The time domain is taken as T : t ≥ 0. The unknowns for the problem become the velocity potential and the free surface elevation η(x, t) : The pressure can hereafter be determined through (3). The unsteady kinematic and dynamic free surface boundary conditions are given by Zakharov [8], where˜signify functions defined only on the free surface. The vertical component of the velocityν = ∂ z φ| z=η is calculated by solving the Laplace problem (2) together with the Zakharov boundary conditions (4) and (5) on the free surface. On the bottom we have the no penetration condition, On the inlet and outlet boundaries (Γ \Γ F S ∪ Γ b ) the gradient of the velocity potential is specified. The computational domain is shown in Fig. 2.

Discretization of Governing Equations
The developed method adopts elements from differential geometry. The unknowns of our system are described by use of differential forms. In a three-dimensional setting we are making use of four types of sub-manifolds: points, curves, surfaces, and volumes, both as inner and outer oriented objects, see an example in Fig. 3. The mimetic spectral element method uses an approach similar to the Galerkin method of the finite element method where the numerical residual is weighted by an arbitrary weight function. In contrast to the traditional finite element method the arbitrary weight functions are taken from the dual space of the function space used by the unknowns.

Basis Functions
For the polynomial representation we use Lagrange polynomials l i (x) and edge polynomials e i (ξ ), see [5]. The Lagrange polynomials are based on a Gauss-Lobatto-Legendre (GLL) point distribution for the nodal values. The Lagrange polynomials and edge polynomials satisfy the properties, and the edge polynomials are explicitly given in terms of the nodal Lagrange basis functions l i (x) as where dl k (ξ ) is the exterior derivative applied to the 0-form l k (ξ ). This definition of the edge polynomial also implies, see [4] and [5],

Mimetic Discretization in 2D
If we let the 0-form φ (0) ∈ Λ 0 (M) be expanded as then we can write φ (0) h as a matrix-vector product where L i,j = l i (ξ j ) and ξ j are the Gauss-Lobatto-Legendre points, GLL points. If we let the 1-form u (1) ∈ Λ 1 (M) be defined as we can expand u ξ and u η using edge polynomials as, The discrete one-form u (1) can also be written as a matrix-vector product, where u is evaluated in the GLL points, where E i,j = e i (ξ j ). The 2-form P (2) ∈ Λ 2 (M) is expanded using only edge polynomials, The Laplace equation can be reformulated using a mixed formulation, see [1], where the equilibrium equation and the constitutive relationship are separated into two equations.
Writing (16) using differential geometry for a 3-D geometry we obtain, where we have utilized the Hodge star operator. The Hodge star operator is a map, which maps p-forms onto (n − p)-forms, where n is the dimension of the domain, Ω. Given a p-form, λ (p) , the hodge star maps as follows: where˜denotes the change of orientation of the new form. The Hodge star is also the coupling between the outer oriented domain and the inner oriented dual space, as seen in Fig. 3. In 2-D, using differential geometry, equations (16) take the form, When the exterior derivative is applied to the balance equation of (21) we obtain, see [5] dq (1) where we have utilized (8). The equilibrium equation, the first equation in (21), is equated to a zero valued 2-form. Expanding the last equation in (21) yields, where f i,j = 0. The basis can then be cancelled and we can rewrite (23) as, where E (2,1) is an incidence matrix, only consisting of 0, 1 and −1. This matrix relates the fluxes of q to the volume integral of the balance equation, see Fig. 4. The first step in developing the discrete system is the formulation of the weak form, where we make use of duality pairing between an arbitrary k-form, α (k) , and The pairing with the (n − k)-form, β (n−k) , takes the role of a weight function in traditional finite element analysis and lives in the dual space and carry the opposite orientation. The result of duality pairing can also be represented as a matrix-vector product, where W contains the Gauss weights and J is the Jacobian matrix. M (k) is a mass matrix of the corresponding discretized kand (n − k)-form pairing, and( ) denotes a matrix of opposite orientation.
Using Stokes generalized theorem [3] and applying integration by parts to the balance equation (the last equation of (21)), we obtain. Ω dq (1) Using duality pairing, an inner product projection for the term with the Hodge star operator, the expansions in (9)-(15), and appropriate boundary conditions we can set up the matrix system for the discrete Laplace operator as shown in (29). ⎡ Using the forward Euler scheme for the temporal term and pairing it with an arbitrary 0-form,α (0) , we can rewrite the Bernoulli's equation as, The density is considered a 2-form, which leaves (30) Hodge invariant. The interior product i is defined in [7]. The discrete version of (30) takes the form, is derived from the interior product of the two 1-forms in the convective term, and contains information of φ from the previous time step and consequently has to be updated at each new time step.
The simulation is initialized by first solving the Laplace equation with the prescribed boundary conditions. The initial velocity potential φ on the free surface, is set toφ(x, t = 0) = x, and the free surface height is set to η(x, t = 0) = 0. At the following time steps, the Zakharov free surface equations are solved to obtain new values ofφ as well as the free surface elevation, η.

Numerical Results
The method is first applied to a non-temporal problem without a free surface. The geometry sketched in Fig. 5 contains a cylinder in the middle of a square. On the horizontal walls of the square and the cylinder wall the no penetration condition is applied. On the left vertical boundary a fully developed velocity profile is specified and on the right vertical boundary a constant velocity potential is defined. The velocity potential φ and streamlines are shown in the middle section of Fig. 5. In the right part of Fig. 5 the pressure field is shown. Figure 6 shows that we obtain spectral convergence for both unknowns.
Furthermore, the balance equation ∇ · u = 0 (conservation of mass) is satisfied both globally and point-wise independent of polynomial order as shown in Fig. 7.
Next we apply the method to a temporal and free surface problem where we have included a bump on the bottom boundary. The Zakharov free surface equations are applied on the top horizontal boundary. In Fig. 8 the pressure field and the free surface are plotted at t = 1,100,200.

Discussion and Conclusion
Using an isoparametric, multi-element formulation the solution of the discretized Laplace equation shows spectral convergence. In addition, we observe that mass is conserved both globally and locally. In (31), the discretized Bernoulli equation was kept Hodge-invariant, leaving the equation metric free. This suggests that the fundamental invariant of the equation is conserved. The Bernoulli equation conserves the total energy of the system. However, in Fig. 9 it is observed that a small amount of energy is gained and lost in a periodic manner. It is also observed that the mean energy is constant. It was possible to time integrate over very long time periods without noticing any degradation of data and we conclude that energy is conserved over long time periods even though fluctuations were observed for short time periods. In the future we plan on using a mimetic time integration scheme, which was used in [6], as well as the mimetic spatial discretization that was used in the present work. The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.