A new stabilized nodal integration scheme is proposed and implemented. In this work, focus is on the natural neighbor meshless interpolation schemes. The approach is a modification of the stabilized conforming nodal integration (SCNI) scheme and is shown to perform well in several benchmark problems

Chen, J S

Puso, M

Zywicz, E

UNT Digital Library

UCRL-PROC-218895A New Stabilized NodalIntegration ApproachM. Puso, E. Zywicz, J. S. ChenFebruary 13, 2006Third International Workshop Meshfree Methods for PartialDifferential EquationsBonn Germany, GermanySeptember 12, 2005 through September 15, 2005This 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.  Updated October 14, 2003 A New Stabilized Nodal Integration ApproachMichael Anthony Puso1, Edward Zywicz2, and J.S. Chen31 Lawrence Livermore National Laboratory, Livermore CA 94550 puso@llnl.gov2 Lawrence Livermore National Laboratory, Livermore CA 94550zywicz@llnl.gov3 University of California at Los Angeles, Los Angeles CA 90095jschen@seas.ucla.eduSummary. A new stabilized nodal integration scheme is proposed and imple-mented. In this work, focus is on the natural neighbor meshless interpolationschemes. The approach is a modification of the stabilized conforming nodal integra-tion (SCNI) scheme and is shown to perform well in several benchmark problems.Key words: Nodal integration, natural neighbor, particle method1 IntroductionNodal integration of the Galerkin weak form of the equations of motion isoften desirable due to its efficiency and applicability in large deformationproblems. For example, in the large deformation case, the stress and materialhistory will move with the nodes and so no re-mapping of the quantities isrequired. Numerous formulations and types of meshless shape functions areavailable, here, we will focus on Galerkin formulations where natural neighbor(NN) interpolation is used; often referred to as the natural element method(NEM) [7]. Different nodal integration approaches have been developed suchas stabilized conforming nodal integration (SCNI) [3, 10, 5] and a least squaresstabilization approach [1]. The latter approach was applied using an EFGformulation and suggested the use of quadratic basis functions due to thesecond derivatives in the least squares term. Furthermore, the least squaresformulation in [1] is not directly applicable to natural neighbor schemes sincethey are not smooth at the nodes. Recently we have found that the SCNImethod applied to NN shape functions produces spurious low energy (notzero energy) modes in some problems and these modes did not appear tovanish with refinement. In this work, a modification to the SCNI method ismade that appears to provide stable results as demonstrated in eigenvalueand large deformation examples provided here.Natural neighbor (NN) shape functions have some distinct advantages overmany meshless schemes in that they interpolate the data making it simpler2 Michael Anthony Puso, Edward Zywicz, and J.S. Chento implement essential boundary conditions and the nodal adjacency onlyincludes near neighbors. The NN/NEM approach has been implemented in thecontext of a meshless method by exploiting the concept of alpha shapes [4, 6] todetermine/treat the free surface of cloud of points. In this work, the integrationmethod in [10, 5] is modified by applying an additional stabilization term. Inwhat follows: Section 2 introduces the global weak form, the NN interpolationscheme and the SCNI approach with added stabilization, Section 3 presentsresults demonstrating the necessity of the added stabilization and effectivenessof the proposed approach.2 FormulationThe formulation will be introduced in the context of linear elasticity and thestraightforward extensions to the nonlinear regime will be given at the end ofthe section.2.1 Galerkin MethodConsidering a body occupying the domain Ω, the governing equations of mo-tion are givenρu¨ = ∇σ + b (1)where u is the displacement field, b is the body force and σ = σ(ε) is theCauchy stress in terms of strainε(u) = 1/2(∇u+∇uT ) (2)Employing test function v, the weak form of (1) can then be written∫Ωρv · u¨ dΩ +∫Ωε(v) : σ(ε(v)) dΩ =∫Ωv · b dΩ +∫Γtv · t dΓ (3)where applied tractions t are specified on the boundary Γ t.2.2 Discrete FormFollowing standard procedure, the discrete displacement field is defineduh =N∑I=1φI(x)uI x ∈ Ωh (4)in terms of shape functions φI and nodal displacements uI over all nodesI = 1, N on the discretized domain Ωh. Unlike finite elements, the definitionof the domain Ωh in meshless methods is not straightforward and is definedA New Stabilized Nodal Integration Approach 3in what follows. Substituting (4) into the weak form (3) and applying nodalintegration, the following discrete form resultsN∑I=1mIvI · u¨I + VIεI(vh) : σ(εI(uh))− vI · f I = 0 (5)where mI is the nodal mass, fI is a collocated form of the applied loads, εIis the nodal strain and VI is the nodal volume. Since the NN shape functionsinterpolate the data, the nodal mass is simply mI = ρVI . The body forcecontribution to f I is computed similarly. Here, the traction force contributionto f I is computed using the surface mesh but other techniques are outlinedin [1]. The εI could be the symmetric gradient of the displacement at xI orcomputed using SCNI as shown in what follows. In this work, VI represents thevolume of the Voronoi cell ΩI about node I (Fig. 1(a)). The Voronoi diagramof a set of sites P := {x1,x2, . . .xN} is a partition such that all points xwithin a Voronoi cell ΩI are closer to the generating site I than to any othersite in P . Considering three dimensional space, this can be mathematicallystatedΩI = {x ∈ R3 : d(x,xI) < d(x,xJ) ∀ J 6= I} (6)To handle arbitrary non-convex geometries, the alpha − shapes approach isused to compute a surface mesh from the cloud of points. The Voronoi diagramis then ”clipped” by the surface in an approach similar to that in [4] and [6]thus producing the discretized domain Ωh. Figure 1(b) illustrates this processfor two dimensions and our implementation performs the equivalent in threedimensions. Of course, this approach may be too computationally intensivefor some applications and one may resort to simpler approaches to computenodal volumes (c.f. [1]).2.3 Natural Neighbor InterpolationThe Voronoi cell for point I is formed by polygonal (usually triangular) cellwalls and the points with common cell walls are the set of natural (or nearest)neighbors N (xI) for point I. The Laplace (or non-Sibson) form of NN inter-polation [8] computes interpolation functions at a point x in the followingway. Consider the Voronoi diagram shown in black in Fig. 2, the interpolationfunction φ at a point x is computed by doing a point insertion where x istreated as an additional vertex in the set P and a new Voronoi cell is formed(shown in blue in Fig. 2(a)). For the two dimensional case, the heights h andcell wall lengths s are used to compute the shape functions evaluated at pointxφJ(x) =sJ/hJ∑K∈N (x) sK/hK∀J ∈ N (x). (7)In three dimensions, the quantity s is a cell wall area. The NN interpolationfunctions are both linear exact and interpolative.4 Michael Anthony Puso, Edward Zywicz, and J.S. Chen 1 2 3 4 5 6 7 (a) 1 2 3 4 5 6 7 (b) Figure 1. (a) Voronoi diagram about cloud of points. (b) Surface mesh and clippedVoronoi cells. 1 2 3 4 5 6 7 x (a) 1 2 3 4 x h4 h3 h2 h1 s1 s4 s2 s3 (b) Figure 2. (a) Point insertion at x. (b) Geometry of Voronoi cell about x2.4 Stabilized Conforming Nodal Integration (SCNI)If smooth shape functions are used, one could merely take derivatives at thenodes to compute the nodal strain εI (although patch test satisfaction maybe violated). Since NN shape functions are not smooth, a different approach isrequired. Here, the nodal strain is computed using the SCNI approach whichis well described in [3, 10, 5]. In short, εI is computed from the volume averageof the strain ε over the Voronoi cell domain ΩI and using divergence theoremresults in the surface integralεI(uh) =1VI∫∂ΩI(uh ⊗ n+ n⊗ uh)/2 dΓ (8)A New Stabilized Nodal Integration Approach 5where n is the unit normal along the contour surface of cell I. Numericalintegration is employed in computing (8) such that one integration point atthe centroid of each Voronoi facet is used. This approach exactly satisfies patchtest and can be used without any additional stabilization in many situations asdemonstrated in [3, 10, 5]. On the other hand, it is shown here that problemsarise such that spurious low energy modes arise in some applications of SCNImotivating the following approach.2.5 Modified SCNIAs mentioned, the SCNI formulation produces spurious modes as seen inFig. 4(b). Although the SCNI approach does not have zero energy modes, astraightforward analysis can demonstrate that the method is not V1 coercive.Referring to (3) and (5), the discrete bilinear form for linear elasticity iswrittena(vh,uh) =N∑I=1εI(vh) : CεI(uh)VI (9)where C is the elasticity tensor and the discrete strain energy results whenvh.= uh. Considering a one dimensional BVP with evenly spaced nodes, theappropriate saw tooth displacement field would cause the nodal strain definedby (8) to be zero everywhere except at the boundaries. This would produce azero energy mode for an infinite domain and a spurious low energy mode fora finite domain. In order to eliminate these modes the Voronoi cells at nodeI are further decomposed into sub-cells ΩcI as shown for the two dimensionalcase in Fig. 3. Here each sub-cell is a triangle/polyhedron formed from thegenerating node and a Voronoi edge/facet in two/three dimensions respec-tively. In practice, we can merge adjacent triangle/polyhedron into reducedset of sub-cells. Examples in this work use eight sub-cells by performing thismerge process although some nodes may actually have less than eight sub-cells since they have less than eight Voronoi facets. The same SCNI approachis used to compute the strain over each sub-cell ΩcIεcI(uh) =1V cI∫∂ΩcI(uh ⊗ n+ n⊗ uh)/2 dΓ c = 1 : N cI (10)where N cI is the number of sub-cells used at node I and VcI is the volume ofthe sub-cell. The modified SCNI weak form is now proposeda(vh,uh) =N∑I=1[εI(vh) : CεI(uh)VI +NcI∑c=1αcI(εI(vh)− εcI(vh)) : C˜(εI(uh)− εcI(uh))V c] (11)6 Michael Anthony Puso, Edward Zywicz, and J.S. Chenwhere αcI is a penalty parameter and C˜ is the elasticity tensor or some modifiedversion. In the proceeding examples, the sub-cell penalty will be uniform i.e.αcI = 1 amounting to what is sometimes referred to as ”physical stabilization”[2, 9]. Here, the stiffness C˜ is chosen in the following way:Isotropic elastic material: Lame´ parameters µ and λ:µ˜ := µ λ˜ := max(λ, 25µ˜) (12)Isotropic plastic material: Lame´ parameter λ and plastic modulus ET :µ˜ := ET /2 λ˜ := max(λ, 25µ˜) (13)For non-isotropic materials, the choice for µ˜ and λ˜ would be less obviousbut still tractable. Using a maximum in (12) and (13) mitigates volumetriclocking. The nodal deformation gradients F I ,F cI are definedF I(uh) =1VI∫∂ΩIuh ⊗ n dΓ F cI(uh) =1V cI∫∂ΩcIuh ⊗ n dΓ e = 1 : N cI(14)and for large deformation kinematics, the internal virtual work (11) is rede-fineda(vh,uh) =N∑I=1[F I(vh) : P (F I(uh))VI +NcI∑e=1αcI(F I(vh)− F cI(vh)) : C˜(F I(uh)− F cI(uh))V c] (15)where P is the first Piola Kirchoff stress and VI and V cI are determined inthe initial reference configuration. The numerical implementation replaces theinternal virtual work (9) with (11) or (15) in the weak form of the equationsof motion (3). 1 2 3 4 6 7 Ω c Figure 3. Subcells in the Voronoi cell formed from xI .A New Stabilized Nodal Integration Approach 72.6 Implementation DetailsIn our implementation, the strain displacement matrices were computed atinitialization and then stored to avoid re-computation of the Voronoi diagramat every time step. This also applies to the nonlinear kinematics case wherethe total Lagrangian approach (15) was employed.The SCNI and the stabilization approach require the computation of shapefunctions at the Voronoi cell walls and polygons separating adjacent sub-cellsfor each node I. To expedite this calculation and to reduce connectivity, thesmallest set of nodes were employed in evaluating shape functions:1. Only the set N (xI) ∪ N (xJ) was used for Voronoi cell wallshape function evaluations where the cell wall was common tothe Voronoi cells of node xI and xJ .2. Only the set N (xI) was used for the interior sub-cell shape func-tion evaluations.The above almost always worked and one could probably even use less nodes.If the approach didn’t work (i.e. the cell wasn’t closed), more nodes wereadded to the set based on adjacency. Unlike changing the smoothing lengthof the weight functions with RKPM or EFG to get the appropriate amountof cover, adding nodes to the set to compute shape functions does not affectany previous computations.The results here used two dimensional triangular finite element interpola-tion everywhere on the exterior boundary surface mesh to do shape functionevaluations along the boundary. We have recently implemented two dimen-sional Laplace shape functions but have yet to compare the difference.3 ResultsThe following results demonstrate the necessity and effectiveness of the addi-tional stabilization. Furthermore, the proposed approach satisfies patch teststo machine precision and does not appear to lock in nearly incompressiblesituations. The constant pressure Q1P0 finite element is used for comparisonsince it performs reasonably well for the nearly incompressible case.EigenanalysisThis example illustrates the spurious modes using the standard SCNI ap-proach and the improved behavior using additional stabilization. An eigen-analysis of a nearly incompressible 1× 1× 1 block (E = 1, ν = 0.499, ρ = 1)was performed and the first eigenmodes using the Q1P0 brick element, theSCNI NN approach and the proposed approach are plotted in Fig. 4. The fre-quency versus eigenmode is plotted in Fig. 5 and it is seen that results fromthe new approach are nearly identical to that of the Q1P0 for the first 100modes. With SCNI, the frequencies for the first one hundred modes are lowdue to the presence of spurious modes.8 Michael Anthony Puso, Edward Zywicz, and J.S. Chen Figure 4. First eigenmode using: (a) brick elements, (b) SCNI (c) modified SCNI.0.20.30.40.50.60.70.80.910 20 40 60 80 100new stabilizationSCNIQ1P0frequency (Hz)eigenmode Figure 5. Frequency versus eigenmode for three approaches.Rubber BilletIn this example, the performance on a large deformation problem is evaluated.A rubber billet 50.8 in length and diameter (Fig. 6(a)) with Neo Hookean ma-terial parameters (E = 1, ν = 0.49) was compressed to 70% strain and onlythe bottom section was modeled due to symmetry. In the model, a unilateralcontact surface was used to vertically constrain the outer surface as it ex-pands (bulges) horizontally. The sequence of deformations from the proposedapproach are shown in Fig. 6(a-c) . The SCNI analysis terminated at 41%strain due to large oscillations (Fig. 6(d)) whereas the stress and deformationmatch nicely for the proposed approach and the Q1P0 hex element results atthis strain threshold (Fig. 6(b,e)). Oscillations do become apparent in the laststate shown in Fig. 6(c) compared to the Q1P0 results Fig. 6(e). The totalA New Stabilized Nodal Integration Approach 9force applied to the billet versus strain is shown plotted in Fig. 7 and it seenthat the proposed approach matches the Q1P0 quite well. (a) (d) (b) (e) (c) (f) Figure 6. Rubber billet: (a-c) sequence of deformation of stabilized NEM (d) un-stabilized NEM (e-f) deformed configurations of Q1P0. The vertical σz stress isindicated in (b-f).Taylor BarA standard benchmark for plasticity is the Taylor bar impact problem. Here,a copper bar impacts a rigid wall at high velocity. The following elastoplasticmaterial properties were used: E = 117 GPa, ν = 0.35, σy = 0.4 GPa, linearhardening modulus ET = 0.1 GPa and density ρ = 8930 kg/m3. The initialbar length was 32.4 mm, the initial radius was 3.2 mm and the initial velocitywas 227 m/s. The final deformed states are shown in Fig. 8 where the colorindicates displacement magnitude. Again, the new stabilization compares wellwith the Q1P0 finite element result but the SCNI shows oscillations.4 DiscussionA new stabilized nodal integration approach was developed and applied to thenatural neighbor meshless method. The approach performed well in eigenvalueand nonlinear benchmark problems. Here, a total Lagrangian approach wasused for large deformation. Extension of the method using an ”Eulerian” type10 Michael Anthony Puso, Edward Zywicz, and J.S. Chen-2000-1500-1000-50000 10 20 30 40 50 60 70new stabilizationSCNIQ1P0force% strain Figure 7. Force versus percent strain for rubber billet (a) (b) (c) Figure 8. Final configuration and plot of displacement magnitude from (a) newstabilized (b) SCNI and (c) Q1P0 approaches.A New Stabilized Nodal Integration Approach 11kernel for extremely large deformations occurring in penetration will be inves-tigated. The method appeared to perform well for the nearly incompressiblecase but would require modifications to treat the fully incompressible case.References1. S. Beissel and T. Belytschko, Nodal integration of the element-free Galerkinmethod, Computer Methods in Applied Mechanics and Engineering 139 (1996),49–74.2. T. Belytschko and Bindeman L.P., Assumed strain stabilization of the eightnode hexahedral, Computer Methods in Applied Mechanics and Engineering105 (1993), 225–260.3. J.S. Chen, C.T. Wu, S. Yoon, and Y. You, A stabilized conforming nodal in-tegration for Galerkin mesh-free methods, International Journal for NumericalMethods in Engineering 50 (2001), 435–466.4. E. Cueto, B. Calvo, and M. Doblare, Modeling three-dimensional piece-wise ho-mogenous domains using the α-shape based Natural Element Method, Interna-tional Journal for Numerical Methods in Engineering 54 (2002), 871–897.5. D. Gonzalez, E. Cueto, M.A. Martinez, and M. Doblare, Numerical integra-tion in Natural neighbor Galerkin methods, International Journal for NumericalMethods in Engineering 60 (2003), 2077–2104.6. S. Idelsohn, E. Onate, N. Calvo, and F. Del Pin, The meshless finite elementmethod, International Journal for Numerical Methods in Engineering 58 (2003),893–912.7. Sukumar N., Moran B., and Belytschko T., The Natural element method insolid mechanics, International Journal for Numerical Methods in Engineering43 (1998), 839–887.8. Sukumar N., Moran B., Yu Semenov Y., and Belikov V.V., Natural neighborGalerkin methods, International Journal for Numerical Methods in Engineering50 (1998), 1–27.9. M.A. Puso, A highly efficient enhanced assumed strain physically stabilized hex-ahedral element, International Journal for Numerical Methods in Engineering49 (2003), 1024–1064.10. J.W. Yoo, B. Moran, and J.S. Chen, Stabilized conforming nodal integration inthe Natural Element Method, International Journal for Numerical Methods inEngineering 60 (2004), 861–890.

A New Stabilized Nodal Integration Approach

https://digital.library.unt.edu/ark:/67531/metadc896093/m2/1/high_res_d/928188.pdf

A New Stabilized Nodal Integration Approach

Abstract

Similar works

Full text

Available Versions

UNT Digital Library