Cut finite element modeling of linear membranes
CCut finite element modeling of linear membranes
Mirza Cenanovicy ∗ , Peter Hansbo † , and Mats G. Larson ‡ Department of Mechanical Engineering, Jönköping University,School of Engineering, SE-55111 Jönköping, Sweden Department of Mathematics and Mathematical Statistics, UmeåUniversity, SE-90187 Umeå, SwedenNovember 7, 2018
Abstract
We construct a cut finite element method for the membrane elasticity problem onan embedded mesh using tangential differential calculus. Both free membranes andmembranes coupled to 3D elasticity are considered. The discretization comes from aGalerkin method using the restriction of 3D basis funtions (linear or trilinear) to thesurface representing the membrane. In the case of coupling to 3D elasticity, we viewthe membrane as giving additional stiffness contributions to the standard stiffnessmatrix resulting from the discretization of the three–dimensional continuum.
In this paper we construct finite element methods for linearly elastic membranes embed-ded in three dimensional space meshed by tetrahedral or hexahedral elements. Thesemeshes do not in general align with the surface of the membrane which instead cutsthrough the elements. For the modeling of the membrane problems we use the tan-gential differential calculus employed by Hansbo and Larson [5] for meshed membranes.We extend this approach following Olshanskii, Reusken, and Grande [7] and constructa Galerkin method by using restrictions of the 3D basis functions defined on the three–dimensional mesh to the surface. This approach can lead to severe ill conditioning, so ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ m a t h . NA ] N ov e adapt a stabilization technique proposed by Burman, Hansbo, and Larson [3] for theLaplace–Beltrami operator to the membrane problem.The main application that we have in mind is the coupling of membranes to 3D elas-ticity. This allows for the modelling of reinforcements, such as shear strengthening andadhesive layers. In the case of adhesives, the method can be further refined using theimperfect bonding approach of Hansbo and Hansbo [4], so that cut 3D elements areused, allowing for relative motion of the continuum on either side of the adhesive. Thisextension is not explored in this paper but has been considered, for adhesives, in a discon-tinuous Galerkin setting in [6]. Here we restrict the method to adding membrane stiffnessto a continuous 3D approximation. The idea of adding stiffness from lower–dimensionalstructures is a classical approach, cf. Zienkiewicz [8, Chapter 7.9], using element sides oredges as lower dimensional entities. Letting the membranes cut through the elements inan arbitrary fashion considerably increases the practical modeling possibilities.The paper is organized as follows: in Section 2 we introduce the membrane modelproblem and the finite element method for membranes and embedded membranes; inSection 3 we describe the implementation details of the method; and in Section 4 wepresent numerical results. In what follows, Γ denotes a an oriented surface, which is embedded in R and equippedwith exterior normal n Γ . The boundary of Γ consists of two parts, ∂ Γ N , where zerotraction boundary conditions are assumed, and ∂ Γ D where zeros Dirichelt boundaryconditions are assumed.We let ρ denote the signed distance function fulfilling ∇ ρ | Γ = n Γ .For a given function u : Γ → R we assume that there exists an extension ¯ u , in someneighborhood of Γ , such that ¯ u | Γ = u . The the tangent gradient ∇ Γ on Γ can be definedby ∇ Γ u = P Γ ∇ u (1)with ∇ the R gradient and P Γ = P Γ ( x ) the orthogonal projection of R onto thetangent plane of Γ at x ∈ Γ given by P Γ = I − n Γ ⊗ n Γ (2)where I is the identity matrix. The tangent gradient defined by (1) is easily shown tobe independent of the extension u . In the following, we shall consequently not makethe distinction between functions on Γ and their extensions when defining differentialoperators. 2he surface gradient has three components, which we shall denote by ∇ Σ u =: (cid:18) ∂u∂x Γ , ∂u∂y Γ , ∂u∂z Γ (cid:19) . For a vector valued function v ( x ) , we define the tangential Jacobian matrix as the trans-pose of the outer product of ∇ Γ and v , ( ∇ Γ ⊗ v ) T := ∂v ∂x Γ ∂v ∂y Γ ∂v ∂z Γ ∂v ∂x Γ ∂v ∂y Γ ∂v ∂z Γ ∂v ∂x Γ ∂v ∂y Γ ∂v ∂z Γ , the surface divergence ∇ Γ · v := tr ∇ Γ ⊗ v , and the in–plane strain tensor ε Γ ( u ) := P Γ ε ( u ) P Γ , where ε ( u ) := 12 (cid:0) ∇ ⊗ u + ( ∇ ⊗ u ) T (cid:1) is the 3D strain tensor. We consider, following [5], the problem of finding u : Γ → R such that −∇ Γ · σ Γ ( u ) = f on Γ , σ Γ = 2 µ ε Γ + λ tr ε Γ P Γ on Γ . (3)where f : Γ → R is a load per unit area and, with Young’s modulus E and Poisson’sratio ν , µ := E ν ) , λ := Eν − ν are the Lamé parameters in plane stress.Splitting the displacement into a normal part u N := u · n Γ and a tangential part u T := u − u N n Γ , the corresponding weak statement takes the form: find u ∈ V := { v : v N ∈ L (Γ) and v T ∈ [ H (Γ)] , v = on Γ D } , such that a ( u , v ) = l ( v ) , ∀ v ∈ V, (4)3here a ( u , v ) = (2 µ ε Γ ( u ) , ε Γ ( v )) Γ + ( λ ∇ Γ · u , ∇ Γ · v ) Γ , l ( v ) = ( f , v ) Γ , and ( v , w ) Γ = (cid:90) Γ v · w d Γ and ( ε Γ ( v ) , ε Γ ( w )) Γ = (cid:90) Γ ε Γ ( v ) : ε Γ ( w ) d Γ are the L inner products. Let (cid:101) T h be a quasi uniform mesh, with mesh parameter < h ≤ h , into shape regulartetrahedra (hexahedra will be briefly discussed in Section 3) of an open and boundeddomain Ω in R completely containing Γ . On (cid:101) T h , let φ be a continuous, piecewise linearapproximation of the distance function ρ and define the discrete surface Γ h as the zerolevel set of φ ; that is Γ h = { x ∈ Ω : φ ( x ) = 0 } (5)We note that Γ h is a polygon with flat faces and we let n h be the piecewise constantexterior unit normal to Γ h . For the mesh (cid:101) T h , we define the active background mesh by T h = { T ∈ (cid:101) T h : T ∩ Γ h (cid:54) = ∅} (6)cf. Figure 2, and its set of interior faces by F h = { F = T + ∩ T − : T + , T − ∈ T h } (7)The face normals n + F and n − F are then given by the unit normal vectors which areperpendicular on F and are pointing exterior to T + and T − , respectively. We observethat the active background mesh T h gives raise to a discrete h –tubular neighborhood of Γ h , which we denote by Ω h = ∪ T ∈T h T (8)with boundary ∂ Ω h consisting of element faces F constituting the trace mesh on T h .Note that for all elements T ∈ T h there is a neighbor T (cid:48) ∈ T h such that T and T (cid:48) sharea face. We denote by ∂ Ω h, D the boundary of Ω h intersected by Γ h, D Finally, let (cid:101) V h denote the space of continuous piecewise linear (or trilinear) polynomialsdefined on (cid:101) T h and V h = (cid:110) v ∈ [ (cid:101) V h | Ω h ] : v = on ∂ Ω h, D (cid:111) (9)be the space of continuous piecewise linear polynomials defined on T h . The finite element4ethod on Γ h takes the form: find u h ∈ V h such that A h ( u h , v ) = l h ( v ) ∀ v ∈ V h (10)Here the bilinear form A h ( · , · ) is defined by A h ( v , w ) = a h ( v , w ) + j h ( v , w ) ∀ v , w ∈ V h (11)with a h ( v , w ) = (2 µ ε Γ h ( v ) , ε Γ h ( w )) Γ h + ( λ ∇ Γ h · v , ∇ Γ h · w ) Γ h (12)and j h ( v, w ) = (cid:88) F ∈F h (cid:90) F τ (cid:2) n + F · ∇ v (cid:3) · (cid:2) n + F · ∇ w (cid:3) ds. (13)Here [ v ] = v + − v − , where w ( x ) ± = lim t → + w ( x ∓ t n + F ) , denotes the jump of v acrossthe face F , and τ is a constant of O (1) . The tangent gradients are defined using thenormal to the discrete surface ∇ Γ h v = P Γ h ∇ v = ( I − n h ⊗ n h ) ∇ v (14)and the right hand side is given by l h ( v ) = ( f , v ) Γ h (15) We next consider the case of a membrane embedded in a surrounding elastic matrix. Thismodel could be used for computation of adhesive interfaces or reinforcements using fibersor fiber plates. The setting is very general and allows for both 2D and 1D models to beadded to the 3D model. Here we only consider adding membrane stiffness to an elastic 3Dmatrix, and we then use the triangulation (cid:101) T h for the discretization of three–dimensionalelasticity. In Ω \ Γ we thus assume there holds −∇ · σ ( u ) = f Ω , σ = 2 µ ε + λ Ω tr ε I , (16)for given body force f Ω , where λ Ω := Eν/ ((1 + ν )(1 − ν )) . We assume for simplicityof presentation that u = on ∂ Ω D , a part of the boundary which is assumed to include ∂ Ω h, D , and that the traction is zero on the rest of the boundary. Our finite elementmethod in the bulk is then based on the finite element space W h = (cid:110) v ∈ [ (cid:101) V h ] : v = on ∂ Ω D (cid:111) , (17)5 - Γ Figure 1: 2D representation of the problem domainand we seek u h ∈ W h such that a Ω ( u h , v ) + a h ( u h , v ) = l Ω ( v ) + l h ( v ) ∀ v ∈ W h (18)where a Ω ( u h , v ) := (2 µ ε ( v ) , ε ( w )) Ω + ( λ ∇ · v , ∇ · w ) Ω , l Ω ( v ) := ( f Ω , v ) Ω . The FEM (18) thus takes into account both the stiffness from the bulk and from themembrane. The bulk stiffness matrix is here established independently of the position ofthe membrane which allows for rapid repositioning of the membrane; this is beneficial forexample for the purpose of optimizing the membrane location. We remark that when themembrane is embedded in a three–dimensional mesh used for elasticity in all of Ω , thenwe can drop the stabilization term (or set τ = 0 ) since the three–dimensional stiffnessmatrix gives stability to the embedded membrane.In case of cohesive zone modelling at the membrane we also need to allow for in-dependent relative motion of the bulk on either side of the membrane, which requirescutting the bulk elements following [4]. In such a case we also need to reintroduce thestabilization term j h ( · , · ) . We will return to this problem in future work. This Section describes the implementational aspects of the embedded membrane modeland provides an algorithm of the implementation.6 h ΓΓ h Figure 2: Surface domain
1. Construct a mesh ˜ T in R d on the domain Ω in which the implicit surface Γ will beembedded. Let x N denote the vector of coordinates in ˜ T .2. Construct the level set function ρ ( x ) either analytically or by the use of surfacereconstruction, see Section(3.2.1) for details.3. Discretize the distance function φ = ρ ( x N ) by evaluating ρ in the nodes of thecomplete underlying mesh ˜ T .4. Find the indices to the elements in the background mesh T h , by using the discretedistance values ( φ > and φ < in some nodes of element T i , see Figure 3.5. Extract the following surface quantities. For every element T i ∈ T h a) Find the zero surface points of the polygons Γ h by looping over all elements T i ∈ T h and interpolating the discrete signed distance function φ linearly forevery element edge that has a difference in function value (see Figure 3 andSection 3.2.2 for details). For simplicity the polygon Γ h,i of element T i is splitinto triangular elements ˆ T .b) Compute the face normal n f of each triangular element, which is used tocompute the Jacobian for the basis functions of element T i . Note that caremust be taken when definining the face normal, such that the orientation ofthe normal field becomes unidirectional. The element normal n φ can be usedto orient n f in the same direction.6. Compute the displacement field u h on T h by solving linear system that results fromthe bilinear equation (11). 7 Γ ih ={ , } K i Γ Γ N N Γ x Γ Γ Γ ϕ N <0 ϕ N >0 ϕ N >0 x xxxxx e e e (a) 2D case N K i Γ Γ N N ϕ N <0 ϕ N >0 ϕ N >0 n n n ϕ =| n n + |n Γ xx x (b) Surface element normal K ihi Γ Γ Γ e e e e e e ϕ N >0 ϕ N <0 ϕ N <0 ϕ N <0 (c) 3D case Figure 3: Surface element Γ ih and parent element K i Γ in 2D and 3D7. Interpolate the solution field u h to u Γ h using the basis functions of the elements in T h . There are a number of ways to generate an implicit surface for analysis. An implicitsurface can be approximated from a CAD surface using surface reconstruction techniques,see Belytschko et al [1]. Another way is to use analytical implicit surfaces descriptions,see Burman et al [2] for an overview. In this paper we use the following analytical surfacedescriptions.Cylinder function: ρ ( x, y, z ) = (cid:112) ( x − x c ) + ( y − y c ) − r (19)where [ x c , y c ] is the center of the cylinder. See Figure 4a.Oblate spherioid: ρ ( x, y, z ) = x + y + (2 z ) − (20)See Figure 4b. The overal procedure is to use linear interpolation on the discrete interface values foreach element edge in order to find the zero level surface point, see Figure 3.1. For every element T i ∈ T h loop over all edges e i .8 a) Cylinder (b) Oblate spherioid Figure 4: Implicit surfacesa) For every edge e i check the sign of the two discrete function values φ | e i todetermine if the edge is cut.b) Linearly interpolate the cut point x Γ ,i along the edge e i using the two vertexcoordinates x me i and x ne i , at nodes m and n (endpoints of e i ) and the functionvalues φ | e i = { φ ( x me i ) , φ ( x ne i ) } .c) Let x e i , = x e i | φ | ei > (the coordinate corresponding to the highest value of φ )and x e i , = x e i | φ | ei < and compute the vector n i = x e i , − x e i , . See Figure3b2. Compute the element vector n φ = (cid:80) n i . Note that n φ points in the generaldirection of ∇ φ and is only used to determine the orientation of the face normals.3. Depending on the number of nodes in element K Γ i and the orientation of the cut,several cut cases must be considered, see Figure 5 for tetrahedral element andFigure 6 for hexahedral. A tessellation into triangles is done for all cases.4. In order to do the tessellation into triangles from an arbitrary polygon, a rotationfrom R into R was performed and then a 2D convex hull algorithm was applied. In this Section we demonstrate one particular application of the elastic membrane. Con-sider a set of membrane surfaces embedded into an elastic material body. We let theembedded membrane stiffen the solution by adding stiffness from the membrane solutionto the bulk solution. This is easily done since the membrane shares the same degrees offreedom as the bulk. 9 a) Triangular (b) Quadliteral
Figure 5: Tetrahedral cut cases (a) Triangular (b) Quadliteral (c) Pentagon (d) Hexagon
Figure 6: Hexahedral cut cases10 .3.1 Algorithm
The algorithm described here is similar to the one in Section 3.1. We allow for severalmembranes (with the possibility of different material properties) to stiffen the bulk.1. Construct a mesh ˜ T in R d on the elastic domain Ω in which the implicit surface Γ will be embedded. Let x N denote the vector of coordinates in ˜ T .2. Create a set of surfaces functions { ρ ( x ) } in the same way as in the previous algo-rithm.3. Discretize the distance functions { φ } = { ρ ( x N ) } by evaluating all functions in theset { ρ } in the nodes of the complete underlying mesh ˜ T .4. Find the sets {T h } that are intersected by the surfaces such that for each { φ } i thereexists a corresponding set of cut elements {T h } i .5. Follow the same approach as described in the previous algorithm to extract thezero surface information for each {T h } i .6. Create a discrete system of equations for the bulk elasticity problem. While as-sembling the bulk stiffness matrix, for each element that is cut by the a membranesurface, compute the membrane element stiffness and add it to the global bulk stiff-ness matrix. Note that no stabilization is needed in this case since the surroundingelements create a well conditioned stiffnessmatrix.7. Solve the linear system. In this Section we show convergence comparison on some geometries for which we cancompute the solutions analytically. We compare the convergence rates of this approachwith a triangulated surface. Numerical examples are given for both tetrahedral andhexahedral elements.The meshsize is defined by: h := 1 √ N N O where
N N O denotes the total number of nodes on the underlying 3D mesh (cid:101) T h , which isuniformly refined. Comparing this approach to the approach previously done by Hansbo and Larson [5],we consider a cylindrical shell of radius r and thickness t , with open ends at x = 0 and11 = L and with fixed axial displacements at x = 0 and radial at x = L , carrying a axialsurface load per unit area f ( x, y, z ) = F πr xL , where F has the unit of force. The resulting axial stress is σ = F (1 − ( x/L ) )4 πrt . We consider the same example as was used in [5] and choose r = 1 , L = 4 , E = 100 , ν = 1 / , t = 10 − and F = 1 .In Figure 7 we show the solution (exaggerated 10 times) on a given tetrahedral meshand in Figure 8 the corresponding solution on a hexahedral mesh with the same meshsize. In Figure 9 we show the L (Σ) error in stress, || σ − σ h || , where σ := σ Γ ( u ) and σ h := σ Γ ( u h ) . See table 1 for convergence rates.Mesh dependent errors occur in a structured tetrahedral mesh case, see Figure 10, theerror becomes less prominent with a finer mesh. This error can be avoided by using anunstructured tetrahedral mesh, see Figure 11. This error is not found in the hexahedralcase, see Figure 12. (a) Underlying linear tetrahedral mesh (b) Zero level iso-surface Figure 7: Displacement field on a tetrahedral mesh (x10 exaggerated)12 a) Underlying linear hexahedral mesh (b) Zero level iso-surface
Figure 8: Displacement field on a hexahedral mesh (x10 exaggerated) meshsize h (log) -5 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 L e rr o r ( l og ) -1.5-1-0.500.511.5 Stress error convergence
Unstructured tetrahedralTriangularHexahedralStructured tetrahedral
Figure 9: Stress convergence comparison for the cylinder.13 a) Displement field on the underlying mesh.(Exaggerated lengths) (b) Displacement field on the extracted linear zerolevel iso-surface.
Figure 10: Mesh dependent errors on low resolution linear tetrahedral mesh. (a) Underlying displacement field (b) Surface displacement field
Figure 11: Displacement fields on unstructered tetrahedral meshFigure 12: Displacement field on the underlying tri-linear hexahedral mesh mesh. (x10exaggerated) 14lementtype Mesh size || σ − σ h || RateTetrahedralunstructured 0.0899 1.7959 -0.0481 0.7733 1.34980.0260 0.4017 1.0650Tetrahedralstructured 0.1330 5.7545 -0.0721 2.1916 1.57490.0376 0.6846 1.79030.0192 0.2550 1.4722Hexahedral 0.1644 2.8832 -0.0899 1.4008 1.19540.0472 0.6698 1.14380.0242 0.3228 1.0924Original 0.1565 0.7580 -0.0783 0.3793 0.99980.0392 0.1897 1.0015Table 1: Convergence of the cylinder
Again we use the same example as in [5] an set the exact solution to be u = ( x, , and compute the stress and then the corrsponding load from (3). We set the parametersE = 1 , ν = 1 / , and t = 1 . The computed displacement field can be seen in Figure13. Compared to the previous work in [5], the superparametric stabilization method isnot needed in this approach since we already stabilize using the term j h ( · , · ) . A stressconvergence comparison can be seen in Figure 14. The convergence rates can be seen inTable 2. As one can see from the number of degree of freedoms and the correspondingmeshsize in each approach (table 2) , the approximation of the embedded approach isdependent on the surface curvature or the surface to volume ratio. If the ratio is smallthen the surface is cutting fewer elements and the resolution of the background mesh canbe kept coarse but still generate a good approximation.15igure 13: Displacement field on the extracted linear zero level iso-surface from a lineartetrahedral mesh.(x10 exaggerated) meshsize h (log) -6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 -2 L e rr o r ( l og ) -5-4.5-4-3.5-3-2.5-2 Stress error convergence
TetrahedralStructured triangular Unstructured triangular
Figure 14: Convergence of stress16lementtype Mesh size Number ofDOFs || σ − σ h || RateTetrahedral 0.0790 2414 0.1493 -0.0474 7218 0.0800 1.22220.256 26961 0.0425 1.0261Structured,triangular 0.0198 2562 0.0242 -0.0099 10242 0.0121 1.00040.0049 40962 0.0061 0.9882Unstructured,triangular 0.0198 2562 0.0288 -0.0099 10242 0.0147 0.97070.0049 40962 0.0078 0.9144Table 2: Convergence of the oblate
A rectangular box (0 ≤ x ≤ and ≤ y, z ≤ ) has a surface load, f =1, applied in thepositive x direction at x = 2 . The bulk material has the following properties, E = 100 ,ν = 0 . . The membrane has E = 1000 , ν = 0 . and t = 0 . . Figure 15 shows adisplacement plot (40x exaggerated) of the elastic bulk material without any embeddedmembrane. Figure 16 shows the same displacement plot but with 8 embedded membranesurfaces. The surfaces are visualized in Figure 17.Finally, in Figures 18 and 19 we how the effect on displacements resulting from insertinga circular membrane into a beam in bending. The material properties are the same as inthe previous example. We note the marked increase in bending stiffness resulting fromthe added membrane stiffness.Figure 15: Elastic beam without embedded elastic membrane17igure 16: Elastic beam with embedded elastic membraneFigure 17: Elastic beam with embedded elastic membraneFigure 18: Elastic beam in bending without embedded elastic membrane18igure 19: Elastic beam in bending with embedded elastic membrane In this paper we have introduced an FE model of curved membranes using higher di-mensional shape functions that are restricted to (the approximation of) the membranesurface. This allows for rapid insertion of arbitrarily shaped membranes into already ex-isting 3D FE models, to be used for example for optimization purposes. We have shownnumerically that the cut element approach to membranes gives errors comparable to tri-angulated membranes, using the same degree of approximation, and we have proposed astabilization method which provides stability to the solution as well as giving the rightconditioning of the discrete system, allowing for arbitrarily small cuts in the 3D mesh.
This research was supported in part by the Swedish Foundation for Strategic ResearchGrant No. AM13-0029, the Swedish Research Council Grants Nos. 2011-4992 and 2013-4708, and the Swedish Research Programme Essence19 eferences [1] T. Belytschko, C. Parimi, N. Moës, N. Sukumar, and S. Usui. Structured extendedfinite element methods for solids defined by implicit surfaces.
Internat. J. Numer.Methods Engrg. , 56(4):609–635, 2003.[2] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing. CutFEM: Discretiz-ing geometry and partial differential equations.
Internat. J. Numer. Methods Engrg. ,104(7):472–501, 2015.[3] E. Burman, P. Hansbo, and M. G. Larson. A stabilized cut finite element method forpartial differential equations on surfaces: the Laplace-Beltrami operator.
Comput.Methods Appl. Mech. Engrg. , 285:188–207, 2015.[4] A. Hansbo and P. Hansbo. A finite element method for the simulation of strongand weak discontinuities in solid mechanics.
Comput. Methods Appl. Mech. Engrg. ,193(33-35):3523–3540, 2004.[5] P. Hansbo and M. G. Larson. Finite element modeling of a linear membrane shellproblem using tangential differential calculus.
Comput. Methods Appl. Mech. Engrg. ,270:1–14, 2014.[6] P. Hansbo and K. Salomonsson. A discontinuous Galerkin method for cohesive zonemodelling.
Finite Elem. Anal. Des. , 102/103:1–6, 2015.[7] M. A. Olshanskii, A. Reusken, and J. Grande. A finite element method for ellipticequations on surfaces.
SIAM J. Numer. Anal. , 47(5):3339–3358, 2009.[8] O. C. Zienkiewicz.
The Finite Element Method . McGraw-Hill, London, 1977. Thethird, expanded and revised, edition of