Skeleton-stabilized IsoGeometric Analysis: High-regularity Interior-Penalty methods for incompressible viscous flow problems

Skeleton-stabilized IsoGeometric Analysis: High-regularity
Interior-Penalty methods for incompressible viscous flow problems

Tuong Hoang Clemens V. Verhoosel Ferdinando Auricchio
 E. Harald van Brummelen
Alessandro Reali Eindhoven University of Technology – Department of Mechanical Engineering,
P.O. Box 513, 5600MB Eindhoven, The Netherlands
IUSS – Istituto Universitario di Studi Superiori Pavia, 27100 Pavia, Italy University of Pavia – Department of Civil Engineering and Architecture, 27100 Pavia, Italy Technische Universität München – Institute for Advanced Study, 85748 Garching, Germany

A Skeleton-stabilized IsoGeometric Analysis (SIGA) technique is proposed for incompressible viscous flow problems with moderate Reynolds number. The proposed method allows utilizing identical finite dimensional spaces (with arbitrary B-splines/NURBS order and regularity) for the approximation of the pressure and velocity components. The key idea is to stabilize the jumps of high-order derivatives of variables over the skeleton of the mesh. For B-splines/NURBS basis functions of degree with -regularity (), only the derivative of order has to be controlled. This stabilization technique thus can be viewed as a high-regularity generalization of the (Continuous) Interior-Penalty Finite Element Method. Numerical experiments are performed for the Stokes and Navier-Stokes equations in two and three dimensions. Oscillation-free solutions and optimal convergence rates are obtained. In terms of the sparsity pattern of the algebraic system, we demonstrate that the block matrix associated with the stabilization term has a considerably smaller bandwidth when using B-splines than when using Lagrange basis functions, even in the case of -continuity. This property makes the proposed isogeometric framework practical from a computational effort point of view.

Isogeometric analysis, Skeleton-stabilized, High-regularity interior-penalty method, Stokes, Navier-Stokes, Stabilization method

1 Introduction

Isogeometric analysis (IGA) was introduced by Hughes et al. HughesCottrellBazilevs:2005 () as a novel analysis paradigm targeting better integration of Computer Aided Design (CAD) and Finite Element Analysis (FEA). The pivotal idea of IGA is that it directly inherits its basis functions from CAD modeling, where Non-uniform Rational B-splines (NURBS) are the industry standard. For analysis-suitable CAD models, geometrically exact analyses can be performed on the coarsest level of the CAD geometry. This contrasts with conventional FEA, which typically uses Lagrange polynomials as basis functions defined on a geometrically approximate mesh. An additional highly appraised property of IGA is that splines allow one to achieve higher-order continuity, in contrast to the -continuity of conventional FEA. We refer to CottrellHughesBazilevs:2009 (); Veiga2014iga () for an overview of established IGA developments.

In the context of viscous flow problems – particularly in the incompressible regime – IGA has been applied very successfully. Within the framework of inf-sup stable spaces for mixed formulations Boffi2013mixed (), a variety of compatible discretizations has been developed, most notably: Taylor-Hood elements bazilevs2006IGAtheory (); BuffaStokes2011 (); BressanStokes2013 (), Nédélec elements BuffaStokes2011 (), subgrid elements ruberg2012NS (); BressanStokes2013 (), and H(div)-conforming elements BuffaStokes2011 (); evans2013Stokes (); evans2013steadNS (); evans2013NS (). The mixed discretization approach leads to a saddle point system where the discrete velocity and pressure spaces are chosen differently in order to satisfy the discrete inf-sup condition. The advantage of this approach is that a stable discrete system is obtained straightforwardly from the continuous weak formulation (without any modifications) if the pair of discrete spaces is chosen appropriately.

In practice, employing the same discrete space for the velocity and pressure fields can provide advantages in term of implementation and computer resources. These advantages become more pronounced in multi-physics problems with many different field variables, for which the derivation of inf-sup stable discrete spaces can be non-trivial. The data structures required to represent the different spaces can make this approach impractical in terms of implementation and computational expenses. Moreover, in the context of IGA, using the same discretization space for all field variables enables direct usage of the CAD basis functions, which is highly beneficial from the vantage point of CAD/FEA integration.

Although there are merits in using the same discrete space for all field variables, without modification this generally leads to an unstable system in the Babuška-Brezzi sense. A common remedy to circumvent this issue is to use stabilization techniques. Various stabilization techniques have been studied in the IGA setting, most notably: Galerkin-least squares and Douglas-Wang stabilization bazilevs2006IGAtheory () and variational multiscale stabilization (VMS) BazilevsCaloCottrell:2007 (). The structure of these approaches is that the stabilization is based on element-by-element residuals. We note that recently a combination of VMS and compatible B-splines is studied in  vanOpstal2017 (). It is also noteworthy that for incompressible elasticity the use of inf-sup stable discretizations can be circumvented by using stream functions auricchio2007stream (), -bar method  elguedj2008 () and -bar method antolin2017 ().

In this contribution we propose a novel skeleton-based stabilization technique for isogeometric analysis of viscous flow problems, like those described by the Stokes equations and incompressible Navier-Stokes equations with moderate Reynold’s numbers. The skeleton-based stabilization allows utilizing identical finite dimensional spaces for the approximation of the pressure and velocity fields. The central idea is to supplement the variational formulation with a consistent penalization term for the jumps of high-order derivatives of the pressure across element interfaces. By taking into account the local continuity at each element interface, the stabilized formulation can be applied to B-splines/NURBS with varying regularities, including the case of multi-patch geometries.

The proposed stabilization technique only controls the -th order derivative in the case of B-splines/NURBS basis functions of degree with -regularity. Therefore it can be regarded as a generalization of the continuous interior penalty finite element method burman2006edge () where Lagrange basis functions are employed. This generalization enables the consideration of a large class of problems in isogeometric analysis for fluid flows. The present work encompasses a detailed study of the effect of the stabilization operator on the sparsity pattern of the mixed matrix – including an analysis of its complexity with respect to the B-splines/NURBS order – from which it is observed that the proposed technique optimally exploits the higher-order continuity properties of isogeometric analysis. We present a series of detailed numerical benchmark simulations to demonstrate the effectivity of the stabilization technique. In particular we show that oscillation-free solutions are attained, and the method yields optimal convergence rates under mesh refinements.

The outline of the paper is as follows. In Section 2 we recall the essential aspects of isogeometric analysis. In particular we introduce the skeleton structure and jump operators, and we discuss the local continuity properties across element interfaces. The skeleton-based isogeometric analysis technique for the Navier-Stokes equations is then introduced in Section 3. In Section 4 we discuss the matrix form and implementation aspects of the method, along with a study of the effect of the skeleton-stabilization operator on the sparsity pattern of the algebraic system. A series of numerical test cases is considered in Section 5 to demonstrate the performance of the proposed method. Conclusions are finally presented in Section 6.

2 Fundamentals of skeleton-based isogeometric analysis

Figure 1: Notations for a parameterization of a multipatch geometry

To provide a setting for the skeleton-based stabilization proposed in Section 3 and to introduce the main notational conventions, we first present multi-patch non-uniform rational B-spline (NURBS) spaces. We consider a domain (with or ) with Lipschitz boundary as exemplified in Figure 1. The domain is parameterized by a, possibly multi-patch (), non-uniform rational B-spline (NURBS) such that


where and are the patch-wise geometric maps and parameter domains, respectively, with the parametric map defined as


where and are the set of NURBS basis functions and the associated set of control points, respectively. The NURBS basis functions are constructed based on a set of non-decreasing knot vectors, , with


such that the number of basis functions per patch is , with the degree of the spline in the direction (). Note that for open B-splines the multiplicity of the first and last knot values is equal to . The regularity of the basis in the parametric directions depends on the order and the multiplicity of the knot value: for . On every patch the knot vectors partition the domain into a parametric mesh . The corresponding partitioning of the domain follows as


The superscript indicates the dependence of the partition on a mesh (resolution) parameter . We associate with the mesh the skeleton111This should not be confused with the topological skeleton concept in geometric modeling.:


Note that since the skeleton-based stabilization technique considered in this work pertains to inter-element continuity properties, the boundary faces are not incorporated in the skeleton. The skeleton (5) can be decomposed in the intra-patch skeleton, , and the inter-patch skeleton, :


It evidently follows from these definitions that and .

Continuity across a patch interface is achieved by matching the knot vectors associated with the two sides of the interface, and by making the corresponding control points on both patches coincident. In terms of the NURBS basis this is equivalent to linking the NURBS basis functions corresponding to the coincident control points. We denote the set of all basis functions over the domain – where interface functions have been linked – by . The space spanned by this basis is denoted by . Let us note that in the general case of a non-conforming muti-patch structure, multi-patch coupling techniques can be used such as the Nitsche’s method Ruess2013 (); nguyen2014 () or the isogeometric mortar method brivadis2015 ().

To define the regularity of the spline space we introduce the plane (or line in the two-dimensional case) in the parameter domain of patch which is perpendicular to the -direction, with its coordinate equal to that of the knot value (see Figure 1):


The regularity of the space across an intra-patch face can then be defined through the unique combination of the patch index , the direction , and the knot index , such that the associated parametric face resides in the plane . In combination with the -continuity condition across patch boundaries, the regularity of the faces is then given by:


For all functions the jumps of its -th normal derivatives across an interface vanish in accordance with


where the jump for some function is defined as , and the superscripts and refer to the traces of on the two opposite sides of .

From (8) it is inferred that in the interior of a patch the regularity per direction is controlled by the knot vector multiplicity, while across patch boundaries merely -continuity of the basis holds. We denote by the spline space with global isotropic degree and per skeleton face regularity in accordance with definition (8). In the special case of a global intra-patch regularity , i.e., we denote the function space by . A special case of this function space is that in which full regularity is achieved, i.e., .

3 Skeleton-stabilized Isogeometric Analysis for the Navier-Stokes equations

In this section we introduce the skeleton-penalty formulation for the Navier-Stokes equations in the context of Isogeometric Analysis. We commence with the formulation of the time-dependent Navier-Stokes equations in Section 3.1. Next, we introduce the discrete skeleton-penalty formulation in Section 3.2.

3.1 The time-dependent Navier-Stokes equations

We consider the unsteady incompressible Navier-Stokes equations on the open bounded domain (with or ). The Lipschitz boundary is split in two complementary open subsets and (such that and ) for Dirichlet and Neumann conditions, respectively. The outward-pointing unit normal vector to is denoted by . For any time instant the Navier-Stokes equations for the velocity field and pressure field read:


Here represents the kinematic viscosity, and the symmetric gradient of the velocity field is denoted by . The exogenous data and , represent the body forces and Neumann conditions, respectively. Without loss of generality we herein assume the Dirichlet data to be homogeneous. The initial conditions in (10) are denoted by .

For any vector space , we denote by a suitable linear space of -valued functions on the time interval . We consider the following weak formulation of (10):


The trilinear, bilinear, and linear forms in this formulation are defined as


where and denote the inner product in and dual product in , respectively. The function spaces in (11) are defined as

In the case of pure Dirichlet boundary conditions, i.e., if coincides with all of , the pressure is determined up to a constant. In that case, the pressure space is subject to the zero average pressure condition:


3.2 The Isogeometric Skeleton-Penalty method with identical discrete spaces of velocity and pressure

In this contribution we study the discretization of (11) by utilizing identical spline discretizations for the velocity and pressure fields. The global isotropic order of the spline space is denoted by and its regularity by (with ; see Section 2):


The semi-discretization in space of the weak form (11) then reads:


The pair of spaces in (14) does not satisfy the inf-sup condition, and hence the discretization in (15) is unstable. To stabilize the system, we propose to supplement the formulation with the skeleton-penalty term,


where is the regularity of the considered spline space at the element interface , is a global stabilization parameter, and is a length scale associated with this element interface. Here we define this length scale as


where and are two elements sharing the interface , and is the -dimensional Hausdorff measure. The stabilized semi-discrete system – to which we refer as the isogeometric skeleton-penalty formulation for the Navier-Stokes equations – then reads:

Remark 1.

The power associated with the interface length in (16) follows from scaling arguments. The global stabilization parameter depends on the utilized spline space . For a sufficiently smooth pressure solution, viz. , the stabilized formulation (18) is variationally consistent with the weak form (11).

Remark 2.

A special case, which is very common for CAD models, is that in which the highest regularity spline space, , is used within each patch of the domain, while -continuity is established between patches. The skeleton-penalty term (16) in this case reads:

Remark 3.

The formulation (18) based on the skeleton-penalty stabilization term (16) can also be applied to Lagrange bases, which is – in terms of function spaces – equivalent with the special case corresponding to regularity . In this case, only the jump of first order derivatives must be stabilized. This case is known as the continuous interior penalty finite element method burman2006edge (). For higher smoothness B-splines, , with regularity , the jump of first order derivatives vanishes, as a consequence of which the formulation in  burman2006edge () cannot be applied. Thus, formulation (16) is the high-regularity generalization of the continuous interior penalty finite element method. Note that although the formulation in  burman2006edge () is equivalent to the special case of , the use of higher-order Bézier elements instead of higher-order Lagrange elements affects the sparsity pattern (see Section 4.3).

Remark 4.

The weak formulation of the steady Stokes problem associated with (11) is given by:


Similar to formulation (18), the isogeometric skeleton-penalty formulation for the Stokes equations reads:


It is well-known that problem (20) is the first-order optimality condition for the saddle point of the Lagrangian functional (see e.g. Boffi2013mixed ())


Analogously, the stabilized discrete system (21) is related to the optimization problem for the modified Lagrangian functional




The stabilized discrete system (21) follows directly from the first-order optimality condition for this modified Lagrangian functional, and the stabilization term (16) appears as the variational derivative of (24). From (24) it is seen that the stabilization term (16) effectively leads to minimization of the jump of high-order derivatives of the pressure over the skeleton in a least-squares sense.

Remark 5.

For quasi-uniform meshes, the length scale can alternatively be defined as


or, even simpler, as


The numerical results presented in Section 5 are based on definition (26).

4 The algebraic form of Skeleton-stabilized Isogeometric Analysis

In this section we discuss various algorithmic aspects of the proposed skeleton-based stabilization framework. In Section 4.1 we briefly discuss the employed solution procedure for the unsteady Navier-Stokes equations, after which the algebraic form of the formulation is introduced in Section 4.2. The effect of the proposed stabilization term on the sparsity pattern of the system matrix is then studied in detail in Section 4.3.

4.1 The unsteady Navier-Stokes solution procedure

We employ a standard solution procedure for the unsteady Navier-Stokes equations. Crank-Nicolson time integration is considered in combination with Picard iterations for solving the nonlinear algebraic problem in each time step. The employed solution strategy is summarized in Algorithm 1. We denote the constant time step size by and the time step index by , such that . The solution at time step is denoted by , and the time-dependence of the non-autonomous linear operator is similarly indicated by a superscript: . The Picard iteration counter is denoted by , and the unresolved solution at iteration by . Note that for the sake of notational brevity we here omit the superscript from the variables.

Input: , , tol   # initial condition, time step, Picard tolerance
# Initialization at
# Time iteration (: Crank-Nicolson)
for in :
        # Picard iteration
        if else for in :
if :
Algorithm 1 Solution procedure for the unsteady Navier-Stokes equations

4.2 The algebraic form

Let and denote two sets of NURBS basis functions for the velocity and pressure fields, respectively. These basis functions span the discrete velocity and pressure spaces


and, accordingly, the approximate velocity field and pressure field can be written as


where and are vectors of degrees of freedom. The corresponding algebraic form of (18) then reads


with the matrix entries given by:


The algebraic form of the solution Algorithm 1 is presented in Algorithm 2.

Input: , , tol   # initial condition vector, time step, Picard tolerance
# Initialization at
# Time iteration (: Crank-Nicolson)
for in :
        # Picard iteration
        if else for in :
               Obtain by solving the linear system:
if :
Algorithm 2 Algebraic form of the solution procedure for the unsteady Navier-Stokes equations

We note that computation of the stabilization matrix requires a data structure related to the skeleton of the mesh . This data structure is constructed such that at each element interface , the jump of high-order derivatives of the basis functions over can be evaluated. It should be noted that this skeleton structure is compatible with the recently proposed efficient row-by-row assembly procedure for IGA calabro2017fastIGA ().

4.3 The k/-complexity of the skeleton-penalty operator on sparsity pattern

Figure 2: Sparsity pattern of the skeleton-penalty matrix, illustrated with univariate cubic spaces: spline space with full regularity (first column), reduced regularity (second column), minimal regularity (third column), and Lagrange space (last column). The top row shows the basis functions, the second row the stabilized (order ) derivatives, and the third row the matrix sparsity pattern of the skeleton-penalty matrix . The bandwidths of in the spline cases are , much smaller than in the Lagrange case .

The skeleton-based stabilization operator (16) affects the sparsity pattern of the discretized Navier-Stokes system due to the fact that the jump operators on the (higher-order) derivatives provide additional connectivity between basis functions. To illustrate this effect we consider the spline space , for which the derivative of order are stabilized. The top row of Figure 2 shows univariate cubic B-spline bases with , , -regularity, and Lagrange (from left to right). The second row plots the stabilized (order ) derivatives for each basis. The third row shows the sparsity pattern of the skeleton-penalty matrix associated with the operator .

The bandwidth222The bandwidth is defined as the smallest non-negative integer such that if . of the skeleton-penalty matrix is equal to , which ranges from for -splines (or at patch interfaces) to a maximum of for splines with full continuity (typical for intra patch interfaces). This observed decrease in bandwidth with decrease in regularity stems from the fact that the number of order derivatives of the basis functions that vanish on the interfaces increases with . This behavior contrasts with classical Lagrange basis functions, for which the bandwidth is equal to (the last column of Figure 2). The resulting increase in bandwidth of the jump stabilization matrix with increase in Lagrange basis order is an important drawback of the interior penalty method compared to element-based stabilization techniques. By construction, B-spline bases ameliorate this issue in the sense that even at full continuity the bandwidth of the skeleton-penalty matrix is considerably smaller than that of the Lagrange basis of equal order. This advantage – which extends to higher dimensions – makes the proposed stabilization technique computationally practical for a wide range of applications.

5 Numerical experiments

In this section we investigate the numerical performance of the Skeleton-stabilized IsoGeometric Analysis framework for a range of numerical test cases for viscous flow problems. These test cases focus on various aspects of the framework, most notably its accuracy and convergence under mesh refinement, its stability, and its robustness with respect to the model parameters.

5.1 Steady Stokes flow in a unit square

We consider the steady two-dimensional Stokes problem – i.e., problem (11) without time-dependent and convective terms – in the unit square domain . The body force is taken in accordance with the manufactured solution BuffaStokes2011 ():


This manufactured solution is visualized in Figure 2(a). Note that homogeneous Dirichlet boundary conditions are imposed on the complete boundary , and that a zero average pressure condition, , is imposed to establish well-posedness. We use a Lagrange multiplier approach to enforce this condition.

In Figure 3 we study the asymptotic -convergence behavior of the proposed method for B-splines of degree with the highest possible regularities, i.e. . The coarsest mesh considered consists of elements, which is uniformly refined until a mesh is obtained. The stabilization parameter is taken as , , . The solution obtained using quadratic splines with elements is shown in Figure 2(a). One can observe that both the pressure and velocity solutions are oscillation-free for all considered cases. Optimal convergence rates are obtained for both the velocity and the pressure field. For the -norm and -norm of the velocity error, Figure 2(b) and 2(c) respectively, asymptotic rates of and are obtained. For the -norm of the pressure shown in Figure 2(d) we observe asymptotic rates of approximately , which is half an order higher than those of the -norm of the velocity error. For inf-sup compatible discretization pairs where the degrees of the pressure and velocity spaces are and respectively, the rate of convergence of the -norm of the pressure error is known to be equal to that of the -norm of the velocity error. We attribute the improved rate for the pressure error using equal order spaces to the fact that compared to the compatible setting the pressure space is one order higher.

(a) Manufactured solution
(b) velocity error
(c) velocity error
(d) pressure error
Figure 3: (a) Solution for the steady Stokes problem in Section 5.1, pressure (color) and velocity (vector field). (b-d) Mesh convergence results for B-splines of order and regularity.

In Figure 4 we study the sensitivity of the computed results with respect to the Skeleton-Penalty stabilization parameter . The -convergence behavior of the solution using -continuous quadratic B-splines is studied for a wide range of stabilization parameters, viz. . We observe that the stabilization parameter does not affect the accuracy of the velocity field in the -norm and -norm, see Figure 3(a) and 3(b), respectively. This is an expected results, as the introduced Skeleton-Penalty term acts only on the pressure field. The pressure solution accuracy is affected by the selection of the stabilization parameter, see Figure 3(c). Choosing too large will lead to ill-conditioning of the system, while taking too small will lead to loss of stability. Figure 4 conveys, however, that the parameter can be selected from a wide range without a significant effect on the accuracy. For the case considered here accuracy deterioration remains very limited in the range . Moreover, for all considered cases we observe the rate of convergence to be independent of the choice of .

(a) velocity error
(b) velocity error
(c) pressure error
Figure 4: Sensitivity of the quadratic spline approximation of the Stokes problem on the unit square with respect to the stabilization parameter .

The performance of the proposed Skeleton-stabilized IsoGeometric Analysis framework is further studied based on the generalized Stokes equations with homogeneous Dirichlet boundary conditions:


This system – for which the body force is selected in accordance with the manufactured solution (31) – is characterized by the Damköhler number


where is the reaction coefficient, and is a characteristic length scale for the problem (in this case the width/height of the unit square). In Figure 5 we study the -convergence behavior of -continuous B-splines for various degrees and . To control the reaction term, we supplement the stabilization term with a contribution from to the scaling ratio, i.e.,

The stabilization parameter is now chosen equal to , , . Note that the non-reactive case of corresponding to resembles the case considered above. For all considered cases we observe the approximation of the velocity solution and pressure solution to be virtually independent of the Damköhler number.

Figure 5: -convergence behavior of -continuous B-spline spaces of degree for various Damköhler numbers.

To understand the effect of reduced regularity – which is of particular importance in the case of multi-patch models – we first study the B-spline discretization of the Stokes problem on the unit square with varying intra-patch regularities. That is, we consider the spline discretizations of order with regularity . A stabilization parameter of – which effectively decreases the penalty parameter with increasing order and regularity – was found to yield an adequate balance between accuracy and stability for the considered simulations. Derivation of a rigorous selection criterion for the penalty parameter is beyond the scope of the current work. Note that because the case of and has already been considered above, we here restrict ourselves to the spline degrees . The -convergence results are collected in Figure 6. Note that we plot the errors versus the square root of the number of degrees of freedom to enable comparison of the various approximations. We observe optimal convergence rates for both the velocity and the pressure approximation for all cases. As anticipated the accuracy per degree of freedom improves with increasing regularity. Note that in the case of – which is equivalent to the Lagrange basis – we observe similar approximation behavior as for the continuous interior-penalty method burman2006edge ().

(a) ,
(b) ,
(c) ,
Figure 6: -convergence results for the Stokes problem on a unit square using B-splines spaces of various degrees and regularities .

5.2 Steady Stokes flow in a quarter annulus ring

To demonstrate the performance of the proposed Skeleton-Penalty stabilization in the context of IsoGeometric Analysis, we consider the steady Stokes problem in the open quarter annulus domain

with inner radius and outer radius . We parametrize this domain using NURBS. Homogeneous Dirichlet boundary conditions are prescribed on the entire boundary and, accordingly, it holds that . The body force is selected in accordance with the manufactured solution Hoang2017mixed (); auricchio2007stream ()