Repository navigation
Fix node 2 rotational spring corrupting the constitutive matrix - #470
Conversation
constitutive_matrix() reduces matrix[1][1] correctly for a node 1 spring, but for node 2 it was touching matrix[2][1] and matrix[1][2] (the moment/rotation coupling terms) instead of matrix[2][2]. A spring at node 2 never actually softened that node's own stiffness, and the off-diagonal terms ended up wrong and asymmetric, which a linear-elastic element's constitutive matrix should never be. The function's own comment above this block already states the correct form (c11 + 1/k1 on one diagonal, c22 + 1/k2 on the other, off-diagonal unchanged), the code just didn't match it for the node 2 case. Fix: matrix[2][2] gets the node 2 spring reduction, off-diagonal terms are left alone for both nodes. Added test_constitutive_matrix.py with a symmetry check, a check that each spring only touches its own diagonal term, and a check for the node-2-only case specifically. Verified all three fail on the old code with the exact wrong values, pass with the fix. Full suite passes, 214 tests (matplotlib needed to be installed for the plotter tests to run at all, otherwise they silently fall back to a null plotter and fail regardless of this change).
|
Gentle ping on this one. With a rotational spring on node 2, |
|
Thank you so much, @Ayberkrk ! This is a great catch, and apparently a 9-year-old bug that meant that rotational springs haven't worked correctly since very early versions of anaStruct... Thank you too for adding the tests alongside. I'll get out a release this weekend with the fix. |
|
Thanks a lot for merging this, and for the kind words! It was fun to dig into, and I had no idea the bug went back that far. I went through the history out of curiosity: it looks like it came in with the "Changed rotational spring order. Needs some checking" commit in August 2017, so I guess the checking finally happened :) While double-checking my fix against closed-form results, I found that this PR only fixed half of the problem, sorry about that. The springs were still being added term by term instead of in series with the element's flexibility. With a soft spring, If you're planning the release this weekend, it might be worth including #487 too, so the stability error from this PR doesn't reach users. Happy to adjust anything you'd like there. Thanks again for maintaining anaStruct! |
anastruct/fem/elements.py, constitutive_matrix(). The function's own comment states the correct flexibility matrix with springs on both sides:
only the diagonal terms change, off-diagonal stays put. The code matches that for node 1 (matrix[1][1]) but for node 2 it modifies matrix[2][1] and matrix[1][2] instead of matrix[2][2]:
So a node 2 spring never softens matrix[2][2] at all, and instead corrupts the moment/rotation coupling terms, leaving the matrix asymmetric, which is not physically possible for a linear-elastic element.
Reproduced directly against a hand-derived matrix (EI=15000, l=4, spring={1:1000, 2:2000}): expected [[937.5, -7500], [-7500, 1764.7]] for the [1:2, 1:2] block, old code returns [[937.5, 2727.27], [731.71, 15000]], asymmetric and node 2's spring stiffness never shows up on the diagonal at all.
Also reproduced the practical effect: a cantilever fixed at node 1 with only a rotational spring at node 2 under a moment load gets flagged as an unstable structure (StabilityError), a false positive caused by the corrupted, singular constitutive matrix.
Fix is switching matrix[2][1]/matrix[1][2] to matrix[2][2] for the node 2 case, matching the node 1 case and the function's own comment.
Added tests/test_constitutive_matrix.py: a symmetry check, a check that each node's spring only touches its own diagonal term against the exact hand-derived values, and a case with only a node 2 spring. Verified all three fail on the old code with the precise wrong values shown above, and pass with the fix.
This is very likely the root cause behind #144 (element rotation wrong with springs at element ends), that report used spring={1: k1, 2: k2} and saw wrong rotation, then suspected geometric nonlinearity, but the linear constitutive matrix itself was already wrong before any nonlinear iteration starts. Left that issue alone since it's not the same fix and didn't want to conflate the two.
Full test suite passes, 214 tests (installed matplotlib in the test venv so the plotter tests actually run instead of silently falling back to a null plotter, this is unrelated to the fix itself).