OpenSees Cloud

OpenSees AMI

Invertible Does Not Mean Stable

Original Post - 11 Dec 2024 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

That you can invert a stiffness matrix does not tell you everything about the numerical stability of a structural model built in OpenSees or any other finite element software. In fact, no finite element software actually forms the inverse of the stiffness matrix, but let’s go with the misleading terminology anyway.

Consider the beam shown below with a pin support at one end and no support at the other end. No loads are applied to the beam, not even self-weight.

Physically, we know this system is unstable. But let’s look at the exact and the numerical stiffness matrix of this model.

Exact Stiffness Matrix

The model has four nodal equilibrium equations but only three basic forces, i.e., the system of equilibrium equations is overdetermined. Undeterred, we can carry the unstable equilibrium through to the exact stiffness matrix:

All diagonal entries are greater than zero and the matrix is symmetric, so this may appear to be a valid, i.e., “invertible”, stiffness matrix at first glance.

But column 3 is linearly dependent on columns 1 and 4–multiply the sum of columns 1 and 4 by -1/L and you will get column 3. This is just a way of the model telling you vertical equilibrium cannot be satisfied at the unsupported end of the beam.

The matrix is rank-deficient. its determinant is zero. The matrix has a zero eigenvalue.

Now matter how you quantify it, this matrix is singular, or “not invertible”.

Numerical Stiffness Matrix

How a finite element analysis software such as OpenSees handles this unstable model can tell you a lot about the software. Here is the unstable beam model with realistic parameter values.

import openseespy.opensees as ops

# Units = kip, inch
L = 48
E = 29000
A = 20
I = 800

ops.wipe()
ops.model('basic','-ndm',2,'-ndf',3)

ops.node(1,0,0); ops.fix(1,1,1,0)
ops.node(2,L,0); #ops.fix(2,1,0,0)

ops.geomTransf('Linear',1)

ops.element('elasticBeamColumn',1,1,2,A,E,I,1)

ops.analysis('Static')
ops.analyze(1)

OpenSees will assemble the following stiffness matrix, obtained using the FullGeneral solver and the printA() command:

Due to its representation with finite precision, the stiffness matrix went from something that should be singular to something that is near-singular. Depending on how each solver in OpenSees factorizes the stiffness matrix (remember, none of the solvers actually form the matrix inverse), a solution may be reached. With no loads, that solution is zero displacements.

Things Get Hairy

Now suppose we placed a one inch long piece of hair on the unsupported end of the beam (and continued to ignore the beam’s self-weight). This piece of hair weighs approximately 2e-5 ounces, or let’s say 1e-9 kip.

P = 1e-9

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
ops.load(2,0,-P,0)

ops.analysis('Static')
ops.analyze(1)

Using the default norm unbalance convergence test with tolerance 1e-6, three orders of magnitude higher than the applied load, a solution is found for some solvers. But the displacement of the unsupported end of the beam depends on which solver is used for the analysis.

Even though the stiffness matrix was “invertible”, you can see that these displacements are meaningless. All displacements are two orders of magnitude higher than the beam length–thank you, linear geometric transformation. And some displacements are positive (up) even though the load was applied downward.

Lowering the convergence tolerance to 1e-9 will force the analysis to take more iterations in finding a solution. But each iteration just racks up the huge displacements until “convergence” of the equilibrium solution algorithm. Try the analysis yourself with the smaller convergence tolerance.