# The Meshless Local Petrov-Galerkin Method for Large Deformation Analysis of Hyperelastic Materials

Authors: De'an Hu; Zhanhua Sun
Journal: ISRN Mechanical Engineering
 (2011)

---

## Abstract

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.


---

## Body

## 1. 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], hp-cloud Method [9], Reproducing Kernel Particle Method (RKPM) [10], Natural Element Method (NEM) [11], Meshless Galerkin Method using Radial Basis Function (RBF) [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–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 trial-function 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–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–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.


## 2. 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(1)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
(2)Fij=∂xi∂Xj=∂ui∂Xj+δij,(3)Eij=12(FkiFkj-δij),(4)Gij=FkiFkj,(5)Sij=∂W∂Eij,
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(6)W(I̅1,I̅2,J)=W̅(I̅1,I̅2)+W̃(J),W̅(I̅1,I̅2)=A10(I̅1-3)+A01(I̅2-3),W̃(J)=K2(J-1)2,
where W̅(I̅1,I̅2) and W̃(J) are the deviatoric and volumetric terms of the strain energy density function, respectively. A10 and A01 are material constants, I̅1 and I̅2 are the reduced invariants, given by(7)I̅1=I1I3-1/3,I̅2=I2I3-2/3,I1, I2, I3 are the first, second, and third invariants of the Green deformation tensor, respectively, and J=I31/2. The second Piola-Kirchhoff stress associated with the modified Mooney-Rivlin strain energy density function defined in (5) is(8)Sij=(2A10I3-1/3+2A01I3-2/3I1)δij-2A01I3-2/3Gij-(23A10I3-1/3I1+43A01I3-2/3I2)Gij-1-PI31/2Gij-1,
where the hydrostatic pressure P is related to W̃(J) by(9)P=-∂W̃∂J=-K(J-1).
The incremental stress-strain relation can be expressed as(10)ΔSij=CijklΔEkl+gijΔP,
where Cijkl is the fourth-order material response tensor and gij is the second-order material response tensor, obtained by(11)Cijkl=C̅ijkl+C̃ijkl=∂2W̅∂Eij∂Ekl+∂(-JGij-1)∂EklP,gij=-I31/2Gij-1.
The explicit expressions of C̅ijkl, C̃ijkl, and gij are given as(12)C̅ijkl=∂2W̅∂Eij∂Ekl=(49A10I3-1/3I1+169A01I3-2/3I2)Gij-1Gkl-1-(43A10I3-1/3+83A01I3-2/3I1)(δijGkl-1+Gij-1δkl)+(23A10I3-1/3I1+43A01I3-2/3I2)(Gik-1Gjl-1+Gil-1Gjk-1)+4A01I3-2/3δijδkl-2A01I3-2/3(δikδjl+δilδjk)+83A01I3-2/3(GijGkl-1+Gij-1Gkl),C̃ijkl=∂(-JGij-1)∂EklP=PJ(Gik-1Gjl-1+Gil-1Gjk-1-Gij-1Gkl-1),gij=-I31/2Gij-1,
where Gij-1 is the ijth component of G-1 not 1/Gij.


## 3. 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 xi, i=1,2,…,N, the approximate function uh(x) of u, for all x∈Ωx,can be defined by [26](13)uh(x)=∑i=1nRi(x)ai+∑j=1mpj(x)bj=RT(x)a+PT(x)b,∀x∈Ωx,x=[x,y,z]T,
where ai are coefficients for the radial basis function Ri(x) and bj are coefficients for the polynomial basis pj(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(14)aT=[a1a2⋯an],bT=[b1b2⋯bm],RT(x)=[R1(x)R2(x)⋯Rn(x)],PT(x)=[p1(x)p2(x)⋯pm(x)].
The radial basis function Ri(x) has many kinds of forms and the multiquadrics (MQ) function is given by(15)Ri(x,y)=(ri2+c2)q=[(x-xi)2+(y-yi)2+c2]q,
where c,q are shape parameters.The polynomial basis has the following monomial terms in the two-dimensional problem(16)pT(x,y)=[1xyx2xyy2⋯].
The coefficients ai and bj 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(17)uk=u(xk,yk)=∑i=1naiRi(xk,yk)+∑j=1mbjPj(xk,yk),k=1,2,…,n,
which cannot be solved because n equations and n+m unknown coefficients are involved, so constraint conditions have to be introduced as the from(18)∑i=1naiPj(xi,yi)=0,j=1,2,…,m.
Equations (17) and (18) can be expressed in matrix form as follows:(19)A{ab}={us0},
where (20)A=[RQPmPmT0],RQ=[R1(x1,y1)R2(x1,y1)⋯Rn(x1,y1)R1(x2,y2)R2(x2,y2)⋯Rn(x2,y2)⋮⋮⋮⋮R1(xn,yn)R2(xn,yn)⋯Rn(xn,yn)]n×n,Pm=[p1(x1,y1)p2(x1,y1)⋯pm(x1,y1)p1(x2,y2)p2(x2,y2)⋯pm(x2,y2)⋮⋮⋮⋮p1(xn,yn)p2(xn,yn)⋯pm(xn,yn)]n×m.
Suppose A-1 exists, (19) can be rewritten in the following form:(21){ab}=A-1{us0}.
Using (21), we have(3)a=Saus,b=Sbus,Sb=[PmTRQ-1Pm]-1PmTRQ-1,Sa=RQ-1[I-PmSb]=RQ-1-RQ-1PmSb,
where PmTRQ-1Pm is termed as a transformed moment matrix.Finally, (13) is expressed as(23)uh(x,xQ)=[RTSa+PTSb]us=Φ(x)us=∑i=1nφi(x)ui,
where the shape functions Φ(x) can be written as follows:(24)Φ(x)=[φ1(x)φ2(x)⋅⋅⋅φn(x)],
in which the shape function φk(x) for the kth node is given by(25)φk(x)=∑i=1nRi(x)Sika+∑j=1mPj(x)Sjkb,
where Sika is the (i,k) element of matrix Sa and Sjkb is the (j,k) element of matrix Sb. The derivatives of φk(x)can be obtained as follows:(26)∂φk∂x=∑i=1n∂Ri∂xSika+∑j=1m∂Pj∂xSjkb,∂φk∂y=∑i=1n∂Ri∂ySika+∑j=1m∂Pj∂ySjkb
The partial derivatives of MQ function can be obtained as follows:(27)∂Ri∂x=2q(x-xi)[(x-xi)2+(y-yi)2+c2]q-1,∂Ri∂y=2q(y-yi)[(x-xi)2+(y-yi)2+c2]q-1.


## 4. 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 Γut, and natural boundary Γtt. 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(28)σij,j+bi=0,inΩt,
where σij is the Cauchy stress and bi is the body force. The boundary conditions are given as
(29)ui=u̅i,on⁡Γut,ti=σijnj=t̅i,on⁡Γtt,
where u̅i and t̅i are the prescribed displacements and tractions on boundary Γut and Γtt, respectively. nj is the unit outward normal to the boundary Γt.A generalized local weak form of (28) and (29) over a local subdomain Ωst can be written as follows:(30)∫tΩs(σij,j+bi)vidΩ=0,
where ui and vi are the trial and test function, respectively.Using(31)σij,jvi=(σijvi),j-σijvi,j,
and divergence theorem in (30) leads to(32)∫Ωst(-σijvi,j+bivi)dΩ+∫∂ΩsttividΓ=0,
where ti=σijnj and nj is a unit outward normal to the boundary ∂Ωst. ∂Ωst is the boundary of subdomain Ωst.It should be mentioned that (32) holds regardless of the size and the shape of Ωst, so the shape of subdomain Ωst can be taken to be a circle in 2D problems without losing generality.Applying the natural boundary condition,ti=σijnj=t̅i on Γstt where Γstt=∂Ωst∩Γtt, we get(33)∫LsttividΓ+∫ΓsuttividΓ+∫Γsttt̅ividΓ+∫Ωst(-σijvi,j+bivi)dΩ=0,
where Γsut and Γstt are the part of the boundary ∂Ωst over which the essential and natural boundary conditions are specified, respectively. In general, ∂Ωst=Γst∪Lst with Γst being a part of the local boundary located on the global boundary and Lst being other part of the local boundary over which no boundary condition is specified.In order to simplify (33), test functions vi are chosen such that they eliminate or simplify the domain integral on Ωst. This can be accomplished by using the Heaviside step function(34)H(x)={ζ,x∈Ωs,0,x∉Ωs,
where ζ is an arbitrary constant (ζ=1 is used in this study). Using this choice, the partial derivatives of the test function vi,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:(35)∫LsttidΓ+∫ΓsuttidΓ=-∫Γsttt̅idΓ-∫ΩstbidΩ.
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(36)t0dΓ0=FjkSikNjdΓ0=σijnjdΓt=tdΓt,
where the boundary Γ0 is referred to the initial configuration and Nj is the initial unit outward normal to the boundary Γ0.Then (35) can be rewritten as(37)∫Ls0FjkSikNjdΓ+∫Γsu0FjkSikNjdΓ=-∫Γst0t̅0idΓ-∫Ωs0b0idΩ,
where t̅0i and b0i 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:(38)ui=ui+Δui,Fij=∂(ui+Δui)∂Xj+δij=∂ui∂Xj+∂Δui∂Xj+δij=Fij+∂Δui∂Xj,Sij=Sij+ΔSij.
Substituting (38) into (37) leads to(39)∫Ls0FjkSikNjdΓ+∫Ls0FjkΔSikNjdΓ+∫Ls0∂Δuj∂XkSikNjdΓ+∫Ls0∂Δuj∂XkΔSikNjdΓ+∫Γsu0FjkSikNjdΓ+∫Γsu0FjkΔSikNjdΓ+∫Γsu0∂Δuj∂XkSikNjdΓ+∫Γsu0∂Δuj∂XkΔSikNjdΓ=-∫Γst0t̅0idΓ-∫Ωs0b0idΩ.
Equation (39) is a nonlinear integral equation that exist; quadratic and cubic terms of Δui after ΔSij are denoted by ΔEij and ΔP using (10) and the Green-Lagrangian strain increment(40)ΔEij=12(Fki∂Δuk∂Xj+Fkj∂Δuk∂Xi+∂Δuk∂Xi∂Δuk∂Xj).
Then linearization must be carried out for solving the integral equation. The linearization of (39) is given by(41)12∫Ls0FjkCikml(Fpm∂Δup∂Xl+Fpl∂Δup∂Xm)NjdΓ+12∫Γsu0FjkCikml(Fpm∂Δup∂Xl+Fpl∂Δup∂Xm)NjdΓ+∫Ls0∂Δuj∂XkSikNjdΓ+∫Γsu0∂Δuj∂XkSikNjdΓ+∫Ls0FjkgikNjΔPdΓ+∫Γsu0FjkgikNjΔPdΓ=-∫Γst0Δt̅0idΓ-∫Ωs0Δb0idΩ-∫Ls0FjkSikNjdΓ-∫Γsu0FjkSikNjdΓ.
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](42)P=λ3(2A10λ3-2A01λ3-3)=2A10λ12λ22-2A01λ12λ22,
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(43)λ32=G33=λ1-2λ2-2=(G11G22-G122)-1.
Then(44)P=2A10(G11G22-G122)-1-2A01(G11G22-G122).
Using (41), the hydrostatic pressure increment is given by(45)ΔP=LΔE,
where(46)LT={-4A10(G11G22-G122)-2G22-4A01G22-4A10(G11G22-G122)-2G11-4A01G118A10(G11G22-G122)-2G12+8A01G12},ΔE=[ΔE11ΔE22ΔE12]T.
For the states of plane strain, we reduce the independent constraint equations in each subdomain Ωs0 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,(47)Θ(P)=‖-K(J-1)-P‖L2(Ωs0)2,
where a constant field is used as the lower-order space, and ∥·∥L2(Ωs0)2 is the L2 norm in the subdomain Ωs0. The minimization of Θ(P) leads to(48)P=-K∫Ωs0(J-1)dΩAs.
Similarly, in each subdomainΩs0, we project the hydrostatic pressure increment onto a constant field by the least squares approximation(49)Θ̂(ΔP)=‖-KJGkl-1ΔEkl-ΔP‖L2(Ωs0)2.
The minimization of Θ̂(ΔP) leads to(50)ΔP=-K∫Ωs0gTBMdΩAsΔu,
where As is the area of subdomain Ωs0. BM is the gradient matrix of Green-Lagrangian strain. If the most severe over-constrained condition is encountered, one-point integration to calculate (48) and (50) can be employed.


## 5. Discretizations

Using the radial basis function coupled with polynomial basis function to construct the trial function,(51)Δuh(x)=∑n=1Nφn(x)Δun.
Substituting (51) into (41) leads to the discrete iterative equation(52)(K)n+1vδu=(f)n+1v,
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(53)KIJ=∫Ls0NF̃CBJMdΓ+∫Γsu0NF̃CBJMdΓ+∫Ls0NTBJGdΓ+∫Γsu0NTBJGdΓ+KIJ**,fI=-∫Γst0t̅dΓ-∫Ωs0bdΩ-∫Ls0NF̃ŜdΓ-∫Γsu0NF̃ŜdΓ.
The term KIJ** results from the hydrostatic pressure increment and can be expressed as(54)KIJ**={∫Ls0NF̃gLBJMdΓ+∫Γsu0NF̃gLBJMdΓfortheplanestress,-KAs-1[(∫Ls0NF̃gdΓ)+(∫Γsu0NF̃gdΓ)]×(∫Ωs0gTBMdΩ)fortheplanestrain.
For the plane problem, the explicit expressions of vectors and matrixes in the above equations are defined by(55)φn=[φ1(x)00φ2(x)]n,N=[N10N200N20N1],C=[C1111C1122C1112C2222C2212symC1212],F̃=[F110F120F22F210F12F11F210F22],T=[S110S1200S220S21S210S2200S120S11],Ŝ={S11S22S12},g={g11g22g12},BJM=[F11∂φ1∂X1F21∂φ2∂X1F12∂φ1∂X2F22∂φ2∂X2F11∂φ1∂X2+F12∂φ1∂X1F21∂φ2∂X2+F22∂φ2∂X1]J,BJG=[∂φ1∂X100∂φ2∂X2∂φ1∂X200∂φ2∂X1]J.
From (53) and (54), it can be found that the domain integral over Ωs0 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.


## 6. 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 ofLs0 and Γs0 for the numerical integration, respectively. Radii of the influence domain rQ is set to αQdi3, di3 is the distance to the third nearest neighboring point from node i. Radii of the subdomain rs is set to αsdi1, di1 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.


### 6.1. Uniaxial Tension-Compression Problem

An8.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(56)e=|(σMLPG-σExat)σExat|,
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 Figure1 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 A10=0.2599 and A01=0.1608 are used. The scaling coefficients are adopted as αs=0.75 and αQ=3.5.Figure 1
Rubber block.Figures2 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.Figure 2
Effect of relative errore by shape parameter (m=0).Figure 3
Effect of relative errore by shape parameter (m=3).Figure 4
Cauchy stress versus deformation ratioλ at axial direction.


### 6.2. Inflation of a Long Rubber Tube

An infinitely long rubber cylinder, with outer radiusr1=12.0 and inner radius r2=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 using279(31×9) regular distributed nodes as shown in Figure 5. In computation, the material constants A10=0.2599, A01=0.1608, and K=105, 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(57)eE=|UMLPG-UExatUExat|,
where UMLPG and UExat are the total strain energies computed by the present MLPG method and analytical method, respectively. The total strain energy can be expressed as(58)U=∫Ω0W(I̅1,I̅2,J)dΩ,
where W is the modified Mooney-Rivlin strain energy density function.Figure 5
Inflation of long rubber tube.Figure6 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.Table 1
Computational efficiency of the MLPG method influenced by the interpolation domain parameterαQ(C=2.0,q=1.03).
ParameterαQ=2.0αQ=3.0αQ=4.0αQ=5.0Normalized time1.02.364.245.92Figure 6
Internal pressurep versus internal radial displacement.Figure 7
Relative errors of energy versus subdomain parameterαs.Figure 8
Relative errors of energy versus interpolation domain parameterαQ.The convergence with different node distribution of the present method is studied for this problem. Three different node densities of126(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.Figure 9
Relative error and convergence Rate of energy norm.


### 6.3. Rubber Beam under Pure Bending

A20.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 A10=1835, A01=146, and K=107 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 10
Rubber beam under pure bending.Figure11 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 Ls0 and Γs0 for the numerical integration, respectively. c*d; c and d are the number of distributed nodes along the thickness and axial direction, respectively.Figure 11
Effect of results by integral points and distributed nodes.The comparison of the MLPG method and FEM with regular and irregular meshes is presented in Figure12. 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.2 using irregular mesh and at M=226.8 using regular mesh because the meshes are distorted. The MLPG method performs without numerical divergences or loss of accuracy at 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.Table 2
Comparison of the MLPG method and the FEM atM=125.
AnalyticalRegularIrregularMLPGFEMMLPGFEM13.9213.9013.2013.8712.26Figure 12
Comparison of the MLPG method and the FEM.


## 6.1. Uniaxial Tension-Compression Problem

An8.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(56)e=|(σMLPG-σExat)σExat|,
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 Figure1 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 A10=0.2599 and A01=0.1608 are used. The scaling coefficients are adopted as αs=0.75 and αQ=3.5.Figure 1
Rubber block.Figures2 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.Figure 2
Effect of relative errore by shape parameter (m=0).Figure 3
Effect of relative errore by shape parameter (m=3).Figure 4
Cauchy stress versus deformation ratioλ at axial direction.


## 6.2. Inflation of a Long Rubber Tube

An infinitely long rubber cylinder, with outer radiusr1=12.0 and inner radius r2=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 using279(31×9) regular distributed nodes as shown in Figure 5. In computation, the material constants A10=0.2599, A01=0.1608, and K=105, 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(57)eE=|UMLPG-UExatUExat|,
where UMLPG and UExat are the total strain energies computed by the present MLPG method and analytical method, respectively. The total strain energy can be expressed as(58)U=∫Ω0W(I̅1,I̅2,J)dΩ,
where W is the modified Mooney-Rivlin strain energy density function.Figure 5
Inflation of long rubber tube.Figure6 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.Table 1
Computational efficiency of the MLPG method influenced by the interpolation domain parameterαQ(C=2.0,q=1.03).
ParameterαQ=2.0αQ=3.0αQ=4.0αQ=5.0Normalized time1.02.364.245.92Figure 6
Internal pressurep versus internal radial displacement.Figure 7
Relative errors of energy versus subdomain parameterαs.Figure 8
Relative errors of energy versus interpolation domain parameterαQ.The convergence with different node distribution of the present method is studied for this problem. Three different node densities of126(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.Figure 9
Relative error and convergence Rate of energy norm.


## 6.3. Rubber Beam under Pure Bending

A20.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 A10=1835, A01=146, and K=107 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 10
Rubber beam under pure bending.Figure11 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 Ls0 and Γs0 for the numerical integration, respectively. c*d; c and d are the number of distributed nodes along the thickness and axial direction, respectively.Figure 11
Effect of results by integral points and distributed nodes.The comparison of the MLPG method and FEM with regular and irregular meshes is presented in Figure12. 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.2 using irregular mesh and at M=226.8 using regular mesh because the meshes are distorted. The MLPG method performs without numerical divergences or loss of accuracy at 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.Table 2
Comparison of the MLPG method and the FEM atM=125.
AnalyticalRegularIrregularMLPGFEMMLPGFEM13.9213.9013.2013.8712.26Figure 12
Comparison of the MLPG method and the FEM.


## 7. 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 use of mesh in the classical FEM is unnecessary. In this paper, Trial functions are constructed using the radial basis function coupled with polynomial basis function, so no additional treatments are required for imposing essential boundary conditions. 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. The Heaviside test function is chosen for eliminating or simplifying domain integrals in the global stiffness matrix. Finally, through the studying of numerical examples, various shape parametersC and q, scaling coefficients αs and αQ, the number of Gauss integration points, and the different node distribution are investigated in the present method for the incompressible 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.



---

*Source: 967512-2011-09-05.xml*