A Higher-order Trace Finite Element Method for Shells
AA Higher-order Trace Finite Element Method for
Shells
D. Schöllhammer , T.P. Fries May 1, 2020 , Institute of Structural AnalysisGraz University of TechnologyLessingstr. 25/II, 8010 Graz, Austria [email protected] , [email protected] Abstract
A higher-order fictitious domain method (FDM) for Reissner-Mindlin shells isproposed which uses a three-dimensional background mesh for the discretization.The midsurface of the shell is immersed into the higher-order background mesh andthe geometry is implied by level-set functions. The mechanical model is based onthe Tangential Differential Calculus (TDC) which extends the classical models basedon curvilinear coordinates to implicit geometries. The shell model is described byPDEs on manifolds and the resulting FDM may typically be called Trace FEM. Thethree standard key aspects of FDMs have to be addressed in the Trace FEM aswell to allow for a higher-order accurate method: (i) numerical integration in thecut background elements, (ii) stabilization of awkward cut situations and eliminationof linear dependencies, and (iii) enforcement of boundary conditions using Nitsche’smethod. The numerical results confirm that higher-order accurate results are enabledby the proposed method provided that the solutions are sufficiently smooth.
Keywords:
Trace FEM; Fictitious domain methods; Tangential differential cal-culus; Shells; Manifolds; Implicit geometries; a r X i v : . [ c s . C E ] A p r CONTENTS
Contents
A shell is a curved, thin-walled structure and due to their high bearing capacity, they occurin many engineering applications such as automotive, aerospace, biomedical- and civil-engineering [12, 22, 58]. In the mechanical modelling of a shell, a dimensional reductionfrom the 3D shell body to its midsurface is a key step. The equilibrium equations resultinto partial differential equations (PDEs) on manifolds embedded in the physical space R . In the classical approach of modelling shells, the midsurface is typically defined by a(piecewise) parametrization which may be called an explicit definition. Then, there existsa map from some two-dimensional parameter space ˆΩ to R . In the context of the classicalSurface FEM, the discrete midsurface is rather defined by an atlas of local element-wisemappings from the reference element to the physical elements. An overview of classical shelltheory is given, e.g., in [15, 55, 56, 3] or in the text books [12, 4, 1, 57, 58]. Alternatively, theshell geometry may also be defined implicitly following the level-set method [25, 26, 47, 54].The main advantage of implicit geometry descriptions is that the application of recent finiteelement techniques such as fictitious domain methods (FDM) for shells is enabled.Herein, we propose a higher-order accurate fictitious domain method for implicitly definedReissner–Mindlin shells. The geometry of the shell is defined implicitly by means of level-set functions. The shell boundary value problem (BVP) is discretized with a higher-orderaccurate Trace FEM approach.The Trace FEM is a fictitious domain method for solving PDEs on manifolds, see e.g.,[48, 43, 44, 31, 29, 53, 5]. This may also be called Cut FEM which, in the last years,became a popular FDM enabling higher-order accuracy [8, 10, 9]. When using the CutFEM for the solution of PDEs on manifolds as done herein, the method becomes analogousto the Trace FEM [11, 13]. This recent finite element technique is fundamentally differentcompared to standard Surface FEM approaches, see e.g., [19, 21, 23, 26]. In the following,the major differences in, (1) geometry definition, (2) mesh and degrees of freedom (DOFs)and, (3) generation of integration points in these two approaches are outlined and visualizedfor the Trace FEM in Fig. 1 and for the Surface FEM in Fig. 2, respectively.1. Firstly, the differences in the geometry definition are emphasized. In Fig. 1(a), animplicitly defined shell by means of level-set functions is presented. In particular,the yellow surface indicates the midsurface of the shell and is defined by the zero-isosurface of a master level-set function. The boundary of the shell is defined by Introduction (a) implicit shell midsurface (b) cut elements (c) integration points on the zeroisosurface
Fig. 1:
Overview of the Trace FEM: (a) Implicit definition of the shell geometry by means of level-set functions, (b) the set of cut three-dimensional background elements are labelled asactive mesh and their nodes imply the DOFs for the numerical simulation, (c) integrationpoints on the zero-isosurface plotted in red for the domain and in blue on the boundary. additional slave level-set functions, i.e., the light green, purple and grey surfaces inFig. 1(a). In contrast to the implicit definition, an explicit representation of the shellmidsurface for this example is shown in Fig. 2(a). Then, the surface is given throughan atlas of element-wise local mappings, as usual in the context of Surface FEM,implying local curvilinear coordinates. Note that a parametrization is only availablein the explicit case and is, in general, not available in implicitly defined geometries.2. Secondly, the mesh generation and the location of DOFs are considered. In theTrace FEM the domain of interest, i.e., the shell midsurface, is embedded in a three-dimensional background mesh. The 3D background mesh may consist of higher-orderLagrange elements with the only requirement that the domain of interest is com-pletely immersed. Neither the shell surface nor the shell boundary have to conformto the background mesh. There is one master level-set function whose zero-isosurfaceimplies the shell surface and additional slave level-set functions imply the boundaries.The set of cut elements is labelled active mesh and is visualized in Fig. 1(b) wherecubic tetrahedral elements are used as an example. For the numerical simulation,the DOFs are located at the nodes of the active mesh, which are clearly not on themidsurface of the shell. The corresponding shape functions are those of the activemesh and are restricted to the trace. This means that for the integration of the weakform, the 3D shape functions are only evaluated on the zero-isosurface of the masterlevel-set function. In contrast, in the Surface FEM a boundary conforming surface (a) explicit shell midsurface (b) surface mesh (c) integration points on the sur-face emsh
Fig. 2:
Overview of the Surface FEM: (a) Explicitly defined shell obtained by an atlas of element-wise local maps, (b) conforming surface mesh consisting of cubic 2D Lagrange elements,(c) integration points for the Surface FEM obtained by standard Gauß integration rules. mesh, see Fig. 2(b), is defined through an atlas of element-wise local mappings andthe DOFs are located at nodes of the surface mesh which are on the discrete mid-surface of the shell. The corresponding shape functions are the 2D shape functionsliving only on the surface mesh. The location of the DOFs and the different dimen-sionality of the shape functions are the most important differences between the twofinite element techniques.3. Lastly, integration points need to be placed on the shell midsurface for the integrationof the weak form. In the case of the Trace FEM, this is not a trivial task especiallywhen a higher-order accurate integration scheme is desired [24, 25, 26, 27, 41, 39],further details are given below. Regarding the Surface FEM standard Gauß integra-tion rules are applicable and mapped from the reference to the surface elements inthe usual manner. In Fig. 1(c) and Fig. 2(c), integration points in the domain (red)and on the boundary (blue) are visualized for the Trace FEM and Surface FEM,respectively.As mentioned above, the Trace FEM approach is a fictitious domain method and the fol-lowing three well-known implementational aspects require special attention: (i) integrationof the weak form, (ii) stabilization and, (iii) enforcement of essential boundary conditions.Regarding the integration of the weak form, suitable integration points with higher-orderaccuracy for multiple level-set functions have to be provided. Herein, the approach, whichnaturally extends to multiple level-set functions as outlined by the authors in [24, 25, 26, 27]
Introduction is employed. Other higher-order integration schemes for implicitly defined surfaces with one level-set function are presented, e.g., in [41, 39]. A stabilization of the stiffness matrixis necessary due to small supports caused by unfavourable cut scenarios and the restrictionof the shape functions to the trace. An overview and analysis of the different stabilizationtechniques in the Trace FEM is presented in [44]. Herein, the “normal derivative volumestabilization” first introduced for scalar-valued problems in [30, 11] and for vector-valuedproblems in [31, 42] is used. The advantage of this particular stabilization technique isthat it is suitable for higher-order accuracy; a straight forward implementation and arather flexible choice of the stabilization parameter is possible. The essential boundaryconditions need to be enforced in a weak manner due to the fact that a strong enforcementby prescribing nodal values does not apply (because the nodes of the background meshare not on the shell boundary). Therefore, a similar approach as presented by the authorsin [27] is employed. In particular, the non-symmetric version of Nitsche’s method is used,see, e.g., [7, 49, 32, 27].The Trace FEM approach is applied in a wide range of applications. In particular, transportproblems are presented, e.g., in [20, 11, 5, 44], flow problems are considered, e.g., in [39,31, 6, 42, 35, 37], moving and evolving manifolds are detailed in [44, 46, 40]. The firstapplication of the Trace FEM to membranes has just recently been achieved for the linearmembrane in [13] and for large deformation membranes by the authors in [27]. This isconsiderably more challenging because in structural mechanics, the governing equationsof membranes and shells have been, until recently, only formulated based on curvilinearcoordinates. This, however, does not directly apply for an implicit shell definition and theuse of the Trace FEM. Therefore, the reformulation of classical (curvilinear) models forshells and membranes during the last years was a crucial preliminary step for the use of theTrace FEM. Most importantly, the use of the Tangential Differential Calculus (TDC) wasfound highly useful and generalizes the mechanical models in the sense that they becomevalid for explicit and implicit geometry definitions [33, 34, 50, 51, 53, 27]. In contrast, inflow and transport applications on curved surfaces, the general coordinate-free definitionof the boundary value problems is a standard for a long time [20, 19, 21, 36, 23], thusenabling the application of the Trace FEM earlier than in structural mechanics as proposedherein. Herein, to the best knowledge of the authors, this is the first time, that a TraceFEM approach is applied to curved Reissner–Mindlin shells. Furthermore, the proposedapproach is also higher-order accurate, thus enabling optimal higher-order convergencerates if the involved fields are sufficiently smooth. Other Trace FEM applications withparticular emphasis on higher-order accuracy are rather scarce. If the geometry is definedby multiple level-set functions we refer to [27] and for one level-set function we refer to[39, 35, 37].The Reissner–Mindlin shell model is suitable to model thin and moderately thick shellsand as mentioned above the employed shell model need to be applicable on implicitlydefined surfaces, where a parametrization of the midsurface is not existent. Therefore, aformulation of the shell equations in the frame of the TDC is employed herein. The recastof the linear shell equations in the frame of the TDC is presented by the authors in [51].The shell formulation is based on a standard difference vector approach. The tangentialityconstraint on the difference vector is accomplished by a projection of a full 3D vector ontothe tangent space combined with a consistent stabilization in normal direction.The paper is organized as follows: In Section 2, the implicit definition of shell geometries isdescribed. Furthermore, the employed differential surface operators and other importantquantities in the frame of the TDC are briefly introduced. In Section 3, the TDC-basedReissner–Mindlin shell model is outlined following [51, 53]. The resulting BVP is formu-lated in strong and weak form including boundary conditions. In Section 4, the higher-orderTrace FEM is introduced in detail. Furthermore, the implementational aspects of the TraceFEM, i.e., (i) integration of the weak form, (ii) stabilization and, (iii) enforcement of es-sential boundary conditions are elaborated. The discrete weak form of the equilibriumis introduced and the discretization of the difference vector is considered. In Section 5,numerical results of the proposed approach are presented with a set of benchmark exam-ples. The results confirm that optimal higher-order convergence rates are achieved whenthe solution is sufficiently smooth. The paper ends in Section 6 with a summary andconclusions.
A shell is a thin-walled, possibly curved structure with thickness t . For the modelling,the 3D shell body Ω may be reduced to its midsurface Γ embedded in the physical space R . In general, the midsurface of a shell can be defined explicitly , see e.g., [3, 1], or implicitly , see e.g., [50, 51, 26]. Herein, the shell geometry is implicitly defined by meansof (multiple) level-set functions. Let there be a master level-set function φ ( x ) : R → R whose zero-isosurface defines the (unbounded) midsurface of the shell in R , see Fig. 3(b). Preliminaries
The boundaries of the shell ∂ Γ i are defined by of additional slave level-set functions ψ i with i = 1 , . . . , n Slaves , see Fig. 3(c), where the blue lines indicates the boundary of the shell ∂ Γ .For the sake of simplicity, all slave level-set functions shall feature the same orientation,i.e., positive inside the domain and negative outside. The bounded midsurface Γ and theboundaries ∂ Γ i of the shell are then defined by Γ := (cid:8) φ ( x ) = 0 ∩ ψ i ( x ) > ∀ i : ∀ x ∈ R (cid:9) , (2.1) ∂ Γ i := (cid:8) φ ( x ) = 0 ∩ ψ i ( x ) = 0 ∀ x ∈ R (cid:9) , (2.2)where the union of all boundaries define ∂ Γ := ∪ i ∂ Γ i . Slave level-set functions are notnecessarily needed when the master-level set function is restricted to some bounded domainof definition Ω ∗ ⊂ R rather than R . Both variants for the implicit definition of thebounded shell geometry work equally well in the method proposed herein.The normal vector of the shell is given through the normalized gradient of the masterlevel-set function φ n Γ ( x ) = ∇ φ ( x ) (cid:107)∇ φ ( x ) (cid:107) . (2.3)Along the boundaries ∂ Γ i , there is an associated tangent vector t ∂ Γ ,i which is defined by thenormalized cross product of the corresponding slave level-set function ψ i and the masterlevel-set function φ t ∂ Γ ,i ( x ) = ψ i ( x ) × φ ( x ) (cid:107) ψ i ( x ) × φ ( x ) (cid:107) , x ∈ ∂ Γ i . (2.4)The orientation of the tangent vector t ∂ Γ ,i is defined with the orientation of the slave level-set function ψ i . Furthermore, the co-normal vector n ∂ Γ at ∂ Γ pointing “outwards” andbeing perpendicular to the boundary yet in the tangent plane T p Γ is n ∂ Γ ,i = t ∂ Γ ,i ( x ) × n Γ ( x ) , x ∈ ∂ Γ i , (2.5)see Fig. 3(d), where the normal vector n Γ and the local triad ( n Γ , n ∂ Γ ,i , t ∂ Γ ,i ) at theboundaries are visualized. .1 Tangential differential calculus (TDC) (a) isosurfaces of φ (b) implicit shell midsurface (c) some bounded shell midsurface n Γ n Γ n Γ n Γ n ∂ Γ , n ∂ Γ , t ∂ Γ , t ∂ Γ , n Γ n Γ n ∂ Γ , n ∂ Γ , t ∂ Γ , t ∂ Γ , (d) normal vector and local triads Fig. 3:
Implicit definition of a (bounded) spherical shell: φ ( x ) = (cid:107) x (cid:107) − r, ψ ( x ) = x, ψ ( x ) = z .(a) The colors are different isosurfaces of φ , (b) implicit shell midsurface defined by thezero-isosurface of φ , (c) definition of boundaries with additional slave level-set functions ψ (purple) and ψ (green), (d) normal (black), co-normal (red) and tangent (green)vectors. The tangential differential calculus provides (TDC) a framework to define surface operatorsindependently of the concrete surface definition. This is an important generalization com-pared to classical models based on curvilinear coordinates which rely on parametrizations,i.e., explicit geometry definitions. In the following, geometrical and differential operatorsamong other important quantities are defined in the frame of the TDC. For a more detailed0
Preliminaries information and derivation we refer to, e.g., [27, 51, 52, 18].
On the manifold Γ , there are two projection operators P ( x ) = I − n Γ ( x ) ⊗ n Γ ( x ) , (2.6) Q ( x ) = I − P ( x ) , (2.7)where P projects arbitrary vectors v (Γ) ∈ R into the tangent space T P Γ of Γ and Q projects an arbitrary vector in normal direction of Γ . The tangential gradient operator ∇ Γ of a differentiable scalar -valued function u : Γ → R on the surface is given by ∇ Γ u ( x ) = P ( x ) · ∇ ˜ u ( x ) , x ∈ Γ (2.8)where ∇ is the standard gradient operator, and ˜ u is an arbitrarily, but sufficiently smoothextension of u in a neighborhood U of the manifold Γ . In the context of the Trace FEM, thefunction ˜ u is naturally available because all quantities are defined in the higher-dimensionalspace, i.e., R or Ω ∗ . Therefore, ˜ u is scalar-valued function in R , i.e., ˜ u : R → R and thefunction u is defined by the restriction of ˜ u to Γ , i.e., u := ˜ u | Γ which means that ˜ u is onlyevaluated on Γ .When surface gradients of vector -valued functions u ( x ) : Γ → R are considered, it isimportant to distinguish directional and covariant gradients: ∇ dir Γ u ( x ) = ∇ dir Γ u ( x ) v ( x ) w ( x ) = ( ∇ Γ u ) T ( ∇ Γ v ) T ( ∇ Γ w ) T = ∇ ˜ u · P , ∇ covΓ u ( x ) = P · ∇ dir Γ u ( x ) = P · ∇ ˜ u · P . Concerning the surface divergence of vector-valued functions u ( x ) : Γ → R and tensor-1valued functions A ( x ) : Γ → R × , there holdsdiv Γ u ( x ) = div Γ ( u, v, w ) T = tr (cid:0) ∇ dir Γ u (cid:1) = tr ( ∇ cov Γ u ) =: ∇ Γ · u , div Γ A ( x ) = div Γ ( A , A , A ) div Γ ( A , A , A ) div Γ ( A , A , A ) =: ∇ Γ · A . The surface gradient of the normal vector n Γ yields the Weingarten map H [16, 36, 51], H = ∇ dir Γ n Γ = ∇ cov Γ n Γ , (2.9)with the property that the directional and covariant gradient are identical. Furthermore,the Weingarten map is a symmetric in-plane tensor, i.e., H = H T , H · n Γ = 0 and its twonon-zero eigenvalues are the principal curvatures κ , = − eig ( H ) . To derive the weak form of the governing equations, the following divergence theorem onmanifolds is needed [17, 18], (cid:90) Γ u · div Γ A d A = − (cid:90) Γ ∇ dir Γ u : A d A + (cid:90) Γ κ · u · A · n Γ d A + (cid:90) ∂ Γ u · A · n ∂ Γ d s, (2.10)where ∇ dir Γ u : A = tr (cid:0) ∇ dir Γ u · A T (cid:1) . For tangential (in-plane) tensor functions with A = P · A · P , the term involving the mean curvature κ = tr ( H ) vanishes and one finds ∇ dir Γ u : A = ∇ cov Γ u : A . The shell is implicitly defined by means of level-set functions. Therefore, it is a crucialrequirement that the shell theory extends to this situation and the obtained boundary valueproblem is well-defined also on implicitly defined geometries, where a parametrization ofthe midsurface does not exist. Herein, the linear Reissner–Mindlin shell theory is recasted2
Governing equations in the frame of the TDC as presented in [51, 53]. The main advantage is that the employedsurface operators are independent of the concrete surface definition. Therefore, the usageof this shell theory in the frame of the TDC is an appealing choice. In the following, theshell equations from [51] are briefly recalled.In [51], the displacement field u Ω ( x ) is decomposed into the displacement of the midsurface u and the rotation of the normal vector is modelled with a difference vector approach w , seeFig. 4. The total rotation of the normal vector is split into two parts: (1) due to bending ofthe midsurface, (2) and due to transverse shear deformations γ . The coordinate in thicknessdirection is labelled ζ and provided that the master level-set function is a signed-distancefunction, ζ = φ ( x ) . Note that the difference vector is a tangential vector. ζ n Γ P ( x Γ , ζ ) Undeformedzero-isosurface Γ ζζ ¯ n Γ uu Ω ¯ P ( x Γ , ζ ) ζ w − ζ ( ∇ dir Γ u ) T · n Γ ζ γ Deformedzero-isosurface ¯Γ yzx x Γ x Γ Fig. 4:
Displacement field u Ω of the Reissner-Mindlin shell. Based on the displacement field u Ω , one may obtain the strain tensor ε Γ by computing thesymmetric part of the surface gradient ∇ dir Γ u Ω . The strain tensor ε Γ is split into in-planestrains ε P Γ and transverse shear strains ε S Γ using the projectors from Section 2.1.1 ε Γ ( x ) = ε P Γ ( x ) + ε S Γ ( x ) , (3.1)with ε P Γ ( x ) = P · ε Γ · P . (3.2) = ε P Γ , Mem + ζ ε P Γ , Bend (3.3) ε S Γ ( x ) = Q · ε Γ + ε Γ · Q . (3.4) .1 Reissner–Mindlin shells ε P Γ , Mem and bendingstrains ε P Γ , Bend resulting in ε P Γ , Mem ( u ) = 12 (cid:2) ∇ cov Γ u + ( ∇ cov Γ u ) T (cid:3) , (3.5) ε P Γ , Bend ( u , w ) = 12 (cid:2) H · ∇ dir Γ u + ( ∇ dir Γ u ) T · H + ∇ cov Γ w + ( ∇ cov Γ w ) T (cid:3) , (3.6) ε S Γ ( u , w ) = 12 (cid:2) Q · ∇ dir Γ u + ( ∇ dir Γ u ) T · Q + n Γ ⊗ w + w ⊗ n Γ (cid:3) . (3.7)For simplicity, we shall use a linear elastic material governed by Hooke’s law with themodified Lamé constants µ = E ν ) , λ = Eν − ν in order to eliminate the normal stress inthickness direction. Based on that, the linear stress tensor yields σ Γ ( x ) = 2 µ ε Γ ( x ) + λ tr [ ε Γ ( x )] I . (3.8)A decomposition of the stress tensor in a similar manner than the strain tensor gives themembrane, bending and transverse shear stresses. With the assumption of a constantshifter in thickness direction, an analytical pre-integration w.r.t. the thickness is possibleand the stress resultants, such as moment tensor m Γ , effective normal force tensor ˜ n Γ andtransverse shear force tensor q Γ are identified as m Γ = (cid:90) t / − t / ζ P · σ Γ · P d ζ = t σ P Γ ( ε P Γ , Bend ) , (3.9) ˜ n Γ = (cid:90) t / − t / P · σ Γ · P d ζ = t σ P Γ ( ε P Γ , Mem ) , (3.10) q Γ = (cid:90) t / − t / Q · σ Γ + σ Γ · Q d ζ = t σ S Γ ( ε S Γ ) . (3.11)The moment and effective normal force tensors are symmetric, in-plane tensors and theirtwo non-zero eigenvalues are the principal moments or forces, respectively. Note that inthe case of curved shells, the physical normal force tensor is n real Γ = ˜ n Γ + H · m Γ and is, ingeneral, not symmetric, but features one zero eigenvalue just as ˜ n Γ .Based on the stress resultants, the force and moment equilibrium for an implicitly defined4 Governing equations
Reissner–Mindlin shell in strong form becomesdiv Γ n real Γ + Q · div Γ q Γ + H · ( q Γ · n Γ ) = − f , (3.12) P · div Γ m Γ − q Γ · n Γ = − c , (3.13)where f is the load vector per area and c is a distributed moment vector on the zero-isosurface Γ . Note that c needs to be in the tangent space of the zero-isosurface, i.e., c · n Γ = 0 . The equilibriums in strong form are a set of second-order surface PDEs andwith suitable boundary conditions, the complete boundary value problem of an implicitly defined shell is defined.For each field u and w there exist two non-overlapping parts at the boundary of the shell ∂ Γ . In particular, the Dirichlet boundary ∂ Γ D,i and the Neumann boundary ∂ Γ N,i , with i = { u , w } . The corresponding boundary conditions are u = ˆ g u on ∂ Γ D , u , n real Γ · n ∂ Γ + ( n Γ · q Γ · n ∂ Γ ) · n Γ = ˆ p on ∂ Γ N , u , w = ˆ g w on ∂ Γ D , w , m Γ · n ∂ Γ = ˆ m ∂ Γ on ∂ Γ N , w , (3.14)where ˆ g u are the displacements, ˆ g w are the rotations of the normal vector, ˆ p are the forcesand ˆ m ∂ Γ are the bending moments at their corresponding part of the boundary. Therotation of the normal vector at the boundary may be written in terms of ( t ∂ Γ , n ∂ Γ ) w = ( w · n ∂ Γ ) (cid:124) (cid:123)(cid:122) (cid:125) ω t ∂ Γ t ∂ Γ + ( w · t ∂ Γ ) (cid:124) (cid:123)(cid:122) (cid:125) ω n ∂ Γ n ∂ Γ . (3.15)In Fig. 5, the possible boundary conditions expressed in terms of the local triad ( t ∂ Γ , n ∂ Γ , n Γ ) and their conjugated forces ( p t ∂ Γ , p n ∂ Γ , p n Γ ) and bending moments ( m t ∂ Γ , m n ∂ Γ ) are vi-sualized. .2 Equilibrium in weak form u t ∂ Γ u n ∂ Γ u n Γ u n Γ ω n ∂ Γ ω t ∂ Γ Γ ∂ Γ yzx (a) displacements and rotations p t ∂ Γ p n ∂ Γ p n Γ p n Γ m n ∂ Γ m t ∂ Γ Γ ∂ Γ yzx (b) forces and moments Fig. 5:
Decomposition of the midsurface displacement u , difference vector w , forces and bendingmoments along the boundary ∂ Γ in terms of t ∂ Γ , n ∂ Γ and n Γ : (a) displacements androtations at the boundary, (b) forces and bending moments at the boundary. In order to convert the equilibrium in strong form to the weak form, we introduce thefollowing function spaces S u = (cid:110) v ∈ (cid:2) H (Γ) (cid:3) : v = ˆ g u on ∂ Γ D , u (cid:111) , (3.16) V u = (cid:110) v ∈ (cid:2) H (Γ) (cid:3) : v = on ∂ Γ D , u (cid:111) , (3.17) S w = (cid:110) v ∈ (cid:2) H (Γ) (cid:3) : v · n Γ = 0 ; v = ˆ g w on ∂ Γ D , w (cid:111) , (3.18) V w = (cid:110) v ∈ (cid:2) H (Γ) (cid:3) : v · n Γ = 0 ; v = 0 on ∂ Γ D , w (cid:111) , (3.19)where H is the space of functions with square integrable first derivatives. Note that thefunctions are 3D functions on the midsurface, i.e, v ( x ) : Γ → R , x ∈ Γ ⊂ R . Lateron for the discrete problem in the frame of the Trace FEM, these functions are definedin the physical space R and then restricted to Γ , with the additional condition that thefunctions need to be in [ H (Γ)] . It is emphasized that the definition of the functions inthe higher-dimensional space and then restricting them to the trace (midsurface) is one ofthe major differences compared to classical Surface FEM.Multiplying the force equilibrium, see Eq. 3.12, with a suitable test function v u and thedivergence theorem for 2D manifolds, see Eq. 2.10, the weak form of the force equilibriumcan be obtained. The task is to find u ∈ S u and w ∈ S w such that for all v u ∈ V u , there6 Trace Finite Element Method holds (cid:90) Γ ∇ dir Γ v u : ˜ n Γ + ( H · ∇ dir Γ v u ) : m Γ + ( Q · ∇ dir Γ v u ) : q Γ d A = (cid:90) Γ v u · f d A + (cid:90) ∂ Γ N , u v u · ˆ p d s . (3.20)The weak form of the moment equilibrium is obtained in a similar manner. The task is tofind u ∈ S u and w ∈ V w such that for all v w ∈ V w , there holds (cid:90) Γ ∇ dir Γ v w : m Γ + v w · ( q Γ · n Γ ) d A = (cid:90) Γ v w · c d A + (cid:90) ∂ Γ N , w v w · ˆ m ∂ Γ d s . (3.21)For further details regarding the derivation of the shell equations in the frame of the TDCand the used identities in order to obtain the equilibrium in weak form, we refer theinterested reader to [51]. The obtained weak form of the equilibrium is discretized with a trace finite element ap-proach (Trace FEM), see e.g., [45, 43, 48, 44, 31, 29, 53]. This approach is a fictitiousdomain method for surface PDEs, where the implicitly defined domain of interest Γ iscompletely immersed in a background domain Ω B ∈ R and a parametrization of Γ is nei-ther available nor needed. Let us first introduce fundamental terms and quantities, whichare required in order to discretize a surface PDE with the Trace FEM. The main ingre-dients for the Trace FEM are the definition of (1) the background mesh, (2) the discretedomain of interest (Γ h , ∂ Γ h ) , and (3) the Trace FEM function space.Firstly, a background mesh Ω B , which completely immerses the domain of interest, i.e., themidsurface of the shell Γ ∈ Ω B , needs to be defined. Without loss of generality, the mesh canbe defined by a set τ B of 3D elements and is not restricted to a certain element type. Herein,the background mesh consists of tetrahedral elements T of complete order k ≥ . Thebackground mesh is then defined by Ω B := ∪ T ∈ τ h T ∈ C . Associated to the reference element ¯ T , there is a fixed set of basis functions { N ki ( r ) } , with i = 1 , . . . , n nodes being the numberof nodes per element. The shape functions are Lagrange basis functions, i.e., N ki ( r j ) = δ ij and N ki ( r ) ∈ P k ( ¯ T ) , where P k ( ¯ T ) is the polynomial basis for 3D tetrahedral Lagrangeelements of complete order k . With an atlas of local, element-wise mappings from the 3D .1 Implementational aspects of the Trace FEM X ( r ) : ¯ T → T , with X ( r ) = N ki ( r ) x i , r ∈ ¯ T , x i ∈ T , where x i are the nodal coordinates, one may define shape functions N ki ( X ) be means of the isoparametric concept. The union of all elements forms a global set of C -continuous basis functions { M kl ( X ) } in Ω B with l = 1 , . . . , n and n being the totalnumber of nodes of the mesh. A general finite element space is then defined as Q k Ω B ,h := (cid:40) v h ∈ C (Ω B ) | v h = n nodes (cid:88) i =1 N ki ( X )ˆ v i | ˆ v i ∈ R (cid:41) ⊂ H (Ω B ) . (4.1)Secondly, the implicitly defined geometry of the shell is defined with a set of level-setfunctions, see Section 2. The continuous level-set functions ( φ, ψ i ) are interpolated withthe shape functions { M kl ( X ) } of the background mesh based on their nodal values, i.e., ˆ φ i = φ ( x i ) , ˆ ψ j,i = ψ j ( x i ) . This means that the level-set data is only needed at the nodesof the background mesh Ω B . The discrete shell midsurface Γ h and normal vector n h Γ isthen implied by φ h and the discrete boundaries of the shell ∂ Γ h may be defined either withthe discrete slave level-set functions ψ hi or the boundary of the background mesh. In thefollowing, the discrete boundary is only defined with additional slave level-set functions,otherwise the overall approach would be limited to a boundary conforming backgroundmesh.The set of elements with a non-empty intersection with Γ h is denoted by τ ΓΩ ,h and definesthe active mesh Ω Γ h := ∪ T ∈ τ ΓΩ ,h T . The definition of the active mesh is a crucial task, becausethe nodes of the active mesh imply the degrees of freedom in the numerical simulation.Lastly, the Trace FEM function space T h is established by the restriction of a “higher-dimensional finite element space” the active mesh Ω Γ h . The finite element space of theactive mesh Q k Ω Γ h is defined in a similar manner as Q k Ω B ,h , but contains only cut backgroundelements. A general Trace FEM function space T h is then defined by T h = (cid:110) v ∈ Q k Ω Γ h : v | Γ h ∈ H (Γ h ) (cid:111) ⊂ H (Ω B ) . (4.2) As mentioned above, the Trace FEM is a fictitious domain method and, compared tostandard Surface finite element approaches, three well-known challenges arise which arefurther outlined below: (i) integration of the weak form on the discrete zero-isosurface8
Trace Finite Element Method and its boundaries, (ii) stabilization of the stiffness matrix due to the restriction of theshape functions to the trace and small supports due to unfavourable cut scenarios and (iii)the enforcement of essential boundary conditions. In the following, these challenges areaddressed in detail enabling a higher-order Trace FEM approach.
The numerical integration of the domain of interest, i.e., on the zero-isosurface of φ , is anon-trivial task, in particular with higher-order accuracy. One approach which is suitablefor higher-order is presented in [39]. The numerical integration is based on a higher-orderaccurate lift of a linear reconstruction of the zero-isosurface. However, the extension tomultiple level-set functions, where the zero-isosurface of the master level-set function isrestricted with additional slave-level-set functions has not been addressed so far.Herein, we employ the integration strategy outlined by the authors in [24, 25, 26, 27]. Theadvantages of this particular approach are an optimal higher-order accurate integrationand a natural extension to multiple level-set functions. In this approach, the placement ofthe integration points on the discrete zero-isosurface is based on a higher-order, recursivereconstruction in the cut reference element. That is, the shell midsurface is reconstructedby higher-order surface elements. It is important that the reconstructed surface element is only used for the generation of integration points in the 3D reference element and may beinterpreted as an integration cell. In Fig. 6, an overview of the procedure is illustrated. Theyellow surface is the implicitly defined zero-isosurface of φ h and is restricted by additionalslave level-set functions ψ hi (green and purple surfaces). In Fig. 6(b), the cut scenarioand integration points in the reference space of the red marked element from Fig. 6(a) arevisualized. In Fig. 6(c), the integration points in the physical domain (red) and on theboundaries (blue) are shown. For further information and details, we refer to [24, 25]. Due to the restriction of the shape functions to the trace, see Eq. 4.2, the shape functionson the manifold only form a frame , which is, in general, not a basis [48, 44]. In Fig. 7,the consequences of the restriction to the trace of the manifold is visualized for a simpleexample.Let us consider a 1D manifold (blue line) embedded in three bi-linear, quadrilateral 2D .1 Implementational aspects of the Trace FEM (a) active background elements (b) integration points in refer-ence element (c) integration points in physicalelements Fig. 6: (a) Active, (black) elements in a background mesh are intersected by the shell midsurface,(b) integration points in one 3D reference element based on a (recursive) decompositionw.r.t. the master level-set function φ h and further restriction to the slave level-set func-tions ψ hi , (c) integration points in the physical domain (red points) and integration pointson the shell boundary (blue points). elements, see Fig. 7(a). Furthermore, a constant function on the manifold (black line) shallbe interpolated based on the nodal values of the active background mesh. As shown inFig. 7(b), the choice of the nodal values in order to interpolate the function on the manifoldis not unique. In particular, three different configurations, which share the same values onthe manifold are visualized. Therefore, a suitable stabilization term need to be added tothe discrete weak form, otherwise the obtained linear system of equations does not havea unique solution w.r.t. the nodal values. In addition to the restriction, depending onunfavourable cut scenarios of cut background elements, unbounded small contributions tothe stiffness matrix may occur, which causes an ill-conditioned system of linear equations.The used stabilization technique addresses both issues and is introduced for scalar-valuedproblems in [30, 11]. This stabilization technique is called “normal derivative volumestabilization” and in [31, 42], the stabilization is applied to vector-valued problems. Thestabilization term added to each unknown field in the discrete weak form is s h ( u h , v h ) := ρ (cid:90) Ω Γ h (cid:16) ∇ u h · n e,h Γ (cid:17) · (cid:16) ∇ v h · n e,h Γ (cid:17) d V , (4.3)where n e,h Γ ( x ) = ∇ φ h ( x ) (cid:107)∇ φ h ( x ) (cid:107) , x ∈ Ω Γ h is a sufficiently smooth extension of the normal vector0 Trace Finite Element Method (a) overview (b) interpolation (c) stabilization
Fig. 7:
Background meshes do not uniquely define a (constant) function on the zero-level set: (a)The black line is a constant, scalar-valued function on the manifold (blue line), which isembedded in three, bi-linear quadrilateral background elements, (b) different possibilitiesfor the interpolation of the function on the manifold (trace), (c) the stabilization addsa constraint in normal direction of the manifold resulting in a unique interpolation (redsurface). n h Γ at the zero-isosurface of φ h . It is noteworthy, that the integral is performed over thewhole active background mesh and not restricted to the trace. However, the integrandis sufficiently smooth and, therefore, a standard numerical integration scheme w.r.t. the active elements is applicable, i.e., a standard 3D Gauss rule. By adding this constraint tothe linear system of equations, the resulting system of equations features a unique solution,see Fig. 7(c). It is recommended in [30] that the stabilization parameter can be chosenwithin the following range h (cid:46) ρ (cid:46) h − , (4.4)where h is the element size of the elements from the active mesh. This stabilization tech-nique is suitable for higher-order shape functions, does not change the sparsity pattern ofthe stiffness matrix, and only first-order derivatives are needed. In addition, the implemen-tation is straightforward and the choice of the stabilization parameter is rather flexible.Other stabilization techniques are presented in [44, 11]. .2 Discretization of the Reissner–Mindlin shell As outlined in [27], the enforcement of essential of boundary conditions is a challengingtask in FDMs due to the fact that it is not possible to directly prescribe nodal values ofthe active background elements. The situation in the case of shells may be quite delicatedue to complex boundary conditions, e.g., membrane support, symmetry support, clampededges, etc., and, therefore, the treatment of boundary conditions requires special attention.The essential (Dirichlet) boundary conditions may, in principle, be enforced in a weakmanner using penalty methods, Lagrange multiplier methods, or Nitsche’s method. Herein,the non-symmetric version of Nitsche’s method is used to enforce the Dirichlet boundaryconditions [7]. The advantage of this method in the context of FDMs is that it doesnot require additional stabilization terms and the discretization of auxiliary fields suchas Lagrange multipliers is not needed. Furthermore, Nitsche’s method is a consistentapproach to enforce essential boundary conditions, which may be an advantage if higher-order convergence rates shall be achieved. For further details about the non-symmetricversion of Nitsche’s method we refer to, e.g., [7, 49, 32, 27]. The additional terms in thediscrete weak form resulting from Nitsche’s method are presented in Section 4.2.2.
In this section, the continuous weak form of the equilibrium, see Eq. 3.20 and Eq. 3.21 isdiscretized with the Trace FEM as described above. The discrete function spaces for thetrial and test functions of the midsurface displacement field are S h u = (cid:8) u h ∈ [ T h ] (cid:9) , (4.5) V h u = (cid:8) v hu ∈ [ T h ] (cid:9) . (4.6)For the discrete difference vector w h , the situation is more complicated due to the kine-matic assumptions as the difference vector needs to be tangential. The discretization oftangent vector fields on implicitly defined manifolds is a non-trivial task and detailed inSection 4.2.1.2 Trace Finite Element Method
Different approaches for the discretization of tangent vector fields are presented in, e.g.,in [51, 42, 37]. Herein, the difference vector w h and its corresponding test function v hw are defined in the general Trace FEM function space without the tangentiality constraint,similar to [42]. The corresponding function spaces are S h w = (cid:8) w h ∈ [ T h ] (cid:9) , (4.7) V h w = (cid:8) v hw ∈ [ T h ] (cid:9) , (4.8) but in the discrete weak form, see Eq. 4.11, only the projected difference vector and testfunction are used, i.e., (cid:101) w h = P · w h , (cid:101) v hw = P · v wh . One may argue that in this approachthe derivatives of the projector P occur, which involves surface derivatives of the normalvector n h Γ . This does not lead to additional computational costs, because in the case ofthe Reissner–Mindlin shell, the Weingarten map H directly appears in the weak form and,therefore, the surface derivatives of the normal vector n h Γ are required independently ofthis approach.As a result of this projection, only the tangential part of w h and v hw is considered in thediscrete weak form and, therefore, the tangentiality constraint is built-in automatically. Itis clear, that without any further measures this would lead to an ill-conditioned systemof equations, because the normal part of w h , i.e., w hn = w h · n h Γ , does not appear in thediscrete weak form and, therefore, w hn is not unique. In order to address this issue, a simplestabilization term, similar to the penalty term in [42] is introduced s w,h := ρ w (cid:90) Γ h (cid:0) w h · n h Γ (cid:1) (cid:0) v hw · n h Γ (cid:1) d A , (4.9)where ρ w is a suitable stabilization parameter. A series of numerical studies regarding thechoice of the stabilization parameter for the Reissner–Mindlin shell has been conducted onflat and curved shell geometries. In detail, the dependency on the (1) material parameter E , (2) thickness t and (3) element size on h w.r.t. the condition number of the stiffnessmatrix and the influence on the results was investigated. Summarizing the outcome of thenumerical studies, the stabilization parameter can be chosen independently of h and for asuitable scaling of the stabilization term, the parameter is set to ρ w = E · t . (4.10) .2 Discretization of the Reissner–Mindlin shell
23A difference of the proposed approach and the method shown in [42] is that only theprojected part of the vector field, i.e., (cid:101) w h = P · w h , is used in the discrete weak formwhich directly enforces the tangentiality constraint. Furthermore, the used stabilizationparameter, which is in [42, 37] a penalty parameter, does not depend on h and only asuitable scaling of the stabilization term may be required. Based on the previous definitions, the discrete weak form with the Trace FEM of theforce equilibrium, see Eq. 3.20, reads as follows: Given material parameters ( E, ν ) ∈ R + ,body forces f ∈ R on Γ h , tractions ˆ p on ∂ Γ h N , u , stabilization parameter ρ ∈ R + find thedisplacement fields ( u h , w h ) ∈ S h u × S h w such that for all test functions ( v hu , v hw ) ∈ V h u × V h w there holds in Γ h (cid:90) Γ h ∇ dir Γ v u : ˜ n Γ ( u h ) + ( H · ∇ dir Γ v u ) : m Γ ( u h , (cid:101) w h ) + ( Q · ∇ dir Γ v u ) : q Γ ( u h , (cid:101) w h ) d A − (cid:90) ∂ Γ h D , u v u · p ( u h , (cid:101) w h ) d s (cid:124) (cid:123)(cid:122) (cid:125) boundary term due to v hu (cid:54) = ∂ Γ D , u + (cid:90) ∂ Γ h D , u u h · p ( v hu , (cid:101) v hw ) d s (cid:124) (cid:123)(cid:122) (cid:125) Nitsche term for displ. on LHS + ρ (cid:90) Ω Γ h (cid:16) ∇ u h · n e,h Γ (cid:17) · (cid:16) ∇ v hu · n e,h Γ (cid:17) d V (cid:124) (cid:123)(cid:122) (cid:125) Trace FEM stabilization, see Section 4.1.2 = (cid:90) Γ h v hu · f d A + (cid:90) ∂ Γ h D , u ˆ g u · p ( v hu , (cid:101) v hw ) d s (cid:124) (cid:123)(cid:122) (cid:125) Nitsche term for displ. on RHS + (cid:90) ∂ Γ h N , u v hu · ˆ p d s , (4.11)where p = n real Γ · n h∂ Γ + (cid:0) n h Γ · q Γ · n h∂ Γ (cid:1) n h Γ , see Eq. 3.14, are the conjugated forces at theDirichlet boundary ∂ Γ h D , u .The discrete weak form of the moment equilibrium, see Eq. 3.21 reads as follows: Givenmaterial parameters ( E, ν ) ∈ R + , distributed moments c ∈ T P Γ h on Γ h , bending mo-ments ˆ m ∂ Γ on ∂ Γ h N , w , stabilization parameters ( ρ, ρ w ) ∈ R + find the displacement fields4 Numerical results ( u h , w h ) ∈ S h u × S h w such that for all test functions ( v hu , v hw ) ∈ V h u × V h w there holds in Γ h (cid:90) Γ h ∇ dir Γ (cid:101) v hw : m Γ ( u h , (cid:101) w h ) + (cid:101) v hw · (cid:104) q Γ ( u h , (cid:101) w h ) · n Γ (cid:105) d A − (cid:90) ∂ Γ D , w (cid:101) v hw · m ∂ Γ ( u h , (cid:101) w h ) d s (cid:124) (cid:123)(cid:122) (cid:125) boundary term due to v hw (cid:54) = ∂ Γ D , w + (cid:90) ∂ Γ D , w (cid:101) w h · m ∂ Γ ( v hu , (cid:101) v hw ) d s (cid:124) (cid:123)(cid:122) (cid:125) Nitsche term for rot. on LHS + ρ (cid:90) Ω Γ h (cid:16) ∇ w h · n e,h Γ (cid:17) · (cid:16) ∇ v hw · n e,h Γ (cid:17) d V (cid:124) (cid:123)(cid:122) (cid:125) Trace FEM stabilization, see Section 4.1.2 + ρ w (cid:90) Γ h (cid:0) w h · n h Γ (cid:1) (cid:0) v hw · n h Γ (cid:1) d A (cid:124) (cid:123)(cid:122) (cid:125) stabilization term for (cid:101) w h see Section 4.2.1 = (cid:90) Γ h v hw · c d A + (cid:90) ∂ Γ h D , w ˆ g w · m ∂ Γ ( v hu , (cid:101) v hw ) d s (cid:124) (cid:123)(cid:122) (cid:125) Nitsche term for rot. on RHS + (cid:90) ∂ Γ h N , w v hw · ˆ m ∂ Γ d s , (4.12)where m ∂ Γ = m Γ · n h∂ Γ , see Eq. 3.14, are the conjugated bending moments at the Dirichletboundary ∂ Γ h D , w .The usual element assembly w.r.t. the active elements yields a linear system of equationsin the following form ( K Stiff + K Nitsche + K Stab ) (cid:124) (cid:123)(cid:122) (cid:125) K · (cid:34) ˆ u ˆ w (cid:35) = ( b Load + b Nitsche ) (cid:124) (cid:123)(cid:122) (cid:125) b , (4.13)with [ ˆ u , ˆ w ] being the sought displacements and rotations of the normal vector at the nodesof the active elements. The matrix K and the load vector b are split into: (1) standardterms for the stiffness matrix, (2) boundary terms, (3) stabilization terms, respectively. In this section, the proposed numerical method for implicitly defined Reissner–Mindlinshells is tested on a set of benchmark examples, consisting of the partly clamped hyperbolicparabolid from [2, 14], the partly clamped gyroid from [28] and a clamped flower-shapedshell inspired by [50, 51]. In the case of the first two examples the shells are rather thin andlocking phenomena can be expected, especially in the case of low ansatz orders. However,when increasing the order p , locking phenomena decrease significantly and, therefore, no .1 Hyperbolic paraboloid u h , w h , and the interpolation of the level-setfunctions, the same order of shape functions are used. The orders are varied as ≤ p ≤ .The element size h is proportional to the factor n which is related to the number of elementsand is varied between ≤ n ≤ .In the presented examples, the stabilization parameter ρ for the normal derivative volumestabilization is set to ρ = / h . In order to achieve a proper scaling of the stiffness matrix,the parameter ρ w , is set to ρ w = Et , as proposed in Section 4.2.1. The first example is the partly clamped hyperbolic paraboloid and is taken from [2, 14]. Theproblem is defined in Fig. 8. The yellow surface is the zero-isosurface of the master level-set function φ and the grey planes are the zero-isosurfaces of the slave level-set functions ψ j , j ∈ [1 , , which define the boundaries of the shell. The blue line is the clamped edgeof the shell. In Fig. 9(a), the active background mesh which contains only cut elements Geometry: Hyperbolic paraboloid φ ( x ) = x − y − zx ∈ [ − . , . y ∈ [ − . , . t = 0 . Material parameters: E = 2 . × ν = 0 . α s = 1 . Load: Gravity load f = [0 , , − · t ] T c = Support: Clamped edge at x = − / Ref. displacement: | u z,i, Ref | = 9 . × − x i = (0 . , , . T Fig. 8:
Definition of the partly clamped hyperbolic paraboloid problem. with p = 4 is shown. In Fig. 9(b), the corresponding integration points are illustrated,those in the domain are plotted in red and those on the boundaries are blue. In Fig. 9(c),the numerical solution of the partly clamped hyperbolic paraboloid is presented. The grey6 Numerical results surface is the undeformed zero-isosurface and the colors on the deformed midsurface of theshell are the Euclidean norm of the displacement field u h . The displacements are magnifiedby a factor of × . (a) active background mesh (b) integration points (c) displacements Fig. 9: (a) Active background mesh, which consists only of cut elements, (b) automatically gen-erated integration points in the domain (red) and on the boundaries (blue), (c) deformedzero-isosurface with scaled displacements u by a factor of × . In the convergence studies, the vertical displacement at x i is compared with the givenreference displacement. Due to the moderate complexity of the master level-set function,the numerical solution converges rather fast and the element factor n is only varied between ≤ n ≤ . In Fig. 10(a), the results of the convergence study are presented. In particular,the normalized displacement u z,i / u z,i, Ref is plotted as a function of the element size h ∼ / n . The behaviour of the convergence is in agreement with the results shown, e.g., in[2, 38, 51]. In particular, the expected locking behaviour is more pronounced for p = 2 anddecreases significantly for higher orders. In Fig. 10(b), the normalized, estimated conditionnumber of the stiffness matrix is plotted as a function of the element size h ∼ / n . Thecondition numbers are obtained with the MATLAB function condest . It can be seen thatthe condition numbers increase with quadratic order as expected for second-order PDEs.The jump between the element orders is well-known in the context of higher-order finiteelement approaches, see e.g., in [26]. .1 Hyperbolic paraboloid (a) convergence (b) condition number Fig. 10: (a) Normalized convergence of reference displacement u z,i, Ref = − . × − atpoint x i = (0 . , , . T , (b) normalized condition numbers, the reference value is . × , which is the condition number at n = p = 2 . Numerical results
The next test case is a partly clamped gyroid and is taken from [28]. The problem isdefined in Fig. 11. Similar as above, the yellow surface is the zero-isosurface of the masterlevel-set function φ , the grey planes are the zero-isosurfaces of the slave level-set functions ψ j , j ∈ [1 , , which bound the master leve-set function. The blue curve is the clampededge of the shell.In contrast to [28], a factor π is inserted into the arguments of the trigonometric functionsof the master level-set function, which is given in [28, Eq. 43]. Otherwise, the obtainedgeometry is not in agreement with the presented geometry in [28, Fig. 12]. In addition,the load is decreased by one order of magnitude in order to decrease the deformations,which shall be significantly smaller than the dimensions of the shell. Therefore, the givenreference displacement needs to be scaled accordingly to . . However, in [28], adifferent shell model (seven-parameter shell model) is used. Herein, the classical Reissner–Mindlin shell, which is often labelled as five-parameter model, is used and, therefore, wecan expect small differences in the displacements. For the Reissner–Mindlin shell model theconverged reference displacement is . , which is a relative error of . compared tothe seven-parameter model, which can be explained with the differences in the kinematicassumptions between the two shell models. This discrepancy could be decreased with asuitable shear correction factor. Geometry: Gyroid φ ( x ) = sin( πx ) cos( πy ) +sin( πy ) cos( πz ) +sin( πz ) cos( πx ) x ∈ [0 , y ∈ [ − . , . z ∈ [ − . , . t = 0 . Material parameters: E = 70 × ν = 0 . α s = 1 . Load: Gravity load f = [0 , , · t ] T c = Support: Clamped edge at x = 0 Ref. displacement: | u z,i, Ref | = 0 . x i = (2 , . , − . T Fig. 11:
Definition of the partly clamped gyroid problem.
Analogously to the example above, in Fig. 12, the active background mesh for p = 4 and .2 Gyroid u h andthe grey surface indicates the undeformed zero-isosurface. (a) active background mesh (b) integration points Fig. 12: (a) Active background mesh, which consists only of cut elements, (b) automaticallygenerated integration points in the domain (red) and on the boundaries (blue).
In the convergence study, the parameter n is varied between ≤ n ≤ . The coarserlevels n = { , } are skipped due to the more complex shape of the shell. In Fig. 13(b),the results of the convergence analyses are presented. Similar to the first test case, thenormalized displacement u z,i / u z,i, Ref is plotted as a function of the element size h ∼ / n . Forthe lower orders p = { , } the expected locking phenomena is more pronounced comparedto the example before. Nevertheless, it is clearly seen that the accuracy for higher ordersincreases significantly and the behaviour of convergence is in agreement with the resultsshown, e.g., in [28]. The condition numbers for this example behave in a similar manneras in Fig. 10(b) and, therefore, the plot is omitted for the sake of brevity.0 Numerical results (a) displacements (b) convergence
Fig. 13: (a) Deformed zero-isosurface, (b) normalized convergence of reference displacement u z,i, Ref = 0 . at point x i = (2 , . , − . T . .3 Flower-shaped shell The geometry of the last example is inspired from [51]. The problem is defined in Fig. 14.Analogously as above, the yellow surface is the zero-isosurface of the master level-set func-tion φ and the intersection with the slave level-set function ψ (grey surface) defines theboundaries of the shell.A characteristic feature of this test case is that smooth solutions in all involved fields can beexpected and, therefore, optimal higher-order convergences rates are enabled. In contrastto the examples before, a reference displacement is not available for the particular example.For the error measurement we use the concept of residual errors in a similar manner asshown in [23, 50, 51]. The residual errors are the evaluation of the equilibrium in strongform, see Section 3.1, integrated over the domain. In particular, the element-wise L -errors Geometry: Flower shell φ ( x ) = / (cid:16) − s (cid:17) + / (cid:16) x − y (cid:17) − zs = r − . . . θ ) r, θ are polar coordinates of x, yψ ( x ) = / (cid:16) x − y (cid:17) − zt = 0 . Material parameters: E = 1 . × ν = 0 . α s = 1 . Load: f = − [1 , , T c = Support: Clamped edgesError measurement Residual errors
Fig. 14:
Definition of shallow flower-shaped shell problem of the force and moment equilibrium are computed in the convergence analyses ε rel,residual,F = τ ΓΩ ,h (cid:88) T =1 (cid:82) T (cid:2) div Γ n real Γ + Q · div Γ q Γ + H · ( q Γ · n Γ ) + f (cid:3) d A (cid:82) T f d A , (5.1) ε residual,M = τ ΓΩ ,h (cid:88) T =1 (cid:90) T [ P · div Γ m Γ − q Γ · n Γ + c ] d A . (5.2)For the computation of the residual errors, second-order surfaces derivatives are requiredwhich implies a theoretical optimal order of convergence O ( p − . The numerical solution2 Conclusions of the problem is visualized in Fig. 15(b) for p = 4 . The displacements are scaled by afactor of × . The colors on the deformed zero-isosurface are the Euclidean norm ofthe displacement field u h . The corresponding integration points are shown in Fig. 15(a) inthe same style than before. (a) integration points (b) displacements Fig. 15: (a) Automatically generated integration points in the domain (red) and on the bound-aries (blue), (b) deformed zero-isosurface with scaled displacements u by a factor of × . In the convergence analyses the parameter n is varied between ≤ n ≤ . The results areplotted in Fig. 16. It is noticeable that the pre-asymptotic range is more pronounced forthe lower ansatz orders p = { , } . Nevertheless, it is clear that higher-order convergencerates are achieved in both residual errors of Eqs. 5.1 and 5.2. In comparison to the originaltest case for the Reissner–Mindlin shell presented in [51], the behaviour of the convergencein the residual errors obtained here with the Trace FEM is in very good agreement to theresults in [51] obtained with isogeometric analysis. A higher-order accurate Trace FEM approach for implicitly defined shells is presented.The shell geometry is implicitly defined by means of multiple level-set functions. Dueto the implicit geometry definition, a parametrization of the midsurface is not availableand the classical shell formulations of the linear Reissner–Mindlin shell are not applicable.Therefore, a more general shell formulation in the frame of the TDC, which extends alsoto implicitly defined shells, is employed.3 (a) force equilibrium (b) moment equilibrium
Fig. 16: (a) residual error of the force equilibrium ε rel,residual,F , (b) residual error of the momentequilibrium ε residual,M . The Trace FEM is a fictitious domain method for PDEs on manifolds enabling higher-orderaccuracy when the following three aspects (well-known for general FDMs) are addressedproperly: (i) numerical integration, (ii) stabilization and (iii) enforcement of essentialboundary conditions. Each aspect is carefully detailed herein, thus proposing an optimalhigher-order accurate FDM for shells for the first time. The employed integration techniqueextends to multiple level-set functions and is based on a recursive reconstruction of the cutreference element into higher-order integration cells. The normal derivative volume stabi-lization also extends to higher-order shape functions without additional measures and thechoice of the stabilization parameter is rather flexible. The essential boundary conditionsare enforced with the non-symmetric version of Nitsche’s method which is a consistentapproach and does neither require additional stabilization terms nor the discretization ofauxiliary fields. In addition, the tangentiality constraint on the difference vector is au-tomatically built-in with an additional projection of a full 3D vector combined with aconsistent stabilization term. In the numerical results, classical and new benchmark ex-amples are presented and optimal higher-order convergence rates are achieved when thephysical fields are sufficiently smooth.4
REFERENCES
References [1] Başar, Y.; Krätzig, W.B.:
Mechanik der Flächentragwerke . Vieweg + Teubner Verlag,Braunschweig, 1985.[2] Bathe, K.J.; Iosilevich, A.; Chapelle, D.: An evaluation of the MITC shell elements.
Computers & Structures , , 1–30, 2000.[3] Bischoff, M.; Ramm, E.; Irslinger, J.: Models and Finite Elements for Thin-WalledStructures. In Encyclopedia of Computational Mechanics (Second Edition) . (Stein, E.;Borst, R.; Hughes, T. J.; Hughes, T. J., Eds.), John Wiley & Sons, 2017.[4] Blaauwendraad, J.; Hoefakker, J.H.:
Structural Shell Analysis , Vol. 200,
Solid Me-chanics and Its Applications . Springer, Berlin, 2014.[5] Bonito, A.; Demlow, A.; Nochetto, R.H.: Chapter 1 - Finite Element Methods forthe Laplace-Beltrami Operator. In
Geometric Partial Differential Equations - Part I .(Bonito, A.; Nochetto, R.H., Eds.), Vol. 21,
Handbook of Numerical Analysis , Elsevier,Amsterdam, 1–103, 2020.[6] Brandner, P.; Reusken, A.: Finite element error analysis of surface Stokes equationsin stream function formulation. arXiv e-prints , 2019. arXiv: 1910.09221.[7] Burman, E.: A penalty-free nonsymmetric Nitsche-type method for the weak imposi-tion of boundary conditions.
SIAM J. Numer. Anal. , , 1959–1981, 2012.[8] Burman, E.; Claus, S.; Hansbo, P.; Larson, M.G.; Massing, A.: CutFEM: Discretizinggeometry and partial differential equations. Internat. J. Numer. Methods Engrg. , ,472–501, 2015.[9] Burman, E.; Elfverson, D.; Hansbo, P.; Larson, M.G.; Larsson, K.: Shape optimiza-tion using the cut finite element method. Comp. Methods Appl. Mech. Engrg. , ,242–261, 2018.[10] Burman, E.; Hansbo, P.; Larson, M.G.: A stabilized cut finite element method forpartial differential equations on surfaces: The Laplace-Beltrami operator. Comp.Methods Appl. Mech. Engrg. , 188–207, 2015.
EFERENCES
ESAIM: Mathematical Modelling and Numerical Analysis , , 2247–2282, 2018.[12] Calladine, C. R.: Theory of Shell Structures . Cambridge University Press, Cambridge,1983.[13] Cenanovic, M.; Hansbo, P.; Larson, M.G.: Cut finite element modeling of linearmembranes.
Comp. Methods Appl. Mech. Engrg. , , 98–111, 2016.[14] Chapelle, D.; Bathe, K.J.: Fundamental considerations for the finite element analysisof shell structures. Computers & Structures , , 19–36, 1998.[15] Chapelle, D.; Bathe, K.J.: The mathematical shell model underlying general shellelements. Internat. J. Numer. Methods Engrg. , , 289–313, 2000.[16] Delfour, M.C.; Zolésio, J.P.: A Boundary Differential Equation for Thin Shells. J.Differential Equations , , 426–449, 1995.[17] Delfour, M.C.; Zolésio, J.P.: Tangential Differential Equations for Dynamical ThinShallow Shells. J. Differential Equations , , 125–167, 1996.[18] Delfour, M.C.; Zolésio, J.P.: Shapes and Geometries: Metrics, Analysis, DifferentialCalculus, and Optimization . SIAM, Philadelphia, 2011.[19] Demlow, A.: Higher-order finite element methods and pointwise error estimates forelliptic problems on surfaces.
SIAM J. Numer. Anal. , , 805–827, 2009.[20] Dziuk, G.: Finite Elements for the Beltrami operator on arbitrary surfaces , 142–155.Springer, Berlin, 1988.[21] Dziuk, G.; Elliott, C.M.: Finite element methods for surface PDEs.
Acta Numerica , , 289–396, 2013.[22] Farshad, M.: Design and Analysis of Shell Structures . Springer, Berlin, 1992.[23] Fries, T.P.: Higher-order surface FEM for incompressible Navier-Stokes flows on man-ifolds.
Int. J. Numer. Methods Fluids , , 55–78, 2018.[24] Fries, T.P.; Omerović, S.: Higher-order accurate integration of implicit geometries. Internat. J. Numer. Methods Engrg. , , 323–371, 2016.6 REFERENCES [25] Fries, T.P.; Omerović, S.; Schöllhammer, D.; Steidl, J.: Higher-order meshing ofimplicit geometries - Part I: Integration and interpolation in cut elements.
Comp.Methods Appl. Mech. Engrg. , , 759–784, 2017.[26] Fries, T.P.; Schöllhammer, D.: Higher-order meshing of implicit geometries - Part II:Approximations on manifolds. Comp. Methods Appl. Mech. Engrg. , , 270–297,2017.[27] Fries, T.P.; Schöllhammer, D.: A unified finite strain theory for membranes and ropes. Comp. Methods Appl. Mech. Engrg. , , 113031, 2020.[28] Gfrerer, M.H.; Schanz, M.: High order exact geometry finite elements for seven-parameter shells with parametric and implicit reference surfaces. Comput. Mech. , ,133–145, 2019.[29] Grande, J.; Lehrenfeld, C.; Reusken, A.: Analysis of a high-order trace finite elementmethod for PDEs on level set surfaces. SIAM J. Numer. Anal. , , 228–255, 2018.[30] Grande, J.; Reusken, A.: A higher order finite element method for partial differentialequations on surfaces. SIAM , , 388–414, 2016.[31] Gross, S.; Jankuhn, T.; Olshanskii, M.A.; Reusken, A.: A trace finite element methodfor vector-laplacians on surfaces. SIAM J. Numer. Anal. , , 2406–2429, 2018.[32] Guo, Y.; Do, H.; Ruess, M.: Isogeometric stability analysis of thin shells: From simplegeometries to engineering models. Internat. J. Numer. Methods Engrg. , , 433–458,2019.[33] Hansbo, P.; Larson, M.G.: Finite element modeling of a linear membrane shell problemusing tangential differential calculus. Comp. Methods Appl. Mech. Engrg. , , 1–14,2014.[34] Hansbo, P.; Larson, M.G.; Larsson, F.: Tangential differential calculus and the finiteelement modeling of a large deformation elastic membrane problem. Comput. Mech. , , 87–95, 2015.[35] Jankuhn, T.; ; Reusken, A.: Higher Order Trace Finite Element Methods for theSurface Stokes Equations. arXiv e-prints , 2019. arXiv: 1909.08327. EFERENCES
Interfaces Free Bound. , , 353–377, 2018.[37] Jankuhn, T.; Olshanskii, M.A.; Reusken, A.; Zhiliako, A.: Error Analysis of HigherOrder Trace Finite Element Methods for the Surface Stokes Equations. arXiv e-prints ,2020. arXiv: 2003.06972.[38] Kiendl, J.; Marino, E.; De Lorenzis, L.: Isogeometric collocation for the Reissner-Mindlin shell problem. Comp. Methods Appl. Mech. Engrg. , , 645–665, 2017.[39] Lehrenfeld, C.: High order unfitted finite element methods on level set domains usingisoparametric mappings. Comp. Methods Appl. Mech. Engrg. , , 716–733, 2016.[40] Lehrenfeld, C.; Olshanskii, M.A.; Xu, X.: A stabilized trace finite element method forpartial differential equations on evolving surfaces. SIAM , , 1643–1672, 2018.[41] Müller, B.; Kummer, F.; Oberlack, M.: Highly accurate surface and volume integra-tion on implicit domains by means of moment-fitting. Internat. J. Numer. MethodsEngrg. , , 512–528, 2013.[42] Olshanskii, M.A.; Quaini, A.; Reusken, A.; Yushutin, V.: A finite element method forthe surface Stokes problem. SIAM J. Sci. Comput. , , A2492–A2518, 2018.[43] Olshanskii, M.A.; Reusken, A.: A finite element method for surface PDEs: Matrixproperties. Numer. Math. , , 491–520, 2009.[44] Olshanskii, M.A.; Reusken, A.: Trace finite element methods for PDEs on surfaces. Lecture Notes in Computational Science and Engineering , , 211–258, 2017.[45] Olshanskii, M.A.; Reusken, A.; Grande, J.: A finite element method for elliptic equa-tions on surfaces. SIAM , , 3339–3358, 2009.[46] Olshanskii, M.A.; Xu, X.: A trace finite element method for PDEs on evolving sur-faces. SIAM , , A1301–A1319, 2017.[47] Osher, S.; Fedkiw, R.P.: Level Set Methods and Dynamic Implicit Surfaces . Springer,Berlin, 2003.8
REFERENCES [48] Reusken, A.: Analysis of trace finite element methods for surface partial differentialequations.
IMA J. Numer. Anal. , , 1568–1590, 2014.[49] Schillinger, D.; Harari, I.; Hsu, M.C.; Kamensky, D.; Stoter, S.K.F.; Yu, Y.; Zhao,Y.: The non-symmetric Nitsche method for the parameter-free imposition of weakboundary and coupling conditions in immersed finite elements. Comp. Methods Appl.Mech. Engrg. , , 625–652, 2016.[50] Schöllhammer, D.; Fries, T.P.: Kirchhoff-Love shell theory based on tangential differ-ential calculus. Comput. Mech. , , 113–131, 2019.[51] Schöllhammer, D.; Fries, T.P.: Reissner–Mindlin shell theory based on tangentialdifferential calculus. Comp. Methods Appl. Mech. Engrg. , , 172–188, 2019.[52] Schöllhammer, D.; Fries, T.P.: Reissner–Mindlin shell theory based on TangentialDifferential Calculus. John Wiley & Sons, Chichester, 2019, PAMM.[53] Schöllhammer, D.; Fries, T.P.: A unified approach for shell analysis on explicitly andimplicitly defined surfaces. Form and Force: Proceedings of the IASS Symposium 2019- Structural Membranes 2019, C. Lázaro, K.U. Bletzinger, E. Oñate (eds.) , Barcelona,750–757, 2019, Form and Force: Proceedings of the IASS Symposium 2019 - StructuralMembranes 2019, C. Lázaro, K.U. Bletzinger, E. Oñate (eds.).[54] Sethian, J.A.:
Level Set Methods and Fast Marching Methods . Cambridge UniversityPress, Cambridge, 2nd edition, 1999.[55] Simo, J.C.; Fox, D.D.: On a stress resultant geometrically exact shell model. Part I:Formulation and optimal parametrization.
Comp. Methods Appl. Mech. Engrg. , ,267–304, 1989.[56] Simo, J.C.; Fox, D.D.; Rifai, M.S.: On a stress resultant geometrically exact shellmodel. Part II: The linear theory; Computational aspects. Comp. Methods Appl.Mech. Engrg. , , 53–92, 1989.[57] Wempner, G.; Talaslidis, D.: Mechanics of Solids and Shells: Theories and Approxi-mations . CRC Press LLC, Florida, 2002.[58] Zingoni, A.: