New Formula for Geometric Stiffness Matrix Calculation

Journal of Applied Mathematics and Physics, 2016, 4, 733-748 Published Online April 2016 in SciRes. http://www.scirp.org/journal/jamp http://dx.doi.org/10.4236//jamp.2016.44084 How to cite this paper: Němec, I., Trcala, M., Ševčík, I. and Štekbauer, H. (2016) New Formula for Geometric Stiffness Matrix Calculation. Journal of Applied Mathematics and Physics, 4, 733-748. http://dx.doi.org/10.4236//jamp.2016.44084 New Formula for Geometric Stiffness Matrix Calculation I. Němec1*, M. Trcala2, I. Ševčík2, H. Štekbauer1 1Faculty of Civil Engineering, Brno University of Technology, Brno, Czech Republic 2FEM Consulting, S.R.O., Brno, Czech Republic Received 4 June 2015; accepted 24 April 2016; published 27 April 2016 Copyright © 2016 by authors and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY). http://creativecommons.org/licenses/by/4.0/ 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 * Corresponding author. I. Němec et al. 734 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 i u of displacement vector u : 1 n i a ia a u Nu = = ∑ (1) 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: [ ] 12 , ,,n NN N=N II I (2) where I 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: = ⋅u Nd (3) 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 , , T , ax a a ay az N NN N  ∂ = ⊗=  ∂  I g II xI (4) and matrix G , which is formed by sub-matrices a g [ ] 12 T ,,,,, an ∂ = = ∂ N G gg g g x  (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: sym. xx xy xz yy yz zz σσσ σσ σ   =⊗=    III I II I Σ σ (6) I. Němec et al. 735 If state of the stress is not negligible, the potential energy of the internal forces should be completed by the following term: ( ) TT TT TT 1 11 dd 2 22 σσ Ω Ω ⌠  ⌡ ∂ ∂ ∂∂ ∏ = ⊗ Ω= Ω = ∂∂ ∂∂ ∫u u NN I d d dKd xx xx Σ σ (7) Then, the following formula for the geometric matrix of the element can be written: T T T dd σ Ω Ω ⌠  ⌡ ∂∂ = Ω= Ω ∂∂ ∫ NN K GG xx ΣΣ (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 σ K relating the element node a to the element node b can also be written simply in matrix notation: ( ) d ab a b NN σ Ω = ∇ ⋅∇ Ω ∫ KσI (9) or in indicial notation: d , 1, 2, 3 ab abij kl ij kl NN K ij xx σ σδ Ω ⌠  ⌡ ∂∂ = Ω= ∂∂ (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 0 Ω (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 ( ) ( ) 0Rx Tx F= −= (11) where ( ) cosTx N α = is the horizontal component of the internal force at node 2 and cos xl α = . () Rx is the residual or out-of-balance force. The horizontal stiffness x K at node 2 is defined simply by the relation 2 2 22 22 22 dd d d d d d ddd d d d d dd 1 cos sin dd x xM x R T Nx N N N l N N N x N K xx x x x l xl l ll x l l l l l Nx N x N N KK ll l l ll σ αα       = = = = += += − +              = + −= + =+   (12) 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: ∂ =∂ r Kd (13) I. Němec et al. 736 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 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 [ ] T 11 2 2 ,, ,uvu v=d (14) 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: 1 0 10 0000 10 1 0 0000 M EA l −    = −   K (15) 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: d d cos d 0N v Tl α −= (16) I. Němec et al. 737 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: 22 d d TN Kvl σ = = (17) 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: 0000 010 1 0000 0 10 1 N l σ   −  =  −  K (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 T K is defined as the sum of the material and geometric stiffness matrix: TM σ = +KK K (19) When applying the general standard algorithm for geometric stiffness matrices to the truss element in question, we obtain: u v  = = ⋅   u Nd (20) [] 12 22 ,NN= N II (21) where 2 I is the identity matrix of order 2 and the base functions i N are defined as follows: 12 1xx NN ll    = − , =       (22) xx ∂∂ = = ∂∂ uN d Gd (23) where 1 0 10 1 0 101 xl −  ∂ = =  − ∂ N G (24) Substituting in the formulae x N A σ = =Σ (25) AlΩ= ⋅ (26) the formula for the geometric stiffness matrix reads: TT 1 0 10 010 1 dd 10 1 0 0 10 1 x l N Nl l σ σ Ω −   −  = Ω= =  −  −  ∫∫ K G G GG (27) I. Němec et al. 738 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 22 1 2 x xx u uv e x xx εη  ∂ ∂∂  =+ +=+    ∂ ∂∂   (28) where 22 1 ,2 xx u uv ex xx η  ∂ ∂∂  = = +    ∂ ∂∂   (29) For truss the principle of virtual work becomes d x x ext SW δε Ω Ω= ∫ (30) where x S is the 2nd Piola-Kirchhoff stress in the x axes at the following calculated time step t + ∆t. Assuming equality tt t x xxxxx S SS S S σ + ≡ = +∆ = +∆  we obtain the incremental expression of the (30) dd d x x x x ext x x S We δε σ δη σ δ ΩΩ Ω ∆ Ω+ Ω= − Ω ∫∫ ∫ (31) and the linearized equation of the principle of virtual work (virtual displacement) simplifies to: dd d x x x x ext x x Ee e W e δ σδη σδ ΩΩ Ω Ω+ Ω= − Ω ∫∫ ∫ (32) Assuming (25) we obtain xx x x EAe el N l F u N el δ δη δ δ +=− (33) where x u ex δ δ ∂ =∂ and x uu vv xx xx δδ δη ∂∂ ∂∂ = + ∂∂ ∂∂ (34) [ ] 11010 x u exl ∂ = = − ∂d (35) [ ] TT 2 1 1 0 10 0 0000 11 1 1010 1 10 1 0 0 0000 xx ee ll l δδ δ −−          =−=    −       d dd d (36) TT T 2 1 0 10 010 1 1 10 1 0 0 10 1 x uu vv xx xx l δδ δη δ δ −   − ∂∂ ∂∂  =+= =  − ∂∂ ∂∂  −  d G Gd d d (37) Then the equation of the principle of virtual work can be written as follows: T T TT M ext int σ δ δ δδ +=− dKd dKd d f d f (38) where 1 0 10 1 0000 0 , 10 1 0 1 0000 0 M int EA N l −−          = =    −       Kf (39) I. Němec et al. 739 T 1 0 10 010 1 d10 1 0 0 10 1 N l σ Ω −   −  = Ω=  −  −  ∫ K GGΣ (40) After transformation into global coordinate system xx yy   =     R , where ( ) ( ) ( ) ( ) cos sin sin cos CS SC αα αα  − −  = =     R (41) =d Td , where  =  R TR 0 0 (42) and after elimination of the vector of virtual displacements we get: T MM =K TKT (43) 22 22 22 22 M C CS C CS CS S CS S EA lC CS C CS CS S CS S  −−  −−  = −−  −−   K (44) T 1 0 10 010 1 10 1 0 0 10 1 N l σσ −   −  = =  −  −  K T KT (45) [ ] T T int in NC SCS= =−−f Tf (46) ( ) M ext int σ +=−K Kd f f (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 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+uu . Using the formula for the body-diagonal of a cuboid with dimensions ddSu+ , dv , dw , we can express the new length of the fibre using the following relation: () 222 d dd d d s Su v w= + ++ (48) Introducing stretch ddsS λ = and considering d dd d d uv w u Sv Sw S xx x δ ∂∂ ∂ = , = , = ∂∂ ∂ (49) we obtain the following relation for stretch of the fibre: 22 2 22 2 d 1 1 12 d x s uvw uuvw S xxx xxxx λε ∂∂∂ ∂∂∂∂    =+==+++ =++++    ∂∂∂ ∂∂∂∂    (50) I. Němec et al. 740 Figure 3. Elongation of fibre dS. Let us consider the binomial theorem: 23 11 2 8 16 AA A A+=+− + + for 21A< (51) and let us take into account only the first two terms. Then we can write: 22 2 1 12 u uvw x xxx λ  ∂ ∂∂∂  =++ + +    ∂ ∂∂∂   (52) and for x ε 22 2 1 12 x u uvw x xxx ελ  ∂ ∂∂∂  = −= + + +    ∂ ∂∂∂   (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: 22 1 12 u vw x xx λ  ∂ ∂∂   =++ +     ∂ ∂∂    (54) and hence 22 1 2 x u vw x xx ε  ∂ ∂∂   =++     ∂ ∂∂    (55) For a 1D problem, therefore, this more accurate expression would be identical to the formula for x ε known from linear mechanics: x u x ε ∂ =∂ (56) Using the more accurate strain measure we obtain: 2 1 2 x xx uv e xx εη ∂∂  =+=+  ∂∂  (57) where x u ex ∂ =∂ , 2 1 2 x v x η ∂  = ∂  I. Němec et al. 741 [ ] TT T T new new 2 0 0000 1 010 1 11 1 0 101 0 0000 1 0 10 1 x vv xx l l l δ δη δ δ δ       −− ∂∂    == = −=    ∂∂    −    dGGdd dd d (58) where is defined a new matrix [] new 10 101 l = − G instead of the standard G . The linearized equation of the principle of virtual work (virtual displacement) modifies to: T T TT 10 10 0000 1 0000 010 1 0 1010 0000 1 0000 0 101 0 ext EA N N ll δ δ δδ −−     −   +=−   −   −   d d d d df d (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: T new new 0000 010 1 d0000 0 10 1 N l σ Ω   −  = Ω=   −  ∫ K GGΣ (60) 22 22 T 22 22 S SC S SC SC C SC C N lS SC S SC SC C SC C σσ  −−  −−  = =  −−  −−   K T KT (61) 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 xE σ <− 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 xE σ = , the unbalanced nodal forces of 1/2 of the load increment value would occur in the first iteration of the last increment. In the 2nd iteration it would be 1/4, and in the i-th iteration the unbalanced force of 12 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 ux∂∂ must be omitted in the standard formula, i.e. the first row of the G 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 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): I. Němec et al. 748 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 List of Variables A Cross section area of a beam C Material tangent moduli E Young modulus , M σ KK Material and geometric tangent stiffness matrix, respectively , M σ KK Material and geometric tangent stiffness matrix in principal axes i N Shape functions N Matrix of shape functions R Rotation tensor S Second Piola-Kirchhoff stress , int ext WW Internal and external virtual work d Vector of nodal displacements ˆ e Infinitesimal strain e Infinitesimal strain in principalaxes int f Internal nodal forces ext f External nodal forces I Unit diagonal matrix l Member length t Time u Displacement field ,,uvw Displacements in the x, y and z directions respectively ,,xyz Spatial (Eulerian) coordinates ,,xyz Coordinates in principal aces ε Modified strain tensor in principalaxes ,ij ηη Quadratic terms of the modified strain tensor in principalaxes σ ∏ Potential energy of geometrical stiffness Σ = ⊗ IΣ σ Σ = ⊗ IΣ σ ,ij σσ Cauchy stress tensor σ Cauchy stress tensor in principal axes 0 , ΩΩ Domain of current (deformed), initial (undeformed)