Numerically Generated Tangent Sti Ness Matrices For Geometrically Non .

4tangent sti ness matrix with results obtained in literature using analytically calculated tangent sti ness matrices or using commercial finite element software. Thepresent solution method was coded in MATLAB in order to solve the examples. Apseudocode for the solution method is provided in Appendix A.Finally, Chapter 6 is a concluding discussion of the results and possible directionsfor future improvements and research.

5Chapter 2SOLUTION METHOD2.1Tangent sti ness matrixThe solution method for a geometrically non-linear problem begins with a linearizationof the non-linear equilibrium equations, such as the general force equilibrium equationsgiven by,Pext Fint (u)(2.1)where Fint is a vector of nodal internal member forces, which are functions ofthe nodal degree of freedom (DOF) displacements1 , u. Pext is a vector of externallyapplied loads.2 Assuming there is a known set of nodal displacement DOFs, u0 , thatsatisfies equation (2.1), the equilibrium equations can be linearized by perturbing theforce about this known solution point. A small perturbation of the externally appliedload corresponds to a perturbation in the nodal DOF displacements and equation(2.1) becomesPext dP Fint (u0 du)(2.2)A first-order taylor series expansion of the right hand side of equation (2.2) resultsin1The term “displacements” is used in this chapter to denote both linear displacements and rotations.2Generally, the external forces may be follower forces and therefore will also be dependant on thenodal displacement degrees of freedom. However, for the purpose of simplicity, follower forces arenot considered in this derivation.

6Pext dP Fint (u0 ) @Fint@udu(2.3)u u0Since u0 is a solution that satisfies equation (2.1), equation (2.3) reduces todP @Fint@udu [KT (u0 )]du(2.4)u u0where266@Fint 6[KT (u)] t,2@um.@Fint,m@um37777775(2.5)is the tangent sti ness matrix, with m being the total number of nodal DOFs inthe structure. From equation (2.5) it can be seen that the tangent sti ness matrix isessentially a matrix of sensitivities. In particular, it is the sensitivities of the internalmember forces to perturbations in the nodal displacement DOFs of the system. Thetangent sti ness matrix contains information regarding both the linear elastic andgeometric sti ness of the structure.2.2Incremental-iterative solution technique: Netwon-Raphson methodAs with any set of non-linear equations, the most e ective way to solve the geometrically non-linear problem is to use an incremental-iterative technique. For the presentwork a Newton-Raphson method is used. The externally applied load is divided intoincrements,P, and at each increment a Newton-Raphson iteration is run in orderto converge on the equilibrium condition, as illustrated in Figure 2.1.

7Figure 2.1: Newton-Raphson incremental-iterative solution methodGiven an initial configuration for the structure, and thus the initial nodal DOFdisplacement vector u0 , the initial tangent sti ness matrix, [KT (u0 )], can be calculated. The externally applied load is then incremented by some small increment,P,and the linearized system of equations from equation (2.4) is solved for the nodalincremental displacement due to the external load:u11 [KT (u0 )]1P [KT (u0 )] 1 Pext,1(2.6)

8Here the subscripts represent the increment number, and the superscripts representthe iteration number.The internal member forces at the updated configuration, Fint (u0 u11 ), arecalculated, and the equilibrium condition from equation (2.1) is checked. That is, forequilibrium to be satisfied,Pext,1Fint (u0 u11 ) 0(2.7)However, since a linear approximation was used to solve for the nodal displacements, it is likely that there is an unbalance in the force equilibrium, otherwise knownas a residual force R. The equilibrium equation is then,Pext,1Fint (u0 u11 ) R11(2.8)and iteration proceeds to converge on a solution for the nodal DOF displacementsthat results in the residual being as close to zero as possible. This is done by calculating the updated tangent sti ness matrix from the updated nodal displacementsand solving for the nodal displacement increments due to the residual force,u21 [KT (u0 u11 )] 1 R11(2.9)The updated internal member forces, Fint (u0 u11 u21 ), are again calculated,and the equilibrium equation becomes,Pext,1Fint (u0 u11 u21 ) R21(2.10)

9If the residual is not close to zero within an acceptable tolerance, then the processis repeated. If the residual is approximately zero, then the equilibrium condition hasbeen achieved and the new updated equilibrium configuration due to the externallyapplied load of Pext,1 is u1 u0 ( u11 u21 · · · ui1 ) u0 u1 , where i isthe number of iterations required to achieve the equilibrium condition.The solution process continues by again incrementing the externally applied loadsuch that Pext,2 Pext,1 P, and proceeding with the iteration on the linearsystem of equations. This is repeated until the full applied load is reached, with thefinal solution being uj u0 u1 u2 · · · uj , where j is the total numberof load increments.While it is one of the simplest iterative methods to implement, the NewtonRaphson method is not the most numerically robust. In particular, for bucklingproblems, or structures with snap-through tendencies, there may be a point at whichthe tangent sti ness matrix becomes singular and the solution ‘blows up’. For thisreason, more sophisticated incremental-iterative methods are often used, includingthe arc length method, the generalized displacement control method as presented in[21], or the BFGS iteration method as presented in [1]. However, for this currentwork, only simple structural problems are analyzed for which the Newton-Raphsonmethod is acceptable.2.3Numerically Generated Tangent Sti ness MatrixCalculating the tangent sti ness matrix is integral to finding a solution to the geometrically non-linear problem. In order to develop a solution method that is generaland easily implementable, the present work focuses on a method for calculating thetangent sti ness matrix numerically.As shown in equation (2.5), the elements of the tangent sti ness matrix are sensitivities of the internal member forces to perturbations in the global nodal DOFdisplacements of the structure. There are many di erent numerical di erentiation

10techniques that could be employed to calculate these sensitivities, most common being any finite di erence method, such as the central di erence method. However,while finite di erence methods are typically easy to implement, step size often becomes an issue with these methods. It is desirable to choose a smaller step size inorder to obtain greater accuracy, but with a finite di erence method a step size thatis too small can result in subtractive cancellation round-o error and numerical instability. This type of error can be completely avoided, however, if the first derivative isestimated using a complex variable di erentiation method.An approximation of the first derivative using complex variable analysis beginsby defining a complex function, f (z) f (x iy) u(x, y) iv(x, y), where u and vare real functions. Combining the Cauchy-Riemann equations derived from complexfunction di erentiation,@u@v @x@y@u@v @y@x(2.11a)(2.11b)and the definition of a derivative, yields@u@vv(x i(y h)) lim@x@y h!0hv(x iy).(2.12)Since we are dealing with real structures, where the forces and degrees of freedomare real values, we restrict ourselves to the real axis, meaning that y 0, so v(x, y) v(x) 0 and f (x iy) f (x) u(x). Equation (2.12) then becomes@f (x)v(x ih)Imag[f (x ih)] lim lim.h!0h!0@xhh(2.13)To approximate the derivative, h is assumed to be a very small finite step size,and equation (2.13) becomes

11@f (x)Imag[f (x ih)] .@xh(2.14)Detailed derivation of the ‘complex-step derivative approximation’ and discussion regarding its numerical robustness compared to other approximation methods ispresented in [10], [11], and [20]. Of particular importance is that complex variabledi erentiation is relatively easy to implement while still allowing for a very small stepsize that translates into very high accuracy without the problem of round-o errors.Equation (2.14) can then be used to approximate the sensitivities of the internalmember forces to perturbations in the global nodal DOFs, and thus used to find thecomponents of the tangent sti ness matrix. The code, instead of being written usingreal variables, is written using complex variables for all entities involved.To achieve this, the first step is to numerically perturb each global nodal DOF.This is done by defining a nodal displacement vector populated with zeros, except forthe component of the vector corresponding to the perturbed DOF, which is updatedfrom zero to ih. This is equivalent to saying that only one particular DOF has beenslightly perturbed while the rest have remained fixed. The resulting vector is thenodal displacement perturbation vector,duperturbed 8 00.9 0 ih . . . : 0 ;(2.15)where h is the chosen step size for the complex variable derivative approximation.With the global nodal perturbation vector established, the global DOFs for the

12structure are updated as follows:uperturbed u duperturbed(2.16)The updated nodal DOFs in the local element axis system can then be calculated.The specific calculation of the local nodal DOF motions from the global nodal DOFmotions is dependant on the particular finite element being used in the structuralproblem, and is typically a known procedure from linear-elastic analysis. It is discussed in further detail for three particular element types in following chapters of thiswork.Next, the internal member forces resulting form the perturbed local nodal DOFmotions are calculated. Since the internal structural forces are a result of elasticdeformation in the member, the linear-elastic equations are used to calculate theforces as follows:[kE ]uperturbed, local fint,perturbed(2.17)where [kE ] is the linear-elastic sti ness matrix for the particular structural element, and fint are the internal elastic forces referred to in the local element axissystem. Since fint are the forces resulting from a complex variable perturbation, thevalues of the force vector will themselves be complex values (the linear-elastic sti ness matrix is of course real valued, but complex algebra is used in the computercode). Note that while large motions are allowed, it is assumed that in their ownlocal coordinates, elements deform only by small deformations. This justifies usinglinear force-deformation relations in those element coordinate systems.The local internal forces are transformed from the element local axis system tothe global axis system using a transformation matrix, [T ],

13Fint,perturbed [T ]T fint,perturbed(2.18)resulting in a complex valued vector representing the nodal member forces in theglobal axes due to a complex variable perturbation of one of the global nodal DOFs.Equation (2.14) is used to numerically approximate the sensitivities of these member forces to the perturbation of the specific DOF as follows:89@Fint,1 @up @Fint,2 @Fint,perturbedImag[Fint,perturbed ]@up . @uph . @F: int,m ;(2.19)@upwhere m is the total number of DOFs in the system, and p is the index corresponding to the particular perturbed DOF.The vector resulting from equation (2.19) becomes the column of the tangentsti ness matrix corresponding to the particular perturbed DOF, specifically the pthcolumn. The above process is repeated for each global nodal DOF perturbation untilall the columns of the tangent sti ness matrix have been generated. This is easilyexecuted in a code by creating a loop that runs through each of the global DOFs andconstructing the corresponding [KT ] column. This process for numerically generating[KT ] is the same for all structures, regardless of element type, resulting in the desiredgenerality. Note that the resulting tangent sti ness matrix is in general not symmetric.

14Chapter 3TRUSS ELEMENTSThe following is a detailed derivation of the present solution method applied to geometrically non-linear three-dimensional structures constructed from truss elements.The derivation is followed by examples with calculated results that are compared withresults from the literature.3.1Solution MethodA simple two-noded truss element is used for analyzing geometrically non-linear threedimensional truss structures. This element carries axial load only, so in the localelement axis has only one degree of freedom at each node. However, when transformedinto the global axis system, each node has three degrees of freedom - displacement inthe X, Y, and Z directions, as illustrated in Figure 3.1.Figure 3.1: Truss element nodal degrees of freedom in local and global axis systems.

15Element DefinitionThe geometry and movement of the truss element in three-dimensional space is completely defined by the location of the two nodes. The node locations are obtainedfrom the location of the nodes in global coordinates. The coordinates are initiallydefined for the undeformed structure, and as the structure deforms the incrementaldisplacements at the nodes are determined by solving the linear equation[KT ] u (PextFint )(3.1)where,u { X1Y1Z1···XmYmZm }Tand m is the total number of nodes in the structure. The global coordinates areupdated from the incremental nodal displacements as illustrated in Figure 3.2.Figure 3.2: Updating of truss element node locations from underformed to deformedconfiguration.

16With the global node locations determined, each truss element in the structure iscompletely defined in space.Force Calculation and Solution ProcessFor a geometrically non-linear structure, finding the increment in internal memberforces as the structure changes geometry is instrumental to the solution process.The process of finding the internal member forces is also an integral component ofthe numerical tangent sti ness generation. The internal forces are a result of linearelastic deformation, which is determined by the motion of the nodes relative to thelocal truss element axis, x̃.The displacements of the nodes relative to the element axis are found by calculating the displacement of node 2 relative to node 1 in the local axis system of eachelement. Therefore, at node 1 the relative displacement isũ1 0(3.2)The displacement of node 2 relative to node 1 is then simply the extension of theelement in the local axis direction, which is the same as the change in element length.Therefore, at node 2 the relative displacement isũ2 LnewLinitial(3.3)where the length of the element, L, is calculated from the node locations in globalcoordinates as follows:L p(X2X1 )2 (Y2Y1 )2 (Z2Z1 ) 2

17Thus, the local elastic deformation vector for the truss element is,uelastic 8 9 ũ 1:ũ ;2and is illustrated in Figure 3.3. 8 :L9 0newLinitial ;(3.4)Figure 3.3: Elastic deformation of truss element in the local axis system.The internal forces in the truss element resulting from the elastic deformation arethen found using the linear-elastic equation[kE ]uelastic felement(3.5)where [kE ] is the local linear-elastic sti ness matrix for the truss element, whichis given as[kE ] 2AE 4 1Linitial11135(3.6)

18where A is the cross-sectional area of the truss element, and E is Young’s Modulusfor the material.The forces calculated in equation (3.5) are due only to the change in geometry ofthe structure and do not include the e ects of any pre-existing stresses in the elements.Therefore, if the structure was prestressed, thus having pre-existing internal elementforces, it is necessary to account for them by adding them to the calculated forces. Thepre-existing internal forces do not change as the structure deforms, but do contributeto the overall nodal forces in the structure.The local element forces, felement , are transformed into the global coordinate system as follows:89 Fx,1 Fy,1 F z,1 Fx,2 Fy,2 :F ;z,2 [T ]T8 9 f 1:f ;2(3.7)elementglobalwhere [T ] is the transformation matrix from the local axis system to the globalaxis system. The transformation matrix for the truss element is made up of directioncosines calculated from the global node coordinates, and typically has the followingform:2[T ] 4mxmymz000000mxmymz35(3.8)where the direction cosines are,mx X2X1L;my Y2Y1L;mz Z2Z1L(3.9)

19The global forces are then assembled into the full internal force vector for thestructure, Fint , and the unbalanced load can then be calculated in order to checkfor equilibrium, as discussed in Section 2.2. The process for calculating the internalmember forces as outlined above is also required for numerically generating the tangent sti ness matrix, which is done using a nodal displacement perturbation vector,as described in Section 2.3.3.2Numeric ResultsThe solution method for the present work involving the numerically generated tangentsti ness matrix is implemented in a MATLAB code, which is then used to obtainresults for two di erent truss structure examples. The same truss structures are solvedby Levy and Spillers in [9], using an analytically derived geometric sti ness matrixwhich is added to the linear-elastic sti ness matrix to create the tangent sti nessmatrix. For both examples, the final updated nodal coordinates of the structure arecompared.3.2.1Biot’s Two-Bar TrussThe Biot two-bar truss is a classical example of a linearly unstable structure which canonly be solved by considering the d

increase in sti ness, known as "stress sti ening" or "geometric sti ening". To solve structures of this nature, it is more accurate to refer the static equilibrium equations to the deformed configuration. The strains are then non-linear functions of the dis-placements, resulting in non-linear equilibrium equations. This kind of .