1 INTRODUCTION

It is surprisingly difficult to design a simple three-node triangular plate finite element based on Kirchhoff plate theory that can represent states of constant curvature and twist (completeness) and give good results when subjected to standard tests (^{Cook et al., 2002}). The DKT (discrete Kirchhoff triangle) element was first published by ^{Stricklin et al. (1969}) as an attempt to bridge this gap. The starting point of this element development is to assume quadratic interpolation for the normal rotations _{
βx
} and _{
βy
} within the element and Hermite cubic interpolation independently for the transverse displacement *w* along the element sides (^{Batoz et al., 1980}). Enforcement of the Kirchhoff constraint at selected points allows converting the approximations for _{
βx
} and _{
βy
} into ones interpolated from the three vertex values of *w*, _{
βx
} and _{
βy
} . The bending stiffness matrix is then obtained from this modified approximations for the rotations using the strain energy of Reissner-Mindlin plates without the contribution of the transverse shear strain. Approximation for *w* within the element is thus neither stated nor required, and the solutions obtained converge to those corresponding to Kirchhoff plates. Elimination of the shear strain contribution also eliminates any possibility of transverse shear locking.

A clear and relatively simple account of the element formulation is given by ^{Batoz et al. (1980}), which concluded that the element is one of the most efficient and reliable three-node triangular elements for Kirchhoff plates. ^{Batoz (1982}) also obtained an explicit expression for the element bending stiffness matrix in local coordinates, which was subsequently coded in global coordinates (^{Jeyachandrabose et al., 1985}). A theoretical study of the element convergence rate can be found in ^{Kikuchi (1975}) and ^{Bernadou (1996}).

Since the insightful work of ^{Batoz et al. (1980}), the DKT element has become attractive and more developments and applications using the element have been reported. It is considered in several papers dealing with linear and nonlinear analysis of general shell structures (^{Bathe and Ho, 1981}; ^{Cook, 1993}; ^{Khosravi et al., 2007}; ^{Zhang et al., 2015}). High-speed component impact analysis (^{Wu et al., 2005}) demonstrates, for instance, that this element performs as good as the main element (^{Belytschko et al., 1984}) used for crashworthiness analyses. Satisfaction of the Kirchhoff constraint at discrete points is adopted in the development of all thin shell elements of the commercial finite element program ABAQUS. In particular, the DKT element is used to formulate the flat shell element STRI3 (^{ABAQUS, 2012}).

As far as the geometric stiffness matrix is concerned, inconsistent forms have been often adopted by assuming ad hoc shape functions for *w* within the element (^{Bathe and Ho, 1981}; ^{Mateus et al., 1997}; ^{Doyle, 2001}; ^{Khosravi et al., 2007}). If the modified approximations obtained for _{
βx
} and _{
βy
} after enforcing the Kirchhoff constraint at selected points are used to generate the geometric stiffness matrix, such a matrix is said to be consistent because their shape functions are the same used in establishing the bending stiffness matrix. Although the consistent version of the matrix can be found in the literature (^{Garnet and Pifko, 1983}), its explicit expression seems to be stated here for the first time. In addition to be consistent, our version exhibits the attractive features of avoiding the expensive exact integration by seven-point Gauss quadrature and of being directly stated in global coordinates. Numerical results are reported on linear buckling of plates to compare the element performance either with the proposed explicit consistent matrix, or with the most commonly used inconsistent matrix obtained by assuming a linear approximation for *w* within the element.

2 LINEARIZED BUCKLING FORMULATION

In the linearized buckling analysis of a Reissner-Mindlin plate, the potential energy ∏ = _{
Ub
} + _{
Us
} + _{
Um
} is quadratic in the variables w, _{
βx
} and _{
βy
} (^{Brush and Almroth, 1975}), where _{
Ub
} , _{
Us
} and _{
Um
} stand for the change due to buckling of the bending strain energy, transverse shear strain energy and membrane strain energy, respectively. The DKT formulation neglects Us contribution and enforces the Kirchhoff constraints at selected points of the element.

2.1 Bending Stiffness Matrix

After enforcement of the Kirchhoff constraints at selected points of the DKT element shown in Figure 1, ^{Batoz et al. (1980}) demonstrate that the normal rotations are effectively approximated by

with the nodal values of *w*, _{
βx = θy
} and _{
βy = θx
} collected in

The shape functions identified by

Have

given in terms of the nondimensional coordinates

The triangle has area *A* and vertices with global coordinates (_{
xi
} , _{
yi
} ). The parameters

are defined as functions of

with *k* = 4, 5, 6 for the triangle sides _{
ij
} = 23, 31, 12, respectively.

Making use of (1), the bending strain energy

reads

where _{
Dij
} are the plate bending rigidities. The element bending stiffness matrix **k** is exactly obtained by ^{Batoz et al. (1980}) using three-point Gauss quadrature and subsequently stated by ^{Batoz (1982}) in explicit form in local coordinates as

where

The matrix **α** is given in Table I of ^{Batoz (1982}).

2.2 Geometric Stiffness Matrix

If the plate is assumed to deform as a membrane prior to buckling with in-plane force distribution _{
Nx
} , _{
Ny
} and _{
Nxy
} , the change in membrane strain energy due to buckling reads (^{Reddy, 2004})

The Kirchhoff constraint is now invoked through the plate so that the rotations _{
βx
} and _{
βy
} can be related to the transverse displacement derivatives by means of

Under this assumption, expression (12) is replaced by

Substitution of (1) into (14) yields

where

is the consistent geometric stiffness matrix and

In order to develop an analytical expression for (16) in closed form, the crucial step is the decomposition of **H** into the product

where

Independent of *ξη* coordinates, **α**
*g* is given by

In view of (18),

where

If **N** is constant over the element,

with

Note that **k**
*g* refers to the global system *xy* and has been exactly obtained in closed form.

The widely used inconsistent geometric stiffness matrix (^{Doyle, 2001})

in which

will be also employed in the buckling analysis for comparison purposes. The matrix is established by adopting the linear displacement field

in order to compute the derivatives ∂w / *∂x* and ∂w / *∂y* required in (12).

The discrete linearized buckling problem, given by *δ*∏= _{
δUb
} + _{
δUm
} = 0 , reduces to the linear eigenvalue problem

Matrix **K** is constructed by assembly of element matrices **k**, whereas **k**
**
g
** is constructed by assembly of element matrices

**k**

*g*(or if the inconsistent geometric stiffness matrix is adopted). Vector

**D**collects the nodal displacements measured from the buckling point and λ is a constant by which the applied load must be multiplied to cause buckling.

3 NUMERICAL TESTS

All the examples refer to square plates with side *a*, as illustrated in Figure 2, subjected to one of the following uniform in-plane loadings: uniaxial compression _{
px
} (_{
py =
}
_{
pxy =
} 0), biaxial compression _{
px = py
} (_{
pxy =
} 0), biaxial compression _{
px
} and tension_{
py = px
} (_{
pxy =
} 0), shear _{
pxy
} (_{
px = py
} = 0. The adopted materials are defined as being

The boundary conditions will be denoted by SSSS (all edges are simply supported), SSCC (edges are simply supported at *x* = 0, *a* and clamped at *y* = 0 , *a*) and CCCC (all edges are clamped). Assuming the in-plane deformation to be unconstrained prior to buckling, **N** has components N_{
x
} = - _{
px
} , N_{
y
} = - _{
py
} and N_{
xy
} = -_{
pxy
} .

Due to symmetry, only one quarter of the plates is defined as the solution domain for uniaxial compression, biaxial compression, biaxial compression and tension. For shear loading, however, the buckled shape does not exhibit any symmetry and thus the whole plate has to be discretized. A typical 4×4 mesh in a quarter plate is depicted in Figure 2.

Tables 1-3 summarize the results for the nondimensional buckling loads

which are obtained for different meshes *N* × *N* (*N* is the number of elements along the plate half side *a*/2 or the plate side *a*, depending on whether the symmetry is adopted or not). Departing from a very coarse mesh 1 × 1, successive meshes are generated by splitting all elements in the current mesh into four smaller ones. For comparison purposes, either exact or approximate reference values of (30) based on Kirchhoff plate theory are provided. While exact values stem from Navier or Lévy solutions (^{Reddy, 2004}), approximate values are accurate results obtained from Ritz solutions employing hierarchic set of functions built from modified Legendre polynomials (^{Bardell, 1991}). No results are tabled for shear loading with 1 × 1 mesh because the introduction of the boundary conditions after discretization of the whole plate does not leave any degree of freedom available. Moreover, 1 × 1 mesh in a quarter may lead to singular geometric stiffness matrices and, consequently, to infinite buckling loads. Both consistent and inconsistent formulations always provide convergence toward exact results with mesh refinement.

^{a} discretization of the whole plate ^{b} consistent matrix ^{c} inconsistent matrix ^{d} exact value ^{e} Ritz solution

^{a} discretization of the whole plate ^{b} consistent matrix ^{c} inconsistent matrix ^{d} exact value ^{e} Ritz solution

The graphical form of the isotropic plate results are shown in Figures 3-6, where more refined meshes (32×32 and 64×64) are included in the log-log plots of Figures 5 and 6 to ensure that the results obtained with the consistent formulation for biaxial compression/tension and shear loads have also achieved the asymptotic range. Defining *h* = *a*/2*N* (or *h* = *a*/*N* if symmetry is not adopted) as a measure of the element size, and _{
ε =|kx/
}*-* 1*|* for _{
pxy =
} 0 or _{
ε = |kxy/
} -1 | for_{
pxy
} ≠ 0 as a measure of the error, the figures also show *ε* plotted against *h* on a log-log scale whose slope of the curves is proportional to the rate of convergence.

For the analyzed isotropic plates, the inconsistent formulation provides an approximately quadratic rate of convergence through the mesh refinement, except for SSCC plates under biaxial compression/tension for which the rate becomes nearly cubic. On the other hand, the rate provided by the consistent formulation may vary with mesh refinement depending on the plate boundary conditions and loading. Such a rate is often higher for coarse meshes, achieving lower values for fine meshes. In the asymptotic range, the consistent formulation also exhibits an approximately quadratic rate of convergence, except for SSSS plates under biaxial compression/tension for which the rate becomes nearly quartic. Although both formulations always provide convergence toward exact results, the buckling load predictions with the consistent formulation are generally better, with higher convergence rate, for coarse meshes and complex modes (plate subjected to shear loading).

The consistent formulation may be a little more time consuming because expression (21) involves more arithmetic operations than (25). In order to confirm our expectation, Figure 7 plots the CPU time to solve CCCC isotropic plates under shear loading. All the calculations were performed in a 3 GH Pentium personal computer with 4 GB of RAM, running a 64-bit Windows 7. It has been observed that the time is practically independent of the loading type and that the change of the boundary conditions from CCCC to SSCC or SSSS has little effect on the time change for the examples considered. Both formulations approximately allows computing the solution in a logarithm of time proportional to *N*.

^{a} discretization of the whole plate ^{b} consistent matrix ^{c} inconsistent matrix ^{d} Ritz solution

4 CONCLUSIONS

To the simple and effective DKT element is offered an explicit form of its consistent geometric stiffness matrix in global coordinates. This proposed form avoids the expensive exact integration by seven-point Gauss quadrature and transformations from local to global coordinates. Based on the examples considered, the following conclusions can be drawn with respect to the performances of the consistent and inconsistent formulations:

Both formulations always provide convergence toward exact results with mesh refinement.

A solution obtained by the consistent formulation is a little more time consuming, and both formulations approximately allows computing the solution in a logarithm of CPU time proportional to

*N*.An approximately quadratic rate of convergence is usually observed for the consistent formulation in the asymptotic range, whereas the inconsistent formulation has such a rate through all the mesh refinement.

For complex buckling modes, as it is the case of plates under shear loading, the consistent formulation gives more precise results for coarse meshes.