TheMeshless Local Petrov-GalerkinMethod for Large Deformation Analysis of Hyperelastic Materials

Nonlinear formulations of the meshless local Petrov-Galerkin method (MLPG) are presented for the large deformation analysis of hyperelastic materials which are considered to be incompressible or nearly incompressible. The MLPG method requires no explicit mesh in computation and therefore avoids mesh distortion difficulties. In this paper, a simple Heaviside test function is chosen for reducing the computational effort by simplifying domain integrals for hyperelasticity problems. Trial functions are constructed using the radial basis function (RBF) coupled with a polynomial basis function. The plane stress hypothesis and a pressure projection method are employed to overcome the incompressibility or nearly incompressibility in the plane stress and plane strain problems, respectively. Effects of the sizes of local subdomain and interpolation domain on the performance of the present MLPG method are investigated. The behaviour of shape parameters of multiquadrics (MQ) function has been studied. Numerical results for several examples show that the present method is effective in dealing with large deformation hyperelastic materials problems.


Introduction
An important class of materials is the hyperelastic materials in engineering.These materials are commonly subjected to large elastic deformations.Due to the highly nonlinear nature of the governing differential equations, analytical solutions exist for only a few problems involving simple geometries and constitutive laws [1].Then a number of numerical schemes have been used to solve this problem.
Despite the finite element method (FEM) is now a standard analysis tool for the finite deformation problem.The application of FEM to large strain problems still has some important drawbacks that cannot be overcome, due to extremely large deformations and the incompressible or nearly incompressible nature of hyperelastic materials.Many FEMs have been developed to handle the volumetric locking and pressure oscillation difficulty resulting from the incompressibility or nearly incompressibility constraint, but the excessive mesh distortion is still reducing the convergence ratio and the overall accuracy of solution [2].Normally this requires remeshing steps and the corresponding projection of the current results in FEM, thus in this case the standard FEM formulations become ineffective or inaccurate without additional improvements [3].Recently, the analysis of hyperelastic materials also has been presented in the boundary element method (BEM) [4,5].However, there exist additional problem of singular integrals occurring in nonlinear boundary integral equations, it is required to regularize the different types of singularities.In addition, the excessive mesh distortion cannot be overcome in the BEM, because the BEM establishing equation is based on the element like FEM.
Another possibility is meshless methods, which have attracted considerable interest over the past decade.Their main advantage is avoiding the use of a mesh due to the development of new shape functions that allow the interpolation of field variables to be accomplished at a global level or a local level, which makes them very adequate for large strain problems.Some of them are Smooth Particle Hydrodynamics (SPH) method [6], Diffuse Element Method (DEM) [7], Element Free Galerkin (EFG) method [8], hpcloud Method [9], Reproducing Kernel Particle Method (RKPM) [10], Natural Element Method (NEM) [11], Meshless Galerkin Method using Radial Basis Function (RBF) ISRN Mechanical Engineering [12], and Meshless Local Petrov-Glerkin Method (MLPG) [13].Most of these methods with the exception of MLPG are not really meshless methods, since they use the background mesh for the numerical integration of the weak form.
The MLPG method is a truly meshless method that computed over a local subdomain based on a weak form.It offers a lot of flexibility to deal with large deformation problems.Remarkable successes of the MLPG method have been reported in solving the potential problems, the convection-diffusion problems and the nonlinear boundary value problems by Atluri et al. [13][14][15]; the fracture mechanics problems by Kim and Atluri [16]; the Navier-Stokes flows by Lin and Atluri [17]; the elasticity problems and plate bending problems by Long [18,19]; static and free vibration analysis of thin plates by Gu and Liu [20]; Crack Tip Fields problems by Ching and Batra [21,22]; anisotropic elasticity problems and crack analysis in 3-D axisymmetric FGM bodies by Sladek et al. [23,24].However, there exist some inconveniences in the MLPG method when the shape functions are obtained from interpolation schemes, such as the Moving Least Squares method (MLS), Partition of Unity Method (PUM), Reproducing Kernel Particle Method (RKPM), or Shepard function.These shape functions lack the Kronecker Delta property for the trialfunction interpolation, and special treatments have to be introduced to impose essential boundary conditions.Most recently, a Heaviside test function and radial basis function coupled with a polynomial basis function are used for developing the MLPG method [25][26][27], which make the MLPG method more flexible in dealing with large deformation problems.
In this paper, trial functions are constructed using a radial basis function coupled with a polynomial basis function.So no additional treatments are required for imposing essential boundary conditions.The Heaviside test function is chosen for eliminating or simplifying domain integrals in the global stiffness matrix.The plane stress hypothesis and a pressure projection method [28][29][30] are developed to overcome the incompressibility or nearly incompressibility in the plane stress and plane strain hyperelastic materials problems, respectively.Further, the effects of the sizes of the subdomain and interpolation domain on the performance of the present method with different shape parameters are illustrated in the large deformation problems.

Governing Equations
Consider a body which initially occupies a region Ω X with boundary Γ X .The deformation of a material particle X ∈ Ω X at time t is described by x(X, t), and the displacement of the particle X is defined by u(X, t) = x(X, t) − X. ( We consider the reference frame to be a rectangular Cartesian system.The deformation gradient F, the Green-Lagrangian strain E, the Green deformation tensor G, and the second Piola-Kirchhoff stress S are where W is the strain energy density function of the hyperelastic materials.Except for the plane stress problems, the incompressibility of hyperelasticity leads to severe numerical difficulties [3].Hence, the nearly incompressible hyperelastisity with a strain energy contribution involving the bulk modulus K is considered.The incompressibility can be recovered by letting K tend to infinity.In the present analysis, the modified Mooney-Rivlin strain energy density function W is used where W(I 1 , I 2 ) and W(J) are the deviatoric and volumetric terms of the strain energy density function, respectively.A 10 and A 01 are material constants, I 1 and I 2 are the reduced invariants, given by I 1 , I 2 , I 3 are the first, second, and third invariants of the Green deformation tensor, respectively, and J = I 1/2 3 .The second Piola-Kirchhoff stress associated with the modified Mooney-Rivlin strain energy density function defined in (5) is where the hydrostatic pressure P is related to W(J) by The incremental stress-strain relation can be expressed as where C i jkl is the fourth-order material response tensor and g i j is the second-order material response tensor, obtained by ISRN Mechanical Engineering 3 The explicit expressions of C i jkl , C i jkl , and g i j are given as where G −1 i j is the i jth component of G −1 not 1/G i j .

Interpolation Approximation
Consider a domain Ω x , the neighborhood of point x, which is located within the problem domain Ω.To approximate the distribution of function u in Ω x over a number of randomly located nodes x i , i = 1, 2, . . ., N, the approximate function u h (x) of u, for all x ∈ Ω x , can be defined by [26] where a i are coefficients for the radial basis function R i (x) and b j are coefficients for the polynomial basis p j (x).The number of the radial basis n is determined by the number of nodes in the support domain Ω x , and the number of the polynomial basis m can be chosen based on the reproduction requirement.A minimum number of terms of the polynomial basis and more terms of the radial basis (m < n) are adopted for better stability.The vectors in (13) are defined as The radial basis function R i (x) has many kinds of forms and the multiquadrics (MQ) function is given by where c, q are shape parameters.
The polynomial basis has the following monomial terms in the two-dimensional problem The coefficients a i and b j in ( 13) are determined by enforcing that the interpolation function passes through all nodes within the support domain.The interpolation at the kth point has the form which cannot be solved because n equations and n + m unknown coefficients are involved, so constraint conditions have to be introduced as the from Equations ( 17) and ( 18) can be expressed in matrix form as follows: where Suppose A −1 exists, ( 19) can be rewritten in the following form: Using ( 21), we have where P T m R −1 Q P m is termed as a transformed moment matrix.Finally, ( 13) is expressed as (23) where the shape functions Φ(x) can be written as follows: in which the shape function ϕ k (x) for the kth node is given by where S a ik is the (i, k) element of matrix S a and S b jk is the ( j, k) element of matrix S b .The derivatives of ϕ k (x)can be obtained as follows: The partial derivatives of MQ function can be obtained as follows:

The MLPG Formulation of Hyperelasticity
Consider a hyperelasitc body initially occupying a domain 0 Ω with boundary 0 Γ.Under the action of external forces, the hyperelastic body deforms into a configuration with domain t Ω, essential boundary t Γ u , and natural boundary t Γ t .During the deformation process, a material point X ∈ 0 Ω is moved to x = X + u ∈ t Ω.For this problem, the static equilibrium equation is described as where σ i j is the Cauchy stress and b i is the body force.The boundary conditions are given as where u i and t i are the prescribed displacements and tractions on boundary t Γ u and t Γ t , respectively.n j is the unit outward normal to the boundary t Γ.
A generalized local weak form of ( 28) and ( 29) over a local subdomain t Ω s can be written as follows: where u i and v i are the trial and test function, respectively.Using and divergence theorem in (30) leads to where t i = σ i j n j and n j is a unit outward normal to the boundary It should be mentioned that (32) holds regardless of the size and the shape of t Ω s , so the shape of subdomain t Ω s can be taken to be a circle in 2D problems without losing generality.
Applying the natural boundary condition, where In order to simplify (33), test functions v i are chosen such that they eliminate or simplify the domain integral on t Ω s .This can be accomplished by using the Heaviside step function where ζ is an arbitrary constant (ζ = 1 is used in this study).Using this choice, the partial derivatives of the test function v i, j are identically zero; hence the corresponding domain integral involved in (33) is identically zero.Using (34) and rearranging (33), we obtain the following local weak form: Due to the current configuration being unknown in (35), the initial configuration is selected as the reference configuration in calculation.The second Piola-Kirchhoff stress and Green-Lagrangian strain are used as stress and strain measures, respectively.The external tractions can be expressed as where the boundary 0 Γ is referred to the initial configuration and N j is the initial unit outward normal to the boundary 0 Γ.Then (35) can be rewritten as where t 0i and b 0i are the prescribed traction and body force measured in the initial configuration, respectively.
To handle geometric and material nonlinearity, we consider an incremental equation.The incremental decomposable terms are expressed as follows: Substituting (38) into (37) leads to Equation ( 39) is a nonlinear integral equation that exist; quadratic and cubic terms of Δu i after ΔS i j are denoted by ΔE i j and ΔP using (10) and the Green-Lagrangian strain increment Then linearization must be carried out for solving the integral equation.The linearization of (39) is given by Many displacement-based methods derived from single field variational principle failed in solving incompressible and nearly incompressibility problems because the volumetric locking and pressure oscillation were encountered [28].
In Section 2, the strain energy density was separated into pure deviatoric and pure volumetric terms.Consequently, the hydrostatic pressure is only related to the volumetric energy.Here the condition of plane stress hypothesis and a pressure projection method are used to overcome the incompressibility or nearly incompressibility in the plane stress and plane strain problems, respectively.For the states of plane stress, the incompressibility condition does not impose any significant difficulties.Working in the principal directions and using the modified Mooney-Rivlin strain energy density function, via the plane stress hypothesis σ 3 = 0, we can obtain [3] where λ 1 , λ 2 , and λ 3 are the principal stretch ratios with respect to both the coordinate system and the strain measure.
In the case of the plane stress, there exists Then Using (41), the hydrostatic pressure increment is given by where

ISRN Mechanical Engineering
For the states of plane strain, we reduce the independent constraint equations in each subdomain 0 Ω s by projecting pressure (9) onto a lower-order space using a pressure projection method proposed by Chen and Pan [28].The lower-order approximation of pressure is performed locally by introducing a least squares method, that is, where a constant field is used as the lower-order space, and is the L 2 norm in the subdomain 0 Ω s .The minimization of Θ(P) leads to Similarly, in each subdomain 0 Ω s , we project the hydrostatic pressure increment onto a constant field by the least squares approximation The minimization of Θ(ΔP) leads to where A s is the area of subdomain 0 Ω s .B M is the gradient matrix of Green-Lagrangian strain.If the most severe overconstrained condition is encountered, one-point integration to calculate (48) and (50) can be employed.

Discretizations
Using the radial basis function coupled with polynomial basis function to construct the trial function, Substituting (51) into (41) leads to the discrete iterative equation where n and v are load steps and iteration counters, respectively.It should be noted that the essential boundary condition can be implemented easily because the shape functions constructed by the radial basis function coupled with polynomial basis function in Section 3 possess the Kronecker Delta function property.Then the stiffness matrix K and the load vector f are defined as The term K * * IJ results from the hydrostatic pressure increment and can be expressed as for the plane strain. (54) For the plane problem, the explicit expressions of vectors and matrixes in the above equations are defined by From ( 53) and ( 54), it can be found that the domain integral over 0 Ω s is altogether avoided except for the terms resulting from projection of the hydrostatic pressure increment in the plane strain problems, which will greatly improve the effectiveness of the present method, especially for solving the large deformation problems.

Numerical Examples
Several numerical examples for hyperelastic material problems are presented to demonstrate the accuracy and efficiency of the present meshless method.In this paper, the modified Mooney-Rivlin strain energy density function is used.In the interpolation approximation, the polynomial basis function and MQ function are employed.Unless specified, otherwise, in computation, 8 and 6 Gauss integration points are used on each section of 0 L s and 0 Γ s for the numerical integration, respectively.Radii of the influence domain r Q is set to α Q d 3 i , d 3 i is the distance to the third nearest neighboring point from node i. Radii of the subdomain r s is set to α s d 1 i , d 1 i is the distance to the nearest neighboring point from node i. α Q and α s are the scaling coefficients for determining the size of interpolation domain and subdomain, respectively.

Uniaxial Tension-Compression
Problem.An 8.0 × 2.0 plane stress rubber block is subjected to tension and compression along the axial direction, with no lateral restraints on two ends.Since the stress-strain relation of this problem is independent of cross-sectional geometry, the analytical solution can be found in reference [30].For the purpose of error estimation, the relative error is defined as where σ MLPG and σ Exat are the axial Cauchy stress computed by the present method and the analytical method, respectively.
The problem domain is discretized as shown in Figure 1 with 85(17 × 5) regular and 90 irregular nodes.The motion of the rubber block is stretched up to 800 percent and compressed down to 20 percent in axial direction, respectively.The material constants A 10 = 0.2599 and A 01 = 0.1608 are used.The scaling coefficients are adopted as α s = 0.75 and α Q = 3.5.
Figures 2 and 3 give the variation of relative errors with shape parameters for MQ function, while the rubber block is stretched up to 800 percent in axial direction.From Figure 2, clearly both 1.03 and 1.99 are optimal shape parameters for q when the polynomial term is not included (m = 0).Figure 3 shows that the present method obtained results have good accuracy with C ranging from 0.5 to 3 and q = −0.5, 0.5, 1.03, 1.99 when a linear polynomial basis (m = 3).Effect of irregular distributed nodes and the polynomial term on axial Cauchy stress are shown in Figure 4 where the shape parameters of MQ function C = 1.42 and q = 1.03 are used.Figure 4 shows that the present method results calculated using regular and irregular with linear and quadratic polynomial basis have excellent agreement with analytical solutions while the rubber block is compressed.In Figure 4, the λ is defined as the deformation ratio that the current length divided by initial length of the rubber block.

Inflation of a Long Rubber
Tube.An infinitely long rubber cylinder, with outer radius r 1 = 12.0 and inner radius Relative error e Shape parameter q r 2 = 8.0, is subjected to an internal pressure p.The analytical solution of this problem can be obtained from [30] for the Mooney-Rivlin strain energy density function.This problem is modeled by a quarter of the plane strain rubber tube using 279 (31 × 9) regular distributed nodes as shown in Figure 5.In computation, the material constants A 10 = 0.2599, A 01 = 0.1608, and K = 10 5 , the shape parameters of MQ function C = 1.0, 2.0, 3.0, and q = 1.03, and the linear polynomial basis (m = 3) are used.Displacement control is used in this analysis and a total of 13 steps are applied to inflate the inner radius of the tube from 8 to 21.For the purpose of parameters and convergence studies, the relative error of energy is defined as where U MLPG and U Exat are the total strain energies computed by the present MLPG method and analytical method, respectively.The total strain energy can be expressed as where W is the modified Mooney-Rivlin strain energy density function.Figure 6 gives the variation of internal pressure p with internal radial displacement while the scaling coefficients are adopted as α s = 0.75, and α Q = 3.5.Comparing the present MLPG results with the analytic solutions it can be found that a good agreement is obtained.The relative errors of energy versus the subdomain scaling coefficients α s and interpolation domain scaling parameters α Q are plotted in Figures 7 and 8, when the internal pressure is equal to 0.262.It can be seen from Figures 7 and 8 that the scaling coefficients α s and α Q affect greatly the accuracy of the present MLPG results.In Figure 7, the excellent accuracy is obtained when α s varies from 0.4 to 1.0 with α Q = 3.5.In Figure 8, the results of the present MLPG method are acceptable only for α Q ranging from 2.5 to 5.0 with α s = 0.8, but not for α Q = 2.0.Table 1 gives computational efficiency of the MLPG method influenced by the size of interpolation domain, where the required computation time is set as 1.0 when α Q = 2.0.From Table 1, it can be found that the required computation time of the present method is increased greatly while enlarging the size of interpolation domain.
The convergence with different node distribution of the present method is studied for this problem.Three different node densities of 126 (21 × 6), 279 (31 × 9), and 451 (41 × 11) are used.The relative errors of energy and convergence rates are shown in Figure 9 with α s = 0.8 and α Q = 3.5.In Figure 9, R is average convergence rates; h is the maximum size of nodal interval.Figure 9 shows that the present method possesses high convergence.

Rubber Beam under Pure
Bending.A 20.0×1.0 cantilever beam as shown in Figure 10 is subjected to a moment on the free end by applying linear normal traction with zero resultant normal force.The material constants A 10 = 1835, A 01 = 146, and K = 10 7 are used.The scaling coefficients are adopted as α s = 0.75 and α Q = 3.5.The shape parameters of MQ function C = 1.42, q = 1.03 and the linear polynomial basis (m = 3) are selected in computation.Unless the computations fail to converge during the loading, the maximum moment M = 500 is applied to the plane strain beam.To generate pure bending deformation, the center node is fixed, while the other nodes at the other end of the beam are only restrained along the axial direction.
Figure 11 shows that effect of results by the integration points and the regular distributed nodes.The results indicate that the node distribution along the thickness direction is critical to solution accuracy.The model with three nodes along thickness direction behaves much stiffer than the analytical solutions [31], while the model with five nodes in thickness direction produces good results.It also can be seen from Figure 11 that the numerical accuracy is slightly variation for the number of Gauss integration points.In Figure 11, GN = a, b; GN is the abbreviation of the number of Gauss integration point; a and b are the number of Gauss integration points that adopted on each section of 0 L s and 0 Γ s for the numerical integration, respectively.c * d; c and d are the number of distributed nodes along the thickness and axial direction, respectively.The comparison of the MLPG method and FEM with regular and irregular meshes is presented in Figure 12.Table 2 compares the MLPG and FEM deflection solutions at M = 125.In FEM, the 9-node Lagrange element with the linear pressure is used that does not have locking or pressure oscillation problem according to the BB condition [32].The FEM calculations diverge at M = 198.2using irregular mesh and at M = 226.8using regular mesh because the meshes are distorted.The MLPG method performs without numerical divergences or loss of accuracy at¡?ehlt?¿ both regular and irregular nodes distribution.It also can be seen from Table 2 that the accuracy of the MLPG solution is better than the FEM.

Conclusions and Discussions
The present meshless method based on the local Petrov-Galerkin formulation is developed for the hyperelastic material problems.Using this meshless method, the structural domain is discretized by a set of nodes, and the   or nearly incompressible large deformation problems.The results of numerical examples indicate that the present method possesses a good performance for the hyperelastic material problems.

Figure 2 :
Figure 2: Effect of relative error e by shape parameter (m = 0).

Figure 8 :
Figure 8: Relative errors of energy versus interpolation domain parameter α Q .

Figure 12 :
Figure 12: Comparison of the MLPG method and the FEM.

Table 1 :
Computational efficiency of the MLPG method influenced by the interpolation domain parameter α Q

Table 2 :
Comparison of the MLPG method and the FEM at M = 125.