[Ayberkrk]

Follow-up to #470.

Root cause

constitutive_matrix() describes the springs correctly in its comment: each rotational spring adds 1/k to the diagonal of the element's rotational flexibility matrix. The code then inverts each stiffness term on its own (k_ii = 1 / (1/k_ii + 1/k)) instead of inverting the flexibility matrix. That is only exact when the off-diagonal coupling is zero. #470 fixed which term was touched and restored symmetry, but the term-by-term inversion remains, and for a spring softer than the element's own end stiffness the rotational block becomes indefinite. For EI = 15000, l = 4 and k = 2000 at one end its eigenvalues are about -1620 and 18385, so solve() raises a stability error. For stiffer springs it solves but gives wrong displacements.

Reproduction

A cantilever (EI = 15000, L = 4, P = 10) with a rotational spring k between the clamped support and the beam has the closed-form tip deflection P L³ / (3 EI) + P L² / k. A second case puts the spring at the end of a first element, followed by a second element of the same length and a tip load.

k	case	before #470	current master	this PR	closed form
2000	spring at node 1	0.051145	stability error	0.094222	0.094222
2000	spring at node 2	0.061356	stability error	0.193778	0.193778
20000	spring at node 1	0.062222	0.033185	0.022222	0.022222
20000	spring at node 2	0.257185	0.189630	0.121778	0.121778
200000	spring at node 1	0.015909	0.015681	0.015022	0.015022
200000	spring at node 2	0.117624	0.119613	0.114578	0.114578

I also checked these cases and examples/ex_7_rotational_spring.py against an independent beam solver that models each spring as a separate rotational degree of freedom. All agree to about 1e-15 relative error.

Fix

The rotational 2x2 block is inverted to a flexibility matrix, 1/k is added on the diagonal for each spring, and the result is inverted back. A hinge at the other end keeps its existing zeroed row and column; the spring then acts on the remaining scalar term. Axial stiffness is unchanged.

Tests

• tests/test_constitutive_matrix.py now checks the block against the inverted flexibility matrix and checks it is symmetric positive definite for soft springs. It also has two structure-level tests against the closed-form deflections above. Four of these five tests fail on current master and all pass with this change. I replaced the assertions I added in Fix node 2 rotational spring corrupting the constitutive matrix #470, because they encoded the term-by-term behaviour.
• The stored solution for example 7 in tests/test_e2e.py is updated to the verified values (node 2 deflection 0.25070634 instead of 0.23558645).

Side note: the it_* checks inside the context_* blocks of tests/test_e2e.py are not collected by pytest-describe with the default describe_prefixes, so they do not run today. That is why the outdated example 7 values went unnoticed. With -o describe_prefixes="describe context", all 45 of those checks pass on this branch except one unrelated system fixture error. I left the configuration alone to keep this PR focused, and I'm happy to open a separate PR for it.

black and pylint (10/10) pass and mypy reports the same errors as master. The full pytest suite passes apart from the six test_plotter.py tests, which fail locally in the same way on master.

[mitchell-tesch]

@Ayberkrk - Nice change. Inverting the whole block makes sense to me for finite springs.

One thing that might be worth a look while you're in here: for a hinge (k = 0), I think the other end may need the condensed 3EI/l rather than 4EI/l. At the moment the released row/column is zeroed, but the other end keeps 4EI/l. For example, test_hinge_with_spring_on_other_node expects L / (4 * EI) where I'd have expected L / (3 * EI).

A quick check: adding spring={2: 0} at a pinned support should change nothing, but here it changes the joint moments from 5.714 / 4.286 to 5.000 / 5.000.

Example and possible fix

ss = SystemElements(EI=1000, EA=1e9)
ss.add_element([[0, 0], [4, 0]])
ss.add_element([[4, 0], [8, 0]], spring={2: 0})  # try with and without the spring
ss.add_support_fixed(1); ss.add_support_roll(2, direction="x"); ss.add_support_hinged(3)
ss.moment_load(2, Ty=10); ss.solve()
print(ss.element_map[1].node_map[2].Tz, ss.element_map[2].node_map[2].Tz)

Condensing before the zeroing in constitutive_matrix():

-    if node_1_hinge or (spring is not None and 1 in spring and spring[1] == 0):
+    hinge_1 = node_1_hinge or (spring is not None and 1 in spring and spring[1] == 0)
+    hinge_2 = node_2_hinge or (spring is not None and 2 in spring and spring[2] == 0)
+    if hinge_1 and not hinge_2:
+        matrix[2][2] = 3 * EI / l
+    if hinge_2 and not hinge_1:
+        matrix[1][1] = 3 * EI / l
+    if hinge_1:
         matrix[1][1] = matrix[1][2] = matrix[2][1] = 0
-
-    if node_2_hinge or (spring is not None and 2 in spring and spring[2] == 0):
+    if hinge_2:
         matrix[1][2] = matrix[2][1] = matrix[2][2] = 0

[Ayberkrk]

@mitchell-tesch thanks, great catch! You're right: releasing one end should condense its rotation out, so the other end gets 3EI/l, and my test had baked in the 4EI/l. I've pushed your suggested change.

Your example now gives 5.714 / 4.286 with and without spring={2: 0}, which matches the 4EI/l to 3EI/l split. I added it as a test, and test_hinge_with_spring_on_other_node now expects L / (3 * EI).

One knock-on effect: the stored solution for example 4 in test_e2e.py changes. In that model node 2 ends up as a hinge, so the column is pinned at both ends and the rotation at node 3 should be M·l / (3EI) = 4.1667e-3. The old value, 3.125e-3, is what you'd get with node 2 fixed, so I updated it. With -o describe_prefixes="describe context" all 45 nested e2e checks now pass (apart from the existing system fixture error).