Sister Mary Kenneth Keller’s dissertation, written in CDC FORTRAN 63, was titled “Inductive Inference on Computer Generated Patterns.”

# FORTRAN 77

# Local Stiffness Matrix for Combined Beam and Spar Element With Axial and Lateral Linear Resistance

Our objective is to develop an element with the following characteristics:

- Two-dimensional, two node element
- Euler-Bernoulli beam theory
- Axial stiffness (“spar” type element)
- “Beam on elastic foundation” characteristic
- Axial elastic resistance

Although it is doubtless possible to start with a single weak-form equation and develop the stiffness matrix, it is more convenient to develop the axial and bending local stiffness matrices separately and then to put them together with superposition.

Both spar and beam elements generally use two nodes, one at each end. For this derivation all of the constants (beam elastic modulus, moment of inertia, cross-sectional area and spring constants) will be assumed to be uniform the full length of the element. If one desires to model non-uniform beams, one can either develop an element with the desired non-uniformity or use more elements, and we see both in finite element analysis.

Let us start with the bending portion. The weak form of the equation for the fourth-order Euler-Bernoulli beam element is

where is the Young’s modulus of the material and is the moment of inertia of the beam. The variable represents the continuous spring constant along the length of the beam relative to the displacement of that beam, the “beam on elastic foundation.” The variable is a uniform load along the beam. The equations were derived using Maple with the idea of the results used on FORTRAN 77, thus the naming convention of some of the variables. An explanation of the weak form, its derivation and the significance of and can be found in a finite element text such as this.

At this point we need to select appropriate weighting functions for the equation. For beam elements we choose weighting functions to satisfy the Hermite interplation of the two primary variables at local nodes 1 and 2, to wit

where are the “displacements” for nodes

1 and 2 and are the first derivative slopes

at these nodes.

This can be expressed in matrix form as follows:

Inverting the matrix, we have

Multiplying the result, we have for the coefficients

The weighting function in its complete form is thus

This breaks down in to weighting functions for each independent variable as follows:

additionally assuming that

If we substitute these weighting functions into the weak form of the governing equations, perform the appropriate substitution, differentiation, integration and algebra, the first term results in the following stiffness matrix:

The second term (for the “elastic foundation”) yields the following stiffness matrix:

The combined local stiffness matrix for bending only is

The FORTRAN 77 code for this is as follows:

K(1,1) = 12/he**3*E1*XI1+13.E0/35.E0*he*g

K(1,2) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g

K(1,3) = -12/he**3*E1*XI1+9.E0/70.E0*he*g

K(1,4) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g

K(2,1) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g

K(2,2) = 4/he*E1*XI1+he**3*g/105

K(2,3) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g

K(2,4) = 2/he*E1*XI1-he**3*g/140

K(3,1) = -12/he**3*E1*XI1+9.E0/70.E0*he*g

K(3,2) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g

K(3,3) = 12/he**3*E1*XI1+13.E0/35.E0*he*g

K(3,4) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g

K(4,1) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g

K(4,2) = 2/he*E1*XI1-he**3*g/140

K(4,3) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g

K(4,4) = 4/he*E1*XI1+he**3*g/105

The vector for the last term is

and the FORTRAN for this is

te(1,1) = he*t/2

te(2,1) = -he**2*t/12

te(3,1) = he*t/2

te(4,1) = he**2*t/12

Now let us turn to the spar element part of the stiffness matrix. The weak form equation for this is

Here is the cross-sectional area of the beam, is a distributed axial spring constant along the spar, and is a distributed axial force along the element. To integrate from to is no different than doing so from to , only the coordinates change.

In this case we select linear weighting functions, to wit

If as before we do the substitutions and integrations, we end up with a local stiffness matrix for the spar element only as follows:

FORTRAN code for this is

K1(1,1) = E1*A/he+c*he/3

K1(1,2) = -E1*A/he+c*he/6

K1(2,1) = -E1*A/he+c*he/6

K1(2,2) = E1*A/he+c*he/3

The right hand side vector is as follows:

and the code for this is

fe1(1,1) = he*q/2

fe1(2,1) = he*q/2

Now we need to combine these. We note that there are three variables:

- x displacement (spar element only)
- y displacement (beam element only)
- rotation (beam element only)

We thus construct a element with the rows and columns in the above order, repeated twice each way for the two nodes. Doing this results in the following local stiffness matrix:

or in code

K2(1,1) = E1*A/he+c*he/3

K2(1,2) = 0

K2(1,3) = 0

K2(1,4) = -E1*A/he+c*he/6

K2(1,5) = 0

K2(1,6) = 0

K2(2,1) = 0

K2(2,2) = 12/he**3*E1*XI1+13.E0/35.E0*he*g

K2(2,3) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g

K2(2,4) = 0

K2(2,5) = -12/he**3*E1*XI1+9.E0/70.E0*he*g

K2(2,6) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g

K2(3,1) = 0

K2(3,2) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g

K2(3,3) = 4/he*E1*XI1+he**3*g/105

K2(3,4) = 0

K2(3,5) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g

K2(3,6) = 2/he*E1*XI1-he**3*g/140

K2(4,1) = -E1*A/he+c*he/6

K2(4,2) = 0

K2(4,3) = 0

K2(4,4) = E1*A/he+c*he/3

K2(4,5) = 0

K2(4,6) = 0

K2(5,1) = 0

K2(5,2) = -12/he**3*E1*XI1+9.E0/70.E0*he*g

K2(5,3) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g

K2(5,4) = 0

K2(5,5) = 12/he**3*E1*XI1+13.E0/35.E0*he*g

K2(5,6) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g

K2(6,1) = 0

K2(6,2) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g

K2(6,3) = 2/he*E1*XI1-he**3*g/140

K2(6,4) = 0

K2(6,5) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g

K2(6,6) = 4/he*E1*XI1+he**3*g/105

As long as all of the elements line up along the x-axis, we are done. But we know that this cannot always be the case. So we need to effect a rotation of the local stiffness matrix. Since each element can be either oriented differently, of different length or both, we need

to rotate the local stiffness matrix before inserting it into the global one. The rotation matrix is

where and are the angles of the elements from the x-axis. To effect a rotation, we need to first premultiply the matrix $K$ by the inverse of and then postmultiply the result by . That process is somewhat simplified by the fact that is orthogonal; thus, its inverse and transpose are identical. Going through that process, the rotated local stiffness matrix is (in code only; we have overwhelmed WordPress’ LaTex conversion capability):

Kglobal(1,1) = cosine**2*(E1*A/he+c*he/3)+sine**2*(12/he**3*E1*XI1

#+13.E0/35.E0*he*g)

Kglobal(1,2) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+

#13.E0/35.E0*he*g)*cosine

Kglobal(1,3) = -sine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)

Kglobal(1,4) = cosine**2*(-E1*A/he+c*he/6)+sine**2*(-12/he**3*E1*X

#I1+9.E0/70.E0*he*g)

Kglobal(1,5) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI

#1+9.E0/70.E0*he*g)*cosine

Kglobal(1,6) = -sine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)

Kglobal(2,1) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+

#13.E0/35.E0*he*g)*cosine

Kglobal(2,2) = sine**2*(E1*A/he+c*he/3)+cosine**2*(12/he**3*E1*XI1

#+13.E0/35.E0*he*g)

Kglobal(2,3) = cosine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)

Kglobal(2,4) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI

#1+9.E0/70.E0*he*g)*cosine

Kglobal(2,5) = sine**2*(-E1*A/he+c*he/6)+cosine**2*(-12/he**3*E1*X

#I1+9.E0/70.E0*he*g)

Kglobal(2,6) = cosine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)

Kglobal(3,1) = -sine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)

Kglobal(3,2) = cosine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)

Kglobal(3,3) = 4/he*E1*XI1+he**3*g/105

Kglobal(3,4) = -sine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)

Kglobal(3,5) = cosine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)

Kglobal(3,6) = 2/he*E1*XI1-he**3*g/140

Kglobal(4,1) = cosine**2*(-E1*A/he+c*he/6)+sine**2*(-12/he**3*E1*X

#I1+9.E0/70.E0*he*g)

Kglobal(4,2) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI

#1+9.E0/70.E0*he*g)*cosine

Kglobal(4,3) = -sine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)

Kglobal(4,4) = cosine**2*(E1*A/he+c*he/3)+sine**2*(12/he**3*E1*XI1

#+13.E0/35.E0*he*g)

Kglobal(4,5) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+

#13.E0/35.E0*he*g)*cosine

Kglobal(4,6) = -sine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)

Kglobal(5,1) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI

#1+9.E0/70.E0*he*g)*cosine

Kglobal(5,2) = sine**2*(-E1*A/he+c*he/6)+cosine**2*(-12/he**3*E1*X

#I1+9.E0/70.E0*he*g)

Kglobal(5,3) = cosine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)

Kglobal(5,4) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+

#13.E0/35.E0*he*g)*cosine

Kglobal(5,5) = sine**2*(E1*A/he+c*he/3)+cosine**2*(12/he**3*E1*XI1

#+13.E0/35.E0*he*g)

Kglobal(5,6) = cosine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)

Kglobal(6,1) = -sine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)

Kglobal(6,2) = cosine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)

Kglobal(6,3) = 2/he*E1*XI1-he**3*g/140

Kglobal(6,4) = -sine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)

Kglobal(6,5) = cosine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)

Kglobal(6,6) = 4/he*E1*XI1+he**3*g/105

The use of “sine” and “cosine” for the trigonometric functions makes it possible to compute these once for each matrix, thus speeding up computations.

One possible application of such a element is with driven piles or deep foundations in soil; the element can be used for both axial and flexural loads. The biggest problem is that the soil response is never linear, so they cannot be used in a “straightforward” fashion, but iteratively.

# Installing OpenWatcom Fortran 77

We offer the OpenWatcom Fortran 77 on our Fortran 77 Resources page. Below is a brief video on how to install it.

The installation is pretty straightforward. If you have a purely 64-bit version of Windows, it won’t install; it doesn’t go beyond 32-bit. The installation in Wine for Linux is similar and it runs fine there too.