A C0 interior penalty finite element method for flexoelectricity
AA C0 interior penalty finite element method forflexoelectricity
J. Ventura ,a , D. Codony ,b , S. Fern´andez-M´endez ,c, ∗ Laboratori de Clcul Numric, E.T.S. de Ingeniera de Caminos,Universitat Politcnica de Catalunya, Barcelona, Spain a Email: [email protected] b Email: [email protected] c Email: [email protected] ∗ Corresponding author
Abstract
We propose a C Interior Penalty Method (C0-IPM) for the computational modelling offlexoelectricity, with application also to strain gradient elasticity, as a simplified case. Standardhigh-order C finite element approximations, with nodal basis, are considered. The proposedC0-IPM formulation involves second derivatives in the interior of the elements, plus integralson the mesh faces (sides in 2D), that impose C continuity of the displacement in weak form.The formulation is stable for large enough interior penalty parameter, which can be estimatedsolving an eigenvalue problem. The applicability and convergence of the method is demon-strated with 2D and 3D numerical examples. Keywords: C finite elements , interior penalty method , strain gradient elas-ticity , flexoelectricity The rising interest on microtechnology evidences the need for mathematical and computationalmodels suitable for small scales, often giving rise to th order Partial Differential Equations (PDEs).In particular flexoelectric effects become relevant, and may be crucial, in the design of smallelectromechanical devices or for the understanding of physical phenomena [1]. The modellingof flexoelectricity involves a two-way coupling between strain gradient and electric field; straingradient elasticity is frequently also included in the model to regularise the problem, leading to asystem of th order PDEs. 1 a r X i v : . [ m a t h . NA ] A ug everal numerical strategies have been recently proposed for the solution of flexoelectricityproblems, based on the use of C approximation spaces or on mixed formulations. Mixed formula-tions split the PDE in two nd order PDEs, allowing the use of C Finite Element (FE) approxima-tions [2, 3]. The approximation spaces, for the primal unknown and for the additional unknowns,must fulfill some conditions for stability that lead to approximation spaces with cumbersome defi-nitions, and difficult extension to 3D or high-order approximations. However, the main drawbackof mixed formulations in the high computational cost due to the additional unknowns.On other hand, C approximations can be directly used for the discretization of the weak formin H , involving nd order derivatives, without additional unknowns. The first successful attemptin this direction considered a meshless method: the maximum entropy method [4]. Unfortunately,the computational cost of accurate meshless methods is high, mainly due to the excessive numberof integration points for accurate solutions, and the large stencils in the discrete matrices. Aimingto improve the efficiency, a solution based on Isogeometric Analysis (IGA) is proposed in [5]. Thesolution is approximated by means of Non-Uniform Rational B-Splines (NURBS). An interestingcritical comparison of IGA, meshless and mixed methods for flexoelectricity can be found at [6],although numerical examples are restricted to simple geometries. The conclusion is that IGA isvery efficient on regular grids, corresponding to a transformation of a rectangle grid, where plainB-Spline approximations can be easily defined. However, in a more general context, defininga NURBS approximation with C continuity in a whole domain with complex shape may notbe straightforward, and the numerical integration of the resulting NURBS may be, again, veryexpensive [7]. An efficient alternative for complex domains is the immersed B-Spline methodproposed in [8]. It considers B-Spline approximations based on a background regular grid, withan embedded domain. The applicability of the proposal is demonstrated with 2D and 3D complexgeometries. The weak points of immersed B-Splines are the usual ones in the context of embeddeddomains: the robust definition of numerical integration in cropped cells (intersected by the domainboundary), which is specially challenging in 3D, and the ill-conditioning problems in the presenceof cells with a small portion in the domain, that can be alleviated with specific techniques [9].Here we propose a C Interior Penalty Method (C0-IPM) for the solution of flexoelectricity.A standard C FE approximation is considered, and C continuity between elements is imposed inweak form by means of the Interior Penalty Method (IPM). The procedure for the derivation ofthe C0-IPM weak form is analogous to the derivation of the IPM in the context of DG methods for nd order PDEs [10], or Nitsche’s method for weak imposition of Dirichlet boundary conditions[11], but now applied to the continuity of normal derivatives on element boundaries. The resultingweak form involves second derivatives in the interior of the elements, plus integrals on the faces(sides in 2D) of the mesh, that impose C continuity of the displacement in weak form. C0-IPMformulations overcome the disadvantages of other methods, because they allow the use of standard C FE approximations. Namely, (i) the computational mesh can be adapted to any geometry, withlocalised refinement were needed, (ii) there is no need to use embedded discretizations, avoidingthe consequent ill-conditioning problems and the definition of special numerical integration forcropped elements, (iii) there are no additional unknowns and (iv) they handle material interfacesin a natural way. In summary, C0-IPM retains the computational efficiency and the versatility that2ake standard FEs the preferred method for many practitioners in the computational mechanicscommunity.In [12], C0-IPM formulations, there referred to as continuous/discontinuous finite elements,are applied to several problems modelled by th order PDEs, including Kirchhoff plates and 1Dstrain gradient elasticity. Numerical experiments show the applicability of the formulation in bothapplications, but convergence studies are limited to 1D examples. The C0-IPM formulation isthen analysed in [13] for the 2D biharmonic equation, with first and second Dirichlet conditions,including a convergence analysis that shows that the method is convergent for p ≥ , but may havesuboptimal convergence depending on the degree and the penalty parameter. An experimentalconvergence study for Kirchhoff plates can be found at [14]. The numerical results demonstratethe applicability of the method for degree greater or equal to , and also show slightly suboptimalconvergence, that slowly deteriorates for larger penalty parameter, in agreement with the analysisin [13]. Variations of C0-IPM have also been applied to strain gradient dependent damage modelsin [15] and to the Cahn-Hilliard equation in [12, 16].This paper develops the C0-IPM method for flexoelectricity and, as a simplified case, for straingradient elasticity, for 2D and 3D computations. Section 2 presents the problem statement andrecalls the weak form in H . The derivation of the C0-IPM method for C approximations, notin H , is presented in section 3. An eigenvalue problem to determine a large enough penaltyparameter, ensuring coercivity of the strain gradient bilinear form, is derived in section 3.1. Finally,in section 4, 2D and 3D numerical experiments demonstrate the applicability of the method andshow, as expected, slightly suboptimal convergence under uniform mesh refinement, but still witha robust high-order convergence for p ≥ . We consider the model in [8], where flexoelectricity is ruled by the following set of PDEs andboundary conditions: ∇ · ( (cid:98) σ ( u , φ ) − ∇ · (cid:101) σ ( u , φ )) + b = in Ω (1a) ∇ · (cid:98) D ( u , φ ) − q = 0 in Ω (1b) u = g on Γ u D (1c) t ( u , φ ) = t n on Γ u N (1d) ∂ u ∂ n = g on Γ u D (1e) r ( u , φ ) = r n on Γ u N (1f) j ( u , φ ) = j ext on C ∂ Ω N (1g) φ = g on Γ φD (1h) w ( u , φ ) = w n on Γ φN (1i)3here Ω ⊂ R n sd is the domain, the displacement u and the electric potential φ are the unknowns, (cid:98) σ and (cid:101) σ are the local and double stress tensors and (cid:98) D is the electric displacement tensor, that is, (cid:98) σ = C : ε − E · e ≡ (cid:98) σ ij = C ijk(cid:96) ε k(cid:96) − E (cid:96) e (cid:96)ij (cid:101) σ = h ... ∇ ε − E · µ ≡ (cid:101) σ ijk = h ijk(cid:96)mn ∂ε (cid:96)m ∂x n − E (cid:96) µ (cid:96)ijk (cid:98) D = κ · E + e : ε + µ ... ∇ ε ≡ (cid:98) D (cid:96) = κ (cid:96)m E m + e (cid:96)ij ε ij + µ (cid:96)ijk ∂ε ij ∂x k ε is the strain tensor, that is ε ij = ( ∂u i /∂x j + ∂u j /∂x i ) / , E = − ∇ φ is the electric field, C is the elasticity tensor, that depends on the Young modulus E and the Poisson ratio ν , h is thestrain-gradient tensor, defined as h ijk(cid:96)mn = l C ij(cid:96)m δ kn with the internal length scale parameter l , e and µ are the tensors of piezoelectric and flexoelectric coefficients and κ contains the dielectricityconstants, see appendix B in [8] for detailed definitions.In the previous equations, and in the rest of the document, Einstein’s notation is assumed. Thatis, repeated indexes sum over the spatial dimensions.The boundary of the domain is split in Dirichlet and Neumann boundaries, for the first andsecond conditions of the mechanical problem and for the electric problem, that is ∂ Ω = Γ u D ∪ Γ u N = Γ u D ∪ Γ u N = Γ φD ∪ Γ φN . Note that all volume and boundary domains are assumed to be open domains, not including theirboundaries.The first mechanical boundary condition, (1c) or (1d), sets the displacement u or the traction t i ( u , φ ) = (cid:32)(cid:98) σ ij − ∂ (cid:101) σ ijk ∂x k − ∇ Sk (cid:101) σ ikj (cid:33) n j + (cid:101) σ ijk (cid:101) N jk , where ∇ Sk (cid:101) σ ikj is the surface divergence of (cid:101) σ ikj , and ˜ N is the second order geometry tensor, see [8]for details. The second mechanical boundary condition, (1e) or (1f), sets the normal derivative ofthe displacement ∂ u /∂ n or the double traction r i ( u , φ ) = (cid:101) σ ijk n j n k . The condition (1g) sets forces on the Neumann boundary edges. That is, the domain boundaryis assumed to be composed of smooth surfaces (curves in 2D) that are joined on sharp boundaryedges (corners in 2D). C ∂ Ω N denotes the union of the boundary edges that are shared by two surfaceswith first Neumann conditions, i.e. the edges in the interior of Γ u N . At the edges shared by at leastone Dirichlet surface the value on the edge is assumed to be the one set on the surface, i.e. u = g for all edges in Γ u D . Line forces (punctual forces in 2D) are defined on boundary edges as j i ( u , φ ) = τ Lj (cid:101) σ Lij(cid:96) n L(cid:96) + τ Rj (cid:101) σ Rij(cid:96) n R(cid:96) , j ( u , φ ) . The superscripts L and R refer to the face (side in 2D) sharing the edge (corner in 2D).being n L and n R the unitary exterior normals on the left and right surfaces sharing the boundaryedge, and τ L and τ R the tangent vectors on each surface pointing outward and perpendicular tothe edge, see an example in figure 1 left. In 2D, τ L and τ R at a corner are just the tangent vectorson each curve sharing the corner and pointing outward, as depicted in figure 1 right.Finally, the electric boundary condition, (1h) or (1i), sets the electric potential φ or the surfacecharge density w ( u , φ ) = − D (cid:96) ( u , φ ) n (cid:96) . Remark 1
For the sake of simplicity, we initially restrict to the case with second Neumann bound-ary conditions in the whole boundary, i.e. Γ u N = ∂ Ω and Γ u D = ∅ . The treatment of secondDirichlet conditions (1e) is commented in Remark 2. It is also worth mentioning that the ra-tionales in this work can also be applied to other models, including converse flexoelectricity orexpressed in terms of polarization, see for instance [6]. If an approximation in H (Ω) can be considered, multiplying (1a) by a weighting vector v ,applying integration by parts twice, and using the symmetries of the stress tensors, leads to (cid:90) Ω v · b d Ω = (cid:90) Ω ε ( v ) : (cid:98) σ ( u , φ ) d Ω + (cid:90) Ω ∇ ε ( v ) ... (cid:101) σ ( u , φ ) d Ω − (cid:90) ∂ Ω v i (cid:32)(cid:98) σ ij − ∂ (cid:101) σ ijk ∂x k (cid:33) n j dS − (cid:90) ∂ Ω ∂v i ∂x j (cid:101) σ ijk n k dS , where ∇ ε ... (cid:101) σ = ∂ ε ij /∂x k (cid:101) σ ijk .Now, to properly treat boundary conditions, the derivative ∂v i /∂x j on the boundary is split innormal and tangential derivatives, and the surface divergence theorem is applied to the term with5he tangential derivative, leading to (cid:90) Ω v · b d Ω = (cid:90) Ω ε ( v ) : (cid:98) σ ( u , φ ) d Ω + (cid:90) Ω ∇ ε ( v ) ... (cid:101) σ ( u , φ ) d Ω − (cid:90) ∂ Ω v · t ( u , φ ) dS − (cid:90) ∂ Ω ∂ v ∂ n · r ( u , φ ) dS − (cid:90) C ∂ Ω v · j ( u , φ ) dl , (2)where C ∂ Ω is the union of all sharp edges of the domain, and the integral on it reduces to a sumevaluating at the boundary corners in 2D.Thus, applying boundary conditions (1c)-(1i), under the assumption Γ u D = ∅ , and adding theweighted residual of the electric potential problem (1b) with (1h) and (1i), the weak form of (1) in H (Ω) is: find u ∈ [ H (Ω)] n sd and φ ∈ H (Ω) such that (1c) and (1h) hold and (cid:90) Ω ε ( v ) : (cid:98) σ ( u , φ ) d Ω + (cid:90) Ω ∇ ε ( v ) ... (cid:101) σ ( u , φ ) d Ω + (cid:90) Ω ∇ ψ · (cid:98) D ( u , φ ) d Ω= s ( v ) − (cid:90) Ω ψq d Ω − (cid:90) Γ φN ψ w n d S (3)for all v ∈ [ H (Ω)] n sd and ψ ∈ H (Ω) such that v = on Γ u D and ψ = 0 on Γ φD where s ( v ) = (cid:90) Γ u N v · t n dS + (cid:90) Γ u N ∂ v ∂ n · r n dS + (cid:90) C ∂ Ω N v · j ext dl + (cid:90) Ω v · b d Ω . (4)This weak form is not suitable when considering C FE approximations, but the same derivationcan be applied in the interior of each element of the mesh, as detailed next. C Interior Penalty Finite Element method
The domain Ω is now split in FEs { Ω e } n el e =1 , and a C piece-wise polynomial approximation isconsidered. That is, the approximation space for the components of the displacement and for thepotential is V h = { v ∈ H (Ω) such that ϕ − e ( v | Ω e ) ∈ P p for e = 1 , . . . , n el } (cid:54)⊂ H (Ω) , where ϕ e is the isoparametric transformation from the reference element to the physical element Ω e , and P p is the space of polynomials of degree less or equal to p for simplexes, and less or equalto p in each direction for quadrilaterals and hexahedra.Since the approximation space in not in H (Ω) , we can not consider the weak form (3). How-ever, the approximation is H (Ω e ) ; thus, considering (2) in each element we have (cid:90) Ω e v · b d Ω = (cid:90) Ω e ε ( v ) : (cid:98) σ ( u , φ ) d Ω + (cid:90) Ω e ∇ ε ( v ) ... (cid:101) σ ( u , φ ) d Ω − (cid:90) ∂ Ω e v · t e ( u , φ ) dS − (cid:90) ∂ Ω e ∂ v ∂ n e · r e ( u , φ ) dS − (cid:90) C ∂ Ω e v · j e ( u , φ ) dl , (5)6here C ∂ Ω e is the union of the edges (corners in 2D) of the element Ω e and j e is the line force(punctual force in 2D) on C ∂ Ω e ; see a representation in figure 2a. The superscripts highlight thatthe surface and line forces, and the normal vector, are from the element Ω e .Summing for all elements, and noting that v is continuous but ∂ v /∂ n is not, we get (cid:90) Ω v · b d Ω = (cid:90) Ω ε ( v ) : (cid:98) σ ( u , φ ) d Ω + (cid:90) (cid:98) Ω ∇ ε ( v ) ... (cid:101) σ ( u , φ ) d Ω − (cid:90) I v · ( t L ( u , φ ) + t R ( u , φ )) dS − (cid:90) I (cid:32) ∂ v L ∂ n L · r L ( u , φ ) + ∂ v R ∂ n R · r R ( u , φ ) (cid:33) dS − n edg (cid:88) k =1 (cid:90) C k v · (cid:88) e ∈ E ( k ) j e ( u , φ ) dl − (cid:90) ∂ Ω v · t ( u , φ ) dS − (cid:90) ∂ Ω ∂ v ∂ n · r ( u , φ ) dS , (6)where I is the union of all the internal element faces and (cid:98) Ω is the union of the interior of theelements, where second derivatives are well-defined, i.e. (cid:98) Ω = n el (cid:91) e =1 Ω e , I = n el (cid:91) e =1 ∂ Ω e \ ∂ Ω , see figure 2b. The supercripts R and L now denote the evaluation from the elements to the leftand right side of the face in I (see figure 2c), and { C k } n edg k =1 are all the edges (corners in 2D) in themesh, being E ( k ) the set of indexes of the elements sharing the edge C k .Now, let us recall that the conditions for interfaces in the domain (also in the case of discon-tinuous material parameters) are the ones corresponding to both continuity for the Dirichlet valuesand equilibrium of Neumann forces. That is, (cid:74) u ⊗ n (cid:75) = (7a) (cid:74) t ( u , φ ) (cid:75) = (7b) (cid:116) ∂ u ∂ n (cid:124) = (7c) (cid:74) r ( u , φ ) ⊗ n (cid:75) = (7d)on the faces in I , where the jump operator is defined as (cid:74) a (cid:75) = a L + a R (a) (b) (c) Figure 2: Example of 2D discretization: (a) interior corner C k (i.e. interior mesh vertex), as seenfrom element Ω e , and representation of the normal and tangent vectors corresponding to the leftand right side of the element sharing the corner, for the definition of the punctual force j e at thecorner, (b) interior faces I in white and broken domain (cid:98) Ω in blue, (c) normal vectors on one faceshared by its left and right elements, for the computation of the jump on a face.and it is used always involving a change of sign due to an odd appearance of the normal vector.In addition, we have to impose equilibrium of forces on the mesh edges. That is, each element Ω e contributes with a force j e on its edges (corners in 2D); and for each edge C k , the sum of theforces for all elements sharing the edge (i.e. for Ω e with e ∈ E ( k ) ) must be zero, or in internalequilibrium with the external forces. That is, (cid:88) e ∈ E ( k ) j e ( u , φ ) = (cid:26) on C k (cid:54)⊂ C ∂ Ω j ext on C k ⊂ C ∂ Ω N (8)where j ext is the force set in (1g). Note that { C k (cid:54)⊂ C ∂ Ω } includes interior edges and also elementedges on ∂ Ω , just excluding the ones in the domain sharp edges. For the edges in Γ u D no valueis set, and the sum of the forces will be in equilibrium with the reaction forces associated to theprescribed displacement (1c).On other hand, using the algebraic identity ( a L n L ) b L + ( a R n R ) b R = { a } (cid:74) b n (cid:75) + (cid:74) a n (cid:75) { b } ,and the equilibrium condition (7d), we can rewrite ∂ v L ∂ n L · r L ( u , φ ) + ∂ v R ∂ n R · r R ( u , φ )= { ∇ v } : (cid:74) r ( u , φ ) ⊗ n (cid:75) + (cid:116) ∂ v ∂ n (cid:124) · { r ( u , φ ) } = (cid:116) ∂ v ∂ n (cid:124) · { r ( u , φ ) } , (9)with the mean operator { a } = ( a L + a R ) / . 8ow, replacing in (6) the identity (9), the Neumann boundary conditions (1d) and (1f), thehomogeneous Dirichlet condition v = on Γ u D related to (1c), the first interface equilibriumcondition (7b) and the equilibrium at interior edges (8), and under the assumption Γ u D = ∅ , (6)simplifies to (cid:90) Ω ε ( v ) : (cid:98) σ ( u , φ ) d Ω + (cid:90) (cid:98) Ω ∇ ε ( v ) ... (cid:101) σ ( u , φ ) d Ω − (cid:90) I (cid:116) ∂ v ∂ n (cid:124) · { r ( u , φ ) } dS = s ( v ) , (10)with s ( v ) defined in (4).The first two integrals in (10) are symmetric and coercive bilinear forms in v and u , as expectedfor the weak form of a strain gradient elasticity operator. However, it is not the case for the integralon the interior faces I .The idea of IPM is adding terms that are analytically zero, thanks to the continuity interfacecondition (7c), to recover symmetry and coercivity of the strain gradient bilinear form. The result-ing weak form for flexoelectricity, under the assumption Γ u D = ∅ , is: find u ∈ [ H (Ω) ∩ H ( (cid:98) Ω)] n sd and φ ∈ H (Ω) such that (1c) and (1h) hold and (cid:90) Ω ε ( v ) : (cid:98) σ ( u , φ ) d Ω + (cid:90) (cid:98) Ω ∇ ε ( v ) ... (cid:101) σ ( u , φ ) d Ω + (cid:90) Ω ∇ ψ · (cid:98) D ( u , φ ) d Ω − (cid:90) I (cid:116) ∂ v ∂ n (cid:124) · { r ( u , φ ) } dS − (cid:90) I { r ( v , ψ ) } · (cid:116) ∂ u ∂ n (cid:124) dS + (cid:90) I β (cid:116) ∂ v ∂ n (cid:124) · (cid:116) ∂ u ∂ n (cid:124) dS = s ( v ) − (cid:90) Ω ψq d Ω − (cid:90) Γ φN ψ w n d S (11)for all v ∈ [ H (Ω) ∩ H ( (cid:98) Ω)] n sd and ψ ∈ H (Ω) such that v = on Γ u D and ψ = 0 on Γ φD .The parameter β is a stabilization parameter that must be taken large enough to ensure coerciv-ity of the strain gradient bilinear form, to get a well-defined saddle point problem [8]. Although itis usually called penalty parameter, thanks to the consistency of the formulation, high-order con-vergence can be achieved with β or order h − . In practice, not very large values are needed foraccurate solutions, avoiding the unaccuracy or ill-conditioning that typically suffer non-consistentpenalty methods [17].The minimum value of the stabilization parameter β can be estimated solving an eigenvalueproblem, as commented in section 3.1. Remark 2 (Second Dirichlet conditions)
If second Dirichlet boundary conditions (1e) are im-posed (i.e. Γ u D (cid:54) = ∅ ), an additional term (cid:82) Γ u D ∂ v /∂ n · r ( u , φ ) dS appears in (10) and, con-sequently, in the weak form (11) . Following the same IPM rationale, two new terms, that are ull thanks to (7c) , are also added in (11) to recover again symmetry and coercivity, namely (cid:82) Γ u D r ( v , ψ ) · ( ∂ u /∂ n − g ) dS + β D (cid:82) Γ u D ∂ v /∂ n · ( ∂ u /∂ n − g ) dS, where β D is a new sta-bilization parameter that can be taken equal to β or tuned separately. Remark 3
It is interesting to note that the C0-IPM weak form (11) reduces to the one for H (Ω) ,i.e. (3) , when a C (Ω) approximation is considered. The C0-IPM formulation keeps the consistencyand is valid for standard C FE approximations, just introducing the proper integrals on the faces I . Also note that the second integral in (11) is in the interior of the elements, (cid:98) Ω , to account for thefact that second derivatives are not defined on I . β In this section we derive an eigenvalue problem to estimate a lower bound for β . The derivation isthe usual one in IPM and Nitsche’s formulations [18, 8, 14].The bilinear form of the strain-gradient elasticity operator is A ( u , v ) = a ( u , v ) − (cid:90) I (cid:116) ∂ v ∂ n (cid:124) · { r sg ( u ) } dS − (cid:90) I { r sg ( v ) } · (cid:116) ∂ u ∂ n (cid:124) dS + β (cid:90) I (cid:116) ∂ v ∂ n (cid:124) · (cid:116) ∂ u ∂ n (cid:124) dSwith a ( u , v ) = (cid:90) (cid:98) Ω ∇ ε ( v ) ... (cid:101) σ sg ( u ) d Ω where r sgi = (cid:101) σ sgijk n j n k and (cid:101) σ sg ( u ) = h ... ∇ ε ( u ) , that is, the mechanical part of the second tractionand the double stress tensor.The bilinear form a is semicoercive (i.e. a ( u , u ) > for any u such that ∇ ε ( u ) (cid:54) = ,and a ( u , u ) = 0 otherwise), leading to well-posed strain gradient elasticity and flexoelectricityproblems for any value of the internal length scale parameter l . However, the addition of theintegrals on the faces I , leads to a bilinear form A that retains semicoercivity only for large enough β . Thus, to ensure well-posedness of the discrete problem for any value of l , we want β such that A ( u , u ) > ∀ u ∈ (cid:101) U h = { u ∈ [ V h ] n sd such that ∇ ε ( u ) (cid:54) = } with A ( u , u ) = a ( u , u ) − (cid:90) I (cid:116) ∂ u ∂ n (cid:124) · { r sg ( u ) } dS + β (cid:90) I (cid:116) ∂ u ∂ n (cid:124) · (cid:116) ∂ u ∂ n (cid:124) dS . − (cid:90) I (cid:116) ∂ u ∂ n (cid:124) · { r sg ( u ) } dS + β (cid:90) I (cid:116) ∂ u ∂ n (cid:124) · (cid:116) ∂ u ∂ n (cid:124) dS ≥ − (cid:107){ r sg ( u ) }(cid:107) L ( I ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:116) ∂ u ∂ n (cid:124) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I ) + β (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:116) ∂ u ∂ n (cid:124) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I ) ≥ − (cid:15) (cid:107){ r sg ( u ) }(cid:107) L ( I ) + ( β − (cid:15) ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:116) ∂ u ∂ n (cid:124) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I ) , for any positive (cid:15) .Thus, considering a positive constant K such that (cid:107){ r sg ( u ) }(cid:107) L ( I ) ≤ Ka ( u , u ) ∀ u ∈ (cid:101) U h , (12)we have A ( u , u ) ≥ (cid:18) − K(cid:15) (cid:19) a ( u , u ) + ( β − (cid:15) ) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:116) ∂ u ∂ n (cid:124) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) L ( I ) and the bilinear form is then positive definite if − K/(cid:15) > and β − (cid:15) > , for any positive (cid:15) .In conclusion, the strain gradient bilinear form A is positive definite in the reduced discretespace if β > K , where K is the constant satisfying (12). This constant can be computed as thelargest eigenvalue of the generalised problem Bx = λ Ax where B and A are the discrete matrices corresponding to the bilinear forms b ( u , v ) and a ( u , v ) in the reduced discrete space (cid:101) U h , with b ( u , v ) = (cid:90) I { r sg ( v ) } · { r sg ( u ) } dS . The computation of the maximum eigenvalue in the reduced space (cid:101) U h can be done from theproblem stated in the complete discrete space [ V h ] n sd setting nodal values to reduce the space orusing the so-called eigenvalue problem deflation [19]. Remark 4
Matrices B and A scale as O ( E l h n sd − ) and O ( El h n sd − ) , respectively, with char-acteristic element size h . Thus, the maximum eigenvalue of (12) scales as O ( El /h ) . Conse-quently, we can consider β = αEl /h, (13)11 ith a large enough constant α , that can be computed solving the eigenvalue problem, or simplytuned, in a coarse mesh with any value of E and l .It is important noting that, differently to non-consistent penalty methods, IPM and Nitsche’smethods provide accurate solutions and high-order convergence with moderate values of β of order O ( h − ) for any degree of approximation, as shown in the numerical examples in section 4 and in[17]. An alternative sufficient condition to have a well-posed discrete problem can be stated includ-ing also the elasticity term in the bilinear form, that is, with a ( u , v ) = (cid:90) Ω ε ( v ) : (cid:98) σ sg ( u ) d Ω + (cid:90) (cid:98) Ω ∇ ε ( v ) ... (cid:101) σ sg ( u ) d Ω , where (cid:98) σ sg ( u ) = C : ε ( u ) is the mechanical part of the local stress tensor. This option leadsto a smaller (sharper) bound for β , specially for small l or large h . However, since the matrixcorresponding to the first elasticity term scales as O ( Eh n sd − ) , the dependency on the mesh sizeand material parameters is not so obvious. The current implementation considers high-order Lagrange nodal basis, with Fekete nodes in thereference element to minimise the condition number of elemental matrices. For degree p ≥ I . To do so, the standard C reference element isextended including second derivatives of the basis functions at the element integration points, thevalue of element basis functions and their derivatives at the integration points of the reference ele-ment faces, a list of the nodes corresponding to each face in the reference element and permutationsfor the integration points of the reference face for flipping .The so-called flipping is a permutation (usually for the nodes in DG methods, but for integrationpoints in our IPM implementation) that has to be applied to the face when seen from the secondelement, to match the orientation of the corresponding face in the first element. In 2D the flippingis the same for any side of the mesh, just using a reverse ordering for the second element sharingthe side. In 3D the possible rotations of the face have to be taken into account to choose the properpermutation for the integration points.A variable storing, for each face, the number of the elements sharing the face, the local num-bering of the face in each one of the two elements and the rotation to be applied for the secondelement, is also computed from the mesh as a preprocess.Dirichlet conditions (1c) are imposed in strong form, just setting the corresponding nodal val-ues, as usual in standard FE computations. Second Dirichlet conditions can be imposed in weakform as commented in Remark 2. 12 Numerical examples
Several numerical examples are included in this section to study the convergence of the C0-IPMformulation in 2D and 3D, and to validate the computational tool by comparison with previousworks. Homogeneous first, second and corner Neumann boundary conditions are assumed whereno boundary condition is specified.
The convergence of the method for the solution of problem (1) is studied in this section. To test themethod with non-regular meshes and curved boundaries, the problem is solved in a square with ahole,
Ω = (0 , \ B ((0 . , . , . . Figure 3 shows the coarsest mesh for nested refinement, withdegree p = 4 .First Dirichlet and second Neumann conditions are imposed on all the boundary. The bodyforce b , the free charge q and the boundary data are set so that the solution is u = [sin (2 π ( x + x )) , cos (2 π ( x + x ))] T ,φ = sin(2 π ( x + x )) + cos(2 π ( x + x )) , (14)the material parameters are E = 2 . , ν = 0 . ,l = 1 . , κ L = 1 . ,e L = 7 . , e T = 1 . , e S = 1 . ,µ L = 1 . , µ T = 1 . , µ S = 5 . , (15)Figure 3: 2D convergence test: initial mesh for the nested refinement, with degree p = 4 .13 log (h) -8-6-4-20 l og ( L2 e rr o r u ) u log (h) -8-6-4-20 l og ( L2 e rr o r ) Figure 4: 2D convergence test for the uncoupled problem ( strain gradient elasticity and potentialequations) with β = 100 El /h : L error under nested mesh refinement (the coarsest mesh isshown in figure 3), for degree of approximation p = 1 . . . . The numbers are the slopes for eachsegment. -2 -1.5 -1 -0.5 log (h) -8-6-4-20 l og ( L2 e rr o r u ) u log (h) -6-4-202 l og ( L2 e rr o r ) -0.5-0.2-0.1 0.91.21.5 3.33.13.1 3.93.94.0 Figure 5: 2D convergence test for the coupled flexoelectricity problem with β = 100 El /h .and the piezoelectric principal direction is x . The definition of the material tensors in terms ofthese parameters can be found, for instance, in appendix B of [8].First we consider the uncoupled problem (with e = and µ = ), that is, an uncoupledsolution of a strain gradient elasticity problem and an electric potential problem. The convergenceplots are shown in figure 4, for penalty parameter β = 100 El /h , and degree p = 1 . . . . For straingradient elasticity, the displacement error behaves in agreement with the results for Kirchhoffplates in [14]. With degree p = 1 , the approximation space is not rich enough to impose C continuity. Moreover, the second derivatives of the displacement in the strain gradient elasticityterms and the flexoelectricity terms cancel out, or are almost zero for curved elements. Thus, the14ethod does not converge for linear approximation. For degree p = 2 much finer meshes would benecessary to reach assymptotic convergence, reducing its practical applicability. Accurate resultswith high-order convergence are obtained for degree p ≥ , with slightly suboptimal convergencefor p = 3 , in agreement with the analysis in [13] for the biharmonic equation. In this particularexample, p = 4 behaves better than expected, exhibiting slightly superoptimal convergence. Theexpected optimal convergence is observed for the uncoupled electric potential problem for anydegree.Figure 5 shows the convergence plots for the flexoelectricity problem, with piezoelectric andflexoelectric coupling. The coupling leads to a reduction in the convergence rate, not relevant forthe displacement, but around one for the potential, for p ≥ . This is probably due to the, smallbut still present, discontinuity of the displacement derivative across element sides, affecting thepotential through the flexoelectricity coupling.The conclusion is then that, even though convergence is suboptimal, the method is able to reachhigh accuracy with high-order convergence for degree p ≥ . The C0-IPM method is thereforepromising for an efficient solution of flexoelectricity.Similar results can be observed with quadrilateral meshes, with better behaviour for the p = 2 approximation thanks to the richer approximation space and the presence of interior nodes in theelement. β The effect of the interior penalty parameter in the accuracy of the numerical solution is studiednext, with the 2D example and meshes of the previous section. Following Remark 4, the parameteris taken as (13), with different orders of magnitude for α , independent of h .Figure 6 shows the convergence plots for the flexoelectricity coupled problem, for the displace-ment u (left) and for the potential φ (right), for degree p = 3 (top) and p = 4 (bottom). The slopesof the segments are shown for the plots with α = 10 for p = 3 , and with α = 100 for p = 4 .We can observe the poor performance of the method for α = 1 , due to the fact that it is not largeenough for a coercive mechanical bilinear form.For degree p = 3 , α = 10 is large enough and provides the best results. Larger values of α ,several orders of magnitude larger, also lead to high-order convergence, proving the robustness ofthe method; but, in agreement with the analysis in [13] the convergence rate slowly decreases forincreasing α .Looking to the results for p = 4 we can observe that, for α = 10 , the bilinear form is coercivefor the first meshes, because the elasticity part dominates in the coefficients of the matrix. This isnot the case for the last mesh, where higher order terms become more relevant. With α ≥ thecondition in Remark 4 is satisfied and convergence is close to p + 1 = 5 for u and around p = 4 for φ , with almost no loss in the accuracy for increasing β .Thus, from this experiment we conclude that C0-IPM with degree p = 4 provides excellentresults, with convergence rates close to p + 1 = 5 for u and around p = 4 for φ , and with littledependency on the particular value of β , for β ≥ El /h .15 log (h) -5-4-3-2-10 l og ( L2 e rr o r u ) u (p=3) log (h) -4-3-2-101 l og ( L2 e rr o r ) (p=3) log (h) -8-6-4-20 l og ( L2 e rr o r u ) u (p=4) =10 -2 -1.5 -1 -0.5 log (h) -6-4-202 l og ( L2 e rr o r ) (p=4) Figure 6: Effect of the β parameter in the solution of the flexoelectricity problem: convergenceplots for the example in section 4.1 with β = αEl /h and different values of α , for degree p = 3 (top) and p = 4 (bottom).The same analysis is performed now for strain gradient elasticity. Figure 7 shows the conver-gence plots for the displacement u for degree p = 3 (left) and degree p = 4 (right), with the sameconclusions. The cantilever beam depicted in figure 8 is considered. The aspect ratio is 20, and the width a variesto show the size-dependent nature of flexoelectricity. The beam is fixed to a wall and grounded onits left end, and it undergoes a punctual force F at the top-right corner. The boundary conditions16 log (h) -5-4-3-2-10 l og ( L2 e rr o r u ) u (p=3) log (h) -8-6-4-20 l og ( L2 e rr o r u ) u (p=4) =10 Figure 7: Effect of the β parameter in the solution of the strain gradient elasticity problem: con-vergence plots for the displacement, for the example in section 4.1, with β = αEl /h and differentvalues of α , for degree p = 3 (left) and p = 4 (right). L x aF Figure 8: Cantilever beam under bending and open-circuit boundary conditions.are thus u = at x = 0 j ( u , φ ) = − F at x = ( L, a/ φ = 0 at x = L, (16)where L = 20 a is the beam length. To reproduce the results obtained in [8] with B-splines, thematerial parameters are E = 100 GPa , , κ = κ = 11 nJ V − m − ,e T = − . − m − , µ T = 1 µ J V − m − ,l = ν = µ L = µ S = e L = e S = 0 , (17)and the piezoelectric principal direction is x . A uniform discretization with × × triangularelements (with characteristic element size h = 0 . a ) of degree p = 4 , and with β = 100 , isconsidered. Since l = 0 , any positive value of β provides good results.17 Figure 9: Cantilever beam: (left) normalised effective piezoelectric constant e (cid:48) as a function of thenormalised beam thickness a (cid:48) , and (right) electric field modulus | E | with a (cid:48) = 1 . for piezoelec-tric (a), pure flexoelectric (b) and flexo-piezoelectric (c) beams.Figure 9 left shows the normalised effective piezoelectric constant, e (cid:48) , versus the normalisedbeam thickness, a (cid:48) , defined as a (cid:48) = − ae T µ − T , e (cid:48) := k eff k eff | µ = , k eff := (cid:115) (cid:82) Ω E · κ · E d Ω (cid:82) Ω ε · C · ε d Ω , where k eff | µ = is the effective piezoelectric constant in the absence of flexoelectric effects, i.e. with µ = .The results are in perfect agreement with the B-spline results in [8], and with the analyticalapproximation in [22]: e (cid:48) | flexo ( a (cid:48) ) (cid:39) (cid:115) a (cid:48) , e (cid:48) | flexo-piezo ( a (cid:48) ) (cid:39) (cid:115) a (cid:48) . The plots in figure 9 also illustrate how flexoelectricity is a size dependent phenomenon, withrelevant and even crucial effect for very small scales.
For further validation of the C0-IPM computational model, we now consider the open and closedcircuit example in [4], where maximum-entropy approximations (LME) were used. The problem18s solved on the same beam with the same FE mesh. The material parameters are now E = 100 GPa , ν = 0 . ,κ = 11nJ V − m − , κ = 12 .
48 nJ V − m − ,e T = − . − m − , µ T = µ L = 1 µ J V − m − ,l = µ S = e S = e L = 0 , and the the piezoelectric principal direction is again x . The mechanical boundary conditions arethe same as in the previous case.For the electrostatic boundary conditions, two different cases are considered: open and closedcircuit. The open circuit is the one considered in the previous example, with grounded right end, L x aF V Figure 10: Cantilever beam under closed-circuit boundary conditions. (a)(b) a' e ' Figure 11: Normalised effective piezoelectric constant as a function of the normalised beam thick-ness, and example of the distribution of electric potential φ in a flexo-piezoelectric beam withclosed (a) and open (b) circuit. Note that the aspect ratio of the beams has been modified to betterobserve the potential distribution along the beam.19hat is φ = 0 at x = L , as shown in figure 8. In the closed circuit, the upper side is grounded andan electrode is placed on the bottom side, that is φ = (cid:26) for x = a/ V for x = − a/ , where V is a free constant value, see figure 10. The electrode condition is enforced setting all po-tential nodal values on the bottom boundary to be equal to the first one, with Lagrange multipliersin our implementation.Figure 11 shows the normalised effective piezoelectric constant e (cid:48) as a function of the nor-malised thickness a (cid:48) . Again, we observe that flexoelectricy becomes relevant for small scales. Forthe open circuit, comparing to the previous results in figure 9, where ν = µ L = 0 , this more gen-eral model gives lower values for the normalised effective piezoelectric constant. On other hand,the open circuit setting leads to larger values of the effective piezoelectric constant. The numericalresults are in perfect agreement with the LME results in [4] demonstrating again the applicabilityof C0-IPM for the study and design of flexoelectric devices. In this section we consider an actuator beam also from [4]. The displacement is fixed on the leftboundary, and a potential difference is applied at the top and bottom sides, leading to a bending ofthe beam. That is, u (0 , x ) = , φ (cid:16) x , a (cid:17) = 0 , φ (cid:16) x , − a (cid:17) = V, on the same beam, i.e. Ω = (0 , a ) × ( − a/ , a/ . The material parameters are (17) and theapplied voltage is V = − a MV .Figure 12 shows the potential on the deformed beam for width a = 2 . µ m . The potentialseems to be smooth, but the section along x = 0 in figure 13 reveals a sharp variation close to theright end. Consequently, the electric field also presents sharp variations close to the right end, asshown in figure 14.Figure 12: Actuator beam with width a = 2 . µ m and l = 0 : deformed beam and potential.20igure 13: Actuator beam with width a = 2 . µ m and l = 0 : normalised potential ( φ/ alongthe x = 0 horizontal mid section. The potential presents a sharp variation close to the rightboundary. Magenta dots correspond to the boundary of the p = 4 quadrilaterals. -0.5 0 0.5 x /a -2024681012 no r m a li z ed E Electric field at the right end (x =20a) E1E2
Figure 14: Actuator beam with l = 0 : detail of the normalised vertical electric field ( E / )close to the right end, and plot of the components of the normalised electric field at the right end( x = 20 a ).Figure 15: Adapted p=4 quadrilateral FE mesh for l = 0 with min( h ) = a/ (top) and for l = 0 . a with min( h ) = a/ (bottom). 21 x /a -5051015 no r m a li z ed E Electric field at the right end (x =20a) E1E2
Figure 16: Actuator beam with l = 0 . a : detail of the normalised vertical electric field E close tothe right end, and plot of the components of the normalised electric field at the right end ( x = 20 a ).These results have been computed on the adapted quadrilateral mesh in figure 15 (top), withdegree p = 4 and β = 1 . The mesh has been refined to capture the sharp variations in thesolution; otherwise, numerical oscillations spoil the solution in the whole domain. It is also worthmentioning that the plot in figure 14 coincides in magnitude and shape with the results in [4] withLME, but getting rid of the smooth oscillations.Sharp variations along the boundary in the solution of flexoelectricity problems can be evenmore pronounced, as can be observed in figure 16. In this case the problem is solved with straingradient elasticity, with l = 0 . a , on the p = 4 adapted mesh in figure 15 with min( h ) = a/ ;with smaller element size along the boundary to capture the high curvatures in the electric field.The stabilization parameter is again taken as β = 100 El / min( h ) , providing stable results. The implementation of periodicity boundary conditions in the C0-IPM method is straightforward,by simply considering the periodicity faces as interior faces and imposing the periodicity con-straints on the boundary nodal values. Considering the periodicity faces as interior faces, that is in I , ensures that C continuity is enforced in weak form and that internal forces are equilibrated alsoon the periodicity boundary. The periodicity conditions for the nodal values can be implemented,for instance, by means of Lagrange multipliers, or reducing the system to the periodic space.As a verification example, figure 17 shows the evolution of the error under nested refinementfor the solution of the flexoelectricity coupled problem (1) in a square domain Ω = (0 , witha regular triangular mesh. First Dirichlet and second Neumann conditions, (1c) and (1d), are seton the top and bottom boundaries, and periodicity is imposed in the x direction. That is, (7) is22 log (h) -6-4-202 l og ( L2 e rr o r u ) u log (h) -4-202 l og ( L2 e rr o r ) -0.8-0.5-0.1 0.31.21.6 2.43.13.3 3.53.94.1 Figure 17: Convergence test for the solution of the flexoelectricity equations in a square domain,with periodicity in the x direction.imposed identifying the left and right boundary as the same boundary and including it in I . Thebody force b , the free charge q , and the data for the boundary conditions on the top and bottomboundaries, are set so that the analytical solution is (14). The stabilization parameter is (13) with α = 100 .The errors exhibit the same behavior as in the convergence analysis in section 4.1. The flexoelectricity equations are now solved in a cube,
Ω = (0 , . , to show the applicability ofthe method also in 3D. The mesh for degree p = 2 and the second level of refinement is shown infigure 18. First Dirichlet and second Neumann boundary conditions are considered in the wholeboundary, and the material parameters are (15). The data is set so that the solution is u = [cos(2 π ( x + 2 x − x )) , sin(2 π ( x + 2 x − x )) , cos(2 π ( x + 2 x − x ))] T ,φ = sin(2 π ( x + 2 x − x )) . Figures 19 and 20 show the convergence plots with β = 100 El /h , for strain gradient elasticity(solving the decoupled problem) and for flexoelectricity, respectively. As in 2D, the method doesnot converge for degree p = 1 ; thus, we show the results for p = 2 , , .Robust high-order convergence is observed in all cases, providing accurate results. Again, inagreement with the analysis in [13], the convergence is suboptimal; but still with order close to p + 1 for the displacement in the strain gradient elasticity problem for p ≥ . Again, we alsoobserve that the flexoelectricy coupling provokes a loss in the convergence rate and the accuracyof the solution; with order close to p for the displacement and the potential in this example.23igure 18: 3D mesh for degree p = 2 and second level of refinement. -1.5 -1 -0.5 0 log (h) -4-3-2-10 l og ( L2 e rr o r u ) u log (h) -5-4-3-2-10 l og ( L2 e rr o r ) Figure 19: 3D solution in a cube: convergence plots for strain gradient elasticity and electricpotential (decoupled problem). -1.5 -1 -0.5 0 log (h) -4-3-2-10 l og ( L2 e rr o r u ) u log (h) -3-2-101 l og ( L2 e rr o r ) Figure 20: 3D solution in a cube: convergence plots for the coupled flexoelectricity problem.24urther numerical experiments show that, with these regular hexahedra meshes, the error hasvery little dependency on the particular value of β ≥ El /h . A novel C0-IPM formulation for strain gradient elasticity and flexoelectricity is proposed. Theweak form involves second derivatives of the displacement in the interior of the elements, plusintegrals on the element faces, weakly imposing continuity of the displacement derivatives, as wellas equilibrium of internal forces across element faces and on interior edges (vertexes in 2D).The formulation is stable, with a symmetric and positive definite matrix for the strain gradientelasticity operator, for large enough interior parameter β . An eigenvalue problem is stated todetermine a bound for β , which leads to a general formula for the parameter: β = αEl /h , withconstant α independent of the element size. Thus, differently to non-consistent penalty methods,and as usual in interior penalty methods, moderate values for β provide stable and accurate results.Standard C FE approximations are considered, retaining the advantages and computationalefficiency of high-order FE. The implementation is based on assembly of elemental matrices, withstandard FE numerical integration and nodal approximation, the discretization can be adapted tothe geometry and locally refined where needed, no additional unknowns are needed, and materialinterfaces can be directly considered just adapting the mesh, as usual in FE computations.The application of C0-IPM to problems with periodicity boundary conditions is straightfor-ward, just considering the periodicity faces as interior faces (thus, imposing C continuity andequilibrium of forces in weak form) and setting the periodicity conditions on the nodal values.Convergence tests, on 2D non-uniform curved triangular meshes and on 3D hexahedra regularmeshes, show high-order convergence of the method for degree p ≥ . A slow continuous loss inthe convergence rate for increasing β is observed for p = 3 , which is in agreement with the analysisfor the biharmonic equation in [13] and the results for Kirchhoff plates in [14]. Fortunately, for p = 4 the convergence shows little dependency on β . In any case, in all examples, the convergencerates are at least close to p for both variables, demonstrating the good behaviour of the method for p ≥ .The computational tool is also validated by comparison with previous works solving realisticactuator and sensor problems on a beam, with perfect agreement. Acknowledgements
This work was supported by the European Research Council (StG-679451 to Irene Arias) andGeneralitat de Catalunya (2017-SGR-1278) 25 eferences [1] Longlong Shu, Renhong Liang, Zhenggang Rao, Linfeng Fei, Shanming Ke, and Yu Wang.Flexoelectric materials and their related applications: A focused review.
Journal of AdvancedCeramics , 8(2):153–173, 2019.[2] Sheng Mao, Prashant Purohit, and N. Aravas. Mixed finite-element formulations in piezo-electricity and flexoelectricity.
Proceedings of the Royal Society A: Mathematical, Physicaland Engineering Science , 472:20150879, 2016.[3] Feng Deng, Qian Deng, Wenshan Yu, and Shengping Shen. Mixed Finite Elements for Flex-oelectric Solids.
Journal of Applied Mechanics , 84(8), 2017.[4] Amir Abdollahi, Christian Peco, Daniel Mill´an, Marino Arroyo, and Irene Arias. Computa-tional evaluation of the flexoelectric effect in dielectric solids.
Journal of Applied Physics ,116:093502–093502, 09 2014.[5] Hamid Ghasemi, Harold S. Park, and Timon Rabczuk. A level-set based IGA formulation fortopology optimization of flexoelectric materials.
Computer Methods in Applied Mechanicsand Engineering , 313:239 – 258, 2017.[6] S.S. Nanthakumar, Xiaoying Zhuang, Harold S. Park, and Timon Rabczuk. Topology opti-mization of flexoelectric structures.
Journal of the Mechanics and Physics of Solids , 105:217– 234, 2017.[7] R. Sevilla and S. Fern´andez-M´endez. Numerical integration over 2D NURBS-shaped do-mains with applications to NURBS-enhanced { FEM } . Finite Elements in Analysis and De-sign , 47(10):1209 – 1220, 2011.[8] D. Codony, O. Marco, S. Fern´andez-M´endez, and I. Arias. An immersed boundary hierar-chical B-spline method for flexoelectricity.
Computer Methods in Applied Mechanics andEngineering , 354:750 – 782, 2019.[9] F. de Prenter, C.V. Verhoosel, G.J. van Zwieten, and E.H. van Brummelen. Condition num-ber analysis and preconditioning of the finite cell method.
Computer Methods in AppliedMechanics and Engineering , 316:297 – 327, 2017. Special Issue on Isogeometric Analysis:Progress and Challenges.[10] Douglas N. Arnold. An interior penalty finite element method with discontinuous elements.
SIAM Journal on Numerical Analysis , 19(4):742–760, 1982.[11] J. Nitsche. ber ein variationsprinzip zur lsung von dirichlet-problemen bei verwendung vonteilrumen, die keinen randbedingungen unterworfen sind.
Abhandlungen aus dem Mathema-tischen Seminar der Universitt Hamburg , 36(1):9–15, 1971.2612] Gerald Engel, Krishna Garikipati, Thomas Hughes, Mats Larson, Luca Mazzei, and R. Tay-lor. Continuous/discontinuous finite element approximations of fourth-order elliptic prob-lems in structural and continuum mechanics with applications to thin beams and plates.
Com-puter Methods in Applied Mechanics and Engineering , 191:3669–3750, 2002.[13] Susanne C. Brenner and Li-Yeng Sung. C0 Interior Penalty Methods for fourth order ellipticboundary value problems on polygonal domains.
Journal of Scientific Computing , 22(1):83–118, 2005.[14] Dani Fojo, David Codony, and Sonia Fern´andez-M´endez. A C0 Interior Penalty Method for4th order PDEs.
Reports@SCM , 2020.[15] Garth Wells, Krishna Garikipati, and Luisa Molari. A discontinuous galerkin method forstrain gradient-dependent damage.
Computer Methods in Applied Mechanics and Engineer-ing , 193:3633–3645, 2003.[16] Susanne C. Brenner, Shiyuan Gu, Thirupathi Gudi, and Li-Yeng Sung. A quadratic C0 in-terior penalty method for linear fourth order boundary value problems with boundary con-ditions of the Cahn-Hilliard type.
SIAM Journal on Numerical Analysis , 50(4):2088–2110,2012.[17] Sonia Fern´andez-M´endez and Antonio Huerta. Imposing essential boundary conditionsin mesh-free methods.
Computer Methods in Applied Mechanics and Engineering ,193(12):1257 – 1275, 2004.[18] Michael Griebel and Marc Schweitzer. A particle-partition of unity method - part v: Bound-ary conditions.
Geometric Analysis and Nonlinear Partial Differential Equations , 41:519–542, 05 2002.[19] Chandrasekhar Annavarapu, Martin Hautefeuille, and John E. Dolbow. A robust Nitsche’sformulation for interface problems.
Computer Methods in Applied Mechanics and Engineer-ing , 225-228:44 – 54, 2012.[20] Qi Chen and Ivo Babu˘ska. Approximate optimal points for polynomial interpolation of realfunctions in an interval and in a triangle.
Computer Methods in Applied Mechanics andEngineering , 128(3):405 – 417, 1995.[21] Eloi Ruiz-Giron´es, Abel Gargallo-Peir´o, Josep Sarrate, and Xevi Roca. Automatically impos-ing incremental boundary displacements for valid mesh morphing and curving.
Computer-Aided Design , 112:47 – 62, 2019.[22] M. Majdoub, Pradeep Sharma, and Tcagin Cagin. Enhanced size-dependent piezoelectricityand elasticity in nanostructures due to the flexoelectric effect.