 codex/beam-support-compliance.py

"""Original linear-beam comparison, writing its PNG in the current directory."""

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from numpy.polynomial import Polynomial


def main():
    simple_left = Polynomial([0, 3 / 48, 0, -4 / 48])
    clamped_left = Polynomial([0, 0, 3 / 48, -4 / 48])
    np.testing.assert_allclose(simple_left(0), 0, atol=1e-14)
    np.testing.assert_allclose(simple_left.deriv(2)(0), 0, atol=1e-14)
    np.testing.assert_allclose(simple_left.deriv()(0.5), 0, atol=1e-14)
    np.testing.assert_allclose(clamped_left(0), 0, atol=1e-14)
    np.testing.assert_allclose(clamped_left.deriv()(0), 0, atol=1e-14)
    np.testing.assert_allclose(clamped_left.deriv()(0.5), 0, atol=1e-14)
    np.testing.assert_allclose(simple_left(0.5), 1 / 48)
    np.testing.assert_allclose(clamped_left(0.5), 1 / 192)
    np.testing.assert_allclose(-2 * simple_left.deriv(3)(0.5), 1)
    np.testing.assert_allclose(-2 * clamped_left.deriv(3)(0.5), 1)

    length = 0.5
    modulus = 200e9
    width = 0.02
    thickness = 0.003
    load = 1.0
    inertia = width * thickness**3 / 12
    rigidity = modulus * inertia
    np.testing.assert_allclose(rigidity, 9)
    scale_mm = load * length**3 / rigidity * 1000
    x = np.linspace(0, 1, 1001)
    folded = np.minimum(x, 1 - x)
    simple = simple_left(folded) * scale_mm
    clamped = clamped_left(folded) * scale_mm
    np.testing.assert_allclose(simple[500] / clamped[500], 4)

    plt.rcParams.update({"font.size": 11})
    fig, ax = plt.subplots(figsize=(8.8, 4.4), dpi=100, facecolor="white")
    fig.subplots_adjust(left=0.11, right=0.98, top=0.78, bottom=0.16)
    fig.suptitle("Support rotation changes beam compliance", fontsize=15, y=0.97)
    fig.text(0.5, 0.875,
             "Original model: L = 0.50 m, E = 200 GPa, b = 20 mm, h = 3 mm; added P = 1 N",
             ha="center", fontsize=10)
    ax.axhline(0, color="#777777", lw=1, ls="--")
    ax.plot(x, simple, color="#0072b2", lw=2.5, label="simply supported")
    ax.plot(x, clamped, color="#d55e00", lw=2.5, label="both ends clamped")
    ax.plot(0.5, simple[500], "o", color="#0072b2")
    ax.plot(0.5, clamped[500], "o", color="#d55e00")
    ax.annotate("0.289 mm", xy=(0.5, simple[500]), xytext=(0.66, 0.285),
                arrowprops={"arrowstyle": "->", "color": "#555555"})
    ax.annotate("0.0723 mm", xy=(0.5, clamped[500]), xytext=(0.66, 0.10),
                arrowprops={"arrowstyle": "->", "color": "#555555"})
    ax.annotate("P", xy=(0.5, 0.04), xytext=(0.5, -0.007), ha="center",
                arrowprops={"arrowstyle": "->", "color": "#333333"})
    ax.set(xlabel="Position along span x/L", ylabel="Downward deflection (mm)",
           xlim=(0, 1), ylim=(0.325, -0.025))
    ax.grid(alpha=0.2)
    ax.set_facecolor("white")
    ax.spines[["top", "right"]].set_visible(False)
    ax.legend(loc="lower left", fontsize=10)
    fig.savefig(Path.cwd() / "beam-support-compliance.png", dpi=100,
                transparent=False, facecolor="white")
    plt.close(fig)


if __name__ == "__main__":
    main()