An Explicit Consistent Geometric Stiffness Matrix for the DKT Element

A large number of references dealing with the geometric stiffness matrix of the DKT finite element exist in the literature, where nearly all of them adopt an inconsistent form. While such a matrix may be part of the element to treat nonlinear problems in general, it is of crucial importance for linearized buckling analysis. The present work seems to be the first to obtain an explicit expression for this matrix in a consistent way. Numerical results on linear buckling of plates assess the element performance either with the proposed explicit consistent matrix, or with the most commonly used inconsistent matrix.

Latin American Journal of Solids and Structures 14 (2017) 613-628 lated from the three vertex values of , and .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 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 within the element (Bathe and Ho, 1981;Mateus et al., 1997;Doyle, 2001;Khosravi et al., 2007).If the modified approximations obtained for and 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 within the element.

LINEARIZED BUCKLING FORMULATION
In the linearized buckling analysis of a Reissner-Mindlin plate, the potential energy Π is quadratic in the variables , and (Brush and Almroth, 1975), where , and 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 contribution and enforces the Kirchhoff constraints at selected points of the element.

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 (1) with the nodal values of , and collected in . (2) The shape functions identified by Have given in terms of the nondimensional coordinates Latin American Journal of Solids and Structures 14 (2017) 613-628 1 2 1 2 . (5) The triangle has area and vertices with global coordinates , .The parameters are defined as functions of (7) with 4, 5, 6 for the triangle sides 23, 31, 12, respectively.Making use of (1), the bending strain energy where are the plate bending rigidities.The element bending stiffness matrix 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).

Geometric Stiffness Matrix
If the plate is assumed to deform as a membrane prior to buckling with in-plane force distribution , and , 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 and 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 into the product ( 18) where 1 0 0 0 0 0 0 . ( Independent of coordinates, is given by In view of (18), 2 Note that refers to the global system 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 / and / required in (12).The discrete linearized buckling problem, given by Π 0, reduces to the linear eigenvalue problem Matrix is constructed by assembly of element matrices , whereas is constructed by assembly of element matrices (or ̅ if the inconsistent geometric stiffness matrix is adopted).Vector collects the nodal displacements measured from the buckling point and λ is a constant by which the applied load must be multiplied to cause buckling.

NUMERICAL TESTS
All the examples refer to square plates with side , as illustrated in Figure 2, subjected to one of the following uniform in-plane loadings: uniaxial compression ( 0), biaxial compression ( 0), biaxial compression and tension ( 0), shear ( 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 0, and clamped at 0, ) and CCCC (all edges are clamped).Assuming the in-plane deformation to be unconstrained prior to buckling, has components , and .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 ( is the number of elements along the plate half side /2 or the plate side , 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.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 /2 (or / if symmetry is not adopted) as a measure of the element size, and / 1 for 0 or / 1 for 0 as a measure of the error, the figures also show plotted against 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 quadrat-Latin American Journal of Solids and Structures 14 (2017) 613-628 ic 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 .

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 . 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.

Figure 1 :
Figure 1: DKT element in the global and nondimensional coordinates.

Figure 2 :
Figure 2: Square plate subjected to uniform in-plane loads , , and 4 4 mesh in a quarter plate.

Table 2 :
Nondimensional buckling load for SSCC square plates:

Table 3 :
Nondimensional buckling load for CCCC square plates: