Large deformation analysis of elastic bodies by nonlinear Petrov

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.

Keywords

Large deformation analysis nonlinear PG-NEM Laplace interpolation function constant strain (CS) basis function linearized total Lagrangian formulation

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 $ℵ$ of N distinct points $x_{I}$ called nodes: $ℵ = {x_{1}, x_{2}, \dots, x_{N} | x_{I} \in ℜ^{2}}$ in $ℜ^{2}$ , a first-order Voronoi diagram $D_{Vo}$ is defined with N Voronoi polygons $ϖ_{I}$ such that $D_{Vo} = {ϖ_{1}, ϖ_{2}, \dots, ϖ_{N} | \bar{ϖ_{1} \cup \cdot \cdot \cdot \cup ϖ_{N}} = ℜ^{2}}$ . Here, the Voronoi polygon $ϖ_{I}$ corresponding to the Ith node is defined by $ϖ_{I} = {x \in ℜ^{2} : d (x, x_{I}) < d (x, x_{J}), \forall J \neq I}$ with $d (x, x_{I})$ being the Euclidean metric in $ℜ^{2}$ , which is a subdomain of the unbounded domain $ϖ \in ℜ^{2}$ with its sides, which perpendicularly bisect the lines connecting $x_{I}$ and the adjacent neighbor nodes in $ℵ$ . 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 $ϖ_{CH} (ℵ)$ of a set $ℵ$ . The convex hull becomes the problem domain: $ϖ_{CH} (ℵ) \equiv Ω \subset ϖ$ 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 $ℑ_{De}$ of Delaunay triangles, $Δ_{JKL} (x)$ , $ℑ_{De} = {Δ_{JKL} (x) : ⋃ \bar{Δ_{JKL}} = ϖ_{CH} (ℵ), J \neq K \neq L}$ , where J, K, and L refer to three distinct vertex nodes $x_{J}, x_{K}$ , and $x_{L}$ of $Δ_{JKL}$ . In other words, the problem domain $Ω$ 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 $cir (Δ_{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 $ℵ$ . It is not hard to realize that the centers of these circumcircles are identical to the Voronoi vertices.

To each node $x_{I} \in ℵ$ , exactly one Laplace interpolation function $φ_{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 $φ_{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 $ϖ_{I}$ , one can construct additional Voronoi polygon $ϖ_{P}$ composed of four subregions, which are divided by the previously defined four Voronoi polygons $ϖ_{I}, ϖ_{J}, ϖ_{K}$ , and $ϖ_{L}$ . Here, four subregions are called the second-order Voronoi cells $ϖ_{PI}$ , and those are mathematically defined by $ϖ_{PI} = {x \in ℜ^{2} : d (x, x_{P}) < d (x, x_{I}) < d (x, x_{J}), \forall J \neq P, I}$ . The number of $ϖ_{PI}$ is not identical with the number of nodes, furthermore the shapes of $ϖ_{P}$ and $ϖ_{PI}$ vary depending on the location of $x_{P}$ .

With the edge lengths $s_{ℓ} = | Γ_{P ℓ} |$ of $ϖ_{P}$ and the relative half distances $h_{ℓ} = d (x_{P}, x_{ℓ}) / 2$ between $x_{P}$ and the neighbor nodes, the weighting factors $α_{ℓ}$ are defined by

α_{ℓ} (x_{P}) = \frac{s_{ℓ} (x_{P})}{h_{ℓ} (x_{P})}, ℓ = I, J, K, L

(1)

Then, the value of Laplace interpolation function $φ_{I}$ and its derivatives $φ_{I, β} (β = x, y)$ at a point $x_{P}$ are calculated by Chinesta et al⁸

\begin{matrix} φ_{I} (x_{P}) = α_{I} (x_{P}) / \sum_{ℓ = 1}^{M} α_{ℓ} (x_{P}), φ_{I, β} (x_{P}) \\ = [α_{I, β} (x_{P}) - φ_{I} (x_{P}) \sum_{ℓ = 1}^{M} α_{ℓ, β} (x_{P})] / \sum_{ℓ = 1}^{M} α_{ℓ} (x_{P}) \end{matrix}

(2)

with $α_{I, β} (x_{P}) = [s_{I, β} (x_{P}) - α_{I} h_{I, β} (x_{P})] / h_{I} (x_{P})$ . And, the Laplace interpolation function $φ_{I}$ corresponding to the I th node $x_{I}$ is defined by replacing $x_{P}$ with an arbitrary point x within the support $supp (φ_{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 $supp (φ_{I})$ of $φ_{I}$ is defined by the intersection of the convex hull $ϖ_{CH} (ℵ)$ and the union $\cup \bar{cir (Δ_{IJK})}$ of Delaunay circumcircles defined by the node $x_{I}$ and its neighbor nodes, $supp (φ_{I} (x)) = ⋃ \bar{cir (Δ_{IJK} (x))} \cap ω_{CH} (ℵ)$ . 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, $Ω^{0} \in ℜ^{2}$ . And, let $Ω^{t + Δ t}$ be the deformed equilibrium with its boundary $Γ^{t + Δ t}$ at the current time, $t + Δ 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 $ξ_{t + Δ t} (x^{0})$ from $Ω^{0}$ to $Ω^{t + Δ t}$ such that

x^{t + Δ t} = ξ_{t + Δ t} (x^{0}) = x^{0} + u^{t + Δ t}

(3)

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

F_{ij}^{t + Δ t} = \frac{\partial ξ_{t + Δ t} (x_{i}^{0})}{\partial x_{j}^{0}} = δ_{ij} + \frac{\partial u_{i}^{t + Δ 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 + Δ t$ is governed by the static equilibrium given by

\frac{\partial σ_{ij}^{t + Δ t}}{\partial x_{j}^{t + Δ t}} + b_{i}^{t + Δ t} = 0 in Ω^{t + Δ t}

(5)

with essential and traction boundary conditions

u_{i}^{t + Δ t} = {\hat{u}}_{ij}^{t + Δ t} on Γ_{D}^{t + Δ t}

(6)

σ_{ij}^{t + Δ t} n_{j}^{t + Δ t} = {\hat{t}}_{i}^{t + Δ t} on Γ_{N}^{t + Δ t}

(7)

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

\begin{matrix} \int_{Ω^{0}} S_{ij}^{t + Δ t} (u^{t + Δ t}) {\tilde{E}}_{ij}^{t + Δ t} (u^{t + Δ t}, v) d Ω^{0} \\ = \int_{Γ_{N}^{0}} {\hat{t}}_{i}^{0} v_{i} d Γ^{0} + \int_{Ω^{0}} b_{i}^{0} v_{i} d Ω^{0} \end{matrix}

(8)

for every virtual displacement v.

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

\begin{matrix} {\tilde{E}}_{ij}^{t + Δ t} (u^{t + Δ t}, v) \\ = \frac{1}{2} [\frac{\partial v_{i}}{\partial x_{j}^{0}} + \frac{\partial v_{j}}{\partial x_{i}^{0}} + \frac{\partial u_{m}^{t + Δ 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 + Δ t}}{\partial x_{j}^{0}}] \end{matrix}

(9)

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

\begin{matrix} E_{ij}^{t + Δ t} & = \frac{1}{2} [F_{mi}^{t + Δ t} F_{mj}^{t + Δ t} - δ_{ij}] \\ = \frac{1}{2} [\frac{\partial u_{i}^{t + Δ t}}{\partial x_{j}^{0}} + \frac{\partial u_{j}^{t + Δ t}}{\partial x_{i}^{0}} + \frac{\partial u_{m}^{t + Δ t}}{\partial x_{i}^{0}} \frac{\partial u_{m}^{t + Δ t}}{\partial x_{j}^{0}}] \end{matrix}

(10)

The weak form equation (8) is nonlinear because $E_{ij}^{t + Δ t}$ includes the displacement $u^{t + Δ 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_{ij}^{t + Δ t} {\tilde{E}}_{ij}^{t + Δ t} \approx S_{ij}^{(k)} {\tilde{E}}_{ij}^{(k)} + \frac{\partial}{\partial u_{m}^{(k)}} (S_{ij}^{(k)} {\tilde{E}}_{ij}^{(k)}) u_{m}^{Δ^{(k + 1)}}

(11)

Referring to Figure 4, two superscripts $(k)$ and $Δ^{(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}^{(k)}} (S_{ij}^{(k)} {\tilde{E}}_{ij}^{(k)}) u_{m}^{Δ^{(k + 1)}} = C_{ijpq}^{(k)} E_{pq}^{Δ^{(k + 1)}} {\tilde{E}}_{ij}^{(k)} + S_{ij}^{(k)} {\tilde{E}}_{ij}^{Δ^{(k + 1)}}

(12)

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

\begin{matrix} E_{ij}^{Δ^{(k + 1)}} (u^{Δ^{(k + 1)}}) \\ = \frac{1}{2} [\frac{\partial u_{i}^{Δ^{(k + 1)}}}{\partial x_{j}^{0}} + \frac{\partial u_{j}^{Δ^{(k + 1)}}}{\partial x_{i}^{0}} + \frac{\partial u_{m}^{Δ^{(k + 1)}}}{\partial x_{i}^{0}} \frac{\partial u_{m}^{(k)}}{\partial x_{j}^{0}} + \frac{\partial u_{m}^{(k)}}{\partial x_{i}^{0}} \frac{\partial u_{m}^{Δ^{(k + 1)}}}{\partial x_{j}^{0}}] \end{matrix}

(13)

{\tilde{E}}_{ij}^{Δ^{(k + 1)}} (u^{Δ^{(k + 1)}}, v) = \frac{1}{2} [\frac{\partial v_{m}}{\partial x_{i}^{0}} \frac{\partial u_{m}^{Δ^{(k + 1)}}}{\partial x_{j}^{0}} + \frac{\partial u_{m}^{Δ^{(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}^{(k)} = A^{(k)} t_{i}^{t + Δ t} (u^{(k)}), b_{i}^{(k)} = J^{(k)} b_{i}^{t + Δ t} (u^{(k)})

(15)

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

\begin{matrix} \int_{Ω^{0}} C_{ijpq}^{(k)} E_{pq}^{Δ^{(k + 1)}} {\tilde{E}}_{ij}^{(k)} d Ω^{0} + \int_{Ω^{0}} S_{ij}^{(k)} {\tilde{E}}_{ij}^{Δ^{(k + 1)}} d Ω^{0} = \int_{Γ_{N}^{0}} t_{i}^{(k)} v_{i} d Γ^{0} \\ + \int_{Ω^{0}} b_{i}^{(k)} v_{i} d Ω^{0} - \int_{Ω^{0}} S_{ij}^{(k)} {\tilde{E}}_{ij}^{(k)} d Ω^{0} \end{matrix}

(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{matrix} W^{(k)} = \frac{1}{2} C_{ijpq}^{(k)} E_{pq}^{(k)} E_{ij}^{(k)}, \\ C_{ijpq}^{(k)} = C_{ijpq} = λ δ_{ij} δ_{pq} + μ (δ_{ip} δ_{jq} + δ_{iq} δ_{jp}) \end{matrix}

(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^{Δ^{(k + 1)}}$ at the (k + 1)th stage and the test displacement v are expanded as

u_{i}^{Δ^{(k + 1)}} (x^{0}) = \sum_{I = 1}^{N} φ_{I} (x^{0}) {\bar{u}}_{Ii}^{Δ^{(k + 1)}}

(18)

v_{i} (x^{0}) = \sum_{J = 1}^{N} ψ_{J} (x^{0}) {\bar{v}}_{Ji}

(19)

Similarly, the displacements at the kth stage (i.e. at time $t + Δ t$ ) and at time stage t are expressed by

u_{i}^{(k)} (x^{0}) = \sum_{I = 1}^{N} φ_{I} (x^{0}) {\bar{u}}_{Ii}^{(k)}

(20)

u_{i}^{t} (x^{0}) = \sum_{I = 1}^{N} φ_{I} (x^{0}) {\bar{u}}_{Ii}^{t}

(21)

for a given natural element grid composed of N nodes. Here, $ϕ_{I}$ are Laplace interpolation functions, while $ψ_{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 $Ω_{IJ}^{int}$ between trial function and test function: (a) $φ_{I}$ and $φ_{J}$ in PG-NEM. (b) $ψ_{I}$ and $φ_{J}$ in BG-NEM.

Three shape functions ${\hat{ψ}}_{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{ψ}}_{i} (x^{0}) = \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}}}{2 A}

(22)

The approximation becomes the BG-NEM when both the trial and test displacements are expanded by Laplace interpolation function, $φ_{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}^{(k)} + {{}^{g}K}_{I}^{(k)}] {\bar{u}}^{Δ^{(k + 1)}} = \sum_{I}^{N} {{{}^{e}F}_{I}^{(k)} - {{}^{i}F}_{I}^{(k)}}

(23)

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

{{}^{m}K}_{I}^{(k)} = \int_{Ω_{v}^{I}} {(B_{v}^{(k)})}^{T} \bar{C} B_{u}^{(k)} d Ω^{0}

(24)

{{}^{g}K}_{I}^{(k)} = \int_{Ω_{v}^{I}} {(H_{v}^{(k)})}^{T} {\hat{S}}^{(k)} H_{u}^{(k)} d Ω^{0}

(25)

{{}^{e}F}_{I}^{(k)} = \int_{Ω_{v}^{I}} Ψ^{T} b^{(k)} d Ω^{0} + \int_{Γ_{N}^{0} \cap Ω_{v}^{I}} Ψ^{T} t^{(k)} d Γ^{0}

(26)

{{}^{i}F}_{I}^{(k)} = \int_{Ω_{v}^{I}} {({\tilde{B}}_{v}^{(k)})}^{T} {\hat{S}}^{(k)} d Ω^{0}

(27)

with $Ω_{v}^{I} = Ω^{0} \cap supp (ψ_{I})$ , ${\overset{\cdot}{S}}^{(k)} = {S_{11}^{(k)}, S_{22}^{(k)}, S_{12}^{(k)}, S_{33}^{(k)}}^{T}$ , the $(3 \times 3)$ material modulus tensor ${\bar{C}}^{(k)}$ , the $(2 \times 2 N)$ basis function matrices $Φ$ and $Ψ$ , $B_{v}^{(k)} = [B_{v}^{1} \dots B_{v}^{I} \dots B_{v}^{N}]$ , $B_{u}^{(k)} = [B_{u}^{1} \dots B_{u}^{J} \dots B_{u}^{N}]$ , $H_{v}^{(k)} = [H_{v}^{1} \dots H_{v}^{I} \dots H_{v}^{N}]$ , and $H_{u}^{(k)} = [H_{u}^{1} \dots H_{u}^{J} \dots 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 $supp (ψ_{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 $ψ_{I}$ is identical with a union of Delaunay triangles, but also the intersection $Ω_{IJ}^{int}$ between the supports of $φ_{J}$ and $ψ_{I}$ is always contained within $supp (ψ_{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

Δ {\bar{u}}_{rel} = \frac{| {\bar{u}}^{Δ^{(k + 1)}} |}{| {\bar{u}}^{(k + 1)} |} \leq α_{u} | {\bar{u}}^{(k + 1)} |

(28)

Δ {\bar{f}}_{rel} = \frac{| {{}^{e}F}^{(k + 1)} - {{}^{i}F}^{(k + 1)} |}{| {{}^{e}F}^{(k + 1)} |} \leq α_{f} | {{}^{e}F}^{(k + 1)} |

(29)

where, $Δ {\bar{u}}_{rel}$ and $Δ {\bar{f}}_{rel}$ mean the magnitudes of the relative displacement increment and the relative force residual, and $α_{u}$ and $α_{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 $θ$ 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 $θ$ becomes smaller than the preset minimum value $θ_{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 mm$ and the section dimensions are set by $t = h = 1.0 mm$ , while the Young’s modulus E and the Poisson’s ratio are set by $400 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 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 $(δ_{h} / L)$ . (b) Vertical $(δ_{v} / L)$ .

Figure 9.

Comparison of the normalized deflections (105 nodes): (a) Horizontal $(δ_{h} / L)$ . (b) Vertical $(δ_{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 mm$ and the Young’s modulus and Poisson’s ratio are set by $4.0 \times GPa$ and 0.25, respectively. The arch domain is uniformly discretized with 205 grid points, and the shear traction is increased up to $12 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 $α$ of $4 . 5^{o}$ . The radius r and the thickness t of the half sphere are 9.0 and $1.0 mm$ , while the Young’s modulus E and Poisson’s ratio $ν$ 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 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 $α$ -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 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 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
Methods	$u_{x} (cm)$	$u_{y} (cm)$	$u_{x} (cm)$	$u_{y} (cm)$
Present	1.7702	−10.1564	8.1052	−7.0460
Calvo et al.²³	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 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 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 cm$ compression. (b) At $14 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.

Footnotes

Appendix 1

Handling Editor: Daxu Zhang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant NRF-2017R1D1A1B03028879).

ORCID iD

JR Cho

References

Belytschko

Element-free Galerkin methods. Int J Numer Methods Eng 1994; 37: 229–256.

Liu

Jun

Zhang

YF.

Reproducing kernel particle methods. Int J Numer Methods Fluid 1995; 20: 1081–1106.

Duarte

Oden

JT.

An H-P adaptive method using clouds. Comput Method Appl M 1996; 139: 237–262.

Melenk

Babuska

The partition of unity finite element method: basic theory and applications. Comput Method Appl M 1996; 139: 289–314.

Krongauz

Belytschko

Enforcement of essential boundary conditions in meshless approximations using finite elements. Comput Method Appl M 1996; 131: 133–145.

Atluri

SN.

A modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin method. Comput Mech 1998; 21: 211–222.

Braun

Sambridge

A numerical method for solving differential equations on highly irregular evolving grids. Nature 1995; 376: 655–660.

Chinesta

Cescotto

Cueto

et al . Natural element method for the simulation of structures and processes. Hoboken, NJ: John Wiley & Sons, 2013.

Cho

Lee

Yoo

WS.

Natural element approximation of Reissner-Mindlin plate for locking-free numerical analysis of plate-like thin elastic structures. Comput Method Appl M 2014; 256: 17–28.

10.

Sukumar

Moran

C1 natural neighbor interpolant for partial differential equations. Numer Meth Partial Diff Eq 1999; 15: 417–447.

11.

Alfaro

Bel

Cueto

et al . Three-dimensional simulation of aluminum extrusion by the shape based natural element method. Comput Method Appl M 2006; 195: 4269–4286.

12.

Nie

Zhou

Han

et al . C1 natural element method for strain gradient linear elasticity and its application to miscrostructures. Acta Mech Sinica 2012; 28: 91–103.

13.

Zhou

Liu

Chen

Upper bound limit analysis of plates utilizing the C1 natural element method. Comput Mech 2012; 50: 543–561.

14.

Zhang

Tan

HP.

Least-square natural element method for radiative heat transfer in graded index medium with semitransparent surfaces. Int J Heat Mass Transfer 2013; 66: 349–354.

15.

Teyssedre

Gilormini

Extension of natural element method to surface tension and wettability for the simulation of polymer flows at the micro and nano scales. J Non-Newtonian Fluid Mech 2013; 200: 9–16.

16.

Marechal

Ramdane

Natural element method applied to electromagnetic problems. IEEE T Magn 2013; 49: 1714–1716.

17.

Strang

Fix

GJ.

An analysis of the finite element method. Upper Saddle River, NJ: Prentice-Hall, 1973.

18.

Cho

Lee

HW.

A Petrov-Galerkin natural element method securing the numerical integration accuracy. J Mech Sci Technol 2006; 20: 94–109.

19.

Wang

Zhou

Shan

The natural neighbour Petrov–Galerkin method for elasto-statics. Int J Numer Method Eng 2005; 63: 1126–1145.

20.

Logan

DL.

A first course in the finite element method. Stanford, CA: Cengage, 2001.

21.

Holzapfel

GA.

Nonlinear solid mechanics: a continuum approach for engineering. Chichester: John Wiley & Sons, 2000.

22.

Doyle

JF.

Nonlinear analysis of thin-walled structures: statics, dynamics and stability. New York: Springer, 2001.

23.

Calvo

Martinez

Doblaré . On solving large strain hyperelastic problems with the natural element method. Int J Numer Methods Eng 2005; 62: 159–185.

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

Abstract

Keywords

Introduction

Laplace interpolation functions

Linearized total Lagrangian formulation

Two-dimensional nonlinear Petrov–Galerkin natural element approximation

Numerical experiments

Conclusion

Footnotes

Appendix 1

Declaration of conflicting interests

Funding

ORCID iD

References