Large deflection and initial instability analysis of anisotropic plates by the generalized finite element method

Abstract

This paper presents investigations laminated plates under moderately large transverse displacements and initial instability, through the Generalized Finite Element Methods - GFEM. The von Kármán plate hypothesis are used along with Kirchhoff and Reissner-Mindlin kinematic plate bending models to approximate transverse displacements and critical buckling loads. The generalized approximation functions are either C ⁰or C ^k continuous functions, with k being arbitrarily large. It is well known that in GFEM, when both the partition of unity (PoU) and the enrichments functions are polynomials, the stiffness matrices are singular or ill conditioned, which causes additional difficulties in applications that requires the solution of algebraic eigenvalues problems, like in the determination of natural frequencies of vibration or the initial buckling loads. Some investigations regarding this problem are presently addressed and some aspects and advantages of using C ^k -GFEM are observed. In addition, comparisons are presented between the classical GFEM and the Stable-GFEM (SGFEM) with regard to the evaluation of the initial critical buckling loads. The numerical experiments use reference values from analytical and numerical results obtained in the open literature.

Keywords:
Large displacements in plates; laminated plate bending; GFEM; SGFEM; arbitrarily continuous approximation functions.

1 INTRODUCTION

Anisotropic laminated plates are amongst the most commonly used types of structural elements in aeronautical and naval industry. They are often subjected to in-plane compressive and shear loads which, when combined with their slenderness, make then prone to geometric instability by buckling, where various parameters such as thickness, layer stack, loading type, boundary conditions and geometry affect its ability to maintain structural integrity. Plate theories, through some kinematic hypotheses, allow the three-dimensional fields to be approximated by two-dimensional ones, thus reducing the complexity and the number of variables required to estimate its mechanical behavior. However, these simplifications lead to important mathematical consequences, such as:

• (a) First order kinematic models such as Kirchhoff and Mindlin, within the linear conditions, allow solutions to be exact with regard to the three-dimensional formulation, in the limit as the transverse displacement becomes small compared to the thickness. A traditional rule of thumb of application of these theories is that the errors involved are admissible in engineering applications if the transverse displacements are limited to some fraction of the thickness, e.g., a half.

• (b) The simplest inclusion of nonlinear kinematic effects is provided by the von Kármán theory. The resultant strain-deformation relation contain a sub-set of the nonlinear terms present in the full relations. As a consequence, the range of applicability of the theory is usually referred to as displacements moderately large. This means it provides better accuracy in the same displacement range of validity of the linear models, but does not allow arbitrarily large displacements. A second aspect of the von Kármán theory is that it allows a first approximation of initial instability of the panel. This is obtained numerically by a standard linear eigenvalue problem which gives, in a relatively inexpensive way, an estimate for the mode shape and critical load of the structure. Even though the load estimate is only an upper bound (Brush and Almroth, 1975; Ventsel and Krauthammer, 2001), it is often very useful in the early stages of the design process.

With regard to finite element (FE) solutions, the plate bending theories available contain different continuity requirements on the function basis used. For instance, the Mindlin kinematic model requires $C^0$ continuous approximation functions and the Kirchhoff model requires $C^1$ functions for the transverse displacements. The Kirchhoff model is interesting by its reduced cost in computational analysis and absence of some of the pathologies which can suffer some of the FE formulations based on the Reissner-Mindlin model, like the shear locking. Frequently, the technologies developed to FE formulations on higher order kinematic models, like the third order shear deformation of Levinson or Reddy (Reddy, 1989, 1984), are inspired in developments made for the Kirchhoff model by using, for example, non-conforming elements to alleviate the $C^1$ requirement on the function basis.

As for the $C^1$ requirement for the Kirchhoff finite element basis, it is virtually unknown in the literature a formulation of arbitrarily shaped plate elements with $C^1$ continuity everywhere in the mesh. Most of them have to be restricted to particular element shapes, or the continuity is restricted to discrete points, only enough to preclude numerical instability. Only recently a few approaches have been developed to fix this limitation. One of them is a variation of the GFEM (Generalized Finite Element Method), based on smooth $C^k$ continuous approximation functions, for arbitrary $k$ (de Barcellos et al., 2009). This is one of the methods which will be tested in this paper, on modeling the nonlinear effects of anisotropic laminated plates by the von Kármán theory with the Kirchhoff and Mindlin kinematic models.

The theoretical and historical basis of the GFEM lies on early developments of meshless methods, particularly on the hp-Clouds method (Duarte and Oden, 1996a, 1996b). These methods share some difficulties, for example on integration and on imposition of boundary conditions. Oden et al. (1998) and, in parallel, Babuška (Babuška and Melenk, 1997; Melenk and Babuška, 1996) proposed the use of clouds (patches of elements around each node) as support of approximation basis functions, with the contours conveniently defined by linear finite element meshes, but keeping from the hp-Cloud method its flexible way of enriching the basis with arbitrary functions defined in global coordinates. After this point the GFEM/XFEM developed exponentially, for example, (Belytschko et al., 2009; Belytschko and Black, 1999; Dolbow et al., 2000; Moës et al., 1999; Sukumar et al., 2000).

Most of GFEM, as well as FEM developments, are formulated with $C^0$ approximation function basis. However, in applications such as the Kirchhoff plate model the required approximation functions must be at least $C^1$ continuous. The $C^k$ version of GFEM has been developed by the use the Edwards proposal (Edwards, 1996), which allows the construction of arbitrarily smooth weighting functions with support in convex clouds. These weighting functions are used to form a Shepard Partition of Unit (Duarte and Oden, 1995; Shepard, 1968). Duarte et al. (Duarte et al., 2006) extends Edwards proposal to non-convex clouds, by the use of a Boolean R-function.

A known characteristic of the GFEM formulation is that its set of approximating functions is linearly dependent, when it is formed by a polynomial partition of unity uniformly enriched with polynomial functions. As a result, the coefficient matrix (stiffness) contains more null eigenvalues than the number of rigid body motions in the model. When different combinations of types of partition of unity and enrichment functions are used, the rank deficiency can be eliminated but the stiffness matrix remains ill conditioned. Babuška (Babuška and Banerjee, 2012, 2011) proposed a modification in the enrichment functions, in which they are substituted by the difference between the original one and a linear interpolating function, defined in each element. This method, referred to as Stable-GFEM (SGFEM), was shown to improve the approximation space and also to alleviate problems due to non-homogeneous enrichment over the mesh. Gupta et al (2013) compared accuracy and conditioning between GFEM and SGFEM in a two-dimensional fracture mechanics problem, and showed that SGFEM improves the accuracy while keeping the conditioning similar to the one of FEM. Later, Gupta et al (2015) extended and improved the SGFEM formulation to the three-dimensional fracture mechanics.

A practical consequence of the rank deficiency in the stiffness matrix of $C^0$ -GFEM is that the number n of null eigenvalues is, a-priori, unknown, and can be relatively large when compared with the order $N$ of the matrix. The standard procedure to evaluate the critical load consists in solving an eigenvalue matrix problem and compute the first, smallest, positive eigenvalue. Since $N$ is large, an adequate method is used to compute not all $N$eigenvalues but only the first few of them. However, the difficulty is posed when one does not know, a-priori, the dimension of the kernel space. This forces the eigenvalue algorithm to work in a vector space of size, at least, $n+1$. On the other hand, the function basis formed by the smooth $C^k$ partition of unity uniformly enriched by polynomials is linearly independent (Mendonça et al., 2011). As a result, after imposition of sufficient Dirichlet boundary conditions, the stiffness matrix has a large condition number, but it is not singular.

The objectives of the present paper are the evaluation of the behavior of both types of GFEM, $C^0$ and $C^k$ in deformation response of Mindlin and Kirchhoff kinematic models under moderately large displacements using the von Kármán theory. Also, the response of both stable and non-stable $C^0$-GFEM in the evaluation of buckling load is analyzed, with regard to the effects of the stiffness matrix rank on the computational effort to obtain the first useful eigenvalue.

The present paper is organized as follows. Section 2 summarizes a few key aspects of the GFEM, in both forms $C^0$and $C^k$, and its stable version; Section 3 shows the basic information about Kirchhoff and Mindlin plate theory with von Kármán hypothesis and its discretization for a generalized finite element procedure; Section 4 details aspects of the initial instability analysis by GFEM and Section 5 presents comparative results for deformation analysis and critical load evaluation in standard problems.

PARTITION OF UNITY AND APPROXIMATION FUNCTIONS

Let us consider a region V, defined by a thickness$H$and a planar middle surface $\Omega$ with boundary region $\Gamma$. The $V$ region may be described, in Cartesian coordinates, as

$V=\left\{(x,y,z)\in R^3|withz\in \left[-\frac{H}{2},\frac{H}{2}\right]andx:=(x,y)\in \Omega ,\Omega \in R^2\right\}$ (2.1)

Let us consider a mesh of elements of domain ${\Omega }_e$, such that $\cup {\Omega }_e=\Omega$ and ${\Omega }_i\cap {\Omega }_j=\varphi$ for $i\ne j$. Define a cloud (a patch of elements) ${\omega }_{\alpha }$ associated to a node $\alpha$ as the union of all element domains that share the node $\alpha$. In the present tests, the mid-surface of a laminated plate is modeled as an open bounded domain, $\Omega \subset$ $R^2$, with a linear triangulation made of $N_{el}$ elements $K_e$, ${\left\{K_e\right\}}_{j=1}^{N_{el}}$, and $N_{node}$ nodes with coordinates ${\left\{x_{\alpha }\right\}}_{\alpha =1}^{N_{node}}$.

A Partition of Unit (PoU) is defined (Oden and Reddy, 1976) as a set of functions ${\varphi }_{\alpha }(x)$, associated with nodes $\alpha =\mathrm{1,}\mathrm{2,...,}{N}_{nos}$, such that: 1) ${{\displaystyle }}_{\alpha =1}^{N_{nos}}{\varphi }_{\alpha }(x)=\mathrm{1,}\forall x\in \Omega$; 2) ${\varphi }_{\alpha }$ is a function at least continuous over the domain, has having its cloud, ${\omega }_{\alpha }$, as a compact support and has the required smoothness. Here, the partition of unity functions have at least $k$ continuous derivative inside the cloud, that is, ${\varphi }_{\alpha }(x)\in C_0^k({\omega }_{\alpha }),k\ge 0$.

Once a PoU is available, the GFEM approximation functions ${\phi }_{\alpha i}(x)$ are obtained by enriching the PoU functions, with enrichment functions $L_{\alpha i},i=\mathrm{1,...,}{n}_e$. That is, for each cloud ${\varphi }_{\alpha }$,

${\phi }_{\alpha i}(x)={\varphi }_{\alpha }(x)L_{\alpha i}(x)$ (2.2)

The set of enrichment functions can be chosen in different ways. The simplest type is the set of monomials of a given degree. Also, it is very common the use a special function extracted from an asymptotic solution around a critical region of the model. In general, the enrichment functions are $C^{\infty }$ continuous over an infinite support. The product with ${\varphi }_{\alpha }$ makes the approximation function ${\phi }_{\alpha i}$ have compact support and inherit the same level of smoothness as the PoU.

The usual GFEM is based on PoU functions which are the traditional FE shape functions, usually the linear tent functions, which naturally constitute a PoU. These functions have low computational cost and are easily integrated by numerical quadrature, although it is limited to $C^0$ continuity.

1.1 Shepard Partition of Unity

For problems requiring higher smoothness of the primal variables, like the transverse displacements in the Kirchhoff and higher order Reddy plate models, the linear tent functions used as PoU in the $C^0$ -GFEM are not adequate and a different approach is needed to generate approximation basis functions with the required level of continuity. Shepard (1968) proposed a procedure to generate smooth PoU by means of smooth weight functions. Let a function $W_{\alpha }:R^2\to R$ with a compact support, ${\omega }_{\alpha }$, belonging to the space $C_0^k({\omega }_{\alpha })$, be the weight function in each cloud ${\omega }_{\alpha }$, of the open coverage $F_{N_{nos}}$ of the domain $\Omega$. Then, the Shepard PoU subordinated to the cover $F_N$ is defined as:

${\varphi }_{\alpha }(x)=\frac{W_{\alpha }}{{\displaystyle \sum \beta \left(x\right)}{W}_{\beta }\left(x\right)}for\beta \left(x\right)\in \left\{\gamma |W_{\gamma }\left(x\right)\ne 0\right\}.$ (2.3)

The regularity of ${\varphi }_{\alpha }(x)$ depends only on the regularity of the weighting functions. Therefore, it remains to define the weight functions with the necessary regularity. In order to do that, a coordinate ${\xi }_j(x)$ normal to the straight edge $j$ of the cloud is defined. It is computed as ${\xi }_j(x)=n_{\alpha ,j}\cdot \left(x-b_{\alpha ,j}\right),$ for $x\in {\omega }_{\alpha }$, where $b_{\alpha ,j}$ is the coordinate of the mid-point of the edge $j$, and $n_{\alpha ,j}$ is the unity vector normal to the edge pointing to the inside of the cloud. Therefore ${\xi }_j(x)$ is the smallest distance between x and the edge $j$. Cloud edge functions with continuity $C^{\infty }$ can be obtained using the following edge function

${\epsilon }_{\alpha ,j}\left[{\xi }_j(x)\right]={\widehat{\epsilon }}_{\alpha ,j}:=\{\begin{array}{ll}AExp\left[-{\frac{{\xi }_j}{\beta }}^{-\gamma }\right]\hfill & if0<{\xi }_j,\hfill \\ 0\hfill & othercase.\hfill \end{array}$ (2.4)

The parameter values are chosen as $\gamma =0.6$ and $\beta =0.3$ according to Duarte et al. (2006). The weighting functions determine the influence of each node within a single element and are defined by

$W_{\alpha }(x)=Exp\left[c_{\alpha }\right]\underset{j=1}{{\displaystyle \prod }_{}^{}^{M_{\alpha }}}{\epsilon }_{\alpha ,j}\left({\xi }_j\right),$ (2.5)

where $c_{\alpha }$ is a scale parameter that guarantees $W_{\alpha }=1$ in the node $x_{\alpha }$ and $M_{\alpha }$ is the number of edges on the boundary of the cloud.

The Shepard PoU functions constructed according to (2.3)-(2.5), as observed by Duarte et al. (2006), have $C^{\infty }$ differentiability for convex clouds. The regularity of the approximation functions is determined by the regularity of the PoU and the enrichment functions used. Since in the present paper only the polynomial enrichment is tested, the resultant approximation functions basis is $C^{\infty }$. If non convex clouds were used, the result would be $C^{\infty }$ everywhere except at the concave nodes, where it would be $C^k$, after using an appropriate Boolean product. In addition to being able to choose a basis that satisfy the differentiability required by any formulation, it is possible to use different regularities for each unknown field through the combined use of $C^0$ and $C^k$ functions (Mendonça et al., 2013).

1.2 Stable-GFEM

The first procedure of stabilization in GFEM was proposed by Babuška (Babuška and Banerjee, 2012, 2011) as a simple local modification in the enrichment functions employed in $C^0$-GFEM. They consist in the difference between the original enrichment function and a linear interpolant function defined in each finite element:

$L_{\alpha i}^{\mathrm{mod}}=L_{\alpha i}-I_{w_{\alpha }}(L_{\alpha i})$ (2.6)

where $L_{\alpha i}^{\mathrm{mod}}$ is the modified SGFEM enrichment function, $L_{\alpha i}$ is the original enrichment function on the cloud ${\omega }_{\alpha }$, $I_{w_{\alpha }}(L_{\alpha i})$ is the interpolant function of the enrichment function defined at the element nodes, obtained with piecewise linear FE tent functions.

The approximation functions are obtained by the product of PoU functions with modified enrichment functions in the same way as in (2.2)

${\phi }_{\alpha i}(x)={\varphi }_{\alpha }(x)L_{\alpha i}^{\mathrm{mod}}(x).$ (2.7)

PLATE MODELS AND GFEM DESCRETIZATION

The kinematic assumptions of Mindlin's model for plates can be summarized as the following, in Cartesian displacement components:

$\begin{array}{c}{u}_x\left(x,z\right)=u_0\left(x\right)+z{\theta }_x\left(x\right)\\ u_y\left(x,z\right)=v_0\left(x\right)+z{\theta }_y\left(x\right)\\ u_z\left(x,z\right)=w\left(x\right).\end{array}$ (3.1)

and the Kirchhoff hypothesis result in ${\theta }_x=-\partial w/\partial x$ and ${\theta }_y=-\partial w/\partial y$. Therefore, Mindlin model involve five generalized displacements and Kirchhoff model three.

For this type of equivalent layer kinematic model, the independent strains are the in-plane strains $\epsilon \left(x,z\right):={\left\{{\epsilon }_x;{\epsilon }_y;{\gamma }_{xy}\right\}}^T$ and the transverse shear strains ${\gamma }_s\left(x,z\right):={\left\{{\gamma }_{yz};{\gamma }_{xz}\right\}}^T$ at an arbitrary point of the plate. The hypothesis (3.1) allows the in-plane deformations to be represented in the separated in the form $\epsilon (x,z)=\epsilon \left(x\right)+z\kappa (x)$, where $\epsilon$ is the membrane strain, and $\kappa$ is the curvature change due to bending. From the von Kármán hypothesis, the membrane deformations $\epsilon$ is divided into linear and non-linear parts ${\epsilon }^0$, ${\epsilon }^{NL}$, as $\epsilon ={\epsilon }^0+{\epsilon }^{NL}$. The complete set of generalized deformations is the following:

${\epsilon }^0=\left\{\begin{array}{c}\frac{\partial u_0}{\partial x}\\ \frac{\partial v_0}{\partial y}\\ \frac{\partial v_0}{\partial x}+\frac{\partial u_0}{\partial y}\end{array}\right\},{\epsilon }^{NL}=\left\{\begin{array}{c}\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\\ \frac{1}{2}{\left(\frac{\partial w}{\partial y}\right)}^2\\ \frac{\partial w}{\partial x}\frac{\partial w}{\partial y}\end{array}\right\},\kappa =\left\{\begin{array}{c}\frac{\partial {\theta }_x}{\partial x}\\ \frac{\partial {\theta }_x}{\partial y}\\ \frac{\partial {\theta }_x}{\partial y}+\frac{\partial {\theta }_y}{\partial x}\end{array}\right\}and{\gamma }^s=\left\{\begin{array}{c}\frac{\partial w}{\partial y}+{\theta }_x\\ \frac{\partial w}{\partial x}+{\theta }_y\end{array}\right\}.$ (3.2)

The transverse shear deformations ${\gamma }^s$ naturally results null in the Kirchhoff model.

The plate is considered an anisotropic laminated body, composed by $N$ layers characterized by the following anisotropic linear stress-strain relations (generalized Hooke's law), such that, for each layer $k_l$

$\sigma \equiv \left\{\begin{array}{l}{\sigma }_x\hfill \\ {\sigma }_x\hfill \\ {\tau }_{xy}\hfill \end{array}\right\}={\overline{Q}}^{k_l}\left(\epsilon +z\kappa \right);{\tau }_s\equiv \left\{\begin{array}{l}{\tau }_{yz}\hfill \\ {\tau }_{xz}\hfill \end{array}\right\}=E^{k_l}{\gamma }^s,$ (3.3)

where ${\overline{Q}}^{k_l}$ and $E^{k_l}$ are the elastic anisotropic material matrices of the layer, $\sigma$ and ${\tau }_s$ are in-plane and transverse shear stresses respectively. For the laminate, the constitutive relations are

$\overline{N}=C\overline{\epsilon }andQ=E^s{\gamma }^s$ (3.4)

where $\overline{N}:={\left\{N;M\right\}}^T$, with $N:={\{N_{xx},N_{yy},N_{xy}\}}^T$ and $M:={\{M_{xx},M_{yy},M_{xy}\}}^T$ being the resultant forces and moments, respectively, defined as $N={{\displaystyle }}_z\sigma dz$ and $M={{\displaystyle }}_z\sigma zdz$. $Q:={\{Q_x,Q_y\}}^T$ are resultant shear forces, defined as $Q={{\displaystyle }}_z{\tau }_{s}dz$. The deformations $\overline{\epsilon }$ are defined as $\overline{\epsilon }\equiv {\{\epsilon ;\kappa \}}^T$.

The bilinear operator for the Mindlin model is

$G(d,\delta d)={{\displaystyle }}_{\Omega }\delta {\overline{\epsilon }}^{T}C\overline{\epsilon }d\Omega +{{\displaystyle }}_{\Omega }\delta {\gamma }^{sT}{E}^s{\gamma }^{s}d\Omega$ (3.5)

where the last integral vanishes in the Kirchhoff case. Separating the linear and nonlinear membrane components, the bilinear form becomes

$G(d,\delta d)={\displaystyle {\int }_{\Omega }(\delta {\overline{\epsilon }}^{0^T}C{\overline{\epsilon }}^0+\delta {\overline{\epsilon }}^{0^T}C{\overline{\epsilon }}^{NL}+\delta {\overline{\epsilon }}^{N{L}^T}C{\overline{\epsilon }}^0+\delta {\overline{\epsilon }}^{N{L}^T}C{\overline{\epsilon }}^{NL})d\Omega }+{\displaystyle {\int }_{\Omega }\delta {\gamma }^{s^T}{E}^s{\gamma }^{s}d\Omega ,}$ (3.6)

where

${\overline{\epsilon }}^{NL}\equiv \left\{\begin{array}{c}{\epsilon }^{NL}\\ 0\end{array}\right\}and{\overline{\epsilon }}^0\equiv \left\{\begin{array}{c}{\epsilon }^0\\ \kappa \end{array}\right\}.$ (3.7)

1.3 Discretization for Kirchhoff and Mindlin plate models

For the Kirchhoff model the generalized displacement fields can be organized as $u(x)={\{u_0,v_0,w\}}^T$ and for Mindlin model as $u(x)={\{u_0,v_0,w,{\theta }_x,{\theta }_y\}}^T$. Each of these components is discretized with the enriched basis of GFEM over an arbitrary element e as

$u_i^e(x)={\displaystyle \sum }_{j=1}^{Nne}{\varphi }_{ij}(x,y)\left\{u_{ij}+{\displaystyle \sum }_{k=1}^{m_{ij}}{L}_{ijk}(x)b_{ijk}\right\},$ (3.8)

where $Nne$ is the number of nodes of element e, ${\varphi }_{ij}$ are the PoU functions associated to the node $j$, $i$ represent a component of $u$ and $k$ is the counter associated to the enrichment monomial functions. The limit $m_{ij}$ is the number of enrichment functions of the displacement $u_i^e$ associated to node $j$, and $L_{jk}$ are the enrichment functions of this same node. The terms $u_{ij}$ are nodal values of $u_i^e$ associated to the PoU functions and $b_{jk}$ are the coefficients of the enriched basis functions ${\varphi }_{ij}{L}_{ijk}$. A compact notation can be obtained collecting the nodal parameters $u_{ij}$ and $b_{jk}$ in a vector $U^e$ and the corresponding basis functions in a matrix $N^e$ of dimensions $q\times N_{dofe}$, where $q=3$ or $q=5$ for Kirchhoff or Mindlin models, respectively, and $N_{dofe}$ is the number of degrees of freedom in the element. Therefore, the discretized generalized displacements are represented in the usual FEM or GFEM way as

$u^e(x)=N^e(x)U^e.$ (3.9)

For future use, the transverse displacement is discretized as

$w(x)=N^w{U}^e.$ (3.10)

Substituting the displacement representation (3.9) in the strain-deformations (3.2), one obtains the generalized linear deformations in discretized form

${\overline{\epsilon }}^0=B{U}^{e}and{\gamma }^s=B^s{U}^e,$ (3.11)

where $B$ is the deformation matrix associated with linear membrane deformation and change of curvature, and $B^s$ is the shear deformation matrix, used in case of Mindlin model.

With regard to the nonlinear generalized deformation ${\overline{\epsilon }}^{NL}$, here we follow the development introduced by Zienkiewicz & Taylor (1991) for the von Kármán theory with Kirchhoff model, and perform the necessary adjustments to GFEM and the Mindlin model. First, one defines the matrices

$A^z={\left[\begin{array}{cccccc}\frac{\partial w}{\partial x}& 0& \frac{\partial w}{\partial y}& 0& 0& 0\\ 0& \frac{\partial w}{\partial y}& \frac{\partial w}{\partial y}& 0& 0& 0\end{array}\right]}^{T}andG=\left[\begin{array}{c}\frac{\partial N^w}{\partial x}\\ \frac{\partial N^w}{\partial y}\end{array}\right].$ (3.12)

Next, a vector $\theta$ of the small rotations of the straight segments normal to the reference surface is defined along with its variation

$\theta =\left\{\begin{array}{c}\frac{\partial w}{\partial x}\\ \frac{\partial w}{\partial y}\end{array}\right\}=GU,and\delta \theta =G\delta U,$ (3.13)

where $U$ is the global vector of displacement coefficients and $\delta$ is the variation operator. The nonlinear strain vector in (3.7) can be decomposed as a product of two linear terms as follows

${\overline{\epsilon }}^{NL}=\frac{1}{2}{A}^z\theta =\frac{1}{2}{A}^{z}GU.$ (3.14)

Its first variation can be written in a more convenient form for computational application after a sequence of operations:

$\delta {\overline{\epsilon }}^{NL}=\frac{1}{2}\delta A^z\theta +\frac{1}{2}{A}^z\delta \theta =A^z\delta \theta =A^{z}G\delta U.$ (3.15)

Defining $B^{NL}=A^{z}G$ in (3.14), we have the discretization of the deformation terms and their variations:

$\begin{array}{c}{\overline{\epsilon }}^{NL}=\frac{1}{2}{B}^{NL}(U)U,\delta {\overline{\epsilon }}^{NL}=B^{NL}(U)\delta U,\\ {\overline{\epsilon }}^0=BU,\delta {\overline{\epsilon }}^0=B\delta U,\end{array}$ (3.16)

supplemented by the transverse shear deformation discretization ${\gamma }^s=B^{s}U$ in case of Mindlin model.

The bilinear operator (3.6), which represents the internal virtual work $\delta W_{\mathrm{int}}$, then becomes

$G(u,\delta u)=\delta U^T{\displaystyle {\int }_{\Omega }\{B^{T}CB+B^{T}C\frac{1}{2}{B}^{NL}+{\left(B^{NL}\right)}^{T}CB+{\left(B^{NL}\right)}^{T}C\frac{1}{2}{B}^{NL}+B^{sT}{E}^s{B}^s\}d\Omega U.}$ (3.17)

Considering $G(u,\delta u)=\delta U^T{F}_{\mathrm{int}}$, where $F_{\mathrm{int}}=K(U)U$ is the nodal vector of internal forces, the stiffness matrix $K$ is

$K(U)={\displaystyle {\int }_{\Omega }[B^{T}CB+\frac{1}{2}{B}^{T}C{B}^{NL}+{\left(B^{NL}\right)}^{T}CB+\frac{1}{2}{\left(B^{NL}\right)}^{T}C{B}^{NL}+B^{sT}{E}^s{B}^s]d\Omega .}$ (3.18)

The equilibrium algebraic system $K(U)U=F_{ext}$ is rewritten in residual form

$R=F_{\mathrm{int}}(U)-F_{ext}.$ (3.19)

Following the usual procedure, the Newton-Raphson iterations are based on:

$K^T\Delta U^k=-R^k,$ (3.20)

where $\Delta U^k$ is a correction such that in $k$-th iteration the solution is approximated by $U^{k+1}=U^k+\Delta U^k$ . $K^T$ is the tangent matrix, defined as $K^T=\partial {\left(KU\right)}^k/\partial U$, which can be represented, in indicial form and using summation rule, as

$K_{ir}^T=\frac{\partial (K_{ij}{U}_j)}{\partial U_r}=\frac{\partial K_{ij}}{\partial U_r}{U}_j+K_{ir}^0,$ (3.21)

and $K_{ir}^0$ is the linear term, given by

$K^0={\displaystyle {\int }_{\Omega }\left[B^{T}CB+B^{sT}{E}^s{B}^s\right]d\Omega .}$ (3.22)

The first term comes from the differentiation of (3.18), which can be obtained taking $B^{NL}=A^{z}G,$ such that

$\frac{\partial B^{NL}}{\partial U}=\frac{\partial A^z}{\partial U}G.$ (3.23)

The derivative of $A^z$ requires the differentiation of w, obtained from (3.10) as $w=N_k^w{U}_k$, such that $\partial w/\partial U_r=N_r^w.$ Taking $A^z$ from (3.12), $\partial A^z/\partial U$ becomes

$\frac{\partial A_{kl}^z}{\partial U_r}={\left[\begin{array}{cccccc}{N}_r^w{,}_x& 0& N_r^w{,}_y& 0& 0& 0\\ 0& N_r^w{,}_y& N_r^w{,}_x& 0& 0& 0\end{array}\right]}^T.$ (3.24)

The symbol $A_{klr}\equiv \partial A_{kl}^z/\partial U_r$ is adopted to compact the notation. With the derivative of $B^{NL}$, the tangential matrix becomes

$K_{ir}^t=K_{ir}^0+A_{mkr}{G}_{ki}{C}_{mn}{B}_{nj}+B_{mi}{C}_{mn}\frac{1}{2}{A}_{nkr}{G}_{kj}+G_{km}{A}_{kir}{C}_{mn}\frac{1}{2}{A}_{nl}^z{G}_{lj}+G_{km}{A}_{ki}^z{C}_{mn}\frac{1}{2}{A}_{nlr}{G}_{lj}.$ (3.25)

The influence of the transverse shear deformations of the Mindlin model appears in the linear stiffness$K_{ir}^0$, the second term in the integral in (3.22), which is null in case of Kirchhoff model.

PLATE MODELS AND GFEM DESCRETIZATION

Here we investigate the behavior of the Generalized Finite Element approach to deal with the initial stability problem of the laminated plate associated with the von Kármán's theory. The weak form is the standard one, given by Brush and Almroth (1975).

${\displaystyle {\int }_{\Omega }\delta {\overline{\epsilon }}^{0T}C{\overline{\epsilon }}^{0}d\Omega }+{\displaystyle {\int }_{\Omega }\delta {\gamma }^{sT}E{\gamma }^{s}d\Omega }-\lambda {\displaystyle {\int }_{\Omega }\delta {\theta }^T\left[\begin{array}{cc}{N}_x& N_{xy}\\ N_{xy}& N_y\end{array}\right]\theta d\Omega =0.}$ (4.1)

Here, $\{N_x,N_y,N_{xy}\}$ form a prescribed profile of in-plane resultant stresses acting in the plate, and one searches for the smallest value $\lambda$ which makes the equality holds for arbitrary admissible weighting functions. Hence, the plate is considered initially undeformed and loaded only by the prescribed in-plane stress resultants. In this way, $\lambda \{N_x,N_y,N_{xy}\}$ gives an estimate for the initial bucking load of the plate, which is proven to be an upper bound of the real bucking load. After discretization, Eq. (4.1) results in the algebraic eigenvalue problem

$\left[K^0-\lambda K_G\right]U=\mathrm{0,}$ (4.2)

where $K^0$ is the linear stiffness matrix, $\lambda$ and $U$ are eigenvalues and eigenvectors and $K_G$ is the geometric matrix, real, symmetric and, in general, singular. One searches for the smallest eigenvalue associated to the critical buckling load. Both in FEM and in GFEM, $K_G$ has non zero entries only on the positions associated to transverse displacement $w$, that is, all rows and columns associated with $u$, $v$ (and ${\theta }_x$, and ${\theta }_y$ in case of the Mindlin model) are zero. In addition, even rows associated with $w$ can be null depending on the prescribed values of $N_x$, $N_y$ and $N_{xy}$. Since only one eigenvalue is needed, the smallest one, an interactive method is usually selected to search only for the first eigenpair, for example the method of subspace iteration or the Lanczos adapted to symmetric matrices. Then, it is sufficient to use a small penalty $\epsilon$ on the zero diagonals of $K_G$. These values generate a set of eigenvalues of order $1/\epsilon$, which are independent from the first ones, and will be the largest ones in the set. In case of FEM, both $K$ and $K_G$ become positive-definite matrices after the imposition of adequate Dirichlet boundary conditions.

In case of the usual C⁰-GFEM, it is proven that the function basis formed by the linear partition of unity uniformly enriched by polynomial functions is linearly dependent. This generates a singular stiffness matrix $K^0$. Besides the number of null eigenvalues associated with the number of rigid body motions of the model, there is an a-priori unknown number of null eigenvalues associated with the linear dependency of the basis. Considering the complete set of $N$ eigenvalues of the algebraic system, the first non-null, the one of interest, will be n -th one in the ordered list. Usually, n can be as large as about a half of $N$. The singularity of $K^0$ is dealt with in this paper by the use of the standard shift technique (solving the eigenproblem $\left[\overline{K}-\overline{\lambda }{K}_G\right]U=0$, with $\overline{K}=K+s{K}_G$ and $\overline{\lambda }=\lambda +s$, for arbitrary $s>0$). Now, $\overline{K}$ becomes positive-definite for C⁰-GFEM. However, the iterative method has to solve for a number of eigenpairs large enough to include the first non-zero one. This defeats the whole purpose of eigenvalue methods like Lanczos or subspace iterations, because the size of the subspace becomes too large to render the method efficient.

On the other hand, function basis formed by the smooth C${}^k$ partition of unity enriched by polynomial monomials is linearly independent (Mendonça et al., 2011). As a result, after imposition of sufficient Dirichlet boundary conditions, the stiffness matrix still has a large condition number but it is nonsingular. In the present paper the shift technique is used in both C⁰-GFEM and C${}^k$-GFEM.

NUMERICAL RESULTS

The present tests with the GFEM formulations were performed modeling the entire plate by a set of uniform meshes, defined by a mesh parameter M, some of which are illustrated in Figure 1. The C${}^k$-GFEM behavior in extremely distorted meshes has already been investigated in (Mendonça et al., 2011) and is not considered here. The same can be stated about benchmark problems with singularities, whose GFEM response have been already intensively investigated in the literature.

Figure 1
Examples of meshes used in the laminated plate numerical tests, with mesh indices M = 1, 2 and 4.

In transverse load cases, the boundaries are simply supported of the type $u_0=v_0=w=0$ along all edges. These constraints are imposed by selectively restricting the enrichment coefficients. For example, for a node on a boundary $x=const.$, the coefficients of functions whose enrichment contains $\overline{x}$ are not constrained, since these functions are naturally zero on that part of the boundary. In this way, the enriched function is preserved throughout the domain. For arbitrarily shaped domains, where the boundary segments are neither straight nor parallel to a given coordinate axis, the procedure proposed in Garcia et al. (2009) can be used, where the partition of unity associated to the boundary node is substituted by it multiplied by a ramp function which is zero at the boundary.

The numerical integration is a special concern when the smooth approximation functions in $C^k$-GFEM is used, because the partition of unity is not polynomial but of rational type. A study about the number of integration points required to integrate adequately the stiffness matrix is described in de Barcellos et al. (2009), using triangular and Gaussian integration rules. In all cases solved in the present paper, an indication of the number of integration points used, $nip$, is shown, and the triangular rule of Wandzura and Xiao (2003) is used.

In comparing results between C${}^0$- and C${}^k$-GFEM, we distinguish the degree $p$ of uniform refinement and the degree $b$ of the least complete polynomial the approximation basis is capable of reproducing. For C${}^0$-GFEM, the PoU is the linear tent and the enriched basis has a polynomial reproducibility of degree $b=p+1$. For C${}^k$-GFEM, the largest polynomial the PoU can reproduce is the constant, and as a consequence, $b=p$.

1.4 Homogeneous, isotropic, thin plate problem

The first evaluation of the behavior of $C^0$-GFEM and $C^k$-GFEM is the simplest case of a homogeneous isotropic square plate, clamped on all edges, modeled by the Mindlin kinematic model, loaded in the nonlinear range. The problem is interesting for a first evaluation because there is an analytic solution available for the thin plate model by Levy (1942). The data used here are the following: plate with sides $a=254$ mm, thickness $H=1.27$ mm (aspect ratio $a/H=200$) and elastic material properties $E=68.947$ GPa and $\nu =0.316$. The transverse distributed load is uniform, $q_0=0.0138$ MPa. Due to the symmetry, only one quarter of the plate is modeled. Considering the Cartesian axes $0xy$ with origin at the plate center, the boundary conditions for the Mindlin model are: at $x=a/2$ and $y=a/2$, fixed conditions $u_0=v_0=w={\theta }_x={\theta }_y=0$, and at $x=0$ and $y=0$, symmetry conditions $u_0={\theta }_x=0$ and $v_0={\theta }_v=0$, respectively.

Figure 2 shows the pointwise relative errors of the transverse displacement at plate center, $w_r=(w-w_{ex})/w_{ex}$, versus the number of dof, for $C^0$-and $C^0$-GFEM, for the different meshes M1 to M32, with uniform polynomial enrichment of degrees $p=\mathrm{0,...,4}$. The results are for the Mindlin model with shear parameter $k_s=1$. The reference value of displacement, from Levy (1942), is $w_{ex}=\mathrm{1,}7272$. In some graphs, the load is parametrized as

$\overline{P}=\frac{q_0{a}^4}{E{h}^4}.$ (5.1)

Figure 2:
Relative error of the transverse displacement at plate center $w_r$ versus dof, for (a) $C^0$-GFEM and (b) $C^k$-GFEM. nip =25, $\overline{P}$=320. Mindlin model with $k_s$=1.

The results in Figure 2 are obtained from $\overline{P}=320$. This load generates maximum transverse displacement of about $w/h\approx \mathrm{1,}7$, which puts the response beyond the usual linear range ($w/h=\mathrm{0,}5$), but still inside the range of validity of the von Kármán model. For the results of $C^k$-GFEM in Figure 2, for lower polynomial enrichments $p=1$and 2, the shape of the basis functions and their derivatives are considered simpler and easier to integrate, therefore we use $nip=25$ integration points. For higher enrichment degrees, $p=3$ and 4, we use $npi=54$ points per element, in order to guarantee a proper integration. The results show a tendency for locking with enrichments $p=0$ and 1, for $C^0$- and $C^k$-GFEM respectively, that is, for a basis of polynomial reproducibility $b=1$. For enrichments $p=2$ and higher, $C^0$-GFEM results show stabilization with mesh refinements, which can be caused by ill conditioning of the stiffness matrix. For the smooth GFEM, the curves show progressive convergence with $h$-refinement. The errors obtained when using mesh M1 are almost always high because this mesh is too crude, particularly in Figure 2 $C^k$-GFEM case.

A more systematic evaluation of the effect of number of integration points ($nip$) in the response of $C^k$-GFEM is described in Figure 3. Here, the relative error of transverse displacement at plate center, $w_r,$ is represented versus $nip$, with enrichments $p=1$ and 2. For each mesh only $nip=\mathrm{7,}23$ and 54 points are included, which are the smallest ones available in the integration rule of Wandzura. It is clear a tendency of convergence at 25 points for the present type of problems.

Figure 3:
Relative error of the transverse displacement at plate center $w_r$ versus the number of integration nip. Mindlin model with $C^k$-GFEM and enrichments (a)$p=1$ and (b) $p=2$.

1.5 Homogeneous, isotropic, thick plate problem

The behavior of $C^0$-GFEM and $C^k$-GFEM approximations for Mindlin model in moderately large displacements in a thick plate is evaluated considering the case of an isotropic square plate with simply supported edges ($u_0=v_0=w=0$) with sides $a=254$ mm and thickness $H=25.47$ mm, such that the aspect ratio is $a/H=10$. The material properties are $E=68.947$ GPa and $\nu =0.316$. The shear correction factor is set to $k_s=5/6$. Only a quarter of the plate is modeled due to the symmetry and the transverse load is uniformly distributed, $q_o$.

Figure 4 shows the pointwise normalized transverse displacement $\overline{w}$ and the normal stress ${\overline{\sigma }}_x$ versus transverse load $\overline{P}$, for $C^0$-GFEM, with polynomial uniform enrichment of degree $p=2$ and $nip=25$ integration points. The displacements are at the plate center and the reference solution is taken from Reddy (Reddy, 2004). This reference uses a finite element solution and the stress is obtained at the point $(x,y,z)=(6.25;6.25;-h/2)$, chosen to coincide with one integration point of his model. That model consists of an uniform mesh of 16 nonconforming quadrilateral elements, Q9. The results obtained from GFEM uses a mesh M8, with polynomial enrichment of degree $p=2$. It must be observed that the results from Reddy (2004) are included here only to illustrate results available in the literature and to show a general trend. Both models have different number of nodes and types of approximation functions. Likely, the present model is more refined and accurate. The values of displacement and normal stresses are listed in Table 1 for the normalized values defined in the Figure 4.

Figure 4:
(a) transverse displacement $\overline{w}$ and (b) normal stress ${\overline{\sigma }}_x$ versus transverse load $\overline{P}$, $C^0$-GFEM, enrichments $p=2$.and nip = 25 Mindlin model with $k_s$=5/6.

Figure 5 shows the evolution of the normalized in-plane stress ${\overline{\sigma }}_x$ along the line $y=0$, from the plate center $x=0$ to the edge, at four different z levels, obtained with $C^0$-GFEM with $p=2$, $\overline{P}=250$ and $nip=25$, and the mesh M8. Stresses obtained from the linear analysis are also shown, multiplied by 0.1. At this load level the nonlinear effect is pronounceable: the membrane stresses causes non zero stresses at the border $x=a/2$, differently from the usual zero value from the linear model; the laminate is symmetric and there is no coupling between bending and membrane behavior in the linear problem, and ${\sigma }_x$ is zero for all z at the border in order to comply with the zero bending moment there. Along the range, from the center to near the border, the stress variation departs from the usual sine like shape of the linear case to an almost constant curve, also because of the nonlinearly generated membrane forces. Most of the transverse load is supported by the transverse component of the membrane stresses instead of the bending moment associated with the transverse gradient of ${\sigma }_x$. This can be seen by the difference between the values of ${\overline{\sigma }}_x$ at the plate center for the linear analysis, ${\overline{\sigma }}_x\approx \pm 76.4$, and nonlinear analysis, ${\overline{\sigma }}_x\approx (1.69;19.6)$ at the top and bottom surfaces respectively.

Table 1
Normalized transverse displacement at plate center ($\overline{w}$) and normalized in-plane stress (${\overline{\sigma }}_x$) for loads ($\overline{P}$) with simply supported boundary condition. Mindlin model with k = 5/6

	Reddy (2004 )		$C^0$ -GFEM		$C^k$ -GFEM
$\overline{P}$	$\overline{w}$	${\overline{\sigma }}_x$	$\overline{w}$	${\overline{\sigma }}_x$	$\overline{w}$	${\overline{\sigma }}_x$
6.25	0.2790	1.861	0.2790	1.872	0.2668	1.800
12.5	0.4630	3.305	0.4630	3.327	0.4487	3.236
25	0.6911	5.319	0.6910	5.360	0.6771	5.264
50	0.9575	8.001	0.9574	8.067	0.9447	7.962
75	1.1333	9.983	1.1332	10.069	1.1211	9.952
100	1.2688	11.634	1.2687	11.736	1.2568	11.606
125	1.3809	13.085	1.3808	13.201	1.3691	13.058
150	1.4774	14.398	1.4773	14.529	1.4658	14.373
175	1.5628	15.608	1.5627	15.756	1.5512	15.585
200	1.6398	16.743	1.6398	16.903	1.6283	16.718
225	1.7102	17.812	1.7102	17.986	1.6987	17.788
250	1.7752	18.831	1.7752	19.017	1.7637	18.805

Figure 5:
Normal stress ${\overline{\sigma }}_x$ along the line $y=0$, at four different $z$ levels at the thickness, for $C^0$-GFEM with $p=2$ and $\overline{P}=250$. (The stresses of linear response are multiplied by 0.1). Mindlin model with $k_s$=5/6.

1.6 Cross-ply laminate problem

The problem considered here consists of a square laminated plate with sides $a=b=200$ mm aligned along $x$ and $y$ axis and thickness $h=2$ mm, with three equal orthotropic layers with orientations $[0^{\circ }/{90}^{\circ }/0^{\circ }]$. Each layer has the following orthotropic properties: $E_1=25E_2$, $G_{12}=G_{13}=0.5E_2$, ${\nu }_{12}={\nu }_{23}=0.25$, $G_{23}=0.2E_2$ with $E_2=7$ GPa. The laminate is thin, with aspect $a/h=100$. The boundary conditions are simply supported: at $x$ and $y=\pm a/2$, $u^0=v^0=w=0$. The tests are made with uniform transverse distributed loads ranging from $0.1$to$1$MPa.

Since the laminate is relatively thin, with aspect $a/h=100$, the behavior of the Kirchhoff model can be evaluated using conform Kirchhoff GFEM formulation with $C^k$ smooth functions. In this way, Figure 6 shows the transverse displacement $\overline{w}$ in the plate center versus the transverse load, using FSDT (Mindlin) with both $C^0$ and $C^k$-GFEM and CPT (Kirchhoff) plate model using $C^k-$GFEM. In all cases, the approximation basis has reproducibility polynomial degree $b=3$.

Figure 6:
Cross-ply laminate. Transverse displacement $\overline{w}$ in the plate center node, using FSDT and CPT plate models, $C^0$- and $C^k$-GFEM with b = 3, a/h = 100.

Figure 7 shows the normal stress ${\overline{\sigma }}_x$ in the central node at top and bottom of the plate, versus transverse load, using FSDT and CPT plate models. The results are compared between the $C^0$ and $C^k$-GFEM with $b=3$. The nonlinear results for FSDT model obtained with $C^0$ and $C^k$-GFEM are very similar to each other, for the stresses at both surfaces. Figure 8 shows the normal stress ${\overline{\sigma }}_x$ across thickness at the center node, computed at two neighboring elements. Because the results are obtained with $C^0$-GFEM, the discontinuity is expected, but it is very small in this case.

Figure 7:
Normal stress ${\overline{\sigma }}_x$ in the central node at top and bottom of the plate, using FSDT and CPT plate models, $C^0$- and $C^k$-GFEM with b = 3. “b” and “t” refer to results at bottom and top surfaces of the laminate, respectively.

Figure 8:
Normal stress ${\overline{\sigma }}_x$ through the thickness at the interfaces elements $e^{1}and{e}^2$ shearing the central node, using Mindlin model, $C^0$-GFEM with $p=3$ and $nip=\mathrm{250,}q=q_{\max }=1.$

The results for displacements and stresses for Mindlin and Kirchhoff are naturally similar because the laminate is thin. Also, the results from $C^0$ and $C^k$-GFEM are similar for the same degree $b$ of the approximation basis. For Mindlin model, a given $b$ implies more degrees of freedom per node in the smooth than $C^0$-GFEM, as seen in Table 2. However, if the laminate is thin, $C^k$-GFEM in Kirchhoff model requires equal or less dof's/node than $C^0$-GFEM to attain a given degree $b$ in the basis, for orders of $b=3$ or higher.

Table 2
Number of degree of freedom per node in Mindlin and Kirchhoff models with $C^0$ and $C^k$-GFEM

Mindlin					Kirchhoff
	$C^0$-GFEM		$C^k$-GFEM		$C^k$-GFEM
	$u_0,v_0,w,{\theta }_x,{\theta }_y$		$u_0,v_0,w,{\theta }_x,{\theta }_y$		$u_0,v_0,w$
$b$	$p$	$dof/node$	$p$	$dof/node$	$p$	$dof/node$
1	0	5x1 = 5	1	5x3 = 15	1	3x3 = 9
2	1	5x3 = 15	2	5x6 = 30	2	3x6 = 18
3	2	5x6 = 30	3	5x10 = 50	3	3x10 = 30
4	3	5x10 = 50	4	5x15 = 75	4	3x15 = 45
5	4	5x15 = 75	5	5x21 = 105	5	3x21 = 63

1.7 Instability analysis

The problem considered here consists of a square laminated plate with equal sides $a=200$ mm aligned along $x$ and $y$ axis, with thickness $H=2$ mm, simply supported at all sides, with four equal orthotropic layers with orientations $[0^{\circ }/{90}^{\circ }/{90}^{\circ }/0^{\circ }]$. Each layer has the following properties in their orthotropic directions: $E_1=25E_2$, $G_{12}=G_{13}=0.5E_2$, ${\nu }_{12}={\nu }_{23}=0.25$, $G_{23}=0.2E_2$ with $E_2=7$ GPa. These geometric and material relations allow analytic solution of the critical load for the Kirchhoff model. Considering simply supported boundary conditions with $N_x=N_y=1$ and $N_{xy}=0$, the first exact eigenvalue is (Whitney, 1987)

$\lambda =\frac{{\pi }^2{D}_{22}}{a^2(1+n^2)}\left[\frac{D_{11}}{D_{22}}+\frac{D}{D_{22}}{n}^2+n^4\right],$ (5.2)

where $D=D_{12}+2D_{66}$. The buckling mode is $\sin \left(n\pi x/a\right)\cos \left(m\pi y/a\right)$, where, $m=1$ and n must be chosen to provide the minimum value of $\lambda$. For these layers properties, the laminate bending elastic matrix has nonzero components $D_{11}=1.029240\cdot {10}^5$ Nmm, $D_{22}=1.871345\cdot {10}^3$ Nmm, $D_{12}=1.169591\cdot {10}^3$ Nmm, $D_{66}=2.333333\cdot {10}^3$ Nmm. Then, $D=5836.24$ Nmm and the critical load happens with mode $n=1$, with value $\lambda =16.446457757$. Where necessary, a reference value for the buckling load associated with the Mindlin model was taken from an uniform mesh of 9 nodes Lagrangian quadrilateral finite elements. It was used a sequence of meshes until of $50\times 50$ elements, with 10201 nodes. The shear correction factor was set to ${\kappa }_s=1.0$. The obtained critical value is $\lambda =\mathrm{16,}36971$.

Figure 9 and Figure 10 display results of relative errors of $\lambda$ versus number of degrees of freedom for Kirchhoff and Mindlin models respectively, with C${}^k$-GFEM. Each curve is associated to a given mesh, as indicated. The dots indicate the enrichment degree $p$ of the basis, starting with $p=1$ at the top, except for the first curve, mesh $M=1$, which starts from $p=2$. In this case, the discretization was too poor to allow a solution. We see the errors and the rates of convergence are similar for both kinematic models. All results were obtained with penalty $\epsilon$ on $K_G$ equal to ${10}^{-2}$ times the smallest non zero value at the diagonal in $K_G$. The shift factor was $s=10$. The eigenvalues were obtained with the subspace iteration method, with error tolerance $tol={10}^{-8}$ for the first 2 eigenvalues. The size of the iteration subspace was set to 10 eigenvectors.

Figure 9:
Relative error of $\lambda$ for Kirchhoff model versus number of degrees of freedom.

Figure 10:
Relative error of $\lambda$ for Mindlin model with $C^k$-GFEM versus number of degrees of freedom.

The abscissa scale in both Figure 9 and Figure 10 are the same. It can be seen a horizontal shift of the curves to the left in the Kirchhoff results, with respect to those of Mindlin model. This is natural because, for the first model, there is only one displacement component (w) and for the first order model there are three ($w,{\theta }_x$ and ${\theta }_y$).

For the Kirchhoff model with enrichment $p=1$ there is no convergence. This can be understood considering that only the transverse displacement is modeled. With $p=1$, the basis is too poor for allowing the adequate derivatives $w{,}_x$ and $w{}_{,y}$ necessary to the model, although the space of approximation is still $C^k$ continuous, with arbitrary $k$. In the case of C${}^0$-GFEM, all displacement components, $w,{\theta }_x$ and ${\theta }_y$, are modeled independently and convergence is observed.

For both, Mindlin and Kirchhoff models, Figure 9 and Figure 10 show very low convergence for enrichment degree $p=4$. This can be imputed to numerical difficulties due to the condition number of the stiffness matrix.

For the tests with C${}^0$-GFEM, as commented in the last section, the stiffness matrix contains an a-priori unknown number of spurious solutions with zero or very small eigenvalues. Figure 11 shows the number $N_s$ of zero eigenvalues obtained for each degree $b$ of polynomial reproducibility of the basis, for Mindlin model. In all cases the $n=N_s+1$-th eigenvalue is an adequate approximation for the reference value of the critical load. The fraction of $N_s$ with the total number of degrees of freedom, $N_s/$dof, grows approximately linearly with the degree $b$ of the basis for all meshes, except for the coarsest one.

Figure 11:
Number of null eigenvalues for Mindlin’s model in $C^0$-GFEM versus degree $b$ of the polynomial basis.

1.7.1 Comparison between GFEM, stable GFEM and C${}^k$-GFEM on the stability analysis

The effects of the use of stabilization in C${}^0$-GFEM is evaluated using several meshes on the entire plate, in the same problem described at the beginning of this section, with $N_x=N_y=1$ and $N_{xy}=0$. Table 3 shows the estimated critical loads for the Mindlin model for polynomial enrichments $p=\mathrm{0,}\mathrm{1,}\mathrm{2,}3$ and 4, for both non-stabilized and stabilized C${}^0$-GFEM. The stabilization procedure (2.7), when applied to polynomial enrichment functions, automatically renders the stabilized linear function identically zero, such that remains in the approximation basis only the partition of unity function and the functions enriched by polynomials of degrees $p=2$ and greater. Except for the very coarse mesh M1 (two elements and four nodes), C${}^0$-GFEM with $p=2$ gives reasonable approximations for the critical load in all meshes, however, the first non-zero eigenvalue is the $n$-th one of the algebraic system. For $p=\mathrm{2,}\mathrm{3,}\cdots ,$ the critical load is the 9th eigenvalue. In the stabilized case, C${}^0$-SGFEM, the critical load is the first eigenvalue for all $p$'s, and the more refined meshes, M4, M6 and M8. For the coarser meshes M1 and M2, the critical load is the first eigenvalue when $p=2$, but can be the 13th one for M2 and $p=4.$ For the refined meshes M4, M6 and M8, and $p=\mathrm{1,}$ the SGFEM critical load is the first eigenvalue but the approximation is too crude to be useful. This is in opposition to the non-stable GFEM, which gives reasonable approximation. For example, for M6 and $p=1$, we have $\lambda =17.306$ and $308.5$ for GFEM and SGFEM respectively. For M2, only with $p=4$ one obtains a reasonable SGFEM approximation, but the algebraic system contains 12 numerically zero eigenvalues. As described in Babuška and Banerjee (Babuška and Banerjee, 2012, 2011), the stabilization procedure aims to improve the condition number and make its rate of growth reduced to the levels of the standard finite element method. It is seems that the stabilization procedure impoverish the approximation basis and does not preclude completely all spurious modes.

Table 3
Critical loads for $C^0$-GFEM and SGFEM versus degree of enrichment p, Mindlin model, nsi is the number of subspace iterations, d1 an d2 are the number of degree of freedom in the model. n is the first useful eigenvalue.

		GFEM					SGFEM
Mesh	$p$	$\lambda$	$n$	d1	nsi	$\lambda$	$n$	d2	nsi	d2/d1
M1	2	91.04	9	72	5	91.04	1	48	2	0.667
	3	17.0403	2	120	3	17.04031	2	96	3	0.8
	4	16.8133	15	180	16	16.81391	2	156	4	0.867
M2	1	27.79	2	81	3
	2	16.56	9	162	9	79.6	1	108	2	0.667
	3	16.48797	25	270	17	24.04	1	216	9	0.8
	4	16.36976	49	405	19	16.36977	13	351	20	0.867
M4	0	646.7	1	75	2
	1	19.063	2	225	5	646.7	1	75	2	0.333
	2	16.43064	9	450	10	47.72	1	300	6	0.667
	3	16.37024	25	750	12	16.74418	1	600	9	0.8
	4	16.37963	49	1125	21	16.37188	1	975	9	0.867
M6	0	308.5	1	147	11
	1	17.3055	3	441	7	308.5	1	147	11	0.333
	2	16.37819	9	882	9	25.18	1	588	6	0.667
	3	16.36966	25	1470	40	16.43857	1	1176	6	0.8
	4	16.36962	49	2205	41	16.37000	1	1911	7	0.867
M8	0	182.8	1	243	14
	1	16.7808	3	729	14	182.8	1	243	14	0.333
	2	16.37153	9	1458	16	19.434	1	972	6	0.667
	3	16.36961	25	2430	16	16.38854	1	1944	8	0.8
	4	16.36962	49	3645	40	16.36968	1	3159	8	0.867

An interesting side effect of the stabilization procedure is the reduction of the number of degrees of freedom for a given polynomial degree of enrichment. For an approximation basis capable of reproducing a polynomial of degree $b=3$, for example, in a two-dimensional space, C${}^0$-GFEM requires 6 functions per node while C${}^0$-SGFEM requires 4. This happens because the PoU function can represent degrees zero and one, and the uniform enrichment with quadratic monomials generates a set of 4 functions per node whose linear combinations are capable of representing all polynomials of degrees 3, 2 and 1. This is summarized in Table 4, where nf is the necessary number of monomials per node. The counting of number of degrees of freedom is given by 3 generalized directions in bending ($w,$ ${\theta }_x,$ ${\theta }_y$), times the number of nodes, times nf. The results are given as d1 and d2 in the Table 3 for GFEM and SGFEM respectively. The difference between both GFEM's gradually reduces with the growth in $p$, but for $p=\mathrm{2,}$ SGFEM needs only 2/3 of the number of dof in GFEM. As for the smooth GFEM, $b=p$, and the number nf of dof's in each generalized direction is the same as in C${}^0$-GFEM. It requires more degrees of freedom to attain a given polynomial reproducibility, but always produces the solution in the first eigenvalue of the algebraic system, without the need of stabilizing the enrichment monomials.

Table 4
Size of the basis nf for GFEM and SGFEM versus degree of enrichment p in two dimensions.

	$C^0$-GFEM		$C^0$-SGFEM		$C^k$-GFEM
	$u_0,v_0,w,{\theta }_x,{\theta }_y$		$u_0,v_0,w,{\theta }_x,{\theta }_y$		$u_0,v_0,w$
$p$	$p$	nf	$p$	nf	$p$	nf
0	1	1	1	1	-	-
1	2	3			1	3
2	3	6	3	4	2	6
3	4	10	4	8	3	10
4	5	15	5	13	4	15

Table 5
Relative errors of critical loads obtained from $C^0$-GFEM, SGFEM and $C^k$-GFEM versus degree of enrichment p, for several meshes. Mindlin model. nsi is the number of subspace iterations, d1, d2 and d3 are the number of degree of freedom in the model.

		$C^0$-GFEM		SGFEM		$C^k$-GFEM
Mesh	$p$	d1	Error	d2	Error	d3	Error
M1	2	72	4.5	48	4.6	72	4.6
	3	120	4.1 E-2	96	4.1 E-2	120	2.0 E-1
	4	180	2.7 E-2	156	2.7 E-2	180	1.2 E-2
M2	1	81	7.0 E-1			81	7.9 E+1
	2	162	1.2 E-2	108	3.9	162	1.1
	3	270	7.2 E-3	216	4.7 E-1	270	2.5 E-2
	4	405	3.2 E-6	351	3.8 E-6	405	5.0 E-4
M4	0	75	3.9 E+1
	1	225	1.6 E-1	75	3.9 E+1	225	1.9 E+1
	2	450	3.7 E-3	300	1.9	450	2.9 E-1
	3	750	3.2 E-5	600	2.3 E-2	750	2.2 E-3
	4	1125	6.1 E-4	975	1.3 E-4	1125	3.3 E-6
M6	0	147	1.8 E+1
	1	441	5.7 E-2	147	1.8 E+1	441	1.2 E+1
	2	882	5.2 E-4	588	5.4 E-1	882	8.4 E-2
	3	1470	3.0 E-6	1176	4.2 E-3	1470	3.3 E-4
	4	2205	5.4 E-6	1911	1.8 E-5	2205	6.9 E-6
M8	0	243	1.0 E+1
	1	729	2.5 E-2	243	1.0 E+1	729	7.9
	2	1458	1.1 E-4	972	1.9 E-1	1458	3.1 E-2
	3	2430	6.0 E-6	1944	1.2 E-3	2430	7.3 E-5
	4	3645	5.4 E-6	3159	1.8 E-6	3645	7.6 E-6

Figure 12:
Relative errors in critical load versus number of degree of freedom for Mindlin model, $C^0$-GFEM and SGFEM and $C^k$-GFEM. Every dot in a curve corresponds to a given degree $p$ of polynomial enrichment, starting from 0, 1 or 2, as indicated.

The results in Table 3 are reorganized in Table 5 to show relative errors in critical load. Here the errors are given versus number of degrees of freedom for Mindlin model, using C${}^0$-GFEM, SGFEM and C${}^k$-GFEM, for several meshes. To easy the interpretations, the results for two meshes, M2 and M8, were plotted in Figure 12. Every dot in a curve corresponds to a given degree $p$ of polynomial enrichment, starting from 0, 1,2,3 or 4, as indicated in the figure. The relation between the results of C${}^0$ and C${}^k$-GFEM follows as expected, because the difference between $b=p+1$ and $b=p$, respectively, generates a vertical translation between the curves, accompanied by a difference in convergence rates. However, it must be considered that the results for C${}^k$-GFEM were obtained computing only the first eigenvalue of the algebraic system, while for C${}^0$-GFEM the quantity $N_s$ of non-physical eigenvalues is unpredictable. Therefore, contrary to the tendencies of the convergence curves alone, C${}^k$-GFEM presents superior behavior in comparison to C${}^0$-GFEM for this type of application.

For C${}^0$-GFEM with mesh M8, it can be seen a lack of convergence between $p=3$ and 4. This is imputed to the degradation of the condition number, which becomes more pronounced than in the mesh M2. The C${}^k$-GFEM results for M8 does not show an effect so intense, because, by construction, its approximation basis is linearly independent, which renders the coefficient matrix still ill conditioned but nonsingular.

The SGFEM aims at a stabilization of the condition number growth, by elimination of the linear dependency in the approximation basis by removing the linear content from the enrichment monomials. In most cases, as seen in n in Table 3, this is enough to renders a numerically nonsingular coefficient matrix, giving the critical load as the first eigenvalue. Even though, for coarse meshes M1 and M2, for $p=4$, it can be seen in the results a list of eigenvalues several orders smaller than the physical one. This happens as a result of the ill conditioning of the matrix, even after the stabilization procedure.

One side effect of the stabilization is the impoverishment of the approximation basis, as seen comparing the curves C${}^0$-GFEM and SGFEM in Figure 12. This has been observed with extensive numerical experimentation by Li (2014) in linear plane elasticity problems. It seems that the linear contents in the enrichment monomials, although generating linear dependency in the basis, provides also a useful enrichment, which is eliminated in the stabilization procedure. Judging by the curves in Figure 12, the resulting stabilized function basis is capable of spanning a space function similar to the one of smooth C${}^k$-GFEM: the curves are very similar in the range $p=2$ and 3.

CONCLUSION

The current paper presents a comparison between two types of GFEM, $C^0$ and $C^k$, in response of Mindlin and Kirchhoff kinematic models under moderately large displacements, modeled under the von Kármán theory. In addition, the response of both stable and non-stable $C^0$-GFEM are evaluated in comparison with the smooth GFEM in the approximation of the buckling load, with regard to the effects of the stiffness matrix rank on the computational effort to obtain the first useful eigenvalue. The results allow the following observations and conclusions:

• 1) For Mindlin model, a given degree $b$ of the polynomial reproducibility of the approximation basis implies more degrees of freedom per node in the smooth than $C^0$-GFEM. However, if the laminate is thin, $C^k$-GFEM allows the use of Kirchhoff model, which requires equal or less dof's/node than $C^0$-GFEM to attain a given degree $b$ in the basis, for orders $b=3$ or higher, thus compensating the natural excess of dof's/node with respect to the $C^0$ version.

• 2) The convergence curves in pointwise relative errors of displacement for Mindlin model with $C^0$-GFEM show stabilization of the curves for all enrichments, after some number of dof's, differently from the curves of smooth GFEM, which keep the convergence for all enrichment degrees, except for approximation basis of degree $b=1$. In fact, both forms of GFEM show locking in the nonlinear response for the Mindlin model with $b=1$. As for the stabilization of the curves for the $C^0$-GFEM, it can be speculated it is related with bad conditioning of the stiffness matrix. The smooth version has the matrix naturally nonsingular, because the basis functions are linearly independent.

• 3) As in the linear problem, the von Kármán response with the smooth GFEM requires a larger effort in the numerical integration than in $C^0$-GFEM. This is due to the highly oscillatory aspect of the smooth enriched functions and their derivatives.

For the standard problem of the approximation of the upper bound of critical instability load by eigenvalue problem with $C^k$-GFEM, we observe the following:

• 4) The convergence curves of the critical load by $C^k$-GFEM and Kirchhoff model show lack of convergence for $p=1$, as expected because only the displacement component $w$ is modeled with linear approximation functions and its first derivatives are too poor. For the Mindlin model, $p=1$ generates a good convergence curve: the approximation basis is of degree $b=2$, applied over the three displacement components, $w$, ${\theta }_x$ and ${\theta }_y$. However, for thin laminates, the Kirchhoff model requires only one third of the dof's in the Mindlin model.

• 5) The behavior of the curves for enrichments $p=2$ and 3 are excellent for both Kirchhoff and Mindlin models. For $p=4$, the response is not convergent, raising the hypothesis of ill conditioning in the coefficient matrix.

Considering the critical load obtained by $C^0$-GFEM and SGFEM, we observe the following:

• 6) The critical loads obtained with $C^0$-GFEM are, in general, as accurate as with smooth GFEM, but it corresponds to the n -th eigenvalue of the problem, where n is a priori unknown. This is different from the$C^k$-GFEM, where all critical loads are the first eigenvalue of the model.

• 7) The stable form of GFEM, for non-coarse meshes, furnishes the critical load at the first eigenvalue, although with inferior accuracy. The processes of stabilization also impoverish the function basis. A positive aspect is that the number of dof's/node is reduced with respect to the non-stable GFEM because, in the process, the enrichment function of degree $p=1$ naturally vanishes from the stabilized basis. However, with the impoverishment of the basis, the SGFEM requires more enrichment than GFEM to attain a given error. The increase in the number of dof's in SGFEM with respect to $C^0$-GFEM, to obtain a given accuracy, is approximately the same as the increase in $C^k$-GFEM, that is, both $C^k$-GFEM and SGFEM show close convergence curves and produce the critical load at the first eigenvalue of the problem.

Acknowledgements

Paulo T. R. Mendonça, Marx Ribeiro and Clovis S. de Barcellos gratefully acknowledge the financial support provided by the Brazilian government agency for scientific research, CNPq (National Council for Scientific and Technological Development) for this research, under research grants 163.461/2012-0, 303.315/2009-1 and 304.698/2013-0, respectively.

References

• Babuška, I., Banerjee, U., 2012. Stable Generalized Finite Element Method (SGFEM). Comput. Methods Appl. Mech. Eng. 201-204, 91-111.
• Babuška, I., Banerjee, U., 2011. Stable Generalized Finite Element Method (SGFEM), ICES Report (11-07), The Institute for Computational Engineering and Sciences, The University of Texas at Austin. Austin, TX,USA.
• Babuška, I., Melenk, J.M., 1997. The Partition of Unity Method. Int. J. Numer. Methods Eng. 40, 727-758.
• Belytschko, T., Black, T., 1999. Elastic crack growth in finite elements with minimal remeshing. Int. J. Numer. Methods Eng. 45, 601-620.
• Belytschko, T., Gracie, R., Ventura, G., 2009. A review of extended/generalized finite element methods for material modeling. Model. Simul. Mater. Sci. Eng. 17(4), 043001 (pp. 24).
• Brush, D.O., Almroth, B.O., 1975. Buckling of Bars, Plates, and Shells. McGraw-Hill Kogakusha. New York, NY, USA.
• de Barcellos, C.S., de Tarso R. Mendonça, P., Duarte, C.A., 2009. A C-k continuous generalized finite element formulation applied to laminated Kirchhoff plate model. Comput. Mech. 44, 377-393.
• Dolbow, J., Moës, N., Belytschko, T., 2000. Discontinuous enrichment in finite elements with a partition of unity method. Finite Elem. Anal. Des. 36, 235-260.
• Duarte, C., Oden, J.T., 1995. A review of some meshless methods to solve partial differential equations, TICAM Report (95-06), The University of Texas at Austin. Texas Institute for Computational and Applied Mathematics Austin, TX.
• Duarte, C.A., Kim, D.J., Quaresma, D.M., 2006. Arbitrarily smooth generalized finite element approximations. Comput. Methods Appl. Mech. Eng. 196, 33-56.
• Duarte, C.A., Oden, J.T., 1996a. H-p clouds - an h-p meshless method. Numer. Methods Partial Differ. Equ. 12, 673-705.
• Duarte, C.A., Oden, J.T., 1996b. An h-p adaptive method using clouds. Comput. Methods Appl. Mech. Eng. 139, 237-262.
• Edwards, H.C., 1996. C-inf Finite Element Basis Functions, TICAM Report (96-45), The University of Texas at Austin.Austin, TX, USA.
• Garcia, O.A., Fancello, E.A., de Tarso R. Mendonça, P., 2009. Developments in the application of the generalized finite element method to thick shell problems. Comput. Mech. 44, 669-682.
• Gupta, V., Duarte, C.A., Babuška, I., Banerjee, U., 2015. Stable GFEM (SGFEM): Improved conditioning and accuracy of GFEM/XFEM for three-dimensional fracture mechanics. Comput. Methods Appl. Mech. Eng. 289, 355-386.
• Gupta, V., Duarte, C.A., Babuška, I., Banerjee, U., 2013. A stable and optimally convergent generalized FEM (SGFEM) for linear elastic fracture mechanics. Comput. Methods Appl. Mech. Eng. 266, 23-39.
• Levy, S., 1942. Square plate with clamped edges under normal pressure producing large deflection, NACA Tech. Rep. 740. Washington, DC, United States.
• Li, H., 2014. Investigation of stability and accuracy of high order generalized finite element methods. University of Illinois at Urbana-Champaign.
• Melenk, J.M., Babuška, I., 1996. The partition of unity finite element method: Basic theory and applications. Comput. Methods Appl. Mech. Eng. 139, 289-314.
• Mendonça, P. de T.R., de Barcellos, C.S., Torres, D.A.F., 2013. Robust generalized FEM approximations for higher-order conformity requirements: Application to Reddy’s HSDT model for anisotropic laminated plates. Compos. Struct. 96, 332-345.
• Mendonça, P.D.T.R., de Barcellos, C.S., Torres, D.A.F., 2011. Analysis of anisotropic Mindlin plate model by continuous and non-continuous GFEM. Finite Elem. Anal. Des. 47, 698-717.
• Moës, N., Dolbow, J., Belytschko, T., 1999. A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 46, 131-150.
• Oden, J.T., Duarte, C.A.M., Zienkiewicz, O.C., 1998. A new cloud-based hp finite element method. Comput. Methods Appl. Mech. Eng. 153, 117-126.
• Oden, J.T., Reddy, J.N., 1976. An introduction to the mathematical theory of finite elements. John Willey & Sons, Inc. New York, NY, USA.
• Reddy, J.N., 2004. An Introduction to Nonlinear Finite Element Analysis: with applications to heat transfer, fluid mechanics, and solid mechanics. OUP Oxford.Oxford, UK.
• Reddy, J.N., 1989. On refined computational models of composite laminates. Int. J. Numer. Methods Eng. 27, 361-382.
• Reddy, J.N., 1984. A refined nonlinear theory of plates with transverse shear deformation. Int. J. Solids Struct. 20, 881-896.
• Shepard, D., 1968. A two-dimensional interpolation function for irregularly-spaced data, in: Proceedings of the 1968 23rd ACM National Conference On -. ACM Press, New York, New York, USA, pp. 517-524.
• Sukumar, N., Moës, N., Moran, B., Belytschko, T., 2000. Extended finite element method for three-dimensional crack modelling. Int. J. Numer. Methods Eng. 48, 1549-1570.
• Ventsel, E., Krauthammer, T., 2001. Thin Plates and Shells: Theory, Analysis, and Applications. Marcel Dekker, Inc., New York.
• Wandzura, S., Xiao, H., 2003. Symmetric quadrature rules on a triangle. Comput. Math. with Appl. 45, 1829-1840.
• Whitney, J.M., 1987. Structural analysis of laminated anisotropic plates. Technomic Publishing Company, Inc., Lancaster, Pennsylvania.
• Zienkiewicz, O.C., Taylor, R.L., 1991. The finite element method. Vol. 2 Solid and fluid mechanics, dynamics and non-linearity, 4th ed. McGraw-Hill.London,UK.

Edited by

Editor:

Rogério José Marczak.

Publication Dates

• Publication in this collection
  28 Oct 2019
• Date of issue
  2019

History

• Received
  19 Nov 2018
• Reviewed
  29 July 2019
• Accepted
  04 Sept 2019
• Published
  11 Sept 2019