Abstract
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.
Keywords:
DKT; geometric stiffness matrix; discrete Kirchhoff triangle; buckling
1 INTRODUCTION
It is surprisingly difficult to design a simple threenode 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., 2002Cook, R.D., Malkus, D.S., Plesha, M.E., Witt R.J. (2002). Concepts and Applications of Finite Element Analysis. John Wiley, New York.). The DKT (discrete Kirchhoff triangle) element was first published by Stricklin et al. (1969Strickling, J.A., Haisler, W.E., Tisdale, P.R., Gunderson, R. (1969). A rapidly converging triangular plate element. AIAA Journal 7(1): 180181.) 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., 1980Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of threenode triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12): 11711812.). 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 ReissnerMindlin 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. (1980Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of threenode triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12): 11711812.), which concluded that the element is one of the most efficient and reliable threenode triangular elements for Kirchhoff plates. Batoz (1982Batoz, J.L. (1982). An explicit formulation for an efficient triangular platebending element. International Journal for Numerical Methods in Engineering 18(7): 10771089.) also obtained an explicit expression for the element bending stiffness matrix in local coordinates, which was subsequently coded in global coordinates (Jeyachandrabose et al., 1985Jeyachandrabose, C., Kirkhope, J., Babu, C. R. (1985). An alternative explicit formulation for the DKT platebending element. International Journal for Numerical Methods in Engineering 21(7): 12891293.). A theoretical study of the element convergence rate can be found in Kikuchi (1975Kikuchi, F. (1975). On a finite element scheme based on the discrete Kirchhoff assumption. Numerische Mathematik 24(3): 211231.) and Bernadou (1996Bernadou, M. (1996). Finite Element Methods for Thin Shell Problems. John Wiley, New York.).
Since the insightful work of Batoz et al. (1980Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of threenode triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12): 11711812.), 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, 1981Bathe, K.J., Ho, L.W. (1981). A simple and effective element for analysis of general shell structures. Computers and Structures 13(56): 673681.; Cook, 1993Cook, R.D. (1993). Further development of a threenode triangular shell element. International Journal for Numerical Methods in Engineering 36(8): 14131425.; Khosravi et al., 2007Khosravi, P., Ganesan, R., Sedaghati, R. (2007). Corotational nonlinear analysis of thin plates and shells using a new shell element. International Journal for Numerical Methods in Engineering 69(4): 859885.; Zhang et al., 2015Zhang, S., Zhang, Y., Zhigao, H., Huamin, Z., Li, J. (2015). The interelement coupling effect of triangular flat shells. Engineering Computations 32(7): 19591980.). Highspeed component impact analysis (Wu et al., 2005Wu, S., Li, G., Belytschko, T. (2005). A DKT shell element for dynamic large deformation analysis. Communications in Numerical Methods in Engineering 21(11): 651674.) demonstrates, for instance, that this element performs as good as the main element (Belytschko et al., 1984Belytschko, T., Lin, J.L., Tsay, C.S. (1984). Explicit algorithms for the nonlinear dynamics of shells. Computer Methods in Applied Mechanics and Engineering 42(2): 225251.) 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, 2012ABAQUS Theory Manual, Version 6.12, ABAQUS, Inc., Dassault Systèmes Simulia Corp., Providence, USA, 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, 1981Bathe, K.J., Ho, L.W. (1981). A simple and effective element for analysis of general shell structures. Computers and Structures 13(56): 673681.; Mateus et al., 1997Mateus, H.C., Soares, C.M.M., Soares, C.A.M. (1997). Buckling sensitivity analysis and optimal design of thin laminated structures. Computers and Structures 64(14): 461472.; Doyle, 2001Doyle, J.F. (2001). Nonlinear Analysis of ThinWalled Structures: Statics, Dynamics, and Stability. Springer, New York.; Khosravi et al., 2007Khosravi, P., Ganesan, R., Sedaghati, R. (2007). Corotational nonlinear analysis of thin plates and shells using a new shell element. International Journal for Numerical Methods in Engineering 69(4): 859885.). 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, 1983Garnet, H., Pifko, A.B. (1983). An efficient triangular plate bending finite element for crash simulation. Computers and Structures 16(14): 371379.), 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 sevenpoint 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 ReissnerMindlin plate, the potential energy ∏ = _{Ub} + _{Us} + _{Um} is quadratic in the variables w, _{βx} and _{βy} (Brush and Almroth, 1975Brush, D.O., Almroth, B.O. (1975). Buckling of Bars, Plates, and Shells. McGrawHill, New York.), 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. (1980Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of threenode triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12): 11711812.) 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. (1980Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of threenode triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12): 11711812.) using threepoint Gauss quadrature and subsequently stated by Batoz (1982Batoz, J.L. (1982). An explicit formulation for an efficient triangular platebending element. International Journal for Numerical Methods in Engineering 18(7): 10771089.) in explicit form in local coordinates as
where
The matrix α is given in Table I of Batoz (1982Batoz, J.L. (1982). An explicit formulation for an efficient triangular platebending element. International Journal for Numerical Methods in Engineering 18(7): 10771089.).
2.2 Geometric Stiffness Matrix
If the plate is assumed to deform as a membrane prior to buckling with inplane force distribution _{Nx} , _{Ny} and _{Nxy} , the change in membrane strain energy due to buckling reads (Reddy, 2004Reddy, J.N. (2004). Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. CRC Press, Boca Raton.)
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, 2001Doyle, J.F. (2001). Nonlinear Analysis of ThinWalled Structures: Statics, Dynamics, and Stability. Springer, New York.)
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 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. if the inconsistent geometric stiffness matrix is adopted). Vector
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 inplane 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
Square plate subjected to uniform inplane loads _{px} , _{py} , _{pxy} and 4x4 mesh in a quarter plate.
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 inplane 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 13 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 Reddy, 2004Reddy, J.N. (2004). Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. CRC Press, Boca Raton.), approximate values are accurate results obtained from Ritz solutions employing hierarchic set of functions built from modified Legendre polynomials (Bardell, 1991Bardell, N.S. (1991). Free vibration analysis of a flat plate using the hierarchical finite element method. Journal of Sound and Vibration 151(2): 263289.). 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. of (30) based on Kirchhoff plate theory are provided. While exact values stem from Navier or Lévy solutions (
The graphical form of the isotropic plate results are shown in Figures 36, where more refined meshes (32×32 and 64×64) are included in the loglog 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/2N (or h = a/N if symmetry is not adopted) as a measure of the element size, and _{ε =kx/} 1 for _{pxy =} 0 or _{ε = kxy/}_{pxy} ≠ 0 as a measure of the error, the figures also show ε plotted against h on a loglog scale whose slope of the curves is proportional to the rate of convergence. 1  for
Uniaxial buckling load (_{px} ≠ 0,_{py = pxy =} 0) for isotropic plates: nondimensional load _{kx} /N and error ε versus element size h. versus
Biaxial compression buckling load (_{px} = _{py} , _{pxy} = 0) for isotropic plates: nondimensional load _{kx} /N and error ε versus element size h. versus
Biaxial compression/tension buckling load ((_{px} = _{py} , _{pxy} = 0)) for isotropic plates: nondimensional load _{kx} /N and error ε versus element size h. versus
Shear buckling load (_{pxy} ≠ 0,_{px = py} = 0) for isotropic plates: nondimensional load _{kxy /}N and error ε versus element size h. versus
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 64bit 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.
CPU time consumption in seconds versus N to solve CCCC isotropic plates under shear loading.
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 sevenpoint 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.
References
 ABAQUS Theory Manual, Version 6.12, ABAQUS, Inc., Dassault Systèmes Simulia Corp., Providence, USA, 2012.
 Bardell, N.S. (1991). Free vibration analysis of a flat plate using the hierarchical finite element method. Journal of Sound and Vibration 151(2): 263289.
 Bathe, K.J., Ho, L.W. (1981). A simple and effective element for analysis of general shell structures. Computers and Structures 13(56): 673681.
 Batoz, J.L. (1982). An explicit formulation for an efficient triangular platebending element. International Journal for Numerical Methods in Engineering 18(7): 10771089.
 Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of threenode triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12): 11711812.
 Belytschko, T., Lin, J.L., Tsay, C.S. (1984). Explicit algorithms for the nonlinear dynamics of shells. Computer Methods in Applied Mechanics and Engineering 42(2): 225251.
 Bernadou, M. (1996). Finite Element Methods for Thin Shell Problems. John Wiley, New York.
 Brush, D.O., Almroth, B.O. (1975). Buckling of Bars, Plates, and Shells. McGrawHill, New York.
 Cook, R.D. (1993). Further development of a threenode triangular shell element. International Journal for Numerical Methods in Engineering 36(8): 14131425.
 Cook, R.D., Malkus, D.S., Plesha, M.E., Witt R.J. (2002). Concepts and Applications of Finite Element Analysis. John Wiley, New York.
 Doyle, J.F. (2001). Nonlinear Analysis of ThinWalled Structures: Statics, Dynamics, and Stability. Springer, New York.
 Garnet, H., Pifko, A.B. (1983). An efficient triangular plate bending finite element for crash simulation. Computers and Structures 16(14): 371379.
 Jeyachandrabose, C., Kirkhope, J., Babu, C. R. (1985). An alternative explicit formulation for the DKT platebending element. International Journal for Numerical Methods in Engineering 21(7): 12891293.
 Khosravi, P., Ganesan, R., Sedaghati, R. (2007). Corotational nonlinear analysis of thin plates and shells using a new shell element. International Journal for Numerical Methods in Engineering 69(4): 859885.
 Kikuchi, F. (1975). On a finite element scheme based on the discrete Kirchhoff assumption. Numerische Mathematik 24(3): 211231.
 Mateus, H.C., Soares, C.M.M., Soares, C.A.M. (1997). Buckling sensitivity analysis and optimal design of thin laminated structures. Computers and Structures 64(14): 461472.
 Reddy, J.N. (2004). Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. CRC Press, Boca Raton.
 Strickling, J.A., Haisler, W.E., Tisdale, P.R., Gunderson, R. (1969). A rapidly converging triangular plate element. AIAA Journal 7(1): 180181.
 Wu, S., Li, G., Belytschko, T. (2005). A DKT shell element for dynamic large deformation analysis. Communications in Numerical Methods in Engineering 21(11): 651674.
 Zhang, S., Zhang, Y., Zhigao, H., Huamin, Z., Li, J. (2015). The interelement coupling effect of triangular flat shells. Engineering Computations 32(7): 19591980.
Publication Dates

Publication in this collection
Mar 2017
History

Received
16 Jan 2016 
Reviewed
19 Dec 2016 
Accepted
29 Jan 2017