Skeletonstabilized IsoGeometric Analysis: Highregularity
InteriorPenalty methods for incompressible viscous flow problems
Abstract
A Skeletonstabilized 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 Bsplines/NURBS order and regularity) for the approximation of the pressure and velocity components. The key idea is to stabilize the jumps of highorder derivatives of variables over the skeleton of the mesh. For Bsplines/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 highregularity generalization of the (Continuous) InteriorPenalty Finite Element Method. Numerical experiments are performed for the Stokes and NavierStokes equations in two and three dimensions. Oscillationfree 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 Bsplines 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.
keywords:
Isogeometric analysis, Skeletonstabilized, Highregularity interiorpenalty method, Stokes, NavierStokes, Stabilization method1 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 Nonuniform Rational Bsplines (NURBS) are the industry standard. For analysissuitable 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 higherorder 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 infsup stable spaces for mixed formulations Boffi2013mixed (), a variety of compatible discretizations has been developed, most notably: TaylorHood 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 infsup 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 multiphysics problems with many different field variables, for which the derivation of infsup stable discrete spaces can be nontrivial. 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škaBrezzi 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: Galerkinleast squares and DouglasWang stabilization bazilevs2006IGAtheory () and variational multiscale stabilization (VMS) BazilevsCaloCottrell:2007 (). The structure of these approaches is that the stabilization is based on elementbyelement residuals. We note that recently a combination of VMS and compatible Bsplines is studied in vanOpstal2017 (). It is also noteworthy that for incompressible elasticity the use of infsup stable discretizations can be circumvented by using stream functions auricchio2007stream (), bar method elguedj2008 () and bar method antolin2017 ().
In this contribution we propose a novel skeletonbased stabilization technique for isogeometric analysis of viscous flow problems, like those described by the Stokes equations and incompressible NavierStokes equations with moderate Reynold’s numbers. The skeletonbased 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 highorder 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 Bsplines/NURBS with varying regularities, including the case of multipatch geometries.
The proposed stabilization technique only controls the th order derivative in the case of Bsplines/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 Bsplines/NURBS order – from which it is observed that the proposed technique optimally exploits the higherorder 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 oscillationfree 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 skeletonbased isogeometric analysis technique for the NavierStokes 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 skeletonstabilization 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 skeletonbased isogeometric analysis
To provide a setting for the skeletonbased stabilization proposed in Section 3 and to introduce the main notational conventions, we first present multipatch nonuniform rational Bspline (NURBS) spaces. We consider a domain (with or ) with Lipschitz boundary as exemplified in Figure 1. The domain is parameterized by a, possibly multipatch (), nonuniform rational Bspline (NURBS) such that
(1) 
where and are the patchwise geometric maps and parameter domains, respectively, with the parametric map defined as
(2) 
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 nondecreasing knot vectors, , with
(3) 
such that the number of basis functions per patch is , with the degree of the spline in the direction (). Note that for open Bsplines 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
(4) 
The superscript indicates the dependence of the partition on a mesh (resolution) parameter . We associate with the mesh the skeleton^{1}^{1}1This should not be confused with the topological skeleton concept in geometric modeling.:
(5) 
Note that since the skeletonbased stabilization technique considered in this work pertains to interelement continuity properties, the boundary faces are not incorporated in the skeleton. The skeleton (5) can be decomposed in the intrapatch skeleton, , and the interpatch skeleton, :
(6a)  
(6b) 
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 nonconforming mutipatch structure, multipatch 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 twodimensional 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):
(7) 
The regularity of the space across an intrapatch 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:
(8) 
For all functions the jumps of its th normal derivatives across an interface vanish in accordance with
(9) 
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 intrapatch 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 Skeletonstabilized Isogeometric Analysis for the NavierStokes equations
In this section we introduce the skeletonpenalty formulation for the NavierStokes equations in the context of Isogeometric Analysis. We commence with the formulation of the timedependent NavierStokes equations in Section 3.1. Next, we introduce the discrete skeletonpenalty formulation in Section 3.2.
3.1 The timedependent NavierStokes equations
We consider the unsteady incompressible NavierStokes 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 outwardpointing unit normal vector to is denoted by . For any time instant the NavierStokes equations for the velocity field and pressure field read:
(10) 
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):
(11) 
The trilinear, bilinear, and linear forms in this formulation are defined as
(12a)  
(12b)  
(12c)  
(12d) 
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:
(13) 
3.2 The Isogeometric SkeletonPenalty 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):
(14) 
The semidiscretization in space of the weak form (11) then reads:
(15) 
The pair of spaces in (14) does not satisfy the infsup condition, and hence the discretization in (15) is unstable. To stabilize the system, we propose to supplement the formulation with the skeletonpenalty term,
(16) 
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
(17) 
where and are two elements sharing the interface , and is the dimensional Hausdorff measure. The stabilized semidiscrete system – to which we refer as the isogeometric skeletonpenalty formulation for the NavierStokes equations – then reads:
(18) 
Remark 1.
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 skeletonpenalty term (16) in this case reads:
(19) 
Remark 3.
The formulation (18) based on the skeletonpenalty 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 Bsplines, , 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 highregularity 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 higherorder Bézier elements instead of higherorder 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:
(20) 
Similar to formulation (18), the isogeometric skeletonpenalty formulation for the Stokes equations reads:
(21) 
It is wellknown that problem (20) is the firstorder optimality condition for the saddle point of the Lagrangian functional (see e.g. Boffi2013mixed ())
(22) 
Analogously, the stabilized discrete system (21) is related to the optimization problem for the modified Lagrangian functional
(23) 
with
(24) 
The stabilized discrete system (21) follows directly from the firstorder 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 highorder derivatives of the pressure over the skeleton in a leastsquares sense.
4 The algebraic form of Skeletonstabilized Isogeometric Analysis
In this section we discuss various algorithmic aspects of the proposed skeletonbased stabilization framework. In Section 4.1 we briefly discuss the employed solution procedure for the unsteady NavierStokes 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 NavierStokes solution procedure
We employ a standard solution procedure for the unsteady NavierStokes equations. CrankNicolson 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 timedependence of the nonautonomous 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.
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
(27) 
and, accordingly, the approximate velocity field and pressure field can be written as
(28) 
where and are vectors of degrees of freedom. The corresponding algebraic form of (18) then reads
(29) 
with the matrix entries given by:
(30a)  
(30b)  
(30c)  
(30d)  
(30e)  
(30f) 
The algebraic form of the solution Algorithm 1 is presented in Algorithm 2.
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 highorder derivatives of the basis functions over can be evaluated. It should be noted that this skeleton structure is compatible with the recently proposed efficient rowbyrow assembly procedure for IGA calabro2017fastIGA ().
4.3 The k/complexity of the skeletonpenalty operator on sparsity pattern
The skeletonbased stabilization operator (16) affects the sparsity pattern of the discretized NavierStokes system due to the fact that the jump operators on the (higherorder) 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 Bspline 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 skeletonpenalty matrix associated with the operator .
The bandwidth^{2}^{2}2The bandwidth is defined as the smallest nonnegative integer such that if . of the skeletonpenalty 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 elementbased stabilization techniques. By construction, Bspline bases ameliorate this issue in the sense that even at full continuity the bandwidth of the skeletonpenalty 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 Skeletonstabilized 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 twodimensional Stokes problem – i.e., problem (11) without timedependent and convective terms – in the unit square domain . The body force is taken in accordance with the manufactured solution BuffaStokes2011 ():
(31a)  
(31b) 
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 wellposedness. We use a Lagrange multiplier approach to enforce this condition.
In Figure 3 we study the asymptotic convergence behavior of the proposed method for Bsplines 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 oscillationfree 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 infsup 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.
In Figure 4 we study the sensitivity of the computed results with respect to the SkeletonPenalty stabilization parameter . The convergence behavior of the solution using continuous quadratic Bsplines 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 SkeletonPenalty 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 illconditioning 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 .
The performance of the proposed Skeletonstabilized IsoGeometric Analysis framework is further studied based on the generalized Stokes equations with homogeneous Dirichlet boundary conditions:
(32) 
This system – for which the body force is selected in accordance with the manufactured solution (31) – is characterized by the Damköhler number
(33) 
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 Bsplines 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 nonreactive 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.



To understand the effect of reduced regularity – which is of particular importance in the case of multipatch models – we first study the Bspline discretization of the Stokes problem on the unit square with varying intrapatch 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 interiorpenalty method burman2006edge ().



5.2 Steady Stokes flow in a quarter annulus ring
To demonstrate the performance of the proposed SkeletonPenalty 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 ()