concrete_section.moment_interaction_diagram

[buddyd16]

concrete-properties/concreteproperties/concrete_section.py

Lines 940 to 987 in 503a3f0

	def moment_interaction_diagram(
	self,
	theta: float = 0,
	control_points: List[Tuple[str, float]] = [
	("kappa0", 0.0),
	("D", 1.0),
	("fy", 1.0),
	("N", 0.0),
	("d_n", 1e-6),
	],
	labels: List[Union[str, None]] = [None],
	n_points: Union[int, List[int]] = [4, 12, 12, 4],
	max_comp: Optional[float] = None,
	max_comp_labels: List[Union[str, None]] = [None, None],
	) -> res.MomentInteractionResults:
	r"""Generates a moment interaction diagram given a neutral axis angle `theta`
	and `n_points` calculation points between the decompression case and the pure
	bending case.
	:param theta: Angle (in radians) the neutral axis makes with the horizontal axis
	(:math:`-\pi \leq \theta \leq \pi`)
	:param control_points: List of control points over which to generate the
	interaction diagram. Each entry in ``control_points`` is a ``Tuple`` with
	the first item the type of control point and the second item defining the
	location of the control point. Acceptable types of control points are
	``"D"`` (ratio of neutral axis depth to section depth), ``"d_n"`` (neutral
	axis depth), ``"fy"`` (yield ratio of the most extreme tensile bar), ``"N"``
	(axial force) and ``"kappa"`` (zero curvature compression - must be at start
	of list, second value in tuple is not used). Control points must be defined
	in an order which results in a decreasing neutral axis depth (decreasing
	axial force). The default control points define an interaction diagram from
	the decompression point to the pure bending point.
	:param labels: List of labels to apply to the ``control_points`` for plotting
	purposes, length must be the same as the length of ``control_points``. If a
	single value is provided, will apply this label to all control points.
	:param n_points: Number of neutral axis depths to compute between each control
	point. Length must be one less than the length of ``control_points``. If an
	integer is provided this will be used between all control points.
	:param max_comp: If provided, limits the maximum compressive force in the moment
	interaction diagram to ``max_comp``
	:param max_comp_labels: Labels to apply to the ``max_comp`` intersection points,
	first value is at zero moment, second value is at the intersection with the
	interaction diagram
	:raises ValueError: If ``control_points``, ``labels`` or ``n_points`` is invalid
	:return: Moment interaction results object
	"""

The docstring for the method notes it takes the neutral axis rotation as the entry parameter but does not appear to vary this value for each interaction point, need to be careful here. Holding the neutral axis rotation constant and solving for the different neutral axis depths to achieve the plot points does not actually give you the proper interaction diagram. If you plot the resulting points in 3D you'll see that they oscillate away from a constant Mx, My plane meaning the returned interaction points would not be in equilibrium with applied loading.

[robbievanleeuwen]

Hi @buddyd16, I see your point here, however the function does what it sets out to do. I don't think there's anything in the docs about keeping the moment ratio constant.

You are correct in noting that for bending about a non-principal axis the ratio between m_x and m_y will vary. The anlysis reflects this and these results are saved in the MomentInteractionResults object, and the user can easily view these results. The plot_diagram() method plots either the m_x, m_y or resultant component of the bending moment (m_x by default).

From a design perspective I can see the value in requiring the ratio between moments is constant rather than the bending axis angle and I think this would be a useful feature to add! Do you have a reference on how this can be done? Seems like a much more difficult problem with two unknowns to be solved (as discussed on engtips forum)?

[buddyd16]

I misunderstood your plots, plotting only one component of the biaxial bending is unfortunate, that would provide a false design utilization ratio.

for design purposes you really need to plot the interaction curves on a constant mx, my plane consistent with the applied loads.

I have two spreadsheets for the problem in the open source group, the approach I took was to nest two root finding methods outer method for the neutral axis rotation and inner method for the neutral axis depth. A constrained optimization approach may be a better fit though minimizing sum of squares for the delta in axial load and deltas in i and j unit vectors for the moments with neutral axis depth constrained to 0 and a depth where the entire cross section is at concrete ultimate strain and the neutral axis angle constrained to be between -pi and pi.

[connorferster]

Maybe I am missing something here but it seems that the kind of plots @buddyd16 is looking for can be generated by the BiaxialBendingResults object where the results are plotted as mx vs my on a constant axial plane. Is this incorrect?

[buddyd16]

@connorferster @robbievanleeuwen hopefully this helps illustrate what I perceive to be an issue:
Cross Section:

Loading:

Mx vs My Slice at applied P:

P v M at applied load angle atan2(mx,my):

3D Interaction Surface showing the above diagrams and what it appears ConcreteProperties is doing:

Top down view showing how a constant neutral angle plot strays from the true interaction curve:

SPCol and other commercial apps are doing the P v M at applied load angle atan2(mx,my) line.

[robbievanleeuwen]

@buddyd16 thanks for the discussion and example, think this has finally clicked for me 😆. The current implementation works for bending about principal axes but as you describe does not give useful (read for design purposes - where m_x / m_y is known) interaction diagrams when there is biaxial bending.

This is showing how out of touch I have become with structural design after a few years out of the game! I've put my designer (rather than academic) hat back on and here's a bit of a brain dump - feel free to comment on this:

How users would currently interact with concreteproperties:

• Uniaxial bending about a principal plane ➡️ set theta equal to 0 or pi / 2 as required and generate interaction diagram - utilisation ratio can be obtained readily from this diagram as there is bending about one axis only ✅
• Biaxial bending present, user wants to generate an interaction diagram for a ratio of m_x to m_y ➡️ there is facility to do this in concreteproperties and user does not know which theta to set (as you describe there is no set theta for a given moment ratio!) - user is confused and cannot generate a moment interaction diagram for this case ❌ (worst case is theta=0 is a non-principal bending axis and user takes these results as gospel ❌❌❌)
• User has a known axial force with biaxial bending ➡️ create biaxial bending diagram for known axial force ✅ (however a utilisation ratio cannot be obtained unless the user has zero axial force as axial force cannot be scaled on a constant axial force plot if it is non-zero ❌)
• Case as described by @buddyd16: User has biaxial bending with an unknown axial force and wants to find the utilisation ratio ➡️ only current option is to painstakingly generate a number of biaxial bending diagrams and somehow interpret these results, if only there was a way to create an interaction diagram for a constant ratio of m_x to m_y!

Solution:
Add an option to generate a moment interaction diagram for a given ratio m_x / m_y rather than a netural axis angle theta (theta is actually useless to the designer, but not to the programmer!) Thinking again about my original comment, would there not actually just be a single unknown that is theta? At least in the way I have implemented my code I can vary the neutral axis depth from zero to infinity, and for each depth just solve for theta that gives the desired moment ratio? I don't think in generating the interaction diagram I would have two unknowns. I can see that in the general non-linear case that given a design axial force and moments, the neutral axis depth and neutral axis angle are unknown, thus the need for something more involved, but we don't have this problem do we? Keen to hear your thoughts on this!

Issues with above solution:
There is currently an option to specify control points in the moment interaction diagram, two of these would be quite difficult with this new implementation.

• Given axial force, adds another unknown and would involve constrained optimisation. I don't see how nested root finding would work as neutral axis depth and neutral axis angle aren't independent variables when it comes to calculating the net axial force (this is a coupled problem).
• Given strain in an extreme bar (used to find balanced point etc.), again adds another unknown to the problem.

Anyway that's a bit of an unstructured rant, hoping it adds to this discussion and makes at least some sense (I guess this is what happens when you try to write a response five times and keep getting interrupted by your 6 week old daughter, how inconsiderate!) Keen to hear your thoughts and look forward to implementing a solution to this!

[robbievanleeuwen]

Thanks for the thoughts @buddyd16, unfortunately I don't have time to explore the biaxial stuff at the moment but hope to jump into it in the near future.

A few things still don't feel quite resolved for me yet, but they may become more apparent once I start to try to tackle the problem. E.g.

Don't want to rely on m_x/m_y ratio as there are instances where m_y can be zero and the neutral axis angle, theta, would also not be 0. I landed on comparing the angle between the applied moments and calculated moments.

I think you might have missed my idea here, I would use the brentq solver to find theta such that the desired m_x/m_y ratio is obtained. So in the example of the L-beam, if the user desired for m_y to be zero, for each neutral axis depth in the interaction diagram, theta would be found to satisfy this moment ratio requirement (theta would be non-zero as you correctly describe).

With a given strain you can solve for the required neutral axis depth directly from the linear strain equation so the only unknown would the neutral axis angle.

If the neutral axis angle theta is unknown, the neutral axis depth cannot be found. This would require some more thought but I feel like there's an answer in there somewhere.

Anyway, I'll be putting this on ice until I have some proper time to dedicate to this. However, I greatly value your experience on this one @buddyd16 and when I get back into it I hope to use these ideas to arrive at a good solution!

[buddyd16]

Thanks for the thoughts @buddyd16, unfortunately I don't have time to explore the biaxial stuff at the moment but hope to jump into it in the near future.

Yeah I've been stopped dead in my tracks on my little side programming endeavor, this is a great tool you've put together so far excited to see how it matures. I'll be peaking into the code when I can and hopefully offer more than just this discussion text.

A few things still don't feel quite resolved for me yet, but they may become more apparent once I start to try to tackle the problem. E.g.

With a given strain you can solve for the required neutral axis depth directly from the linear strain equation so the only unknown would the neutral axis angle.

If the neutral axis angle theta is unknown, the neutral axis depth cannot be found. This would require some more thought but I feel like there's an answer in there somewhere.

Speaking more towards the root finding process, if you go the option of nesting two brentq solvers you would pass a neutral axis and an optional required extreme tensile strain into the outer solver, if the strain was passed in solve for neutral axis depth directly from the linear strain equation if no strain passed in then proceed to the inner brentq solver for the neutral axis depth.

Anyway, I'll be putting this on ice until I have some proper time to dedicate to this. However, I greatly value your experience on this one @buddyd16 and when I get back into it I hope to use these ideas to arrive at a good solution!

Looking forward to being able to help out. You've made incredibly strides on the state of the tool already.

[Agent6-6-6]

To add to the discussion, until the analysis based on load angle is implemented the only true way to see how close you are to the curve is to compare on a biaxial plot.

To show the effect of the 'wander' of the interaction surface based on the constant neutral axis angle, you can see the effect at the base of this interaction surface created by concrete-properties (look at the black lines). The section was a 600 square column with 12-DH20 bars.

Where you can see the vertical black lines all tend to converge at 4 specific orientations, or further up they wander around the surface which means a lot less accuracy as you get near to the surface point especially if you're relying on the M/N curves and a load point lying inside this. You really should only be comparing $M_x$ and $M_y$ load data on the biaxial plane as this is correct

[robbievanleeuwen]

Hi @Agent6-6-6, cool plot! Would you mind sharing the code you used to generate this? It will be helpful for when I get to tackling this one 😄

[Agent6-6-6]

@robbievanleeuwen, sure something along these lines was used.

Just sub in the max compression/tension loads and strength reduction factors from the AS3600 code if using that (as the example code was generated based on #40). Note design_code=NZS3101() in the example.

Might be nice to have a standard method to generate something like this?

Note takes a while to run though!

# determine maximum compression load for the column section
phi = 0.85
max_comp = phi * design_code.max_comp_strength()
max_ten = phi * design_code.max_ten_strength()
steps = 15
# generate axial load list
n_list = np.linspace(-max_ten + 20000, max_comp, steps)
theta_list = np.linspace(0.0, 2 * np.pi, 48, False)

results_bb = []
results_mi = []

for theta in theta_list:
    f_mi_res, _, _ = design_code.moment_interaction_diagram(theta=theta)
    results_mi.append(f_mi_res)
for n in n_list:
    f_bb_res, _ = design_code.biaxial_bending_diagram(n=n)
    results_bb.append(f_bb_res)

# plot all the diagrams on one image
fig = plt.figure(figsize=(12, 10))
ax = plt.axes(projection="3d")
n_scale = 1e-3
m_scale = 1e-6
for mi_res in results_mi:
    n, m_x_list = mi_res.get_results_lists(moment="m_x")
    _, m_y_list = mi_res.get_results_lists(moment="m_y")
    n = np.array(n) * n_scale
    m_x_list = np.array(m_x_list) * m_scale
    m_y_list = np.array(m_y_list) * m_scale

    ax.plot3D(m_x_list, m_y_list, n, "-k", linewidth=0.5)

for i, bb_res in enumerate(results_bb):
    m_x_list, m_y_list = bb_res.get_results_lists()
    m_x_list = np.array(m_x_list) * m_scale
    m_y_list = np.array(m_y_list) * m_scale

    n = np.array(n_list[i]) * n_scale
    ax.plot3D(m_x_list, m_y_list, n, "-r", linewidth=0.5)

ax.set_xlabel("Bending Moment $M_x$")
ax.set_ylabel("Bending Moment $M_y$")
ax.set_zlabel("Axial Force $N$")

plt.show()