Plot Beam Diagrams

PlotBeamDiagrams.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 any form of warranty.
# Please see the Programming note for how to get started, and notice
# that you must copy certain functions into the file terjesfunctions.py
# ------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
import terjesfunctions as fcns


# ------------------------------------------------------------------------
# INPUT
# ------------------------------------------------------------------------

E = 200000.0
I = 80900000.0
A = 9100.0
L = 6000.0

# Distributed load [negative means downwards]
q = -10.0

# Local displacement vector from structural analysis
#               u1  u2     u3     u4   u5    u6
ul = np.array([0.0, 3.2, 0.0002, 0.0, 5.2, 0.0003])


# ------------------------------------------------------------------------
# CALCULATE RESULTS USING MATRIX STRUCTURAL ANALYSIS
# ------------------------------------------------------------------------

# Get local stiffness matrix
[Kl, Kgeometric, Tlg] = fcns.getFrameStiffness2D(E, 0, I, A, 0, L, L, 0, 0)

# End forces from displacement vector, without fixed-end moments from distributed load
Fl = Kl.dot(ul)

# Add contributions from distributed loads [signs to match plots]
print('\n'"Moment at the left end from matrix analysis [kNm]:", -(Fl[2] + q * L**2 / 12.0)/1e6)
print("Moment at the right end from matrix analysis [kNm]:", (Fl[5] - q * L**2 / 12.0)/1e6)

print('\n'"Shear force at the left end from matrix analysis [kN]:", (Fl[1] - q * L / 2)/1e3)
print("Shear force at the right end from matrix analysis [kN]:", (-Fl[4] + q * L / 2)/1e3)


# ------------------------------------------------------------------------
# CALCULATE RESULTS USING THE DIFFERENTIAL EQUATION
# ------------------------------------------------------------------------
xData = []
wData = []
thetaData = []
MData = []
VData = []
beamLine = []
num = 20
xi = np.linspace(0,1,num)
for i in range(len(xi)):

    # x-value corresponding to xi
    x = xi[i] * L
    xData.append(x)

    # Reference line for the plots
    beamLine.append(0)

    # Bending stiffness and axial stiffness
    EI = E * I
    EA = E * A

    # Rename the displacements for clarity
    u1 = ul[0]
    u2 = ul[1]
    u3 = ul[2]
    u4 = ul[3]
    u5 = ul[4]
    u6 = ul[5]

    # Coefficients in the solution to the differential equation
    C1 = -q*L/(12.0*EI) - 2.0*u5/L**3 - u3/L**2 + 2.0*u2/L**3 - u6/L**2
    C2 = q*L*L/(24.0*EI) + 3.0*u5/L**2 + 2.0*u3/L - 3.0*u2/L**2 + u6/L
    C3 = -u3
    C4 = u2

    # Displacement
    wData.append(q * x**4/(24.0 * EI) + C1 * x**3 + C2 * x**2 + C3 * x + C4)

    # Rotation [clockwise positive]
    thetaData.append(q * x**3/(6.0*EI) + 3.0*C1 * x**2 + 2.0 * C2 * x + C3)

    # Bending moment [on the tension side]
    MData.append(-(q * x**2/2.0 + EI * 6.0 * C1 * x + 2.0 * EI * C2)/1e6)

    # Shear force [clockwise positive]
    VData.append((q * x + EI * 6.0 * C1)/1000.0)

print('\n'"Moment at the left end from diff. eq. [kNm]:", MData[0])
print("Moment at the right end from diff. eq. [kNm]:", MData[num-1])

print('\n'"Shear force at the left end from diff. eq. [kN]:", VData[0])
print("Shear force at the right end from diff. eq. [kN]:", VData[num-1])


# ------------------------------------------------------------------------
# PLOT DIAGRAMS
# ------------------------------------------------------------------------
fig = plt.figure()
plt.axis('off')
plt.title('Deformations and Forces in Beam')

axes1 = fig.add_subplot(411)
axes1.get_xaxis().set_visible(False)
plt.plot(xData, wData, 'k-', linewidth=1.0)
plt.plot(xData, beamLine, 'k-', linewidth=2.0)
plt.plot([xData[0], xData[0]], [beamLine[0], wData[0]], 'k-', linewidth=1.0)
plt.plot([xData[num-1], xData[num-1]], [beamLine[num-1], wData[num-1]], 'k-', linewidth=1.0)
plt.ylabel('w [mm]')

axes2 = fig.add_subplot(412)
axes2.get_xaxis().set_visible(False)
plt.plot(xData, thetaData, 'k-', linewidth=1.0)
plt.plot(xData, beamLine, 'k-', linewidth=2.0)
plt.plot([xData[0], xData[0]], [beamLine[0], thetaData[0]], 'k-', linewidth=1.0)
plt.plot([xData[num-1], xData[num-1]], [beamLine[num-1], thetaData[num-1]], 'k-', linewidth=1.0)
plt.ylabel(r'$\theta$ [rad]')

axes3 = fig.add_subplot(413)
axes3.get_xaxis().set_visible(False)
plt.plot(xData, MData, 'k-', linewidth=1.0)
plt.plot(xData, beamLine, 'k-', linewidth=2.0)
plt.plot([xData[0], xData[0]], [beamLine[0], MData[0]], 'k-', linewidth=1.0)
plt.plot([xData[num-1], xData[num-1]], [beamLine[num-1], MData[num-1]], 'k-', linewidth=1.0)
plt.ylabel('M [kNm]')

axes4 = fig.add_subplot(414)
plt.plot(xData, VData, 'k-', linewidth=1.0)
plt.plot(xData, beamLine, 'k-', linewidth=2.0)
plt.plot([xData[0], xData[0]], [beamLine[0], VData[0]], 'k-', linewidth=1.0)
plt.plot([xData[num-1], xData[num-1]], [beamLine[num-1], VData[num-1]], 'k-', linewidth=1.0)
plt.ylabel('V [kN]')
plt.xlabel('x')

plt.show()