fib-international / structuralcodes

A Python library for structural engineering calculations
https://fib-international.github.io/structuralcodes/
Apache License 2.0
83 stars 22 forks source link

N-M calculation discrepancies in cross section #142

Open DanielGMorenaFhecor opened 2 months ago

DanielGMorenaFhecor commented 2 months ago

Contrasting the results from this library with proven software, the N-M diagram results are inconsistent. When the center of gravity is at the origin, the results are more accurate. It would also be beneficial to obtain the full N-M diagram, not just the lower range of M.

See Point 10 that @MestreCarlos prepared in this notebook for further information and illustrations.

talledodiego commented 1 month ago

For the first part, see also #140 and #141 that I think explains the apparent discrepancies.

For the second part: I remember we have spent some time thinking about it with @mortenengen (there is still a notebook somewhere with my comment saying # Morten: two calls? Or a single call? 🤣) Since in the API for calculate_nm_interaction_domain we included the angle theta of the neutral axis, it makes sense that only half is drawn sinche with theta equal to pi the other part is returned.

Something like this:

# Morten: two calls? Or a single call?
res_1 = section.section_analyzer.calculate_nm_interaction_domain(theta=0)
res_2 = section.section_analyzer.calculate_nm_interaction_domain(theta=0+np.pi)

fig_mn, ax_mn = plt.subplots()
ax_mn.plot(res_1.n / 1000, res_1.m_y *1e-6)
ax_mn.plot(res_2.n / 1000, res_2.m_y *1e-6)
ax_mn.set_xlabel('N [kN]')
ax_mn.set_ylabel('M [kNm]')
ax_mn.axvline(0.0, linestyle='--', color ='k')
ax_mn.axhline(0.0, linestyle='--', color ='k')
ax_mn.set_xlim([-6500, 1500])
ax_mn.xaxis.set_inverted(True)

i_begin = 0
colors = colormaps['tab20']
j = 2
for i in (3,10,40,10,40,40):
    i_end = i_begin + i + 1

    ax_mn.plot(res_1.n[i_begin:i_end] / 1000,res_1.m_y[i_begin:i_end] *1e-6,'-o',markersize=2,color=colors(j))

    ax_mn.plot(res_2.n[i_begin:i_end] / 1000,res_2.m_y[i_begin:i_end] *1e-6,'-o',markersize=2,color=colors(j+1))
    ax_mn.plot([0, res_2.n[i_begin:i_end][-1]*1e-3], [0, res_2.m_y[i_begin:i_end][-1]*1e-6], '-k')
    ax_mn.plot([0, res_1.n[i_begin:i_end][-1]*1e-3], [0, res_1.m_y[i_begin:i_end][-1]*1e-6], '-k')
    i_begin = i_end - 1
    j+=2

image

mortenengen commented 3 weeks ago

I see the point of returning the full domain in one call, but can we do the following: