Add biaxial bending diagram

Author: robbievanleeuwen

README.md

@@ -61,9 +61,8 @@ documentation can found at [https://robbievanleeuwen.github.io/concrete-properti
 - [x] Squash load
 - [x] Tensile load
 - [x] Moment interaction diagrams
-  - [x] Basic M-N curve
-  - [ ] Combined M-N curve (sagging & hogging)
-  - [ ] Biaxial bending curve
+  - [x] M-N curves
+  - [x] Biaxial bending curve
 
 ### Stress Analysis
 - [ ] Uncracked stresses

concreteproperties/concrete_section.py

@@ -632,7 +632,9 @@ def calculate_service_section_actions(
         self,
         d_n: float,
         kappa: float,
-        moment_curvature: res.MomentCurvatureResults = res.MomentCurvatureResults(),
+        moment_curvature: res.MomentCurvatureResults = res.MomentCurvatureResults(
+            theta=0
+        ),
     ) -> res.MomentCurvatureResults:
         """Given a neutral axis depth `d_n` and curvature `kappa`, calculates the
         resultant axial force and bending moment.
@@ -782,7 +784,7 @@ def ultimate_bending_capacity(
         b = d_t  # neutral axis at extreme tensile fibre
 
         # initialise ultimate bending results
-        ultimate_results = res.UltimateBendingResults()
+        ultimate_results = res.UltimateBendingResults(theta=theta)
 
         # find neutral axis that gives convergence of the axial force
         (d_n, r) = brentq(
@@ -805,8 +807,7 @@ def ultimate_normal_force_convergence(
         ultimate_results: results.UltimateBendingResults,
     ) -> float:
         """Given a neutral axis depth `d_n` and neutral axis angle `theta`, calculates
-        the difference between the target net axial force `n` and the axial force
-        given `d_n` & `theta`.
+        the difference between the target net axial force `n` and the axial force.
 
         :param float d_n: Depth of the neutral axis from the extreme compression fibre
         :param float n: Net axial force
@@ -829,7 +830,9 @@ def ultimate_normal_force_convergence(
     def calculate_ultimate_section_actions(
         self,
         d_n: float,
-        ultimate_results: results.UltimateBendingResults,
+        ultimate_results: res.UltimateBendingResults = res.UltimateBendingResults(
+            theta=0
+        ),
     ) -> results.UltimateBendingResults:
         """Given a neutral axis depth `d_n` and neutral axis angle `theta`, calculates
         the resultant bending moments `mx`, `my`, `mv` and the net axial force `n`.
@@ -962,13 +965,16 @@ def calculate_ultimate_section_actions(
     def moment_interaction_diagram(
         self,
         theta: float = 0,
+        m_neg: bool = False,
         n_points: int = 24,
     ) -> res.MomentInteractionResults:
         """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 float theta: Angle the neutral axis makes with the horizontal axis
+        :param bool m_neg: If set to true, also calculates the moment interaction for
+            theta = theta + pi, i.e. positive and negative bending
         :param int n_points: Number of calculation points between the decompression
             case and the pure bending case.
 
@@ -994,9 +1000,14 @@ def moment_interaction_diagram(
 
         # create progress bar
         with utils.create_known_progress() as progress:
+            progress_length = n_points
+
+            if m_neg:
+                progress_length *= 2
+
             task = progress.add_task(
                 description="[red]Generating M-N diagram",
-                total=n_points,
+                total=progress_length,
             )
 
             for d_n in d_n_list:
@@ -1007,16 +1018,100 @@ def moment_interaction_diagram(
                 mi_results.m.append(ultimate_results.mv)
                 progress.update(task, advance=1)
 
+            if not m_neg:
+                progress.update(
+                    task,
+                    description="[bold green]:white_check_mark: M-N diagram generated",
+                )
+
+            # add tensile load
+            mi_results.n.append(self.gross_properties.tensile_load)
+            mi_results.m.append(0)
+
+            # if calculating negative bending
+            if m_neg:
+                theta += np.pi
+                # compute extreme tensile fibre
+                _, d_t = utils.calculate_extreme_fibre(
+                    points=self.geometry.points, theta=theta
+                )
+
+                # compute neutral axis depth for pure bending case
+                ultimate_results = self.ultimate_bending_capacity(theta=theta, n=0)
+
+                # generate list of neutral axes
+                d_n_list = np.linspace(
+                    start=d_t, stop=ultimate_results.d_n, num=n_points
+                )
+
+                for d_n in reversed(d_n_list):
+                    ultimate_results = self.calculate_ultimate_section_actions(
+                        d_n=d_n, ultimate_results=ultimate_results
+                    )
+                    mi_results.n.append(ultimate_results.n)
+                    mi_results.m.append(-ultimate_results.mv)
+                    progress.update(task, advance=1)
+
+                progress.update(
+                    task,
+                    description="[bold green]:white_check_mark: M-N diagram generated",
+                )
+
+                # add squash load
+                mi_results.n.append(self.gross_properties.squash_load)
+                mi_results.m.append(0)
+
+        return mi_results
+
+    def biaxial_bending_diagram(
+        self,
+        n: float = 0,
+        n_points: int = 48,
+    ) -> res.BiaxialBendingResults:
+        """Generates a biaxial bending diagram given a net axial force `n` and
+        `n_points` calculation points.
+
+        :param float n: Net axial force
+        :param int n_points: Number of calculation points between the decompression
+            case and the pure bending case.
+
+        :return: Biaxial bending results
+        :rtype: :class:`~concreteproperties.results.BiaxialBendingResults`
+        """
+
+        # initialise results
+        bb_results = res.BiaxialBendingResults(n=n)
+
+        # calculate d_theta
+        d_theta = 2*np.pi / n_points
+
+        # generate list of thetas
+        theta_list = np.linspace(start=0, stop=2*np.pi-d_theta, num=n_points)
+
+        # create progress bar
+        with utils.create_known_progress() as progress:
+            task = progress.add_task(
+                description="[red]Generating biaxial bending diagram",
+                total=n_points,
+            )
+
+            # loop through thetas
+            for theta in theta_list:
+                ultimate_results = self.ultimate_bending_capacity(theta=theta, n=n)
+                bb_results.mx.append(ultimate_results.mx)
+                bb_results.my.append(ultimate_results.my)
+                progress.update(task, advance=1)
+
+            # add first result to end of list top
+            bb_results.mx.append(bb_results.mx[0])
+            bb_results.my.append(bb_results.my[0])
+
             progress.update(
                 task,
-                description="[bold green]:white_check_mark: M-N diagram generated",
+                description="[bold green]:white_check_mark: Biaxial bending diagram generated",
             )
 
-        # add tensile load
-        mi_results.n.append(self.gross_properties.tensile_load)
-        mi_results.m.append(0)
-
-        return mi_results
+        return bb_results
 
     def get_c_local(
         self,

concreteproperties/results.py

@@ -184,10 +184,13 @@ def calculate_transformed_properties(
 
 @dataclass
 class MomentCurvatureResults:
-    """Class for storing moment curvature results."""
+    """Class for storing moment curvature results.
+
+    :param float theta: Angle the neutral axis makes with the horizontal axis
+    """
 
     # results
-    theta: float = 0
+    theta: float
     kappa: List[float] = field(default_factory=list)
     moment: List[float] = field(default_factory=list)
 
@@ -234,9 +237,12 @@ def plot_results(
 
 @dataclass
 class UltimateBendingResults:
-    """Class for storing ultimate bending results."""
+    """Class for storing ultimate bending results.
+
+    :param float theta: Angle the neutral axis makes with the horizontal axis
+    """
 
-    theta: float = 0
+    theta: float
     d_n: float = 0
     n: float = 0
     mx: float = 0
@@ -329,3 +335,45 @@ def plot_multiple_diagrams(
                 ax.legend(loc="center left", bbox_to_anchor=(1, 0.5))
 
         return ax
+
+@dataclass
+class BiaxialBendingResults:
+    """Class for storing biaxial bending results.
+
+    :param float n: Net axial force
+    """
+
+    n: float
+    mx: List[float] = field(default_factory=list)
+    my: List[float] = field(default_factory=list)
+
+    def plot_diagram(
+        self,
+        m_scale: float = 1e-6,
+        **kwargs,
+    ) -> matplotlib.axes._subplots.AxesSubplot:
+        """Plots a biaxial bending diagram.
+
+        :param float m_scale: Scaling factor to apply to bending moment
+        :param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
+
+        :return: Matplotlib axes object
+        :rtype: :class:`matplotlib.axes._subplots.AxesSubplot`
+        """
+
+        # create plot and setup the plot
+        with plotting_context(title=f"Biaxial Bending Diagram, $N = {self.n:.3e}$", **kwargs) as (
+            fig,
+            ax,
+        ):
+            # scale results
+            mx = np.array(self.mx) * m_scale
+            my = np.array(self.my) * m_scale
+
+            ax.plot(mx, my, "o-", markersize=3)
+
+            plt.xlabel("Bending Moment $M_x$")
+            plt.ylabel("Bending Moment $M_y$")
+            plt.grid(True)
+
+        return ax

examples/ex1.py

@@ -41,7 +41,7 @@
     d=600,
     n=32,
     dia=20,
-    n_bar=6,
+    n_bar=7,
     n_circle=4,
     area_conc=np.pi * 600 * 600 / 4,
     area_bar=310,
@@ -51,6 +51,6 @@
 )
 
 conc_sec = ConcreteSection(geometry)
-
 print(conc_sec.ultimate_bending_capacity())
 conc_sec.moment_interaction_diagram().plot_diagram()
+conc_sec.biaxial_bending_diagram(n=4000e3).plot_diagram()

examples/ex2.py

@@ -43,13 +43,13 @@
     b=400,
     d=600,
     dia_top=16,
-    n_top=6,
-    dia_bot=16,
-    n_bot=6,
+    n_top=3,
+    dia_bot=24,
+    n_bot=3,
     n_circle=4,
-    cover=66,
+    cover=30,
     area_top=200,
-    area_bot=200,
+    area_bot=450,
     conc_mat=concrete,
     steel_mat=steel,
 )
@@ -58,7 +58,7 @@
 
 pprint(conc_sec.gross_properties)
 pprint(conc_sec.get_transformed_gross_properties(elastic_modulus=32.8e3))
-mi_res = conc_sec.moment_interaction_diagram()
+mi_res = conc_sec.moment_interaction_diagram(m_neg=True)
 pprint(mi_res)
 mi_res.plot_diagram()
 cracked_res = conc_sec.calculate_cracked_properties()