Large deformation analysis of elastic bodies by nonlinear Petrov–Galerkin natural element method

Abstract

A two-dimensional nonlinear Petrov–Galerkin natural element method is presented for the large deformation analysis of elastic structures. The large deformation problem is formulated according to the linearized total Lagrangian method based on Taylor series expansion. The displacement increment is approximated with the Voronoi polygon-based Laplace interpolation functions, while the admissible virtual displacement is expanded by the constant strain basis functions that are supported on Delaunay triangles. The iterative computation is performed by Newton–Raphson method, and the numerical integration is carried out by applying the conventional Gauss quadrature rule to Delaunay triangles. The proposed nonlinear Petrov–Galerkin natural element method is illustrated through the numerical experiments, and its numerical accuracy is compared with MSC/Marc, constant strain finite element method, and Bubnov–Galerkin natural element method and the symmetry preservation in the very large strain is presented.

Introduction

The finite element method (FEM) has been widely and successfully employed during several decades to numerically solve various engineering problems using nonoverlapping finite elements that can systematically and effectively discretize a problem domain. But, at the same time, it possesses the inherent drawbacks such as the element connectivity preservation, the numerical quality deterioration stemming from the excessive mesh distortion and the painstaking mesh adaptation. In this context, several kinds of mesh-free methods were introduced^1–4 to resolve the drawbacks of FEM by substituting finite elements with grid points. These mesh-free methods provide highly smooth approximation solutions, but the subdomains are not perfectly separated and overlapped so that their basis functions do not obey the Kronecker delta property. Thus, additional techniques are required to impose the essential boundary conditions.^5,6 Furthermore, these mesh-free methods need to generate additional background cells for the numerical integration using the conventional Gauss quadrature rule. What is worse, the numerical integration accuracy is influenced by how to construct the background cell.

Differing from these grid-point-based mesh-free methods, the natural element method (NEM) introduced by Braun and Sambridge ⁷ employs the basis functions called Laplace interpolation functions that are defined in terms of Voronoi polygons and Delaunay triangles. Since Laplace interpolation functions are defined based on the geometric relationship of Voronoi polygons, those strictly obey the Kronecker delta property.^8,9 Hence, the troublesome treatment of essential boundary conditions in the conventional mesh-free methods can be completely resolved. Another important advantage of NEM is no extra effort is needed to generate background cell for the numerical integration, because Delaunay triangles that were generated for defining Laplace interpolation functions automatically serve as a background cell. In addition, like the basis functions in conventional mesh-free methods, Laplace interpolation functions exhibit the sufficiently smooth $C^1$ -interpolant characteristic. ¹⁰

Thanks to the abovementioned merits, the NEM has been rapidly extended to solve linear and nonlinear problems in various engineering fields.^11–16 Meanwhile, the NEM may be classified into Bubanov–Galerkin and Petrov–Galerkin, depending on whether the trial and test basis functions are identical or not. In the former case, the intersection of trial and test basis function supports does not coincide with a union of Delaunay triangles, which gives rise to the numerical integration inaccuracy.^17,18 This problem owing to the inconsistency between the integrand function support and the background cell can be effectively resolved by employing the test basis function that is different from Laplace interpolation function. Furthermore, in Petrov–Galerkin, the uniform Delaunay triangularization always provides the symmetric integration meshes. This symmetry at the integration mesh level guarantees the preservation of deformation symmetry, even at very large strain. Wang et al. ¹⁹ introduced a sort of Petrov–Galerkin natural element method (PG-NEM) in which the support of test basis function is identical with a Delaunay triangle. The PG-NEM leads to nonsymmetric matrices, which is an indispensable drawback to secure the numerical integration accuracy. However, nowadays it becomes a less critical issue thanks to the advances in computational facilities and technologies.

In this context, the purpose of this article is to extend a two-dimensional PG-NEM to the geometrically large deformation analysis of elastic structures. The geometric nonlinear problem is formulated according to the linearized total Lagrangian method based on the Taylor series expansion. The displacement increment is approximated by Laplace interpolation functions, while the virtual displacement is expanded by the constant strain (CS) basis functions ²⁰ that are supported on Delaunay triangles. The numerical integration is performed by applying the conventional Gauss quadrature rule to Delaunay triangles, and the iterative computation is carried out by Newton–Raphson method. The proposed nonlinear PG-NEM is illustrated through the numerical experiments, and its numerical accuracy is verified through the comparison with MSC/Marc, constant strain FEM (CS-FEM), and Bubnov–Galerkin natural element method (BG-NEM). As well, the symmetry preservation in the deformed configuration is presented at very large strain.

Laplace interpolation functions

Figure 1(a) represents the Voronoi diagram and the Delaunay triangulation, which lay down a geometric framework for constructing a natural element grid and defining Laplace interpolation functions. For a given set $\aleph$ of N distinct points $x_I$ called nodes: $\aleph =\{x_1,x_2,\dots ,x_N|x_I\in {\Re }^2\}$ in ${\Re }^2$ , a first-order Voronoi diagram $D_{\mathrm{Vo}}$ is defined with N Voronoi polygons ${\varpi }_I$ such that $D_{\mathrm{Vo}}=\{{\varpi }_1,{\varpi }_2,\dots ,{\varpi }_N|\overline{{\varpi }_1\cup ···\cup {\varpi }_N}={\Re }^2\}$ . Here, the Voronoi polygon ${\varpi }_I$ corresponding to the Ith node is defined by ${\varpi }_I=\{x\in {\Re }^2:d(x,x_I)<d(x,x_J),\forall J\ne I\}$ with $d(x,x_I)$ being the Euclidean metric in ${\Re }^2$ , which is a subdomain of the unbounded domain $\varpi \in {\Re }^2$ with its sides, which perpendicularly bisect the lines connecting $x_I$ and the adjacent neighbor nodes in $\aleph$ . The Voronoi polygons are uniquely defined to each node with the vertex points called Voronoi vertices, and those become to be unbounded when their nodes are located on the boundary of the convex hull ${\varpi }_{\mathrm{CH}}(\aleph )$ of a set $\aleph$ . The convex hull becomes the problem domain: ${\varpi }_{\mathrm{CH}}(\aleph )\equiv \mathrm{\Omega }\subset \varpi$ in case of convex domain, but does not for nonconvex domain.

Figure 1.

(a) Voronoi diagram and Delaunay triangulation. (b) Second-order Voronoi cells.

As a geometric dual of the Voronoi diagram, the Delaunay triangulation generates a set ${\Im }_{\mathrm{De}}$ of Delaunay triangles, ${\mathrm{\Delta }}_{\mathrm{JKL}}(x)$ , ${\Im }_{\mathrm{De}}=\{{\mathrm{\Delta }}_{\mathrm{JKL}}(x):\bigcup \overline{{\mathrm{\Delta }}_{\mathrm{JKL}}}={\varpi }_{\mathrm{CH}}(\aleph ),J\ne K\ne L\}$ , where J, K, and L refer to three distinct vertex nodes $x_J,x_K$ , and $x_L$ of ${\mathrm{\Delta }}_{\mathrm{JKL}}$ . In other words, the problem domain $\mathrm{\Omega }$ is discretized with a finite number of nonoverlapping Delaunay triangles. These three vertex nodes should be chosen such that their Voronoi polygons share common edges, and a circumcircle defined by these three nodes is called the Delaunay circumcircle $\mathrm{cir}({\mathrm{\Delta }}_{\mathrm{JKL}})$ . An important property of the Delaunay triangulation is the empty circumcircle criterion implying that each Delaunay circumcircle should not contain any node $x_I$ in $\aleph$ . It is not hard to realize that the centers of these circumcircles are identical to the Voronoi vertices.

To each node $x_I\in \aleph$ , exactly one Laplace interpolation function ${\phi }_I(x)$ is defined, where the subscript I designates the I th node in a natural element grid. In order to define the I th Laplace interpolation function ${\phi }_I(x)$ , let us consider five Delaunay circumcircles as shown in Figure 1(b), which are generated by node $x_I$ and its five neighbor nodes. By adding a point $x_P$ into the I th Voronoi polygon ${\varpi }_I$ , one can construct additional Voronoi polygon ${\varpi }_P$ composed of four subregions, which are divided by the previously defined four Voronoi polygons ${\varpi }_I,{\varpi }_J,{\varpi }_K$ , and ${\varpi }_L$ . Here, four subregions are called the second-order Voronoi cells ${\varpi }_{\mathrm{PI}}$ , and those are mathematically defined by ${\varpi }_{\mathrm{PI}}=\{x\in {\Re }^2:d(x,x_P)<d(x,x_I)<d(x,x_J),\forall J\ne P,I\}$ . The number of ${\varpi }_{\mathrm{PI}}$ is not identical with the number of nodes, furthermore the shapes of ${\varpi }_P$ and ${\varpi }_{\mathrm{PI}}$ vary depending on the location of $x_P$ .

With the edge lengths $s_{\ell }=\left|{\mathrm{\Gamma }}_{P\ell }\right|$ of ${\varpi }_P$ and the relative half distances $h_{\ell }=d(x_P,x_{\ell })/2$ between $x_P$ and the neighbor nodes, the weighting factors ${\alpha }_{\ell }$ are defined by

$${\alpha }_{\ell }\left(x_P\right)=\frac{s_{\ell }\left(x_P\right)}{h_{\ell }\left(x_P\right)},\ell =I,J,K,L$$

(1)

Then, the value of Laplace interpolation function ${\phi }_I$ and its derivatives ${\phi }_{I,\beta }(\beta =x,y)$ at a point $x_P$ are calculated by Chinesta et al ⁸

$$\begin{array}{c}{\phi }_I\left(x_P\right)={\alpha }_I\left(x_P\right)/\sum _{\ell =1}^M{\alpha }_{\ell }\left(x_P\right),{\phi }_{I,\beta }\left(x_P\right)\hfill \\ =[{\alpha }_{I,\beta }\left(x_P\right)-{\phi }_I\left(x_P\right)\sum _{\ell =1}^M{\alpha }_{\ell ,\beta }\left(x_P\right)]/\sum _{\ell =1}^M{\alpha }_{\ell }\left(x_P\right)\hfill \end{array}$$

(2)

with ${\alpha }_{I,\beta }(x_P)=[s_{I,\beta }(x_P)-{\alpha }_I{h}_{I,\beta }(x_P)]/h_I(x_P)$ . And, the Laplace interpolation function ${\phi }_I$ corresponding to the I th node $x_I$ is defined by replacing $x_P$ with an arbitrary point x within the support $\mathrm{supp}({\phi }_I)$ shown in Figure 2(a). The supports and the shapes of Laplace interpolation functions corresponding to three representative nodes are given in Figure 2. Here, the support $\mathrm{supp}({\phi }_I)$ of ${\phi }_I$ is defined by the intersection of the convex hull ${\varpi }_{\mathrm{CH}}\left(\aleph \right)$ and the union $\cup \overline{\mathrm{cir}({\mathrm{\Delta }}_{\mathrm{IJK}})}$ of Delaunay circumcircles defined by the node $x_I$ and its neighbor nodes, $\mathrm{supp}({\phi }_{\mathrm{I}}(x))=\bigcup \overline{\mathrm{cir}({\mathrm{\Delta }}_{\mathrm{IJK}}(x))}\cap {\omega }_{\mathrm{CH}}(\aleph )$ . The reader may refer to literature^8–10 for the detailed properties of Laplace interpolation functions.

Figure 2.

Laplace interpolation functions: (a) Supports. (b) Shapes.

Linearized total Lagrangian formulation

Referring to Figure 3, let us consider the geometrically nonlinear problem of a two-dimensional large deformable elastic body occupying initially an open bounded domain, ${\mathrm{\Omega }}^0\in {\Re }^2$ . And, let ${\mathrm{\Omega }}^{t+\mathrm{\Delta }t}$ be the deformed equilibrium with its boundary ${\mathrm{\Gamma }}^{t+\mathrm{\Delta }t}$ at the current time, $t+\mathrm{\Delta }t$ . Note that the time variable does not represent real time, but the parameter represents the load step. We assume that there exists a unique smooth mapping function ${\xi }_{t+\mathrm{\Delta }t}(x^0)$ from ${\mathrm{\Omega }}^0$ to ${\mathrm{\Omega }}^{t+\mathrm{\Delta }t}$ such that

$$x^{t+\mathrm{\Delta }t}={\xi }_{t+\mathrm{\Delta }t}\left(x^0\right)=x^0+u^{t+\mathrm{\Delta }t}$$

(3)

between the initial and current coordinates $x^0$ and $x^{t+\mathrm{\Delta }t}$ . Then, the deformation gradient is defined by

$$F_{\mathrm{ij}}^{t+\mathrm{\Delta }t}=\frac{\partial {\xi }_{t+\mathrm{\Delta }t}\left(x_i^0\right)}{\partial x_j^0}={\delta }_{\mathrm{ij}}+\frac{\partial u_i^{t+\mathrm{\Delta }t}}{\partial x_j^0}$$

(4)

Figure 3.

Initial and deformed configurations of a large deformable elastic body.

The deformed configuration of the body at the current time $t+\mathrm{\Delta }t$ is governed by the static equilibrium given by

$$\frac{\partial {\sigma }_{\mathrm{ij}}^{t+\mathrm{\Delta }t}}{\partial x_j^{t+\mathrm{\Delta }t}}+b_i^{t+\mathrm{\Delta }t}=0\mathrm{in}{\mathrm{\Omega }}^{t+\mathrm{\Delta }t}$$

(5)

with essential and traction boundary conditions

$$u_i^{t+\mathrm{\Delta }t}={\hat{u}}_{\mathrm{ij}}^{t+\mathrm{\Delta }t}\mathrm{on}{\mathrm{\Gamma }}_D^{t+\mathrm{\Delta }t}$$

(6)

$${\sigma }_{\mathrm{ij}}^{t+\mathrm{\Delta }t}{n}_j^{t+\mathrm{\Delta }t}={\hat{t}}_i^{t+\mathrm{\Delta }t}\mathrm{on}{\mathrm{\Gamma }}_N^{t+\mathrm{\Delta }t}$$

(7)

The undeformed initial configuration is taken as a reference frame using the second Piola–Kirchhoff stress $S_{\mathrm{ij}}^{t+\mathrm{\Delta }t}$ and the Jacobian $J^{t+\mathrm{\Delta }t}(=\left|F_{\mathrm{ij}}^{t+\mathrm{\Delta }t}\right|)$ of the deformation gradient. Then, the total Lagrangian formulation of the above boundary value problem in equations (5) to (7) becomes: Find $u^{t+\mathrm{\Delta }t}$ such that

$$\begin{array}{c}\underset{{\mathrm{\Omega }}^0}{\int }{S}_{\mathrm{ij}}^{t+\mathrm{\Delta }t}\left(u^{t+\mathrm{\Delta }t}\right){\tilde{E}}_{\mathrm{ij}}^{t+\mathrm{\Delta }t}(u^{t+\mathrm{\Delta }t},v)d{\mathrm{\Omega }}^0\hfill \\ =\underset{{\mathrm{\Gamma }}_N^0}{\int }{\hat{t}}_i^0{v}_{i}d{\mathrm{\Gamma }}^0+\underset{{\mathrm{\Omega }}^0}{\int }{b}_i^0{v}_{i}d{\mathrm{\Omega }}^0\hfill \end{array}$$

(8)

for every virtual displacement v.

The strain-like tensor ${\tilde{E}}_{\mathrm{ij}}^{t+\mathrm{\Delta }t}$ is defined by

$$\begin{array}{c}{\tilde{E}}_{\mathrm{ij}}^{t+\mathrm{\Delta }t}(u^{t+\mathrm{\Delta }t},v)\hfill \\ =\frac{1}{2}[\frac{\partial v_i}{\partial x_j^0}+\frac{\partial v_j}{\partial x_i^0}+\frac{\partial u_m^{t+\mathrm{\Delta }t}}{\partial x_i^0}\frac{\partial v_m}{\partial x_j^0}+\frac{\partial v_m}{\partial x_i^0}\frac{\partial u_m^{t+\mathrm{\Delta }t}}{\partial x_j^0}]\hfill \end{array}$$

(9)

in terms of test and trial functions v and u, and it represents the variation of the Green–Lagrange (GL) strain $E_{\mathrm{ij}}^{t+\mathrm{\Delta }t}$ defined by

$$\begin{array}{cc}\multicolumn{1}{c}{E_{\mathrm{ij}}^{t+\mathrm{\Delta }t}}& =\frac{1}{2}[F_{\mathrm{mi}}^{t+\mathrm{\Delta }t}{F}_{\mathrm{mj}}^{t+\mathrm{\Delta }t}-{\delta }_{\mathrm{ij}}]\hfill \\ \multicolumn{1}{c}{}& =\frac{1}{2}[\frac{\partial u_i^{t+\mathrm{\Delta }t}}{\partial x_j^0}+\frac{\partial u_j^{t+\mathrm{\Delta }t}}{\partial x_i^0}+\frac{\partial u_m^{t+\mathrm{\Delta }t}}{\partial x_i^0}\frac{\partial u_m^{t+\mathrm{\Delta }t}}{\partial x_j^0}]\hfill \end{array}$$

(10)

The weak form equation (8) is nonlinear because $E_{\mathrm{ij}}^{t+\mathrm{\Delta }t}$ includes the displacement $u^{t+\mathrm{\Delta }t}$ to be solved, so it is needed to be linearized according to the Taylor series expansion for the iterative numerical computation. The first term in the left-hand side of equation (8) is linearized at the kth equilibrium stage such that

$$S_{\mathrm{ij}}^{t+\mathrm{\Delta }t}{\tilde{E}}_{\mathrm{ij}}^{t+\mathrm{\Delta }t}\approx S_{\mathrm{ij}}^{\left(k\right)}{\tilde{E}}_{\mathrm{ij}}^{\left(k\right)}+\frac{\partial }{\partial u_m^{\left(k\right)}}\left(S_{\mathrm{ij}}^{\left(k\right)}{\tilde{E}}_{\mathrm{ij}}^{\left(k\right)}\right)u_m^{{\mathrm{\Delta }}^{(k+1)}}$$

(11)

Referring to Figure 4, two superscripts $(k)$ and ${\mathrm{\Delta }}^{(k+1)}$ indicate the values at the kth iteration stage and the increments in the (k + 1)-th iteration stage, respectively. In addition, the second term in the right-hand side of equation (11) is to be expanded as

$$\frac{\partial }{\partial u_m^{\left(k\right)}}\left(S_{\mathrm{ij}}^{\left(k\right)}{\tilde{E}}_{\mathrm{ij}}^{\left(k\right)}\right)u_m^{{\mathrm{\Delta }}^{(k+1)}}=C_{\mathrm{ijpq}}^{\left(k\right)}{E}_{\mathrm{pq}}^{{\mathrm{\Delta }}^{(k+1)}}{\tilde{E}}_{\mathrm{ij}}^{\left(k\right)}+S_{\mathrm{ij}}^{\left(k\right)}{\tilde{E}}_{\mathrm{ij}}^{{\mathrm{\Delta }}^{(k+1)}}$$

(12)

with $C^{(k)}$ being a fourth-order elastic modulus tensor at the kth equilibrium stage. And, $E^{{\mathrm{\Delta }}^{(k+1)}}$ and ${\tilde{E}}^{{\mathrm{\Delta }}^{(k+1)}}$ are interpreted as the increments of the Green–Lagrangian strain tensor and its variation at the (k + 1)th stage, respectively

$$\begin{array}{c}{E}_{\mathrm{ij}}^{{\mathrm{\Delta }}^{(k+1)}}\left(u^{{\mathrm{\Delta }}^{(k+1)}}\right)\hfill \\ =\frac{1}{2}[\frac{\partial u_i^{{\mathrm{\Delta }}^{(k+1)}}}{\partial x_j^0}+\frac{\partial u_j^{{\mathrm{\Delta }}^{(k+1)}}}{\partial x_i^0}+\frac{\partial u_m^{{\mathrm{\Delta }}^{(k+1)}}}{\partial x_i^0}\frac{\partial u_m^{\left(k\right)}}{\partial x_j^0}+\frac{\partial u_m^{\left(k\right)}}{\partial x_i^0}\frac{\partial u_m^{{\mathrm{\Delta }}^{(k+1)}}}{\partial x_j^0}]\hfill \end{array}$$

(13)

$${\tilde{E}}_{\mathrm{ij}}^{{\mathrm{\Delta }}^{(k+1)}}(u^{{\mathrm{\Delta }}^{(k+1)}},v)=\frac{1}{2}[\frac{\partial v_m}{\partial x_i^0}\frac{\partial u_m^{{\mathrm{\Delta }}^{(k+1)}}}{\partial x_j^0}+\frac{\partial u_m^{{\mathrm{\Delta }}^{(k+1)}}}{\partial x_i^0}\frac{\partial v_m}{\partial x_j^0}]$$

(14)

Figure 4.

Schematic representation of the iterative computation scheme.

Meanwhile, according to the Taylor series expansion, the surface traction and the body force are linearized as

$$t_i^{\left(k\right)}=A^{\left(k\right)}{t}_i^{t+\mathrm{\Delta }t}\left(u^{\left(k\right)}\right),b_i^{\left(k\right)}=J^{\left(k\right)}{b}_i^{t+\mathrm{\Delta }t}\left(u^{\left(k\right)}\right)$$

(15)

with $A^{(k)}=J^{(k)}\left|{[F_{\mathrm{ij}}^{(k)}]}^{-1}{n}_j^0\right|$ . Then, finally the linearization of weak form equation (8) ends up with

$$\begin{array}{c}\underset{{\mathrm{\Omega }}^0}{\int }{C}_{\mathrm{ijpq}}^{\left(k\right)}{E}_{\mathrm{pq}}^{{\mathrm{\Delta }}^{(k+1)}}{\tilde{E}}_{\mathrm{ij}}^{\left(k\right)}\mathrm{d}{\mathrm{\Omega }}^0+\underset{{\mathrm{\Omega }}^0}{\int }{S}_{\mathrm{ij}}^{\left(k\right)}{\tilde{E}}_{\mathrm{ij}}^{{\mathrm{\Delta }}^{(k+1)}}\mathrm{d}{\mathrm{\Omega }}^0=\underset{{\mathrm{\Gamma }}_N^0}{\int }{t}_i^{\left(k\right)}{v}_i\mathrm{d}{\mathrm{\Gamma }}^0\hfill \\ +\underset{{\mathrm{\Omega }}^0}{\int }{b}_i^{\left(k\right)}{v}_i\mathrm{d}{\mathrm{\Omega }}^0-\underset{{\mathrm{\Omega }}^0}{\int }{S}_{\mathrm{ij}}^{\left(k\right)}{\tilde{E}}_{\mathrm{ij}}^{\left(k\right)}\mathrm{d}{\mathrm{\Omega }}^0\hfill \end{array}$$

(16)

Elastic structures are modeled as a Kirchhoff material, which is path-independent and possess an elastic strain energy potential. The strain energy density functional and the elastic modulus tensor are as follows ²¹

$$\begin{array}{c}{W}^{\left(k\right)}=\frac{1}{2}{C}_{\mathrm{ijpq}}^{\left(k\right)}{E}_{\mathrm{pq}}^{\left(k\right)}{E}_{\mathrm{ij}}^{\left(k\right)},\hfill \\ C_{\mathrm{ijpq}}^{\left(k\right)}=C_{\mathrm{ijpq}}=\lambda {\delta }_{\mathrm{ij}}{\delta }_{\mathrm{pq}}+\mu ({\delta }_{\mathrm{ip}}{\delta }_{\mathrm{jq}}+{\delta }_{\mathrm{iq}}{\delta }_{\mathrm{jp}})\hfill \end{array}$$

(17)

Two-dimensional nonlinear Petrov–Galerkin natural element approximation

For the PG-NE approximation of the linearized weak form equation (16), the displacement increment $u^{{\mathrm{\Delta }}^{(k+1)}}$ at the (k + 1)th stage and the test displacement v are expanded as

$$u_i^{{\mathrm{\Delta }}^{(k+1)}}\left(x^0\right)=\sum _{I=1}^N{\phi }_I\left(x^0\right){\overline{u}}_{\mathrm{Ii}}^{{\mathrm{\Delta }}^{(k+1)}}$$

(18)

$$v_i\left(x^0\right)=\sum _{J=1}^N{\psi }_J\left(x^0\right){\overline{v}}_{\mathrm{Ji}}$$

(19)

Similarly, the displacements at the kth stage (i.e. at time $t+\mathit{\Delta }t$ ) and at time stage t are expressed by

$$u_i^{\left(k\right)}\left(x^0\right)=\sum _{I=1}^N{\phi }_I\left(x^0\right){\overline{u}}_{\mathrm{Ii}}^{\left(k\right)}$$

(20)

$$u_i^t\left(x^0\right)=\sum _{I=1}^N{\phi }_I\left(x^0\right){\overline{u}}_{\mathrm{Ii}}^t$$

(21)

for a given natural element grid composed of N nodes. Here, ${\varphi }_I$ are Laplace interpolation functions, while ${\psi }_J$ are the constant-strain basis functions that are supported on a union of Delaunay triangles, as shown in Figure 5(a).

Figure 5.

Intersection ${\mathrm{\Omega }}_{\mathrm{IJ}}^{\mathrm{int}}$ between trial function and test function: (a) ${\phi }_I$ and ${\phi }_J$ in PG-NEM. (b) ${\psi }_I$ and ${\phi }_J$ in BG-NEM.

Three shape functions ${\hat{\psi }}_i(x^0)(i=1,2,3)$ defined within a Delaunay triangle with three vertex points $x_1^0,x_2^0$ and $x_3^0$ and the area A are given by

$${\hat{\psi }}_i\left(x^0\right)=\frac{\{(x_j^0{y}_k^0-x_k^0{y}_j^0)+(y_j^0-y_k^0)x^0+(x_k^0-x_j^0)y^0\}}{2A}$$

(22)

The approximation becomes the BG-NEM when both the trial and test displacements are expanded by Laplace interpolation function, ${\phi }_I$ . Meanwhile, the finite element approximation using the constant-strain basis functions as for both trial and test functions is called by CS-FEM in the current study.

Substituting equations (18) to (21) into equation (16) leads to the matrix equations for geometrically nonlinear problems, which are given by

$$\sum _I^N[{{}^{m}K}_I^{\left(k\right)}+{{}^{g}K}_I^{\left(k\right)}]{\overline{u}}^{{\mathrm{\Delta }}^{(k+1)}}=\sum _I^N\{{{}^{e}F}_I^{\left(k\right)}-{{}^{i}F}_I^{\left(k\right)}\}$$

(23)

Here, ${{}^m\mathbf{K}}^{(k)}$ and ${{}^g\mathbf{K}}^{(k)}$ are the $(2N\times 2N)$ material tangent stiffness matrix and geometric tangent stiffness matrix at the kth iteration step at time $t+\mathrm{\Delta }t$ , respectively. ${{}^e\mathbf{F}}^{(k)}$ and ${{}^i\mathbf{F}}^{(k)}$ represent the $(2N\times 1)$ external force vector and internal force vector at the kth iteration step. These node-wise matrices and vectors are defined by

$${{}^m\mathbf{K}}_I^{\left(k\right)}=\underset{{\mathrm{\Omega }}_v^I}{\int }{\left(B_v^{\left(k\right)}\right)}^T\overline{C}{B}_u^{\left(k\right)}\mathrm{d}{\mathrm{\Omega }}^0$$

(24)

$${{}^g\mathbf{K}}_I^{\left(k\right)}=\underset{{\mathrm{\Omega }}_v^I}{\int }{\left(H_v^{\left(k\right)}\right)}^T{\hat{S}}^{\left(k\right)}{H}_u^{\left(k\right)}\mathrm{d}{\mathrm{\Omega }}^0$$

(25)

$${{}^e\mathbf{F}}_I^{\left(k\right)}=\underset{{\mathrm{\Omega }}_v^I}{\int }{\mathrm{\Psi }}^T{b}^{\left(k\right)}\mathrm{d}{\mathrm{\Omega }}^0+\underset{{\mathrm{\Gamma }}_N^0\cap {\mathrm{\Omega }}_v^I}{\int }{\mathrm{\Psi }}^T{t}^{\left(k\right)}\mathrm{d}{\mathrm{\Gamma }}^0$$

(26)

$${{}^i\mathbf{F}}_I^{\left(k\right)}=\underset{{\mathrm{\Omega }}_v^I}{\int }{\left({\tilde{B}}_v^{\left(k\right)}\right)}^T{\hat{S}}^{\left(k\right)}\mathrm{d}{\mathrm{\Omega }}^0$$

(27)

with ${\mathrm{\Omega }}_v^I={\mathrm{\Omega }}^0\cap \mathrm{supp}({\psi }_{\mathrm{I}})$ , ${\stackrel{·}{S}}^{(k)}={\{S_{11}^{(k)},S_{22}^{(k)},S_{12}^{(k)},S_{33}^{(k)}\}}^T$ , the $(3\times 3)$ material modulus tensor ${\overline{C}}^{\left(k\right)}$ , the $(2\times 2N)$ basis function matrices $\Phi$ and $\mathrm{\Psi }$ , ${\mathbf{B}}_{\mathbf{v}}^{(k)}=[{\mathbf{B}}_v^1\dots {\mathbf{B}}_v^I\dots {\mathbf{B}}_v^N]$ , ${\mathbf{B}}_u^{(k)}=[{\mathbf{B}}_u^1\dots {\mathbf{B}}_u^J\dots {\mathbf{B}}_u^N]$ , ${\mathbf{H}}_{\mathbf{v}}^{(k)}=[{\mathbf{H}}_v^1\dots {\mathbf{H}}_v^I\dots {\mathbf{H}}_v^N]$ , and ${\mathbf{H}}_u^{(k)}=[{\mathbf{H}}_u^1\dots {\mathbf{H}}_u^J\dots {\mathbf{H}}_u^N]$ . The submatrices included in equations (24) to (27) are defined in Appendix 1. These node-wise matrices and load vectors are to be integrated within individual supports $\mathrm{supp}({\psi }_I)$ of the test basis functions using the conventional Gauss quadrature rule. Differing from the inconsistency shown in Figure 5(b) for BG-NEM, not only the support of ${\psi }_I$ is identical with a union of Delaunay triangles, but also the intersection ${\mathrm{\Omega }}_{\mathrm{IJ}}^{\mathrm{int}}$ between the supports of ${\phi }_J$ and ${\psi }_I$ is always contained within $\mathrm{supp}({\psi }_I)$ . Thus, the numerical integration of these matrices and load vectors can be accurately and easily performed as in the FEM.

Figure 6 represents a flowchart for the large deformation analysis by the current PG-NEM, which is composed of the outer load increment loop and the inner nonlinear iteration loop. The inner iterative computation is performed by the Newton–Raphson method until the following two convergence criteria are satisfied

$$\mathrm{\Delta }{\overline{u}}_{\mathrm{rel}}=\frac{\left|{\overline{u}}^{{\mathrm{\Delta }}^{(k+1)}}\right|}{\left|{\overline{u}}^{(k+1)}\right|}\le {\alpha }_u\left|{\overline{u}}^{(k+1)}\right|$$

(28)

$$\mathrm{\Delta }{\overline{f}}_{\mathrm{rel}}=\frac{|{{}^e\mathbf{F}}^{(k+1)}-{{}^i\mathbf{F}}^{(k+1)}|}{\left|{{}^e\mathbf{F}}^{(k+1)}\right|}\le {\alpha }_f\left|{{}^e\mathbf{F}}^{(k+1)}\right|$$

(29)

where, $\mathrm{\Delta }{\overline{u}}_{\mathrm{rel}}$ and $\mathrm{\Delta }{\overline{f}}_{\mathrm{rel}}$ mean the magnitudes of the relative displacement increment and the relative force residual, and ${\alpha }_u$ and ${\alpha }_f$ are set by $1\times {10}^{-4}$ for the current study. When the convergence is satisfied, then the structure geometry and Delaunay triangles are updated and the applied load is increased. And, the convergence iteration restarts and finally the analysis is completed successfully when the load increment reaches the applied load. On the contrary, if the convergence iteration number is larger than the preset maximum iteration number, $k^{max}$ , then the control parameter $\theta$ of load increment is reduced by half and the convergence iteration is repeated. Still the convergence is not satisfied until the load increment control parameter $\theta$ becomes smaller than the preset minimum value ${\theta }_{min}$ , then the analysis is terminated unsuccessfully, and the analysis should be retried by adjusting all the input data.

Figure 6.

Flowchart of the large deformation analysis.

Numerical experiments

A test program of nonlinear PG-NEM was coded according to the numerical formula given in the previous sections, and its validity and accuracy were illustrated and investigated through four representative large deformation problems. Figure 7(a) shows a cantilever beam subject to an upward vertical tip load P, where the beam is assumed to be in the plane stress condition. The length L is $20\mathrm{mm}$ and the section dimensions are set by $t=h=1.0\mathrm{mm}$ , while the Young’s modulus E and the Poisson’s ratio are set by $400\mathrm{MPa}$ and 0, respectively. To examine the effect of NEM grid density, two uniform NEM grids are taken. The nonlinear analysis was carried out by increasing the vertical tip load P up to $400\mathrm{MPa}$ in 10 uniform load steps. The analytic solution of this problem is available in a book by Doyle. ²²

Figure 7.

(a) A cantilever beam-like structure. (b) Final deformed configurations obtained by three different methods.

The problem was also solved by BG-NEM and CS-FEM for the comparison purpose, and the total number of finite element nodes is kept the same for CS-FEM. The numerical results obtained by PG-NEM, BG-NEM, and CS-FEM are compared with the analytic solution. Figure 7(b) compares the final deformed shapes of cantilever beam, where PG-NEM and BG-NEM show almost the similar shape. But, CS-FEM leads to the remarkably different shape such that the upward deformation is much lower than that of PG-NEM. It is because the CS basis function used for CS-FEM leads to the stiffness matrix that is relatively stiffer, when compared with Laplace interpolation function. Figures 8 and 9 comparatively represent the load-deflection curves for two different NEM grids. It is clearly observed that both NEMs using Laplace interpolation function provide much better result than CS-FEM using the CS basis function. When compared with CS-FEM, NEMs provide relatively accurate results even at coarse NEM grids. Meanwhile, when compared with PG-NEM, BG-NEM leads to the slightly poor result, in particular at coarse grid. It demonstrates that the numerical integration inaccuracy of BG-NEM does also deteriorate the numerical accuracy of the large deformation analysis. However, one can find from Figure 9 that the effect of numerical integration inaccuracy becomes insignificant as the grid density increases.

Figure 8.

Comparison of the normalized deflections (33 nodes): (a) Horizontal $({\delta }_h/L)$ . (b) Vertical $({\delta }_v/L)$ .

Figure 9.

Comparison of the normalized deflections (105 nodes): (a) Horizontal $({\delta }_h/L)$ . (b) Vertical $({\delta }_v/L)$ .

The next example is an arch shown in Figure 10(a) that is subject to the tip shear traction T, where the radius r and the thickness t are 9.0 and $1.0\mathrm{mm}$ and the Young’s modulus and Poisson’s ratio are set by $4.0\times \mathrm{GPa}$ and 0.25, respectively. The arch domain is uniformly discretized with 205 grid points, and the shear traction is increased up to $12\mathrm{MPa}$ in 10 uniform load steps. The problem was also solved by commercial nonlinear FEM code MSC/Marc using 1613 eight-node cubic elements. Figure 10(b) compares the load step-wise deformed shapes, where PG-NEM shows an excellent agreement with MSC/Marc up to the final load step. On the contrary, BG-NEM shows the gradual discrepancy in proportional to the load step. It is clearly observed from Figure 11 that comparatively represents the load-deflection curves in the horizontal and vertical directions, respectively.

Figure 10.

(a) An arch subject to the tip shear traction. (b) Comparison of deformed configurations (205 nodes).

Figure 11.

Comparison of load-deflection curves (205 nodes): (a) Horizontal. (b) Vertical.

Figure 12(a) depicts a hollow half sphere subject to a uniform line traction concentrated at the top center. The uniform traction is applied on a circle, which has the conical angle $\alpha$ of $4.5^o$ . The radius r and the thickness t of the half sphere are 9.0 and $1.0\mathrm{mm}$ , while the Young’s modulus E and Poisson’s ratio $\nu$ are the same with the previous arch problem shown in Figure 10(a). Two uniform NEM grids with different grid densities are taken to investigate the effect of grid density. As is the previous problems, the uniform traction is loaded up to $400\mathrm{MPa}$ in 10 uniform load steps. For comparison purpose, the problem was also solved by CS-FEM and MSC/Marc. For the sake of the reference solution, MSC/Marc was solved using a fine mesh that was generated with the total of 1613 eight-node quadratic finite elements.

Figure 12.

(a) A hollow half-sphere subject to circular concentrated load. (b) Natural element grids.

Figure 13(a) and (b) compares the pressure–deflection curves at the top that are obtained with 205 and 405 grid points. The comparison confirms that the PG-NEM provides the most accurate solution that is almost the same with MSC/Marc for both coarse and fine grids. However, similar to the previous cantilever problem, CS-FEM shows the remarkably poor accuracy compared with the other two methods. The deformed shapes obtained by PG-NEM at three different load steps are represented in Figure 14(a), and the final deformed shapes of all the methods are compared in Figure 14(b). One can clearly notice that PG-NEM shows an excellent agreement with MSC/Marc, but the other two methods, particularly CS-FEM, are fairly different from MSC/Marc. Thus, the accuracy of PG-NEM has been justified through three model problems, which is thanks to the fact that Laplace interpolation function provides high smoothness and the Delaunay triangle-based test function completely resolves the numerical integration error.

Figure 13.

Comparison of vertical deflections: (a) 205 nodes. (b) 405 nodes.

Figure 14.

Deformed configurations: (a) PG-NEM at different steps. (b) Comparison at the final step.

Figure 15(a) represents a two-dimensional rectangular rubber component in plane strain, which is bonded by two rigid plates on its top and bottom. This example was dealt by Calvo et al. ²³ to investigate the hyperelastic neo-Hookean behavior for the rubber material by the $\alpha$ -shape-based NEM. The bottom plate is fixed while the top plated is allowed to move vertically, and the neo-Hookean material constants $C_1$ and $1/D$ are set by $0.035472\mathrm{MPa}$ and $3.8095\times {10}^{-9}P{a}^{-1}$ , respectively. For the compression simulation, a uniform $41\times 41$ NEM grid shown in Figure 15(b) is used. For the numerical integration for both the stiffness matrices, seven Gauss points were used.

Figure 15.

Constrained compression: (a) Rectangular rubber (unit, cm). (b) Uniform NEM grid (1681 points).

Table 1 compares the displacements at two Points A and B when the top plate is compressed by $14\mathrm{cm}.$ The compression was terminated in 34 uniform load steps, while the total iteration numbers for each load step were not uniform. It is found from the table that the present method leads to the displacement values close to those of Calvo et al. such that the relative errors are 3.251 ∼ 8.550%. The maximum relative error occurs for $u_y$ at Point B.

Table 1.

Comparison of displacements at two points (34 load steps).

Methods	Point A		Point B
	$u_x(\mathrm{cm})$	$u_y(\mathrm{cm})$	$u_x(\mathrm{cm})$	$u_y(\mathrm{cm})$
Present	1.7702	−10.1564	8.1052	−7.0460
Calvo et al. 23	1.8297	−10.0970	7.6732	−6.4910

Figure 16(a) and (b) comparatively represents the deformed configuration of rubber component at the vertical compression of 10 and $14\mathrm{cm}$ , respectively. In the former case, the compression simulation was terminated in 10 load steps. From the figure, both methods commonly show the lateral expansion in proportional to the vertical compression. In the compression of $10\mathrm{cm}$ , both methods produce almost the similar deformed shapes to each other, showing the clear symmetry with respect to the center vertical line. But, a noticeable difference is observed from Figure 16(b) for very large strain deformations. Calvo et al. shows a nonsymmetric deformation caused by a buckling instability effect, which is attributed to the nonsymmetry in the integration cells. In Calvo et al., ²³ it is reported that this effect does also occur at the half model, but not at a quarter model. However, it is clearly observed that this phenomenon completely disappears in the present method. Thus, the present method does not suffer from the buckling instability effect at the very large strain compression. It is because uniform Delaunay triangularization in our method always preserves the symmetry in integration meshes.

Figure 16.

Comparison of deformations at selected time steps (1681 points; unit, cm): (a) At $10\mathrm{cm}$ compression. (b) At $14\mathrm{cm}$ compression.

Figure 17 represents the displacement-reaction force curves, where the reaction forces for both nonsymmetric mesh and half model are less rigid than those of the other two, owing to the buckling effect.

Figure 17.

Comparison of displacement-reaction force curves (1681 points).

The present method predicts less reaction forces at small compression, while it approaches the quarter model without buckling as the compression increases. The difference between the present method and the quarter model is not large with the maximum relative difference of 12%.

Conclusion

A 2-D nonlinear PG-NEM based on Voronoi polygons and Delaunay triangles has been introduced for the geometrically large deformation analysis of elastic structures. The displacement increments were formulated according to the linearized total Lagrangian method incorporated with Taylor series expansion. A combined use of Laplace interpolation functions and the Delaunay triangle–supported CS test functions provides the high smoothness and easily and accurately deals with the essential boundary condition and the numerical integration. Through the numerical experiments, the validity and accuracy of the proposed method have been verified. In particular, the proposed method provided the high interpolation accuracy even for coarse grids. On the contrary, CS-FEM and BG-NEM suffer from the numerical inaccuracy owing to the high stiffness matrix and the numerical integration error. Furthermore, the proposed method exactly preserves the symmetry in deformed configuration at very large strain. Consequently, the potential of PG-NEM for analyzing the large deformation problems has been sufficiently justified.