A hierarchical approach to the a posteriori error estimation of isogeometric Kirchhoff plates and Kirchhoff-Love shells

TL;DR

This paper introduces a hierarchical a posteriori error estimator for isogeometric Kirchhoff plates and shells, enabling adaptive refinement with optimal convergence rates and accurate error approximation in complex structural simulations.

Contribution

It presents a novel hierarchical spline-based error estimator for fourth-order PDEs, specifically tailored for Kirchhoff plates and shells, improving adaptive simulation accuracy.

Findings

Achieves optimal convergence rates in energy norm

Provides accurate true error approximation

Effective on problems with smooth and singular solutions

Abstract

This work focuses on the development of a posteriori error estimates for fourth-order, elliptic, partial differential equations. In particular, we propose a novel algorithm to steer an adaptive simulation in the context of Kirchhoff plates and Kirchhoff-Love shells by exploiting the local refinement capabilities of hierarchical B-splines. The method is based on the solution of an auxiliary residual-like variational problem, formulated by means of a space of localized spline functions. This space is characterized by $C^1$ continuous B-splines with compact support on each active element of the hierarchical mesh. We demonstrate the applicability of the proposed estimator to Kirchhoff plates and Kirchhoff-Love shells by studying several benchmark problems which exhibit both smooth and singular solutions. In all cases, we obtain optimal asymptotic rates of convergence for the error measured in the energy norm and an excellent approximation of the true error.

Equations

$${\mathcal{B}}_{\mathbf{i},\mathbf{p}}=\prod _{r=1}^{d_r}{b}_{i_r,p_r}^r$$

$$\mathbf{S}=\sum _{\mathbf{i}}{\mathcal{B}}_{\mathbf{i},\mathbf{p}}{\mathbf{P}}_{\mathbf{i}},$$

$${\mathcal{Q}}^l=\left\{{\mathcal{Q}}_i^l|{\mathcal{Q}}_i^l=({\hat{\xi }}_i^l,{\hat{\xi }}_{i+1}^l),i=0,\dots ,{\hat{n}}^l-1\right\}$$

$${\mathcal{B}}_{-}^l$$

$${\mathcal{B}}_{+}^l$$

$${\mathcal{B}}_{=}^l$$

$$\mathcal{H}\mathcal{B}$$

$${\mathcal{H}\mathcal{B}}^l$$

$$D{\Delta }^{2}u$$

$$u$$

$$-\nabla u\cdot \mathbf{d}$$

$$\nu D\Delta u+(1-\nu )D\mathbf{d}\cdot (\nabla \nabla u)\mathbf{d}$$

$$D(\nabla (\Delta u)+(1-\nu )\Psi (u))\cdot \mathbf{d}$$

$$D=\frac{E{t}^3}{12(1-{\nu }^2)},$$

$$a(u,v)=F(v)\forall v\in V,$$

$$a(u,v)$$

$$F(v)$$

$$V_h=\text{span }\{b\circ {\mathbf{S}}^{-1}|b\in \mathcal{H}\mathcal{B}\}.$$

$$a(u_h,v_h)=F(v_h)\forall v_h\in V_h,$$

$$u_h=\sum _{b\in \mathcal{H}\mathcal{B}}b{c}_b,$$

$$a(\mathbf{u},\mathbf{v})=F(\mathbf{v})\forall \mathbf{v}\in V,$$

$$a(\mathbf{u},\mathbf{v})$$

$$F(\mathbf{v})$$

$$V_h$$

$$a\left({\mathbf{u}}_h,{\mathbf{v}}_h\right)=F\left({\mathbf{v}}_h\right)\forall {\mathbf{v}}_h\in V_h.$$

$${\tilde{V}}_h^p=V_h^p\oplus W_h^{p+1},$$

$$W_h^{p+1}=\bigcup _{l=0}^N{W}_{h,l}^{p+1},$$

$$W_{h,l}^{p+1}=\text{span }\{b\circ {\mathbf{S}}^{-1}|b\in {\mathcal{B}}_{a,p+1}^l\cap C^1(\overline{{\hat{\Omega }}^l})\text{ and }\exists ϵ\subseteq {\hat{\Omega }}^l:\text{ supp}(b)=\overline{ϵ}\}.$$

$$a(e_h,b_h)=F(b_h)-a(u_h,b_h)\forall b_h\in W_h^{p+1},$$

$${\eta }_{ϵ}=C_a{\Vert e_h\Vert }_{E(ϵ)}\forall ϵ\in Q,$$

$${\eta }_{ϵ}>\gamma {\tilde{\eta }}_{\epsilon },\text{where}{\tilde{\eta }}_{\epsilon }=\underset{\epsilon \in Q}{\max }{\eta }_{\epsilon },$$

$$u$$

$$\nu D\Delta u+(1-\nu )D\mathbf{d}\cdot (\nabla \nabla u)\mathbf{d}$$

$$g=64{\pi }^4\sin (2\pi x)\sin (2\pi y).$$

$${\Vert u_{ex}-u_h\Vert }_{H^s(\Omega )}\le C{h}^{p+1-s}{\Vert u_{ex}\Vert }_{H^{p+1}(\Omega )}p+1>s,$$

Full text

A hierarchical approach to the a posteriori error estimation of isogeometric Kirchhoff plates and Kirchhoff-Love shells

Pablo Antolin

Institute of Mathematics, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Annalisa Buffa

Institute of Mathematics, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Istituto di Matematica Applicata e Tecnologie Informatiche ‘E. Magenes’ (CNR), Pavia, Italy.

Luca Coradello luca.coradello@epfl.ch, Corresponding Author Institute of Mathematics, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

Abstract

This work focuses on the development of a posteriori error estimates for fourth-order, elliptic, partial differential equations. In particular, we propose a novel algorithm to steer an adaptive simulation in the context of Kirchhoff plates and Kirchhoff-Love shells by exploiting the local refinement capabilities of hierarchical B-splines. The method is based on the solution of an auxiliary residual-like variational problem, formulated by means of a space of localized spline functions. This space is characterized by $C^{1}$ continuous B-splines with compact support on each active element of the hierarchical mesh. We demonstrate the applicability of the proposed estimator to Kirchhoff plates and Kirchhoff-Love shells by studying several benchmark problems which exhibit both smooth and singular solutions. In all cases, we obtain optimal asymptotic rates of convergence for the error measured in the energy norm and an excellent approximation of the true error.

Keywords: isogeometric analysis, a posteriori error estimator, adaptivity, hierarchical B-splines, shells, plates.

1 Introduction

IsoGeometric Analysis (IGA) has been a thriving area of research since the first pioneering work [32] was published in 2005. Employing smooth B-splines, NURBS or variances thereof as basis functions for the solution field has shown excellent properties in many mathematical and engineering applications. For instance, in many cases higher global continuity yields a better accuracy per degree-of-freedom (dof) and it allows to discretize higher-order variational problems in their primal form, e.g. Kirchhoff plates and Kirchhoff-Love shells [33] or Cahn-Hilliard problems [28]. For a detailed review of the method and its applications, the reader is referred to [32, 19, 1].

However, in order to make IGA competitive for industrial applications, a crucial feature is the capability of building locally refined bases suitable for analysis in an efficient and robust manner. Indeed, we remark that the tensor-product nature of B-splines is broken by local refinement. To tackle these issues, several alternatives have been proposed such as hierarchical B-splines (HB) [23, 29, 35] and their recent variant denoted by truncated hierarchical B-splines (THB) [26, 27], T-splines [8, 44, 9] and locally-refinable splines (LR-splines) [21, 11]. Also the inverse operation of refinement (known as coarsening), which plays an important role for the computational efficiency of dynamic problems, proves to be a non-trivial task and another booming area of research, see for instance [38, 25, 31, 16]. In the framework of this work, we employ hierarchical splines as they provide the most straightforward construction and an easy implementation in the context of isogeometric analysis.

In order to automate the process of local refinement, a posteriori error estimators have been successfully used to steer adaptive simulations for standard finite elements, where the literature is rich and several families of estimators have been thoroughly studied (see for instance [5, 2, 48] and references therein). We remark that a good estimator should properly capture and resolve local features of the solution that cause the error to be much larger than in the rest of the computational domain. In the realm of isogeometric analysis, the field of error estimation is still in its early stages of development. Adaptivity using T-splines has been presented in the pioneering work [22]. More results on recovery-based error estimators have been published in [36, 37] where refinement is achieved by using LR-splines. Additionally, functional-type error estimates have been studied in [34] but not within a real adaptive loop. More recently, also recovery-based methods for polynomial splines over T-meshes have been proposed [3]. In the context of (truncated) hierarchical B-splines, we mention [13, 24] where residual-based element-wise estimators have been analyzed for linear second-order elliptic problems. These concepts have been extended in [12] for function-based estimators.

The main goal of this paper is to introduce, building upon the work of [6, 49] on multi-level estimators, a simple a posteriori error estimator for linear fourth-order elliptic partial differential equations (PDEs) of engineering relevance, such as Kirchhoff plates and Kirchhoff-Love shells, by making use of hierarchical B-splines. The proposed method is based on the solution of an additional residual-like problem which is formulated on a so-called bubble space. Given the locality of the functions in the aforementioned space, the resulting linear system is in general small and easy to solve. We remark that by exploiting this technique we bypass completely the computational burden stemming from the evaluation of higher order derivatives and their jumps across edges, which would be required in classical residual-based error indicators. This is particular advantageous for Kirchhoff-Love shells, where the terms associated to higher derivatives involve the cumbersome and error-prone task of taking (nested) covariant derivatives of quantities of interest, see for instance [39]. Indeed, the use of residual-based error estimators for Kirchhoff-Love shells is computationally expensive and an efficient implementation is far from being trivial. Moreover, contrary to residual-based techniques, we do not require the computation of additional integrals which take into account the jumps of derivatives across element interfaces. These features play again in favor of the computational efficiency and simplicity of the proposed methodology.

A known drawback of this family of estimators is that their reliability and efficiency is subjected to the saturation assumption onto the underlying augmented space. As noted in [6], this assumption is, in general, problem dependent and it can, potentially, become crucial. However, we highlight that our focus is to find a good indicator for steering adaptive simulations in the scope of structural mechanics and all our numerical experiments confirm that the proposed estimator performs well in this task. Furthermore, it is worth noting that in all our tests the estimator provides also an excellent approximation of the true error.

The structure of the paper is as follows. Section 2 introduces the fundamental concepts needed for hierarchical B-splines. Section 3 briefly recalls the Kirchhoff plate and Kirchhoff-Love formulations. Then, Section 4 discusses the proposed bubble error estimator and a possible algorithm. Section 5 presents several numerical examples. Finally, some conclusions are drawn in Section 6.

2 IsoGeometric Analysis

In this Section, the notation and basic concepts and definitions related to B-splines and NURBS are reviewed, following closely the derivation in [20]. For further details, the reader is referred to [40, 19], and references therein.

A B-spline basis function of degree $p$ is generated starting from a non-decreasing sequence of real values that corresponds to a set of coordinates in the parameter space referred to as knot vector, denoted in the following as $\Xi$. Given a knot vector, a univariate B-spline basis function $b_{i,p}$ can be defined recursively using the Cox-de Boor algorithm [40]. It is worth mentioning that the smoothness of the obtained B-spline basis is $C^{p-k}$ at every knot, where $k$ denotes the multiplicity of the considered knot, while it is $C^{\infty}$ everywhere else.

The definition of multivariate B-splines is achieved in a straight-forward manner using the tensor product of univariate B-splines. Indeed, the multivariate B-splines $\mathcal{B}_{\mathbf{i},\mathbf{p}}$ are obtained as:

[TABLE]

with $r=1,...,d_{r}$, where $d_{r}$ is the dimension of the parameter space and $p_{r}$ denotes the polynomial degree in the parametric direction $r$, respectively. Additionally, the multi-index $\mathbf{i}=\left\{i_{1},...,i_{d_{r}}\right\}$ denotes the position in the tensor product structure and $\mathbf{p}=\left\{p_{1},...,p_{r}\right\}$ indicates the vector of polynomial degrees.

Lastly, let us define a B-spline parametrization $\mathbf{S}$ as a linear combination of multivariate B-spline basis functions and control points as follows:

[TABLE]

where the coefficients $\mathbf{P}_{\mathbf{i}}\in\mathbb{R}^{d}$ of the linear combination are the control points and $d$ denotes the dimension of the physical space. We remark that all the concepts summarized here can be readily transferred to NURBS, for further details we refer to [19]. In the rest of the paper, without loss of generality, the degree $p$ will be considered equal in each parametric direction and will be omitted from the notation.

2.1 Hierarchical B-splines

In this Subsection, the concept of hierarchical B-spline basis, denoted by $\mathcal{HB}$, is introduced. This allows us to build a basis that is locally refinable and therefore to overcome the limitations intrinsic to the tensor-product nature of B-splines and NURBS.

Let $V^{0}\subset V^{1}\subset\cdots\subset V^{N}$ be a sequence of nested spaces of splines defined on a parametric domain $\hat{\Omega}$, where, to keep the notation simple, we consider the one-dimensional case. Each space $V^{l}$, $l=0,\ldots,N$ is spanned by the B-spline basis $\mathcal{B}^{l}$ of degree $p$, associated to level $l$ and the corresponding knot vector $\Xi^{l}$.

Our task now is to identify a set of analysis-suitable functions $\mathcal{N}\subset\bigcup_{l}\mathcal{B}^{l}$. To this end, let us define as isogeometric elements $Q$ a partition of $\hat{\Omega}$. In particular, given $\hat{n}^{l}\in\mathbb{N}$, let $\hat{\Xi}^{l}=(\hat{\xi}^{l}_{0},\ldots,\hat{\xi}^{l}_{\hat{n}^{l}})$ be the knot vector composed of non-decreasing knots of $\Xi^{l}$ without repetition, and let:

[TABLE]

be the set of open intervals constituting the non-empty knot spans of $\Xi^{l}$. The elements of the multi-level mesh can be any non-overlapping partition $Q\subset\bigcup_{l}\mathcal{Q}^{l}$ of $\hat{\Omega}$, such that $\bigcup_{\epsilon\in Q}\bar{\epsilon}=\hat{\Omega}$. Moreover, let $Q^{l}=Q\cap\mathcal{Q}^{l}$ be the elements of level $l$, and let $\hat{\Omega}^{l}=\bigcup_{\epsilon\in Q^{l}}\bar{\epsilon}$ denote their domain. Furthermore, let $\hat{\Omega}^{l}_{+}=\bigcup_{l^{*}=l+1}^{N}\hat{\Omega}^{l^{*}}$ be the refined domain with respect to level $l$ and, analogously, let $\hat{\Omega}^{l}_{-}=\bigcup_{l^{*}=0}^{l-1}\hat{\Omega}^{l^{*}}$ be the coarser domain. Given a set of elements $Q$, we still need to define a set of basis functions $\mathcal{N}\subset\bigcup_{l}\mathcal{B}^{l}$ suitable for analysis. To this end, we consider the set of functions with support on the elements of level $l$ as the set of active functions $\mathcal{B}^{l}_{a}=\left\{b\;|b\in\mathcal{B}^{l},\,\textnormal{supp}(b)\cap\hat{\Omega}^{l}\neq\emptyset\right\}\subset\mathcal{B}^{l}$. Among these, a subset of linearly independent functions has to be chosen. Following [30, 20], we partition $\mathcal{B}^{l}_{a}$ into:

[TABLE]

Then, the hierarchical B-spline basis $\mathcal{HB}$ [49, 35] is defined as follows:

[TABLE]

Namely, $\mathcal{HB}$ is the set of B-splines of each level $l$ whose support intersects only elements of level $\tilde{l}\geq l$ and at least one element of level $l$. It was proven in [49] that this set is composed of linearly independent functions and therefore suitable for the analysis.

3 Formulations

In this section, the Kirchhoff plate and the Kirchhoff-Love shell formulations are stated in the context of IGA. The reader is referred to [42, 33, 18, 39] for a more comprehensive review of the formulations.

3.1 The Kirchhoff plate formulation

Here, following the derivation in [42], we introduce the Kirchhoff plate problem, governed by the bi-Laplace differential operator. Let us define an open set $\Omega\subset\mathbb{R}^{2}$ with a sufficiently smooth boundary $\partial\Omega$, such that the normalized normal vector $\boldsymbol{d}$ is well-defined (almost) everywhere. Additionally, we assume that the boundary $\Gamma=\partial\Omega$ can be decomposed into $\Gamma=\mkern 0.1mu\overline{\mkern-0.1mu\Gamma_{w}\cup\Gamma_{Q}\mkern-0.1mu}\mkern 0.1mu$ and $\Gamma=\mkern 0.1mu\overline{\mkern-0.1mu\Gamma_{\phi}\cup\Gamma_{M}\mkern-0.1mu}\mkern 0.1mu$, such that $\Gamma_{w}\cap\Gamma_{Q}=\varnothing$ and $\Gamma_{\phi}\cap\Gamma_{M}=\varnothing$, respectively. We formulate the strong form of the problem as follows:

[TABLE]

where $u$ represents the deflection of the plate, $D$ its bending stiffness, $\nu$ is the Poisson ratio, $g$ is the load per unit area in the thickness direction, $u_{\Gamma}$, $\phi_{\Gamma}$, $M_{\Gamma}$ and $Q_{\Gamma}$ are the prescribed deflection, rotation, bending moments and effective shear, respectively. The bending stiffness $D$ is defined as:

[TABLE]

where $E$ is the Young modulus and $t$ denotes the thickness of the plate, which without loss of generality we suppose to be a constant in $\Omega$. Using variational calculus, the weak form of problem (2) reads: find $u\in\>V$ such that

[TABLE]

where $V\subset H^{2}(\Omega)$ denotes an infinite-dimensional space that depends in general on the boundary conditions of the problem at hand, for further details we refer to [17]. For simplicity and without loss of generality, in the following we assume the Poisson ratio $\nu$ to be zero and we choose $E$ and $t$ such that $D=1$. Consequently, the bilinear form $a:V\times V\rightarrow\mathbb{R}$ can be expanded as follows:

[TABLE]

and similarly the linear functional $F:V\rightarrow\mathbb{R}$ reads:

[TABLE]

It can be shown that (3) satisfies the conditions of the Lax-Milgram theorem, therefore $u\in V$ is the unique solution of the variational problem (3).

Finally, we discretize the weak form (3) employing the Bubnov-Galerkin method and using as basis the hierarchical B-splines $\mathcal{HB}$. Let $V_{h}\subset V$ be a finite dimensional subspace defined as:

[TABLE]

where $\mathbf{S}$ is the spline parametrization defined in Section 2. We remark that due to the requirements on the discrete admissible space for the deflection, the basis functions must be at least $C^{1}$ continuous globally such that the bilinear form is well-defined. This requirement is easily fulfilled within one patch by the use of (hierarchical) B-splines of maximum continuity (of degree $p\geq 2$). Then, the discrete weak form can be written as: find $u_{h}\in\>V_{h}$ such that

[TABLE]

where $u_{h}$ represents the discrete solution defined as:

[TABLE]

with $c_{b}\in\mathbb{R}$ denoting the unknown control variables and $b\in\mathcal{HB}$ the corresponding hierarchical basis introduced in Section 2.

3.2 The Kirchhoff-Love shell formulation

In the following, the Kirchhoff-Love formulation in its weak form is briefly introduced, following closely the derivation and the notation used in [33, 18].

3.2.1 The weak formulation

Similarly to what has been described in Section 3.1 for the Kirchhoff plate, starting from the strong formulation of the Kirchhoff-Love problem (for further details see [39]) and applying variational calculus, the weak form reads: find $\boldsymbol{u}\in\>V$ such that

[TABLE]

where the choice of space $V\subset H^{2}(\Omega)$ depends on the boundary conditions of the problem at hand. Moreover, $a:V\times V\rightarrow\mathbb{R}$ is a continuous and strongly coercive bilinear form and $F:V\rightarrow\mathbb{R}$ is a continuous linear functional. They can be expanded, respectively, as follows:

[TABLE]

where $\boldsymbol{\varepsilon}$, $\boldsymbol{\kappa}$ denote the membrane and bending strain tensors, respectively, and $\boldsymbol{n}$, $\boldsymbol{m}$ are their energetically conjugate stress resultants. Additionally, $\boldsymbol{b}$ is the applied body load, $\boldsymbol{p}$ and $\boldsymbol{r}$ represent the applied traction force and traction moment, respectively. For a detailed derivation, we refer to [33, 4]. Note that the weak form (5) satisfies the conditions of the Lax-Milgram theorem and therefore the solution $\boldsymbol{u}\in V$ exists and it is unique.

Employing the IGA paradigm, we utilize the same (hierarchical) B-splines and NURBS basis functions for the geometry description and for the field approximation. Let us recall the definition of the discrete space:

[TABLE]

by employing the Bubnov-Galerkin discretization of (5), we get the following discrete weak formulation of the Kirchhoff-Love shell problem: find $\boldsymbol{u}_{h}\in V_{h}$ such that

[TABLE]

Again, it is worth noting that due to the requirements on the discrete admissible space for the displacement field, the basis functions must be at least $C^{1}$ continuous globally such that the bending operator is well-defined.

Finally, we remark that in most cases, due to efficiency, the parametrization $\mathbf{S}$ can be expressed in terms of tensor product B-splines of the coarsest level $\mathcal{B}^{0}$ instead of the hierarchical basis $\mathcal{HB}$ [25].

4 A posteriori error estimator for Kirchhoff plates and Kirchhoff-Love shells

In the following Section we introduce a variant of the error estimator analyzed in [6] and we extend its isogeometric version, proposed in [22], to fourth-order PDEs of mechanical relevance, such as Kirchhoff plates and Kirchhoff-Love shells. We choose these models since, due to the Kirchhoff hypotheses, we avoid a priori the negative effects related to shear locking. Then, later in the Section, we show a possible implementation of the proposed indicator in the context of hierarchical B-splines.

4.1 The bubble error estimator

This family of a posteriori estimators was introduced in [6] in the context of the p-version of the Finite Element Method and has been successfully applied to T-splines in [22] for second-order boundary value problems. To introduce the concept in the context of hierarchical IGA, let us define the finite dimensional solution space $V^{p}_{h}$ as the span of hierarchical B-splines basis functions of order $p$, where we restrict our analysis to the case $p>2$ (see Remark 2 for the case $p=2$). Then, let us recall the discrete solution $u_{h}\in V^{p}_{h}$, where without loss of generality we make no distinction between the plate or the shell problem. We now assume there exists a larger space $V^{p}_{h}\subset\widetilde{V}^{p}_{h}\subset V$ such that the following decomposition holds:

[TABLE]

where $W^{p+1}_{h}$ is the space in which we seek a good estimation of the error $e_{h}\approx e=\left\lVert u-u_{h}\right\rVert$, in a suitable norm $\left\lVert\cdot\right\rVert$. In particular, recalling some definitions provided in Section 2, let us characterize $W^{p+1}_{h}$ as follows:

[TABLE]

where

[TABLE]

In words, $W^{p+1}_{h,l}$ is the space spanned by active B-splines (or better in this case, Bernstein polynomials) of level $l$ and degree $p+1$ such that their support is compact and overlaps exactly with one active element $\epsilon$ of level $l$, e.g. see Figure 1 where one level is depicted. We are now ready to define the a posteriori error estimate as: find $e_{h}\in\>W^{p+1}_{h}$ such that

[TABLE]

where $b_{h}$ is referred to in the following as bubble functions. We remark that we require these functions to be at least $C^{1}$ on their support such that the corresponding bilinear form $a$ is well-defined. Notice that due to the compact support of $b_{h}$, the linear system corresponding to the discrete error weak form (8) is block diagonal and in general easy to solve. Moreover, we highlight that the evaluation of the bilinear form $a$ appearing in (8) follows the same steps as the original plate/shell problem, where different test functions $b_{h}$ are used. Hence, setting up the bubble problem requires only minor changes into an existing isogeometric code.

Once we have found the unknown coefficients $e_{h}$, we define the element-wise error estimator $\eta_{\epsilon}$ as follows:

[TABLE]

where we recall that $\epsilon$ is an isogeometric element of the hierarchical mesh $Q$ and $\left\lVert\cdot\right\rVert_{E(\epsilon)}$ denotes the energy norm restricted to element $\epsilon$. Additionally, we introduce $C_{a}$ as an empirical constant which, considering all our numerical experiments, we claim is independent from the chosen degree and from the problem at hand. In the following, we set $C_{a}=3$ once and for all in order to shift the estimator above the true error.

At this point, it is worth highlighting that the proposed error estimator is computationally inexpensive, easy to implement and is embarrassingly parallelizable due to the choice of the disjoint bubble space $W^{p+1}_{h}$. A possible algorithm to estimate the error is summarized in Algorithm 1, where we design the code to solve the additional problem (8) level-wise over the hierarchical mesh $Q$.

Remark 1** (Inhomogeneous boundary conditions).**

On one hand, we highlight that whenever inhomogeneous boundary conditions of Neumann-type are applied to the problem at hand, problem (8) must also contain the corresponding additional terms in the right-hand-side. In these cases, we augment the bubble space $W^{p+1}_{h}$ with suitable boundary bubbles to properly capture the error on the imposition of natural boundary conditions. An example is depicted in Figure 1 for bubbles of degree $p=4,5$. On the other hand, we make the assumption that the error on the imposition of inhomogeneous Dirichlet-type boundary conditions is negligible and therefore no additional shape functions are needed for these cases.

Remark 2** (The case $p=2$).**

We remark that isogeometric Kirchhoff-Love elements notoriously suffer from membrane locking phenomena, specially for low degree discretization [33], and it can be considered common practice to use the coarsest (and lowest degree) description of the geometry that is still exact while employing a k-refined space for the solution field to mitigate these effects. Nevertheless, for the sake of completeness, we should highlight that the construction of the bubble space $W^{p+1}_{h}$ is not straightforward in the case where the underlying B-splines are quadratic. We propose as possible remedies:

Construction of an additional sub-grid obtained as dyadic refinement of the original hierarchical mesh and definition of the bubble functions on the corresponding space defined by the new grid. 2. 2.

Introduction of suitable edge functions in the bubble space.

Notice that, although both approaches could be feasible, they add a non-negligible complexity to the method. In the first case, a new, finer hierarchical grid has to be created at every iteration of the adaptive loop and numerical integration of (8) has to be performed over it. In the second case, the definition of which edge bubbles have to be considered in the hierarchical case is not unique and considerably worsen the simplicity of the method. Indeed, the corresponding linear system would become more difficult to solve since the element-wise locality of the bubble functions would be lost. Therefore, in light of these considerations, we will not address further the case $p=2$ in the scope of this work.

Remark 3** (Jumps in the residual).**

In the last remark, we notice that, although the residual exhibits jumps across element edges for the discretization of degree $p=3$, in all our computations the estimator behaves still optimally even if the corresponding jump terms are not accounted for in (8).

4.2 Marking and refine strategies

Here, we briefly recall that once an estimate of the error element-wise is available, there exist several strategies proposed in the literature to mark elements for refinement. In all our numerical examples we use the so-called maximum strategy, which can be summarized as follows. Let $\gamma\in(0,1)$ be a user-defined threshold, all elements such that

[TABLE]

are marked for refinement. In the following, if not stated otherwise, we will always use $\gamma=0.5$. Additionally, since for basis functions of order $p$ their support usually spans up to $p+1$ elements, we try to reduce the number of isolated marked elements. Indeed, in general, isolated refined elements do not improve the solution space $V^{p}_{h}$ (they do not add any additional degree of freedom) but only the associated numerical integration is refined. To avoid that, we slightly modify the maximum strategy marking algorithm such that once an element $\epsilon$ of level $l$ is marked for refinement, the algorithm also marks all the neighbors of $\epsilon\in Q^{l}$, where $Q^{l}$ is the set of active elements of level $l$. Finally, we remark that our marking algorithm is also designed to preserve the admissibility (as defined in [15]) of the hierarchical mesh between consecutive iterations of the adaptive loop. We recall that admissible meshes are guaranteed to have a bounded number of basis functions acting on any element of the mesh, such that interaction between very fine and very coarse functions is avoided. If not stated otherwise, we will set the class of admissibility $m$ to be $m=p-1$ in all our numerical examples. For the sake of brevity we do not discuss the details here, the reader is referred to [15, 14, 10] for a comprehensive review of the concept of admissibility.

5 Numerical Examples

All the numerical experiments presented in the following section have been implemented in the open-source and free Octave/Matlab package GeoPDEs [47], which is a software suite for the solution of partial differential equations specifically designed for isogeometric analysis.

5.1 Kirchhoff plate

In this Section we present several adaptive computations in the context of Kirchhoff plates, thought to give a first assessment of the performance of the proposed error estimator. We demonstrate the applicability of our method both for problems which exhibit smooth and singular solutions. In all cases presented here, the proposed error estimator shows excellent performance in steering the adaptive simulation, yielding the expected optimal rates of convergence in the asymptotic regime.

5.1.1 Smooth solution on a square plate

In the first example, we analyzed the behavior on the bubble error estimator compared to a classical residual-based type error estimator. We define the computational domain to be the unit square $\Omega=[0,1]^{2}$ and we impose the following homogeneous boundary conditions:

[TABLE]

Additionally, the applied load $g$ is constructed such that it fulfills the following manufactured solution $u_{ex}=\sin(2\pi x)\sin(2\pi y)$. Namely, $g$ reads as follows:

[TABLE]

The structure, physical parameters and the exact solution are depicted in Figure 2.

Notice that, as predicted by classical a priori error estimates (see for instance [17, 41] or in the context of IGA [7]) and since $u_{ex}$ is sufficiently regular, we expect the following convergence rate:

[TABLE]

where $p$ denotes the polynomial degree of the hierarchical basis and $h$ represents the maximum mesh size, namely:

[TABLE]

Additionally, in all our numerical examples, we choose $s=2$. Moreover, for simplicity and without loss of generality, we measure the error in the $H^{2}$ semi-norm, which corresponds to the energy norm of the problem at hand. It is worth noting that in most cases we will use the number of degrees-of-freedom (dofs) instead of the mesh size $h$, since we can relate these two quantities via $h\approx\text{dofs}^{-1/d}$, where $d$ denotes the dimensionality of the problem.

For the sake of comparison, in the case of Kirchhoff plates only, we drive the adaptive method by means of both the residual-based and the bubble-based estimators. To test our implementation, we perform at first uniform refinement for different degrees $p=3,4,5$ and check the convergence rate of the error in the energy norm and both estimators against the element size $h$, as depicted in Figure 3 for the bubble and residual-based, respectively. In all the presented cases an optimal asymptotic rate of convergence is observed, both for the error and the estimator. Furthermore, we remark that the bubble estimator is much closer to the real error compared to the residual-based.

Finally, we run the same example letting the bubble and residual-based estimators drive automatically the adaptive simulation. The results are reported in Figure 4, where the error in energy norm and the estimator are plotted against the square root of the number of dofs. We observe analogous results to the uniform refinement case, where for $p=3,4,5$ the expected asymptotic rates of convergence are achieved. Then, let us define the effectivity index $\theta$ as:

[TABLE]

where the optimal value of $\theta$ is $1$. Also in the adaptive example, the bubble estimator provides a better estimate of the real error compared to a classical residual-based estimator, yielding an effectivity index much closer to the optimal value, as shown in Figure 5. This statement holds in an analogous way for all the presented examples and shows numerically the efficiency of the proposed method.

5.1.2 Singular solution on a square plate

In the next example, we consider again the computational domain to be a unit square $\Omega=[0,1]^{2}$, see Figure 6.

However, this time the manufactured solution $u_{ex}=x^{\alpha}y^{\beta}$ with $\alpha=\beta=2.8$ is constructed such that a singularity is present along the bottom and the left edges of the plate. In particular, it holds $u_{ex}\in H^{3}(\Omega)\setminus H^{4}(\Omega)$. The applied load $g$ is again constructed such that it fulfills the manufactured solution $u_{ex}$ and it is given as follows:

[TABLE]

The boundary conditions are also constructed from the exact solution as:

[TABLE]

where $M_{ex}$ denotes the exact bending moment. The reduction in regularity of the solution limits the rate of convergence in case of uniform refinement, even for increasing $p$. From [17, 41], the result on a priori convergence for uniform refinement is in this case a more general version of (10) and reads:

[TABLE]

where $r$ represents the regularity of the exact solution (for this example we have $r\in\left]3,4\right[$) and $l$ is defined as $l=\min{(r,p+1)}.$ This effect can be clearly seen in Figure 7 for $p=3,4$. Figure 7 also shows the results for an adaptive simulation driven by the bubble estimator, where optimal rates of convergence are recovered by our method and a significant increase in accuracy per degree-of-freedom is achieved. Additionally, the obtained mesh at different iterations of the adaptive algorithm are depicted in Figure 8, where it can be seen that the singularities are accurately detected and resolved by the estimator.

5.1.3 Point load on a rectangular plate

In this example, we test the performance of the proposed error estimator when a point load is applied as an external force. We perform again our computation on a unit square $\Omega=[0,1]^{2}$ and we suppose that the plate is simply supported on the entire boundary $\partial\Omega$. Additionally, we apply the external point load at the center of the structure. There exists an analytical solution for the deflection under the load, given as [43]:

[TABLE]

where $g$ represents the applied external force, $L$ is the length of the plate and $D$ denotes its flexural stiffness. We set $g=-1$, $L=1$ and the physical parameters such that $D=1$. Computing (11) with an adequate number of terms for the double Fourier series yields a reference value of $u_{ex}=-0.011600839735872\ldots$ for the deflection. In Figure 9(a) the convergence of the normalized deflection is plotted against the number of degrees-of-freedom, where the former is defined as $|1-u_{h}/u_{ex}|$. Moreover, the obtained solution and hierarchical mesh obtained after $k=9$ steps of the adaptive algorithm are depicted in Figure 9(b). Here, the advantages of using an adaptive scheme are clearly highlighted, in terms of efficiency and accuracy per degree-of-freedom. For instance, at around 2500 dofs the adaptive strategy is already two orders of magnitude more accurate compare to the results obtained by uniform refinement.

Remark 4** (Regularization of the Dirac delta).**

It is worth noting that a point load is modeled as a Dirac delta, which rigorously speaking is a distribution. Therefore, for the evaluation of the residual in a strong sense, needed in the residual-based estimator, we cannot directly use it inside our computations but instead we must regularize it. In our example, this is achieved with a steep Gaussian function. A considerable amount of literature has been written on how to regularize and correctly integrate the Dirac delta according to the corresponding application (e.g. we refer to [46] and references therein) and although our approach might not be optimal we feel it suits the purpose of our test case while at the same time it proves to be easily implementable.

5.2 Kirchhoff-Love shell

In this Section we demonstrate the applicability of the bubble error estimator to the Kirchhoff-Love formulation. Again, in all our numerical experiments the estimator performs well in driving adaptive computations. Once more, it is worth highlighting that with the proposed method we avoid the evaluation of higher order derivatives (this time defined on a manifold) which turns out to be particularly challenging and computationally involved in the framework of Kirchhoff-Love shells, see for instance [39].

5.2.1 Gravity load on a pinned roof

In this example, we study the behavior of a structure subjected to its self-weight where pinned supports are applied to the entire boundary. The geometry, boundary conditions and physical parameters are given in Figure 10, where the problem setup (except for the boundary conditions) is taken from the classical Scordelis-Lo benchmark [45]. The boundary conditions are modified in order to avoid numerical issues during the computation of the reference solution and corresponding error, since the original problem definition of the Scordelis-Lo roof is not well-posed, see for instance [33].

To the best of the authors’ knowledge, there is no closed-form global solution available in the literature for the problem at hand. Therefore, all the convergence studies presented in the following assess the behavior against a reference solution $\boldsymbol{u}_{h}^{\text{ref}}$, which was computed using B-splines of degree $p=8$ on a fine uniform mesh with $\simeq 200\,000$ degrees-of-freedom. In particular, we define an approximation of the error in the energy norm as:

[TABLE]

where the numerical integration of this quantity is performed on the fine mesh.

In Figure 11 the convergence behavior of $||\tilde{e}||$ is depicted for the case of uniform refinement and adaptive refinement driven by the bubble estimator, against the number of dofs. We notice that since the true solution of the problem is regular enough, we obtained the optimal asymptotic rate of convergence, both for the error in energy norm and the estimator.

This example is thought as a first assessment of the performance of the proposed estimator on shell geometries. We remark that in standard shell benchmarks adopted in the literature, convergence is usually tested against a reference value for the displacement in some points of interest of the structure. While this information is relevant in many engineering applications, we feel that a deeper look into the behavior of derived global quantities, like for instance the error in energy norm used here, can be useful and mathematically more rigorous.

5.2.2 Point load on the Scordelis-Lo roof

The next example is meant to demonstrate once again the higher accuracy per-degree-of-freedom achievable using local refinement. The geometrical setup and physical parameters are taken as defined in the Scordelis-Lo benchmark and are given in Figure 12(a), where we change the boundary conditions. Indeed, as external loading, we apply a point load of magnitude $10^{5}$ in the middle of the structure, directed in the negative vertical direction.

The convergence behavior of the displacement under the point load is depicted in Figure 12(b), where a reference value of $-0.206794699852788\ldots$ has been obtained with an overkill solution, which was computed using B-splines of degree $p=4$ on a fine uniform mesh with $\simeq 200\,000$ degrees-of-freedom. It can be seen that the solution obtained with an adaptive simulation based on the bubble estimator is several orders of magnitude more accurate compared to performing uniform refinement of the patch, for the same number of dofs. Finally, in Figure 13, the obtained displacement in the z-direction, hierarchical mesh and Von Mises stress distribution are presented at different steps of the adaptive loop where we remark that once again the estimator properly captures and resolves the sharp features of the solution.

6 Conclusions

We have introduced a novel a posteriori error estimator in the context of linear fourth-order partial differential equations, namely we have successfully applied it to adaptive simulation, based on the use of hierarchical B-splines, of isogeometric Kirchhoff plates and Kirchhoff-Love shells. This can be seen as an extension of similar approaches present in the literature, e.g. we refer to [6, 49]. The evaluation of the estimator is based on the solution of an additional, residual-like variational problem formulated in the so-called bubble space, which is composed of Bernstein polynomials defined locally on active elements. We remark that, thanks to this choice of the aforementioned space, the resulting linear system is in general small, block-diagonal and easily-invertible. Moreover, this method is suitable for parallelization and straightforward to implement into existing isogeometric codes. More importantly, it is computationally cheap compared to classical residual-type estimators since it avoids the computation of the residual in a strong sense. On one hand, this is a major advantage particularly for Kirchhoff-Love shells, since the evaluation of covariant derivatives is a tedious task, becomes very quickly computationally expensive and, for all practical purposes, almost intractable from a numerical standpoint. On the other hand, with this technique we also avoid the computation of integral terms involving the evaluation of the jump of the derivatives across element boundaries. Finally, we have also observed that in all our experiments the proposed estimator acts as an excellent approximation of the true error.

To conclude, we have numerically demonstrated the applicability and robustness of the proposed error estimator to steer an adaptive simulation of a wide range of problems of engineering relevance, where superior efficiency and accuracy per-degree-of-freedom have been achieved thanks to the local refinement capabilities of hierarchical B-splines.

Acknowledgements

The authors would like to thank Prof.Alessandro Reali and Dr.Rafael Vázquez for the fruitful discussions on the subject of this paper and Dr.Rafael Vázquez for his help with GeoPDEs. The authors also gratefully acknowledge the support of the European Research Council, via the ERC AdG project CHANGE n.694515.