Stabilization techniques for immersogeometric analysis of plate and shell problems in explicit dynamics

TL;DR

This paper develops stabilization techniques for immersogeometric analysis of plates and shells in explicit dynamics, enabling stable simulations with lumped mass matrices despite challenges posed by high-order PDEs and cut elements.

Contribution

The authors extend previous methods to achieve stable immersogeometric analysis of plates and shells with lumped mass matrices, using polynomial extensions to improve accuracy.

Findings

Stable immersogeometric analysis achieved for plates and shells

Lumped mass matrices can be used without spurious oscillations

Accuracy comparable to boundary-fitted discretizations

Abstract

Finite element plate and shell formulations are ubiquitous in structural analysis for modeling all kinds of slender structures, both for static and dynamic analyses. The latter are particularly challenging as the high order nature of the underlying partial differential equations and the slenderness of the structures all impose a stringent constraint on the critical time step in explicit dynamics. Unfortunately, badly cut elements in immersed finite element discretizations further aggravate the issue. While lumping the mass matrix often increases the critical time step, it might also trigger spurious oscillations in the approximate solution thereby compromising the numerical solution. In this article, we extend our previous work in \cite{voet2025stabilization} to allow stable immersogeometric analysis of plate and shell problems with lumped mass matrices. This technique is based on polynomial extensions and restores a level of accuracy comparable to boundary-fitted discretizations.

Equations

$$\Delta t\le \Delta t_c=\frac{C}{{\omega }_n}$$

$${\mathbf{a}}_{\alpha }={\mathcal{F}}_{,\alpha }$$

$${\mathbf{a}}_3=\frac{{\mathbf{a}}_1\times {\mathbf{a}}_2}{\Vert {\mathbf{a}}_1\times {\mathbf{a}}_2\Vert }$$

$${\mathbf{a}}^{\alpha }\cdot {\mathbf{a}}_{\beta }={\delta }_{\alpha \beta }$$

$$[a^{\alpha \beta }]=[a_{\alpha \beta }{]}^{-1}$$

$${\mathbf{a}}^{\alpha }=a^{\alpha \beta }{\mathbf{a}}_{\beta }.$$

$$b_{\alpha \beta }={\mathbf{a}}_{\alpha ,\beta }\cdot {\mathbf{a}}_3.$$

$$\text{I}(\mathbf{u},\mathbf{v})=a_{\alpha \beta }{u}^{\alpha }{v}^{\beta }\text{and}\text{mbox}II(\mathbf{u},\mathbf{v})=b_{\alpha \beta }{u}^{\alpha }{v}^{\beta }$$

$$\kappa (\mathbf{u})=\frac{\text{mbox}II(\mathbf{u},\mathbf{u})}{\text{I}(\mathbf{u},\mathbf{u})}.$$

$${\kappa }_{\min }\le \kappa (\mathbf{u})\le {\kappa }_{\max }$$

$$R_{\min }=\frac{1}{\max |\kappa (\mathbf{u})|}.$$

$$\Omega =S\times T.$$

$$\mathcal{G}({\xi }_1,{\xi }_2,{\xi }_3)=\mathcal{F}({\xi }_1,{\xi }_2)+{\xi }_3{\mathbf{a}}_3({\xi }_1,{\xi }_2),$$

$${\mathbf{g}}_{\alpha }$$

$${\mathbf{g}}_3$$

$$\eta =\frac{R_{\min }}{\tau }\text{and}\zeta =\frac{L}{\tau }$$

$${\int }_{\Omega }\rho \ddot{\mathbf{u}}\cdot \mathbf{v}\sqrt{g}\mathrm{d}\Omega +{\int }_{\Omega }\mathbf{ϵ}(\mathbf{u}):\mathbf{C}:\mathbf{ϵ}(\mathbf{v})\sqrt{g}\mathrm{d}\Omega ={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}\sqrt{g}\mathrm{d}\Omega +{\int }_{\partial {\Omega }_N}\mathbf{h}\cdot \mathbf{v}\sqrt{g}\mathrm{d}\partial \Omega \forall \mathbf{v}\in V$$

$$C^{ijkl}=\lambda g^{ij}{g}^{kl}+\mu (g^{ik}{g}^{jl}+g^{il}{g}^{jk}),$$

$$\mathbf{u}({\xi }_1,{\xi }_2,{\xi }_3,t)=\mathbf{u}({\xi }_1,{\xi }_2,t)+{\xi }_3\mathbf{\theta }({\xi }_1,{\xi }_2,t)$$

$$[\mathbf{u}{]}_{\mathcal{E}}=(u_1,u_2,u_3)\text{and}[\mathbf{\theta }{]}_{\mathcal{A}}=({\theta }_1,{\theta }_2).$$

$${ϵ}_{ij}=\frac{1}{2}({\mathbf{g}}_i\cdot {\mathbf{u}}_{,j}+{\mathbf{g}}_j\cdot {\mathbf{u}}_{,i}).$$

$${\mathbf{u}}_{,\alpha }$$

$${\mathbf{u}}_{,3}$$

$${ϵ}_{\alpha \beta }$$

$${ϵ}_{3\alpha }$$

$${ϵ}_{33}$$

$$2{\epsilon }_{\alpha \beta }$$

$$2{\kappa }_{\alpha \beta }$$

$$2{\chi }_{\alpha \beta }$$

$${\gamma }_{\alpha }$$

$$0=({\mathbf{a}}_3\cdot \mathbf{\theta }{)}_{,\alpha }$$

$${\theta }_{\alpha ,\beta }=({\mathbf{a}}_{\alpha }\cdot \mathbf{\theta }{)}_{\beta }$$

$$\mathsf{V}=\{\mathsf{v}:[\mathsf{v}{]}_{\mathsf{E}}\in [H^1{]}^3\times [H^1{]}^2,\mathsf{v}{|}_{\partial S_D}=0\}$$

$$b(\mathsf{u},\mathsf{v})+a(\mathsf{u},\mathsf{v})=F(\mathsf{v})\forall \mathsf{v}\in \mathsf{V}$$

Full text

Keywords: Plate and shell formulations, Trimming, Isogeometric analysis, Explicit dynamics, Mass lumping.

1 Introduction

Plate and shell structural elements are practically encountered in every modern structure, ranging from buildings, bridges, ships to aircraft fuselage. Engineering design commonly relies on advanced finite element software built from classical plate and shell mathematical models that were established more than a century ago. The intrinsically small thickness characterizing plates and shells often allows simplifying general 3D elasticity equations by assuming the displacement and/or stress fields satisfy special properties. Therefore, complex three-dimensional bodies are essentially described through their mid-surface as two-dimensional ones. Contrary to standard elasticity, the models derived either feature different primary unknowns (e.g., displacements and rotations) or high order derivatives that naturally impose greater regularity on the solution field. Unfortunately, building smooth approximation spaces from classical $C^{0}$ finite elements is rather cumbersome and despite many attempts, none were fully satisfactory [2]. Isogeometric analysis, on the contrary, naturally fulfills the smoothness requirement and is ideally suited for high order formulations.

Isogeometric analysis (IGA) is a spline-based finite element method, which was originally designed with the intention of bridging the gap between design and analysis by using a single model for both [3, 4]. The method came along with many other advantages. Most notably, spline functions, which are the primary technology in Computer-Aided-Design (CAD), may completely eliminate the geometry discretization error that characterizes standard finite element discretizations. Moreover, smooth spline functions have better approximation properties than classical $C^{0}$ Lagrange functions [5, 6, 7], which translates into greater accuracy per degree of freedom. Finally, there exist positive, compactly supported and stable bases, ideal for numerical computations. Therefore, isogeometric analysis is an outstanding candidate for high order finite element discretizations and held its promises, as testified by its numerous successful applications [8, 9, 10], particularly for plate and shell problems [11, 12].

However, for transient processes as in structural dynamics, a time discretization must complement the spatial one. In principle, isogeometric analysis could be used for both and is actively being investigated as an instance of space-time discretizations [13]. In practice though, practitioners still mostly resort to standard explicit time integration techniques based on finite differences. One of the reasons for favoring explicit methods over implicit ones stems from the physical constraint on the step size. Indeed, highly oscillatory (nonlinear) processes already inherently limit the step size, independently of the stability constraint. However, other less oscillatory processes in the realm of linear structural dynamics could benefit from greater step sizes and the stability constraint might then get in the way. This stability constraint, commonly referred to as the Courant–Friedrichs–Lewy (CFL) condition, is inversely proportional to the largest discrete frequency of the system ($\omega_{n}$) and states that

[TABLE]

where $\Delta t_{c}$ is the critical time step and $C$ is a constant that depends on the actual method. Thus, inaccurate high-frequency optical branches characterizing standard finite element discretizations often impose a stringent constraint on the step size. Although maximally smooth isogeometric analysis removes most of the optical branches, some outlier frequencies persist [8] and grow larger as the mesh gets finer and the degree increases [14]. Dampening the optical branches and removing outliers is a central topic in finite element and isogeometric analysis. While algebraic mass scaling techniques have been developed for standard finite element methods [15, 16, 17], isogeometric analysis (and tensorized finite element discretizations more generally) enjoy significantly more structure, which often forms the cornerstone of tailored outlier removal strategies. In this context, common techniques include spline reparameterizations [8], weak [18] or strong [19, 20] imposition of additional boundary conditions (also related to $n$-width optimal spline spaces) and eigenvalue deflation techniques [21]. Most of these methods are both applicable to second and fourth order problems.

Another great reason for employing explicit methods comes from mass lumping. Indeed, explicit time integration of undamped dynamical systems leads to solving a linear system with the mass matrix at each time step. Instead of solving them “exactly”, practitioners commonly substitute the (consistent) mass matrix with an easily invertible (preferably diagonal) matrix, a process that became known as mass lumping. Its most successful instance is undoubtedly for spectral element discretizations that consist in choosing as finite element nodes the Gauss-Lobatto quadrature nodes. Integrating the mass matrix with the same quadrature rule then naturally produces accurate diagonal lumped mass matrices, a method nowadays commonly known as nodal quadrature [22, 23]. In this context, it may also coincide with the row-sum [2] and diagonal scaling (or HRZ) techniques [24], two other well-known algebraic lumping strategies [25]. Moreover, mass lumping is often praised for increasing the critical time step, a property formally proved for the row-sum and extended to classes of (block) banded matrices in isogeometric analysis [26, 21]. Unfortunately, the nodal quadrature method is built around interpolatory bases and tensor product elements. Choosing a suitable lumping technique for more general bases and elements is far less obvious: while the row-sum technique may produce indefinite lumped mass matrices, the diagonal scaling method is often rather inaccurate [25, 27]. Although the positivity of the B-spline basis in isogeometric analysis is often seen as an advantage for guaranteeing positive definite lumped mass matrices, the row-sum technique in this context is still only second order accurate, independently of the spline order [8].

Intensive research is currently underway to improve its accuracy and several paths are being explored. Most notably, several authors have suggested applying the row-sum technique in a Petrov-Galerkin framework by choosing (approximate) dual functions as test functions. Although the idea, whose origins date back to Cottrell [4], initially did not attract much attention, it was taken up again recently in [28] with promising results. Since then, there has been a surge of interest in approximate dual functions [29, 30, 31, 32]. Unfortunately, such functions come with all sorts of difficulties. Firstly, (approximate) duality only holds in the parametric domain. Although the rescaling suggested in [29] ensures this property also holds in the physical domain, it complicates the assembly of the stiffness matrix. Secondly, as for any discretization involving dual functions, the treatment of essential boundary conditions requires ad hoc techniques [30]. Finally, the method has only been tested in rather simple academic settings and its extension to trimmed geometries is not in sight.

In another line of research, other authors [33, 34] have built interpolatory spline bases and combined them with the row-sum technique in an attempt to mimic the nodal quadrature method. The method is seamlessly applicable to wave type equations [33] as well as plate models [34] but although it improves the convergence rate, it remains sub-optimal. Moreover, several critical issues remain unattended. In particular, interpolatory spline functions are globally supported, leading to dense system matrices unsuited for large scale simulations. Additionally, quadrature rules with positive weights are essential for constructing positive definite lumped mass matrices. Yet, the Greville quadrature, as proposed in [33, 34], does not fulfill this condition and as a matter of fact, negative weights are encountered for non-uniform knot spans or spline spaces of intermediate regularity [35]. This last issue completely undermines the applicability of the method and has been overlooked by the authors.

Due to the aforementioned shortcomings, the standard row-sum technique often remains the method of choice for complex industrial applications, particularly for trimmed isogeometric geometries [36]. Trimming is ubiquitous in industrial CAD models and consists in discarding or merging portions of the domain delimited by spline curves or surfaces [37, 38]. Although it grants virtually unlimited flexibility in geometric design, small trimmed elements often cause great difficulties in analysis, including instabilities [39], ill-conditioning [40] and arbitrarily small step sizes (for consistent mass formulations) [41, 42]. The latter prohibits explicit time integration altogether unless some form of stabilization is employed [41] or alternative time integration schemes are sought [43]. Mass lumping is sometimes a remarkably simple workaround, which has been widely promoted in the engineering community [44, 45]. Indeed, the critical time step in this case may no longer depend on how badly trimmed elements are, provided the discretization is sufficiently smooth [36]. This surprising finding was later confirmed in several independent studies [44, 46, 45, 47] and recently unraveled theoretically [41, 42]. Nevertheless, there is another side to the story: while the largest discrete eigenvalues constraining the critical time step may remain bounded, the smallest ones instead converge to zero and bring in spurious high frequency modes in the low-frequency spectrum [47, 46]. As recently shown in [1], if these modes are unluckily activated, they may completely ruin the solution by triggering spurious high-frequency oscillations. As shown later in this contribution, those issues are not confined to wave type problems and may also arise for plate and shell models. Fortunately, a solution for the former may also be a solution for the latter and we will confirm it in this work. Tentative mitigation strategies were already listed in [47], including incomplete lumping strategies. However, although they might restore some accuracy, conveniently inverting the partially lumped mass matrix is less obvious, particularly in dimension $d\geq 2$, when trimmed elements may constitute a significant portion of the boundary. Instead, in [1], the authors stabilized the approximation space by extending polynomial segments of large elements into small neighboring elements. This technique both removed trimming-related outlier frequencies from the consistent mass and prevented spurious modes from surfacing in the low-frequency spectrum when it is lumped. After stabilization, the accuracy and CFL condition were in all aspects similar to boundary-fitted discretizations. Moreover, the stabilization method may be combined with advanced mass lumping techniques [26, 21] for further improving the accuracy of the row-sum. In this work, we extend the stabilization technique to plate and shell problems to produce similar results.

The outline for the rest of our article is as follows: After reviewing in Section 2 some elementary principles from differential geometry and shell theory, the variational problem is discretized in Section 3 with immersed isogeometric methods. The stabilization technique is then recalled in Section 4 and numerical experiments follow in Section 5. Although plate and shell models are prone to spurious oscillations for lumped mass formulations, they are cured with the same stabilization technique as for standard wave equations. Finally, Section 6 summarizes our findings and concludes our article.

2 Shell theory

In this section, we first recall some fundamental notions of differential geometry at the core of shell theories before formulating the problem in its weak (or variational) form. Throughout this section, we adopt some common conventions of shell theory whereby Greek indices take values $1,2$ and Latin indices takes values $1,2,3$. Additionally, we employ in this section (and this section only) Einstein’s convention where summations are carried out over repeated indices. Several key references have guided the construction of this section, including [48, 49, 50, 51], which the reader may consult for further details.

2.1 Shell geometry

The deformation of shells is always described with respect to their mid-surface, defined by a parameterized surface; i.e. a smooth injective map (or chart) $\bm{\mathcal{F}}\in C^{2}(S;\ \mathbb{R}^{3}$), which maps a point $\bm{\xi}=(\xi_{1},\xi_{2})$ in the parametric domain $S\subset\mathbb{R}^{2}$ to a point $\bm{x}=\bm{\mathcal{F}}(\bm{\xi})$ in the physical space $\mathbb{R}^{3}$. The parametric domain $S$ is assumed open, bounded and connected with a Lipschitz boundary $\partial S$. The mid-surface is then defined as $\widehat{S}=\bm{\mathcal{F}}(S)$, as shown in Figure 2.1.

For complex shell structures, several maps may be necessary for defining the mid-surface but for the sake of the presentation, we only consider a single one. To each point $\bm{x}=\bm{\mathcal{F}}(\xi_{1},\xi_{2})$ on the mid-surface is associated a local basis (called the covariant basis) whose vectors are defined as

[TABLE]

where $\bm{\mathcal{F}}_{,\alpha}$ denotes differentiation of $\bm{\mathcal{F}}$ with respect to $\xi_{\alpha}$. The covariant basis vectors $\bm{a}_{\alpha}$ are tangent to the mid-surface and span the tangent plane (see Figure 2.1.b). A unit normal vector (or director) is computed from the covariant basis as

[TABLE]

where $\times$ denotes the cross product and $\|\bullet\|$ the standard Euclidean norm. One of the most prominent objects in differential geometry is the metric tensor (or first fundamental form), which is nothing but a standard Gram matrix for the Euclidean inner product (denoted $\cdot$) on the tangent space. Its components in the covariant basis are $a_{\alpha\beta}=\bm{a}_{\alpha}\cdot\bm{a}_{\beta}$ and we let $a=\det(\nabla\bm{\mathcal{F}}^{T}\nabla\bm{\mathcal{F}})$ denote its determinant.

Shell structures may also experience deformations relative to its mid-surface. For this reason, one often introduces another more convenient basis, called the contravariant basis. Its basis vectors $\bm{a}^{\alpha}$ are dual to the covariant basis; i.e. they satisfy

[TABLE]

where $\delta_{\alpha\beta}$ is the Kronecker delta symbol ($\delta_{\alpha\beta}=1$ if $\alpha=\beta$ and [math] otherwise). The components of the metric tensor in the covariant and contravariant bases are related by

[TABLE]

and allow computing the vectors of the local contravariant basis as

[TABLE]

Shell theory distinguishes thick and thin shells based on their curvature and thickness. The components of the curvature tensor (or second fundamental form) are defined as

[TABLE]

The first and second fundamental forms induce the symmetric bilinear forms

[TABLE]

where $\bm{u}=u^{\alpha}\bm{a}_{\alpha}$ and $\bm{v}=v^{\beta}\bm{a}_{\beta}$ are vectors in the tangent plane and $u^{\alpha}$ and $v^{\beta}$ are their contravariant components. Note that the components in the covariant basis are actually called “contravariant”, while the components in the contravariant basis are called “covariant”. From the first and second fundamental forms, we define the curvature along the cross section normal to $\bm{u}$ by the Rayleigh quotient

[TABLE]

Since I and II are symmetric and I is positive definite,

[TABLE]

where $\kappa_{\min}$ and $\kappa_{\max}$ are called the principal curvatures. The curvature may take positive or negative values but in shell theory, one is mostly interested in its largest magnitude or equivalently, in the smallest radius of curvature

[TABLE]

The radius of curvature is one of the key quantities for defining the slenderness of shell structures. The other important geometrical quantity is the shell’s thickness $\tau$ that now shapes it as a fully three-dimensional body. Denoting $T=\left(-\frac{\tau}{2},\ \frac{\tau}{2}\right)$, we define the 3D parametric domain

[TABLE]

In principle, the shell’s thickness may vary from one point to another (i.e., $\tau=\tau(\xi_{1},\xi_{2})$) but for simplicity we do not consider such cases here. The reference (or initial) configuration of the full three-dimensional body of the shell is then defined as $\widehat{\Omega}=\bm{\mathcal{G}}(\Omega)$, where $\bm{\mathcal{G}}\colon\Omega\to\mathbb{R}^{3}$ is the 3D chart

[TABLE]

also sometimes called the canonical extension of $\bm{\mathcal{F}}$. Similarly to the surface, the covariant basis vectors for the volume $\bm{g}_{i}=\bm{\mathcal{G}}_{,i}$ are given by

[TABLE]

and the corresponding contravariant basis vectors $\bm{g}^{i}$ and determinant $g=\det(\nabla\bm{\mathcal{G}}^{T}\nabla\bm{\mathcal{G}})$ are defined analogously as before. Given some characteristic length $L$ for the in-plane dimension of $\Omega$, the ratios

[TABLE]

are called slenderness ratios and play an important role in (heuristically) classifying shell or plate theories. The former is commonly used for shells while the latter is customary for plates. To ensure all quantities are well-defined, one must assume that $\eta>0.5$ [50, Fig. 2.6]. Depending on how large $\eta$ or $\zeta$ is, various assumptions on the strain and/or stress fields may apply and underpin different theories. Since plates are actually a particular case of shells, we will focus on the latter, which is much more general. Among the most common and well-studied theories are the Kirchhoff-Love and Reissner-Mindlin theories. The Kirchhoff-Love theory assumes that the vectors normal to the mid-surface remain normal after deformation and therefore neglects transverse shear deformations. This assumption applies to thin shells (e.g., $\eta\geq 20$ [11]) and relates the deformation within the shell to the deformation of the mid-surface only. Although the model only features displacement unknowns, it also involves high order derivatives and imposes greater regularity requirements. As a matter of fact, finite element discretizations based on the Kirchhoff-Love theory require $C^{1}$ approximation spaces, a condition that is not easily fulfilled with classical (Lagrange) finite elements. This was historically one of the main reasons for adopting alternative shell theories in finite element codes [11].

One very popular alternative is the Reissner-Mindlin theory, where transverse shear deformations are properly accounted for and is applicable to moderately thick shells (e.g., $\eta<20$ [11]). However, the model requires introducing additional unknowns (rotations) for relating the displacement of the mid-surface to the displacement through the thickness. Those additional unknowns avoid introducing high order derivatives, which has tangible consequences: contrary to the Kirchhoff-Love theory, finite element discretizations based on the Reissner-Mindlin theory only require classical $C^{0}$ approximation spaces and are therefore more widespread [11]. For this reason, we adopt the Reissner-Mindlin theory in the rest of the article, while keeping in mind that the methods presented herein virtually apply to any shell formulation, including Kirchhoff-Love theory.

2.2 Shell kinematics

The variational formulation of shell problems originally stems from simplifications of the variational formulation of three-dimensional elasticity, albeit stated in curvilinear coordinates. Curvilinear coordinates are very natural in this setting and allow separating the in-plane contributions from the out-of-plane (or transverse) ones, which is essential for reducing the three-dimensional problem to a simpler two-dimensional one along the mid-surface. We merely state below the result of this transformation and refer interested readers to [49, Chapter 1] for its derivation. In curvilinear coordinates, the three-dimensional variational form reads: find $\mathbf{u}(t)\in V=\{\mathbf{v}\in[H^{1}]^{3}\colon\mathbf{v}|_{\partial\Omega_{D}}=0\}$ such that

[TABLE]

where the Dirichlet and Neumann boundaries ($\partial\Omega_{D}$ and $\partial\Omega_{N}$, respectively) form a disjoint partition of the boundary ($\partial\Omega=\partial\Omega_{D}\cup\partial\Omega_{N}$ and $\partial\Omega_{D}\cap\partial\Omega_{N}=\emptyset$) and we have assumed homogeneous Dirichlet boundary conditions for simplicity. In (2.7), $\rho$ is the density, $\bm{\epsilon}$ is the linearized strain tensor (defined later), $\mathbf{f}$ and $\mathbf{h}$ are the volume force and “surface traction”, respectively, and $\mathbf{v}$ is a test function. We must stress that the “surface traction” $\mathbf{h}$ differs from the expression one is accustomed to in standard elasticity and accounts for additional terms stemming from the change of variables. Its expression is derived in [49, Chapter 1] but is not needed in the sequel. Finally, $\bm{C}$ is the fourth order elasticity tensor, whose contravariant components for a linear elastic homogeneous and isotropic material are given by

[TABLE]

where $\lambda$ and $\mu$ are the Lamé constants. In this contribution, we limit ourselves to isotropic materials that cover the majority of industrial applications. For the constitutive relationship of laminates, which are also relevant in the context of shell structures, the interested reader may consult [52, 53].

Shell theories assume that the displacement field takes a specific form, which is another cause for diverging theories. Nevertheless, all shell theories share one common goal: when combined with ad hoc approximations, they reduce the problem to a two-dimensional one along the mid-surface. Our work focuses on classical shell theories that assume a unit vector normal to the mid-surface remains straight and unstretched after deformation. This assumption translates into expressing the displacement field $\mathbf{u}$ as

[TABLE]

where $\bm{u}(\xi_{1},\xi_{2},t)$ and $\bm{\theta}(\xi_{1},\xi_{2},t)$ denote the displacement and rotation, respectively, at a point $\bm{x}=\bm{\mathcal{F}}(\xi_{1},\xi_{2})$ of the mid-surface at time $t$. We have deliberately chosen the same letter for the displacement field in the shell and along the mid-surface but different font styles. This convention is systematically employed in analogous expressions.

The Reissner-Mindlin kinematical assumption [50, p. 95] states that $\bm{a}_{3}\cdot\bm{\theta}=0$, which implies that $\bm{\theta}$ lies in the tangent plane (i.e., it does not have any component along $\bm{a}_{3}$). Contrary to the Kirchhoff-Love theory, where the rotations are related to the displacements (derivatives), the Reissner-Mindlin theory features both displacements and rotations as primary unknowns. The displacements pertain to the mid-surface while the rotations relate the displacement of the mid-surface to the displacement through the thickness. Using the covariant coordinates $\theta_{\alpha}$ of the rotation vector is advantageous in this context since it reduces the number of unknowns. However, the lack of assumptions on the mid-surface’s displacement field does not suggest a particularly attractive basis and we simply stick to the Cartesian basis; i.e., $\bm{u}=u_{i}\bm{e}_{i}$ and $\bm{\theta}=\theta_{\alpha}\bm{a}^{\alpha}$. For a given vector field $\bm{v}$ and a fixed basis $\mathcal{B}$, we will denote $[\bm{v}]_{\mathcal{B}}$ the component vector of $\bm{v}$ in the basis $\mathcal{B}$. Hence, denoting $\mathcal{E}=\{\bm{e}_{1},\;\bm{e}_{2},\;\bm{e}_{3}\}$ the Cartesian basis and $\mathcal{A}=\{\bm{a}^{1},\;\bm{a}^{2}\}$ the surface’s contravariant basis,

[TABLE]

Remark 2.1**.**

Multiple other forms of the displacement field have been suggested in the literature. One rather natural generalization of (2.8) consists in adding high order terms in $\xi_{3}$ to better model the stress distribution through the thickness. The interested reader may consult [54, 55, 56] for an example of high order shell theories. There exist numerous other shell theories, which arise from combining various geometrical and mechanical assumptions.

As we will see, the special form of the displacement field (2.8) allows to reduce the general three-dimensional elasticity problem (2.7) to a simpler two-dimensional one. For understanding how, we must simplify the expression of the linearized strain tensor, whose covariant components are given by

[TABLE]

The derivation of the foregoing expression (rarely detailed in most articles) is deferred to Appendix A. Given the derivatives of the displacement field

[TABLE]

and the covariant basis vectors (2.5), the covariant components of the strain tensor (2.9) reduce to

[TABLE]

where

[TABLE]

These expressions (called generalized strains) merely follow from substituting (2.5) and (2.10) into (2.9) and exploiting the fact that

[TABLE]

to simplify and rewrite the relations. In particular, combining the special form of the displacement field with the kinematic assumption already enforces the place strain assumption ($\epsilon_{33}=0$). The conflicting plane stress assumption $(\sigma_{33}=0)$ is later enforced in the constitutive equations (see Appendix B). The quadratic dependency on $\xi_{3}$ in (2.11a) also later justifies neglecting $\chi_{\alpha\beta}$ for thin shells. The remaining quantities $\varepsilon_{\alpha\beta}$, $\kappa_{\alpha\beta}$ and $\gamma_{\alpha}$ are called membrane, bending and shear strains, respectively. Their actual expression is not very important in the context of this work; one must only note that they are independent of the transverse variable $\xi_{3}$.

2.3 Two-dimensional variational problem

The two-dimensional variational problem along the mid-surface stems from its full three-dimensional counterpart (2.7), after substituting the expressions for the displacement and strain fields (2.8) and (2.11) and combining them with several geometrical approximations. The final result is reported in this section while its derivation is postponed to Appendix B. We have not found a satisfactory reference presenting a sufficiently concise and comprehensive treatment of the Reissner-Mindlin model. For completeness, we have summarized and regrouped the main steps in the appendix and included pointers to relevant literature for a more detailed discussion.

Denoting $\bm{\mathsf{u}}=(\bm{u},\bm{\theta})$, $\bm{\mathsf{v}}=(\bm{v},\bm{\phi})$ and $\mathsf{E}=\mathcal{E}\times\{\bm{0}\}\cup\{\bm{0}\}\times\mathcal{A}$ the basis, we seek $\bm{\mathsf{u}}(t)\in\mathsf{V}$, where

[TABLE]

such that

[TABLE]

where

[TABLE]

where the Dirichlet and Neumann boundaries ($\partial S_{D}$ and $\partial S_{N}$, respectively) form a disjoint partition of the shell’s boundary ($\partial S=\partial S_{D}\cup\partial S_{N}$ and $\partial S_{D}\cap\partial S_{N}=\emptyset$). In (2.14b), $W=\operatorname{\operatorname{diag}}(\rho_{u}I_{3},\;\rho_{\theta}I_{3})$ with $\rho_{u}=\tau\rho$ and $\rho_{\theta}=\frac{\tau^{3}}{12}\rho$ and in (2.14c)

[TABLE]

are the shell’s “external forces” and “boundary tractions”, whose components are defined in Appendix B. The tensors $\bm{N}$, $\bm{M}$ and $\bm{Q}$ appearing in (2.14a), called the generalized membrane, bending and shear stresses, respectively, are defined as

[TABLE]

where $\alpha_{s}=\nicefrac{{5}}{{6}}$ is the shear correction factor and the fourth order tensors $E_{u}^{\alpha\beta\gamma\delta}=\tau E^{\alpha\beta\gamma\delta}$ and $E_{\theta}^{\alpha\beta\gamma\delta}=\frac{\tau^{3}}{12}E^{\alpha\beta\gamma\delta}$ are related to the shell elasticity tensor, whose contravariant components are

[TABLE]

Verifying the well-posedness of this reduced variational problem falls outside the scope of this contribution. The interested reader may nevertheless consult [49] for useful hints on the matter.

3 Isogeometric discretization and mass lumping

In finite element analysis, the components of an approximate solution to the variational problem (2.13) are sought in a finite dimensional subspace. In isogeometric analysis, this subspace is a spline space; i.e., a space of piecewise polynomials $\mathcal{S}_{h}$, parameterized by a mesh size $h$, spline order $p$ and continuity $0\leq k\leq p-1$. In principle, different components of the solution may be approximated in different spaces and the discretization parameters of a given space may depend on the parametric direction but for simplifying the notation, we assume that the spaces are identical for each component and that the discretization parameters are uniform. Extending the analysis to mixed approximation spaces is straightforward, albeit more cumbersome. Furthermore, the solution of time-dependent PDEs is often approximated in two stages, by first discretizing in space and then in time. The approximate solution is called semi-discrete when it is only discretized in space while it is called fully discrete when it is discretized in both space and time. Although space-time discretization techniques (simultaneously discretizing both in space and in time) are gaining momentum, the former is still the most popular technique and the one employed in this work. Hence, denoting $m=\dim(\mathcal{S}_{h})$ the dimension of the spline space and $\bm{\xi}=(\xi_{1},\xi_{2})$ the two-dimensional curvilinear coordinates, for each time $t$, $\bm{\mathsf{u}}(t)\approx\bm{\mathsf{u}}_{h}(t)=(\bm{u}_{h}(t),\bm{\theta}_{h}(t))$ where $\bm{\mathsf{u}}_{h}(t)\in\mathsf{V}_{h}=\{\bm{\mathsf{v}}_{h}\colon[\bm{\mathsf{v}}_{h}]_{\mathsf{E}}\in[\mathcal{S}_{h}]^{3}\times[\mathcal{S}_{h}]^{2},\,\bm{\mathsf{v}}_{h}|_{\partial S_{D}}=\bm{0}\}$ such that

[TABLE]

where $B_{j}(\bm{\xi})$ are the B-spline basis functions, $\bm{u}_{j}(t)\in\mathbb{R}^{3}$ and $\bm{\theta}_{j}(t)\in\mathbb{R}^{2}$ are unknown coefficient vectors and $Q=[\bm{a}^{1},\;\bm{a}^{2}]\in\mathbb{R}^{3\times 2}$ is the matrix whose columns are the (surface) contravariant basis vectors. This matrix accounts for the fact that $\bm{\theta}_{h}$ is expressed in the contravariant basis. For trimmed isogeometric shells, the B-spline basis is typically constructed over an entire fictitious (non-trimmed) parametric domain $S_{0}\supseteq S$ before retaining the functions whose support intersects the parametric domain $S$, defining the so-called active basis. This procedure is exemplarily illustrated in Figure 3.1 for the parameterized surface of Figure 2.1. The construction of the B-spline basis over the fictitious domain $S_{0}$ is fairly standard and the reader may consult [3, 4] for the details.

A standard Galerkin approximation of the variational problem now leads to solving the following second order system of ordinary differential equations (ODEs) [2, 57]

[TABLE]

where $\bm{\textsf{u}}(t)$ is the unknown coefficient vector containing all the displacement $\bm{u}_{j}(t)$ and rotation $\bm{\theta}_{j}(t)$ degrees of freedom, sorted componentwise. Hence, the mass matrix $M$ and stiffness matrix $K$ are $5\times 5$ block matrices. While the stiffness matrix $K$ is generally fully block-coupled and its expression is rather complicated (and not needed here), the mass matrix $M$ is only partially coupled and is of much greater interest in the context of this work. Indeed, since $W$ is diagonal in (2.14b), $M$ is actually block diagonal. Its structure can be conveniently described through the direct sum of matrices. Given square matrices $A_{i}$, their direct sum is the block diagonal matrix

[TABLE]

With this definition, the consistent mass matrix is expressed as $M=M_{u}\oplus M_{\theta}$, where

[TABLE]

are the mass matrices related to the displacement and rotation degrees of freedom, respectively. Their blocks $\mathrm{M}_{u}$ and $\mathrm{M}^{\alpha\beta}_{\theta}$ (denoted with straight letters) are defined as

[TABLE]

Clearly, $M$ (and its diagonal blocks $M_{u}$ and $M_{\theta}$) are symmetric positive definite. Those properties still hold in more general cases with non-diagonal matrices $W$, although $M$ would then have some block-coupling between displacement and rotation degrees of freedom, just as the stiffness matrix. However, although the B-spline basis functions are nonnegative, the mass matrix (and more specifically the off-diagonal blocks of $M_{\theta}$) may feature negative entries, due to the coefficients $a^{\alpha\beta}$ in (3.2). This property sharply contracts with classical elasticity or plate models and might jeopardize standard mass lumping techniques. Indeed, when explicit time-stepping schemes are employed for approximately solving (3.1), the consistent mass $M$ is oftentimes substituted with an ad hoc (diagonal) lumped mass matrix. Various lumping techniques are known but a very popular choice is the row-sum technique [2]. For a given $M\in\mathbb{R}^{n\times n}$, this technique is defined algebraically as

[TABLE]

where $d_{i}=\sum_{j=1}^{n}m_{ij}$ for $i=1,\dots,n$. Practitioners directly substitute in (3.1) the consistent mass $M$ by its lumped mass approximation $\mathcal{L}(M)$, which in our context becomes

[TABLE]

However, since $M_{\theta}$ is not nonnegative, the lumped mass matrix $\mathcal{L}(M_{\theta})$ might end up being indefinite or even singular, both of which preclude simulations in explicit dynamics. In the next lemma, we present a sufficient condition guaranteeing positive definite lumped mass matrices. This condition depends on the angle between the covariant basis vectors, denoted $\theta(\bm{a}_{1},\bm{a}_{2})$.

Lemma 3.1**.**

If

[TABLE]

then $\mathcal{L}(M)$ is positive definite.

Proof.

First note that

[TABLE]

where $\mathrm{L}^{\alpha}=\sum_{\beta=1}^{2}\mathcal{L}(\mathrm{M}^{\alpha\beta}_{\theta})$. Since $\mathrm{M}_{u}$ is nonnegative, our only concern is ensuring that $\mathrm{L}^{\alpha}$ is positive definite. Yet,

[TABLE]

Condition (3.4) guarantees that $\sum_{\beta=1}^{2}a^{\alpha\beta}>0$ and therefore $(\mathrm{L}^{\alpha})_{ii}>0$. Indeed,

[TABLE]

Thus, recalling (2.3),

[TABLE]

since $a$ (the determinant of the surface metric tensor) is positive. Consequently, $\mathcal{L}(M)$ is positive definite. ∎

Remark 3.2**.**

Although sufficient, the condition stated in Lemma 3.1 is not necessary. Indeed, the integral in (3.5) may be positive even though the integrand becomes negative pointwise. Anyway, the condition is not very restrictive and essentially means that the norm of the covariant basis vectors should be balanced and the angle between them should not be too small. In the sequel, we will assume this condition is always satisfied.

Although more advanced mass lumping strategies have been developed for isogeometric analysis [28, 29, 33], most of them are only applicable to boundary-fitted discretizations. To the best of our knowledge, the only methods supporting single-patch, multi-patch as well as trimmed geometries, are based on algebraic generalizations of the row-sum technique that effectively enlarge its bandwidth [26, 21]. However, the resulting lumped mass matrices are no longer diagonal and we will stick to the classical row-sum technique for simplicity. Apart from facilitating the solution of linear systems, mass lumping is also praised for increasing the critical time step of explicit time integration schemes, which is inversely proportional to the largest discrete frequency of the system (see (1.1)). The constraint on the step size is often referred to as the Courant–Friedrichs–Lewy (CFL) condition.

For shell problems, mass lumping is often combined with a form of rotational mass scaling that artificially increases the inertia of rotational degrees of freedom [12, 17]. Although it increases the critical time step, it might also heavily degrade the accuracy of the solution. Thus, we have not incorporated any form of mass scaling in this work.

Mass lumping has always played a central role in explicit dynamics and became even more prominent after Leidinger [36] showed that the largest generalized eigenvalues for the row-sum lumped mass remained bounded regardless of how badly cut the elements were, provided the spline space was sufficiently smooth. This finding has since then been thoroughly investigated, both numerically [44, 45, 46, 47] and more theoretically [41, 42] and is by now relatively well understood. It turns out this “blessing of smoothness” is more of a double-edged sword. Indeed, while the largest eigenvalues may no longer depend on the cut elements, the smallest ones do and converge to zero at a rate that depends on the polynomial degree and cut configuration [42]. These eigenvalues bring in spurious oscillatory modes in the low-frequency spectrum [46, 47] and if those modes are unluckily activated, they may cause spurious oscillations in the solution itself. Their disastrous effect is evident from the recent counter-examples reported in [1] for 1D and 2D wave problems. Thus, mass lumping is not a miracle and merely trades one problem for another. Once spurious eigenvalues and modes are lying in the low-frequency spectrum, filtering them out is not so obvious. Instead, the most convenient way of eliminating them is to prevent them from surfacing in the first place.

Fortunately, spurious high-frequency modes for the consistent mass and spurious low-frequency modes for the lumped mass have exactly the same origin and any solution for the former is also a solution for the latter. Among them are stabilization techniques based on local polynomial extensions [39, 58, 1] and extended B-splines [59, 60, 61]. In [1], the authors used the former for stabilizing the discrete formulation and subsequently formed the row-sum lumped mass from the stabilized consistent mass. This technique removes the large spurious (trimming-related) eigenvalues for the consistent mass and prevents the corresponding modes from surfacing in the low-frequency spectrum when the mass matrix is lumped. Although this solution has only been studied for a simple wave equation, we believe it also applies to plate and shell models, where similar issues might arise. We will demonstrate it in Section 5 but before doing so, the stabilization technique is recalled in the next section.

4 Stabilization

In this section, we explain how to build a stable approximation space that is later used for approximating both the displacement and rotation fields. This introduction closely follows the presentation in [39, 58, 1], where a similar construction was employed for the wave equation and the Poisson problem.

We consider the B-spline basis $\mathcal{B}=\{B_{i}\}_{i=1}^{n}$ for a spline space in the trimmed parametric domain $S$. As explained earlier, for trimmed geometries, a standard tensor product basis is first constructed over a fictitious domain $S_{0}$ containing the parametric domain $S$ and $\mathcal{B}$ then collects all basis functions whose support intersects the parametric domain. Similarly, $\mathcal{T}_{h}$ collects all active mesh elements. Unfortunately, the support of some of the functions in $\mathcal{B}$ may only intersect the domain over arbitrarily small elements and those functions typically cause great stability and conditioning issues [40, 38]. The idea presented in [39] is to locally substitute those functions with the extension of polynomial segments from large neighboring elements. For a fixed tolerance $\gamma\in[0,1]$, an element $T\in\mathcal{T}_{h}$ is labeled as large (or good) if

[TABLE]

and is small (or bad) otherwise. The sets of large and small elements, denoted $\mathcal{T}_{h}^{L}$ and $\mathcal{T}_{h}^{S}$ respectively, partition $\mathcal{T}_{h}$ into disjoint sets such that $\mathcal{T}_{h}=\mathcal{T}_{h}^{L}\cup\mathcal{T}_{h}^{S}$. Given a small element $T\in\mathcal{T}_{h}^{S}$, there always exist a large neighboring element $T^{\prime}\in\mathcal{T}_{h}^{L}$ provided the boundary is sufficiently regular and the mesh sufficiently fine [62]. This large neighboring element may either be chosen such that it maximizes $|T^{\prime}\cap S|$ among all neighbors [39] or minimizes the distance between the centers of mass of $|T\cap S|$ and $|T^{\prime}\cap S|$ [58]. The B-splines on the small element $T$ are then locally substituted with their extension from the good neighbor $T^{\prime}$, as in Figure 4.1. In 1D (Figure 1(a)), this technique coincides with extended B-splines [60, 61]. However, in dimension $d\geq 2$ (Figure 1(b)), independent extensions may create discontinuities between neighboring elements. Consequently, the approximation space is generally non-conforming and is not a subspace of the original spline space, which precludes forming the system matrices from their non-stabilized counterpart. Instead, the assembly routine of element matrices must be modified, as described in [1, Algorithm 5.1]. Adapting the algorithm presented therein to other operators and right-hand sides is rather straightforward. Alternatively, one may restore conformity by combining the extension with a projection step [58], but this is not necessary from an approximation point of view. Another possible remedy are extended B-splines [60, 61], although the increased overlap of extended basis functions may slightly reduce the sparsity of system matrices. On the contrary, the local extension technique ensures the sparsity pattern of stabilized system matrices is contained in their non-stabilized counterpart [1, Lemma 5.2], which also allows the deployment of advanced mass lumping techniques [21], if need be.

In any case, the extension procedure generally does not preserve the positivity of the B-spline basis (see Figure 1(a)). Consequently, the stabilized consistent mass matrix may feature negative entries and again the row-sum technique no longer guarantees positive definite lumped mass matrices. Fortunately, the magnitude of such entries is typically very small unless large elements are stabilized. Thus, indefinite lumped mass matrices are very unlikely provided $\gamma$ is sufficiently small. In all our experiments, we chose $\gamma=0.1$ and never encountered any issue. More generally, values of $\gamma$ between $0.05$ and $0.15$ are equally good choices. For a more complete overview of polynomial extension techniques and implementation aspects, we refer to [59, 58, 1] and the references therein. It is important to realize that none of these techniques are restricted to smooth spline discretizations and also apply to classical $C^{0}$ finite elements.

5 Numerical experiments

In this section, we present a few numerical experiments showing on the one hand the alarming effect mass lumping may have on the numerical solution and on the other hand how the polynomial extension technique might resolve the issue. All our experiments are carried out in GeoPDEs [63], an open-source MATLAB/Octave software package for isogeometric analysis. Our experiments range from academic examples to more realistic industrial ones, featuring irregularly shaped cut elements. In all our experiments, we chose $\gamma=0.1$ for classifying large and small elements.

Example 5.1** (Plate externally trimmed).**

Our simplest example is a trimmed square plate of side length $L=1$ m (see Figure 5.1). In this case, the shell model of Section 2 reduces to a plate one. For illustrating how mass lumping may ruin the simulation, we consider a manufactured solution whose displacement and rotation fields are

[TABLE]

where the spatial part $w\colon S\to\mathbb{R}$ and temporal part $\phi\colon[0,\;+\infty)\to\mathbb{R}$ are defined as

[TABLE]

with the functions

[TABLE]

with parameter values $a_{0}=10$ cm, $b_{0}=0.3$, $n_{0}=5$, $\omega=\frac{1}{10L}\sqrt{\frac{E}{\rho}}$, unit Young modulus $E$ and density $\rho$ and Poisson ratio $\nu=0.25$.

Neumann (or type 2) boundary conditions are prescribed over the entire boundary and are computed from the exact solution, just as the source term and initial conditions. Snapshots of the vertical component of the manufactured displacement field $(u_{3})$ are shown in Figure 5.2.e1-e3 for three different times $t_{1}$-$t_{3}$ as indicated in Figure 5.2.a. The plate is discretized with maximally smooth B-splines of degree $p$ over a background grid of mesh size $h$. The left and right (resp. top and bottom) boundaries are at a distance of $\delta$ from the nearest vertical (resp. horizontal) grid line. For future reference, we denote $\epsilon=\delta/h$ the relative trimming parameter. The numerical solutions are computed with (1) the consistent mass, (2) the lumped mass and (3) the stabilized lumped mass for various spline orders $p$, trimming parameters $\epsilon$ and thicknesses $\tau=\{1,\;5,\;10\}$ cm. For the stabilized lumped mass, the approximation space is first stabilized with the polynomial extension technique recalled in Section 4 before lumping the mass matrix with the standard row-sum technique. Stabilizing the mass matrix prior to mass lumping filters out the spurious (trimming-related) outlier frequencies and modes, thereby preventing them from moving toward the low-frequency part of the spectrum when the mass matrix is lumped. Note that computing the numerical solution for the (non-stabilized) consistent mass with an explicit scheme is practically infeasible, due to the stringent constraint on the step size. Hence, we have used the Central Difference method (an explicit Newmark method) for the stabilized and non-stabilized lumped mass solutions, while relying on an implicit unconditionally stable Newmark method for the consistent mass solution, which was computed with the same step size as the explicit solutions to ensure a comparable error in time.

Snapshots of the solutions are provided in Figure 5.2 for the consistent (c1-c3), lumped (l1-l3) and stabilized lumped (s1-s3) solutions at three different times $t_{1}$-$t_{3}$. As shown in Figure 5.2.s1-s3, the stabilization technique drastically improves the accuracy of the lumped mass solution. Indeed, while the solution for the consistent mass (Figure 5.2.c1-c3) closely resembles the exact solution (Figure 5.2.e1-e3), the solution for the lumped mass is utterly inaccurate and features spurious oscillations nears the corners of the plate (Figure 5.2.l1-l3), even when the exact solution vanishes. As a matter of fact, the numerical solution in this case is completely worthless and further deteriorates as the spline order and slenderness increase or as the trimming parameter decreases (see Figure 5.3). The error is much more pronounced in the $L^{\infty}$ norm that measures pointwise inaccuracies (Figure 5.4). Although refining the mesh (i.e., decreasing $h/L$) improves the lumped mass solution (see Figure 5.4), an acceptable accuracy might only be achieved for infeasibly small mesh sizes. Therefore, refining the mesh might not be a viable solution due to computer limitations. Remarkably, stabilizing the space prior to lumping the mass completely eliminates the oscillations. As shown in Figure 5.5, the error for the stabilized lumped mass now perfectly overlaps with the error for the consistent mass. We must stress that the improvement is not attributed to an improvement of the lumping technique per se, since both the stabilized and non-stabilized lumped mass solutions are computed with the row-sum technique. However, its accuracy is space (and even basis) dependent, which again hints at modifying the space prior to lumping the mass.

Simulations in explicit dynamics are a fine balance of accuracy and efficiency. On the one hand, the consistent mass oftentimes yields the best accuracy one can hope for, but is prohibitively expensive. On the other hand, the lumped mass drastically improves the efficiency if one is willing to give up on some accuracy. Admirably, the stabilized lumped mass retains the efficiency of the (non-stabilized) lumped mass while improving its accuracy. Indeed, as shown in Figure 5.6, the critical time step of neither the lumped mass nor its stabilized counterpart depend on the trimming parameter, contrary to the (non-stabilized) consistent mass. For the sake of completeness, we have also included results for the stabilized consistent mass. While stabilizing the discrete formulation eliminates the deleterious effect trimming has on the spectrum, all other dependencies on the discretization parameters remain unchanged, particularly for the largest eigenvalue indirectly shown in Figure 5.6. Indeed, as for boundary-fitted discretizations, the critical time step still decreases as the spline order increases or as the thickness and the mesh size decrease.

Example 5.2** (Rotated plate).**

As we will show, relatively simple examples are not exempt from developing spurious oscillations. A square plate of thickness $\tau=5$ cm initially centered at the origin is trimmed along its bottom and right edges before being rotated $45\degree$ counterclockwise (see Figure 5.7). Neumann (or type 2) boundary conditions are prescribed along its two trimmed edges while Dirichlet (or type 1) boundary conditions are enforced on the remaining two edges. Contrary to the previous example, we prescribe boundary data and initial conditions instead of computing them from a manufactured solution. The shell’s “external forces” $\bm{\mathsf{f}}=(\bm{f},\bm{m})$ are

[TABLE]

and the “boundary tractions” $\bm{\mathsf{h}}=(\bm{h},\bm{n})$ are

[TABLE]

where $\bm{r}$ denotes the outward unit normal vector and $\phi(t)=\sin(\omega t)$. The parameter values are $f_{0}=0.1$ $\mathrm{N/m^{2}}$, $m_{0}=0.1$ $\mathrm{N/m}$, $h_{0}=0.1$ $\mathrm{N/m^{2}}$, $b_{0}=6$, $n_{0}=28$ and $\omega=\frac{1}{10L}\sqrt{\frac{E}{\rho}}$ with unit elastic modulus $E$ and density $\rho$ and Poisson ratio $\nu=0.25$. Finally, the initial conditions are

[TABLE]

where $v_{0}=7.74\times 10^{-5}$ $\mathrm{m/s}$. The background mesh discretizing the geometry aligns with the plate’s boundaries, whose trimmed edges are only $\delta$ away from the nearest grid lines. The ratio $\epsilon=\delta/h$, where $h$ denotes the mesh size, measures the severity of trimming. In this experiment, we have set $\epsilon=10^{-9}$ and $L=1$ m. Figure 5.8.e1-e2 shows a reference solution computed over a fine mesh at two distinct times $t_{2}$ and $t_{3}$, as indicated in Figure 5.8.a. The numerical solutions are computed with exactly the same strategy as in the previous example. Namely, the Central Difference method for the lumped mass solutions and an unconditionally stable Newmark method for the consistent mass with the same step size. As shown in Figure 5.8.c1-c2 and Figure 5.8.s1-s2, the solutions for the consistent and stabilized lumped mass closely resemble the reference solution. On the contrary, the solution for the (non-stabilized) lumped mass once again features spurious oscillations that are particularly pronounced at the corner of the two trimmed edges. This example proves that these oscillations can occur even for relatively simple expressions of the prescribed data and, therefore, they could a priori arise in practical industrial applications. As shown in Figure 5.8.s1-s2, stabilizing the discrete space prior to lumping the mass is a very convenient way of ruling out this possibility.

Example 5.3** (Plate with cut-out).**

This third example slightly modifies the first one by creating a square cut-out in the middle of the plate of side length $2a=0.4$ m (Figure 5.9). Moreover, the inner boundary is trimmed instead of the outer one. The manufactured solution is modified accordingly as

[TABLE]

where the spatial part $w\colon S\to\mathbb{R}$ and temporal part $\phi\colon[0,\;+\infty)\to\mathbb{R}$ are defined as

[TABLE]

with the functions

[TABLE]

and

[TABLE]

with parameter values $a_{0}=10$ cm, $n_{0}=4$, $n=6$, Poisson ratio $\nu=0.25$ and $\omega=\frac{1}{2L}\sqrt{\frac{E}{\rho}}$ with unit Young modulus $E$ and density $\rho$.

Dirichlet (or type 1) boundary conditions are prescribed along the outer boundary while Neumann (or type 2) boundary conditions are prescribed along the inner boundary. The exact solution (5.1) is depicted in Figure 5.10.e1-e3 at three distinct times $t_{1}$-$t_{3}$. The plate is discretized with maximally smooth B-splines over a grid of mesh size $h$. The inner boundaries are $\delta$ away from the nearest grid lines and $\epsilon=\delta/h$ again denotes the relative trimming parameter. The solutions for the consistent mass and (stabilized) lumped mass, computed with the same methods as in the first example, are shown in Figure 5.10 for times $t_{1}$-$t_{3}$. Shortly after the beginning of the simulation (at time t1), all solutions look reasonably accurate. However, oscillations gradually develop in the solution for the lumped mass and amplify over time. Nevertheless, we do not expect these oscillations to indefinitely grow over time. Based on the results in [1], the error should still behave periodically, albeit with an amplitude and period that depends on both $\epsilon$ and $p$. Indeed, the pattern shown in Figure 5.11 corroborates this finding. For relatively large values of $\epsilon$, the error behaves periodically and we expect it still holds over larger time scales as $\epsilon$ becomes smaller. Moreover, by comparing Figure 5.11.a1-a3 to Figure 5.11.b1-b3, increasing the degree exacerbates this trend. As explained in [1], this finding correlates with the behavior of the smallest discrete eigenvalues, which were analyzed in [42]. As shown in Figure 5.12, stabilizing the discrete formulation (nearly) fixes the issue. However, contrary to the first example, the error for the stabilized lumped mass does not perfectly overlap with the consistent mass for thick plates. We believe this is due to the inherent inaccuracy of the row-sum technique (independently of trimming) and advanced mass lumping strategies [21] might eliminate or mitigate this discrepancy.

Example 5.4** (Fuselage window).**

Our fourth example is one step closer to industrial type examples and models an aircraft fuselage window. This is our first encounter with an actual shell structure, whose mid-surface is defined by the parameterization

[TABLE]

where $R=1$ m is the radius of curvature. A window, consisting in a rectangular opening with rounded corners, is cut out from a fuselage panel of side length $L=1$ m (Figure 5.13.b). The width and height of the opening are $2a=0.3$ m and $2b=0.4$ m, respectively, and its corners form circular arcs of radius $r=0.08$ m (Figure 5.13.a). As in the previous examples, we consider a manufactured solution

[TABLE]

where the spatial part $w\colon S\to\mathbb{R}$ and temporal part $\phi\colon[0,\;+\infty)\to\mathbb{R}$ are defined as

[TABLE]

with the functions

[TABLE]

and

[TABLE]

with parameter values $a_{0}=10$ cm, $\beta=10$, $n_{0}=3$, $n=6$, Poisson ratio $\nu=0.25$ and $\omega=\frac{1}{10L}\sqrt{\frac{E}{\rho}}$ with unit material parameters.

Dirichlet (or type 1) boundary conditions are prescribed along the outer sides of the panel while Neumann (or type 2) boundary conditions are prescribed all along the window’s boundaries. The boundary data (as well as the source term and initial conditions) are computed from the expressions in (5.3). The exact solution is shown in Figure 5.16.e1-e3 for times t1-t3. The problem is discretized with maximally smooth B-splines of order $p$ over a grid of mesh size $h$. The vertical sides of the window are at a distance $\delta$ from the nearest grid line and as usual $\epsilon=\delta/h$ denotes the relative trimming parameter. The numerical solutions for the consistent, lumped and stabilized lumped solutions are computed with the same schemes as before (i.e., the Central Difference method is used for both solutions with a lumped mass and an unconditionally stable Newmark method is used for the solution with the consistent mass). Snapshots of the solutions are shown in Figure 5.14 for the consistent (c1-c3), lumped (l1-l3) and stabilized lumped (s1-s3) solutions. Once again, spurious oscillations are visible along the trimmed boundary for the lumped mass solution, while the other two solutions are visually indistinguishable from the exact one. Those oscillations are not confined to the displacement field and may also affect the stress field, as shown in Figure 5.15.c for the first component of the bending stress ($M^{11}$) at the final time t3. Even without an exact or reference solution, the solution for the lumped mass would raise suspicions as it contradicts physical intuition. However, the error is not in the implementation but in the lumping itself. The error for the consistent and lumped mass solutions is reported in Figure 5.16 and follows a trend in all aspects similar to previous examples. While the accuracy of the solution for the consistent mass is not affected by small trimmed elements, the solution for the lumped mass rapidly deteriorates as $\epsilon$ becomes smaller or $p$ becomes larger (note the different range of values for $\epsilon$ in Figure 5.16.a1-a3 and Figure 5.16.b1-b3). Increasing the slenderness also negatively affects the solution, although to a milder extent. As expected, stabilizing the solution completely resolves most of these issues (Figure 5.17). However, the quality of the solution still depends on the slenderness, as it does for boundary-fitted discretizations. Hence, the stabilization only resolves trimming-related issues and we may envision combining it with different techniques for solving common shell-related issues (e.g., locking).

Example 5.5** (B-Pillar).**

Our final example, inspired from [64], treats the B-pillar of a car. This vertical part, located between the front and rear doors, is an essential component to the vehicle’s safety. The irregularly shaped object of length $L=1.4$ m, Young modulus $E=210\cdot 10^{9}$ Pa and Poisson ratio $\nu=0.3$ is clamped on one side in all directions (type 1 boundary conditions) and is subjected to a displacement field $u_{3}=-0.025$ m on its opposite side. The other components of the displacement field on this same side are left unconstrained by prescribing Neumann boundary conditions instead. We refer to such hybrid cases as type 3 boundary conditions. Finally, to add additional complexity to the model, an L-shaped cut-out is created on the base of the B-pillar and Neumann (or type 2) boundary conditions are prescribed along its inner boundaries. All geometric details and boundary conditions are depicted in Figure 5.18.a and a close-up of the cut-out is shown in Figure 5.18.b together with the underlying mesh. Once again, some of the grid lines are exceedingly close to the inner boundaries, only at a distance $\delta=5\cdot 10^{-8}$. External forces and moments are prescribed on the triangular portion of the domain between the L-shaped cut-out the dashed blue segment $P_{B}-P_{A}=(0.09,\;0.09)$ m (see Figure 5.18.b). We denote $\bm{\eta}$ the downward pointing normal to this segment and $\bm{r}$ the usual outward pointing normal to the boundary. The shell’s “external forces” $\bm{\mathsf{f}}=(\bm{f},\bm{m})$ are

[TABLE]

and the “boundary tractions” $\bm{\mathsf{h}}=(\bm{h},\bm{n})$ are

[TABLE]

where $\eta$ denotes the coordinate along the normal vector $\bm{\eta}$ and $\phi(t)=\sin(\omega t)$. The parameter values are $f_{0}=1.7\cdot 10^{7}$ Pa/m2, $m_{0}=1.7\cdot 10^{7}$ Pa/m, $h_{0}=1.7\cdot 10^{7}$ Pa/m, $b_{0}=0.1$, $n_{0}=90$ and $\omega=10^{4}$ rad/s. The initial velocity conditions in the local triangular portion are

[TABLE]

where $v_{0}=5\cdot 10^{2}$ m/s. Finally, $\bm{u}(\bm{x},0)$ and $\bm{\theta}(\bm{x},0)$ are given by the solution of the static problem that only accounts for the boundary conditions of type 1 and 3. With such settings, this example aims at simulating a local perturbation to a structural component otherwise in equilibrium. The exact solution to this problem in unknown and we must content ourselves with a reference solution computed over a fine mesh. The numerical solutions for the consistent and (stabilized) lumped mass are computed with the same schemes as before. At a first glance, all solutions reported in Figure 5.19 at final time appear satisfactory. The magnitude of the displacement field is maximal at the top of B-pillar subjected to the prescribed displacement and gradually decreases to zero, matching the homogeneous Dirichlet boundary conditions at its base. However, a closer inspection reveals oscillations near one of the corners of the cut-out, which are confirmed by the close-up figures in the bottom row of Figure 5.19. Although the results are not as clear-cut as in the previous examples, the oscillations for the (non-stabilized) lumped mass are somewhat sporadic and unexpected (Figure 5.19.l2). Such patterns do not appear in the solutions for the consistent and stabilized lumped mass (Figure 5.19.c2 and Figure 5.19.s2, respectively), suggesting they are a pure artifact of mass lumping and trimming. This finding provides yet another example of the deleterious effects mass lumping may have for trimmed geometries, even in complex industrial-like settings.

6 Conclusion

Finite element shell models are ubiquitous in structural engineering and, as for other structural models, they commonly rely on mass lumping techniques for fast explicit solutions. However, for trimmed isogeometric discretizations, mass lumping may introduce spurious frequencies and modes in the low-frequency spectrum. These modes are not an issue as long as the structural response is orthogonal to them. However, if they are unluckily activated, they may trigger spurious oscillations in the solution thereby compromising simulation results. This fact was recently uncovered for standard wave equations in [1] and the extensive numerical experiments within this article indicate that the same issue arises for plate and shell models. Since our counter-examples were obtained not only for manufactured solutions but also for prescribed data, we cannot a priori exclude real engineering applications from suffering the same fate. Unfortunately, once lodged in the low-frequency spectrum, filtering out spurious eigenvalues and modes is difficult. Hence, we suggest modifying the discrete formulation prior to lumping the mass matrix. The modification consists in replacing polynomial segments of basis functions on small trimmed elements by extending those from large neighboring elements. This stabilization technique prevents trimming-related spurious modes from surfacing in the low-frequency spectrum when the mass matrix is lumped. As shown in our experiments, the stabilization restores in most cases a level of accuracy comparable to boundary-fitted discretizations. Nevertheless, the stabilization only resolves trimming-related issues and combining it with advanced mass lumping techniques or locking-free strategies might be necessary for further improving the accuracy of the solution.

Acknowledgments

The authors kindly acknowledge the support of the Swiss National Science Foundation through the project “PDE tools for analysis-aware geometry processing in simulation science” n. 200021 215099, and through the project “Fast simulation tools for lattice structures” n. 200021 214987.

Appendix A Covariant components of the linearized strain tensor

The covariant components $\epsilon_{ij}$ of the strain tensor are its components in the contravariant basis; i.e.

[TABLE]

where $\{\bm{g}^{1},\;\bm{g}^{2},\;\bm{g}^{3}\}$ is the contravariant basis pertaining to the volume mapping. Thanks to the Kronecker delta property relating the covariant and contravariant bases, the $(i,j)$th component of the strain tensor is given by

[TABLE]

All we need for computing it is an expression for the strain tensor in curvilinear coordinates, which we compute from its expression in Cartesian coordinates. Indeed, the displacement field as a function of the Cartesian and curvilinear coordinates, denoted $\hat{\mathbf{u}}$ and $\mathbf{u}$, respectively, are related through the volume mapping $\bm{\mathcal{G}}$

[TABLE]

Thus, given the familiar expression for the strain tensor in Cartesian coordinates

[TABLE]

we obtain its expression in curvilinear coordinates as

[TABLE]

where the subscript for the gradient indicates the differentiation variable. Finally, recalling that the columns of $\nabla_{\bm{\xi}}\bm{\mathcal{G}}=[\bm{g}_{1},\;\bm{g}_{2},\;\bm{g}_{3}]$ are the covariant basis vectors, we obtain

[TABLE]

which is the expression in (2.9).

Appendix B From volume to surface integration

In order to ensure that the article is self-contained, this appendix outlines how the general three-dimensional variational problem integrated over the volume may be reduced to a simpler two-dimensional one integrated over the mid-surface, as is customary for shell formulations [65, 66]. For the Reissner-Mindlin model, this process is not so obvious and requires several ad hoc assumptions and approximations, which are somewhat scattered in the literature [48, 49, 50, 51]. Hence, we summarize the main steps in this appendix.

We recall that the transverse coordinate $\xi_{3}$ varies in the interval $T=\left(-\frac{\tau}{2},\ \frac{\tau}{2}\right)$ and note that the Neumann boundary is composed of three parts: $\partial\Omega_{N}=\partial\Omega_{N}^{-}\cup\partial\Omega_{N}^{+}\cup\partial\Omega_{N}^{0}$ (see Figure B.1) where

[TABLE]

The assumption on the displacement field (2.8) allows to considerably simplify the general three-dimensional variational problem: find $\mathbf{u}(t)\in V=\{\mathbf{v}\in[H^{1}]^{3}\colon\mathbf{v}|_{\partial\Omega_{D}}=0\}$ such that

[TABLE]

The steps for reducing (B.1) to an analogous (but approximate) two-dimensional problem are rather nontrivial. The main difficulty lies in isolating the dependency on $\xi_{3}$ from the other two curvilinear coordinates in order to integrate along the thickness (assumed constant). Our derivation is quite informal and follows an approach similar to [49, Chapter 3] whereby high order terms in $\xi_{3}$ are dropped out. We will first investigate how $g$ (the determinant of the volume metric tensor) is approximated before analyzing each term of (B.1) individually from left to right.

B.1 Volume metric

Let us begin with the determinant of the volume metric tensor $g=\det(\nabla\bm{\mathcal{G}}^{T}\nabla\bm{\mathcal{G}})=\det(\nabla\bm{\mathcal{G}})^{2}$. From the multilinearity of the determinant and the expression for the covariant basis vectors (2.5), we obtain

[TABLE]

Moreover, since $\bm{a}_{i}\cdot\bm{a}_{3}=\delta_{i3}$,

[TABLE]

Hence, a first order approximation is given by

[TABLE]

This approximation justifies substituting $\sqrt{g}$ with $\sqrt{a}$ in (B.1) for thin-walled structures.

B.2 Inertia forces

Now we simplify each term from left to right in (B.1) after replacing the trial and test functions by their expression according to the Reissner-Mindlin assumption

[TABLE]

The first term then becomes

[TABLE]

If $\rho$ is assumed independent of $\xi_{3}$, the above expression reduces to

[TABLE]

with $\bm{\mathsf{u}}=(\bm{u},\bm{\theta})$, $\bm{\mathsf{v}}=(\bm{v},\bm{\phi})$ and $W=\operatorname{\operatorname{diag}}(\rho_{u}I_{3},\;\rho_{\theta}I_{3})$ where

[TABLE]

Remark B.1**.**

If $\rho$ depends on $\xi_{3}$, $W$ is non-diagonal (but symmetric).

B.3 Internal forces

The derivations for the second term (related to the internal forces) is quite lengthy. We first express it componentwise as

[TABLE]

where $C^{ijkl}$ are the contravariant components of $\bm{C}$, which, for a linear elastic homogeneous and isotropic material, are given by

[TABLE]

and $\epsilon_{ij}$ are the covariant components of the strain tensor, listed in (2.11) and recalled below:

[TABLE]

The expressions for the membrane, bending and shear strains ($\bm{\varepsilon}$, $\bm{\kappa}$ and $\bm{\chi}$, respectively) are given in (2.12) but are not needed here; one must only bear in mind that they do not depend on $\xi_{3}$. Since the components of the elasticity tensor $\bm{C}$ do depend on $\xi_{3}$, we will first substitute them with a first order approximation. For this purpose, recall that the contravariant components of the volume metric tensor $g^{ij}$ are the components of $(\nabla\bm{\mathcal{G}}^{T}\nabla\bm{\mathcal{G}})^{-1}$. Let us decompose $\nabla\bm{\mathcal{G}}$ as

[TABLE]

where

[TABLE]

Since $\bm{A}$ is invertible, $\nabla\bm{\mathcal{G}}^{-1}=\bm{A}^{-1}(\bm{I}+\xi_{3}\bm{X}\bm{A}^{-1})^{-1}$ and if $\xi_{3}$ is small enough (such that $\|\xi_{3}\bm{X}\bm{A}^{-1}\|<1$ in any suitable matrix norm), a Neumann series expansion yields (see [67, Corollary 5.6.16] and the subsequent exercise therein)

[TABLE]

Hence, $\bm{A}^{-1}$ is a first order approximation of $\nabla\bm{\mathcal{G}}^{-1}$ and therefore

[TABLE]

The foregoing expression justifies approximating the contravariant components of the volume metric tensor as $g^{ij}\approx a^{ij}$. Hence,

[TABLE]

This approximation allows reducing the three-dimensional constitutive model to a two-dimensional one. For future reference, we note that $A$ inherits the major and minor symmetries of $C$ (i.e. $A^{ijkl}=A^{klij}$, $A^{ijkl}=A^{jikl}$ and $A^{ijkl}=A^{ijlk}$) and yields several useful identities:

[TABLE]

Imposing the plane stress assumption ($\sigma_{33}=0$) allows to eliminate $\epsilon_{33}$ from the equations and the Lamé constants in this context are given by

[TABLE]

Note that, within the Reissner-Mindlin theory, the plane stress assumption is incompatible with the condition $\epsilon_{33}=0$, which is a direct consequence of the kinematic hypothesis on the displacement field. However, the accuracy of the theory is not affected by this inconsistency. From the three-dimensional constitutive relation (with $C$ substituted for $A$) and the identities in (B.5),

[TABLE]

Now we simplify the constitutive relation by exploiting (B.5) and substituting $\epsilon_{33}$ by its expression in (B.6)

[TABLE]

We denote

[TABLE]

the shell elasticity tensor. Substituting $\epsilon_{\alpha\beta}$ with its expression in (B.4a), neglecting high order terms and integrating over the thickness yields

[TABLE]

where

[TABLE]

are called the generalized membrane, bending and shear stresses, respectively, and $E_{u}^{\alpha\beta\gamma\delta}=\tau E^{\alpha\beta\gamma\delta}$ and $E_{\theta}^{\alpha\beta\gamma\delta}=\frac{\tau^{3}}{12}E^{\alpha\beta\gamma\delta}$ are simply rescaled versions of the shell elasticity tensor. In practice, the contribution from shear stresses is scaled with an ad hoc shear correction factor $\alpha_{s}$ such that $Q^{\alpha}=\mu\alpha_{s}\tau a^{\alpha\beta}\gamma_{\beta}$. This factor ensures consistency with classical bending theory [2, Chapter 6]. Finally, putting all the pieces back together, we obtain

[TABLE]

B.4 Volume forces

For the volume forces (i.e., the third term), we obtain

[TABLE]

where $\bm{\mathsf{f}}_{1}=(\bm{f}_{1},\bm{m}_{1})$ with

[TABLE]

representing distributed forces and moments applied to the mid-surface $S$, respectively.

B.5 Surface tractions

Since $\partial\Omega_{N}=\partial\Omega_{N}^{-}\cup\partial\Omega_{N}^{+}\cup\partial\Omega_{N}^{0}$, the surface tractions give rise to three distinct contributions: those from the “bottom” $\partial\Omega_{N}^{-}=S\times\left\{-\frac{\tau}{2}\right\}$ and “top” $\partial\Omega_{N}^{+}=S\times\left\{\frac{\tau}{2}\right\}$ surfaces of the shell as well as from its “sides” $\partial\Omega_{N}^{0}=\partial S_{N}\times T$. Denoting $\mathbf{h}^{-}$, $\mathbf{h}^{+}$ and $\mathbf{h}^{0}$ these respective contributions, the fourth and last term in (B.1) becomes

[TABLE]

where $\bm{\mathsf{f}}_{2}=(\bm{f}_{2},\bm{m}_{2})$ and $\bm{\mathsf{h}}=(\bm{h},\bm{n})$ whose components are

[TABLE]

Finally, combining (B.3), (B.8), (B.9) and (B.10), we obtain the two-dimensional shell variational formulation

[TABLE]

where $\bm{\mathsf{f}}=\bm{\mathsf{f}}_{1}+\bm{\mathsf{f}}_{2}=(\bm{f},\bm{m})$.