New Formula for Geometric Stiffness Matrix Calculation

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.


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.

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 i u of displacement vector u : where ia u is the value of displacement i u in node a and n is the number of element nodes.Let us define matrix N as follows: [ ] where I is the unit diagonal matrix of the order 3 × 3, where 3 is the dimension of the problem.Then, the fol- lowing relation can be written for the displacement vector: where d is the vector of deformation parameters of the element containing all the components ia u in such an arrangement that for each node a all components i u are listed.Let us define matrix a g containing the first derivatives of base functions for node a with respect to spatial coordinates and matrix G , which is formed by sub-matrices a g [ ] 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: sym.
If state of the stress is not negligible, the potential energy of the internal forces should be completed by the following term: ( ) Then, the following formula for the geometric matrix of the element can be written: Integration is carried out on the deformed body Ω (in the current configuration) and the derivatives are per- formed with respect to the spatial coordinates.
The component of the matrix σ K relating the element node a to the element node b can also be written simply in matrix notation: or in indicial notation: 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 0 Ω (in the orig- inal configuration) while the derivatives are performed with respect to the material coordinates.

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 where ( ) is the horizontal component of the internal force at node 2 and cos x l α = .( ) R x is the residual or out-of-balance force.The horizontal stiffness x K at node 2 is defined simply by the relation This formula is independent of any strain measure or pertinent constitutive relations.It can be seen that stiffness x K consists of two parts.The first part, xM K , represents the material stiffness and depends on the strain measure and constitutive relations.The second part, x K σ , 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 d be a vector of the nodal displacements of an element, and let r be a vector of residual forces; the stiffness matrix of the element can then be defined as follows:  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 x with the x axis aligned with the axis of the rod and corresponding displacement vector u and let's restrict the deformation to the xy plane.
Let the vector of the nodal displacements of the element be where u and v are the displacement components in the direction of the x and y axis, respectively.The well known material stiffness matrix of the truss element in 2D is then defined by the following relation: Note that the truss element has no lateral material stiffness.
In general, arbitrary term of a stiffness matrix ij K is defined as the derivative of an unbalanced force i r with respect to the deformation parameter j d 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 dv in node 1 is sufficient to obtain the transversal diagonal stiffness term 22 K : The moment equilibrium condition can be written as follows: For the infinitesimal angle dα it can be assumed that cos d 1 α = , and the following term for the stiffness term 22 K σ can be derived: When introducing a displacement du 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 24 K σ , 42 K σ and 44 K σ .The remaining coefficients of the matrix are zeros.The geometric stiffness matrix then has the following form: 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 T K is defined as the sum of the material and geometric stiffness matrix: When applying the general standard algorithm for geometric stiffness matrices to the truss element in question, we obtain: [ ] where 2 I is the identity matrix of order 2 and the base functions i N are defined as follows: where Substituting in the formulae the formula for the geometric stiffness matrix reads: 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 xy plane, the Green-Lagrange strain tensor is defined 2 2 1 2 where For truss the principle of virtual work becomes where x S is the 2nd Piola-Kirchhoff stress in the x axes at the following calculated time step t + ∆t.Assuming equality we obtain the incremental expression of the (30) and the linearized equation of the principle of virtual work (virtual displacement) simplifies to: Assuming (25) we obtain where [ ] Then the equation of the principle of virtual work can be written as follows: where After transformation into global coordinate system , where and after elimination of the vector of virtual displacements we get: 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 dS .Without any loss of generalization, let us introduce a system of coordinates x 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 ds 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 d + u u.Using the formula for the body-diagonal of a cuboid with dimensions d d S u + , dv , dw , we can express the new length of the fibre using the following relation: we obtain the following relation for stretch of the fibre: Let us consider the binomial theorem: and let us take into account only the first two terms.Then we can write: and for x ε 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: and hence For a 1D problem, therefore, this more accurate expression would be identical to the formula for x ε known from linear mechanics: Using the more accurate strain measure we obtain: where is defined a new matrix The linearized equation of the principle of virtual work (virtual displacement) modifies to: 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: Resulting stiffness matrix σ K 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 x E σ = − , zero normal tangent stiffness would be obtained for the truss element, although there is no physical reason for this.For x E σ < − 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 x E σ = , the unbalanced nodal forces of 1/2 of the load increment val- ue 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 1 2 i 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 u x ∂ ∂ must be omitted in the standard formula, i.e. the first row of the G matrix must be filled in with zeros.

An Improved Formula for a Geometric Stiffness Matrix
Introducing a fibre of constant cross section area A in the direction x 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): where x σ 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 d d d d x y z Ω = , x being the coordinate system in the principal stress axes given by the relation It was earlier shown that for a rod (see formula 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: where ê being the infinitesimal strain and ( ) The linearized equation of the principle of virtual work (virtual displacement) for 2D or 3D continuum is similar to (32) and reads: ( yielding its following form in terms of finite element matrices: where The difference from the standard formula (8) lies in the fact that in the Gd expression the members 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 must be omitted in the standard formula.There is no reason why a normal stress component should influence the stiffness in the same 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 x .Then, for the updated Lagrangian formulation, the following formulae for 3D solid elements hold: , , , , , where I being the diagonal unit matrix of the order 2 × 2; σ is the stress tensor in the principal stress axes x with x σ , y σ , z σ being the principal stresses.Then, the geometric stiffness matrix in the axes of principal stresses is defined by the following formula: The component of the matrix σ K relating the element node a to the element node b can also be written in indicial notation: To obtain the geometric stiffness matrix in the global coordinate system the following transformation must be performed: The transformation from the global coordinate system x to the principal stress coordinate system x can be defined as follows: = x Rx (75) Then, the relations between the first derivatives of the base functions and stresses are the following:

Illustration on a Quadrilateral Plane Stress Element
For a quadrilateral plane stress element the following can be obtained: 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 0 Ω (in the original configuration), and the derivatives are performed by the material coordinates.

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   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 1 e − or 1 1 e − (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.

Convergence of the Standard and New Approach-A Comparison
Let the sequence { }

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

The New Approach
Figure 8 shows that the new approach yields the quadratic convergence even for very large strain.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.

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.

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

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

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

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

Figure 6 .
Figure 6.Magnitude of displacement due touniform normal load of the value-E acting on two opposite sides.
then p is called the order of convergence of the sequence.The constant λ is called the asymptotic error.If p is large, the sequence { } n u converges rapidly to c u .If 1 p = , 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.