Direct Stiffness Method of structural analysis is a matrix method of structural analysis. Being a matrix method and able to analyse a wide class of structures using the same set of steps, it is an ideal problem to solve using a software such as Scilab. Since matrix is a built-in data type in Scilab with a rich set of operators and functions that work with matrices, it allows the user to focus on understanding the procedure, direct stiffness method in this case, without first having to code and test a library of functions to represent and operate on matrices.
Here is the set of steps to analyse a Plane Frame by the Direct Stiffness Method:
Degree of freedom numbering: Based on the type of structure (and hence the maximum number of degrees of freedom per node), count the number of degrees of freedom after leaving out the restrained degrees of freedom at supports. The sequence in which the degrees of freedom are numbered is not very important and can be left to the choice of the algorithm. But remember that the band-width of the stiffness matrix depends on this numbering. To minimize the bandwidth, you can optimize the numbering sequence (Reverse Cuthill-McKee algorithm). More...
Local Stiffness Matrix of a Plane Frame Element
Rotation Matrix of a Plane Frame Element
Global Stiffness Matrix of a Plane Frame Element
Assembling Structure Stiffness Matrix by superposing Global Stiffness Matrices of individual Plane Frame Elements
Assembling Global Load Vector from Nodal Loads
Assembling Global Load Vector from Member Loads: Member loads must be converted to equivalent joint loads and then superposed on the existing Global Load Vector. Member loads may be input by the user in either global coordinates or local coordinates. Equivalent joint loads are the negative values of fixed end reactions of the plane frame element for the given member load. Since there are a number of member load types (point load, uniformly distributed load over part or entire length of the element, uniformly varying load, concentrated moment, uniformly distributed moment etc.), this step may become fairly involved. The computed equivalent joint loads must be in global coordinate system.
Solving the Stiffness Equation to compute displacements for each degree of freedom: This step will depend on how the structure stiffness matrix is stored. If it is a dense matrix (a square matrix with as many rows and columns as the number of degrees of freedom), it could be solved using Gauss Elimination or one its variations. Cholesky method is ideal since structure stiffness matrix is symmetric and positive definite. However, it is common practice to store the structure stiffness matrix in other forms, usually taking advantage of its banded nature or skyline storage taking advantage of the fact that zeros in a column above a certain row remain zero throughout. In such cases, the solver will be dependent on how the structure stiffness matrix is stored. Further, it may also depend on whether the solver stores the entire structure stiffness matrix in memory or only a part of it in memory.
Extracting member end forces from nodal displacements: Member end forces are expressed in local coordinate system and can be easily obtained by multiplying the local stiffness matrix of the member with its nodal displacements in local coordinate system. Since the nodal displacements are in global coordinate system, it is necessary to transform the displacements to local coordinate system before obtaining the member end forces.
Extracting support reactions from member end forces: Computing support reactions can be approached in one of two ways. First, by determining the list of members meeting at the support, then obtaining the member end forces for those members at the end joining the support, transforming these member end forces to global coordinate system, superposing these transformed member end forces. Second approach is the compute that sub-matrix of the structure stiffness matrix corresponding to the zero displacement boundary conditions and then multiplying this sub-matrix with the global nodal displacements.
You must understand the following limitations of the approach described above:
It can analyse only Plane Frames, because it has only one type of elements, namely, Plane Frame elements.
It computes results only for loads applied either at nodes or on members. It cannot handle problems involving settlement of supports (non-zero displacement boundary conditions), thermal forces.
You can download a tutorial with an implementation as described above from my Lulu page.