Add plot_warping_function(), enforce Lagrangian Multiplier constraint

[TLCFEM]

Theory

Using $\sum\omega_i=0$ over all nodes leads to $c=[1,1,1,\cdots,1]$ for the boarding matrix. This is only an approximation of the exact constraint $\int\omega~dA=0$, and the error is acceptable if the mesh density is uniform and it is only used to compute torsional constant. But if this warping function is further used in other analysis, this error is not acceptable.

The exact constraint $\int\omega~dA=0$ should be enforced, and this gives $c^T\omega_i=0$, where $c$ is the assembly of $\int{}N_{el}~dA$.

Test Code

See the difference of the solution of the following section. The symmetric section needs to have extrema of approx. identical maginitudes.

from matplotlib import pyplot as plt

from sectionproperties.analysis.section import Section
from sectionproperties.pre.library.primitive_sections import rectangular_section

rec1 = rectangular_section(b=20, d=10).shift_section(-20, -10)
rec2 = rectangular_section(b=20, d=10).shift_section(-20, 0)
rec3 = rectangular_section(b=20, d=10).shift_section(0, -10)
rec4 = rectangular_section(b=20, d=10).shift_section(0, 0)
rec = rec1 + rec2 + rec3 + rec4
rec.create_mesh(mesh_sizes=[1000, 1000, 1000, 0.1])
rectangle_section = Section(geometry=rec, time_info=True)
rectangle_section.calculate_geometric_properties()
rectangle_section.calculate_warping_properties()
rectangle_section.plot_mesh()

loc = rectangle_section.mesh["vertices"]
omega = rectangle_section.section_props.omega

print(sum(omega), min(omega), max(omega))

fig = plt.figure()
ax = fig.add_subplot(111)
ax.tricontourf(loc[:, 0], loc[:, 1], omega, cmap="viridis")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.tricontour(loc[:, 0], loc[:, 1], omega, colors="k", levels=20)
plt.show()

[TLCFEM]

You may also want to update the relevant part in the JOSS paper before submission.

[robbievanleeuwen]

Hi @TLCFEM, thanks for implementing this. This is a stricter enforcement of the Lagrangian Multiplier method that I describe.

Interestingly these changes do not change any results obtained by sectionproperties. I remembered reading about this briefly in Pilkey (not sure if you have access). On page 176:

It can be shown that the warping function ω(y, z) can be determined only up to a constant. This is proven by uniqueness studies in such books as Little (1973). However, the shear stresses of Eq. (5.20) are found to be unique since they involve derivatives of ω. This non-uniqueness property of ω(y, z) must be taken into account if ω(y, z) is to be employed in a study.

So it seems the non-uniqueness does not affect the results in sectionproperties, however may affect your use in suanPan? Either way, this is a good and necessary addition.

I also like the addition of the plot_warping_function() method - thanks! Only suggestion is to document the three arguments, level, cmap and with_lines.

I will make the changes in the JOSS paper - thanks for the heads up!

[TLCFEM]

So it seems the non-uniqueness does not affect the results in sectionproperties.

The shape of the warping function is uniquely determined by the PDE, but it's location is not, as the matrix is singular, it can be shifted up/down (like rigid body motion) while still satisfying the PDE. A unique location can be determined by applying a constraint, such as the net warping shall be zero over the section.

This does not affect torsional constant, as $\omega\dot{}P$ is not affected by the shift, as $P$ (shear force) is self-balancing.

This shift does not affect the partial derivatives as well (imaging rigid body motion does not affect strain).

However, the warping function value will be used to compute bimoment, a shifted value affects bimoment, which further affect the interaction with other sectional forces via higher-order interactions.

[TLCFEM]

Also, the section is always shifted to the barycentre to perform the warping analysis, which physically means that the barycentre is the axis of twist (no warping, as $\int{}\omega~dA=0$).

In a general setup, sometimes the section is deliberately defined and placed with an offset, shifting to the barycentre ignores this consideration. What implications it may bring is not clear to me ATM. Any thoughts?

[robbievanleeuwen]

In a general setup, sometimes the section is deliberately defined and placed with an offset, shifting to the barycentre ignores this consideration. What implications it may bring is not clear to me ATM. Any thoughts?

It's not super clear to me either. I did a quick test of your above example with the rectangle shifted away from the origin and the shifting in section.py disabled (lines 353-355). The warping function still solves ok and seems to give the same values however I cannot solve the shear functions as the error becomes too high (4% for my chosen shift).

Unfortunately I don't recall where I got the idea that I had to shift the problem to the barycentre - I can't seem to find it in Pilkey. It is interesting however, that the values of the warping function do not change after the shift (with your new code in this PR). I obviously haven't tested extensively as I've only looked at this case, but it may be something worth investigating.

[TLCFEM]

It is interesting however, that the values of the warping function do not change after the shift (with your new code in this PR).

This is expected.

The warping function still solves ok and seems to give the same values however I cannot solve the shear functions as the error becomes too high (4% for my chosen shift).

Yes. The warping function is not affected. But when computing shear stress, the shear force applied through origin now creates not only transverse shear but also torsion on the section --- if the section is not shifted to the barycentre.

The quite important fact is that, shifting the section to its barycentre invalidates different layouts that may be intentional. Clearly, a transverse shear force applied through barycentre or with a leverage would leads to different shear stresses on the section.

[robbievanleeuwen]

If this is only related to stress calculation surely this could be accounted for by varying the shear force and torsion moment accordingly. At the moment it's not clearly documented that the shear forces and torsion moment are acting at and about the shear centre, this could be improved.

[TLCFEM]

If this is only related to stress calculation surely this could be accounted for by varying the shear force and torsion moment accordingly. At the moment it's not clearly documented that the shear forces and torsion moment are acting at and about the shear centre, this could be improved.

Shear centre is another point, not the barycentre, right? And shear centre is not unique.

[robbievanleeuwen]

My understanding is that a shear force acting at the shear centre will not generate a twisting moment on the section. In fact this defines the shear centre. I'm not sure I follow "and shear centre is not unique".

The current implementation is such that when you give calculate_stress() a shear force there is no twisting, i.e. the integration of the shear stresses results solely in a force acting through the shear centre with no moment. (I think this is a separate issue to the below?)

Also I'm a little bit confused with this (apologies it's been a very long time since originally working on this and my brain is in another field now!):

However, the warping function value will be used to compute bimoment, a shifted value affects bimoment, which further affect the interaction with other sectional forces via higher-order interactions.

If we are not expecting a change in the values for the warping function ("This is expected"), what values are you expecting would change if we do not first shift the section to the centroid?

[TLCFEM]

shear centre is not unique

Because there are many definitions of shear centre.

If we are not expecting a change in the values for the warping function ("This is expected"), what values are you expecting would change if we do not first shift the section to the centroid?

I guess we are on different things.
Without the constraint, the warping function is not unique as mentioned earlier. For the constraint itself, whether to use the approximation $c=[1,1,1,\cdots,1]$ or the exact constraint $\int{}\omega~dA=0$ will affect the warping function values, this difference would affect computation of other section resultants. The constraint is actually shifting $\omega(y,z)$ along $x$ axis to meet the assigned condition.

With the proper constraint (say, for example, whichever), shifting the section around (via Geometry.shift_section(), that is shifting $\omega(y,z)$ in the $y-z$ plane) does not change warping function. But it shall change stress distribution, as forces are applied differently.

The book does assume the transverse forces are applied via shear centre (section 6.1.2, last paragraph).

[TLCFEM]

In both issues I used shift, which may be very confusing. Actually I am talking about two separate things.

[robbievanleeuwen]

Thanks for the clarification.

1. Shifting along longitudinal axis - I am happy with the new constraint that you introduced ($\int{}\omega~dA=0$) as it gives more appropriate values for the warping function.
2. Shifting in the plane of the cross-section - let me know if you are interested in investigating this further. The challenge will be to figure out how to solve for the shear functions $\Psi$ and $\Phi$ and what the resulting impact on the shear stress distribution will be. As you said before, we are currently constrained to applying the shear forces through the shear centre. The warping section properties all appear to be functions of the derivatives of $\omega$, $\Psi$ and $\Phi$, so these will not be affected.

[TLCFEM]

$\psi$ and $\phi$ involve coordinates $y$ and $z$, and they may be affected. If the section is not shifted to its barycentre, then the same shape functions used in gemetry computation can be reused. I'll find some time to look into it.