G2 Example 8

G2Example8.py

# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any kind.

# The code implemented in this file is underpinned by axial strain and axial
# stress, which are related to both axial force (P) and bending moment (M).
# The capacity of a cross-section to carry axial force and bending moment
# depends on several factors: its geometry, the materials it is made up of,
# and also the criteria used to define failure, i.e., the capacity limits.
# Importantly, the presence of axial force affects the bending moment capacity,
# and vice versa. This is why the objective in this document is to create a
# P-M interaction diagram.

# A rectangular solid steel cross-section is considered in this file, together
# with an elastic-perfectly-plastic uniaxial material model. That means we can
# apply different deformation scenarios to the cross-section, ramping up the
# deformations until all material fibres yield, and simply print the resulting
# axial force and bending moment as the relevant capacity values. (In this
# particular situation, if the P-M diagram was the only objective, a rigid-
# plastic model without an elastic region would be quicker and better.) The
# resulting capacity values can be compared with the limits Pu and Mu obtained
# from the "lower-bound theorem" of simple plastic capacity analysis. In fact,
# the axes of the P-M interaction diagram are scaled by Pu and Mu, which makes
# the graph independent of the cross-section geometry. The plot is therefore
# identical to that presented in Figure 6.21 in Filippou & Fenves' 2004 chapter
# entitled Methods of Analysis for Earthquake-Resistant Structures in a book
# edited by Bozorgnia and Bertero.

# P-M interaction diagrams for reinforced concrete cross-sections differ in
# several ways, including failure criteria. Typically, a maximum strain values
# in the steel material and the concrete material are used. One important point
# on the P-M interaction diagram for a reinforced concrete cross-section is the
# "balanced point" at which the concrete is at its maximum compressive strain
# and the steel is at its maximum tensile strain. That point is at or near the
# maximum moment capacity, and is associated with moderate axial compressive
# force. Many good references on creating P-M interaction diagram are available
# online, including Professor Michael Scott's Portwood Digital blog:
# https://portwooddigital.com/2022/06/12/p-m-interaction-by-the-book/

import numpy as np
import matplotlib.pyplot as plt
from G2SectionRectangular import *

# Material properties (N, mm)
E = 200e3
fy = 350
alpha = 0.0

# Cross-section geometry (mm)
h = 100
b = 50

# Number of fibres (number of integration points)
n = 20

# Independent ultimate capacities
Mu = fy * b * (h/2)**2
Pu = fy * b * h
print('\n'"Independent ultimate capacities: Mu=%.0fkNm, Pu=%.0fkN" % (Mu/1e6, Pu/1e3))

# Yield strain and curvature
epsy = fy/E
kappay = epsy / (h/2)

# Storage for final moment and axial force values
finalMstorage = []
finalPstorage = []

# Curvature increments
kappaIncrements = kappay / 10

# Let deformation scenarios be defined by location of zero strain
for i in range(11):

    # Location with zero strain
    location = i / 10 * h

    # Initialize section matrices
    us = np.zeros((2, 3))
    Ks = np.eye(2)

    # Create and recreate the cross-section
    theSection = rectangularSection(h, b, n, ['Bilinear', E, fy, alpha])

    # Initialize storage for forces and deformations
    epsStorage = []
    kappaStorage = []
    Pstorage = []
    Mstorage = []

    # Swiftly deal with locations 0 and h, i.e., full compression or tension
    if location == 0 or location == h:

        # Set incremental section deformations
        us[0, 1] = 0
        if location == 0:
            us[1, 1] = epsy
        else:
            us[1, 1] = -epsy

        # Set total section deformations
        us[0, 0] += us[0, 1]
        us[1, 0] += us[1, 1]

        # State determination
        Fs, Ks = theSection.state(us)

        # Commit history variables
        theSection.commit()

        # Store forces and deformations
        kappaStorage.append(us[0, 0])
        epsStorage.append(us[1, 0])
        Mstorage.append(Fs[0])
        Pstorage.append(Fs[1])

    else:

        # Add curvature gradually, stop when cross-section has zero stiffness
        while np.linalg.det(Ks) > 0:

            # Set incremental section deformations
            us[0, 1] = kappaIncrements
            us[1, 1] = kappaIncrements * (h/2 - location)

            # Set total section deformations
            us[0, 0] += us[0, 1]
            us[1, 0] += us[1, 1]

            # State determination
            Fs, Ks = theSection.state(us)

            # Commit history variables
            theSection.commit()

            # Store forces and deformations
            kappaStorage.append(us[0, 0])
            epsStorage.append(us[1, 0])
            Mstorage.append(Fs[0])
            Pstorage.append(Fs[1])

    # Store final moment and axial force values for interaction diagram
    finalMstorage.append(Fs[0])
    finalPstorage.append(Fs[1])

# Obtain points on the other side of the interaction diagram, i.e., for opposite sign of moment
for i in range(9):

    # Location with zero strain
    location = (i+1) / 10 * h

    # Initialize section matrices
    us = np.zeros((2, 3))
    Ks = np.eye(2)

    # Create and recreate the cross-section
    theSection = rectangularSection(h, b, n, ['Bilinear', E, fy, alpha])

    # Initialize storage for forces and deformations
    epsStorage = []
    kappaStorage = []
    Pstorage = []
    Mstorage = []

    # Add curvature gradually, stop when cross-section has zero stiffness
    while np.linalg.det(Ks) > 0:

        # Set incremental section deformations
        us[0, 1] = -kappaIncrements
        us[1, 1] = -kappaIncrements * (h/2 - location)

        # Set total section deformations
        us[0, 0] += us[0, 1]
        us[1, 0] += us[1, 1]

        # State determination
        Fs, Ks = theSection.state(us)

        # Commit history variables
        theSection.commit()

        # Store forces and deformations
        kappaStorage.append(us[0, 0])
        epsStorage.append(us[1, 0])
        Mstorage.append(Fs[0])
        Pstorage.append(Fs[1])

    # Store final moment and axial force values for interaction diagram
    finalMstorage.append(Fs[0])
    finalPstorage.append(Fs[1])

# Connect the last point to the first, for the sake of the plot
finalMstorage.append(finalMstorage[0])
finalPstorage.append(finalPstorage[0])

# Analytical solution
M = []
P = []
granularity = 100
for widthP in np.linspace(h, 0, granularity):
    widthM = (h-widthP)/2
    M.append(widthM*b*fy*(widthP+widthM))
    P.append(widthP*b*fy)

for widthP in np.linspace(h/granularity, h, granularity):
    widthM = (h-widthP)/2
    M.append(widthM*b*fy*(widthP+widthM))
    P.append(-widthP*b*fy)

for widthP in np.linspace(h-h/granularity, 0, granularity):
    widthM = (h-widthP)/2
    M.append(-widthM*b*fy*(widthP + widthM))
    P.append(-widthP*b*fy)

for widthP in np.linspace(h/granularity, h, granularity):
    widthM = (h-widthP)/2
    M.append(-widthM*b*fy*(widthP + widthM))
    P.append(widthP*b*fy)

# Plot
plt.ion()
plt.figure()
plt.title('P-M Interaction Diagram')
plt.grid(True)
plt.autoscale(True)
plt.ylabel('$P/P_u$')
plt.xlabel('$M/M_u$')
plt.plot(np.divide(finalMstorage, Mu), np.divide(finalPstorage, Pu), 'bo-', linewidth=1.0, markersize=5, label="From force-deformation curves")
plt.plot(np.divide(M, Mu), np.divide(P, Pu), 'r-', linewidth=1.0, label="Analytical solution")
plt.legend(loc='center')
print('\n'"Click somewhere in the plot to continue...")
plt.waitforbuttonpress()