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

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.

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 Ω_𝑋 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𝑊𝐼₁,𝐼₂,𝐽=𝑊𝐼₁,𝐼₂+𝑊(𝐽),𝑊𝐼₁,𝐼₂𝑊=𝐴(𝐽₁₀)=𝐼₁𝐾−3+𝐴₀₁𝐼₂−3,2(𝐽−1)²,(6) where 𝑊(𝐼₁,𝐼₂) and 𝑊(𝐽) are the deviatoric and volumetric terms of the strain energy density function, respectively. 𝐴₁₀ and 𝐴₀₁ are material constants, 𝐼₁ and 𝐼₂ are the reduced invariants, given by𝐼₁=𝐼₁𝐼^3−1/3,𝐼₂=𝐼₂𝐼^3−2/3,(7)𝐼₁, 𝐼₂, 𝐼₃ 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𝐴₁₀𝐼^3−1/3+2𝐴₀₁𝐼^3−2/3𝐼₁𝛿_𝑖𝑗−2𝐴₀₁𝐼^3−2/3𝐺_𝑖𝑗−23𝐴¹⁰𝐼^3−1/3𝐼₁+43𝐴⁰¹𝐼^3−2/3𝐼₂𝐺^−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𝐶_{𝑖𝑗𝑘𝑙}=𝐶_{𝑖𝑗𝑘𝑙}+𝐶_{𝑖𝑗𝑘𝑙}=𝜕²𝑊𝜕𝐸_𝑖𝑗𝜕𝐸_𝑘𝑙+𝜕−𝐽𝐺^−1𝑖𝑗𝜕𝐸_𝑘𝑙𝑔𝑃,_𝑖𝑗=−𝐼^31/2𝐺^−1𝑖𝑗.(11) The explicit expressions of 𝐶_{𝑖𝑗𝑘𝑙}, 𝐶_{𝑖𝑗𝑘𝑙}, and 𝑔_𝑖𝑗 are given as𝐶_{𝑖𝑗𝑘𝑙}=𝜕²𝑊𝜕𝐸_𝑖𝑗𝜕𝐸_𝑘𝑙=49𝐴¹⁰𝐼^3−1/3𝐼₁+169𝐴⁰¹𝐼^3−2/3𝐼₂𝐺^−1𝑖𝑗𝐺^−1𝑘𝑙−43𝐴¹⁰𝐼^3−1/3+83𝐴⁰¹+2𝐼^3−2/3𝐼₁𝛿^𝑖𝑗𝐺^−1𝑘𝑙+𝐺^−1𝑖𝑗𝛿_𝑘𝑙3𝐴¹⁰𝐼^3−1/3𝐼₁+43𝐴⁰¹+4𝐴₀₁𝐼^3−2/3𝐼₂𝐺^−1𝑖𝑘𝐺^−1𝑗𝑙+𝐺^−1𝑖𝑙𝐺^−1𝑗𝑘+8𝐼^3−2/3𝛿_𝑖𝑗𝛿_𝑘𝑙−2𝐴₀₁𝐼^3−2/3𝛿_𝑖𝑘𝛿_𝑗𝑙+𝛿_𝑖𝑙𝛿_𝑗𝑘3𝐴⁰¹𝐶_{𝑖𝑗𝑘𝑙}𝐼^3−2/3𝐺_𝑖𝑗𝐺^−1𝑘𝑙+𝐺^−1𝑖𝑗𝐺_𝑘𝑙,=𝜕−𝐽𝐺^−1𝑖𝑗𝜕𝐸_𝑘𝑙𝑔_𝑖𝑗𝑃=𝑃𝐽𝐺^−1𝑖𝑘𝐺^−1𝑗𝑙+𝐺^−1𝑖𝑙𝐺^−1𝑗𝑘−𝐺^−1𝑖𝑗𝐺^−1𝑘𝑙,=−𝐼^31/2𝐺^−1𝑖𝑗,(12) where 𝐺^−1𝑖𝑗 is the 𝑖𝑗th component of 𝐆⁻¹ 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𝐛𝐚^𝑇_𝑇=𝑎₁𝑎₂⋯𝑎_𝑛,𝐑_𝑇=𝑏₁𝑏₂⋯𝑏_𝑚,𝐏(𝐱)=𝑅₁(𝐱)𝑅₂(𝐱)⋯𝑅_𝑛(𝐱)_𝑇,(𝐱)=𝑝₁(𝐱)𝑝₂(𝐱)⋯𝑝_𝑚(𝐱).(14) The radial basis function 𝑅_𝑖(𝐱) has many kinds of forms and the multiquadrics (MQ) function is given by𝑅_𝑖(𝑥,𝑦)=𝑟^2𝑖+𝑐²^𝑞=𝑥−𝑥_𝑖²+𝑦−𝑦_𝑖²+𝑐²^𝑞,(15) where 𝑐,𝑞 are shape parameters.

The polynomial basis has the following monomial terms in the two-dimensional problem𝐩^𝑇(𝑥,𝑦)=1𝑥𝑦𝑥²𝑥𝑦𝑦²⋯.(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 𝐀=𝐑^𝑄𝐏_𝑚𝐑_𝑄𝐏^𝑇𝑚𝟎,=⎡⎢⎢⎢⎢⎢⎣𝑅¹𝑅₁𝑥₁,𝑦₁𝑅₂𝑥₁,𝑦₁⋯𝑅_𝑛𝑥₁,𝑦₁𝑅𝑥⋮₂,𝑦₂𝑅₂𝑥⋮₂,𝑦₂⋯⋮𝑅_𝑛𝑥⋮₂₁,𝑦₂𝑥_𝑛,𝑦_𝑛𝑅₂𝑥_𝑛,𝑦_𝑛⋯𝑅_𝑛𝑥_𝑛,𝑦_𝑛⎤⎥⎥⎥⎥⎥⎦_𝑛×𝑛𝐏_𝑚,=⎡⎢⎢⎢⎢⎢⎣𝑝¹𝑝₁𝑥₁,𝑦₁𝑝₂𝑥₁,𝑦₁⋯𝑝_𝑚𝑥₁,𝑦₁𝑝𝑥⋮₂,𝑦₂𝑝₂𝑥⋮₂,𝑦₂⋯⋮𝑝_𝑚𝑥⋮₂₁,𝑦₂𝑥_𝑛,𝑦_𝑛𝑝₂𝑥_𝑛,𝑦_𝑛⋯𝑝_𝑚𝑥_𝑛,𝑦_𝑛⎤⎥⎥⎥⎥⎥⎦_𝑛×𝑚.(20) Suppose 𝐀⁻¹ exists, (19) can be rewritten in the following form:𝐛𝐚=𝐀^−𝟏𝐮𝟎^𝐬.(21) Using (21), we have𝐛𝐚==𝑆𝑆_𝑎𝑏𝐮_𝑠,𝑆_𝑏𝐮_𝑠,𝑆_𝑎=𝐏^𝑇𝑚𝐑^𝑄−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:𝚽(𝐱)=𝜑₁(𝐱)𝜑₂(𝐱)⋅⋅⋅𝜑_𝑛(𝐱),(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𝑞𝑥−𝑥^𝑖𝑥−𝑥_𝑖²𝜕𝑅_𝑖+𝑦−𝑦_𝑖²+𝑐²^𝑞−1,𝜕𝑦=2𝑞𝑦−𝑦^𝑖𝑥−𝑥_𝑖²+𝑦−𝑦_𝑖²+𝑐²^𝑞−1.(27)

4. The MLPG Formulation of Hyperelasticity

Consider a hyperelasitc body initially occupying a domain ⁰Ω with boundary ⁰Γ. 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 𝐗∈⁰Ω 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𝐭₀d⁰Γ=𝐹_𝑗𝑘𝑆_𝑖𝑘𝑁_𝑗d⁰Γ=𝜎_𝑖𝑗𝑛_𝑗d^𝑡Γ=𝐭d^𝑡Γ,(36) where the boundary ⁰Γ is referred to the initial configuration and 𝑁_𝑗 is the initial unit outward normal to the boundary ⁰Γ.

Then (35) can be rewritten as⁰L_𝑠𝐹_𝑗𝑘𝑆_𝑖𝑘𝑁_𝑗dΓ+₀Γ_𝑠𝑢=−₀𝐹_𝑗𝑘𝑆_𝑖𝑘𝑁_𝑗dΓΓ_𝑠𝑡𝑡_0𝑖dΓ−₀Ω_𝑠𝑏_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 to⁰L_𝑠𝐹_𝑗𝑘𝑆_𝑖𝑘𝑁_𝑗dΓ+₀L_𝑠𝐹_𝑗𝑘Δ𝑆_𝑖𝑘𝑁_𝑗dΓ+₀L_𝑠𝜕Δ𝑢_𝑗𝜕𝑋_𝑘+𝑆_𝑖𝑘𝑁_𝑗dΓ₀L_𝑠𝜕Δ𝑢_𝑗𝜕𝑋_𝑘Δ𝑆_𝑖𝑘𝑁_𝑗dΓ+₀Γ_𝑠𝑢+𝐹_𝑗𝑘𝑆_𝑖𝑘𝑁_𝑗dΓ₀Γ_𝑠𝑢+𝐹_𝑗𝑘Δ𝑆_𝑖𝑘𝑁_𝑗dΓ₀Γ_𝑠𝑢𝜕Δ𝑢_𝑗𝜕𝑋_𝑘𝑆_𝑖𝑘𝑁_𝑗dΓ+₀Γ_𝑠𝑢𝜕Δ𝑢_𝑗𝜕𝑋_𝑘=−₀Δ𝑆_𝑖𝑘𝑁_𝑗dΓΓ_𝑠𝑡𝑡_0𝑖dΓ−₀Ω_𝑠𝑏_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 by12⁰L_𝑠𝐹_𝑗𝑘𝐶_{𝑖𝑘𝑚𝑙}𝐹^𝑝𝑚𝜕Δ𝑢_𝑝𝜕𝑋_𝑙+𝐹_𝑝𝑙𝜕Δ𝑢_𝑝𝜕𝑋_𝑚𝑁^𝑗+1dΓ2⁰Γ_𝑠𝑢𝐹_𝑗𝑘𝐶_{𝑖𝑘𝑚𝑙}𝐹^𝑝𝑚𝜕Δ𝑢_𝑝𝜕𝑋_𝑙+𝐹_𝑝𝑙𝜕Δ𝑢_𝑝𝜕𝑋_𝑚𝑁^𝑗+dΓ₀L_𝑠𝜕Δ𝑢_𝑗𝜕𝑋_𝑘𝑆_𝑖𝑘𝑁_𝑗dΓ+₀Γ_𝑠𝑢𝜕Δ𝑢_𝑗𝜕𝑋_𝑘+𝑆_𝑖𝑘𝑁_𝑗dΓ₀L_𝑠𝐹_𝑗𝑘𝑔_𝑖𝑘𝑁_𝑗Δ𝑃dΓ+₀Γ_𝑠𝑢=−₀𝐹_𝑗𝑘𝑔_𝑖𝑘𝑁_𝑗Δ𝑃dΓΓ_𝑠𝑡Δ𝑡_0𝑖dΓ−₀Ω_𝑠−Δ𝑏_0𝑖dΩ₀L_𝑠𝐹_𝑗𝑘𝑆_𝑖𝑘𝑁_𝑗dΓ−₀Γ_𝑠𝑢𝐹_𝑗𝑘𝑆_𝑖𝑘𝑁_𝑗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 𝜎₃=0, we can obtain [3]𝑃=𝜆₃2𝐴₁₀𝜆₃−2𝐴₀₁𝜆³⁻³=2𝐴¹⁰𝜆²¹𝜆²²−2𝐴₀₁𝜆²¹𝜆²²,(42) where 𝜆₁, 𝜆₂, and 𝜆₃ 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) Then𝑃=2𝐴₁₀𝐺₁₁𝐺₂₂−𝐺²¹²⁻¹−2𝐴₀₁𝐺₁₁𝐺₂₂−𝐺²¹².(44) Using (41), the hydrostatic pressure increment is given byΔ𝑃=𝐋Δ𝐄,(45) where𝐋^𝑇=⎧⎪⎨⎪⎩−4𝐴¹⁰𝐺₁₁𝐺₂₂−𝐺²¹²⁻²𝐺₂₂−4𝐴₀₁𝐺₂₂−4𝐴₁₀𝐺₁₁𝐺₂₂−𝐺²¹²⁻²𝐺₁₁−4𝐴₀₁𝐺₁₁8𝐴₁₀𝐺₁₁𝐺₂₂−𝐺²¹²⁻²𝐺₁₂+8𝐴₀₁𝐺₁₂Δ𝐸=Δ𝐸₁₁⎫⎪⎬⎪⎭,Δ𝐸₂₂Δ𝐸₁₂^T.(46) For the states of plane strain, we reduce the independent constraint equations in each subdomain ⁰Ω_𝑠 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𝐿₂₍₀Ω_𝑠),(47) where a constant field is used as the lower-order space, and ‖⋅‖^2𝐿₂₍₀Ω_𝑠) is the 𝐿₂ norm in the subdomain ⁰Ω_𝑠. The minimization of Θ(𝑃) leads to𝑃=−𝐾∫⁰Ω_𝑠(𝐽−1)dΩ𝐴_𝑠.(48) Similarly, in each subdomain⁰Ω_𝑠, we project the hydrostatic pressure increment onto a constant field by the least squares approximationΘ(Δ𝑃)=‖‖−𝐾𝐽𝐺^−1𝑘𝑙Δ𝐸_𝑘𝑙−Δ𝑃‖‖^2𝐿₂₍₀Ω_𝑠).(49) The minimization of Θ(Δ𝑃) leads toΔ𝑃=−𝐾∫⁰Ω_𝑠𝐠^𝑇𝐁^𝑀dΩ𝐴_𝑠Δ𝐮,(50) where 𝐴_𝑠 is the area of subdomain ⁰Ω_𝑠. 𝐁^𝑀 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𝐊_𝐼𝐽=₀L_𝑠𝐍𝐅𝐂𝐁^𝑀𝐽dΓ+₀Γ_𝑠𝑢+𝐍𝐅𝐂𝐁^𝑀𝐽dΓ₀L_𝑠𝐍𝐓𝐁^𝐺𝐽dΓ+₀Γ_𝑠𝑢𝐟_𝐼𝐍𝐓𝐁^𝐺𝐽dΓ+𝐊^{∗∗𝐼𝐽},=−₀Γ_𝑠𝑡𝐭dΓ−₀Ω_𝑠−𝐛dΩ₀L_𝑠𝐍𝐅𝐒̂dΓ−₀Γ_𝑠𝑢𝐍𝐅𝐒̂dΓ.(53) The term 𝐊^{∗∗𝐼𝐽} results from the hydrostatic pressure increment and can be expressed as𝐊^{∗∗𝐼𝐽}=⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩⁰L_𝑠𝐍𝐅𝐠𝐋𝐁^𝑀𝐽dΓ+₀Γ_𝑠𝑢−𝐾𝐴_𝑠−1forthepl𝐍an𝐅e𝐠𝐋st𝐁re^𝑀𝐽ssd,Γ⁰L_𝑠𝐍𝐅𝐠dΓ+₀Γ_𝑠𝑢×𝐍𝐅𝐠dΓ₀Ω_𝑠𝐠^T𝐁^𝑀dΩfortheplanestrain.(54) For the plane problem, the explicit expressions of vectors and matrixes in the above equations are defined by𝝋_𝑛=𝜑¹0(𝐱)𝜑₂0(𝐱)_𝑛,𝐍=𝑁¹0𝑁0₂𝑁₂0𝐂=⎡⎢⎢⎢⎣𝐶₁₁₁₁0𝑁₁,𝐶₁₁₂₂𝐶₁₁₁₂𝐶₂₂₂₂𝐶₂₂₁₂sym𝐶₁₂₁₂𝐅=⎡⎢⎢⎢⎢⎢⎣𝐹₁₁⎤⎥⎥⎥⎦,0𝐹₁₂0𝐹₂₂𝐹₂₁0𝐹₁₂𝐹₁₁𝐹₂₁0𝐹₂₂⎤⎥⎥⎥⎥⎥⎦,𝐓=⎡⎢⎢⎢⎢⎢⎣𝑆¹¹0𝑆0₂₂𝑆₁₂00𝑆₂₁𝑆0₂₁𝑆0₁₂𝑆₂₂00𝑆₁₁𝐒̂=⎧⎪⎨⎪⎩𝑆₁₁⎤⎥⎥⎥⎥⎥⎦,𝑆₂₂𝑆₁₂⎫⎪⎬⎪⎭,𝐠=⎧⎪⎨⎪⎩𝑔¹¹𝑔₂₂𝑔₁₂𝐁_𝑀𝐽⎫⎪⎬⎪⎭,=⎡⎢⎢⎢⎢⎢⎢⎣𝐹¹¹𝜕𝜑₁𝜕𝑋₁𝐹₂₁𝜕𝜑₂𝜕𝑋₁𝐹₁₂𝜕𝜑₁𝜕𝑋₂𝐹₂₂𝜕𝜑₂𝜕𝑋₂𝐹₁₁𝜕𝜑₁𝜕𝑋₂+𝐹₁₂𝜕𝜑₁𝜕𝑋₁𝐹₂₁𝜕𝜑₂𝜕𝑋₂+𝐹₂₂𝜕𝜑₂𝜕𝑋₁⎤⎥⎥⎥⎥⎥⎥⎦_𝐽𝐁_𝐺𝐽,=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝜕𝜑¹𝜕𝑋₁0𝜕𝜑0₂𝜕𝑋₂𝜕𝜑₁𝜕𝑋₂0𝜕𝜑0₂𝜕𝑋₁⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦_𝐽.(55) From (53) and (54), it can be found that the domain integral over ⁰Ω_𝑠 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 ⁰L_𝑠 and ⁰Γ_𝑠 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 𝐴₁₀=0.2599 and 𝐴₀₁=0.1608 are used. The scaling coefficients are adopted as 𝛼_𝑠=0.75 and 𝛼_𝑄=3.5.

Figure 1 

Rubber block.

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.

Figure 2 

Effect of relative error 𝑒 by shape parameter (𝑚=0).

Figure 3 

Effect of relative error 𝑒 by shape parameter (𝑚=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 radius 𝑟₁=12.0 and inner radius 𝑟₂=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 𝐴₁₀=0.2599, 𝐴₀₁=0.1608, and 𝐾=10⁵, 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𝑈=₀Ω𝑊𝐼₁,𝐼₂,𝐽dΩ,(58) where 𝑊 is the modified Mooney-Rivlin strain energy density function.

Figure 5 

Inflation of long rubber tube.

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.

Table 1 

Computational efficiency of the MLPG method influenced by the interpolation domain parameter 𝛼_𝑄(𝐶=2.0,𝑞=1.03).

Figure 6 

Internal pressure 𝑝 versus internal radial displacement.

Figure 7 

Relative errors of energy versus subdomain parameter 𝛼𝑠.

Figure 8 

Relative errors of energy versus interpolation domain parameter 𝛼𝑄.

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.

Figure 9 

Relative error and convergence Rate of energy norm.

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 𝐴₁₀=1835, 𝐴₀₁=146, and 𝐾=10⁷ 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 10 

Rubber beam under pure bending.

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 ⁰L_𝑠 and ⁰Γ_𝑠 for the numerical integration, respectively. 𝑐∗𝑑; 𝑐 and 𝑑 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 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.

Table 2 

Comparison of the MLPG method and the FEM at 𝑀=125.

Figure 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 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.

Acknowledgment

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