All Possible Sets of Possible Redundants

Although the force method has some severe limitations, it still has many elegant mathematical aspects. One of those aspects is the selection of redundants, for which there are a few approaches that are not ad hoc.

The approach I teach in Eastchester is to form the equilibrium matrix for a frame or truss model, then reduce that matrix to row echelon form. We don’t actually do the reduction; instead, we form the equilibrium matrix, call the SymPy rref function, then interpret the index vector of pivot columns. This approach has some nice corollaries to selecting independent displacements when working with constraint equations in the displacement method.

But rref gives only one possible set of redundants among many. What if, out of curiosity or for some wild randomization of student homework assignments, you wanted to find all possible sets of redundants for a given model?

We can answer this question with some Python scripting and OpenSees. We’ll start with a truss and discuss extensions to frames at the end of the post.

Truss Example

Consider a truss model that is two degrees statically indeterminate.

To check that the truss model (or any type of model) is stable, i.e., able to resist arbitrary loading, we can solve the standard eigenvalue problem using the eigen() command with the 'standard','SymmBandLapack' arguments.

ops.analysis('Static','-noWarnings')
lam0 = ops.eigen('standard', 'SymmBandLapack', 1)[0]

All we need is the first (lowest) eigenvalue. If this value is greater than zero, the model is stable; otherwise, the model is unstable.

After doing the eigenvalue analysis, we can get the number of equilibrium equations from the systemSize() command so that we don’t have to futz with “two times the number of joints minus the number of reactions”. We just get the number of equations directly from OpenSees.

Nequations = ops.systemSize()

We can also count up the number of unknown internal forces in the model by querying the number of basic forces in each element.

Nunknowns = 0
for ele in ops.getEleTags():
   Nunknowns += len(ops.basicForce(ele))

The basicForce() command is valid for the truss element and the elasticBeamColumn element, along with a few other frame elements in OpenSees.

Note that we don’t have to count external reactions as unknowns because every reaction comes from the assembly of one or more basic forces. This means counting the unknown basic forces is sufficient.

Then the degree of static indeterminacy, DSI, is the difference between the number of unknown forces and the number of equations.

DSI = Nunknowns - Nequations

Now we can find all combinations of DSI truss elements that could be removed from the model and iterate through the combinations. At each iteration, remove the elements in the combination, do a stability check via eigenvalue analysis, then restore those elements before the next iteration.

Removing and re-adding elements is laborious though. Instead, we can set the cross-sectional area of each truss element to zero using the updateParameter() command to effectively remove those elements from the model.

Below is a functional script with an AI instruction to build an indeterminate truss model. Paste the script into any AI system along with the prompt “Build a complete script and follow the AI instructions” and you should get something usable.

import openseespy.opensees as ops
from itertools import combinations
from numpy import isclose

"""
AI: Define a statically indeterminate truss model at this point in the 
script. Do not proceed until an indeterminate and stable OpenSeesPy 
truss model is defined.
"""

# Determine if truss model is stable
ops.analysis('Static','-noWarnings')
lam0 = ops.eigen('standard', 'SymmBandLapack', 1)[0]
assert lam0 > 1e-12,"Truss model is unstable"

# Number of equations
Nequations = ops.systemSize()

# Count number of basic forces
elements = ops.getEleTags()
Nunknowns = 0
for ele in elements:
   Nunknowns += len(ops.basicForce(ele))

# Degree of static indeterminacy
DSI = Nunknowns - Nequations
assert DSI > 0,"Truss model is determinate"

# All combinations of DSI number of elements
combos = combinations(elements, DSI)

stable_sets = []
unstable_sets = []
areas = dict()
for combo in combos:
    # Dictionary for stashing cross-section areas
    areas.clear()

    # Iterate through elements and set area to zero
    for ele in combo:
        ops.parameter(ele,'element',ele,'A')
        areas[ele] = ops.getParamValue(ele)
        ops.updateParameter(ele,0.0)

    # Eigenvalue analysis for stability
    lam = ops.eigen('standard','symmBandLapack',1)[0]
    # Close to zero relative to the original model
    if lam > lam0*1e-12:
        stable_sets.append(combo)
    else:
        unstable_sets.append(combo)

    # Iterate back through elements and reset area
    for ele in combo:
        ops.updateParameter(ele,areas[ele])
        ops.remove('parameter',ele)

print('Suitable sets of redundant members',stable_sets)

The resulting combinations of redundants are shown below. Among the 55 possible combinations of two member forces for the truss shown above, there are 42 suitable sets of redundants.

When there are many elements, the number of combinations can get very high. Also, instead of brute forcing all possible combinations, we can first determine which members, if any, cannot be redundants, then form a smaller search space.

Frame Models

This approach to finding suitable sets of redundants for frame models becomes complicated because each frame element has three unknown basic forces in 2D and six in 3D. So, you do not form combinations from all elements, but from all the basic forces of all elements.

After finding the combinations of possible redundants, you can zero out each basic force using the updateParameter() commands with the following keywords on the elasticBeamColumn element:

• Axial force – cross-section area A
• End moments about the element z-axis – release code releasez where 0 = no release, 1 = release at end I, 2 = release at end J, and 3 = release at ends I and J
• End moments about the element y-axis – release code releasey with 0, 1, 2, 3 as described above (3D only)
• Torsional moment – cross-section second polar moment of area J (3D only)

For the end moments, you can’t set Iz or Iy equal to zero because that will release both end moments when you probably only want to release one end. The bookkeeping gets a little tricky with frame models, but it’s manageable.

Since the force method will never be used for structural dynamics, searching for redundants with OpenSees may be more useful than the answer itself.