SciELO - Scientific Electronic Library Online

 
vol.14 issue4Dynamic Response of Dam-Reservoir Systems: Review and a Semi-Analytical ProposalInvestigation of High-Speed Cryogenic Machining Based on Finite Element Approach author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

Share


Latin American Journal of Solids and Structures

Print version ISSN 1679-7817On-line version ISSN 1679-7825

Lat. Am. j. solids struct. vol.14 no.4 Rio de Janeiro Mar. 2017

http://dx.doi.org/10.1590/1679-78252799 

Articles

An Explicit Consistent Geometric Stiffness Matrix for the DKT Element

Eliseu Lucena Netoa  * 

Marcus Antônio Ferreira Araripea 

Francisco Alex Correia Monteiroa 

José Antônio Hernandesa 

a Instituto Tecnológico de Aeronáutica, 12228-900 São José dos Campos, SP, Brazil. 1 mafararipe@gmail.com, 2 facm@ita.br, 3 hernandes@ita.br

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

(1)

with the nodal values of w, βx = θy and βy = θx collected in

(2)

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

The shape functions identified by

(3)

Have

(4)

given in terms of the nondimensional coordinates

(5)

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

(6)

are defined as functions of

(7)

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

Making use of (1), the bending strain energy

(8)

reads

(9)

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

(10)

where

(11)

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)

(12)

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

(13)

Under this assumption, expression (12) is replaced by

(14)

Substitution of (1) into (14) yields

(15)

where

(16)

is the consistent geometric stiffness matrix and

(17)

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

(18)

where

(19)

Independent of ξη coordinates, α g is given by

(20)

In view of (18),

(21)

where

(22)

If N is constant over the element,

(23)

with

(24)

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)

(25)

in which

(26)

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

(27)

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

(28)

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

(29)

Figure 2: Square plate subjected to uniform in-plane 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 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

(30)

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.

Table 1: Nondimensional buckling load for SSSS square plates: kx = px a 2 / π2 D 22 ; kxy = pxy a 2 / π2 D 22 and is the reference value. 

kx / ( pxy = 0) kxy / a
mesh py = 0 py = px py = - px px = py = 0
isotropic 1×1 1.119b 1.121 1.231 -
0.903c 0.903 -
2×2 1.016 1.020 1.003 1.546
0.992 0.992 1.167 0.422
4×4 1.003 1.004 0.994 1.007
0.998 0.998 1.041 0.813
8×8 1.001 1.001 0.997 0.982
0.999 0.999 1.010 0.954
16×16 1.000 1.000 0.999 0.994
1.000 1.000 1.003 0.988
4.000d 2.000d 8.333d 9.325e
orthotropic 1×1 1.085 1.097 1.264 -
0.837 0.837 -
2×2 1.016 1.012 0.996 1.649
0.973 0.969 1.245 0.449
4×4 1.003 1.001 0.989 1.029
0.991 0.991 1.056 0.826
8×8 1.001 1.000 0.995 0.986
0.998 0.997 1.013 0.950
16×16 1.000 1.000 0.999 0.999
0.999 0.999 1.003 0.986
13.000d 6.500d -11.333d 28.066e

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

Table 2: Nondimensional buckling load for SSCC square plates: kx = px a 2 / π2 D 22 ; kxy = pxy a 2 / π2 D 22 and is the reference value. 

kx / ( pxy = 0) kxy / a
mesh py = 0 py = px py = - px px = py = 0
isotropic 1×1 1.092b 1.140 0.951 -
c 0.743 -
2×2 1.004 1.055 0.913 1.385
1.054 1.032 1.071 0.487
4×4 0.995 1.014 0.958 0.966
1.014 1.013 1.010 0.854
8×8 0.998 1.003 0.988 0.976
1.003 1.003 1.001 0.964
16×16 1.000 1.001 0.997 0.993
1.001 1.001 1.000 0.991
7.691d 3.830d 10.788d 12.565e
orthotropic 1×1 1.338 1.416 1.072 -
0.739 0.846 -
2×2 1.182 1.104 1.395 1.858
1.089 1.077 2.120 0.497
4×4 1.038 1.020 1.044 1.114
1.020 1.023 1.179 0.864
8×8 1.008 1.004 1.008 1.016
1.004 1.004 1.039 0.958
16×16 1.002 1.001 1.002 1.000
1.001 1.001 1.009 0.988
17.604d 7.693d 15.169d 29.567e

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

Figure 3: Uniaxial buckling load ( px ≠ 0, py = pxy = 0) for isotropic plates: nondimensional load kx / versus N and error ε versus element size h

Figure 4: Biaxial compression buckling load ( px = py , pxy = 0) for isotropic plates: nondimensional load kx / versus N and error ε versus element size h

Figure 5: Biaxial compression/tension buckling load (( px = - py , pxy = 0)) for isotropic plates: nondimensional load kx / versus N and error ε versus element size h

Figure 6: Shear buckling load ( pxy ≠ 0, px = py = 0) for isotropic plates: nondimensional load kxy / versus N and error ε versus element size h

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.

Figure 7: CPU time consumption in seconds versus N to solve CCCC isotropic plates under shear loading. 

Table 3: Nondimensional buckling load for CCCC square plates: kx = px a 2 / π2 D 22 ; kxy = pxy a 2 / π2 D 22 and is the reference value. 

kx / ( pxy = 0) kxy / a
mesh py = 0 py = px py = - px px = py = 0
isotropic 1×1 2.002b 1.920 1.087 -
0.809c 0.768 -
2×2 1.055 1.099 1.096 1.326
1.043 1.062 1.340 0.583
4×4 1.010 1.019 0.989 1.000
1.018 1.019 1.080 0.861
8×8 1.002 1.004 0.994 0.980
1.005 1.005 1.018 0.962
16×16 1.000 1.001 0.998 0.993
1.001 1.001 1.004 0.989
10.074d 5.304d 14.966d 14.642d
orthotropic 1×1 1.622 0.975 -
0.649 -
2×2 1.094 1.247 1.160 1.257
0.992 1.433 1.237 0.768
4×4 1.016 1.016 1.015 1.137
1.002 1.060 1.191 0.909
8×8 1.004 1.002 0.999 1.014
1.000 1.012 1.042 0.947
16×16 1.001 1.000 0.999 0.994
1.000 1.003 1.010 0.983
46.289d 16.681d 24.066d 46.726d

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.

References

ABAQUS Theory Manual, Version 6.12, ABAQUS, Inc., Dassault Systèmes Simulia Corp., Providence, USA, 2012. [ Links ]

Bardell, N.S. (1991). Free vibration analysis of a flat plate using the hierarchical finite element method. Journal of Sound and Vibration 151(2): 263-289. [ Links ]

Bathe, K.J., Ho, L.W. (1981). A simple and effective element for analysis of general shell structures. Computers and Structures 13(5-6): 673-681. [ Links ]

Batoz, J.L. (1982). An explicit formulation for an efficient triangular plate-bending element. International Journal for Numerical Methods in Engineering 18(7): 1077-1089. [ Links ]

Batoz, J.L., Bathe, K.J., Ho, L.W. (1980). A study of three-node triangular plate bending elements. International Journal for Numerical Methods in Engineering 15(12): 1171-1812. [ Links ]

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): 225-251. [ Links ]

Bernadou, M. (1996). Finite Element Methods for Thin Shell Problems. John Wiley, New York. [ Links ]

Brush, D.O., Almroth, B.O. (1975). Buckling of Bars, Plates, and Shells. McGraw-Hill, New York. [ Links ]

Cook, R.D. (1993). Further development of a three-node triangular shell element. International Journal for Numerical Methods in Engineering 36(8): 1413-1425. [ Links ]

Cook, R.D., Malkus, D.S., Plesha, M.E., Witt R.J. (2002). Concepts and Applications of Finite Element Analysis. John Wiley, New York. [ Links ]

Doyle, J.F. (2001). Nonlinear Analysis of Thin-Walled Structures: Statics, Dynamics, and Stability. Springer, New York. [ Links ]

Garnet, H., Pifko, A.B. (1983). An efficient triangular plate bending finite element for crash simulation. Computers and Structures 16(1-4): 371-379. [ Links ]

Jeyachandrabose, C., Kirkhope, J., Babu, C. R. (1985). An alternative explicit formulation for the DKT plate-bending element. International Journal for Numerical Methods in Engineering 21(7): 1289-1293. [ Links ]

Khosravi, P., Ganesan, R., Sedaghati, R. (2007). Corotational non-linear analysis of thin plates and shells using a new shell element. International Journal for Numerical Methods in Engineering 69(4): 859-885. [ Links ]

Kikuchi, F. (1975). On a finite element scheme based on the discrete Kirchhoff assumption. Numerische Mathematik 24(3): 211-231. [ Links ]

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(1-4): 461-472. [ Links ]

Reddy, J.N. (2004). Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. CRC Press, Boca Raton. [ Links ]

Strickling, J.A., Haisler, W.E., Tisdale, P.R., Gunderson, R. (1969). A rapidly converging triangular plate element. AIAA Journal 7(1): 180-181. [ Links ]

Wu, S., Li, G., Belytschko, T. (2005). A DKT shell element for dynamic large deformation analysis. Communications in Numerical Methods in Engineering 21(11): 651-674. [ Links ]

Zhang, S., Zhang, Y., Zhigao, H., Huamin, Z., Li, J. (2015). The inter-element coupling effect of triangular flat shells. Engineering Computations 32(7): 1959-1980. [ Links ]

Received: January 16, 2016; Revised: December 19, 2016; Accepted: January 29, 2017

* Corresponding author: eliseu@ita.br

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License