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.

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. [1315]; 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 [2527], 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 [2830] 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 Ω𝑋 with boundary Γ𝑋. The deformation of a material particle 𝐗Ω𝑋 at time 𝑡 is described by 𝐱(𝐗,𝑡), and the displacement of the particle 𝐗 is defined by𝐮(𝐗,𝑡)=𝐱(𝐗,𝑡)𝐗.(1) We consider the reference frame to be a rectangular Cartesian system. The deformation gradient 𝐹, the Green-Lagrangian strain 𝐸, the Green deformation tensor 𝐺, and the second Piola-Kirchhoff stress 𝑆 are 𝐹𝑖𝑗=𝜕𝑥𝑖𝜕𝑋𝑗=𝜕𝑢𝑖𝜕𝑋𝑗+𝛿𝑖𝑗,𝐸(2)𝑖𝑗=12𝐹𝑘𝑖𝐹𝑘𝑗𝛿𝑖𝑗𝐺,(3)𝑖𝑗=𝐹𝑘𝑖𝐹𝑘𝑗𝑆,(4)𝑖𝑗=𝜕𝑊𝜕𝐸𝑖𝑗,(5) where 𝑊 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 𝐾 is considered. The incompressibility can be recovered by letting 𝐾 tend to infinity. In the present analysis, the modified Mooney-Rivlin strain energy density function 𝑊 is used𝑊𝐼1,𝐼2=,𝐽𝑊𝐼1,𝐼2+𝑊(𝐽),𝑊𝐼1,𝐼2=𝐴10𝐼13+𝐴01𝐼2,𝑊𝐾3(𝐽)=2(𝐽1)2,(6) where 𝑊(𝐼1,𝐼2) and 𝑊(𝐽) are the deviatoric and volumetric terms of the strain energy density function, respectively. 𝐴10 and 𝐴01 are material constants, 𝐼1 and 𝐼2 are the reduced invariants, given by𝐼1=𝐼1𝐼31/3,𝐼2=𝐼2𝐼32/3,(7)𝐼1, 𝐼2, 𝐼3 are the first, second, and third invariants of the Green deformation tensor, respectively, and 𝐽=𝐼31/2. The second Piola-Kirchhoff stress associated with the modified Mooney-Rivlin strain energy density function defined in (5) is𝑆𝑖𝑗=2𝐴10𝐼31/3+2𝐴01𝐼32/3𝐼1𝛿𝑖𝑗2𝐴01𝐼32/3𝐺𝑖𝑗23𝐴10𝐼31/3𝐼1+43𝐴01𝐼32/3𝐼2𝐺1𝑖𝑗𝑃𝐼31/2𝐺1𝑖𝑗,(8) where the hydrostatic pressure 𝑃 is related to 𝑊(𝐽) by𝜕𝑊𝑃=𝜕𝐽=𝐾(𝐽1).(9) The incremental stress-strain relation can be expressed asΔ𝑆𝑖𝑗=𝐶𝑖𝑗𝑘𝑙Δ𝐸𝑘𝑙+𝑔𝑖𝑗Δ𝑃,(10) where 𝐶𝑖𝑗𝑘𝑙 is the fourth-order material response tensor and 𝑔𝑖𝑗 is the second-order material response tensor, obtained by𝐶𝑖𝑗𝑘𝑙=𝐶𝑖𝑗𝑘𝑙+𝐶𝑖𝑗𝑘𝑙=𝜕2𝑊𝜕𝐸𝑖𝑗𝜕𝐸𝑘𝑙+𝜕𝐽𝐺1𝑖𝑗𝜕𝐸𝑘𝑙𝑔𝑃,𝑖𝑗=𝐼31/2𝐺1𝑖𝑗.(11) The explicit expressions of 𝐶𝑖𝑗𝑘𝑙, 𝐶𝑖𝑗𝑘𝑙, and 𝑔𝑖𝑗 are given as𝐶𝑖𝑗𝑘𝑙=𝜕2𝑊𝜕𝐸𝑖𝑗𝜕𝐸𝑘𝑙=49𝐴10𝐼31/3𝐼1+169𝐴01𝐼32/3𝐼2𝐺1𝑖𝑗𝐺1𝑘𝑙43𝐴10𝐼31/3+83𝐴01𝐼32/3𝐼1𝛿𝑖𝑗𝐺1𝑘𝑙+𝐺1𝑖𝑗𝛿𝑘𝑙+23𝐴10𝐼31/3𝐼1+43𝐴01𝐼32/3𝐼2𝐺1𝑖𝑘𝐺1𝑗𝑙+𝐺1𝑖𝑙𝐺1𝑗𝑘+4𝐴01𝐼32/3𝛿𝑖𝑗𝛿𝑘𝑙2𝐴01𝐼32/3𝛿𝑖𝑘𝛿𝑗𝑙+𝛿𝑖𝑙𝛿𝑗𝑘+83𝐴01𝐼32/3𝐺𝑖𝑗𝐺1𝑘𝑙+𝐺1𝑖𝑗𝐺𝑘𝑙,𝐶𝑖𝑗𝑘𝑙=𝜕𝐽𝐺1𝑖𝑗𝜕𝐸𝑘𝑙𝐺𝑃=𝑃𝐽1𝑖𝑘𝐺1𝑗𝑙+𝐺1𝑖𝑙𝐺1𝑗𝑘𝐺1𝑖𝑗𝐺1𝑘𝑙,𝑔𝑖𝑗=𝐼31/2𝐺1𝑖𝑗,(12) where 𝐺1𝑖𝑗 is the 𝑖𝑗th component of 𝐆1 not 1/𝐺𝑖𝑗.

3. Interpolation Approximation

Consider a domain Ω𝑥, the neighborhood of point 𝐱, which is located within the problem domain Ω. To approximate the distribution of function 𝑢 in Ω𝑥 over a number of randomly located nodes 𝑥𝑖, 𝑖=1,2,,𝑁, the approximate function 𝑢(𝐱) of 𝑢, for all 𝐱Ω𝑥,can be defined by [26]𝑢(𝐱)=𝑛𝑖=1𝑅𝑖(𝐱)𝑎𝑖+𝑚𝑗=1𝑝𝑗(𝐱)𝑏𝑗=𝐑𝑇(𝐱)𝐚+𝐏𝑇(𝐱)𝐛,𝐱Ω𝑥[],𝐱=𝑥,𝑦,𝑧𝑇,(13) where 𝑎𝑖 are coefficients for the radial basis function 𝑅𝑖(𝐱) and 𝑏𝑗 are coefficients for the polynomial basis 𝑝𝑗(𝐱). The number of the radial basis 𝑛 is determined by the number of nodes in the support domain Ω𝑥, and the number of the polynomial basis 𝑚 can be chosen based on the reproduction requirement. A minimum number of terms of the polynomial basis and more terms of the radial basis (𝑚<𝑛) are adopted for better stability. The vectors in (13) are defined as𝐚𝑇=𝑎1𝑎2𝑎𝑛,𝐛𝑇=𝑏1𝑏2𝑏𝑚,𝐑𝑇𝑅(𝐱)=1(𝐱)𝑅2(𝐱)𝑅𝑛(,𝐏𝐱)𝑇(𝑝𝐱)=1(𝐱)𝑝2(𝐱)𝑝𝑚.(𝐱)(14) The radial basis function 𝑅𝑖(𝐱) has many kinds of forms and the multiquadrics (MQ) function is given by𝑅𝑖𝑟(𝑥,𝑦)=2𝑖+𝑐2𝑞=𝑥𝑥𝑖2+𝑦𝑦𝑖2+𝑐2𝑞,(15) where 𝑐,𝑞 are shape parameters.

The polynomial basis has the following monomial terms in the two-dimensional problem𝐩𝑇(𝑥,𝑦)=1𝑥𝑦𝑥2𝑥𝑦𝑦2.(16) The coefficients 𝑎𝑖 and 𝑏𝑗 in (13) are determined by enforcing that the interpolation function passes through all nodes within the support domain. The interpolation at the 𝑘th point has the form𝑢𝑘𝑥=𝑢𝑘,𝑦𝑘=𝑛𝑖=1𝑎𝑖𝑅𝑖𝑥𝑘,𝑦𝑘+𝑚𝑗=1𝑏𝑗𝑃𝑗𝑥𝑘,𝑦𝑘,𝑘=1,2,,𝑛,(17) which cannot be solved because 𝑛 equations and 𝑛+𝑚 unknown coefficients are involved, so constraint conditions have to be introduced as the from𝑛𝑖=1𝑎𝑖𝑃𝑗𝑥𝑖,𝑦𝑖=0,𝑗=1,2,,𝑚.(18) Equations (17) and (18) can be expressed in matrix form as follows:𝐀𝐚𝐛=𝐮𝐬𝟎,(19) where 𝐑𝐀=𝑄𝐏𝑚𝐏𝑇𝑚𝟎,𝐑𝑄=𝑅1𝑥1,𝑦1𝑅2𝑥1,𝑦1𝑅𝑛𝑥1,𝑦1𝑅1𝑥2,𝑦2𝑅2𝑥2,𝑦2𝑅𝑛𝑥2,𝑦2𝑅1𝑥𝑛,𝑦𝑛𝑅2𝑥𝑛,𝑦𝑛𝑅𝑛𝑥𝑛,𝑦𝑛𝑛×𝑛,𝐏𝑚=𝑝1𝑥1,𝑦1𝑝2𝑥1,𝑦1𝑝𝑚𝑥1,𝑦1𝑝1𝑥2,𝑦2𝑝2𝑥2,𝑦2𝑝𝑚𝑥2,𝑦2𝑝1𝑥𝑛,𝑦𝑛𝑝2𝑥𝑛,𝑦𝑛𝑝𝑚𝑥𝑛,𝑦𝑛𝑛×𝑚.(20) Suppose 𝐀1 exists, (19) can be rewritten in the following form:𝐚𝐛=𝐀𝟏𝐮𝐬𝟎.(21) Using (21), we have𝐚=𝑆𝑎𝐮𝑠,𝐛=𝑆𝑏𝐮𝑠,𝑆𝑏=𝐏𝑇𝑚𝐑𝑄1𝐏𝑚1𝐏𝑇𝑚𝐑𝑄1,𝑆𝑎=𝐑𝑄1𝐈𝐏𝑚𝑆𝑏=𝐑𝑄1𝐑𝑄1𝐏𝑚𝑆𝑏,(3) where 𝐏𝑇𝑚𝐑𝑄1𝐏𝑚 is termed as a transformed moment matrix.

Finally, (13) is expressed as𝑢𝐱,𝐱𝑄=𝐑𝑇𝑆𝑎+𝐏𝑇𝑆𝑏𝐮𝑠=𝚽(𝐱)𝐮𝑠=𝑛𝑖=1𝜑𝑖(𝐱)𝑢𝑖,(23) where the shape functions Φ(𝐱) can be written as follows:𝚽𝜑(𝐱)=1(𝐱)𝜑2(𝐱)𝜑𝑛(𝐱),(24) in which the shape function 𝜑𝑘(𝐱) for the 𝑘th node is given by𝜑𝑘(𝐱)=𝑛𝑖=1𝑅𝑖(𝐱)𝑆𝑎𝑖𝑘+𝑚𝑗=1𝑃𝑗(𝐱)𝑆𝑏𝑗𝑘,(25) where 𝑆𝑎𝑖𝑘 is the (𝑖,𝑘) element of matrix 𝑆𝑎 and 𝑆𝑏𝑗𝑘 is the (𝑗,𝑘) element of matrix 𝑆𝑏. The derivatives of 𝜑𝑘(𝐱)can be obtained as follows:𝜕𝜑𝑘=𝜕𝑥𝑛𝑖=1𝜕𝑅𝑖𝑆𝜕𝑥𝑎𝑖𝑘+𝑚𝑗=1𝜕𝑃𝑗𝑆𝜕𝑥𝑏𝑗𝑘,𝜕𝜑𝑘=𝜕𝑦𝑛𝑖=1𝜕𝑅𝑖𝑆𝜕𝑦𝑎𝑖𝑘+𝑚𝑗=1𝜕𝑃𝑗𝑆𝜕𝑦𝑏𝑗𝑘(26) The partial derivatives of MQ function can be obtained as follows:𝜕𝑅𝑖𝜕𝑥=2𝑞𝑥𝑥𝑖𝑥𝑥𝑖2+𝑦𝑦𝑖2+𝑐2𝑞1,𝜕𝑅𝑖𝜕𝑦=2𝑞𝑦𝑦𝑖𝑥𝑥𝑖2+𝑦𝑦𝑖2+𝑐2𝑞1.(27)

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 𝑡Ω, essential boundary 𝑡Γ𝑢, and natural boundary 𝑡Γ𝑡. During the deformation process, a material point 𝐗0Ω is moved to 𝐱=𝐗+𝐮𝑡Ω. For this problem, the static equilibrium equation is described as𝜎𝑖𝑗,𝑗+𝑏𝑖=0,in𝑡Ω,(28) where 𝜎𝑖𝑗 is the Cauchy stress and 𝑏𝑖 is the body force. The boundary conditions are given as 𝑢𝑖=𝑢𝑖,on𝑡Γ𝑢,𝑡𝑖=𝜎𝑖𝑗𝑛𝑗=𝑡𝑖,on𝑡Γ𝑡,(29) where 𝑢𝑖 and 𝑡𝑖 are the prescribed displacements and tractions on boundary 𝑡Γ𝑢 and 𝑡Γ𝑡, respectively. 𝑛𝑗 is the unit outward normal to the boundary 𝑡Γ.

A generalized local weak form of (28) and (29) over a local subdomain 𝑡Ω𝑠 can be written as follows:𝑡Ω𝑠𝜎𝑖𝑗,𝑗+𝑏𝑖𝑣𝑖dΩ=0,(30) where 𝑢𝑖 and 𝑣𝑖 are the trial and test function, respectively.

Using𝜎𝑖𝑗,𝑗𝑣𝑖=𝜎𝑖𝑗𝑣𝑖,𝑗𝜎𝑖𝑗𝑣𝑖,𝑗,(31) and divergence theorem in (30) leads to𝑡Ω𝑠𝜎𝑖𝑗𝑣𝑖,𝑗+𝑏𝑖𝑣𝑖dΩ+𝜕𝑡Ω𝑠𝑡𝑖𝑣𝑖dΓ=0,(32) where 𝑡𝑖=𝜎𝑖𝑗𝑛𝑗 and 𝑛𝑗 is a unit outward normal to the boundary 𝜕𝑡Ω𝑠. 𝜕𝑡Ω𝑠 is the boundary of subdomain 𝑡Ω𝑠.

It should be mentioned that (32) holds regardless of the size and the shape of 𝑡Ω𝑠, so the shape of subdomain 𝑡Ω𝑠 can be taken to be a circle in 2D problems without losing generality.

Applying the natural boundary condition, 𝑡𝑖=𝜎𝑖𝑗𝑛𝑗=𝑡𝑖 on 𝑡Γ𝑠𝑡 where 𝑡Γ𝑠𝑡=𝜕𝑡Ω𝑠𝑡Γ𝑡, we get𝑡𝐿𝑠𝑡𝑖𝑣𝑖dΓ+𝑡Γ𝑠𝑢𝑡𝑖𝑣𝑖dΓ+𝑡Γ𝑠𝑡𝑡𝑖𝑣𝑖+dΓ𝑡Ω𝑠𝜎𝑖𝑗𝑣𝑖,𝑗+𝑏𝑖𝑣𝑖dΩ=0,(33) where 𝑡Γ𝑠𝑢 and 𝑡Γ𝑠𝑡 are the part of the boundary 𝜕𝑡Ω𝑠 over which the essential and natural boundary conditions are specified, respectively. In general, 𝜕𝑡Ω𝑠=𝑡Γ𝑠𝑡L𝑠 with 𝑡Γ𝑠 being a part of the local boundary located on the global boundary and 𝑡L𝑠 being other part of the local boundary over which no boundary condition is specified.

In order to simplify (33), test functions 𝑣𝑖 are chosen such that they eliminate or simplify the domain integral on 𝑡Ω𝑠. This can be accomplished by using the Heaviside step function𝐻(𝐱)=𝜁,𝐱Ω𝑠,0,𝐱Ω𝑠,(34) where 𝜁 is an arbitrary constant (𝜁=1 is used in this study). Using this choice, the partial derivatives of the test function 𝑣𝑖,𝑗 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:𝑡L𝑠𝑡𝑖dΓ+𝑡Γ𝑠𝑢𝑡𝑖dΓ=𝑡Γ𝑠𝑡𝑡𝑖dΓ𝑡Ω𝑠𝑏𝑖dΩ.(35) 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𝐭0d0Γ=𝐹𝑗𝑘𝑆𝑖𝑘𝑁𝑗d0Γ=𝜎𝑖𝑗𝑛𝑗d𝑡Γ=𝐭d𝑡Γ,(36) where the boundary 0Γ is referred to the initial configuration and 𝑁𝑗 is the initial unit outward normal to the boundary 0Γ.

Then (35) can be rewritten as0L𝑠𝐹𝑗𝑘𝑆𝑖𝑘𝑁𝑗dΓ+0Γ𝑠𝑢𝐹𝑗𝑘𝑆𝑖𝑘𝑁𝑗dΓ=0Γ𝑠𝑡𝑡0𝑖dΓ0Ω𝑠𝑏0𝑖dΩ,(37) where 𝑡0𝑖 and 𝑏0𝑖 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) Substituting (38) into (37) leads to0L𝑠𝐹𝑗𝑘𝑆𝑖𝑘𝑁𝑗dΓ+0L𝑠𝐹𝑗𝑘Δ𝑆𝑖𝑘𝑁𝑗dΓ+0L𝑠𝜕Δ𝑢𝑗𝜕𝑋𝑘𝑆𝑖𝑘𝑁𝑗+dΓ0L𝑠𝜕Δ𝑢𝑗𝜕𝑋𝑘Δ𝑆𝑖𝑘𝑁𝑗dΓ+0Γ𝑠𝑢𝐹𝑗𝑘𝑆𝑖𝑘𝑁𝑗+dΓ0Γ𝑠𝑢𝐹𝑗𝑘Δ𝑆𝑖𝑘𝑁𝑗+dΓ0Γ𝑠𝑢𝜕Δ𝑢𝑗𝜕𝑋𝑘𝑆𝑖𝑘𝑁𝑗dΓ+0Γ𝑠𝑢𝜕Δ𝑢𝑗𝜕𝑋𝑘Δ𝑆𝑖𝑘𝑁𝑗dΓ=0Γ𝑠𝑡𝑡0𝑖dΓ0Ω𝑠𝑏0𝑖dΩ.(39) Equation (39) is a nonlinear integral equation that exist; quadratic and cubic terms of Δ𝑢𝑖 after Δ𝑆𝑖𝑗 are denoted by Δ𝐸𝑖𝑗 and Δ𝑃 using (10) and the Green-Lagrangian strain incrementΔ𝐸𝑖𝑗=12𝐹𝑘𝑖𝜕Δ𝑢𝑘𝜕𝑋𝑗+𝐹𝑘𝑗𝜕Δ𝑢𝑘𝜕𝑋𝑖+𝜕Δ𝑢𝑘𝜕𝑋𝑖𝜕Δ𝑢𝑘𝜕𝑋𝑗.(40) Then linearization must be carried out for solving the integral equation. The linearization of (39) is given by120L𝑠𝐹𝑗𝑘𝐶𝑖𝑘𝑚𝑙𝐹𝑝𝑚𝜕Δ𝑢𝑝𝜕𝑋𝑙+𝐹𝑝𝑙𝜕Δ𝑢𝑝𝜕𝑋𝑚𝑁𝑗+1dΓ20Γ𝑠𝑢𝐹𝑗𝑘𝐶𝑖𝑘𝑚𝑙𝐹𝑝𝑚𝜕Δ𝑢𝑝𝜕𝑋𝑙+𝐹𝑝𝑙𝜕Δ𝑢𝑝𝜕𝑋𝑚𝑁𝑗+dΓ0L𝑠𝜕Δ𝑢𝑗𝜕𝑋𝑘𝑆𝑖𝑘𝑁𝑗dΓ+0Γ𝑠𝑢𝜕Δ𝑢𝑗𝜕𝑋𝑘𝑆𝑖𝑘𝑁𝑗+dΓ0L𝑠𝐹𝑗𝑘𝑔𝑖𝑘𝑁𝑗Δ𝑃dΓ+0Γ𝑠𝑢𝐹𝑗𝑘𝑔𝑖𝑘𝑁𝑗Δ𝑃dΓ=0Γ𝑠𝑡Δ𝑡0𝑖dΓ0Ω𝑠Δ𝑏0𝑖dΩ0L𝑠𝐹𝑗𝑘𝑆𝑖𝑘𝑁𝑗dΓ0Γ𝑠𝑢𝐹𝑗𝑘𝑆𝑖𝑘𝑁𝑗dΓ.(41) 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]𝑃=𝜆32𝐴10𝜆32𝐴01𝜆33=2𝐴10𝜆21𝜆222𝐴01𝜆21𝜆22,(42) 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𝜆23=𝐺33=𝜆12𝜆22=𝐺11𝐺22𝐺2121.(43) Then𝑃=2𝐴10𝐺11𝐺22𝐺21212𝐴01𝐺11𝐺22𝐺212.(44) Using (41), the hydrostatic pressure increment is given byΔ𝑃=𝐋Δ𝐄,(45) where𝐋𝑇=4𝐴10𝐺11𝐺22𝐺2122𝐺224𝐴01𝐺224𝐴10𝐺11𝐺22𝐺2122𝐺114𝐴01𝐺118𝐴10𝐺11𝐺22𝐺2122𝐺12+8𝐴01𝐺12,Δ𝐸=Δ𝐸11Δ𝐸22Δ𝐸12T.(46) For the states of plane strain, we reduce the independent constraint equations in each subdomain 0Ω𝑠 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,Θ(𝑃)=𝐾(𝐽1)𝑃2𝐿2(0Ω𝑠),(47) where a constant field is used as the lower-order space, and 2𝐿2(0Ω𝑠) is the 𝐿2 norm in the subdomain 0Ω𝑠. The minimization of Θ(𝑃) leads to𝑃=𝐾0Ω𝑠(𝐽1)dΩ𝐴𝑠.(48) Similarly, in each subdomain0Ω𝑠, we project the hydrostatic pressure increment onto a constant field by the least squares approximationΘ(Δ𝑃)=𝐾𝐽𝐺1𝑘𝑙Δ𝐸𝑘𝑙Δ𝑃2𝐿2(0Ω𝑠).(49) The minimization of Θ(Δ𝑃) leads toΔ𝑃=𝐾0Ω𝑠𝐠𝑇𝐁𝑀dΩ𝐴𝑠Δ𝐮,(50) where 𝐴𝑠 is the area of subdomain 0Ω𝑠. 𝐁𝑀 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,Δ𝐮(𝐱)=𝑁𝑛=1𝜑𝑛(𝐱)Δ𝐮𝑛.(51) Substituting (51) into (41) leads to the discrete iterative equation(𝐊)𝑣𝑛+1𝛿𝐮=(𝐟)𝑣𝑛+1,(52) where 𝑛 and 𝑣 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 𝐊 and the load vector 𝐟 are defined as𝐊𝐼𝐽=0L𝑠𝐍𝐅𝐂𝐁𝑀𝐽dΓ+0Γ𝑠𝑢𝐍𝐅𝐂𝐁𝑀𝐽+dΓ0L𝑠𝐍𝐓𝐁𝐺𝐽dΓ+0Γ𝑠𝑢𝐍𝐓𝐁𝐺𝐽dΓ+𝐊𝐼𝐽,𝐟𝐼=0Γ𝑠𝑡𝐭dΓ0Ω𝑠𝐛dΩ0L𝑠𝐍𝐅̂𝐒dΓ0Γ𝑠𝑢𝐍𝐅̂𝐒dΓ.(53) The term 𝐊𝐼𝐽 results from the hydrostatic pressure increment and can be expressed as𝐊𝐼𝐽=0L𝑠𝐍𝐅𝐠𝐋𝐁𝑀𝐽dΓ+0Γ𝑠𝑢𝐍𝐅𝐠𝐋𝐁𝑀𝐽dΓfortheplanestress,𝐾𝐴𝑠10L𝑠𝐍+𝐅𝐠dΓ0Γ𝑠𝑢𝐍×𝐅𝐠dΓ0Ω𝑠𝐠T𝐁𝑀dΩfortheplanestrain.(54) For the plane problem, the explicit expressions of vectors and matrixes in the above equations are defined by𝝋𝑛=𝜑1(𝐱)00𝜑2(𝐱)𝑛𝑁,𝐍=10𝑁200𝑁20𝑁1,𝐶𝐂=1111𝐶1122𝐶1112𝐶2222𝐶2212sym𝐶1212,𝐹𝐅=110𝐹120𝐹22𝐹210𝐹12𝐹11𝐹210𝐹22𝑆,𝐓=110𝑆1200𝑆220𝑆21𝑆210𝑆2200𝑆120𝑆11,̂𝑆𝐒=11𝑆22𝑆12𝑔,𝐠=11𝑔22𝑔12,𝐁𝑀𝐽=𝐹11𝜕𝜑1𝜕𝑋1𝐹21𝜕𝜑2𝜕𝑋1𝐹12𝜕𝜑1𝜕𝑋2𝐹22𝜕𝜑2𝜕𝑋2𝐹11𝜕𝜑1𝜕𝑋2+𝐹12𝜕𝜑1𝜕𝑋1𝐹21𝜕𝜑2𝜕𝑋2+𝐹22𝜕𝜑2𝜕𝑋1𝐽,𝐁𝐺𝐽=𝜕𝜑1𝜕𝑋100𝜕𝜑2𝜕𝑋2𝜕𝜑1𝜕𝑋200𝜕𝜑2𝜕𝑋1𝐽.(55) From (53) and (54), it can be found that the domain integral over 0Ω𝑠 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 of 0L𝑠 and 0Γ𝑠 for the numerical integration, respectively. Radii of the influence domain 𝑟𝑄 is set to 𝛼𝑄𝑑3𝑖, 𝑑3𝑖 is the distance to the third nearest neighboring point from node 𝑖. Radii of the subdomain 𝑟s is set to 𝛼𝑠𝑑1𝑖, 𝑑1𝑖 is the distance to the nearest neighboring point from node 𝑖. 𝛼𝑄 and 𝛼𝑠 are the scaling coefficients for determining the size of interpolation domain and subdomain, respectively.

6.1. 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||||𝜎𝑒=MLPG𝜎Exat𝜎Exat||||,(56) 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 𝐴10=0.2599 and 𝐴01=0.1608 are used. The scaling coefficients are adopted as 𝛼𝑠=0.75 and 𝛼𝑄=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 𝑞 when the polynomial term is not included (𝑚=0). Figure 3 shows that the present method obtained results have good accuracy with 𝐶 ranging from 0.5 to 3 and 𝑞=0.5, 0.5, 1.03, 1.99 when a linear polynomial basis (𝑚=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 𝐶=1.42 and 𝑞=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.

6.2. Inflation of a Long Rubber Tube

An infinitely long rubber cylinder, with outer radius 𝑟1=12.0 and inner radius 𝑟2=8.0, is subjected to an internal pressure 𝑝. 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 𝐴10=0.2599, 𝐴01=0.1608, and 𝐾=105, the shape parameters of MQ function 𝐶=1.0, 2.0, 3.0, and 𝑞=1.03, and the linear polynomial basis (𝑚=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𝑒E=||||𝑈MLPG𝑈Exat𝑈Exat||||,(57) where 𝑈MLPG and 𝑈Exat are the total strain energies computed by the present MLPG method and analytical method, respectively. The total strain energy can be expressed as𝑈=0Ω𝑊𝐼1,𝐼2,𝐽dΩ,(58) where 𝑊 is the modified Mooney-Rivlin strain energy density function.

Figure 6 gives the variation of internal pressure 𝑝 with internal radial displacement while the scaling coefficients are adopted as 𝛼𝑠=0.75, and 𝛼𝑄=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 𝛼𝑠 and interpolation domain scaling parameters 𝛼𝑄 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 𝛼𝑠 and 𝛼𝑄 affect greatly the accuracy of the present MLPG results. In Figure 7, the excellent accuracy is obtained when 𝛼𝑠 varies from 0.4 to 1.0 with 𝛼𝑄=3.5. In Figure 8, the results of the present MLPG method are acceptable only for 𝛼𝑄 ranging from 2.5 to 5.0 with 𝛼𝑠=0.8, but not for 𝛼𝑄=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 𝛼𝑄=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 𝛼𝑠=0.8 and 𝛼𝑄=3.5. In Figure 9, 𝑅 is average convergence rates; is the maximum size of nodal interval. Figure 9 shows that the present method possesses high convergence.

6.3. 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 𝐴10=1835, 𝐴01=146, and 𝐾=107 are used. The scaling coefficients are adopted as 𝛼𝑠=0.75 and 𝛼𝑄=3.5. The shape parameters of MQ function 𝐶=1.42, 𝑞=1.03 and the linear polynomial basis (𝑚=3) are selected in computation. Unless the computations fail to converge during the loading, the maximum moment 𝑀=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=𝑎,𝑏; GN is the abbreviation of the number of Gauss integration point; 𝑎 and 𝑏 are the number of Gauss integration points that adopted on each section of 0L𝑠 and 0Γ𝑠 for the numerical integration, respectively. 𝑐𝑑; 𝑐 and 𝑑 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 𝑀=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 𝑀=198.2 using irregular mesh and at 𝑀=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.

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 parameters 𝐶 and 𝑞, scaling coefficients 𝛼𝑠 and 𝛼𝑄, 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.


The financial support from National Natural Science Foundation of China (10902038) is gratefully acknowledged.