A stabilized, nodally integrated linear tetrahedral is formulated and analyzed. It is well known that linear tetrahedral elements perform poorly in problems with plasticity, nearly incompressible materials, and acute bending. For a variety of reasons, linear tetrahedral elements are preferable to quadratic tetrahedral elements in most nonlinear problems. Whereas, mixed methods work well for linear hexahedral elements, they don't for linear tetrahedrals. On the other hand, automatic mesh generation is typically not feasible for building many 3D hexahedral meshes. A stabilized, nodally integrated linear tetrahedral is developed and shown to perform very well in problems with plasticity, nearly incompressible materials and acute bending. Furthermore, the formulation is analytically and numerically shown to be stable and optimally convergent. The element is demonstrated to perform well in several standard linear and nonlinear benchmarks

Puso, M

UNT Digital Library

UCRL-CONF-211697An Improved Linear TetrahedralElement for PlasticityMichael PusoApril 26, 20052005 Joint ASME/ASCE/SES Conference on Mechanics andMaterialsBaton Rouge, LA, United StatesJune 1, 2005 through June 3, 2005Disclaimer   This document was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor the University of California nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or the University of California. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or the University of California, and shall not be used for advertising or product endorsement purposes.  Proceedings of McMat 20052005 Joint ASME/ASCE/SES Conference on Mechanics and MaterialsJune 1-3, 2005, Baton Rouge, Louisiana, USA743AN IMPROVED LINEAR TETRAHEDRAL ELEMENT FOR PLASTICITY1Michael PusoLawrence Livermore National LaboratoryDepartment of Mechanical EngineeringUniversity of CaliforniaLivermore, California 94551Email: puso@llnl.govABSTRACTA stabilized, nodally integrated linear tetrahedral is formu-lated and analyzed. It is well known that linear tetrahedral ele-ments perform poorly in problems with plasticity, nearly incom-pressible materials, and acute bending. For a variety of reasons,linear tetrahedral elements are preferable to quadratic tetrahe-dral elements in most nonlinear problems. Whereas, mixed meth-ods work well for linear hexahedral elements, they don’t for lin-ear tetrahedrals. On the other hand, automatic mesh genera-tion is typically not feasible for building many 3D hexahedralmeshes. A stabilized, nodally integrated linear tetrahedral is de-veloped and shown to perform very well in problems with plas-ticity, nearly incompressible materials and acute bending. Fur-thermore, the formulation is analytically and numerically shownto be stable and optimally convergent. The element is demon-strated to perform well in several standard linear and nonlinearbenchmarks.INTRODUCTION1 A nodally integrated tetrahedral element was first formu-lated in [1] and further analyzed in [2]. The element was shownto perform well in several 2D bending problems but it was notedthat the formulation was prone to spurious low energy modes. Inthis work, a stabilization of the nodally integrated tetrahedral isproposed. It is then shown analytically that the proposed schemeis stable and consistent for linear elasticity. As it turns out, the1Work performed under the auspices of the U.S. Department of Energy byLawrence Livermore National Laboratory under Contract W-7405-Eng-48.nodal integration provides reduced constraints such that coarsemesh accuracy can be achieved for nearly incompressible mate-rials and plasticity. Finally, example problems demonstrate theimproved numerical response of the proposed tetrahedral formu-lation.FORMULATIONIn what follows, superscript indices in lower case (typicallye) refer to element quantities and superscript indices in uppercase (typically I) refer to nodal quantities. Here, the standardweak form of linear elasticity is considereda(u,w) = f (w) u,w ∈ (H1(Ω))3 (1)for the domain Ω ⊂ R3. The discrete problem readsah(uh,wh) = f (wh) (2)where uh wh are the discrete trial and test functions respectivelybased on linear tetrahedral interpolation and the discrete bilinearoperator ah(uh,wh) is based on average nodal integration. Thedisplacement gradient on the domain of element e (i.e. Ωe) willbe written ∇ueh = {∇uh : x ∈ Ωe} and is observed to be constanton Ωe. Furthermore, strains based on the trial and test functionsare respectively written εe = ∇sueh and δεe = ∇swe where ∇s isthe symmetric gradient. In order to present the average nodal1strain formulation, the set SI is defined to be the group of ele-ments common to node I and Ne and Nn are respectively the totalnumber of elements and nodes defined on Ω. The average nodalgradient ∇uI and strain εI at node I are defined∇uIh =1V I∑e∈SIV e4∇uehεI =12[∇uIh +(∇uIh)T ] or εI =1V I∑e∈SIV e4εe (3)where V e = vol(Ωe) and the average nodal volume is givenV I =∑e∈SIV e4(4)”Virtual” nodal strain quantities analogous to (3) based on thetest functions are also defined∇wIh =1V I∑e∈SIV e4∇weδεI = 12[∇wIh +(∇wIh)T ] or δεI =1V I∑e∈SIV e4δεe (5)Considering the standard tetrahedral finite element discretizationfor the equations of motion and applying the preceeding notation,the bilinear form can be written as a sum over elements or a sumover nodes as followsa(wh,uh) =Ne∑e=1V eδεe : Cεe (6)=Nn∑I=1∑e∈SIV e4δεe : Cεe (7)where, again, SI is the set of elements that are common to node Iand C is the material stiffness.The following modified, nodal strain definition of the bilin-ear form is proposedah(wh,uh)=Nn∑I=1V IδεI :CεI +Nn∑I=1∑e∈SIαe,IV e4(δεI−δεe) : ˜C(εI−εe)(8)where αe,I is a stabilization parameter that can potentially de-pend on element e and node I and ˜C could be an alternate ma-terial stiffness. Of course C = ˜C would be the natural choicebut given a nearly incompressible material with Poisson’s rationν = 0.4999, better results may be obtained by letting Poisson’sratio be ν = 0.4 in ˜C. The average nodal strain formulation pro-posed in [1] and analyzed in [2] is recovered when αe,I = 0. Infact, results from aforementioned references are quite good de-spite the absence of the stabilization term. Nonetheless, the eige-nalysis in the RESULTS Section demonstrate the necessity of thestabilization term.STABILITYThe strain energy energy from (8) is writtenah(uh,uh)=12Nn∑I=1V IεI :CεI + 12Nn∑I=1∑e∈SIαe,IV e4(εI−εe) : ˜C(εI−εe)(9)With αe,I = 0, a spurious zero energy mode can arise in the eventof an infinite domain with a regular lattice of points where thedisplacement oscillates such that strains εe are positive to the leftof every node and negative to the right of every node thus pro-ducing an average strain εI = 0 at every node and thus a zeroenergy mode. For a finite domain, the zero energy mode can-not propagate since the nodes on the boundaries have strain con-tributions from only one side thus precluding the zero energymode. Nonetheless, the energy of this above described mode canbe very small and stability in the H1 norm is not ensured. Withthe stabilization, the energy in (9) reduces to the following whenεI = 0ah(uh,uh) =12Nn∑I=1∑e∈SIαe,IV e4εe : Cεe (10)From Korn’s lemma (10) is always non zero for any displacementfield modulo rigid body modes such thatah(uh,uh)≥ α‖uh‖21 (11)where α is a mesh independent constant. Consequently, the form(9) has no spurious zero energy modes and is thus V coercive. Inthe event that αe,I = 0, the discrete bilinear form may be greaterthan zero for non-trivial displacements but α in (11) will dependon h precluding it from the following convergence analysis.CONVERGENCEIn this section, optimal convergence in the H1 (energy norm) isproved analytically for the proposed discrete bilinear form in (8).Repeated indices will imply summation unless otherwise speci-fied and the following notation will be used throughout for norms2and semi-norms‖y‖1 =(∫Ω(yiyi + yi, jyi, j)dΩ)1/2|y|1 =(∫Ωyi, jyi, j dΩ)1/2(12)and the Frobenius norms will be used for matrices where for ex-ampley = Rm,n ‖y‖= (yi jyi j)1/2 i = 1 : m, j = 1 : nStrang’s first Lemma is defined‖u−uh‖1 ≤C infvh∈Vh[‖u− vh‖1 + supwh∈Vha(vh,wh)−ah(vh,wh)‖wh‖1 ](13)where u is the exact solution to the boundary value problem anduh is the solution to the approximate weak form ah(wh,uh) =f (wh) Because (8) is coercive (i.e. stable) (13) can be used toestablish its rate of convergence. From hereon, letvh = Πhu (14)be the interpolant of u thus from approximation theory [3] ‖u−vh‖1 ≤Ch. Hence, it will suffice to show thatsupwh∈Vha(vh,wh)−ah(vh,wh)‖wh‖1 <Ch (15)Exploiting the symmetries of the elasticity tensor C the exact bi-linear form (7) can be written in terms of the gradientsa(wh,vh) =Nn∑I=1∑e∈SIV e4∇we : C∇ve (16)Furthermore, the modified bilinear form (8) can be rewrittenah(wh,vh) = anh(wh,vh)+ash(wh,vh) (17)whereanh(wh,vh) =Nn∑I=1∑e∈SIV e4∇we : C∇vIh (18)ash(wh,vh) =Nn∑I=1∑e∈SIαe,IV e4(∇wI −∇we) : ˜C(∇vI −∇ve) (19)and where (5) was used to eliminate δεI in (18). Our first task isto show thatsupwh∈Vha(vh,wh)−anh(vh,wh)‖wh‖1 <Ch (20)Equations (16) and (18) provide the following differencea(vh,wh)−anh(vh,wh) =Nn∑I=1∑e∈SIV e4∇we : C(∇ve−∇vIh) (21)and by Cauchy Schwarz, the following inequality can be estab-lisheda(vh,wh)−anh(vh,wh)≤ MNn∑I=1∑e∈SIV e4‖∇we‖‖∇ve−∇vIh‖(22)where ‖ · ‖ is the norm ‖∇v‖ = (vi, jvi, j)1/2 as per (12) and M isthe norm of the elasticity tensor M = (Ci jklCi jkl)1/2. In order tobound ‖∇ve−∇vIh‖ in (22), it is required that the exact solution uis sufficiently smooth (u ∈C2) such that following Taylor seriesexpansion can be made about point xe∇u(x) = ∇u(xe)+∇2u(xe∗)(x− xe) (23)where xe is the coordinates for the centroid of the element eand xe∗ ∈ Ωe. For a 3D tetrahedral e, xe corresponds to ξ =(1/3,1/3,1/3) in the parent coordinates of the tetrahedral. Next,the following difference is computed using (3) and (4)∇ve−∇vI = ∇ve−∇u(xI)− 1V I∑e¯∈SIV e¯4(∇ve¯−∇u(xI)) (24)where xI is the coordinates of node I. Expanding ∇u(xI) in termsof (23) and substituting into (24) yields∇ve−∇vI = ∇ve−∇u(xe)−∇2u(xe∗)(xI − xe)− 1V I∑e¯∈SIV e¯4(∇ve¯−∇u(xe¯))+∑e¯∈SIV e¯4(∇2u(xe¯∗)(xI − xe¯) (25)Taking the norm of the left and right sides of (25) and applyingthe triangle inequality yields the inequality‖∇ve−∇vI‖ ≤ ‖∇ve−∇u(xe)‖+‖∇2u(xe∗)‖h+1V I∑e¯∈SIV e¯4‖(∇ve¯−∇u(xe¯))‖+ 1V I∑e¯∈SIV e¯4‖∇2u(xe¯∗)‖h (26)3where h is the maximum element diameter and thus ‖xI−xe‖≤ hfor all e. It turns out to be most convenient to bound (26) by theL∞ norm; consequently, for some vector u ∈ R3|u|2,∞,Ωe = maxi jksupx∈Ωe∣∣∣∣ ∂2ui∂x j∂xk∣∣∣∣ (27)Based on (14)and approximation theory [3,4], the following errorbound is exploited∣∣∣∣ ∂ui∂x j (x)−∂vi∂x j∣∣∣∣x∈Ωe≤Ch|u|2,∞,Ωe ∀ i = 1 : 3 (28)Using (3) and substituting‖∇2u(xe)‖=∑i jk( ∂2ui∂x j∂xk(xe))21/2≤ 3√3|u|2,∞,Ωe (29)and (28) into (26) yields‖∇ve−∇vI‖ ≤ (Ce +3√3)h |u|2,∞,Ωe+1V I∑e¯∈SIV e¯4(Ce¯ +3√3)h |u|2,∞,Ωe¯≤C h |u|2,∞,Ω (30)since |u|2,∞,Ω| ≥ |u|2,∞,Ωe ∀e. Substituting the result (30) into(22) and applying the Cauchy Schwartz inequality yieldsNn∑I=1∑e∈SIV e4‖∇we‖‖∇ve−∇vI‖ ≤ Nn∑I=1∑e∈SIV e4‖we‖21/2× Nn∑I=1∑e∈SIV e4C h2 |u|22,∞,Ω1/2≤( Ne∑e=1V e ‖we‖2)1/2( Ne∑e=1V e C h2 |u|22,∞,Ω)1/2≤ C |wh|1,Ω |u|2,∞,ΩV h(31)where summations over nodes I and sets SI were reorganizedas summations over elements e as in (7). and the semi-norm iscomputed|wh|1,Ω =( Ne∑e=1V e ‖∇we‖2)1/2(32)since the ∇we is constant on Ωe. The desired bound in (20) is aconsequence of (31).The next task is to show thatsupwh∈Vhash(vh,wh)‖wh‖1 <Ch (33)Applying the results from (30),(31) and the Cauchy Schwarz andtriangle inequalities to (19) provides the following boundash(wh,vh)≤ MNn∑I=1∑e∈SIαe,IV e4‖∇wI −∇we‖‖∇ve−∇vI‖≤ αˆMNn∑I=1∑e∈SIV e4(‖∇wI‖+‖∇we‖)‖∇ve−∇vI‖≤ αˆMNn∑I=1∑e∈SIV e4‖∇wI‖(C1 |u|2,∞,Ω h)+ αˆMC2|w|1,Ω |u|2,∞,ΩV h(34)where αˆ = maxe,I αe,I . At this point it remains to bound the re-maining summation from the last inequality in (34). Using thedefinition for nodal volume (4) and nodal strain (3) givesNn∑I=1∑e∈SIV e4‖∇wI‖ ≤Nn∑I=1∑e∈SIV e4‖ 1V I∑eˆ∈SIV eˆ4∇weˆ‖≤Nn∑I=1∑e∈SIV e4V I∑eˆ∈SIV eˆ4‖∇weˆ‖≤Nn∑I=1∑eˆ∈SIV eˆ4‖∇weˆ‖≤ Nn∑I=1∑eˆ∈SIV eˆ4‖∇weˆ‖21/2 Nn∑I=1∑eˆ∈SIV eˆ41/2≤ |wh|1,ΩV (35)The desired bound in (33) is a arrived at by substituting the resultof (35) into (34). Substitution of (20) and (33) into (13) providesthe final result‖u−uh‖1 ≤Ch (36)provided sufficient smoothness of the exact solution u ∈ C2 ∩(H1)3.4Table 1. EIGENFREQUENCIES FROM TET AND HEX MESHESModes α = 0.05 α = 0 hexMode 1 0.258 0.209 0.258Mode 10 0.424 0.236 0.404Mode 16 0.452 0.248 0.482Figure 1. FIRST EIGENMODE: α = 0.05 (LEFT) AND α = 0 (RIGHT)RESULTSThe following example problems demonstrate the necessityof the stabilization added to the nodal integration scheme and thegood convergence characteristics of the proposed approach. Auniform stabilization parameter was employed throughout suchthat αe,I = α with α = 0.05EigenanalysisAn eigenanalysis reveals the spurious modes of the nodal in-tegration approach (i.e. α = 0). The eigenvalues and first eigen-mode of a 1x1x1 block with ρ = 1, E = 1 and ν = 0.499 areshown in Table 1 and Fig. 1 respectively. The eigenfreqenciesare compared to a mesh composed of 512 incompatible modeshexahedral elements. With α = 0, the spurious modes are notzero energy modes since the method is still stable in the L2 normbut not the H1 norm. On the other hand, it is seen that furthermesh refinement yields convergence to the wrong eigenfrequen-cies with α = 0. Whereas, the stabilized formulation convergesto the correct eigenfrequencies.Asymptotic ErrorA standard benchmark is a cantilever beam ν= 0.499 loadedin shear as described in [5]. The energy of the discretization erroris plotted in Fig. 2 for the standard linear triangle element and thenodally integrated triangle with α = 0.05. The nodally integrated-3.5-3-2.5-2-1.5-1-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6linearnodal (0.05)log(e)log( h )11Figure 2. DISCRETIZATION ERROR FOR CANTILEVER BEAMelement actually converges at a rate of 1.2 and is clearly superiorto the linear triangle.PlasticityThe Cook membrane Fig. 3 (see [6] for dimensions) withVon Mises plasticity (E = 70, ν = 0.333, σy = 0.243 and a lin-ear hardening modulus Et = 0.15) is considered. Meshes withn= 8,16,32,64 elements along the edges are considered and Fig.3 shows the plastic strain for the n = 16 case. The stabilizationstiffness ˜C (8) was chosen as follows: µ˜ = αEt/(2(1+ ν)) and˜λ = αλ where λ = E/((1+ν)(1− 2ν)) with α = 0.05. The tipdisplacement versus n is shown plotted in Fig. 4 for the nodallyintegrated triangle, linear triangle and the QM6 incompatiblemodes quadrilateral. The linear tetrahedral is very stiff whereasthe nodally integrated triangle performs as well or better than theincompatible modes quadrilateral.DISCUSSIONA stabilized nodally integrated tetrahedral element formula-tion was developed. It was shown analytically and numericallythat the proposed formulation was stable and optimally conver-gent. Although the element is not shown to be LBB stable, it per-forms well in some cases where nearly incompressible or plasticmaterial were used. More studies will be done for more generalmaterial cases such as nonlinear hardening and large deforma-tions.5Plastic strain3.60.00.18Figure 3. COOK MEMBRANE: EFFECTIVE PLASTIC STRAIN n = 16REFERENCES[1] Dohrmann, C., Heinstein, M., Jung, J., Key, S., andWitkowski, W., 2000. “Node-based uniform strain elementsfor three-node triangular and four-node tetrahedral meshes”.International Journal for Numerical Methods in Engineer-ing, 47, pp. 1549–1568.[2] Bonet, J., Marriot, H., and Hassan, O., 1992. “An aver-aged nodal deformation gradient linear tetrahedral elementfor large strain explicit dynamic applications”. Communi-cations in Numerical Methods in Engineering, 17, pp. 551–561.[3] Ciarlet, P., 2002. The finite element method for elliptic prob-lems. North-Holland, Amsterdam.[4] Johnson, C., 1987. Numerical solution of partial differen-tial equations by the finite element method. North-Holland,Amsterdam.[5] Hughes, T., 1987. The Finite Element Method: Linear Staticand Dynamic Finite Element Analysis. Prentice-Hall, Engle-wood Cliffs, NJ.[6] Puso, M., 2003. “A highly efficient enhanced assumed strainphysicall stabilized hexahedral element”. International Jour-nal for Numerical Methods in Engineering, 49, pp. 1024–1064.0.511.522.50 1 0 2 0 3 0 4 0 5 0 6 0 7 0nodallinearQM 6tip displacementnFigure 4. TIP DISPLACEMENT VERSUS ELEMENTS PER SIDE6

An Improved Linear Tetrahedral Element for Plasticity

https://digital.library.unt.edu/ark:/67531/metadc885686/m2/1/high_res_d/918409.pdf

An Improved Linear Tetrahedral Element for Plasticity

Abstract

Similar works

Full text

Available Versions

UNT Digital Library