- Email: [email protected]
A stable, meshfree, nodal integration method for nearly incompressible solids
Finite Elements in Analysis and Design 51 (2012) 81–85
Contents lists available at SciVerse ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
A stable, meshfree, nodal integration method for nearly incompressible solids William Elmer c,n, J.S. Chen a, Mike Puso b, Ertugrul Taciroglu a a
Civil and Environment Engineering Department, University of California, Los Angeles, CA 90095, United States Methods Development Group, Lawrence Livermore National Laboratory, Livermore, CA 94550, United States c Defense Systems Analysis Group, Lawrence Livermore National Laboratory, 7000 East Ave., L-126, Livermore, CA 94550, United States b
a r t i c l e i n f o
abstract
Article history: Received 2 March 2010 Received in revised form 11 May 2011 Accepted 5 November 2011 Available online 23 December 2011
An improved nodal integration method for nearly incompressible materials is proposed. Although nodal integration methods can avoid displacement type locking with nearly incompressible materials, they often encounter pressure oscillations. A mixed approach to nodal integration is exploited to avoid the oscillation. The method is applied to nodal integration using meshless shape functions and appears to work well in simple two dimensional benchmark tests. The method currently relies on a structured discretization but generalizations are proposed as the basis for future work. & 2011 Elsevier B.V. All rights reserved.
Keywords: Pressure modes Nodal integration Incompressibility
1. Introduction Nodal integration is an attractive approach when implementing meshless methods because it reduces the number of stress point calculations and can be applied to large deformation problems in a natural way. However, issues of stability arise with meshfree methods that need to be addressed. The under-integration of nodally integrated meshfree methods has been dealt with in many different works [5,2,18,16]. In most of these studies, nearly incompressible materials were treated and locking was typically not encountered. On the other hand, pressure oscillations have been encountered with nodal integration of more constrained problems, as demonstrated, for example, in [17,1]. A simple constraint count has also been used to identify these issues in nearly incompressible problems [4]. As is well known, nodal integration methods do not satisfy the patch test in general; unless they satisfy the so-called ‘‘integration constraint’’ as demonstrated by Chen and co-workers in [5]. Chen’s method is general, but relies on Voronoi cells to satisfy the patch test. Convergent behavior still occurs for stabilized nodal integration when the nodal discretization is structured, as demonstrated in [16]. Improved behavior can be achieved by adding stress points [3,9–11], but the stress-point approach is less natural for meshless methods. Nodal integration of meshfree
n
Corresponding author. Tel.: þ 1 925 423 8403. E-mail address: [email protected] (W. Elmer).
0168-874X/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2011.11.001
methods is particularly attractive in very large deformation situations [7]. In the work to follow, the discretization is structured (e.g., nodes are placed on a rectilinear grid) and a linear elastic solution is adopted. A number of meshless treatments for nearly incompressible materials have been presented in the literature; including those by Dolbow and Belytschko [8] and Chen et al. [6]. Both of these works focus on implementations with background grids, and propose ways to reduce the ‘‘constraint count’’ associated with incompressibility. Extending [6] to nodal integration would be straightforward. A mixed method could be adopted that projects integration point stresses onto local polynomials in order to provide a more stable pressure. The approach taken in this work applies a mixed method to nodal integration; and no projection is made. Certain nodes are defined as ‘‘master’’ nodes where pressure is calculated directly from the strain at the node. These pressures are then interpolated over a local domain. Nodes deemed as ‘‘slaves’’ use only the interpolated pressure. The master nodes define the pressure over the domain and the resulting pressure interpolation is continuous. All nodes use a linear basis for the MLS kernel. The proposed method is similar in spirit to the composite element of Guo et al. [12]. These authors use the LBB-stable 10node tetrahedral element with continuous pressure interpolation as their starting point. Then, in lieu of using Lagrange (quadratic) shape functions for displacement fields, they subdivide the 10node tetrahedron into six linear tetrahedra. The pressure field is still linearly interpolated over the element domain and globally continuous. The pressure on the six linear sub-elements is determined by this linear pressure interpolation.
82
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
In what follows, the formulation of the pressure stabilization technique and its implementation for meshfree particles on a structured grid is presented first (Section 2). Extension of this method to more general problems and results from simple benchmark problems are provided in Sections 3 and 4, respectively. Finally, a discussion of the results and the conclusions are given in Section 5.
p2 ¼ 1=2ðp1 þp3 Þ, p6 ¼ 1=2ðp3 þp9 Þ,
2. Formulation The internal virtual work of elasticity for the nodal integration scheme is given by Z X r : de dO deA : rA V A þoptional stabilization terms ð1Þ O
A ¼ 1,n
where n is the number of nodes, rA and eA are the stress and strain at node A, and VA is the nodal volume. The stabilization terms may include some gradient based stabilization or additional stress points. Meshless methods typically use moving least squares (MLS) or the related reproducing kernel (RK) type shape functions to represent the displacement field. The simplest form of nodal integration is applied for this work. Here, MLS/RK shape functions are used to discretize the displacement field and the smoothed strain at the node A with coordinate xA is computed from the gradient of shape functions NB X eA ¼ rsym NB ðxA ÞuB ð2Þ B A SA
where uB is the displacement coefficient at node B, and S A is the connectivity at node A defined by the overlapping supports of the surrounding meshless shape functions. For appropriately structured sets of points, the definition for nodal strain (2) satisfies the patch test exactly due to symmetry. To treat the general case of scattered, unstructured nodes, the nodal strain defined by the SCNI method [5] can be used in place of (2) so as to satisfy the integration constraint, and thus, the patch test. The gradient based stabilization presented in [16] is used in this work. Here, an additional strain gradient energy stabilization is used in (1) such that Z X 1 X r : de dO deA : rA V A þ ½ð@x deA Þ : Cdev : ð@xi eA ÞV A 3 i ¼ 1;3 i O A ¼ 1,n ð3Þ where the deviatoric stiffness Cdev is used for the stabilization energy in place of the full stiffness to avoid volumetric locking for nearly incompressible behavior. Smoothed nodal strain eA , and the gradient over a small domain near the node, are used in the following developments. The use of (1) and (2) appears to avoid displacement locking due to incompressibility but can produce pressure oscillations in some problems as will be demonstrated later (cf., Fig. 5). The following simple approach is proposed to eliminate these oscillations: Given a nodal discretization, specific nodes are chosen as master nodes and the remaining are chosen as slave nodes. The pressure at master nodes is computed from the associated nodal volumetric strain as is standard, and the pressure at slave nodes is interpolated from the neighboring master nodes. This pressure interpolation is made using shape functions fM centered at P master nodes that satisfy partition of unity ð M fM ðxÞ ¼ 1Þ, as in X fM ðxS ÞpM ð4Þ pS ¼ M A SS
where S A fset of slave nodesg with support at node Sg.
The simplest possible case considers a square structured nodal discretization in two dimensions where there are nine nodes as shown in Fig. 1. Here nodes {1, 3, 7, 9} are the master nodes and {2, 4, 5, 6, 8} are the slaves. Stress and pressure at the master (corner) nodes are computed from strains at the corners. Given the appropriate symmetry and support sizes of the pressure interpolation functions fM , the pressure at the slaves boundary edge nodes {2,4,6,8} is computed as in
and
p4 ¼ 1=2ðp1 þ p7 Þ p8 ¼ 1=2ðp7 þ p9 Þ
ð5Þ
and pressure at interior node 5 is computed as p5 ¼ 1=4ðp1 þp3 þ p7 þ p9 Þ
ð6Þ
Following this structured pattern for larger sets of nodes, one would generate the discretization of master nodes and then place slave nodes between them so that all interior nodes have four master nodes, and all edge nodes have two masters. Minor deviations from rectilinear distributions may be handled by assigning more weight to the master node closer to a given slave, as would a partition of unity shape function. An extension to nodal distributions without a grid structured nodal neighborhood is proposed in Section 3. A B (B-bar) approach can be defined using the above method such that strain at node 5 is computed via
e 5ij ¼ e5ij 1=3½e5kk 1=4ðe1kk þ e3kk þ e7kk þ e9kk Þ
ð7Þ
The strain displacement matrix BA computes the strain at any node A in the domain via the spatial gradient operator in (2) such that
eA ¼ BA u
ð8Þ
where u is the vector of all nodal displacements. Similar procedures can be performed to calculate strains at other ‘‘side’’ of slave nodes. All points A are integration points, but different operations are performed on a portion of the strain displacement data. Assuming the structure of BA takes the following form for plane strain: 8 9 2 3 ðe Þ ðBA Þ1 > < A 11 > = 7 ðeA Þ22 ¼ 6 ð9Þ 4 ðBA Þ2 5u > : ðe Þ > ; ðBA Þ3 A 12 the B-bar strain displacement matrix at slave node S, for the very specific conditions outlined, is defined as 2 3 2 3 ðBS Þ1 þ ðBS Þ2 ðBM Þ1 þðBM Þ2 16 1 X 7 6 ð10Þ B S ¼ BS 4 ðBS Þ1 þ ðBS Þ2 5 þ f ðxS Þ4 ðBM Þ1 þðBM Þ2 7 5 3 3 MAS M T S 0 0T It should be noted that the connectivity of slave nodes is increased due to the influence of the master nodes. Because the proposed approach will likely be applied to an explicit scheme for large scale deformations and damage, an unsymmetric formulation
2
1
4
5
7
8
3
6
9
S S ¼ fset of master nodes Fig. 1. Meshfree template, master nodes 1, 3, 7 and 9.
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
is applied for the sake of efficiency. Stress r will be computed from B-bar strain e and force will be computed using the standard B. The implementation is ordered such that all master node stresses, strains and nodal forces are computed first. Slave node quantity calculations then follow. Here, master node volumetric strains saved from the first pass will be used to compute e S as in (7). Using these B-bar strains, slave node stress is computed, and the nodal forces are obtained using the (non-averaged) strain displacement B in lieu of B to minimize the amount of scatter. Here, the stiffness matrix for the unsymmetric formulation is written for slave nodes as in K S ¼ BTS CB S V S
Slave nodes that are placed on the facet will always be halfway in between two master nodes by construction. The weights from each node would be equal and available a priori to the averaging scheme. The Voronoi cells are then updated to recognize the new nodes. This procedure could potentially be done locally. The method would effectively produce a pressure distribution similar to the LBB-stable six-node quadratic displacement element with linear pressure distribution. Although quadratic bases are not used, stability is still achievable in the same way as the composite triangles proposed by Guo et al. in [14], which is shown with a dotted line connecting three nodes of Fig. 2.
ð11Þ
Another reason to use the unsymmetric formulation is satisfaction of the integration constraint [5]. This constraint now needs to be calculated with the averaged pressures, which nullifies the symmetry argument that previously held for smoothed nodal integration on a regular grid. Since the integration constraint fails to hold in a symmetric construction, the Petrov–Galerkin variational form is chosen and the averaged stresses are filtered by the regular strain, i.e., Z X 1 X r : de dO d eA : C : e A V A þ ð@x deA Þ : C : ð@xi e A ÞV A 3 i ¼ 1;3 i O A ¼ 1,n
4. Example problems 4.1. Cantilever beam with tip load A commonly used benchmark problem to investigate locking is the tip-loaded—and not necessarily slender—cantilever beam under flexure [14]. The analytical solution to this problem is due to Timoshenko [19], and can be compared to numerical
ð12Þ We note here that the same unsymmetric formulation is also used to compute the stabilization forces for expedience. The implementation of the stabilization term above exploits the same averaging scheme in Eq. (3). That is, the average smoothed strain gradients are computed at the master nodes and interpolated to the slave nodes. Also note that, only the deviatoric portion of the material stiffness Cdev was used for stabilization to avoid volumetric locking for the implementation in [16], as shown earlier in Eq. (3). This is not the case for the proposed B-bar averaging scheme where the full material stiffness is used in the stabilization energy of Eq. (12) without the effect of locking.
3. Generalization The techniques proposed so far can be generalized through the use of stabilized conforming nodal integration (SCNI). Extending the method to SCNI would allow the method to pass the patch test on grossly unstructured point distributions. For this transition to take place, the Delaunay tessellation would be performed and Voronoi cells constructed as shown in Fig. 2. The generalized method would then be applied by implanting slave nodes at the midway between nearest neighbors.
Fig. 2. Slave nodes would be placed on locations marked ‘‘ ’’.
83
Fig. 3. Energy norm versus mesh size. Slope 1.8.
Fig. 4. Pressure error norm versus mesh size. Slope 1.9.
84
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
Fig. 5. Pressure oscillation without stabilization (left). Gradient stabilization fails to control the oscillation (right).
Fig. 6. The finite element solution (left). Single point stresses are shown rendered. Elimination of the pressure oscillation using the proposed method (right). Slave nodes are indicated with blue ’s, masters with red J’s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
simulation results. Here, Poisson’s ratio is chosen to be n ¼ 0:4999 to represent an incompressible elastic material, and n ¼ 0:3 for compressibility. A parabolic distributed shear load is applied to the right side of the cantilever. The distributed forces of opposing shear force and bending moment resultants are applied to the left side. The nodal spacing in the axial and transverse directions may differ slightly, but the elements are not skewed. The method proposed here passes the linear patch test on a structured grid. The energy error norm and pressure error norm are shown in Figs. 3 and 4, respectively. The pressure (9n9l2 ) norm is calculated by integrating the absolute value of the difference between analytical and computational values at each node, weighted by site volume, multiplying by 1/2 and taking the square root. The error norm follows the same procedure, but on the norm of the difference between qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R computed and analytical stress ( 12 O ðsa sc Þ : ðsa sc Þ dO). Both norms display convergent behavior. The slopes of the energy and pressure error norms are approximately 1.8 and 1.9, respectively.1 The cantilever benchmark problem shows the pressure averaging scheme does nothing to disturb the convergence behavior of the method, as anticipated by analytical work on FEM [13].
1 This type of—‘‘better than optimal’’—convergence behavior is not unexpected from nodally integrated meshfree methods.
4.2. The incompressible punch problem Whereas the nodally integrated meshless method performs very well in some nearly incompressible problems, it does experience severe pressure oscillations for highly constrained problems. The problem considered here—namely, the incompressible punch problem—demonstrates the superior behavior of the mixed formulation described herein. The problem involves an elastic block that is highly constrained in three dimensions. The block has a thickness that is equal to 1/10th of its unit width/ height. The shape is constrained in the normal direction on the front and back, and on three sides. The left, right and bottom side set of nodes are all constrained from motion in the sides’ normal direction, but are free to move tangentially. The middle third of the top side is pressed down under displacement control, up to a total of 2% of the block’s height. The material is nearly incompressible with a Young’s modulus of 2 106 and Poisson’s ratio of 0.49999. An analytical solution is not available to this problem, and therefore, the meshfree results will be compared to those obtained using the linear, under-integrated finite element. The model has 13 nodes per side in the plane, and only two sheets of nodes in the out-of-plane direction. In this, and any other graded mesh configuration, the ellipsoidal supports [16] are able to tightly cover their neighboring nodes. The nodes on the boundary fully express the entire field’s value at that location, and
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
85
problem with an analytical solution and is found to preserve strong convergence. The method is also able to eliminate the pressure oscillation on another benchmark problem in which most boundary nodes are subject to constraint conditions. The method described here focused on structured point sets, extensions to the general case of unstructured point sets have been proposed, and will be the basis of future work.
Acknowledgments The research reported herein was supported by the National Science Foundation under Grant No. CMS-0547670, and in part, by funding from the Faculty Grants Program of UCLA Academic Senate Council on Research. This support is gratefully acknowledged. The opinions expressed in this paper are those of the authors and do not necessarily reflect those of the sponsors. References
Fig. 7. Improper displacements when stabilization is not used.
boundary conditions are directly imposed as the desired displacement values. Results using the stabilized non-conforming nodal integration approach, per Eq. (1), are shown in Fig. 5. The displacements in this, and in all subsequent figures, are scaled by a factor of 10. The subplots on the left and right show results without and with the gradient stabilization, respectively. Gradient stabilization reduces some oscillations, but is not adequate. The maximum pressures are an order of magnitude higher than those predicted by the finite element results (cf., subplot on the left of Fig. 6). Furthermore, highly negative pressures exist near the upper corners. The pressure oscillation remediation achieved by the proposed method is seen in Fig. 6 with the finite element solution for comparison. Pressures in the meshfree solution are reported at the nodes. The highest meshfree pressure is experienced by the nodes in the corners of the indented zone. By contrast, pressures reported by the finite element simulation are from the Gausspoints at the element center, and are smoothed for rendering.2 The pressure of the field at the center of the element may naturally be lower than at the nodes (7% difference observed). The fringe pattern of pressure agrees well with the finite element simulation. Stabilization is required to achieve correct displacements (cf. Fig. 7) when using the pressure oscillation remediation technique.
5. Conclusions The pressure oscillation encountered in a meshfree representation of a nearly incompressible continuum is eliminated while preserving nodal integration. The connectivity of the problem is increased, but the increase is not dramatic, and is limited for all problems. The method is tested on a common benchmark
2 The patch function in Matlab [15] is used to visualize the color data, so a skewed rendering may be visible as an artifact of the interpolation on a triangular facet.
[1] F.M. Andrade Pires, E.A. de Souza Neto, J.L. de la Cuesta Padilla, An assessment of the average nodal volume formulation for the analysis of nearly incompressible solids under finite strains, Commun. Numer. Methods Eng. 20 (7) (2004). [2] S. Beissel, T. Belytschko, Nodal integration of the element-free Galerkin method, Comput. Methods Appl. Mech. Eng. 139 (1–4) (1996) 49–74. [3] T. Belytschko, S. Xiao, Stability analysis of particle methods with corrected derivatives, Comput. Math. Appl. 43 (3–5) (2002) 329–350. [4] J.S. Chen, W. Hu, M.A. Puso, Y. Wu, X. Zhang, Strain smoothing for stabilization and regularization of Galerkin meshfree methods, in: Lecture Notes in Computational Science and Engineering, vol. 57, 2006, p. 57. [5] J.S. Chen, C.T. Wu, S. Yoon, Y. You, A stabilized conforming nodal integration for Galerkin mesh-free methods, Int. J. Numer. Methods Eng. 50 (2) (2001) 435–466. [6] J.S. Chen, S. Yoon, H.P. Wang, W.K. Liu, An improved reproducing kernel particle method for nearly incompressible finite elasticity, Comput. Methods Appl. Mech. Eng. 181 (1–3) (2000) 117–145. [7] J.S. Chen, S. Yoon, C.T. Wu, Nonlinear version of stabilized conforming nodal integration for Galerkin meshfree methods, Int. J. Numer. Methods Eng. 52 (2002) 2587–2615. [8] J. Dolbow, T. Belytschko, Volumetric locking in the element free Galerkin method, Int. J. Numer. Methods Eng. 46 (6) (1999). [9] Q. Duan, T. Belytschko, On the stabilization of stress-point integration in the element free Galerkin method, in: Meshfree Methods for Partial Differential Equations IV, 2008, p. 47. [10] Q. Duan, T. Belytschko, Gradient and dilatational stabilizations for stresspoint integration in the element-free Galerkin method, Int. J. Numer. Methods Eng. 77 (6) (2009). [11] T.P. Fries, T. Belytschko, Convergence and stabilization of stress-point integration in mesh-free and particle methods, Int. J. Numer. Methods Eng. 74 (7) (2008). [12] Y. Guo, M. Ortiz, T. Belytschko, E.A. Repetto, Triangular composite finite elements, Int. J. Numer. Methods Eng. 47 (2000) 287–316. [13] B.P. Lamichhane, From the Hu–Washizu formulation to the average nodal strain formulation, Comput. Methods Appl. Mech. Eng. 198 (49–52) (2009) 3957–3961. [14] W.K. Liu, Y. Guo, S. Tang, T. Belytschko, A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis, Comput. Methods Appl. Mech. Eng. 154 (1–2) (1998) 69–132. [15] MathWorks, MATLAB Release Notes. The MathWorks, Inc., Natick, MA, version 7.9 (r2009b) edition, 2009. [16] M.A. Puso, J.S. Chen, E. Zywicz, W. Elmer, Meshfree and finite element nodal integration methods, Int. J. Numer. Methods Eng. 74 (3) (2008) 416. [17] M.A. Puso, J. Solberg, A formulation and analysis of a stabilized nodally integrated tetrahedral, Int. J. Numer. Methods Eng. 67 (2006) 841–867. [18] T. Rabczuk, T. Belytschko, S.P. Xiao, Stable particle methods based on Lagrangian kernels, Comput. Methods Appl. Mech. Eng. 193 (12–14) (2004) 1035–1063. [19] S.P Timoshenko, J.N. Goodier, Theory of Elasticity, vol. 880, McGraw-Hill, 1990, p. 866.
Contents lists available at SciVerse ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
A stable, meshfree, nodal integration method for nearly incompressible solids William Elmer c,n, J.S. Chen a, Mike Puso b, Ertugrul Taciroglu a a
Civil and Environment Engineering Department, University of California, Los Angeles, CA 90095, United States Methods Development Group, Lawrence Livermore National Laboratory, Livermore, CA 94550, United States c Defense Systems Analysis Group, Lawrence Livermore National Laboratory, 7000 East Ave., L-126, Livermore, CA 94550, United States b
a r t i c l e i n f o
abstract
Article history: Received 2 March 2010 Received in revised form 11 May 2011 Accepted 5 November 2011 Available online 23 December 2011
An improved nodal integration method for nearly incompressible materials is proposed. Although nodal integration methods can avoid displacement type locking with nearly incompressible materials, they often encounter pressure oscillations. A mixed approach to nodal integration is exploited to avoid the oscillation. The method is applied to nodal integration using meshless shape functions and appears to work well in simple two dimensional benchmark tests. The method currently relies on a structured discretization but generalizations are proposed as the basis for future work. & 2011 Elsevier B.V. All rights reserved.
Keywords: Pressure modes Nodal integration Incompressibility
1. Introduction Nodal integration is an attractive approach when implementing meshless methods because it reduces the number of stress point calculations and can be applied to large deformation problems in a natural way. However, issues of stability arise with meshfree methods that need to be addressed. The under-integration of nodally integrated meshfree methods has been dealt with in many different works [5,2,18,16]. In most of these studies, nearly incompressible materials were treated and locking was typically not encountered. On the other hand, pressure oscillations have been encountered with nodal integration of more constrained problems, as demonstrated, for example, in [17,1]. A simple constraint count has also been used to identify these issues in nearly incompressible problems [4]. As is well known, nodal integration methods do not satisfy the patch test in general; unless they satisfy the so-called ‘‘integration constraint’’ as demonstrated by Chen and co-workers in [5]. Chen’s method is general, but relies on Voronoi cells to satisfy the patch test. Convergent behavior still occurs for stabilized nodal integration when the nodal discretization is structured, as demonstrated in [16]. Improved behavior can be achieved by adding stress points [3,9–11], but the stress-point approach is less natural for meshless methods. Nodal integration of meshfree
n
Corresponding author. Tel.: þ 1 925 423 8403. E-mail address: [email protected] (W. Elmer).
0168-874X/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2011.11.001
methods is particularly attractive in very large deformation situations [7]. In the work to follow, the discretization is structured (e.g., nodes are placed on a rectilinear grid) and a linear elastic solution is adopted. A number of meshless treatments for nearly incompressible materials have been presented in the literature; including those by Dolbow and Belytschko [8] and Chen et al. [6]. Both of these works focus on implementations with background grids, and propose ways to reduce the ‘‘constraint count’’ associated with incompressibility. Extending [6] to nodal integration would be straightforward. A mixed method could be adopted that projects integration point stresses onto local polynomials in order to provide a more stable pressure. The approach taken in this work applies a mixed method to nodal integration; and no projection is made. Certain nodes are defined as ‘‘master’’ nodes where pressure is calculated directly from the strain at the node. These pressures are then interpolated over a local domain. Nodes deemed as ‘‘slaves’’ use only the interpolated pressure. The master nodes define the pressure over the domain and the resulting pressure interpolation is continuous. All nodes use a linear basis for the MLS kernel. The proposed method is similar in spirit to the composite element of Guo et al. [12]. These authors use the LBB-stable 10node tetrahedral element with continuous pressure interpolation as their starting point. Then, in lieu of using Lagrange (quadratic) shape functions for displacement fields, they subdivide the 10node tetrahedron into six linear tetrahedra. The pressure field is still linearly interpolated over the element domain and globally continuous. The pressure on the six linear sub-elements is determined by this linear pressure interpolation.
82
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
In what follows, the formulation of the pressure stabilization technique and its implementation for meshfree particles on a structured grid is presented first (Section 2). Extension of this method to more general problems and results from simple benchmark problems are provided in Sections 3 and 4, respectively. Finally, a discussion of the results and the conclusions are given in Section 5.
p2 ¼ 1=2ðp1 þp3 Þ, p6 ¼ 1=2ðp3 þp9 Þ,
2. Formulation The internal virtual work of elasticity for the nodal integration scheme is given by Z X r : de dO deA : rA V A þoptional stabilization terms ð1Þ O
A ¼ 1,n
where n is the number of nodes, rA and eA are the stress and strain at node A, and VA is the nodal volume. The stabilization terms may include some gradient based stabilization or additional stress points. Meshless methods typically use moving least squares (MLS) or the related reproducing kernel (RK) type shape functions to represent the displacement field. The simplest form of nodal integration is applied for this work. Here, MLS/RK shape functions are used to discretize the displacement field and the smoothed strain at the node A with coordinate xA is computed from the gradient of shape functions NB X eA ¼ rsym NB ðxA ÞuB ð2Þ B A SA
where uB is the displacement coefficient at node B, and S A is the connectivity at node A defined by the overlapping supports of the surrounding meshless shape functions. For appropriately structured sets of points, the definition for nodal strain (2) satisfies the patch test exactly due to symmetry. To treat the general case of scattered, unstructured nodes, the nodal strain defined by the SCNI method [5] can be used in place of (2) so as to satisfy the integration constraint, and thus, the patch test. The gradient based stabilization presented in [16] is used in this work. Here, an additional strain gradient energy stabilization is used in (1) such that Z X 1 X r : de dO deA : rA V A þ ½ð@x deA Þ : Cdev : ð@xi eA ÞV A 3 i ¼ 1;3 i O A ¼ 1,n ð3Þ where the deviatoric stiffness Cdev is used for the stabilization energy in place of the full stiffness to avoid volumetric locking for nearly incompressible behavior. Smoothed nodal strain eA , and the gradient over a small domain near the node, are used in the following developments. The use of (1) and (2) appears to avoid displacement locking due to incompressibility but can produce pressure oscillations in some problems as will be demonstrated later (cf., Fig. 5). The following simple approach is proposed to eliminate these oscillations: Given a nodal discretization, specific nodes are chosen as master nodes and the remaining are chosen as slave nodes. The pressure at master nodes is computed from the associated nodal volumetric strain as is standard, and the pressure at slave nodes is interpolated from the neighboring master nodes. This pressure interpolation is made using shape functions fM centered at P master nodes that satisfy partition of unity ð M fM ðxÞ ¼ 1Þ, as in X fM ðxS ÞpM ð4Þ pS ¼ M A SS
where S A fset of slave nodesg with support at node Sg.
The simplest possible case considers a square structured nodal discretization in two dimensions where there are nine nodes as shown in Fig. 1. Here nodes {1, 3, 7, 9} are the master nodes and {2, 4, 5, 6, 8} are the slaves. Stress and pressure at the master (corner) nodes are computed from strains at the corners. Given the appropriate symmetry and support sizes of the pressure interpolation functions fM , the pressure at the slaves boundary edge nodes {2,4,6,8} is computed as in
and
p4 ¼ 1=2ðp1 þ p7 Þ p8 ¼ 1=2ðp7 þ p9 Þ
ð5Þ
and pressure at interior node 5 is computed as p5 ¼ 1=4ðp1 þp3 þ p7 þ p9 Þ
ð6Þ
Following this structured pattern for larger sets of nodes, one would generate the discretization of master nodes and then place slave nodes between them so that all interior nodes have four master nodes, and all edge nodes have two masters. Minor deviations from rectilinear distributions may be handled by assigning more weight to the master node closer to a given slave, as would a partition of unity shape function. An extension to nodal distributions without a grid structured nodal neighborhood is proposed in Section 3. A B (B-bar) approach can be defined using the above method such that strain at node 5 is computed via
e 5ij ¼ e5ij 1=3½e5kk 1=4ðe1kk þ e3kk þ e7kk þ e9kk Þ
ð7Þ
The strain displacement matrix BA computes the strain at any node A in the domain via the spatial gradient operator in (2) such that
eA ¼ BA u
ð8Þ
where u is the vector of all nodal displacements. Similar procedures can be performed to calculate strains at other ‘‘side’’ of slave nodes. All points A are integration points, but different operations are performed on a portion of the strain displacement data. Assuming the structure of BA takes the following form for plane strain: 8 9 2 3 ðe Þ ðBA Þ1 > < A 11 > = 7 ðeA Þ22 ¼ 6 ð9Þ 4 ðBA Þ2 5u > : ðe Þ > ; ðBA Þ3 A 12 the B-bar strain displacement matrix at slave node S, for the very specific conditions outlined, is defined as 2 3 2 3 ðBS Þ1 þ ðBS Þ2 ðBM Þ1 þðBM Þ2 16 1 X 7 6 ð10Þ B S ¼ BS 4 ðBS Þ1 þ ðBS Þ2 5 þ f ðxS Þ4 ðBM Þ1 þðBM Þ2 7 5 3 3 MAS M T S 0 0T It should be noted that the connectivity of slave nodes is increased due to the influence of the master nodes. Because the proposed approach will likely be applied to an explicit scheme for large scale deformations and damage, an unsymmetric formulation
2
1
4
5
7
8
3
6
9
S S ¼ fset of master nodes Fig. 1. Meshfree template, master nodes 1, 3, 7 and 9.
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
is applied for the sake of efficiency. Stress r will be computed from B-bar strain e and force will be computed using the standard B. The implementation is ordered such that all master node stresses, strains and nodal forces are computed first. Slave node quantity calculations then follow. Here, master node volumetric strains saved from the first pass will be used to compute e S as in (7). Using these B-bar strains, slave node stress is computed, and the nodal forces are obtained using the (non-averaged) strain displacement B in lieu of B to minimize the amount of scatter. Here, the stiffness matrix for the unsymmetric formulation is written for slave nodes as in K S ¼ BTS CB S V S
Slave nodes that are placed on the facet will always be halfway in between two master nodes by construction. The weights from each node would be equal and available a priori to the averaging scheme. The Voronoi cells are then updated to recognize the new nodes. This procedure could potentially be done locally. The method would effectively produce a pressure distribution similar to the LBB-stable six-node quadratic displacement element with linear pressure distribution. Although quadratic bases are not used, stability is still achievable in the same way as the composite triangles proposed by Guo et al. in [14], which is shown with a dotted line connecting three nodes of Fig. 2.
ð11Þ
Another reason to use the unsymmetric formulation is satisfaction of the integration constraint [5]. This constraint now needs to be calculated with the averaged pressures, which nullifies the symmetry argument that previously held for smoothed nodal integration on a regular grid. Since the integration constraint fails to hold in a symmetric construction, the Petrov–Galerkin variational form is chosen and the averaged stresses are filtered by the regular strain, i.e., Z X 1 X r : de dO d eA : C : e A V A þ ð@x deA Þ : C : ð@xi e A ÞV A 3 i ¼ 1;3 i O A ¼ 1,n
4. Example problems 4.1. Cantilever beam with tip load A commonly used benchmark problem to investigate locking is the tip-loaded—and not necessarily slender—cantilever beam under flexure [14]. The analytical solution to this problem is due to Timoshenko [19], and can be compared to numerical
ð12Þ We note here that the same unsymmetric formulation is also used to compute the stabilization forces for expedience. The implementation of the stabilization term above exploits the same averaging scheme in Eq. (3). That is, the average smoothed strain gradients are computed at the master nodes and interpolated to the slave nodes. Also note that, only the deviatoric portion of the material stiffness Cdev was used for stabilization to avoid volumetric locking for the implementation in [16], as shown earlier in Eq. (3). This is not the case for the proposed B-bar averaging scheme where the full material stiffness is used in the stabilization energy of Eq. (12) without the effect of locking.
3. Generalization The techniques proposed so far can be generalized through the use of stabilized conforming nodal integration (SCNI). Extending the method to SCNI would allow the method to pass the patch test on grossly unstructured point distributions. For this transition to take place, the Delaunay tessellation would be performed and Voronoi cells constructed as shown in Fig. 2. The generalized method would then be applied by implanting slave nodes at the midway between nearest neighbors.
Fig. 2. Slave nodes would be placed on locations marked ‘‘ ’’.
83
Fig. 3. Energy norm versus mesh size. Slope 1.8.
Fig. 4. Pressure error norm versus mesh size. Slope 1.9.
84
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
Fig. 5. Pressure oscillation without stabilization (left). Gradient stabilization fails to control the oscillation (right).
Fig. 6. The finite element solution (left). Single point stresses are shown rendered. Elimination of the pressure oscillation using the proposed method (right). Slave nodes are indicated with blue ’s, masters with red J’s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
simulation results. Here, Poisson’s ratio is chosen to be n ¼ 0:4999 to represent an incompressible elastic material, and n ¼ 0:3 for compressibility. A parabolic distributed shear load is applied to the right side of the cantilever. The distributed forces of opposing shear force and bending moment resultants are applied to the left side. The nodal spacing in the axial and transverse directions may differ slightly, but the elements are not skewed. The method proposed here passes the linear patch test on a structured grid. The energy error norm and pressure error norm are shown in Figs. 3 and 4, respectively. The pressure (9n9l2 ) norm is calculated by integrating the absolute value of the difference between analytical and computational values at each node, weighted by site volume, multiplying by 1/2 and taking the square root. The error norm follows the same procedure, but on the norm of the difference between qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R computed and analytical stress ( 12 O ðsa sc Þ : ðsa sc Þ dO). Both norms display convergent behavior. The slopes of the energy and pressure error norms are approximately 1.8 and 1.9, respectively.1 The cantilever benchmark problem shows the pressure averaging scheme does nothing to disturb the convergence behavior of the method, as anticipated by analytical work on FEM [13].
1 This type of—‘‘better than optimal’’—convergence behavior is not unexpected from nodally integrated meshfree methods.
4.2. The incompressible punch problem Whereas the nodally integrated meshless method performs very well in some nearly incompressible problems, it does experience severe pressure oscillations for highly constrained problems. The problem considered here—namely, the incompressible punch problem—demonstrates the superior behavior of the mixed formulation described herein. The problem involves an elastic block that is highly constrained in three dimensions. The block has a thickness that is equal to 1/10th of its unit width/ height. The shape is constrained in the normal direction on the front and back, and on three sides. The left, right and bottom side set of nodes are all constrained from motion in the sides’ normal direction, but are free to move tangentially. The middle third of the top side is pressed down under displacement control, up to a total of 2% of the block’s height. The material is nearly incompressible with a Young’s modulus of 2 106 and Poisson’s ratio of 0.49999. An analytical solution is not available to this problem, and therefore, the meshfree results will be compared to those obtained using the linear, under-integrated finite element. The model has 13 nodes per side in the plane, and only two sheets of nodes in the out-of-plane direction. In this, and any other graded mesh configuration, the ellipsoidal supports [16] are able to tightly cover their neighboring nodes. The nodes on the boundary fully express the entire field’s value at that location, and
W. Elmer et al. / Finite Elements in Analysis and Design 51 (2012) 81–85
85
problem with an analytical solution and is found to preserve strong convergence. The method is also able to eliminate the pressure oscillation on another benchmark problem in which most boundary nodes are subject to constraint conditions. The method described here focused on structured point sets, extensions to the general case of unstructured point sets have been proposed, and will be the basis of future work.
Acknowledgments The research reported herein was supported by the National Science Foundation under Grant No. CMS-0547670, and in part, by funding from the Faculty Grants Program of UCLA Academic Senate Council on Research. This support is gratefully acknowledged. The opinions expressed in this paper are those of the authors and do not necessarily reflect those of the sponsors. References
Fig. 7. Improper displacements when stabilization is not used.
boundary conditions are directly imposed as the desired displacement values. Results using the stabilized non-conforming nodal integration approach, per Eq. (1), are shown in Fig. 5. The displacements in this, and in all subsequent figures, are scaled by a factor of 10. The subplots on the left and right show results without and with the gradient stabilization, respectively. Gradient stabilization reduces some oscillations, but is not adequate. The maximum pressures are an order of magnitude higher than those predicted by the finite element results (cf., subplot on the left of Fig. 6). Furthermore, highly negative pressures exist near the upper corners. The pressure oscillation remediation achieved by the proposed method is seen in Fig. 6 with the finite element solution for comparison. Pressures in the meshfree solution are reported at the nodes. The highest meshfree pressure is experienced by the nodes in the corners of the indented zone. By contrast, pressures reported by the finite element simulation are from the Gausspoints at the element center, and are smoothed for rendering.2 The pressure of the field at the center of the element may naturally be lower than at the nodes (7% difference observed). The fringe pattern of pressure agrees well with the finite element simulation. Stabilization is required to achieve correct displacements (cf. Fig. 7) when using the pressure oscillation remediation technique.
5. Conclusions The pressure oscillation encountered in a meshfree representation of a nearly incompressible continuum is eliminated while preserving nodal integration. The connectivity of the problem is increased, but the increase is not dramatic, and is limited for all problems. The method is tested on a common benchmark
2 The patch function in Matlab [15] is used to visualize the color data, so a skewed rendering may be visible as an artifact of the interpolation on a triangular facet.
[1] F.M. Andrade Pires, E.A. de Souza Neto, J.L. de la Cuesta Padilla, An assessment of the average nodal volume formulation for the analysis of nearly incompressible solids under finite strains, Commun. Numer. Methods Eng. 20 (7) (2004). [2] S. Beissel, T. Belytschko, Nodal integration of the element-free Galerkin method, Comput. Methods Appl. Mech. Eng. 139 (1–4) (1996) 49–74. [3] T. Belytschko, S. Xiao, Stability analysis of particle methods with corrected derivatives, Comput. Math. Appl. 43 (3–5) (2002) 329–350. [4] J.S. Chen, W. Hu, M.A. Puso, Y. Wu, X. Zhang, Strain smoothing for stabilization and regularization of Galerkin meshfree methods, in: Lecture Notes in Computational Science and Engineering, vol. 57, 2006, p. 57. [5] J.S. Chen, C.T. Wu, S. Yoon, Y. You, A stabilized conforming nodal integration for Galerkin mesh-free methods, Int. J. Numer. Methods Eng. 50 (2) (2001) 435–466. [6] J.S. Chen, S. Yoon, H.P. Wang, W.K. Liu, An improved reproducing kernel particle method for nearly incompressible finite elasticity, Comput. Methods Appl. Mech. Eng. 181 (1–3) (2000) 117–145. [7] J.S. Chen, S. Yoon, C.T. Wu, Nonlinear version of stabilized conforming nodal integration for Galerkin meshfree methods, Int. J. Numer. Methods Eng. 52 (2002) 2587–2615. [8] J. Dolbow, T. Belytschko, Volumetric locking in the element free Galerkin method, Int. J. Numer. Methods Eng. 46 (6) (1999). [9] Q. Duan, T. Belytschko, On the stabilization of stress-point integration in the element free Galerkin method, in: Meshfree Methods for Partial Differential Equations IV, 2008, p. 47. [10] Q. Duan, T. Belytschko, Gradient and dilatational stabilizations for stresspoint integration in the element-free Galerkin method, Int. J. Numer. Methods Eng. 77 (6) (2009). [11] T.P. Fries, T. Belytschko, Convergence and stabilization of stress-point integration in mesh-free and particle methods, Int. J. Numer. Methods Eng. 74 (7) (2008). [12] Y. Guo, M. Ortiz, T. Belytschko, E.A. Repetto, Triangular composite finite elements, Int. J. Numer. Methods Eng. 47 (2000) 287–316. [13] B.P. Lamichhane, From the Hu–Washizu formulation to the average nodal strain formulation, Comput. Methods Appl. Mech. Eng. 198 (49–52) (2009) 3957–3961. [14] W.K. Liu, Y. Guo, S. Tang, T. Belytschko, A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis, Comput. Methods Appl. Mech. Eng. 154 (1–2) (1998) 69–132. [15] MathWorks, MATLAB Release Notes. The MathWorks, Inc., Natick, MA, version 7.9 (r2009b) edition, 2009. [16] M.A. Puso, J.S. Chen, E. Zywicz, W. Elmer, Meshfree and finite element nodal integration methods, Int. J. Numer. Methods Eng. 74 (3) (2008) 416. [17] M.A. Puso, J. Solberg, A formulation and analysis of a stabilized nodally integrated tetrahedral, Int. J. Numer. Methods Eng. 67 (2006) 841–867. [18] T. Rabczuk, T. Belytschko, S.P. Xiao, Stable particle methods based on Lagrangian kernels, Comput. Methods Appl. Mech. Eng. 193 (12–14) (2004) 1035–1063. [19] S.P Timoshenko, J.N. Goodier, Theory of Elasticity, vol. 880, McGraw-Hill, 1990, p. 866.