A Stable Mixed FE Method for Nearly Incompressible Linear Elastostatics
RReceived: Added at production Revised: Added at production Accepted: Added at productionDOI: xxx/xxxx
RESEARCH ARTICLE
A Stable Mixed FE Method for Nearly Incompressible LinearElastostatics
Eirik Valseth* | Albert Romkes | Austin R. Kaul | Clint Dawson Oden Institute for ComputationalEngineering and Sciences, The Universityof Texas at Austin, Texas, USA Department of Mechanical Engineering,South Dakota School of Mines &Technology, South Dakota, USA
Correspondence *Eirik Valseth, Oden Institute forComputational Engineering and Sciences,The University of Texas at Austin, Austin,TX 78712, USA. Email: [email protected]
Summary
We present a new, stable, mixed ๏ฌnite element (FE) method for linear elastostat-ics of nearly incompressible solids. The method is the automatic variationally stableFE (AVS-FE) method of Calo, Romkes and Valseth, in which we consider a Petrov-Galerkin weak formulation where the stress and displacement variables are in thespace ๐ป ( div , ฮฉ) and ๐ป (ฮฉ) , respectively. This allows us to employ a fully conform-ing FE discretization for any elastic solid using classical FE subspaces of ๐ป ( div , ฮฉ) and ๐ป (ฮฉ) . Hence, the resulting FE approximation yields both continuous stressesand displacements.To ensure stability of the method, we employ the philosophy of the discontinuousPetrov-Galerkin (DPG) method of Demkowicz and Gopalakrishnan and use optimaltest spaces. Thus, the resulting FE discretization is stable even as the Poissonโs ratio ๐ โ . , and the system of linear algebraic equations is symmetric and positivede๏ฌnite. Our method also comes with a built-in a posteriori error estimator as wellas well as indicators which are used to drive mesh adaptive re๏ฌnements. We presentseveral numerical veri๏ฌcations of our method including comparisons to existing FEtechnologies. KEYWORDS: discontinuous Petrov-Galerkin method, a priori error estimation, adaptive mesh re๏ฌnement, nearlyincompressible elasticity, composite materials
Linear elastostatics is arguably the most successful area of application of the classical (Galerkin) FE method. For homoge-neous isotropic engineering materials, such as steel and aluminum, the Galerkin method is stable, satis๏ฌes a best approximation a r X i v : . [ m a t h . NA ] J a n EIRIK VALSETH
ET AL property in terms of elastic strain energy and is computationally e๏ฌcient. However, for commonly employed modern engineer-ing materials such as rubbers and soft plastics, i.e., nearly incompressible materials, the Galerkin method su๏ฌers from locking and loss of discrete stability (see, e.g., ). In , Phillips and Wheeler investigate this phenomenon in great detail and highlightthat error estimates for the classical FE method depend on the factor ๐ ) , which clearly tends to in๏ฌnity as ๐ โ . .Mixed FE methods provide functional settings in which certain FE discretizations are stable for mixed forms of the elasto-statics problem as well as for the Stokes equations which and can be shown to be equivalent to the equations of linear elastostaticswhen ๐ = 0 . . Other mixed FE methods based on the consideration of the underlying equations of elastostatics using the com-pliance tensor do not su๏ฌer from this loss of stability, but often require additional constraints to ensure symmetric stresses ,called Hellinger-Reissner formulations, and are typically more computationally costly than the primal Galerkin formulation.In , Arnold and Winther present an element for the Hellinger-Reissner formulation which was one of the ๏ฌrst stable elementsutilizing polynomial bases for both stress and displacement. The di๏ฌculty in building conforming approximation spaces for thestress has also resulted in several nonconforming mixed methods, see e.g., where the FE approximation of the stress does notreside in the space dictated by the weak formulation. The task of establishing approximation spaces for these mixed FE methodsis certainly not trivial and is an active area of investigation and recent publications include .Stabilized FE methods that adjust the functionals of the weak formulation can be used to ensure discrete stability . This typeof stabilization is performed for both the mixed and classical FE methods, see the work of Nakshatrala et al. as well as Masud et al. . However, stabilized methods generally require arduous analyses to establish a proper choice of penalization/stabilizationparameters. Reduced integration methods are also commonly used when approaching the incompressible limit . The discontin-uous Galerkin (DG) method also remains a popular choice for nearly incompressible elastostatics . Generally, these achievestability by adjusting the inter-element jump or average terms by weights in a manner similar to the stabilized FE methods.Stable FE methods such as the least squares FE methods (LSFEMs) (see, e.g., text by Bochev and Gunzberger ) or thediscontinuous Petrov-Galerkin (DPG) Method of Demkowicz and Gopalakrishnan can be employed to resolve the stabilityissue. The LSFEM has been applied to linear elastostatics in and in , a weighted ๏ฌrst-order system least squares is appliedsuccessfully to nearly incompressible materials. Gopalakrishnan and Qiu provide an analysis of the well-posedness of the DPGmethod applied to linear elastostatics in . In , Bramwell et al. consider two distinct DPG methods for the nearly incompressibleelastostatics problem that are locking free and present numerical veri๏ฌcations highlighting capabilities as ๐ โ . . The DPGhas also been successfully applied to this problem in several works, including the fully incompressible case in employing thecompliance tensor to avoid locking in that case. In , the DPG method is applied to the problem of linear elastostatics andseveral variational formulations are considered including for the case of nearly incompressible materials. In particular, in , theidea of coupling multiple weak formulations throughout the computational domain is explored in great detail. IRIK VALSETH
ET AL In the classical FE method, the approximations of displacements of the equivalent weak form of the underlying partial dif-ferential equation (PDE) of static equilibrium are sought in ๐ถ continuous polynomial spaces and stress approximations areestablished by computing gradients of the displacements, i.e., the stresses are piecewise discontinuous. On the other hand, mixedFE methods for the linear elastostatics problem consider an equivalent ๏ฌrst-order system of the underlying PDE. This ๏ฌrst-order system description can lead to weak forms which allow stresses to be in ๐ป ( div , ฮฉ) and displacements that are in ๐ฟ (ฮฉ) (see Section 2.4 of for a thorough discussion on other options). Hence, in the FE approximations the displacements must besought in piecewise discontinuous polynomial function spaces. The theory of distributions ensures that optimally convergentFE solutions can be established for both classical and mixed FE methods, as well as for their properly stabilized counterparts if ๐ is close to . . However, the resulting numerical approximations are not physical, as we know that both the displacement andcertain components of the stress ๏ฌelds are continuous. We know of two remedies to establish both continuous displacementsand stresses, ๐ ) the isogeometric FE methods of Hughes et al. which uses higher order bases for the discrete FE approxima-tion, i.e., ๐ถ ๐ continuity and ๐๐ ) the ๐ โ version FE method of Surana et al. which employs higher order bases as well as a leastsquares approach. The popularity of the isogeometric FE method has grown signi๏ฌcantly over the last decade, but the stabilityissue of nearly incompressible materials still persist. In , Taylor introduces a mixed version of the isogeometric FE methodsfor incompressible solids where discontinuous stress approximations are sought.The automatic variationally stable ๏ฌnite element (AVS-FE) method introduced by Calo, Romkes and Valseth in providesa framework, much like the DPG of Demkowicz and Gopalakrishnan , to establish stable FE approximations for any PDE.However, the AVS-FE di๏ฌers in its choice of trial spaces while employing the DPG concept of optimal discontinuous testfunctions. In addition to the approximation of the trial variables, the AVS-FE also comes with a "built-in" error estimator anderror indicators that can be employed to drive mesh adaptive strategies. The stability property of the AVS-FE allows us toderive Petrov-Galerkin weak formulations that are posed with trial functions that are in classical Hilbert spaces, e.g., ๐ป (ฮฉ) and ๐ป ( div , ฮฉ) . Hence, the corresponding FE approximations are to be sought in classical continuous FE approximation spacesyielding continuous
FE approximations for all trial variables. The LSFEMs presented in also pose weak formulations inHilbert spaces as the AVS-FE but considers alternative formulations for the elasticity problem and considers nonconformingapproximations for the displacement.In this paper, we build upon the preliminary investigation of Valseth in for the AVS-FE method applied to linear elastostaticsof nearly incompressible media. We introduce our model problem and notations in Section 2.1, The weak formulation and itscorresponding FE discretzation are presented in Section 2 in conjunction with a brief review of the AVS-FE methodology. InSection 2.4, we present optimal a priori error estimates. Several numerical veri๏ฌcations are presented in Section 3 highlightingthe stability of our method as ๐ โ . , including an asymptotic convergence study with comparisons to existing FE methods.We draw conclusions and discuss future works in Section 4. EIRIK VALSETH
ET AL
The AVS-FE method provides a functional setting to analyze linear boundary value problems (BVPs) in which the underlyingdi๏ฌerential operator is non self-adjoint or leads to unstable FE discretizations. In this section we introduce our model problem,and brie๏ฌy review the AVS-FE method. A thorough introduction can be found in . Let ฮฉ โ โ be a bounded open domain, containing a linearly elastic, nearly incompressible, and possible heterogeneous solid.The boundary ๐ ฮฉ is partitioned into two open and disjoint segments ฮ ๐ก and ฮ ๐ฎ , such that ๐ ฮฉ = ฮ ๐ก โช ฮ ๐ฎ . As depicted in Figure 1,the body is in static equilibrium under the action of external body loads ๐ โ [ ๐ฟ (ฮฉ)] in ฮฉ , surface tractions ๐ญ โ ๐ป โ (ฮ ๐ก ) on ฮ ๐ก ,as well prescribed displacements ๐ฎ โ ๐ป (ฮ ๐ฎ ) on ฮ ๐ฎ . Since the solid is assumed to be linearly elastic, its constitutive behavioris governed by Generalized Hookeโs Law, i.e.: ๐ = ๐ ๐บ , (1)where ๐ denotes the (2D) Cauchy stress tensor, ๐บ the (2D) Green strain tensor, and ๐ the fourth order (Riemann) elasticitytensor, with elliptic and symmetric Riemann coe๏ฌcients ๐ธ ๐๐๐๐ โ ๐ฟ โ (ฮฉ) . In this work, we limit our focus to problems in whichthe deformations in the material remain small and therefore the kinematic relation between the strain tensor ๐บ and displacement๏ฌeld ๐ฎ is linear and governed by: ๐บ = 12 [ ๐๐ฎ + ( ๐๐ฎ ) ๐ ] . (2)With these notations and relations in force, the equilibrium state of the solid is represented by the following BVP, governing thedisplacement ๏ฌeld ๐ฎ : Find ๐ฎ โ ๐ป (ฮฉ) such that: โ ๐ โ ( ๐๐๐ฎ ) = ๐ , in ฮฉ , ๐๐๐ฎ ๐ง = ๐ญ , on ฮ ๐ก , ๐ฎ = , on ฮ ๐ฎ , (3)where ๐ง denotes the outward unit normal vector to ๐ ฮฉ . In this paper, we consider the speci๏ฌc scenario in which the solid iscomprised of one or more constituents with nearly incompressible material properties. Hence, the Riemann coe๏ฌcients ๐ธ ๐๐๐๐ can involve values of the Poisson Ratio ๐ that are very close to, but still less than, . .In the following, we shall use the following notations:โข inner products between vector valued functions are denoted with the single dot symbol โ , and inner products betweentensor valued functions are denoted by the colon or double dot symbol โถ . IRIK VALSETH
ET AL ๐ญ โชข ฮฉฮ ๐ฎ ๐ง๐ FIGURE 1
The model problem.โข โ ๐ is the diameter of element ๐พ ๐ .โข in weak formulations, we present edge integrals using trace the operators: ๐ ) ๐พ ๐ โถ ๐ป ( ๐พ ๐ ) โถ โโ ๐ป ( ๐๐พ ๐ ) as the localzeroth order trace operator and ๐๐ ) ๐พ ๐ ๐ง ๐ โถ ๐ป ( div , ๐พ ๐ ) โโ ๐ป โ1โ2 ( ๐๐พ ๐ ) denote the local normal trace operators where ๐ง ๐ isthe outward unit normal vector to the element boundary ๐๐พ ๐ (e.g., see ). AVS-FE weak formulations are established using techniques that are similar to DG and DPG methods by considering element-wise weak formulations that are subsequently summed throughout the FE mesh partition to yield global weak formulations. Wemention only key points here and omit the full derivation here for brevity but refer to for detailed derivations.To establish AVS-FE weak formulations, we ๏ฌrst require a partition ๎ผ โ of ฮฉ into convex elements ๐พ ๐ , such that: ฮฉ = int ( โ ๐พ ๐ โ ๎ผ โ ๐พ ๐ ) , ๐พ ๐ โฉ ๐พ ๐ , ๐ โ ๐. The partition ๎ผ โ is such that any discontinuities in ๐ธ ๐๐๐๐ are restricted to the boundaries of each element ๐๐พ ๐ . The BVP (3) isrecast as a ๏ฌrst-order system by using the stress tensor from the constitutive law (1), i.e.,Find ( ๐ฎ , ๐ ) โ ๐ป (ฮฉ) ร ๐ป ( div , ฮฉ) such that: ๐ โ ๐ ๐บ = in ฮฉ , โ ๐ โ ๐ = ๐ , in ฮฉ , ๐ ๐ง = ๐ญ , on ฮ ๐ก , ๐ฎ = , on ฮ ๐ฎ , (4) EIRIK VALSETH
ET AL where ๐ป ( div , ฮฉ) is the Hilbert space of tensor-valued functions which divergence is weakly continuous and ๐บ = ๐บ ( ๐ฎ ) denotesthe gradient operator in (2). Note that this ๏ฌrst-order system, or mixed form, BVP is standard for mixed FE methods for lineareleastostatics.Next, the ๏ฌrst-order system is multiplied by test functions ( ๐ฏ , ๐ฐ ) โ ๐ฟ ( ๐พ ๐ ) and enforced weakly on each individual element ๐พ ๐ โ ๎ผ โ . We then apply integration by parts locally on each element ๐พ ๐ to the term involving the divergence of the stress ๏ฌeld ๐ โ ๐ to enable weak applications of the Neumann boundary condition (BC). A subsequent summation of all elements in ๎ผ โ ,weak enforcement of Neumann and strong enforcement of Dirichlet BCs leads us to the AVS-FE weak formulation:Find ( ๐ฎ , ๐ ) โ ๐ (ฮฉ) such that: ๐ต (( ๐ฎ , ๐ ); ( ๐ฏ , ๐ฐ )) = ๐น ( ๐ฏ ) , โ( ๐ฏ , ๐ฐ ) โ ๐ ( ๎ผ โ ) , (5)where the test space ๐ ( ๎ผ โ ) is broken, the bilinear form, ๐ต โถ ๐ (ฮฉ) ร ๐ ( ๎ผ โ ) โโ โ , and linear functional, ๐น โถ ๐ ( ๎ผ โ ) โโ โ arede๏ฌned: ๐ต (( ๐ฎ , ๐ ); ( ๐ฏ , ๐ฐ )) def = โ ๐พ ๐ โ ๎ผ โ { โซ ๐พ ๐ [ ( ๐ โ ๐ ๐บ ( ๐ฎ )) โ ๐ฐ ๐ + ๐ โถ ๐๐ฏ ๐ ] d ๐ฑ โ โฎ ๐๐พ ๐ โงต ฮ ๐ก โชฮ ๐ฎ ๐พ ๐ ๐ง ( ๐ ) ๐พ ๐ ( ๐ฏ ๐ ) d ๐ } ,๐น ( ๐ฏ ) def = โ ๐พ ๐ โ ๎ผ โ โงโชโจโชโฉ โซ ๐พ ๐ ๐ โ ๐ฏ ๐ d ๐ฑ + โฎ ๐๐พ ๐ โฉฮ ๐ก ๐ญ ๐พ ๐ ( ๐ฏ ๐ ) d ๐ โซโชโฌโชโญ , (6)and the function spaces are de๏ฌned: ๐ (ฮฉ) def = { ( ๐ฎ , ๐ ) โ [ ๐ป (ฮฉ)] ร ๐ป ( div , ฮฉ) โถ ๐พ ๐ ( ๐ฎ ) | ๐๐พ ๐ โฉฮ ๐ฎ = , โ ๐พ ๐ โ ๎ผ โ } ,๐ ( ๎ผ โ ) def = { ( ๐ฏ , ๐ฐ ) โ [ ๐ป ( ๎ผ โ )] ร [ ๐ฟ (ฮฉ)] โถ ๐พ ๐ ( ๐ฏ ๐ ) | ๐๐พ ๐ โฉฮ ๐ฎ = , โ ๐พ ๐ โ ๎ผ โ } , (7)with norms โ โ โ ๐ (ฮฉ) โถ ๐ (ฮฉ) โโ [0 , โ) and โ โ โ ๐ ( ๎ผ โ ) โถ ๐ ( ๎ผ โ ) โโ [0 , โ) de๏ฌned as: โ ( ๐ฎ , ๐ ) โ ๐ (ฮฉ) def = โโโโ โซ ฮฉ [ ๐๐ฎ โถ ๐๐ฎ + ๐ฎ โ ๐ฎ + ( ๐ โ ๐ ) + ๐ โ ๐ ] d ๐ฑ , โ ( ๐ฏ , ๐ฐ ) โ ๐ ( ๎ผ โ ) def = โโโโโ โ ๐พ ๐ โ ๎ผ โ โซ ๐พ ๐ [ โ ๐ ๐๐ฏ ๐ โถ ๐๐ฏ ๐ + ๐ฏ ๐ โ ๐ฏ ๐ + ๐ฐ ๐ โ ๐ฐ ๐ ] d ๐ฑ . (8)The norm โ โ โ ๐ ( ๎ผ โ ) is equivalent to the more classical broken norm: โ ( ๐ฃ, ๐ฐ ) โ ๐ def = โโโโโ โ ๐พ ๐ โ ๎ผ โ โซ ๐พ ๐ [ ๐ ๐ฃ ๐ โ ๐๐ฏ ๐ + ๐ฃ ๐ โ ๐ฃ ๐ + ๐ฐ ๐ โ ๐ฐ ๐ ] d ๐ฑ . (9)Note that the edge integrals in (6) are to be interpreted as duality pairings in ๐ป ( ๐๐พ ๐ ) ร ๐ป โ1โ2 ( ๐๐พ ๐ ) , but we employ notationthat is engineering convention here using the integral representation. Most importantly, since ( ๐ฎ , ๐ ) โ ๐ป (ฮฉ) ร ๐ป ( div , ฮฉ) , IRIK VALSETH
ET AL these integrals are well de๏ฌned and our trial space is (weakly) continuous. As the trial and test spaces are of di๏ฌerent regularitywe are in a Petrov-Galerkin setting, particularly a DPG setting, since the test space is broken. However, since the trial space is ๐ป (ฮฉ) ร ๐ป ( div , ฮฉ) our functional setting di๏ฌers from that of DPG methods in which the regularity of the trial space is reducedby introducing variables on the edge of each element.In the spirit of the DPG method, we now introduce an equivalent norm on the trial space, the energy norm โ โ โ B โถ ๐ (ฮฉ) โโ [0 , โ) : โ ( ๐ฎ , ๐ ) โ B def = sup ( ๐ฏ , ๐ฐ )โ ๐ ( ๎ผ โ ) โงต {( , )} | ๐ต ( ๐ฎ , ๐ ); ( ๐ฏ , ๐ฐ )) |โ ( ๐ฏ , ๐ฐ ) โ ๐ ( ๎ผ โ ) , (10)and a Riesz representation problem for ( ๐ฉ , ๐ซ ) , the optimal test functions: ( ( ๐ฉ , ๐ซ ); ( ๐ฏ , ๐ฐ ) ) ๐ ( ๎ผ โ ) = ๐ต ( ( ๐ฎ , ๐ ); ( ๐ฏ , ๐ฐ ) ) , โ( ๐ฃ, ๐ฐ ) โ ๐ ( ๎ผ โ ) . (11)The Riesz representation problem is well posed with unique solutions due to the inner product in the left hand side (LHS) andguarantees the stability of DPG methods. The Riesz representation problem also leads to the following norm equivalence: โ ( ๐ฎ , ๐ ) โ B = โ ( ๐ฉ , ๐ซ ) โ ๐ ( ๎ผ โ ) , (12)which will be employed extensively in the following. For details on optimal test functions and proof of the norm equivalence,we refer to . Due to the energy norm, the bilinear form (6) has continuity and inf-sup constants equal to one and the loadfunctional is also continuous which can be shown using classical techniques. Hence, we have established a well posed weakformulation of the linear elastostatics BVP using continuous trial spaces for both displacement and stress ๏ฌelds, i.e., ๐ป (ฮฉ) ร ๐ป ( div , ฮฉ) , in terms of the energy norm (10).Well-posedness in terms of the energy norm is essentially an assumption of DPG methods as we de๏ฌne a norm that ensure inf-sup and continuity conditions of the bilinear form. For completeness, we also provide the following lemma of well-posednessin standard Sobolev norms by ๏ฌrst stating two important results: Proposition ( ๐ฎ D , ๐ D ) โ [ ๐ป ๐ฎ (ฮฉ)] ร [ ๐ฟ (ฮฉ)] be the solution of the dual mixed formulation:Find ( ๐ฎ D , ๐ D ) โ [ ๐ป ๐ฎ (ฮฉ)] ร [ ๐ฟ (ฮฉ)] such that: โซ ฮฉ [( ๐ D โ ๐ ๐บ ( ๐ฎ D ) ) โ ๐ฐ + ๐ D โถ ๐๐ฏ ] d ๐ฑ โโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโโ ๐ D ( ๐ฎ D , ๐ D ) , ( ๐ฏ , ๐ฐ )) = โซ ฮฉ ๐ โ ๐ฏ d ๐ฑ , โ( ๐ฏ , ๐ฐ ) โ [ ๐ป ๐ฎ (ฮฉ)] ร [ ๐ฟ (ฮฉ)] , (13)which is well posed. Hence, the bilinear form satis๏ฌes the inf-sup condition: โ ๐พ > ( ๐ฏ , ๐ฐ )โ ๐ (ฮฉ) | ๐ D ( ๐ฎ D , ๐ D ) , ( ๐ฏ , ๐ฐ )) |โ ( ๐ฏ , ๐ฐ ) โ ๐ (ฮฉ) โฅ ๐พ โ ( ๐ฎ D , ๐ D ) โ ๐ ๐ท (ฮฉ) , (14) EIRIK VALSETH
ET AL where ๐ ๐ท (ฮฉ) = ๐ (ฮฉ) = [ ๐ป ๐ฎ (ฮฉ)] ร [ ๐ฟ (ฮฉ)] and ๐ป ๐ฎ (ฮฉ) is the space of ๐ป functions that satisfy homogeneous Dirichletconditions on ฮ ๐ฎ = ๐ ฮฉ . Proof : see Theorem 2.1 in .We write the bilinear form (6) as: ๐ต (( ๐ฎ , ๐ ) , ( ๐ฏ , ๐ฐ )) = ๐ D (( ๐ฎ , ๐ ) , ( ๐ฏ , ๐ฐ )) + โจ ๐พ ๐ ๐ง ( ๐ ) , ๐พ ๐ ( ๐ฏ ๐ ) โฉ ฮ โ , (15)where โจ ๐พ ๐ ๐ง ( ๐ ) , ๐พ ๐ ( ๐ฏ ๐ ) โฉ ฮ โ def = โ ๐พ ๐ โ ๎ผ โ โฎ ๐๐พ ๐ { ๐พ ๐ ๐ง ( ๐ ) ๐พ ๐ ( ๐ฏ ๐ ) } d ๐ . Proposition ๐ โ ๐ป ( div , ๎ผ โ ) and ๐ฃ ๐ โ ๐ป ( ๎ผ โ ) . Then: โ ๐พ ๐ > ๐ฃ โ ๐ป ( ๎ผ โ ) |โจ ๐พ ๐ ๐ง ( ๐ ) , ๐พ ๐ ( ๐ฏ ๐ ) โฉ ฮ โ |โ ๐ฃ โ ๐ป ( ๎ผ โ ) โฅ ๐พ ๐ โ ๐ โ ฬ๐ (ฮ โ ) , (16)where ๐ป ( div , ๎ผ โ ) denotes the broken ๐ป ( ๐๐๐ฃ ) space and โ ๐ โ ฬ๐ (ฮ โ ) is the minimum energy extension norm: โ ๐ โ ฬ๐ (ฮ โ ) def = sup ๐ฃ โ ๐ป ( ๎ผ โ ) |โจ ๐พ ๐ ๐ง ( ๐ ) , ๐พ ๐ ( ๐ฏ ๐ ) โฉ ฮ โ |โ ๐ฃ โ ๐ป ( ๎ผ โ ) = inf โ ๐ โ H( div , ฮฉ) . (17)Additionally, if ๐ โ ๐ป ( div , ฮฉ) : โจ ๐พ ๐ ๐ง ( ๐ ) , ๐พ ๐ ( ๐ฏ ๐ ) โฉ ฮ โ = 0 . (18) Proof : see Theorem 2.3 in . Lemma ( ๐ฎ , ๐ ) โ ๐ (ฮฉ) and ( ๐ฏ , ๐ฐ ) โ ๐ ( ๎ผ โ ) . Then, the AVS-FE weak formulation (5) satis๏ฌes all conditions of the Babuลกka Lax-Milgram Theorem and is well posed. Proof : The load functional and bilinear form (6) are continuous in terms of the norms โ โ โ ๐ and โ โ โ ๐ (ฮฉ) (see (9) and (8)) dueto the Cauchy-Schwarz inequality. The following inf-sup condition: โ ๐ถ > ( ๐ฏ , ๐ฐ )โ ๐ ( ๎ผ โ ) | ๐ต ( ๐ฎ , ๐ ) , ( ๐ฏ , ๐ฐ )) |โ ( ๐ฏ , ๐ฐ ) โ ๐ โฅ ๐พ โ ( ๐ฎ , ๐ ) โ ๐ (ฮฉ) , (19)is satis๏ฌed due to Theorem 3.3 in due to Propositions 2.2 and 2.1. Essentially, Theorem 3.3 in states that the stability of abroken variational formulation is inherited from its unbroken analogue. IRIK VALSETH
ET AL Remark Keith et al. consider several possible weak formulationsfor the elastostatics problem and the DPG method and perform a rigorous analysis showing their well posedness.
To establish FE discretizations of the weak formulation (5), the AVS-FE takes the approach of classical FE methods and seekscontinuous polynomial approximations that are in conforming subspaces of the ๐ป ( div , ฮฉ) ร ๐ป (ฮฉ) trial spaces. Hence, for thedisplacement ๏ฌeld we use classical ๐ถ (ฮฉ) continuous Lagrange polynomials. Generally, in mixed FE methods this choice leadsto unstable and inconsistent FE discretizations and is avoided and the stress ๏ฌeld is sought in a Raviart-Thomas (RT) or Brezzi-Douglas-Marini (BDM) space . Due to the stability properties of the AVS-FE method, it is often convenient to employ the same ๐ถ (ฮฉ) polynomials for the stress variable. As reported in , for convex domains and smooth solutions, the ๐ถ (ฮฉ) are superior.In Section 3 we present numerical veri๏ฌcations comparing the approximations from these spaces. Furthermore, we present averi๏ฌcation where we again compare the classical RT spaces with the ๐ถ (ฮฉ) polynomials for a physical application in whichthe stress ๏ฌeld is such that it is discontinuous in the tangential direction.Hence, we seek numerical approximations ( ๐ฎ โ , ๐ โ ) of ( ๐ฎ , ๐ ) of the weak formulation (5) and represent the approximations aslinear combinations of the trial bases ( ๐ ๐ ( ๐ฑ ) , ๐ ๐ ( ๐ฑ )) โ ๐ โ (ฮฉ) and their corresponding degrees of freedom: ๐ฎ โ ( ๐ฑ ) = ๐ โ ๐ =1 ๐ฎ โ๐ ๐ ๐ ( ๐ฑ ) , ๐ โ ( ๐ฑ ) = ๐ โ ๐ =1 ๐ โ,๐ ๐ ๐ ( ๐ฑ ) . (20)Now, the test space, which is discontinuous, is to be constructed by the DPG philosophy using optimal test functions de๏ฌnedby the discrete equivalent of the Riesz representation problem (11). Thus, the optimal test space is spanned by functions thatare solutions of the global weak problems, e.g., for a trial basis function ๐ ๐ ( ๐ฑ ) for the displacement variable, its correspondingoptimal test function ( ฬ๐ ๐ , ฬ๐ ๐ ) is de๏ฌned by: ( ( ๐ซ , ๐ณ ); ( ฬ๐ ๐ , ฬ๐ ๐ ) ) ๐ ( ๎ผ โ ) = ๐ต ( ( ๐ ๐ , ); ( ๐ซ , ๐ณ ) ) , โ( ๐ซ , โ ๐ ( ๎ผ โ ) , ๐ = 1 , โฆ , ๐. (21)Inspection of (21) reveals the ingenuity of the DPG philosophy since the test space is broken, we do not need to solve thisproblem globally but rather element-wise local analogues which can be solved in a completely decoupled local fashion. Hence, EIRIK VALSETH
ET AL the AVS-FE approximation is governed by:Find ( ๐ฎ โ , ๐ โ ) โ ๐ โ (ฮฉ) such that: ๐ต (( ๐ฎ โ , ๐ โ ); ( ๐ฏ โ , ๐ฐ โ )) = ๐น ( ๐ฏ โ ) , โ( ๐ฏ โ , ๐ฐ โ ) โ ๐ โ ( ๎ผ โ ) , (22)where the test space ๐ โ ( ๎ผ โ ) is spanned by the approximated optimal test functions computed from local equivalents of (21).The choice we have made of fully continuous trial spaces has several important consequences: ๐ ) as the bilinear form (6) issuch that information is transferred from element-to-element by the continuity of trial functions alone, the optimal test functionsin the AVS-FE have the same support as its trial basis functions. Hence, the resulting global sti๏ฌness matrix has the samesparseness as mixed FE methods. ๐๐ ) the local optimal test function problems can be solved by using the same polynomialdegree of approximation as the trial functions which de๏ฌne each problem. Thus, the cost incurred to establish the optimal testfunctions is kept as low as possible (see Remark 2.2). ๐๐๐ ) ๏ฌnally, the AVS-FE optimal test functions can be implemented inlegacy FE software in which continuous polynomials are the only available basis functions by rede๏ฌning the element sti๏ฌnessmatrix assembly process. Remark ๐ฟ , i.e., the LSFEM, this is not possible inpractical computations and we consider an approximation of these functions . Thus, there is a potential loss of discrete stabilityif the optimal test functions are computed without su๏ฌcient accuracy. Su๏ฌcient accuracy is ensured by the existence of (local)Fortin operators . The construction of such operators for the DPG method is studied in great detail in , and its analysis wasrecently further re๏ฌned in . For second order PDEs, a Fortin operatorโs existence and thus discrete stability is ensured if thelocal Riesz representation problems are solved using polynomials of order ๐ = ๐ + ฮ ๐ , where ๐ is the degree of the trial spacediscretization and ฮ ๐ = ๐ the space dimension. However, while this enrichment degree ensures the existence of the requiredFortin operator, numerical evidence suggest that in most cases ฮ ๐ = 1 is typically su๏ฌcient . Alternative test spaces for theDPG method for singular perturbation problems are investigated in , even for the case of ฮ ๐ = 0 .In the AVS-FE method, numerical evidence suggest that ๐ = ๐ is su๏ฌcient for convection-di๏ฌusion PDEs as well asextensive numerical experimentation for the linear elastostatics PDE. Since the test functions are sought in a discontinuouspolynomial space, using ๐ = ๐ still result in a larger space than the trial as the discontinuous spaces contain additional degreesof freedom. Furthermore, in the limit โ โ the space ๐ ( ๎ผ โ ) is essentially ๐ฟ , i.e, any polynomial degree above constants isinherently an enrichment of the test space. IRIK VALSETH
ET AL The approach to establishing AVS-FE approximations described until this point can be established in FE software with relativeease. However, there are other alternative interpretations of DPG methods which are even more straightforward in terms ofimplementation aspects. The inventors of the DPG, Demkowicz and Gopalakrishnan refer to this as di๏ฌerent "hats" of DPGmethods and the one we have explained here is that of a Petrov-Galerkin method with optimal test functions which leadsto (22). As for the DPG method, we are also going to consider another "hat", in which a saddle point, or mixed, interpretationof the AVS-FE method which allows us to employ high level FE solvers such as FEniCS . Hence, let us introduce a newunknown function ( ๐ , ๐ ) , the error representation function. This function derives its name since it is a Riesz representer of theapproximation error induced by the AVS-FE approximation (22) of the weak formulation (5):Find ( ๐ , ๐ ) โ ๐ ( ๎ผ โ ) such that: ( ( ๐ , ๐ ); ( ๐ฏ , ๐ฐ ) ) ๐ ( ๎ผ โ ) = ๐น ( ๐ฏ ) โ ๐ต (( ๐ฎ โ , ๐ โ ); ( ๐ฏ , ๐ฐ )) , โ( ๐ฏ , ๐ฐ ) โ ๐ ( ๎ผ โ ) . (23)Again, the broken nature of the test space allows this function to be approximated on each element ๐พ ๐ โ ๎ผ โ a posteriori to thesolution of (22) to be used as an error estimate and error indicator. Using basic arguments (see, e.g., ) the following saddlepoint system can be established:Find ( ๐ฎ , ๐ ) โ ๐ โ (ฮฉ) โง ( ๐ , ๐ ) โ ๐ ( ๎ผ โ ) such that: ( ( ๐ , ๐ ); ( ๐ฏ , ๐ฐ ) ) ๐ ( ๎ผ โ ) โ ๐ต (( ๐ฎ โ , ๐ โ ); ( ๐ฏ , ๐ฐ )) = โ ๐น ( ๐ฏ ) , โ( ๐ฏ , ๐ฐ ) โ ๐ ( ๎ผ โ ) ,๐ต (( ๐ฒ โ , ๐ณ โ ); ( ๐ , ๐ )) = 0 , โ( ๐ฒ โ , ๐ณ โ ) โ ๐ โ (ฮฉ) . (24)The discretization of ( ๐ , ๐ ) in (24) follows standard FE methodology and the space ๐ ( ๎ผ โ ) is to be discretized with discon-tinuous polynomials. Clearly, the global computational cost of this saddle point system is larger than that of computing optimaltest functions to establish (22). However, this cost is justi๏ฌed as the error representation function is to be used as an a pos-teriori error estimator as well as element-wise error indicators to be used in mesh adaptive strategies. Additionally, the e๏ฌortin implementation into high level FE solvers with well established documentation and capabilities is embarrassingly low. Thenorm equivalence (12) between the energy norm and the norm on ๐ ( ๎ผ โ ) of Riesz representers leads to the follwoing identityfor the error representation function: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B = โ ( ๐ , ๐ ) โ ๐ ( ๎ผ โ ) , (25) EIRIK VALSETH
ET AL which allows us to approximate the approximation error in the energy norm, which is not directly computable due to thesupremum, as well as element-wise error indicators: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B โ โ ( ๐ โ , ๐ โ ) โ ๐ ( ๎ผ โ ) ,๐ = โ ( ๐ โ , ๐ โ ) โ ๐ ( ๐พ ๐ ) , (26)where ( ๐ โ , ๐ โ ) is the approximation of ( ๐ , ๐ ) computed from the discretization of the saddle point system (24). Remark
In this section, we establish a priori error estimates for the AVS-FE method. While we here assume that the components ๐ โ are discretized with continuous polynomials in ๐ถ (ฮฉ) , the analysis can be performed with minor modi๏ฌcations using, e.g., RTor BDM discretizations. Furthermore, we assume that the optimal test functions, to be computed in the approach of (22) aresought in local polynomial spaces of the same degree as the trial functions. Equivalently, we assume that the discrete errorrepresentation function is in discontinuous polynomial spaces of the same degree as the trial functions. We shall use the arbitraryconstant ๐ถ to denote generic mesh independent constants.The starting point of our analysis is the best approximation property of the AVS-FE and DPG methods in terms of the energynorm . Hence, let ( ๐ฎ , ๐ ) โ ๐ (ฮฉ) be the exact solution of the weak formulation (5) and ( ๐ฎ โ , ๐ โ ) โ ๐ โ (ฮฉ) its approximationcomputed from the AVS-FE discretization (22) or equivalently from the saddle point system (24). Then, the energy norm of theapproximation error satis๏ฌes: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B โค โ ( ๐ฎ โ ๐ฏ โ , ๐ โ ๐ฐ โ ) โ B , (27)where ( ๐ฏ โ , ๐ฐ โ ) are arbitrary functions in ๐ โ (ฮฉ) . The proof of this inequality is established using classical techniques fromfunctional analysis and both inf-sup and continuity constants being unity. Additionally, since the energy norm is an equivalentnorm to โ โ โ ๐ (ฮฉ) on the trial space, we consequently have the following quasi-best approximation property: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ ๐ (ฮฉ) โค ๐ถ โ ( ๐ฎ โ ๐ฏ โ , ๐ โ ๐ฐ โ ) โ ๐ (ฮฉ) , (28)where the mesh independent norm equivalence constant ๐ถ depends on the continuity constants of of the bilinear form and aFortin operator . Another key component in the following analysis is the convergence of polynomial interpolating functions. IRIK VALSETH
ET AL Hence, there exist a global polynomial interpolation operator ฮ โ๐ : ฮ โ๐ โถ ๐ โ ๐ โ๐ . (29)Thus, ฮ โ๐ ( ๐ข ) represents an interpolant of ๐ข consisting of piecewise continuous polynomials, then : Theorem ๐ข โ ๐ป ๐ (ฮฉ) and ฮ โ๐ ( ๐ข ) โ ๐ โ๐ be the interpolant of ๐ข (29). Then, there exists ๐ถ > such that the interpolation error can bebounded as follows: โ ๐ข โ ฮ โ๐ ( ๐ข ) โ ๐ป ๐ (ฮฉ) โค ๐ถ โ ๐ โ ๐ ๐ ๐ โ ๐ โ ๐ข โ ๐ป ๐ (ฮฉ) , (30)where โ is the maximum element diameter, ๐ the minimum polynomial degree of interpolants in the mesh, ๐ โค ๐ , and ๐ = min( ๐ + 1 , ๐ ) .To establish error estimates in terms of the energy norm, we ๏ฌrst establish a bound on the Riesz representers of the trialfunctions, i.e., the optimal test functions. Lemma ( ๐ฎ , ๐ ) โ ๐ (ฮฉ) be the exact solution of the AVS-FE weak formulation (5) and ( ๐ฎ โ , ๐ โ ) โ ๐ โ (ฮฉ) its corresponding AVS-FEapproximation from (22). Then: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B โค ๐ถ โ ๐ โ1 ๐ ๐ ๐ฉ โ1 ๐ฎ , (31)where โ is the maximum element diameter, ๐ = min ( ๐ ๐ฎ + 1 , ๐ ) , ๐ ๐ฎ the minimum polynomial degree of approximation of ๐ฎ โ inthe mesh, and ๐ ๐ฉ the regularity of the solution ๐ฉ of the distributional PDE underlying the Riesz representation problem (11) Proof:
The RHS of (27) can be bounded by the error in the Riesz representers of the exact and approximate AVS-FE trialfunctions by the energy norm equivalence in (12), and the map induced by the Riesz representation problem (11) to yield: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B โค โ ( ๐ฉ โ ๐ฉ โ , ๐ซ โ ๐ซ โ ) โ ๐ ( ๎ผ โ ) , (32)where ( ๐ฉ , ๐ซ ) โ ๐ ( ๎ผ โ ) are the exact Riesz representers of ( ๐ฎ , ๐ ) through (11), and ( ๐ฉ โ , ๐ซ โ ) โ ๐ โ ( ๎ผ โ ) are the approximate Rieszrepresenters of ( ๐ฎ โ , ๐ โ ) through a FE discretization of (11). The de๏ฌnition of the norm โ โ โ ๐ and its equivalence to โ โ โ ๐ ( ๎ผ โ ) thengives: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B โค โ ๐ฉ โ ๐ฉ โ โ ๐ป ( ๎ผ โ ) + โ ๐ซ โ ๐ซ โ โ ๐ฟ (ฮฉ) . Since we pick trial functions that are polynomial interpolants for the Riesz representers ( ๐ฉ โ , ๐ซ โ ) of the same degree ๐ ๐ฎ , we bound โ ๐ซ โ ๐ซ โ โ ๐ฟ (ฮฉ) using Theorem 2.1 and โ ๐ฉ โ ๐ฉ โ โ ๐ป ( ๎ผ โ ) by an extension of this theorem for broken Hilbert spaces as introduced EIRIK VALSETH
ET AL by Riviรฉre et al. in to get: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B โค ๐ถ โ ๐ โ1 ๐ ๐ ๐ฉ โ1 ๐ฎ + ๐ถ โ ๐ ๐ ๐ ๐ซ ๐ฎ , (33)where ๐ = min ( ๐ ๐ฎ + 1 , ๐ ๐ฉ ) , ๐ = min ( ๐ ๐ฎ + 1 , ๐ ๐ซ ) , ๐ ๐ฉ the regularity of the Riesz representer of the PDE underlying (11). Sincethe ๏ฌrst term in the RHS of (33) is larger than the second, the proof is completed. โก Next, with the the quasi-best approximation property (28) at hand, we can readily introduce a priori bounds in classicalSobolev norms. First, the bound in terms of the norm โ โ โ ๐ (ฮฉ) is governed by the following lemma: Lemma ( ๐ฎ , ๐ ) โ ๐ (ฮฉ) be the exact solution of the AVS-FE weak formulation (5) and ( ๐ฎ โ , ๐ โ ) โ ๐ โ (ฮฉ) its corresponding AVS-FEapproximation through (22). Then: โ ๐ถ > โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ ๐ (ฮฉ) โค ๐ถ ( โ ๐ โ1 ๐ ๐ ๐ฎ โ1 ๐ฎ + โ ๐ โ1 ๐ ๐ ๐ โ1 ๐ฎ ) , (34)where โ is the maximum element diameter, ๐ = min ( ๐ ๐ฎ + 1 , ๐ ๐ฎ ) , ๐ ๐ฎ the minimum polynomial degree of approximation of ๐ฎ โ , ๐ = min ( ๐ ๐ฎ + 1 , ๐ ๐ ) , in the mesh, ๐ ๐ฎ the regularity of the solution ๐ฎ of the governing PDE (3), and ๐ ๐ the regularity of thesolution ๐ of the governing ๏ฌrst order system PDE (4). Proof:
By the quasi-best approximation property: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ ๐ (ฮฉ) โค ๐ถ โ ( ๐ฎ โ ๐ฏ โ , ๐ โ ๐ฐ โ ) โ ๐ (ฮฉ) , the de๏ฌnition of the norm on ๐ (ฮฉ) (8) leads to: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ ๐ (ฮฉ) โค ๐ถ { โ ๐ฎ โ ๐ฏ โ โ ๐ป (ฮฉ) + โ ๐ โ ๐ฐ โ โ ๐ป ( div , ฮฉ) } , by choosing trial functions ( ๐ฎ โ , ๐ โ ) that are polynomial interpolants and note that โ ๐ โ ๐ฐ โ โ ๐ป ( div , ฮฉ) โค โ ๐ โ ๐ฐ โ โ ๐ป (ฮฉ) , theapproximation property in Theorem 2.1 with ๐ = 1 gives: โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ ๐ (ฮฉ) โค ( ๐ถ โ ๐ โ1 ๐ ๐ ๐ฎ โ1 ๐ฎ + ๐ถ โ ๐ โ1 ๐ ๐ ๐ โ1 ๐ ) , where we use the de๏ฌnitions of Lemma 2.3 for the ๐ โs and ๐ โs. Finally, combining the constants ๐ถ , ๐ถ and noting that in theAVS-FE method we always pick ๐ ๐ฎ = ๐ ๐ the desired error bound is established. โก IRIK VALSETH
ET AL TABLE 1
Material data for the nearly incompressible problem.Property Symbol ValueYoungโs modulus ๐ธ MPaPoissonโs ratio ๐ . Remark for details. To assess the performance of our method, we ๏ฌrst consider a problem with a smooth solution which allows us to asses theconvergence properties of the AVS-FE method to verify the a priori bounds of Section 2.4. Then, we consider an exampleproblem considered by Brenner in , with a manufactured exact solution that is dependent on the Poissonโs ratio which is usedin a comparison between the LSFEM and the Bubnov-Galerkin FE method. As ๏ฌnal numerical veri๏ฌcations, we present twoengineering applications, ๐ ) an example in which a commonly applied engineering structure, a cantilever beam, is consideredand ๐๐ ) the deformation of a composite structure.In the numerical veri๏ฌcations presented in this section we use the FE solvers Firedrake and FEniCS . In particular, forall presented veri๏ฌcations using uniform meshes we use Firedrake, whereas in the case of mesh-adaptive re๏ฌnements, we useFEniCS. In all cases, we use the linear solvert MUMPS to perform the inversion of the resulting sti๏ฌness matrices To present the convergence properties of the AVS-FE for nearly incompressible elastostatics, we consider a 2D model problemwith a smooth exact solution which ensures the stress regularity is ๐ป (ฮฉ) , i.e., we use polynomial approximations for bothvariables. The domain is the unit square, i.e., ฮฉ = (0 ,
1) ร (0 , โ โ , consisting of a material that is nearly incompressiblewith physical properties listed in Table 1, and a sinusoidal exact solution. ๐ฎ ๐๐ฅ ( ๐ฑ ) = โงโชโจโชโฉ ๐ข ๐๐ฅ๐ฅ ( ๐ฑ ) ๐ข ๐๐ฅ๐ฆ ( ๐ฑ ) โซโชโฌโชโญ = โงโชโจโชโฉ sin ( ๐๐ฅ ) sin ( ๐๐ฆ ) sin ( ๐๐ฅ ) sin ( ๐๐ฆ ) โซโชโฌโชโญ . (35)where ๐ is the Lamรฉ parameter: Inspection of (35) reveals the proper boundary conditions are homogeneous Dirichlet conditionson the entire boundary ๐ ฮฉ and the source ๐ is chosen such that it is the di๏ฌerential operator of the PDE (3) acting on ๐ฎ ๐๐ฅ ( ๐ฑ ) . EIRIK VALSETH
ET AL
As an initial veri๏ฌcation, we consider uniform mesh re๏ฌnements starting from a mesh consisting of two triangle elements,and we use ๐ถ (ฮฉ) polynomials of increasing order. With this exact solution, the a priori error estimates from Section 2.4 are: โ ( ๐ข โ ๐ข โ , ๐ โ ๐ โ ) โ ๐ (ฮฉ) โค ๐ถ โ ๐ , โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B โค ๐ถ โ ๐ . (36)In Figures 2 and 3, we present the convergence history for the norms indicated in (36) with the exception of the โ ๐ โ ๐ โ โ ๐ป ( div , ฮฉ) as this is a component of โ ( ๐ข โ ๐ข โ , ๐ช โ ๐ช โ ) โ ๐ (ฮฉ) . In these ๏ฌgures we also show the errors in the individual norms โ ๐ฎ โ ๐ฎ โ โ ๐ฟ (ฮฉ) and โ ๐ฎ โ ๐ฎ โ โ ๐ป (ฮฉ) . The rates of convergence are as predicted in (36) for the energy norm and the norm on ๐ (ฮฉ) . For theindividual norms we observe the expected rates of convergence for polynomial FE approximations with the exception of thecases ๐ = 4 and ๐ = 2 where the observed rates are one order higher in โ ๐ฎ โ ๐ฎ โ โ ๐ฟ (ฮฉ) and โ ๐ฎ โ ๐ฎ โ โ ๐ป (ฮฉ) . This can be seen inFigures 2(a) and 2(b) where the slopes for ๐ = 4 and ๐ = 2 are equal to the slopes for ๐ = 5 and ๐ = 3 , respectively. This problem -7 -6 -5 -4 -3 -2 -1 E rr o r no r m dofs p = 1p = 2p = 3p = 4p = 5 (a) โ ๐ฎ โ ๐ฎ โ โ ๐ฟ (ฮฉ) . -5 -4 -3 -2 -1 E rr o r no r m dofs p = 1p = 2p = 3p = 4p = 5 (b) โ ๐ฎ โ ๐ฎ โ โ ๐ป (ฮฉ) . FIGURE 2
Asymptotic convergence results.was also considered in for the DPG method and we point out that the convergence rates of the AVS-FE approximations matchthe presented rates in . IRIK VALSETH
ET AL E rr o r no r m dofs p = 1p = 2p = 3p = 4p = 5 (a) โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ ๐ (ฮฉ) . -1 E rr o r no r m dofs p = 1p = 2p = 3p = 4p = 5 (b) โ ( ๐ฎ โ ๐ฎ โ , ๐ โ ๐ โ ) โ B . FIGURE 3
Asymptotic convergence results.
To compare the AVS-FE method for nearly incompressible elasticity to other FE methods, we consider a model problem of whichthe solution is a function of the Poissonโs ratio from . The exact solution is given in (37), we use the same Youngโs modulus asin the preceding example. However, we increase the Poissonโs ratio to . to make the problem more challenging witha Lamรฉ parameter ๐ of the order . ๐ฎ ๐๐ฅ ( ๐ฑ ) = โงโชโจโชโฉ ๐ข ๐๐ฅ๐ฅ ( ๐ฑ ) ๐ข ๐๐ฅ๐ฆ ( ๐ฑ ) โซโชโฌโชโญ = โงโชโจโชโฉ sin (2 ๐๐ฆ ) [โ1 + cos (2 ๐๐ฅ )] + sin ( ๐๐ฅ ) sin ( ๐๐ฆ )1 + ๐ sin (2 ๐๐ฅ ) [1 โ cos (2 ๐๐ฆ )] + sin ( ๐๐ฅ ) sin ( ๐๐ฆ )1 + ๐ โซโชโฌโชโญ , (37)The methods we consider are the Galerkin and ๏ฌrst-order system least squares FE methods in addition to the AVS-FE method.For the Galerkin method, we approximate ๐ฎ using ๐ถ (ฮฉ) linear polynomials. Likewise, for the ๏ฌrst-order system least squaresmethod we use linear polynomials for both variables ๐ฎ and ๐ , and in the AVS-FE approximation we use linear polynomialsand ๏ฌrst order Raviart-Thomas bases for ๐ฎ and ๐ , respectively. The LSFEM is implemented in the same fashion as the AVS-FE method, i.e., we consider a saddle point system of the type (24). The initial mesh we consider is uniform, consisting of 8triangular elements, and we proceed to perform uniform mesh re๏ฌnements. EIRIK VALSETH
ET AL
In Figure 4 we compare the convergence history of these three methods in terms of the ๐ฟ (ฮฉ) , and ๐ป (ฮฉ) norms of the erroron ๐ฎ โ ๐ฎ โ . The Galerkin method performs well initially, but upon continued mesh re๏ฌnements becomes unstable, as indicatedby the diverging errors. While the least squares method remains stable initially in the re๏ฌnement process, the error never reachesthe asymptotic range of convergence and begins to increase as the mesh gets very ๏ฌne. Conversely, the AVS-FE method doesnot su๏ฌer from these issues, and retains the optimal rates of convergence in both ๐ฟ (ฮฉ) and ๐ป (ฮฉ) norms. The increasing errorsof the Galerkin and least squares FE methods is likely due to ill conditioning of the resulting sti๏ฌness matrices. By loweringthe Poisson ratio this e๏ฌect of ill conditioning is negated, as expected. Thereby limiting the applicability of these methods fornearly incompressible materials. -4 -3 -2 -1 E rr o r no r m dofs AVS-FEMLSFEMGalerkin (a) โ ๐ฎ โ ๐ฎ โ โ ๐ฟ (ฮฉ) . -1 E rr o r no r m dofs AVS-FEMLSFEMGalerkin (b) โ ๐ฎ โ ๐ฎ โ โ ๐ป (ฮฉ) . FIGURE 4
Comparison between AVS-FE method, ๏ฌrst-order system least squares FE method and Galerkin FE method.
Modern engineering materials often consist of multiple constituents, i.e., composites. Material inclusions often lead to stressconcentrations and need to be accurately resolved to ensure con๏ฌdence in the engineering design. Resolution of such featuresis typically achieved by carefully constructed FE meshes which requires a signi๏ฌcant e๏ฌort from analysts, or mesh re๏ฌnements.Until this point, we have presented numerical veri๏ฌcations in which the mesh partitions are uniform and their re๏ฌnements are
IRIK VALSETH
ET AL ๐ธ ๐ , ๐ ๐ ๐ธ ๐ , ๐ ๐ ๐ฅ๐ฆ . ๐ญ FIGURE 5
Linear elastic problem with inclusion.uniform as well, since the solutions are known to be rather well behaved, a priori . While these certainly give us con๏ฌdence inthe AVS-FE approximations, the computational cost of uniform mesh re๏ฌnements becomes very large as โ โ . To this end,we consider adaptive mesh re๏ฌnements that are guided by the built-in error indicators (26). As an adaptive strategy, we choosethe marking strategy and re๏ฌnement criteria of Dรถr๏ฌer using a ๏ฌxed parameter ๐ = 0 . .As an example of a simpli๏ฌed composite material we consider a unit square domain with an inclusion in the center at ๐ฅ = ๐ฆ = 0 . , as shown in Figure 5 The matrix material is nearly incompressible and the inclusion is a sti๏ฌ material with Youngโsmodulus several orders of magnitude larger than the bulk material. In Table 2 the properties of the materials are listed. Materialswith these properties are, e.g., Silicone based rubber for the matrix and an epoxy based inclusion. The boundary conditions are TABLE 2
Material properties.Physical property Youngโs Modulus Poissonโs Ratio.Matrix 1,500 MPa 0.49Inclusion 10,000 MPa 0.3as shown in Figure 5 i.e., a surface traction ๐ญ = {100 MPa , ๐ on the right side and a fully clamped boundary on a portionof the left side. Hence, we expect signi๏ฌcant concentration of stresses near the inclusion. We employ the error representationfunction (26) as an a posteriori error estimate. The goal of the adaptive algorithm is then to minimize this error based on itslocal error indicators.To establish an initial mesh taking into account the circular geometry with reasonable accuracy, we employ the built-in FEn-iCS mesh generation functionality. In particular, we use the tool "mshr" to create an inclusion consisting of 250 line segmentsto accurately represent the circumference of the circle. The initial mesh is shown in Figure 6. For the AVS-FE approximations, EIRIK VALSETH
ET AL
FIGURE 6
Initial mesh for the inclusion problem.we consider second order approximations for all variables in both cases of RT and ๐ถ (ฮฉ) stress approximations.The normal stresses ๐ ๐ฅ๐ฅ and ๐ ๐ฆ๐ฆ are presented in Figure 7 where the stresses, as expected, are concentrated near the sti๏ฌinclusion. Furthermore, we see in this ๏ฌgure that the stresses perpendicular to the loading, i.e., ๐ ๐ฆ๐ฆ are dominated by the Poissone๏ฌect which is another indication that the AVS-FE approximation of this problem is consistent with the expected physics. Asa sanity check, we also note that the far ๏ฌeld stresses are equal to the applied surface traction ๐ญ , see Figure 7(a). The shear (a) ๐ โ๐ฅ๐ฅ ( ๐ฑ ) . (b) ๐ โ๐ฆ๐ฆ ( ๐ฑ ) . FIGURE 7
AVS-FE approximate normal stress components (MPa).stress component ๐ ๐ฅ๐ฆ is shown in Figure 8 along with the adapted mesh after 12 re๏ฌnements. The ๏ฌnal adapted mesh shown in๏ฌgure 8(b) highlights that the built-in error indicators as well as the marking strategy of Dรถr๏ฌer lead to mesh re๏ฌnements inlocations critical to proper resolution of physical features. IRIK VALSETH
ET AL (a) ๐ โ๐ฅ๐ฅ ( ๐ฑ ) . (b) Final adapted mesh using RT stress approximations. FIGURE 8
AVS-FE approximate shear stress component (MPa) and ๏ฌnal mesh.Computing these results in the AVS-FE method is straightforward due to its built-in error estimate and discrete stability.Hence, the e๏ฌort in implementation of the adaptive algorithm is minimal and requires only a FE solver with mesh re๏ฌnementcapabilities. However, it is worth mentioning that the required e๏ฌort to establish similar results in the Galerkin FE methodwould be signi๏ฌcant. While several a posteriori error estimation techniques exists, the stability issue that arises with the jumpin material coe๏ฌcients will likely require signi๏ฌcant e๏ฌorts in analysis to ensure stability of the error estimation or prohibitivelylarge mesh generation e๏ฌorts.Finally, we compare ๐ถ (ฮฉ) and RT stress approximations. The inclusion leads to a stress ๏ฌeld that has continuous normalcomponents across the interface between the two materials, whereas the tangential components are discontinuous. While thisdomain is still convex, we therefore expect that the generally advocated ๐ถ (ฮฉ) stress approximations in the AVS-FE methodwill be less accurate than the standard mixed FE choice of RT discretization. Thus, as a veri๏ฌcation we consider both types ofstress approximations in this experiment and measure the di๏ฌerence between the two by the approximate energy norm. To thisend, we employ the same adaptive algorithm for the case of ๐ถ (ฮฉ) stresses and in Figure 9, the convergence histories of theapproximate error in the energy norm through (25) for both cases are presented. As expected, in this case the RT approximationslead to lower energy norm errors due to the pollution e๏ฌect of enforcing continuity of tangential stress components for the ๐ถ (ฮฉ) case. Finally, in Figure 10 the ๏ฌnal adapted mesh for the case of ๐ถ (ฮฉ) stresses are presented. The overly restrictive stressapproximations lead to mesh re๏ฌnements that do not resolve the stress ๏ฌeld near the inclusion in the same was as the RT case asevident from comparison of Figures 8(b) and 10. Hence, the ๐ถ (ฮฉ) stress approximations are not applicable for this problem. EIRIK VALSETH
ET AL -1 E rr o r no r m dofs RTC FIGURE 9
Convergence of the energy norm comparing ๐ถ (ฮฉ) and RT stress approximations. FIGURE 10
Final adapted mesh using ๐ถ (ฮฉ) stress approximations. To further illustrate the AVS-FE method in its application to FE analyses of nearly incompressible solids, we present a commonengineering application of non-uniform bending of beams. The 2D problem concerns a beam with a length ๐ฟ = 2 m and slender-ness ratio ๐ฟ โ ๐ป = 10 , and consists of a homogeneous linearly elastic isotropic material with a Youngโs Modulus equal to thatof a nitrile based rubber, i.e., ๐ธ = 1 . MPa. The Poisson ratio ๐ is chosen such that โ ๐ = 10 โ9 and therefore the material isnearly incompressible. The beam is subject to kinematic constraints along its left edge, where material particles are prohibitedfrom moving in the ๐ฅ direction but free to move in the ๐ฆ direction. To eliminate the rigid body translation in the ๐ฆ direction, the IRIK VALSETH
ET AL ๐๐ฟ ๐ฅ๐ฆ๐ป FIGURE 11
Linear elastic beam problem.bottom left corner point is kept ๏ฌxed. In terms of loading, the beam is subject to a downward uniform distributed force ๐ , of . N/m as depicted in Figure 11).For our numerical experiment of the AVS-FE method, we start with an initial uniform mesh of four triangular elements (twoelements in the length, versus one in the height) and subsequently conduct uniform โ -re๏ฌnements in which each element ispartitioned into four new elements. In all computations, we employ quadratic (i.e., ๐ = 2 ) ๐ถ trial functions for the displacementsand second order Raviart-Thomas functions functions for the stress variable. The evolution of the elastic strain energy of theAVS-FE solutions throughout the โ -re๏ฌnement process are shown in orange in Figure 12. In comparison, the correspondingresults for the classical FE, or Bubnov-Galerkin, method are also shown in this graph in blue. These were established by usingthe same mesh partitions as for the AVS-FE method but employing the classical ๐ถ quadratic Lagrangian trial functions for thedisplacements and a standard displacement based weak formulation, as is common in classical Bubnov-Galerkin analyses.Figure 12 shows that, as expected for this level of near incompressibility, the classical FE solutions fail to converge and startto exhibit spurious behavior as โ โโ with greatly changing values for the elastic strain energy between successive solutions.The AVS-FE method, on the other hand, maintains numerical stability from the onset and exhibits convergence upon continueduniform mesh re๏ฌnements. The converged AVS-FE solutions are physically valid, as can be seen in Figure 13, in which contourplots of the distributions of the normal stress ๐ ๐ฅ๐ฅ (Figure 13(a)) and shear stress ๐ ๐ฅ๐ฆ (Figure 13(b)) are shown of the AVS-FEsolution for the mesh consisting of
22 ร 11 elements. Since Raviart-Thomas trial functions have been used to compute the stressvariables, minor discontinuities across inter-element edges can be observed, but these attenuate as the mesh is further re๏ฌned.Hence, the AVS-FE method shows less sensitivity to the nearly incompressible constitutive behavior of the material than theclassical FE method. EIRIK VALSETH
ET AL
FIGURE 12
Non-uniform bending of a nearly incompressible beam โ elastic strain energy evolution. (a) Normal Stress ๐ ๐ฅ๐ฅ (kPa).(b) Shear stress ๐ ๐ฅ๐ฆ (kPa). FIGURE 13
AVS-FE stress distributions of beam in bending with โ ๐ = 10 โ9 . IRIK VALSETH
ET AL We have presented a mixed FE method that results in continuous FE approximations of both the displacement and (normal)stress ๏ฌelds in linear elasticity. In particular, we considered nearly incompressible materials as classical FE methods su๏ฌerfrom loss of discrete stability for such materials. The DPG philosophy of optimal test spaces leads to stable FE approximationswithout reformulation of the problem using the compliance tensor for the case for the nearly incompressible case. The DPGmethod we present distinguishes itself from the DPG method by only breaking the test space while keeping the trial spaceglobally con๏ฌrming. Hence, in the corresponding FE approximation we use classical FE bases such as Lagrange polynomialsand Raviart-Thomas functions.We present a priori error bounds in terms of norms of the numerical approximation errors of both displacement and stress trialvariables. These bounds are are all optimal in the sense that the approximation errors converge at rates equal to their underlyinginterpolating functions. Numerical veri๏ฌcations of the asymptotic convergence properties con๏ฌrm the established error boundin all appropriate norms. A convergence study comparing our method to existing FE methods, i.e., Galerkin FE method and leastsquares FE method, show that the AVS-FE method is superior for the case of nearly incompressible materials for the presentedveri๏ฌcation. In the veri๏ฌcation presented, the LSFEM reached a point in which the presented error norms began to diverge. TheAVS-FE method did not exhibit this behavior, nor have we observed such behavior for other veri๏ฌcations we have performed.We attribute this to the scaling term โ ๐ of the ๐ป seminorm portion of โ โ โ ๐ ( ๎ผ โ ) (8) as it ensures entries in the resulting sti๏ฌnessmatrix are of similar order of magnitude. However, we cannot rule out such behavior for the AVS-FE method for very ๏ฌne meshesas this is an issue of numerical linear algebra and not the LSFEM or AVS-FE method. Similar behavior has been observed byStorn in for the DPG method, where it is attributed to the inversion of the Gram matrix in the computation of optimal testfunctions.By considering a global saddle point form of the AVS-FE method, we establish both approximations of the displacement anstress ๏ฌelds as well as an approximation of an error representation function which measures the global energy error of the AVS-FE approximation. This error representation leads to a posteriori error estimate and error indicators which we employ in a meshadaptive strategy. We present a numerical veri๏ฌcation of a challenging physical application of a composite material where thebuilt-in error indicator is used to drive adaptive mesh re๏ฌnements to resolve the stress ๏ฌeld in the composite.While successful for the presented composite material, the built-in error indicator is a local indication of the residual of theAVS-FE approximation (see (25)). However, in certain applications, localized solution features may be of higher importancethan the global energy error. Hence, goal-oriented error estimates and error indicators based on local quantities of interest canprovide alternative mesh re๏ฌnement strategies. As shown in , alternative AVS-FE goal-oriented error estimates are capableof accurately estimating these errors and driving goal-oriented mesh re๏ฌnements. While we have considered a single AVS-FE EIRIK VALSETH
ET AL weak formulation here (5), as mentioned in Remark 2.1, this is not a unique choice. In , multiple weak formulations for linearelasticity are considered for the DPG method, all of which provide slightly di๏ฌerent FE approximations. Such investigation forthe AVS-FE method and comparison of the DPG and AVS-FE methods for linear elasticity are postponed to future works. ACKNOWLEDGEMENTS
Authors Dawson and Valseth have been supported by the United States National Science Foundation - NSF PREEVENTS Track2 Program, under NSF Grant Number 1855047. Authors Kaul, Romkes, and Valseth have been supported by the United StatesNational Science Foundation - NSF CBET Program, under NSF Grant Number 1805550.
References
1. Babuลกka I, Suri M. Locking e๏ฌects in the ๏ฌnite element approximation of elasticity problems.
Numerische Mathematik
SIAM Journal on Numerical Analysis
An introduction to the mathematical theory of ๏ฌnite elements . Dover Publications . 2012.4. Phillips PJ, Wheeler MF. Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach.
Computational Geosciences
Mixed and Hybrid Finite Element Methods . 15. Springer-Verlag . 1991.6. Falk RS. Finite element methods for linear elasticity. In: Springer. 2008 (pp. 159โ194).7. Arnold DN, Winther R. Mixed ๏ฌnite elements for elasticity.
Numerische Mathematik
Mathematical models and methods in appliedsciences
IMA Journal of Numerical Analysis
Computer Methods in Applied Mechanics and Engineering
IRIK VALSETH
ET AL
11. Ambartsumyan I, Khattatov E, Nordbotten JM, Yotov I. A multipoint stress mixed ๏ฌnite element method for elasticity onsimplicial grids.
SIAM Journal on Numerical Analysis
Computer methods in applied mechanics and engineering
Computational Mechanics
Journal of AppliedMechanics, Transactions ASME
Computer Methods in Applied Mechanics and Engineering
Computer methods in applied mechanics and engineering
Computers & structures
Least-Squares Finite Element Methods . 166. Springer Science & Business Media . 2009.19. Demkowicz L, Gopalakrishnan J. A Class of Discontinuous Petrov-Galerkin Methods. II. Optimal Test Functions.
Numerical Methods for Partial Di๏ฌerential Equations
SIAM journal onnumerical analysis
Numerical Methods for Partial Di๏ฌerential Equations: An International Journal
Mathematics of Computation
Numerische Mathematik
Journal of Computational Physics EIRIK VALSETH
ET AL
25. Keith B, Fuentes F, Demkowicz L. The DPG methodology applied to di๏ฌerent variational formulations of linear elasticity.
Computer Methods in Applied Mechanics and Engineering
Computer methods in applied mechanics and engineering โ, ๐, ๐
Mathematical and Computational Finite Element Framework for Boundary Valueand Initial Value Problems.
Acta Mechanica Solida Sinica
International Journal for Numerical Methods in Engineer-ing
Automatic Variationally Stable Analysis for Finite Element Computations . PhD thesis. South Dakota School ofMines and Technology, Rapid City, South Dakota; 2019.31. Girault V, Raviart PA. Finite Element Methods for Navier-Stokes Equations; Theory and Algorithms. In: . 5. Springer-Verlag. 1986.32. Carstensen C, Demkowicz L, Gopalakrishnan J. A Posteriori Error Control for DPG Methods.
SIAM Journal on NumericalAnalysis
Computers & Mathematics with Applications
Numerische Mathematik
Computers and Mathematics with Applications
Mixed ๏ฌnite element methods and applications . 44. Springer . 2013.37. Nagaraj S, Petrides S, Demkowicz LF. Construction of DPG Fortin operators for second order problems.
Computers &Mathematics with Applications
Computers & Mathematics with Applications
IRIK VALSETH
ET AL
39. Salazar J, Mora J, Demkowicz L. Alternative enriched test spaces in the DPG method for singular perturbation problems.
Computational Methods in Applied Mathematics
ICES Report
Archive of Numerical Software โ๐ version of the ๏ฌnite element method with quasiuniform meshes. ESAIM: Mathematical Modellingand Numerical Analysis-Modรฉlisation Mathรฉmatique et Analyse Numรฉrique
Computational Geosciences
SIAMJournal on Numerical Analysis
ACMTransactions on Mathematical Software (TOMS)
SIAM Journal on Matrix Analysis and Applications
SIAM Journal on Numerical Analysis
Computers & Mathematicswith Applications