Errors with cantilever beam

[sitadrost]

Hi Ritchie,

I just installed anaStruct, and so far I like it a lot, so thanks for sharing! However, when playing around with a simple cantilever beam case, I'm running into a number of errors that I don't quite understand, so I was hoping you could give some explanation.

So here's my toy problem:

from anastruct import SystemElements

Ne = 10     # number of elements
F_y = 10    # point load

ss = SystemElements(EI=5e3, EA=1e5)
ss.add_multiple_elements([[0, 0], [10, 0]], Ne)
# clamp first node
ss.add_support_fixed(1)
# point load on last node
ss.point_load(ss.id_last_node, Fy = -F_y)

# solve, linear
ss.solve()

# solve with geometrical non-linearity
ss.solve(geometrical_non_linear = True)

With the values of Ne and F_y used above, everything goes fine with the linear solver (I'll get to the geometrical_non_linear later on). However, if I increase Ne to, say, 100, ss.solve() gives me this error:

Traceback (most recent call last):

  File "<ipython-input-16-a7e1225dab88>", line 1, in <module>
    ss.solve()

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system.py", line 532, in solve
    "The eigenvalues of the stiffness matrix are non zero, "

FEMException: ('StabilityError', 'The eigenvalues of the stiffness matrix are non zero, which indicates a instable structure. Check your support conditions')

Solving the same problem (i.e. Ne = 10, Fy = 10) with geometrical_non_linear set to True results in the following error (independent of whether I add extra discretisation, e.g. discretize_kwargs=dict(n=20)):

Traceback (most recent call last):

  File "<ipython-input-47-3c8b776ea164>", line 1, in <module>
    ss.solve(geometrical_non_linear = True)

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system.py", line 554, in solve
    self, verbosity, discretize_kwargs=discretize_kwargs

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system_components/solver.py", line 132, in geometrically_non_linear
    buckling_factor = det_linear_buckling(buckling_system)

  File "/Applications/miniconda3/lib/python3.6/site-packages/anastruct/fem/system_components/solver.py", line 104, in det_linear_buckling
    kg = system.reduced_system_matrix - k0

ValueError: operands could not be broadcast together with shapes (20,20) (30,30) 

Can you explain what's causing these errors?

Many thanks, Sita

P.S. I just noticed my code might be a bit unclear/misleading, so to clarify: I'm using either ss.solve() or ss.solve(geometrical_non_linear = True), not both in one go.

[ritchie46]

Hmm.. this is a strange bug indeed. I noticed that at 43/ 44 elements the eigenvalues of the reduced system matrix suddenly become negative. I don't really understand how / why this happens? I would need some research or help on that part.

[sitadrost]

Thanks for getting back, sorry, I nearly forgot about this. I'd be happy to help if I can, what exactly do you need?

[ritchie46]

Need some help to understand why the eigenvalues become negative. Is this a bug in the code, or is this a property of a system matrix with too many nodes?

[sitadrost]

Ah, I see. I'll try to have a look when I have time

[smith120bh]

As one possible path to pursue, I've spotted this error before when elements are too small as an absolute length number. For example, note in the above script, if you replace: ss.add_multiple_elements([[0, 0], [10, 0]], Ne) with ss.add_multiple_elements([[0, 0], [1000, 0]], Ne) and run it again with 100 elements, it will run perfectly happily.

My wrapper script for running anaStruct limits the minimum size of elements, and we stopped encountering this. I'd suspect that there's a numerical instability, perhaps related to floating point roundoff error, occurring somewhere in the code? I have no idea where though - I can try to have a poke at some point as well.

[ritchie46]

Has it to do with the ratio of the element size and the force magnitudes? I can image that thos influence the final eigenvalues of the stiffness matrix. When it is a linear calculation, we can probably do a calculation with a unit load and upscale the results afterwards.

[smith120bh]

Ah - got it! It's specifically a bug with the add_multiple_elements() function, and it is indeed hitting a floating point roundoff issue. That is that, with small float values, n*(length/n) == length is not necessarily true.

Running sitadrost's original example with Ne=100, there were actually 101 elements created in which the following occurred:

[el.l for el in ss.element_map.values()]
# [0.10000000149011612, ..., 0.10000000149011612, 1.9073486328125e-06]

That relatively much smaller element at the end there was causing the global stiffness matrix's condition number to explode (into the 1e17+ range), and the eigenvalues to go crazy.

@ritchie46 : I'll PR a fix to this momentarily. It also might be worth adding a warning if the condition number of the reduced system matrix becomes too large and/or if the relative lengths of elements varies too much.