An improved meshless method for finite deformation problem in compressible hyperelastic media

φ(Ω)
σ : ∇δudV −
∫
φ(Ω)
δu · b¯dV −
∫
φ(Γt)
δu · t¯dΓ = 0, (6)
where σ is the Cauchy stress and δu is the vector of test functions.
In this study, the meshless RPIM that is developed basing on the Galerkin weak form
and discretization procedure [21,22]. The displacement field u is approximated by the use
of RPIM shape functions Φ and nodal displacement vector usp in a support domain. The
30 Nha Thanh Nguyen, Minh Ngoc Nguyen, Thien Tich Truong, Tinh Quoc Bui
Fig. 1. Hyperelastic model with boundary conditions
approximated displacement, its variational and the gradient of test functions are written
as
u = Φusp; δu = Φδusp; ∇δu = Bδusp, (7)
where B is the B-operator, this matrix contents the derivatives of shape functions with
respect to the current configuration.
Substituting Eq. (7) into Eq. (6), the weak form can be rewritten as follow∫
φ(Ω)
BTσdV −
∫
φ(Ω)
ΦTb¯dV −
∫
φ(Γt)
ΦT t¯dΓ = 0, (8)
or R (u) = Fint − Fext = 0, (9)
where R (u) denotes the residual vector that is the function of displacement field.
Apply the linearization procedure to Eq. (8) or Eq. (9), obtain the two nonlinear parts
for large deformation hyperelastic problem [23, 24]
∆R (u) =
∫
φ(Ω)
∆BTσdV +
∫
φ(Ω)
BT∆σdV = Ktan∆u. (10)
The first part in Eq. (10) relates to the stress state of the current configuration, so it
represents for the geometric nonlinear of the problem. The second one shows the varia-
tion of the stress and it stands for the material nonlinear of the model. Finally, the tangent
stiffness matrix in Eq. (10) can be written as
Ktan =
∫
φ(Ω)
GTMGdV +
∫
φ(Ω)
BT∆CeBdV, (11)
where the matrix M and matrix G at node i are defined as
M = σ ⊗ I2×2; GI =

∂φi
∂x
I2×2
∂φi
∂y
I2×2
 . (12)
where I2×2 is the (2× 2) identity matrix.
An improved meshless method for finite deformation problem in compressible hyperelastic media 31
3.2. CTM-based meshless RPIM method
In this study, the RPIM is used for the definition of meshfree shape functions. The
procedure for the construction of RPIM shap functions is described in [21, 22]. Similar
to the conventional element free Galerkin (EFG) method, the meshless RPIM approach
is also based on the global weak form. The most advantage of RPIM shape function
is its satisfaction of Kronecker delta property, this leads to convenience in treatment of
boundary conditions. However, if the traditional RPIM method is used, the integrals in
Eq. (6) or Eq. (11) are still conducted with the help of background cells. In this proposed
approach, the CTM is chosen to alternate the Gaussian integration. Suppose that a 2D
domain integral is to be calculated
I =
∫
Ω
f (x, y) dΩ, (13)
where f (x, y) is an arbitrary function. Let ΩA be the auxiliary domain that contains the
domain Ω and the new function g (x, y) is defined over the domain ΩA as
g(x, y) =
{
f (x, y), (x, y) ∈ Ω
0, (x, y) /∈ Ω (14)
Fig. 2. An auxiliary domain ΩA containing domain Ω
In this study, a rectangle is used as an auxiliary domain (see Fig. 2), the domain
integral in Eq. (13) can be rewritten as a double boundary integral
I =
∫ y2
y1
[∫ b
a
g (x, y) dx
]
dy =
∫ y2
y1
h (y) dy, (15)
where h (y) =
∫ b
a
g (x, y) dx.
To calculate the integral in Eq. (15), the boundaries Γ1 and Γ3 are divided into n inte-
gration intervals and the Gaussian quadrature method is used. The integral I in Eq. (15)
can be computed as
I =
n
∑
i=1
(∫ yi+1
yi
h(y)dy
)
=
n
∑
i=1
(∫ 1
−1
h(η)Jdη
)
=
n
∑
i=1
m
∑
j=1
h(ηj)Jwj, (16)
32 Nha Thanh Nguyen, Minh Ngoc Nguyen, Thien Tich Truong, Tinh Quoc Bui
where ηj is Gaussian points, wj is weight numbers and J = dy/dη = (yi+1 − yi)/2 de-
notes the Jacobian of the 1D transformation.
In numerical practice, to perform the integrals, a system of CTM integration points
is generated on the rectangular auxiliary domain ΩA and a simple algorithm is applied
to detect if an integration point belongs to the domain or not. In CTM-based meshless
method, there is no 2D integration cell created. This is the advantage of CTM approach
because meshing is a time-consuming task and CTM can make the RPIM method become
a truly meshless method. The detail procedure for CTM technique can be found in papers
[15–18].
4. NUMERICAL EXAMPLES
4.1. Cook’s membrane problem
The first numerical example for finite deformation analysis deals with a popular
benchmark problem, Cook’s membrane, which aims to show the accuracy of the devel-
oped CTM-based RPIM. As shown in Fig. 3(a), the model has a trapezoidal shape that is
clamped at the left edge and loaded by a traction force along the right edge. The com-
pressible neo-Hookean model is used for the hyperelastic material of the problem, the
shear modulus is µ = 80.194 N/mm2 and bulk modulus is κ = 120.291 N/mm2. Various
values of traction force are chosen for the test (i.e., F = 8, 12, 16, 24 N/mm2). There are
five nodal distributions (10× 10, 15× 15, 20× 20, 25× 25 and 30× 30) selected to inves-
tigate the accuracy of the proposed method in the nonlinear analysis. Fig. 3(b) shows the
distribution of 10× 10 scattered nodes and a system of CTM integration points is also
displayed (5196 points).
Fig. 3. Cook’s membrane problem: geometry and scattered model
The obtained results are verified with the solutions given by variational differential
quadrature (VDQ) method [19] that used 23× 23 nodes for the simulation. The present
CTM-based RPIM offers acceptable results as the computed vertical displacements at
point A (see Fig. 3) are in good agreement with reference solutions [19], which are clearly
shown in Fig. 4. The larger the traction force is, the higher value of vertical displacement
is obtained. Fig. 5 displays the first Piola Kirchhoff stresses (P = JσF−T) distributions in
An improved meshless method for finite deformation problem in compressible hyperelastic media 33
Cook’s membrane problem in cases F = 16 N/mm2 and F = 24 N/mm2, 20× 20 nodes
are used for the simulation.
Fig. 4. Vertical displacement of point A with five nodal distributions given by CTM-based RPIM,
comparing with VDQ (23× 23 nodes) [19]
Fig. 5. First Piola Kirchhoff stress distributions in Cook’s membrane problem
The computational time estimated for this example by both methods, the conven-
tional Gaussian and CTM techniques, whose computed results are given in Fig. 6, show-
ing the advantage of our development. In this comparison, several nodal densities are
used for the Cook’s membrane problem with shear force F = 8 N/mm2 that is divided in
34 Nha Thanh Nguyen, Minh Ngoc Nguyen, Thien Tich Truong, Tinh Quoc Bui
to 5 load steps. Because Gaussian-based RPIM method uses a system of background cells
for integration, so the number of integration points increases significantly when a large
number of nodes are used for the interpolation and this essentially requires more time in
the computation. As a consequence, the implementation for large deformation problems
is effective with the developed CTM-based RPIM. In other words, the computational cost
for the present model can be reduced by using the CTM technique.
Fig. 6. Computational time comparison between Gaussian-based and CTM-based RPIM methods
Tab. 1 shows the comparison of the vertical displacement of point A between
Gaussian-based and CTM-based RPIM methods. A model of 25 × 25 nodes are used
for both methods and several cases of number of integration points are investigated for
CTM approach. The data shows the stable of CTM method and the obtained results from
CTM fit well with those from conventional RPIM method although CTM-based RPIM
uses less integration points than Gaussian-based method.
Table 1. Vertical displacement of point A given by Gaussian-based and CTM-based RPIM
(25× 25 nodes are used for both methods)
q (N/mm2)
Vertical displacement of point A (mm)
Gaussian-based CTM-based
9216 points 5196 points 7112 points 8604 points
8 10.52 10.49 10.47 10.46
12 13.80 13.74 13.70 13.72
16 16.37 16.29 16.26 16.28
24 20.29 20.21 20.25 20.26
An improved meshless method for finite deformation problem in compressible hyperelastic media 35
4.2. Inhomogenous compression problem
The second numerical example is devoted to numerical simulation for a rectangular
plate under compression using our CTM-based RPIM. This is to further show the ap-
plicability and effectiveness of the developed approach in modeling large deformation
problem. The material of plate is compressible hyperelastic, the neo-Hookean model is
used and the material parameters are the same with the previous problem. It is assumed
that the horizontal displacement of the top edge and vertical displacement of the bottom
edge are set zero. A constant distributed load applies on a part of the top edge. The
geometry, dimensions and boundary conditions are presented in Fig. 7(a). Because of
the symmetry, only a half of the plate is considered in the scattered model, as shown in
Fig. 7(b).
Fig. 7. Inhomogenous compression problem: geometry and scattered model
Similar to the previous example, to investigate the accuracy and nodal independence
of the proposed method, five nodal distributions are selected. The percent of compres-
sion is computed at point I with four cases of distributed force F = 50, 100, 150 and
200 N/mm2. Graphs in Fig. 8 (a) show the compression percentages at point I that is
obtained by the proposed method and compared with those of VDQ method [19]. The
higher values of compression force give larger percents of compression. It is clearly to
see that the CTM-based RPIM can give accurate results with less degrees of freedom than
the reference method. Fig. 9 depicts the non-linear relationship between the vertical dis-
placement of point I and pressure force F, three cases of load steps are also investigated
to verify the accuracy of the proposed approach.
The charts in Fig. 10 display the comparison of computational time between the con-
ventional Gaussian and CTM techniques. The distributed force F = 200 N/mm2 with
five load steps is chosen and five cases of nodal distribution are used for the comparison.
Once again, the CTM-based RPIM shows its effectiveness in implementation for non-
linear problem. Finally, Fig. 11 represents the deformations in which the color indicates
the distribution of norm of the first Piola Kirchhoff stress
(
‖P‖ =
√
PijPij
)
on the models
with different values of pressure force.
36 Nha Thanh Nguyen, Minh Ngoc Nguyen, Thien Tich Truong, Tinh Quoc Bui
Fig. 8. Percent of compression at point I with five nodal distributions given by CTM-based RPIM,
comparing with VDQ (45× 45 nodes) [19]
Fig. 9. Force-Displacement diagrams for vertical displacement of point I
with different numbers of load steps
An improved meshless method for finite deformation problem in compressible hyperelastic media 37
Fig. 10. Computational time comparison between Gaussian-based and CTM-based RPIM
methods for inhomogenous compression problem
Fig. 11. Norm of first Piola Kirchhoff stress distributions for inhomogenous compression problem
(grey grid indicates the initial undeformed configuration of the model)
4.3. Curved beam under bending
In the last example, a model of curved beam under bending load is considered. The
geometry and boundary conditions are given in Fig. 12. The dimensions are Rin = 9 mm
Rout = 10 mm. The beam is fixed at the lowest edge while its right edge is subject to
a distributed force q with angle α = 45◦. The material is compressible hyperelastic, the
38 Nha Thanh Nguyen, Minh Ngoc Nguyen, Thien Tich Truong, Tinh Quoc Bui
properties are the same with the previous examples. Various values of distributed force
q = 0.2, 0.3, 0.4 and 0.5 N/mm2 are considered in the simulations. A distribution of
9× 101 nodes and twenty load steps are used for the problem.
Fig. 12. Curved beam under bending problem: geometry and scattered model
Fig. 13. Force - total displacement diagrams given by CTM-based RPIM,
Gaussian-based RPIM and FEM
The charts in Fig. 13 show the non-linear relationship between the external force and
total displacement at point A. The obtained results from both CTM-based and Gaussian-
based RPIM methods match well with FEM solutions given by ANSYS program. Finally,
the distributions of total displacement in the curved beam in four cases of distributed
force q are displayed in Fig. 14.
An improved meshless method for finite deformation problem in compressible hyperelastic media 39
Fig. 14. Distributions of total displacement in the curved beam
5. CONCLUSIONS
In this paper, we have developed an enhanced and effective CTM-based RPIM for
finite deformation analysis of compressible hyperelastic materials. The CTM integra-
tion technique is used as alternative to traditional Gaussian quadrature method and this
makes the meshless RPIM method become independent from the 2D integration cells.
There is no 2D meshing algorithm is required to generate the background cells. Three
numerical examples are performed and the obtained results are verified with reference
solution. The proposed method has shown good performance in both accuracy and com-
putation cost for large deformation problem. It is promising for further investigation
on behavior of hyperelastic materials such as dynamic problems and fracture analysis.
However, in this study, only regular nodal distributions are considered so the variations
of nodal density would be taken in account in the future works. Moreover, because CTM
method is being developed, there is still not any general algorithm to generate CTM inte-
gration points for every type of geometry. These drawbacks should be improved in next
studies.
ACKNOWLEDGMENT
This research is funded by Vietnam National Foundation for Science and Technology
Development (NAFOSTED) under grant number 107.02-2019.327.
40 Nha Thanh Nguyen, Minh Ngoc Nguyen, Thien Tich Truong, Tinh Quoc Bui
REFERENCES
[1] A. M. Maniatty, Y. Liu, O. Klaas, and M. S. Shephard. Higher order stabilized finite element
method for hyperelastic finite deformation. Computer Methods in Applied Mechanics and Engi-
neering, 191, (13-14), (2002), pp. 1491–1503. https://doi.org/10.1016/s0045-7825(01)00335-8.
[2] A. A. Ramabathiran and S. Gopalakrishnan. Automatic finite element formulation and as-
sembly of hyperelastic higher order structural models. Applied Mathematical Modelling, 38,
(11-12), (2014), pp. 2867–2883. https://doi.org/10.1016/j.apm.2013.11.021.
[3] A. Nomoto, H. Yasutaka, S. Oketani, and A. Matsuda. 2-dimensional homogenization
FEM analysis of hyperelastic foamed rubber. Procedia Engineering, 147, (2016), pp. 431–436.
https://doi.org/10.1016/j.proeng.2016.06.335.
[4] A. Angoshtari, M. F. Shojaei, and A. Yavari. Compatible-strain mixed finite element methods
for 2D compressible nonlinear elasticity. Computer Methods in Applied Mechanics and Engineer-
ing, 313, (2017), pp. 596–631. https://doi.org/10.1016/j.cma.2016.09.047.
[5] T. Belytschko, Y. Y. Lu, and L. Gu. Element-free Galerkin methods. Interna-
tional Journal for Numerical Methods in Engineering, 37, (2), (1994), pp. 229–256.
https://doi.org/10.1002/nme.1620370205.
[6] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl. Meshless methods: An
overview and recent developments. Computer Methods in Applied Mechanics and Engineering,
139, (1), (1996), pp. 3–47. https://doi.org/10.1016/s0045-7825(96)01078-x.
[7] N. S. Atluri and T. Zhu. A new meshless local Petrov-Galerkin (MLPG) approach
in computational mechanics. Computational Mechanics, 22, (2), (1998), pp. 117–127.
https://doi.org/10.1007/s004660050346.
[8] J. G. Wang and G. R. Liu. A point interpolation meshless method based on radial basis func-
tions. International Journal for Numerical Methods in Engineering, 54, (11), (2002), pp. 1623–1648.
https://doi.org/10.1002/nme.489.
[9] L. Gu. Moving kriging interpolation and element-free Galerkin method. Inter-
national Journal for Numerical Methods in Engineering, 56, (1), (2002), pp. 1–11.
https://doi.org/10.1002/nme.553.
[10] Y. T. Gu, Q. X. Wang, and K. Y. Lam. A meshless local Kriging method for large deformation
analyses. Computer Methods in Applied Mechanics and Engineering, 196, (9-12), (2007), pp. 1673–
1684. https://doi.org/10.1016/j.cma.2006.09.017.
[11] D. A. Hu, S. Y. Long, X. Han, and G. Y. Li. A meshless local Petrov–Galerkin method for large
deformation contact analysis of elastomers. Engineering Analysis with Boundary Elements, 31,
(7), (2007), pp. 657–666. https://doi.org/10.1016/j.enganabound.2006.11.005.
[12] D. Hu, Z. Sun, C. Liang, and X. Han. A mesh-free algorithm for dynamic impact
analysis of hyperelasticity. Acta Mechanica Solida Sinica, 26, (4), (2013), pp. 362–372.
https://doi.org/10.1016/s0894-9166(13)60033-6.
[13] Y. Zhang, W. Ge, Y. Zhang, Z. Zhao, and J. Zhang. Topology optimization of
hyperelastic structure based on a directly coupled finite element and element-
free Galerkin method. Advances in Engineering Software, 123, (2018), pp. 25–37.
https://doi.org/10.1016/j.advengsoft.2018.05.006.
[14] E. Khosrowpour, M. R. Hematiyan, and M. Hajhashemkhani. A strong-form meshfree
method for stress analysis of hyperelastic materials. Engineering Analysis with Boundary El-
ements, 109, (2019), pp. 32–42. https://doi.org/10.1016/j.enganabound.2019.09.013.
[15] A. Khosravifard and M. R. Hematiyan. A new method for meshless integration in 2D and
3D Galerkin meshfree methods. Engineering Analysis with Boundary Elements, 34, (1), (2010),
pp. 30–40. https://doi.org/10.1016/j.enganabound.2009.07.008.
An improved meshless method for finite deformation problem in compressible hyperelastic media 41
[16] A. Khosravifard, M. R. Hematiyan, and L. Marin. Nonlinear transient heat conduction anal-
ysis of functionally graded materials in the presence of heat sources using an improved
meshless radial point interpolation method. Applied Mathematical Modelling, 35, (9), (2011),
pp. 4157–4174. https://doi.org/10.1016/j.apm.2011.02.039.
[17] T. Q. Bui, A. Khosravifard, C. Zhang, M. R. Hematiyan, and M. V. Golub. Dy-
namic analysis of sandwich beams with functionally graded core using a truly mesh-
free radial point interpolation method. Engineering Structures, 47, (2013), pp. 90–104.
https://doi.org/10.1016/j.engstruct.2012.03.041.
[18] N. T. Nguyen, T. Q. Bui, M. N. Nguyen, and T. T. Truong. Meshfree thermomechanical crack
growth simulations with new numerical integration scheme. Engineering Fracture Mechanics,
(2020), p. 107121. https://doi.org/10.1016/j.engfracmech.2020.107121.
[19] R. Hassani, R. Ansari, and H. Rouhi. Large deformation analysis of 2D hyper-
elastic bodies based on the compressible nonlinear elasticity: A numerical varia-
tional method. International Journal of Non-Linear Mechanics, 116, (2019), pp. 39–54.
https://doi.org/10.1016/j.ijnonlinmec.2019.05.003.
[20] J. P. Pascon. Large deformation analysis of plane-stress hyperelastic problems via triangular
membrane finite elements. International Journal of Advanced Structural Engineering, 11, (3),
(2019), pp. 331–350. https://doi.org/10.1007/s40091-019-00234-w.
[21] G. R. Liu, G. Y. Zhang, Y. T. Gu, and Y. Y. Wang. A meshfree radial point interpolation method
(RPIM) for three-dimensional solids. Computational Mechanics, 36, (6), (2005), pp. 421–430.
https://doi.org/10.1007/s00466-005-0657-6.
[22] G. R. Liu. Mesh free methods - Moving beyon the finite element method. CRC Press, (2009).
[23] P. Wriggers. Nonlinear Finite Element Methods. Springer Berlin Heidelberg, (2008).
https://doi.org/10.1007/978-3-540-71001-1.
[24] H. D. Huynh, P. Tran, X. Zhuang, and H. Nguyen-Xuan. An extended polygonal finite el-
ement method for large deformation fracture analysis. Engineering Fracture Mechanics, 209,
(2019), pp. 344–368. https://doi.org/10.1016/j.engfracmech.2019.01.024.