New Formula for Geometric Stiffness Matrix Calculation

I. Němec^1*, M. Trcala², I. Ševčík², H. Štekbauer¹

¹Faculty of Civil Engineering, Brno University of Technology, Brno, Czech Republic.
²FEM Consulting, S.R.O., Brno, Czech Republic.

DOI: 10.4236/jamp.2016.44084                

Abstract

The standard formula for geometric stiffness matrix calculation, which is convenient for most engineering applications, is seen to be unsatisfactory for large strains because of poor accuracy, low convergence rate, and stability. For very large compressions, the tangent stiffness in the direction of the compression can even become negative, which can be regarded as physical nonsense. So in many cases rubber materials exposed to great compression cannot be analyzed, or the analysis could lead to very poor convergence. Problems with the standard geometric stiffness matrix can even occur with a small strain in the case of plastic yielding, which eventuates even greater practical problems. The authors demonstrate that amore precisional approach would not lead to such strange and theoretically unjustified results. An improved formula that would eliminate the disadvantages mentioned above and leads to higher convergence rate and more robust computations is suggested in this paper. The new formula can be derived from the principle of virtual work using a modified Green-Lagrange strain tensor, or from equilibrium conditions where in the choice of a specific strain measure is not needed for the geometric stiffness derivation (which can also be used for derivation of geometric stiffness of a rigid truss member). The new formula has been verified in practice with many calculations and implemented in the RFEM and SCIA Engineer programs. The advantages of the new formula in comparison with the standard formula are shown using several examples.

Keywords

Geometric Stiffness, Stress Stiffness, Initial Stress Stiffness, Tangent Stiffness Matrix, Finite Element Method, Principle of Virtual Work, Strain Measure

Received 4 June 2015; accepted 24 April 2016; published 27 April 2016

1. Introduction

Stress stiffening is an important source of stiffness and must be taken into account when analyzing structures. The standard formula for geometric stiffness matrices is introduced by a number of authors, such as Zienkiewicz, Bathe, Cook, Belytschko, Simo, Hughes, Bonet, de Souza Neto and others [1] - [10] . The standard formula has been shown to be satisfactory in a large amount of cases, though certain difficulties such as low accuracy, poor convergence rate and poor solution stability were discovered when solving problems that included the evaluation of extreme stress and strain states. Some authors, e.g. Cook [4] , have suggested an improvement for bars and some authors dealt with nonlinear models describing large (finite) deformation (strain) behavior of materials and structures [11] - [21] . However, as far as the authors know, no general solution to the problem has been suggested for a 2D or 3D continuum. Upon this ascertainment, thoughts arose concerning the physical essence of geometric (or stress) stiffness and the formula for evaluating geometric stiffness matrices. As a result, a new formula for geometric stiffness matrix calculation is suggested. The presentation of this new formula, which should substantially improve analysis of structures exposed to large strain, is the subject matter of this paper. In Section 2, the standard formula for geometric stiffness matrices is presented. Section 3 shows the physical background of geometric stiffness based on equilibrium. In Section 4, the new, improved formula for geometric matrices is introduced. The advantages of the new formula, including a substantially improved rate of convergence and stability, are demonstrated by examples in Section5. Conclusions are presented in Section 6.

2. The Standard Formula for Stress-Stiffness Matrices

Let us show the general calculation algorithm for the geometric stiffness matrix (sometimes also called the stress stiffness matrix or initial stress matrix) of an element in an updated Lagrangian formulation.

Let the following hold for each component of displacement vector:

(1)

where is the value of displacement in node a and n is the number of element nodes.

Let us define matrix as follows:

(2)

where is the unit diagonal matrix of the order 3 × 3, where 3 is the dimension of the problem. Then, the following relation can be written for the displacement vector:

(3)

where is the vector of deformation parameters of the element containing all the components in such an arrangement that for each node a all components are listed.

Let us define matrix containing the first derivatives of base functions for node a with respect to spatial coordinates

(4)

and matrix, which is formed by sub-matrices

(5)

The operator denotes the tensor (Kronecker) matrix product.

Further, let us define matrix by multiplying each component of the Cauchy stress tensor by the unit diagonal matrix:

(6)

If state of the stress is not negligible, the potential energy of the internal forces should be completed by the following term:

(7)

Then, the following formula for the geometric matrix of the element can be written:

(8)

Integration is carried out on the deformed body (in the current configuration) and the derivatives are performed with respect to the spatial coordinates.

The component of the matrix relating the element node a to the element node b can also be written simply in matrix notation:

(9)

or in indicial notation:

(10)

Similar formulae also hold for a total Lagrangian formulation, but the second Piola-Kirchhoff stress tensor is then used instead of the Cauchy stress, and integration is carried out on the undeformed body (in the original configuration) while the derivatives are performed with respect to the material coordinates.

3. The Source of Geometric Stiffness―The Physical Background

Let us consider the truss member shown in Figure 1. Node 2 is loaded by the force F parallel to the x axis and sliding in the same direction. The equilibrium equation in the x direction in node 2 can be written as follows

(11)

where is the horizontal component of the internal force at node 2 and. is the residual or out-of-balance force. The horizontal stiffness at node 2 is defined simply by the relation

(12)

This formula is independent of any strain measure or pertinent constitutive relations. It can be seen that stiffness consists of two parts. The first part, , represents the material stiffness and depends on the strain measure and constitutive relations. The second part, , which does not depend on the material or the strain and stress measures chosen, but only on the geometry and the normal force, represents so-called geometric stiffness. It can be seen that if the angle is zero, no geometric stiffness will occur regardless of the normal force value.

Let us show a derivation of a formula for geometric stiffness matrix of a truss member (see Figure 2) in a finite element formulation and let us start with a simple derivation based on equilibrium conditions.

Let be a vector of the nodal displacements of an element, and let be a vector of residual forces; the stiffness matrix of the element can then be defined as follows:

(13)

Figure 1. Truss member in an arbitrary position in 2D.

Figure 2. Truss member: the x axis is the axis of the rod in its original position.

A geometric (stress) stiffness matrix can be obtained by an equilibrium condition when only the initial stress state and pertinent infinitesimal nodal displacement for each row of the matrix is taken into account. Such a definition of a geometric stiffness matrix is independent of the strain tensor chosen.

To simplify the following derivations let’s introduce both, the coordinates with the axis aligned with the axis of the rod and corresponding displacement vector and let’s restrict the deformation to the plane.

Let the vector of the nodal displacements of the element be

(14)

where and are the displacement components in the direction of the and axis, respectively. The well known material stiffness matrix of the truss element in 2D is then defined by the following relation:

(15)

Note that the truss element has no lateral material stiffness.

In general, arbitrary term of a stiffness matrix is defined as the derivative of an unbalanced force with respect to the deformation parameter as is defined by (13). Based on this definition, the geometric stiffness matrix of the truss element subjected to tensile force N can be easily derived. The moment equilibrium condition for the truss member in the configuration with the lateral displacement in node 1 is sufficient to obtain the transversal diagonal stiffness term:

The moment equilibrium condition can be written as follows:

(16)

For the infinitesimal angle it can be assumed that, and the following term for the stiffness term can be derived:

(17)

When introducing a displacement in the direction of the axis of the member, the end forces are in equilibrium and no additional force and therefore no geometric stiffness will occur in this direction.

From equilibrium equations and symmetry of the stiffness matrix it is easy to determine the other coefficients of the geometric stiffness matrix, particularly, and. The remaining coefficients of the matrix are zeros. The geometric stiffness matrix then has the following form:

(18)

The same formula corresponds with Formula (12) and is presented also by Cook in [4] , the same as many other authors. The geometric stiffness matrix for a truss member can also be derived from the principle of virtual work, which will be described later. Then a strain measure and constitutive law must be introduced, which is not applicable for a rigid truss, where geometric stiffness also exists.

The resulting tangent stiffness matrix is defined as the sum of the material and geometric stiffness matrix:

(19)

When applying the general standard algorithm for geometric stiffness matrices to the truss element in question, we obtain:

(20)

(21)

where is the identity matrix of order 2 and the base functions are defined as follows:

(22)

(23)

where

(24)

Substituting in the formulae

(25)

(26)

the formula for the geometric stiffness matrix reads:

(27)

This geometric stiffness matrix differs from that in Formula (18) and introduces also an axial stiffening. But no reason was found by the authors for concluding that normal force had led to a change in the axial stiffness of the element. So let us derive the geometric stiffness matrix of a truss element in a more undisputable way based on the principle of virtual work.

With deformation restricted to the plane, the Green-Lagrange strain tensor is defined

(28)

where

(29)

For truss the principle of virtual work becomes

(30)

where is the 2nd Piola-Kirchhoff stress in the axes at the following calculated time step t + ∆t. Assuming equality we obtain the incremental expression of the (30)

(31)

and the linearized equation of the principle of virtual work (virtual displacement) simplifies to:

(32)

Assuming we obtain

(33)

where and (34)

(35)

(36)

(37)

Then the equation of the principle of virtual work can be written as follows:

(38)

where

(39)

(40)

After transformation into global coordinate system

, where (41)

, where (42)

and after elimination of the vector of virtual displacements we get:

(43)

(44)

(45)

(46)

(47)

The geometric stiffness matrix (45) is the same as that obtained by use the standard Formula (27) and the first row of the matrix does not correspond with Formula (12). Let us try to derive the geometric stiffness matrix of a truss element using a more accurate strain measure.

The approximate nature of the linear relation between the deformation and displacement can be shown on a fibre of initial length. Without any loss of generalization, let us introduce a system of coordinates with the origin at the starting point of the fibre and with the x axis oriented in the original direction of the fibre. Let us denote by the length of the fibre in the deformed body (Figure 3).

Let us denote by the vector of displacement of the starting point of the fibre. The end-point of the fibre will be displaced by vector.

Using the formula for the body-diagonal of a cuboid with dimensions, , , we can express the new length of the fibre using the following relation:

(48)

Introducing stretch and considering

(49)

we obtain the following relation for stretch of the fibre:

(50)

Figure 3. Elongation of fibre dS.

Let us consider the binomial theorem:

for (51)

and let us take into account only the first two terms. Then we can write:

(52)

and for

(53)

If we want to be more accurate and take into account three terms of the binomial expansion, and if we neglect the third and higher powers of the derivatives of the displacement components, we get a more accurate expression for the stretch:

(54)

and hence

(55)

For a 1D problem, therefore, this more accurate expression would be identical to the formula for known from linear mechanics:

(56)

Using the more accurate strain measure we obtain:

(57)

where,

(58)

where is defined a new matrix instead of the standard.

The linearized equation of the principle of virtual work (virtual displacement) modifies to:

(59)

After transformation into global coordinate system and elimination of the vector of virtual displacements we get different geometric stiffness matrix in the rotated and thus also in global coordinate system:

(60)

(61)

Resulting stiffness matrix derived from the principle of virtual work, using the more accurate strain measure, is the same as that derived from equilibrium conditions (18) and corresponds with Formula (12).

It can be seen that the standard formula has produced a different geometric matrix for the 2D truss element (27) than Formulae (18), (12) and (61) derived earlier and theoretically unjustified geometric axial stiffness was also produced. This formula would lead to a poor convergence rate, inaccuracy and even, in the case of extreme compression, to singularity. E.g. for, zero normal tangent stiffness would be obtained for the truss element, although there is no physical reason for this. For the normal tangent stiffness would even be negative, which would be absurd. In the case of tension no stability problem would occur, but the low convergence problem is still present. E.g. when, the unbalanced nodal forces of 1/2 of the load increment value would occur in the first iteration of the last increment. In the 2^nd iteration it would be 1/4, and in the i-th iteration the unbalanced force of of the load increment value would still occur. These problems are known, and therefore for the geometric stiffness of truss elements Formula (18) is widely used instead of Formula (27), which is derived from the general Formula (8) or (9). Then, in many computer programs different rates of convergence are obtained for a rod modeled by a truss element than in the case of a truss modeled by solid elements.

To obtain the same geometric stiffness matrix for the 2D truss element (18) as was derived above from the equilibrium, the influence of the member must be omitted in the standard formula, i.e. the first row of the matrix must be filled in with zeros.

4. An Improved Formula for a Geometric Stiffness Matrix

Introducing a fibre of constant cross section area A in the direction of principal stress in a 2D or 3D continuum instead of a rod, and assuming only nonzero strain in the direction of the fibre, and that all the other components of the strain tensor are zero, we can write a similar formula to (12):

(62)

where represents the a principal stress.

To evaluate the first part of the expression, a strain measure and pertinent constitutive relation must be chosen. This part represents material stiffness. The second part, which is the matter of our interest, represents geometric stiffness.

The contributions of the two remaining principal stresses to the stiffness could be derived in a similar manner.

Let us introduce the infinitesimal volume element of continuum, being the coordinate system in the principal stress axes given by the relation.

It was earlier shown that for a rod (see formula (12)) the first derivative of a displacement component with respect to the same direction does not generate geometric stiffness. For the 2D or 3D continuum a similar formula to (60) can be derived in a similar way as in the case of a rod.

New measure of deformation defined in the principal axes can be, similarly as in the case of rod, divided into the linear and nonlinear parts as follows:

(63)

where

(64)

being the infinitesimal strain and

(65)

The linearized equation of the principle of virtual work (virtual displacement) for 2D or 3D continuum is similar to (32) and reads:

(66)

yielding its following form in terms of finite element matrices:

(67)

(68)

where

(69)

The difference from the standard formula (8) lies in the fact that in the expression the members, and are omitted.

A particular case where the standard formula was applied to a 2D truss element producing an unintentional change in the axial stiffness was presented earlier. This phenomenon can also be generally observed when the standard formula is used. It is clear that the uniaxial stress state will provide the same result regardless of the way it is modeled, i.e. a truss member modeled as a 3D solid should provide the same result as one modeled by a truss member or by shell elements. To guarantee this and to improve the influence of the stress state on stiffness,

the members, and must be omitted in the standard formula. There is no reason why a normal

stress component should influence the stiffness in the same direction. To ensure objectivity (independence from any arbitrary coordinate system) of the geometric stiffness matrix, the omission of the above mentioned terms must be evaluated in the principal stress axes. Then, for the updated Lagrangian formulation, the following formulae for 3D solid elements hold:

(70)

(71)

where

(72)

being the diagonal unit matrix of the order 2 × 2; is the stress tensor in the principal stress axes with, , being the principal stresses. Then, the geometric stiffness matrix in the axes of principal stresses is defined by the following formula:

(73)

The component of the matrix relating the element node a to the element node b can also be written in indicial notation:

(74)

To obtain the geometric stiffness matrix in the global coordinate system the following transformation must be performed:

The transformation from the global coordinate system to the principal stress coordinate system can be defined as follows:

(75)

Then, the relations between the first derivatives of the base functions and stresses are the following:

(76)

(77)

Illustration on a Quadrilateral Plane Stress Element

For a quadrilateral plane stress element the following can be obtained:

(78)

(79)

(80)

(81)

(82)

where C = cos(α); S = sin(α); α is the angle between principal and global directions; t is the element thickness and A is the area of the element.

A similar formula also holds for the total Lagrangian formulation for such an element, but the second Piola- Kirchhoff stress tensor is then used instead of the Cauchy stress, integration is carried out on the undeformed body (in the original configuration), and the derivatives are performed by the material coordinates.

5. Examples

An application of the new formula for the geometric stiffness matrix for large strain was demonstrated on the example of a unit cube represented by different computational models (rod, shell, solid elements) with different orientations in space (see Figure 4) assuming isotropic hyperelastic material with linear relation between the

Figure 4. Different computational models (rod, shell, solid elements) with different orientations in space.

Figure 5. Magnitude of displacement due touniform normal load of the value E acting on two opposite sides.

Figure 6. Magnitude of displacement due touniform normal load of the value-E acting on two opposite sides.

logarithmic strain and Cauchy stress tensors. Let E be the Young modulus and for simplicity let us assume zero Poisson ratio. The cube was exposed to uniform stress of the magnitude E or −E normal on two opposite sides. A logarithmic strain value of 1 or −1 and the prolongation or shortening of the value of or (e being the base of the natural logarithm) should be obtained. Different computational models of the cube were tested, utilizing rod, shell and solid elements. Calculations were performed for three orientations in space (the basic configuration, a rotation of 30 degrees and a rotation of 45 degrees). Several orientations in space were also tested. In Figure 5 and Figure 6, which are graphical outputs from the RFEM program, it is shown that practically exact results were obtained for all computational models and orientations.

5.1. Convergence of the Standard and New Approach―A Comparison

Let the sequence. converge to. If such a positive number and constant exists that:

(83)

then p is called the order of convergence of the sequence. The constant is called the asymptotic error.

If p is large, the sequence converges rapidly to. If, the convergence is said to be linear. If p = 2, then convergence is quadratic, and if p = 3, it is cubic, etc. Most sequences converge linearly or quadratically. Quadratic convergence is sufficient for computationally efficient numerical methods.

5.2. The Standard Approach

Figure 7 shows that the standard approach provides only linear convergence for large strain.

Figure 7. Investigation of linear (, upper panel) and quadratic (, lower panel) convergence, where.

5.3. The New Approach

Figure 8 shows that the new approach yields the quadratic convergence even for very large strain.

Figure 8. Investigation of linear (, upper panel), quadratic (, middle panel) and cubic (, lower panel) convergence, where.

The numerical solution of the presented example has shown that to reach a sufficiently good result using the standard formula (ANSYS etc.) 15 iterations were needed whereas using the improved approach presented in this paper (RFEM) only 5 iterations were needed to obtain the same precision.

6. Conclusions

The present formula for a geometric stiffness matrix, which has been published in many books, is widely utilized, objective and simply defined. However, stability and convergence problems occur when analyzing large strains, or, what is more important in practice, in a case of yielding. If the yield criterion is satisfied, then the

material stiffness decreases substantially. The stress state remains high and in case of compression the tangent stiffness in the direction of the compression can become negative even with a small strain.

This is caused by a theoretically unsupported change in pressure stiffness in the direction of compression produced by the standard formula. This results in a correspondingly high nodal force unbalance, poor convergence and possibly also instability. The origin of the problem arises from the approximation of strain, in which only the first two terms of the binomial series are applied.

The presented algorithm is slightly more complicated, but remains objective and provides a solution with increased stability, a higher rate of convergence in the case of a large strain, or plastic yielding, and improved accuracy over the present formula. In case of very large strain, the number of iterations needed could be several times less using the new formula comparing to the standard formula. In many cases the new formula can even provide solutions in cases where the standard formula has failed. This new formula for a geometric matrix has been implemented in the RFEM program and has been demonstrated to be much more stable and faster than the standard formula.

Acknowledgements

This outcome has been achieved with the financial support of the Czech Science Foundation (GACR) project 14-25320S “Aspects of the use of complex non linear material models”.

List of Variables

Cross section area of a beam

Material tangent moduli

Young modulus

Material and geometric tangent stiffness matrix, respectively

Material and geometric tangent stiffness matrix in principal axes

Shape functions

Matrix of shape functions

Rotation tensor

Second Piola-Kirchhoff stress

Internal and external virtual work

Vector of nodal displacements

Infinitesimal strain

Infinitesimal strain in principalaxes

Internal nodal forces

External nodal forces

Unit diagonal matrix

Member length

Time

Displacement field

Displacements in the x, y and z directions respectively

Spatial (Eulerian) coordinates

Coordinates in principal aces

Modified strain tensor in principalaxes

Quadratic terms of the modified strain tensor in principalaxes

Potential energy of geometrical stiffness

Cauchy stress tensor

Cauchy stress tensor in principal axes

Domain of current (deformed), initial (undeformed)

NOTES

^*Corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Bathe, K.-J. (1996) Finite Element Procedures. Prentice Hall, Upper Saddle River.
[2]	Belytschko, T., Liu, W.K., Moran, B. and Elkhodary, K. (2000) Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, Hoboken.
[3]	Bonet, J. and Wood, R.D. (2008) Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press, Cambridge.
[4]	Cook, R.D., Malkus, D.S., Plesha, M.E. and Witt, R.J. (2002) Concepts and Applications of Finite Element Analysis. John Wiley & Sons, Hoboken.
[5]	Němec I., et al. (2010) Finite Element Analysis of Structures: Principles and Praxis. Shaker-Verlag, Aachen.
[6]	Reddy, J.N. (2004) Nonlinear Finite Element Analysis. Oxford University Press, Oxford.
[7]	Zienkiewicz, O.C. and Taylor, R.L. (2000) The Finite Element Method. Butterworth Heinemann, Oxford.
[8]	de Souza Neto, E.A., Periæ, D. and Owen, D.R.J. (2008) Computational Methods for Plasticity: Theory and Applications. John Wiley & Sons, Hoboken.
[9]	Simo, J.C. & Hughes, T.J.R. (2008) ComputationalInelasticity. Springer, New York, 392 p.
[10]	Wriggers, P. (2008) Nonlinear Finite Element Methods. Springer-Verlag, New York.
[11]	Lacarbonara, W. and Pacitti, A. (2008) Nonlinear Modeling of Cables with Flexural Stiffness. Mathematical Problems in Engineering, 2008, Article ID: 370767. http://dx.doi.org/10.1155/2008/370767
[12]	Xiao, H. and Chen, L.S. (2002) Hencky’s Elasticity Model and Linear Stress-Strain Relations in Isotropic Finite Hyperelasticity. Acta Mechanica, 157, 51-60. http://dx.doi.org/10.1007/BF01182154
[13]	Xiao, H., Bruhns, O.T. and Meyers, A. (1998) Objective Corotational Rates and Unified Work-Conjugacy Relation between Eulerian and Lagrangean Strain and Stress Measures. Archives of Mechanics, 50, 1015-1045.
[14]	Meyers, A. (1999) On the Consistency of Some Eulerian Strain Rates. Zeitschrift fur Angewandte Mathematik und Mechanik, 79, 171-177. http://dx.doi.org/10.1002/(SICI)1521-4001(199903)79:3<171::AID-ZAMM171>3.0.CO;2-6
[15]	Simo, J.C. and Pister, K.S. (1984) Remarks on Rate Constitutive Equations for Finite Deformation Problem: Computational Implications. Computer Methods in Applied Mechanics and Engineering, 46, 201-215. http://dx.doi.org/10.1016/0045-7825(84)90062-8
[16]	Curnier, A. and Rakotomanana, L. (1991) Generalized Strain and Stress Measures: Critical Survey and New Results. Engineering Transactions, 39, 461-538.
[17]	Chiskis, A. and Parnes, R. (2000) Linear Stress-Strain Relations in Nonlinear Elasticity. Acta Mechanica, 146, 109-113. http://dx.doi.org/10.1007/BF01178798
[18]	Farahani, K. and Naghdabadi, R. (2000) Conjugate Stresses of the Seth-Hill Strain Tensors. International Journal of Solids and Structures, 37, 5247-5255. http://dx.doi.org/10.1016/S0020-7683(99)00209-7
[19]	Darijani, H. and Naghdabadi, R. (2010) Constitutive Modeling of Solids at Finitede Formation Using a Second-Order Stress-Strain Relation. International Journal of Engineering Science, 48, 223-236. http://dx.doi.org/10.1016/j.ijengsci.2009.08.006
[20]	Hill, R. (1978) Aspects of Invariance in Solid Mechanics. Advances in Applied Mechanics, 18, 1-75. http://dx.doi.org/10.1016/S0065-2156(08)70264-3
[21]	Farahani, K. and Bahai, H. (2004) Hyper-Elastic Constitutive Equations of Conjugate Stresses and Strain Tensors for the Seth-Hill Strain Measures. International Journal of Engineering Science, 42, 29-41. http://dx.doi.org/10.1016/S0020-7225(03)00241-6