Methods of truss, the TI-89 Program

Introduction

This document briefly covers the methods used to do the calculations in the TI-89 graphing calculator program truss.  It does not go into significant detail about the methods of matrix structural analysis or finite element methods, but you should be able to follow along reasonably well with a basic knowledge of linear algebra and systems of linear equations.  I will assume that you, like me, are interested in writing your own finite elements programs for your calculator or computer.

My programs are the result of my having taken a beginning level "finite elements" class, Matrix Structural Analysis, at my engineering school.  One of the texts is a deceptively small and inexpensive paperback, Matrix Analysis of Structures, by Robert E. Sennet, ISBN 1-57766-143-5.  It is an easy book and a gentle but very useful introduction into what can otherwise be a very tough and inaccessible subject.  If you are interested in studying the subject on your own, I would recommend buying this book as a starter instead of spending $100 or more for something that will be put away after a few weeks of frustration.  Once Sennet's book is well taken in and the basics are understood, any heavy-duty finite elements text should be a much easier read.  My calculator programs use no more theory than that presented in that text.

I use very few of the routines in that text because of the particular capabilities of the TI-89 graphing calculator.  The programs are written in TI-basic, which is very slow with things such as programmed routines, but gives access to some fast built-in capabilities.  I have stayed away from the more powerful calculator languages (assembly and C with TIGCC) specifically because I don't want to spend a lot of time learning techniques for a calculator that will likely be replaced with something newer in a few years.  I find TI-basic to be a pleasantly easy language with enough capability to do most things I neeed with enough speed.

This document was created with Mathematica 4.  All equations and expressions use Mathematica's conventions.  These sometimes differ from representations in text, e.g. matrices in text will use square brackets while those by themselves will not.   Vectors will be shwn in parentheses () in text.

General Concepts

The basic idea with matrix analysis (stiffness method) is that you have a stiffness matrix [k] for a structure, a vector of applied forces (f), and a vector of nodal displacements (u) that are related in the following way:

The above is the familiar statement of the "system of linear equation" problem [M](x)=(b) for the unknown vector (x) and the knwon vector (b).  One difficulty with structural analysis is that some forces, but not all, are known, and likewise for displacements.  You will recall from linear algebra that the above system may be partitioned to handle this discrepancy in the following manner:

Here, ()is the part of the applied force vector that is known, ()is the unknown part, ()represents the unknown displacements, and ()contains known displacements.  Subscrips p and s designate those items at free nodes and supports, respectively  The system may then be solved completely by separating the partitions.  ()gets solved for first, then ().

Notice that all components of [k] are known, and is known (the displacements of the supports, either specified or zero), so the latter part of the both equations will always reduce to a vector of forces. (A stiffness matrix multiplied by a displacement vector is always a force vector.)  If is a zero vector, there is no need to combine anything, since the resultant term will also be a zero vector.  By combining that resulting vector with the vector of known applied forces into a new vector (v), the problem can be brought back to the simplicity of [M](x)=(b).

can then be determined and substituted back into an appropriate equation to find .  The "normal" way to find is to use the inverse of :

As luck would have it, the matrix inversion routine in the TI-89 is excruciatingly slow (~9 min for a 25x25 matrix).  The calc does have an LU decomposition, and is is much faster (~40-50 sec for the same problem).  Solutions using LU are as fast as I know how to make them in basic, but there is a lot of room for improvement.  HP calculators have excellent structural analysis programs written by Mr. Caspar Lugtmeier.  I can't beat FEM48/49's speed, even though the HP calculators use slower hardware (~4Mhz vs. ~12Mhz).  I believe Mr. Lugtmeier uses a Cholesky method programmed in assembly.  This method takes advantage of the fact that is symmetric and has positive values on its diagonal.  If you want to take advantage of the TI-89's speed and use memory much more effectively, try doing a Cholesky routine in C and compiling with TIGCC.

The TI-89 LU is not a complete solution to the problem - it only decomposes into upper and lower triangular matrices and supplies a permutation matrix to "straighten" things out.  I used a simple forward and backward substitution routine from the book Numerical Methods for Engineers, 3rd ed. by Chapra and Canale, ISBN 0-07-010938-9, to finish the job.  If you don't want to fuss with LU, another built-in methd, simult(), is almost as fast.

Structural Stiffness Matrix

The subject of generating the structural stiffness matrix is outside the scope of this document at the lowest level.  See the Sennett book for particulars there.  Suffice it to say that the structural stiffness matrix is made by combining smaller stiffness matrices for each structural element, i.e. the things that connect the nodes.  Each type of structural element has its own general formula for the elemental stiffness matrix.  In global coordinates, the elemental matrix for a truss member is:

where A is the member cross-sectional area, E is Young's modulus, l is the member length, c is the cosine of the member's angle from its beginning (i node) to its end (j node), i.e. Δx/l which might be negative, and s is the sine of that angle.  The matrix is 4x4 because each member has two degrees of freedom (x and y translation) at each end.  The member's i node x-translation corresponds to the first row/column, i node y-translation to the second, j node x-translation to the third, and j node y-translation to the fourth row/column.  Also, each node has "coordinates" for its two degrees of freedom.  Given node number n, the degree of freedom "coordinates" are 2n-1 for x-translation, and 2n for y-translation.  The actual location in the structural matrix depends on the way in which the nodes are numbered.  Wit  hout worrying too much about the details, understand that the entries of each elemental matrix get spattered all over the structural stiffness matrix, which is a composite of all the elemental matrices.  It doesn't take too much imagination to figure out that entries corresponding to knowns and unknowns get placed in quasi-random positions, which makes partitioning a difficult task.

A relatively simple bookkeeping method that solves this problem is to "renumber" the degrees of freedom such that knowns and unknowns are guaranteed to be placed together in the structural matrix equation.  The original numbering scheme (2n-1, 2n) needs to be restored later, in order to fetch the displacement and external force for each node and display it.

The sorting method I use for this is to iterate through the nodes on order and assign their degrees of freedom to the front matrix "coordinates" if unrestrained, and to the end if restrained.  For example, a 5 node structure with 3 restraints (node 3 in the x and y directions; node 4 in the y direction) might end up with the following:

node	direction
1	x	1	no	1
1	y	2	no	2
2	x	3	no	3
2	y	4	no	4
3	x	5	yes	10
3	y	6	yes	9
4	x	7	no	5
4	y	8	yes	8
5	x	9	no	6
5	y	10	no	7

From the above chart one can see how easy it is to retrieve, say, the y displacement for node 2 in old "coordinates".  Since it is a y displacement, the formula is 2(node number) or matrix coordinate 4.  New coordinates of the right-hand column do not have such regularity, but do guarantee that known and uknown terms are each grouped with their like terms.  To derive the new coordinates, one assigns the degrees of freedom in the order encountered according to the old coordinates to either the front or back of the list, depending on restraint conditions.  The counter controlling unrestrained nodes counts forward; the counter for restrained nodes starts at the end and counts backward.  It can also be seen that the structural stiffness matrix will be 10x10 (rows x columns) and the associated partitions will be: 7x7, 7x3, 3x7, 3x3, 7x1, 3x1, 7x1, 3x1.

If the new matrix coordinates are , , , and then the entries in each element matrix can be seen as destined for the structural matrix entry given by the 's as follows:

According to this scheme, the elemental matrix entry would be placed in the structural matrix in entry , where the 's are determined by the matrix coordinates for the degrees of freedom corresponding to the i node ( , ) and the j node ( , ).  For the 5 node example above, say some member runs from node 3 to node 5.  Its 's would then be

so for this elemental matrix gets stored in and so on.

There are other ways to keep track of which structural matrix entries do what.  I chose this method because it maintains symmetry and allows easy partitioning with the facilities of the TI-89 calculator, which does not offer a Take[] or slice() function.  A C++ or Fortran programmer may want to do things differently.

When the time comes to actually insert the values in the structural matrix, I have elected to skip the creation of the elemental matrix entirely.  The values are calculated and stored directly in the structural matrix.  This is made substantially easier by the matrices symmetry and the fact that only a few values are reused in the elemental matrix, namely , , and .  Each one of those three is calculated in turn and it or its negative is stored in one half of the structural matrix.  If a value is destined for the diagonal, it is halved.  The resulting half matrix is then added to its transpose and multiplied by to arrive at the full structural matrix, ready for partitioning.   Again, this procedure is intended to minimize the use of slow TI-basic code and take advantage of faster built-in features. A 6x6 (3 node) structural matrix is generalized below.

Member Forces

Once the requisite partitioning and solving is done,  all the entries in the (u) vector are completely defined, i.e. all the displacements of the unrestrained nodes are now known.  Two distinct types of forces are then desired.  The first is the set of reactions at all the supports, which is obtained by simply premultiplying the displacement vector with the stiffness matrix, i.e. [k](u)=(f) or its partitioned variant.

Member forces require a different strategy since only part of the vector sum of forces at any one node is desired and this part is desired in local, or member, coordinates instread of glabal (x, y) coordinates.  Without going into the derivation (see Sennett's book), the following matrix equation gives the forces acting on the end of each member in its local coordinates:

It turns out that the only force that is really needed is one of the axial forces, or .  Notice that and are guaranteed to be zero due to the first 4x4 matrix.  This matrix was derived in such a way because truss members effectively carry only axial force, due to their "pinned" ends.  Because of the sign conventions used, turns out to be negative if the member is in compression, positive if in tension.  Therefore the expression is evaluated and only that part which gives is included in the program:

A list of all the member 's is made and kept for later display.  All unkowns have now been solved for, and the analysis is complete.  Please see the source code for details about data storage, plotting, and other peripheral issues.