Divergence-Conforming Isogeometric Collocation Methods for the Incompressible Navier-Stokes Equations

The de Rham complex is a cochain complex that is often used as a starting point for developing mixed finite element methods which preserve topological properties of the continuous problem and are typically more stable in practice [24]. In 3D, it is typically written as:

In the context of fluid flow, these infinite dimensional spaces can be interpreted as the spaces of scalar potential fields ($\Phi$), vector potential fields ($\boldsymbol{\Psi}$), velocity fields ($\boldsymbol{\mathcal{V}}$), and pressure fields ($\mathcal{Q}$). This complex is exact for simply connected domains, meaning that the range of each map is the same as the null space of the following map.

For completeness, we also state the rotated 2D de Rham complex:

Not only is the topological structure of the incompressible Navier-Stokes equations embedded in the de Rham complex, but stability conditions result from the complex as well. By creating approximation spaces for the unknown variables that conform to discrete analogs of Equations (21) and (26) we can generate numerical methods which inherit these properties. More concretely, for 3D problems let the space $\Phi_{h}$ contain the discrete scalar potentials, $\boldsymbol{\Psi}_{h}$ contain the discrete vector potentials as well as the discrete vorticity, $\boldsymbol{\mathcal{V}}_{h}$ contain the discrete velocity, and $Q_{h}$ contain the discrete pressure. Then if there exist projection operators $\Pi_{\Phi}:\Phi\rightarrow\Phi_{h}$, $\Pi_{\boldsymbol{\Psi}}:\boldsymbol{\Psi}\rightarrow\boldsymbol{\Psi}_{h}$, $\Pi_{\boldsymbol{\mathcal{V}}}:\boldsymbol{\mathcal{V}}\rightarrow\boldsymbol{\mathcal{V}}_{h}$, and $\Pi_{\mathcal{Q}}:\mathcal{Q}\rightarrow\mathcal{Q}_{h}$ such that the following commuting diagram holds

The same holds in 2D, where we instead let $\Psi_{h}$ define the discrete space to which the vorticity belongs (as well as the streamfunction), $\boldsymbol{\mathcal{V}}_{h}$ define the discrete velocity space, and $\mathcal{Q}_{h}$ define the discrete pressure space. The required commuting diagram in this case is

Of course, we have yet to define the specifics of how to construct discrete spaces such that these discrete complexes hold. For the purposes of this paper we will use compatible B-spline spaces, and the following section is devoted to introducing the basics of B-spline basis functions.

The construction of a B-spline basis in one dimension requires two objects: the degree of the basis (denoted $k$) and a series of numbers called the knot vector (denoted $\Xi=\{\xi_{1},...\xi_{n+k+1}\}$). The knots $\xi_{i}$ are non-decreasing and denote the locations in parametric space where the parametrization can change, similar to element boundaries in standard FEA. The number $n$ in the previous relation represents the total number of functions in the basis. The basis functions themselves are defined through the Cox-de Boor recursion: The $k=0$ basis functions are built as

In higher dimensions (two or three for the purposes of this paper), B-spline basis functions are constructed by simply taking the tensor product of one dimensional B-spline bases in each parametric direction. Note that different polynomial degrees and knot vectors can be used in each direction.

B-spline basis functions possess a number of useful properties for numerical method development. In particular, the smoothness of the global basis at the knot locations is controlled by the repetition of the knot value in $\Xi$. The basis at these locations is $C^{k-r}$, where $r$ is the multiplicity of the knot. Compared to a standard finite element basis, this basis has improved global continuity, which enables the use of collocation as more classical derivatives of the functions are well defined. Note that if the first and last entries in the knot vector are repeated $k+1$ times, the spline basis will become interpolatory at those locations. Such knot vectors are referred to as open knot vectors, and allow for easy specification of Dirichlet boundary conditions. For the results within this paper, all other entries in the knot vector have multiplicity one, meaning we are utilizing a B-spline basis with the highest possible order of continuity. Further, let us define the so-called regularity vector $\boldsymbol{\alpha}$ for a basis. The size of this vector is equal to the number of distinct knots, with entries equal to the polynomial degree of the basis minus the multiplicity of the corresponding knot. In terms of the regularity vector, the global basis functions are $C^{\alpha_{j}}$-continuous across the $\alpha_{j}$ unique knot. To simplify notation, the space of functions spanned by a 1D B-spline basis of degree $k$ and a provided knot vector is denoted as:

Using the basics of B-spine functions above allows us to develop discrete approximation spaces for the vorticity, velocity, and pressure. These results are built upon work in the area of isogeometric discrete differential forms [36, 37], which we will not fully develop here. The construction of these types of spaces in the context of Galerkin approximation for the Navier-Stokes equations can also be found in [22, 23]. For 3D problems, the B-spline spaces used to discretize our unknown fields are given by:

In practice we usually define $k_{1}=k_{2}=k_{3}$, and thus we can define the polynomial degree of the spline bases constructed in the above manner using a single number $k^{\prime}=k_{1}-1=k_{2}-1=k_{3}-1$. This indicates that the pressure space $Q_{h}$ will have degree equal to $k^{\prime}$ in each direction. Then, according to the above, each velocity component will have degree $k^{\prime}+1$ in one direction and degree $k^{\prime}$ in the other two. Similarly, the vorticity components will have degree $k^{\prime}+1$ in two directions and degree $k$ in the last.

In 2D, we define the following spline spaces:

We start by reviewing the form of the divergence-conforming isogeometric Galerkin methods which inspired our collocation schemes. Let us consider a problem with Dirichlet boundary conditions on the velocity for concreteness. Then we define the discrete test and trial function spaces for velocity as $\boldsymbol{\mathcal{V}}_{h,0}$ and $\boldsymbol{\mathcal{V}}_{h,\mathbf{g}}$, which are defined as the same $\boldsymbol{\mathcal{V}}_{h}$ from Equation (40) with either no penetration boundary conditions strongly enforced (for the test space $\boldsymbol{\mathcal{V}}_{h,0}$) or with the normal velocity prescribed as given by the boundary data $\mathbf{g}$ at specified collocation points (for the trial space $\boldsymbol{\mathcal{V}}_{h,\mathbf{g}}$). Similarly, define the test and trial space for pressure as $\mathcal{Q}_{h,0}$, where $\mathcal{Q}_{h}$ is the same space as in Equation (41) but with the added condition that the pressure must have zero integral. Then the Galerkin formulation for the velocity-pressure form would read

Note that in the above we have used a Nitsche approach to enforce the tangential boundary conditions in the momentum equations. This Galerkin formulation is valid if and only if the minimum entry in the regularity vectors satisfy $\min\{\boldsymbol{\alpha_{1}}-1\}\geq 0$ and $\min\{\boldsymbol{\alpha_{2}}-1\}\geq 0$, where $\boldsymbol{\alpha_{1}}$ and $\boldsymbol{\alpha_{2}}$ are the regularity vectors from Equations (39) - (41). We can write a similar Galerkin form of the vorticity-velocity-pressure form of the Navier-Stokes equations, which would yield

However, we wish to pursue a collocation method inspired by the divergence-free mixed finite element form, which we describe below. Collocation imposes additional regularity requirements on the spaces of unknowns, as the unknown fields and their derivatives will be evaluated at points rather than integrated over the domain. The spaces developed above are not only divergence-conforming, meaning discrete velocity approximations will be pointwise divergence-free, but are also regular enough to use in collocation provided the polynomial degree is sufficiently large and the global basis is sufficiently regular.

Similarly to the Galerkin setting, each discrete unknown in the collocation schemes is assumed to lie in the corresponding space from above: The discrete velocity $\mathbf{u}^{h}\in\boldsymbol{\mathcal{V}}_{h,\mathbf{g}}$, the discrete pressure $p^{h}\in\mathcal{Q}_{h,0}$, and, when applicable, the discrete vorticity $\omega^{h}\in\Psi_{h}$. For now we ignore any boundary conditions; these will be discussed in the following section. To generate the system of equations needed to solve for the coefficients for each basis function, we define sets of collocation points of the form $\boldsymbol{\tau}_{i}$ for $i=1,...,N$ for each of the governing equations. Note that the total number of collocation points should be equal to the total number of degrees of freedom in the discretization. The full discrete system is formed by requiring the strong form of the governing equations to hold at each of the collocation points.

We choose the Greville abscissae of a B-spline space as collocation points. In one dimension, the Greville abscissae are defined by

As each of the discrete unknowns lie in a different B-spline space, each of the governing equations will be collocated at a different set of Greville points. In particular, for both the velocity-pressure formulation and the vorticity-velocity-pressure formulation we use the Greville abscissae associated with the basis of the $x$-velocity ($S^{k_{1},k_{2}-1}_{\boldsymbol{\alpha_{1}},\boldsymbol{\alpha_{2}-1}}$) as collocation points for the $x$-momentum equation, the Greville abscissae associated with the basis of the $y$-velocity ($S^{k_{1}-1,k_{2}}_{\boldsymbol{\alpha_{1}-1},\boldsymbol{\alpha_{2}}}$) as collocation points for the $y$-momentum equation, and the Greville abscissae associated with the basis of the pressure ($S^{k_{1}-1,k_{2}-1}_{\boldsymbol{\alpha_{1}}-1,\boldsymbol{\alpha_{2}}-1}$) as collocation points for the continuity equation. The constitutive relation within the three field formulation is collocated at the Greville abscissae for the vorticity basis ($S^{k_{1},k_{2}}_{\boldsymbol{\alpha_{1}},\boldsymbol{\alpha_{2}}}$). The left of Figure 1 details an example of this construction for the velocity-pressure scheme, while the left of Figure 2 shows example grids for the vorticity-velocity-pressure scheme.

The last unspecified aspect of the method is enforcement of the Dirichlet boundary conditions. The enforcement of the normal boundary condition is done strongly and collocation points along a boundary for the velocity component orthogonal to that boundary are removed, as the boundary condition specifies the value of the solution at these points. This is shown on the right of Figure 1 and Figure 2, which depict the same scenarios as their counterparts but with normal boundary conditions enforced.

Enforcement of the tangential boundary condition is slightly more subtle. Recall that in Equation (42) we utilized Nitsche’s method to enforce this boundary condition. This motivates the enforcement in the velocity-pressure collocation scheme. Indeed if we take this equation and undo the integration by parts, the consistency term vanishes by construction and we are left with just the penalty terms. If we approximate the integral of the test function as done in [41] the collocated momentum equations will be of the form

In the vorticity-velocity-pressure scheme we do not use the same Nitsche-based approach. The method utilized here is directly inspired by the Enhanced Collocation method for enforcing Neumann boundary conditions in isogeometric collocation schemes [41].

We start by considering the weak form of the constitutive relation given by Equation (47). The final term on the left hand side is the boundary term which would be used to enforce natural boundary conditions in a Galerkin method. In a similar vein to the Enhanced Collocation approach, we can undo the integration by parts to arrive at

By approximating the integrals of the test functions as done in [41, 38] we arrive at a modified strong form statement of the constitutive relation which can be collocated along the boundaries

Finally, the results of the previous sections are collected and we present the final form of the discrete equations used to solve for the discrete unknowns. The velocity-pressure scheme is considered first. Let us define $\boldsymbol{\tau}^{u_{x}}_{\ell}$ for $\ell=1,...,M^{u_{x}}$ to be the set of Greville points for $S^{k_{1},k_{2}-1}_{\boldsymbol{\alpha_{1}},\boldsymbol{\alpha_{2}-1}}$ with the points corresponding to no-penetration boundaries removed as discussed previously. Define in a similar manner $\boldsymbol{\tau}^{u_{y}}_{\ell}$ for $\ell=1,...,M^{u_{y}}$, which are the Greville points of $S^{k_{1}-1,k_{2}}_{\boldsymbol{\alpha_{1}-1},\boldsymbol{\alpha_{2}}}$ with no-penetration boundary points removed. Lastly, $\boldsymbol{\tau}^{p}_{\ell}$ for $k=1,...,N^{p}$ are the Greville points of $\mathcal{Q}_{h}$. Then the discrete 2D problem reads:

In the above we have split the momentum equations into expressions valid on the interior collocation points (Equations (52) and (54)) and expressions valid on the remaining boundary collocation points (Equations (53) and (55)).

Similarly, for the vorticity-velocity-pressure scheme we also define $\boldsymbol{\tau}^{\omega}_{\ell}$ for $\ell=1,...,N^{\omega}$ as the Greville points of $\Psi_{h}$. With this scheme the discrete 2D problem reads:

In the three field formulation we split the constitutive law into an expression for the interior collocation points (Equation (60)) and another expression for boundary collocation points (Equation (61)).

Resulting from these equations are nonlinear systems of equations which we can use to solve for the unknown coefficients of velocity, pressure, and vorticity using a Newton-Raphson method.

From the commuting diagrams our spaces form, shown in Equations (30) and (31), it is simple to show that both of the resulting collocation methods return an exactly pointwise divergence free velocity field. The commuting diagrams reveal that the divergence of the discrete velocity lies within the discrete pressure space, $\nabla\cdot\mathbf{u}^{h}\in\mathcal{Q}_{h}$. We can thus equivalently write the divergence of the velocity as a linear combination of the pressure basis functions:

If the choice of collocation points yields a set of linearly independent equations, that is to say $\mathbf{M}$ is invertible, then we know that the solution to Equation (63) is $\mathbf{c}=\mathbf{0}$, and thus the velocity field is exactly divergence free pointwise.

As a first numerical experiment, we consider a vortex problem on the unit square constructed using the method of manufactured solutions. This solution was originally developed in [21] and employs the following velocity and pressure fields:

For the velocity-pressure scheme we define the forcing term to be

To understand the accuracy of this collocation method, we test the convergence rates on a variety of grids and with different degree B-spline bases. For this test we set the Reynolds number to $Re=\frac{1}{\nu}=1$. We measure error using the $L^{2}$ norm as well as the $H^{1}$ semi-norm. Figure 3 details the convergence rates of velocity and pressure when using the two field formulation. In this case the errors in both velocity and pressure converge at a rate of $k^{\prime}$ when $k^{\prime}$ is even and $k^{\prime}-1$ for odd $k^{\prime}$. These are the standard, expected rates that have been seen in other studies of isogeometric collocation, and are one and two orders suboptimal in $L^{2}$ for odd and even $k^{\prime}$.

Figure 4 details the convergence of velocity, kinematic pressure, and vorticity as we refine the meshes in the three field scheme. Using this scheme, all of the unknowns converge in the $L^{2}$ norm at a rate of approximately $k^{\prime}$ for even values of $k^{\prime}$ and at a rate of $k^{\prime}+1$ for odd values of $k^{\prime}$. These rates match the rates achieved using even $k^{\prime}$ in the two field formulation, and these rates are two orders faster for odd $k^{\prime}$. In fact, this formulation recovers optimal convergence rates in the $L^{2}$ norm for odd $k^{\prime}$.