{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Approximation of Torsion Constant for Trapezoidal Sections\n\nTrapezoidal elements or components of a cross-section are\nquite common in bridge structures, either concrete or steel\ncomposite construction. However, it's common to determine \nthe torsion constant of the trapezoidal section by using a\nrectangular approximation. For example, this is done in\nthe Autodesk Structural Bridge Design software when there\nis a haunch in a [Steel Composite Beam](https://knowledge.autodesk.com/support/structural-bridge-design/learn-explore/caas/CloudHelp/cloudhelp/ENU/ASBD-InProdAU/files/structure/tech-info/ti-torsion/ASBD-InProdAU-structure-tech-info-ti-torsion-Torsion-property-for-SAM-composite-beams-html-html.html)\n\nThe question then arises, when is it appropriate to make\nthe rectangular approximation to a trapezoidal section,\nand what might the expected error be?\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Define the Imports\nHere we bring in the primative section shapes and also the more generic\nShapely `Polygon` object.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom shapely import Polygon\nimport sectionproperties.pre.geometry as geometry\nimport sectionproperties.pre.pre as pre\nimport sectionproperties.pre.library.primitive_sections as sections\nfrom sectionproperties.analysis.section import Section"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Define the Calculation Engine\nIt's better to collect the relevant section property calculation in a\nsingle function. We are only interested in the torsion constant, so this\nis straightforward enough. `geom` is the Section Property Geometry object;\n`ms` is the mesh size.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def get_section_j(geom, ms, plot_geom=False):\n    geom.create_mesh(mesh_sizes=[ms])\n    section = Section(geom)\n    if plot_geom:\n        section.plot_mesh()\n    section.calculate_geometric_properties()\n    section.calculate_warping_properties()\n    return section.get_j()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Define the Mesh Density\nThe number of elements per unit area is an important input to the calculations\neven though we are only examining ratios of the results. A nominal value of 100\nis reasonable.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n = 100  # mesh density"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Create and Analyse the Section\nThis function accepts the width `b` and a slope `S` to create the trapezoid.\nSince we are only interested in relative results, the nominal dimensions are\nimmaterial. There are a few ways to parametrize the problem, but it has been\nfound that setting the middle height of trapezoid (i.e. the average height)\nto a unit value works fine.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def do_section(b, S, d_mid=1, plot_geom=False):\n    delta = S * d_mid\n    d1 = d_mid - delta\n    d2 = d_mid + delta\n\n    # compute mesh size\n    ms = d_mid * b / n\n\n    points = []\n    points.append([0, 0])\n    points.append([0, d1])\n    points.append([b, d2])\n    points.append([b, 0])\n    if S < 1.0:\n        trap_geom = geometry.Geometry(Polygon(points))\n    else:\n        trap_geom = sections.triangular_section(h=d2, b=b)\n    jt = get_section_j(trap_geom, ms, plot_geom)\n\n    rect_geom = sections.rectangular_section(d=(d1 + d2) / 2, b=b)\n    jr = get_section_j(rect_geom, ms, plot_geom)\n    return jt, jr, d1, d2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Example Section\nThe analysis for a particular section looks as follows:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "b, S = 4.0, 0.3\njt, jr, d1, d2 = do_section(b, S, plot_geom=True)\nprint(f\"{b=}; {S=}; {jr=}; {jt=}; {jr/jt}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Create Loop Variables\nThe slope `S` is 0 for a rectangle, and 1 for a triangle and is\ndefined per the plot below. A range of `S`, between 0.0 and 1.0 and\na range of `b` are considered, between 1 and 10 here (but can be extended)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "b_list = np.logspace(0, np.log10(10), 10)\nS_list = np.linspace(0.0, 1.0, 10)\nj_rect = np.zeros((len(b_list), len(S_list)))\nj_trap = np.zeros((len(b_list), len(S_list)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The Main Loop\nExecute the double loop to get the ratios for combinations of `S` and `b`.\n\nAn optional deugging line is left in for development but commented out.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for i, b in enumerate(b_list):\n    for j, S in enumerate(S_list):\n        jt, jr, d1, d2 = do_section(b, S)\n        j_trap[i][j] = jt\n        j_rect[i][j] = jr\n\n        # print(f\"{b=:.3}; {S=:.3}; {d1=:.5}; {d2=:.5}; J rect = {j_rect[i][j]:.5e}; J trap = {j_trap[i][j]:.5e}; J ratio = {j_rect[i][j]/j_trap[i][j]:.5}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Calculate the Ratios\nCourtesy of numpy, this is easy:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "j_ratio = j_rect / j_trap"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot the Results\nHere we highlight a few of the contours to illustrate the accuracy and behaviour\nof the approximation.\n\nAs expected, when the section is rectangular, the error is small, but as it increases\ntowards a triangle the accuracy generally reduces. However, there is an interesting\nline at an aspect ratio of about 2.7 where the rectangular approximation is always\nequal to the trapezoid's torsion constant.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "levels = np.arange(start=0.5, stop=1.5, step=0.05)\nplt.figure(figsize=(12, 6))\ncs = plt.contour(\n    S_list,\n    b_list,\n    j_ratio,\n    levels=[0.95, 0.99, 1.00, 1.01, 1.05],\n    colors=(\"k\",),\n    linestyles=(\":\",),\n    linewidths=(1.2,),\n)\nplt.clabel(cs, colors=\"k\", fontsize=10)\nplt.contourf(S_list, b_list, j_ratio, 25, cmap=\"Wistia\", levels=levels)\n# plt.yscale('log')\nminor_x_ticks = np.linspace(min(S_list), max(S_list), 10)\n# plt.xticks(minor_x_ticks,minor=True)\nplt.minorticks_on()\nplt.grid(which=\"both\", ls=\":\")\nplt.xlabel(r\"Slope $S = (d_2-d_1)/(d_2+d_1); d_2\\geq d_1, d_1\\geq 0$\")\nplt.ylabel(\"Aspect $b/d_{ave}; d_{ave} = (d_1 + d_2)/2$\")\nplt.colorbar()\nplt.title(\n    r\"Accuracy of rectangular approximation to trapezoid torsion constant $J_{rect}\\, /\\, J_{trapz}$\",\n    multialignment=\"center\",\n)\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}