1.0 Introduction
In this era of information technology, civil engineers rely heavily on software to perform their design tasks. Unfortunately, most commercial structural analysis packages are closed-source, which means that the operations that the program performs cannot be inspected by the user.
Moreover, such software packages are invariably very pricey, and, hence, are generally not affordable for students and smaller engineering firms.
The objective of this design project was to design a structural analysis program that would be free of charge and available to all. This computer program was to be open source and well commented, so that its users could comprehend the operations performed in the analysis of a given structure.
To accomplish these objectives, the generalized stiffness method of structural analysis was implemented into a computer algorithm. This algorithm, called “TrussT Structural Analysis”, is a collection of visual basic modules embedded in a Microsoft Excel document using Visual Basic for Applications (VBA). This design report outlines the theory behind TrussT Structural
Analysis, as well as the methods by which that theory was implemented into computer algorithms. The first two sections of this report present the theory of the generalized stiffness method of structural analysis and its implementation into a computer algorithm. The following sections present the procedures by which the stiffness method was modified to incorporate the analysis of structure with special characteristics such as member applied loads, member releases or support settlements. A computer implementation of the Euler method of analysis is described to account for the geometric non-linearity of structures. Finally, algorithms that can generate member force diagrams and moment diagrams are presented.
2.0 Generalized Stiffness Method
The term stiffness refers to a body’s ability to resist imposed displacements by generating internal forces. Conversely, the term flexibility refers to a body’s ability to deflect when subjected to
6
applied forces. These two terms are analogous to two methods of structural analysis: the flexibility method, which generates compatibility equations to solve for forces, and the stiffness method, which generates equilibrium equations to solve for displacements. The preliminary report done in Phase I of this project outlined the reasons why the stiffness method is more suitable to computer implementation than its counterpart. The computer algorithm presented in this report was designed based on the stiffness method of structural analysis.
The stiffness method is based on the fundamental equilibrium equation shown in Equation 1. (1)
This equation equates force (P) to the product of stiffness (KK) and displacement (d). In the simplest terms, this equation represents the equilibrium between external forces and internal forces. In terms of a structure, the applied joint forces (P) and joint displacements (d) are vectors, whereas the structure stiffness is expressed as a matrix of n-by-n terms, n being the degrees-of-freedom of that structure. Generally, both the applied forces (P) and the structure’s stiffness (KK) are known, which leads to a system of n equations with n unknown joint displacements (d).
3.0 Computer Implementation of the Stiffness Method
The analysis engine driving TrussT Structural Analysis is a computer implementation of the stiffness method. This section describes the basic approach that was taken to implement the stiffness method in a computer program. That implementation as utilized in TrussT Structural
Analysis is illustrated in Figure 1.
7
Obtain Input Parameters
Organize Inputs in Arrays
(Record Keeping)
Calculate Member Stiffnesses (s)
Assemble Global Stiffness Matrix (K)
Solve for Global Displacements (d)
Transform Displacements to
Local Coordinates (u)
Calculate Member Forces (Q)
Transform Member Forces to
Global Coordinates (F)
Assemble Reaction Forces (R)
User Inputs
Joint
Displacements
Reaction
Forces
Figure 1: Analysis Procedure Flow Chart
3.1 Input
TrussT Structural Analysis requires information about the structure to be inputted by a user in three parts: joint inputs, member inputs, and load inputs. Appendix A of this report presents a user’s manual, which outlines how to use the program and input the necessary variables.
The variables required by the program to analyse a structure are described in following subsections. Consistent units must be used when inputting information in the user interface.
Although the user interface suggests the metric units of mm, kN and GPa, any set of consistent units will produce correct results (e.g. in, kips and ksi).
3.1.1 Joint Inputs
For each joint, the program requires:
 Joint Coordinates - numerical values for the three dimensional coordinates (x, y, z).
8
 Joint Boundary Conditions – indicating whether translations and rotations are free or restrained along the x, y, and z axes. Alternatively, the user can specify whether a joint is fixed, pinned or free.
3.1.2 Member Inputs
For each member, the program requires:
 Member Joints - the two joints the member spans.
 Member Orientation (ψ) – defined as either (1) for general members or (2) for vertical members. (1) The angle, measured clockwise positive, when looking in the negative x direction, through which the local xyz coordinate system must be rotated around its x-axis so that the xy plane is vertical with the y-axis pointing upwards.
(2) The angle, measured clockwise positive, when looking in the negative x direction, through which the local xyz coordinate system must be rotated around its x-axis so that the z axis is parallel to, and points to global z.
 Member Section Properties – modulus of elasticity E, area A, shear modulus G, moment of inertia about the strong axis Iz, moment of inertia about the weak axis Iy, and the torsion constant J.
3.1.3 Load Inputs
For each load the program requires the:
 Load Location – the joint on which an external load is applied.
 Load Magnitude – numerical value of the load (force and moment) in each axis.
3.2 Record Keeping
In order to organize adequately all the input and output information in a set of matrices, the program assigns record-keeping indices to every joint translation and rotation. These indices, later referred to as codes, represent the position of a certain joint displacement or reaction in the
9
global vectors and matrices presented in Equation 1. In three dimensions, each joint has six possible movements (translations and rotations about each of the x, y, and z axes). The computer program first assigns codes to the joint movements that are free, that is, joint movements that are not boundary conditions. The total number of such movements is referred to as the degrees-offreedom of the structure. Once all the degrees-of-freedom have been assigned codes, the computer program then assigns codes to the joint movements with boundary conditions.
3.3 Member Stiffness
Based on the inputted member properties, the program generates a 12 x 12 member stiffness matrix (s ss) and transformation matrix (T TT) for each member. Equation 2 illustrates the local stiffness matrix for a space frame element based on the member properties, while Equation 3 shows the transformation matrix required to transform the member’s stiffness matrix into the structure’s coordinate system.

(2)

(3)
…where (rrr) is given by Equation 4 for general members and by Equation 5 for vertical members.
10

(4)

(5)
The member stiffness matrix (s ss) for each member is converted from the Local Coordinate
Systems (LCS) into the Global Coordinate System (GCS) through Equation 6, which yields the member stiffness matrix in the GCS (k kk). k (6)
3.4 Global Stiffness
A structure’s global stiffness (KK) consists of a n-by-n matrix of coefficients, where n represents the degrees-of-freedom of the structure. To assemble (KK), the program uses the codes assigned in the record keeping phase to assemble the individual member stiffness matrices to the global stiffness matrix.
3.5 Determining Displacements
The joint displacements at the free degrees-of-freedom (d) are computed using the joint equilibrium equation, Equation 1, which relates these displacements to the known applied joint loads vector (P) and to the known structure’s global stiffness (KK). A Gauss-Jordan solver is used to solve for these member displacements.
11
3.6 Converting Displacements to Local Coordinates
Once the displacements in the GCS (d) are determined, they are stored as the member end displacements (v) using the record keeping codes. These global displacements are then converted to the LCS through the member transformation matrices, as shown in Equation 7. (7)
3.7 Calculating Member Forces
Once the joint local member displacements (u) are known, the algorithm calculates the member forces (Q) based on the member’s local stiffness matrix (s ss) with Equation 8. Equation 8 is the member-space equivalent of the basic equilibrium equation, Equation 1. (8)
3.8 Transformation of Member Forces
The member forces in LCS (Q) are transformed back to the GCS using Equation 9, at which point they are assembled based on the structure’s configuration (using the record keeping codes) to obtain the global reactions. (9)
3.9 Computer Algorithm
The following pseudo-code outlines the stiffness method presented in this section to solve for a framed structure under joint loads. The nomenclature of the variables and the procedures used are described following the code.
'Assembly of Applied Loads Vector (P)
Loop i = 1 To numNodes
M1() = GetCodes(Node(i))
M2() = GetForces(Node(i))
Loop j = 1 To 6
P(M1(j)) = M2(j)
End Loop
End Loop
12
‘Assembly of Member Stiffness Matrix (S) and Transformation Matrix (T); Performed in Member Class
Module for each member
S_local() = MemberStiffness(E, A, L, Iy, Iz, J)
T() = Transformation(x1, x2, y1, y2, z1, z2)
S_global() = Transpose(T()) * S_local() * T()
'Assembly of Structure’s Stiffness Matrix (K)
Loop i = 1 To numMem
M1() = GetCodes(Member(i))
M2() = GlobalStiffness(Member(i))
Loop j = 1 To 12
If M1(j) minDOF Then
P(M1(j)) = P(M1(j)) + M2(j)
End If
End Loop
End Loop
Description of Variables and Procedures
M1() - Temporary matrix
M2() - Temporary matrix numNodes – Number of nodes in the system numMem – Number of members in the system
S_local() - Local member stiffness
S_global() - Global member stiffness matrix
K() - Structure’s overall stiffness matrix
Kff() - Structure’s partitioned stiffness matrix for free nodes
E, A, Iy, Iz, J – Member properties specified by the user
L – Member length
T() – Member transformation matrix x1, x2, y1, y2, z1, z2 – Member coordinates minDOF – Number of degrees-of-freedom in the system
P() - Joint forces vector
Q() - Member forces vector in local coordinates
F() - Member forces vector in global coordinates u() - Nodal displacements vector in local coordinates
D() - Nodal displacements vector in global coordinates
Node(i) - Class module for node i
13
Member(i) - Class module for member i
Getcodes - Invokes code number stored in the class module.
GetForces - Invokes user-specified load stored in the class module.
MemberStiffness - Calculates the member stiffness matrix using Equation 2.
Transformation - Calculates the member transformation matrix using Equation 3.
MemberStiffness - Invokes the member stiffness in the global coordinate system stored in the class module.
MemberForce - Invokes the member forces stored in the class module.
Transpose - Performs transposed matrix operation.
MatInv - Performs inverse matrix operation.
4.0 Equivalent Joint Loads
As the basic stiffness method considers only loads applied at the joints, the concept of Equivalent
Joint Loads (EJL) can be used to expand the scope of the stiffness method to include member loads. To determine the EJL of a structure, the EJL of each member must be determined and assembled using the record-keeping codes. Kassimali [1] and others have shown that these EJLs are equivalent to the negative of the fixed-end forces generated by a set of loads on a fully fixed member. Fixed-end forces for a member subjected to any set of loads can be found based on the already derived formulae found in many structural analysis textbooks, such as Kassimali [1]. The equations used in TrussT Structural Analysis to solve for the fixed-end forces are included in
Appendix B. TrussT Structural Analysis has the capability of converting member point loads and linearly variable distributed loads to equivalent joint loads. No formula could be found in available literature that described the fixed-end forces caused by a linearly variable axial load.
Therefore, a formula was derived for linearly variable axial loads based on the axial point load formula presented by Kassimali [1]. The derivation of this can be found in Appendix B.
The fixed-end forces are assembled from local coordinates (Qf) to global coordinates (Pf) using the same coding mechanism used to assemble the structure’s stiffness matrix.
The equivalent joint loads are incorporated in the procedure described in section 3 by subtracting the fixed-end forces (adding the EJL) from the applied joint loads. This operation is integrated into the analysis by modifying Equations 1 and 8 to Equations 10 and 11. (10) (11)
14
The following pseudo-code illustrates the conversion of member loads to equivalent joint loads in
TrussT Structural Analysis.
If bDist = False Then
R() = PointLoad(L, Pos, sForces(3))
EJL(3) = EJL(3) + R(1)
EJL(5) = EJL(5) + R(2)
EJL(9) = EJL(9) + R(3)
EJL(11) = EJL(11) + R(4)
R() = PointLoad(L, Pos, sForces(2))
EJL(2) = EJL(2) + R(1)
EJL(6) = EJL(6) + R(2)
EJL(8) = EJL(8) + R(3)
EJL(12) = EJL(12) + R(4)
R() = AxialLoad(L, Pos, sForces(1))
EJL(1) = EJL(1) + R(1)
EJL(7) = EJL(7) + R(2)
R() = Moment(L, Pos, sForces(6))
EJL(2) = EJL(2) + R(1)
EJL(6) = EJL(6) + R(2)
EJL(8) = EJL(8) + R(3)
EJL(12) = EJL(12) + R(4)
R() = Moment(L, Pos, sForces(5))
EJL(3) = EJL(3) + R(1)
EJL(5) = EJL(5) + R(2)
EJL(9) = EJL(9) + R(3)
EJL(11) = EJL(11) + R(4)
R() = Torque(L, Pos, sForces(4))
EJL(4) = EJL(4) + R(1)
EJL(10) = EJL(10) + R(2)
Else
R() = Distributed(L, Pos1,Pos2, sForces(5), sForces(6))
EJL(3) = EJL(3) + R(1)
EJL(5) = EJL(5) + R(2)
EJL(9) = EJL(9) + R(3)
EJL(11) = EJL(11) + R(4)
R() = Distributed(L, Pos1,Pos2, sForces(3), sForces(4))
EJL(2) = EJL(2) + R(1)
EJL(6) = EJL(6) + R(2)
EJL(8) = EJL(8) + R(3)
EJL(12) = EJL(12) + R(4)
R() = AxialDistributed(L, Pos1, Pos2, sForces(1), sForces(2))
EJL(1) = EJL(1) + R(1)
EJL(7) = EJL(7) + R(2)
End If
Description of Variables and Procedures bDist - Identifier for distributed loads sForces() - Member load vector
( sForces(1) = Fx; sForces(4) = Mx for point loads )
( sForces(1) = Fx1; sForces(2) = Fx2 for distributed loads )
L - Member length
Pos - Location of point load
Pos1 – Location of near-end distributed load
Pos2 – Location of far-end distributed load
R() – Temporary matrix to store computed equivalent joint load for the member load being considered.
EJL() – Equivalent joint load vector
15
PointLoad - Converts point member load to equivalent joint load.
AxialLoad - Converts point axial load to equivalent joint load.
Moment - Converts applied moment to equivalent joint load.
Torque - Converts torque to equivalent joint load.
Distributed - Converts distributed load to equivalent joint load.
AxialDistributed - Converts distributed axial load to equivalent joint load.
*Equations used for the above sub-function are included in Appendix B
5.0 Member Releases
TrussT Structural Analysis has the ability to account for member releases in a structure. A member is said to be released at a certain location when its ability to develop an internal force at that location has been removed. Physical examples of member releases are hinges (rotational releases) and slotted connections (translational releases). Figure 2 depicts the possible locations of member releases on a space frame member.
Figure 2: Possible Locations of Member Releases
Releases at intermediate locations along the member can also be modeled using the approach outlined below by utilizing two smaller members and releasing the joining ends.
To account for a member release, one must modify the basic equilibrium equations to ensure that the force at the location being released is zero. That is, if one is to release the translation at #3, then the internal force Q3 in Equation 11 must be set to zero. Setting Q3 equal to zero implies that the right-hand side of the Q3 equilibrium equation is also zero, as shown in Equation 12. (12)
Equation 12 shows that there is inter-dependence between the set of member displacements.
Equation 12 can therefore be re-written as Equation 13. … / (13)
16
The stiffness matrix needs to be modified to reflect the physical reality of the member release.
Given that Q3 = 0, any set of displacements (u) will cause Q3 to remain 0. Therefore, the third row of the member stiffness matrix should be a row of zeros. As Equation 14 shows, the equilibrium conditions now consist of 11 equations containing 12 unknown displacements.

,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,

(14)
Given that there is one too many unknowns for this system of equations to be solvable, the displacement u3 must be distributed to the other displacements according to Equation 13.
Additionally, as Equation 13 shows, the released displacement u3 also depends on Qf3, the member local fixed-end force at location #3. Therefore as u3 is being substituted into the other rows of the stiffness matrix, the fixed-end forces corresponding to this row must also be modified according to Equation 12.
The final element stiffness matrix will be an 11 x 11 matrix. This released matrix can also be thought of as a 12 x 12 matrix with a row and column of zeros. The row of zeros signifies that no force can be generated at that location regardless of the member displacements, and the column of zeros is necessary to have a solvable system of equations. The effect of that displacement has been expressed as a sum of other displacements and fixed-end forces, and the stiffness and member fixed-end forces coefficients have been modified accordingly.
The following procedure is analogous to how TrussT Structural Analysis implements the member release procedure.
1. Determine location to be released (r)
2. Express ur as the sum of other displacements:
17

3. For every ij element of the stiffness matrix, subtract the following stiffness according to
Equation 12:

4. Similarly, for the each i fixed-end force, subtract the following force according to
Equation 12:

5. Set all ri elements of the stiffness matrix equal to zero: 0
6. Set all ir elements of the stiffness matrix equal to zero: 0
The complete algorithm for this process can be found in the cMember class module of the program code, under the sub-procedures ReleaseMatrix and ReleasedQF. TrussT Structural
Analysis uses the above procedure to transform the space frame stiffness matrix into various released matrices. Instead of having different algorithms to solve for the displacements of trusses, beams or other simplified structures, TrussT uses member releases to transform the space frame stiffness matrix to the desired simplified matrix. As Equation 15 shows, a space frame stiffness matrix released 10 times at r = 2, 3, 4, 5, 6, 8, 9, 10, 11, and 12 is equivalent to a 2 x 2 axial member stiffness matrix.
18

(15)
6.0 Support Settlement
The alternative formulation of the stiffness method, as presented by Kassimali [1], is particularly well suited to the incorporation of a support settlement into the analysis of a structure. This method involves assembling a global stiffness matrix (KK* **) of size m-by-m, where m represents not only the quantity of degrees-of-freedom, but the total quantity of degrees in the structure.
This method also involves assembling an m sized applied load vector (P*), fixed-end forces vector (Pf*) and global displacement vector (d*). The resulting equilibrium is shown in Equation
16.

(16)
Equation 16 is particular in the way that it has unknowns on both sides of the equation. To solve such a system, Kassimali’s proposed partitioning the KK* ** matrix into four sub-matrices and partitioning the P*, Pf* and d vectors into two sub-vectors each as follows:

(17)

(18)

(19)

(20)
19
…where dR is the vector of known support settlement and PR is the vector of unknown applied loads (support reactions).
The KK and KKF FFR RR sub-matrices are used to compute the joint displacements (d) as per Equation 21, whereas the KKR RRF FF and KKR RRR RR sub-matrices are used to compute the unknown reactions forces (PR) as per Equation 22. (21) (22)
Although this formulation requires larger matrices, it allows for the unknown global reactions to be solved directly, without using the assembly process described by Equation 9.
One drawback of this method of analysis arises when conducting second-order analyses (as discussed in Section 7). Since second-order analyses required the computation of the member forces, it was found not to be computationally advantageous to compute reaction forces as per
Equation 22. Second-order analyses require that the systems of equations be solved numerous times, therefore the reduction of the number of operations within the loop of iterations will drastically reduce the overall computation time.
For that reason, a hybrid between Kassimali’s alternative formulation and the original formulation was adopted by TrussT Structural Analysis, whereby the global displacements are computed by Equation 21, and the support reactions are computed as before, by assembling the member forces (F). Therefore, Equation 22 is not utilized in the TrussT Structural Analysis algorithm. It should be noted that a structure with no support settlement has a (dR) displacement vector of 0; therefore, Equation 21 reverts to its original form in Equation 10.
20
7.0 Second-Order Analysis
In first-order analysis, it is assumed that the load-displacement relationships of structural elements are linear, such that the linearly-assumed displacements under an increasing applied load (P) will deviate from the actual displacement, as shown in Figure 3.
Figure 3: Equilibrium Paths of Linear and Nonlinear Analyses
As Figure 3 illustrates, the discrepancy between the assumed displacements and the actual displacements grows enormously as the applied load increases.
Second-order analysis involves continuously or intermittently changing the stiffness of a structure so that it responds to applied loads based on its current configuration, rather than based on its initial configuration. This procedure more accurately models the displacement behaviour of an element. TrussT Structural Analysis has a function that incorporates geometric nonlinearity into the analysis. The algorithm is based on the Euler method, also known as the incremental single-step method. The Euler method is the most elementary second-order analysis method; however, it is suitable for solving a system with moderate nonlinearity, which applies to most structural engineering applications. Figure 4 illustrates an element equilibrium path using second-order analysis based on the Euler method alongside the actual equilibrium path.
21
Figure 4: Euler method second-order iteration
As Figure 4 illustrates, the actual equilibrium path of an element is more accurately modelled with the Euler method than with first-order analysis. Furthermore, as the number of iterations increases, the linear segments outlining the equilibrium path will become shorter and the accuracy of the model will increase.
In the program, once the user chooses to perform a second-order analysis, the user may also specify the number of iterations (n). The prescribed loads are then discretized into n incremental loads (ΔP). The program performs n linear analyses that are analogous to the first-order analysis, described in sections 2 and 3, to calculate the incremental displacements (Δu) and the member forces (ΔQ). The second-order analysis is different from the first-order analysis in the sense that it uses the tangent stiffness matrix (Kttt), which is the sum of the elastic stiffness matrix (K) and the geometric stiffness matrix (Kggg). The member geometric stiffness matrix is given in Equation
23.
22
(23)
As shown in Equation 23, the member geometric stiffness matrix is a function of the member axial force (Fx2), hence, the geometric stiffness (and the tangent stiffness) will change as the nodal displacements change along the iteration. Once the iteration is complete, the calculated incremental displacements and member forces