Compute section properties of a circle by using boundary method

ComputeSectionCircleValues.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing the section values of a circle\n",
    "\n",
    "This notebook computes the caracteristic values of a circular section.\n",
    "We use [this example](https://sectionproperties.readthedocs.io/en/latest/sphinx_gallery_examples/00-basic/01_simple.html#sphx-glr-sphinx-gallery-examples-00-basic-01-simple-py) to compare the values.\n",
    "\n",
    "The idea is from [this video](www.youtube.com/watch?v=Y2cBpBLJlfU): instead of computing the integral on the area (like in the example with many triangles), we transform the integral over the domain into a integral over the boundary by using Green's theorem:\n",
    "\n",
    "$$\\oint_{\\partial D} P \\ dx + Q \\ dy = \\iint_{D} \\left(\\dfrac{\\partial Q}{\\partial x} - \\dfrac{\\partial P}{\\partial y}\\right) \\ dx \\ dy$$\n",
    "\n",
    "All the equations were gotten from the [theorical documentation](https://sectionproperties.readthedocs.io/en/latest/rst/theory.html)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial configuration\n",
    "\n",
    "<img src=\"https://sectionproperties.readthedocs.io/en/latest/_images/arbitrarySection.png\" width=\"400\"></img>\n",
    "\n",
    "For an arbitrary cross-section, we can parametrize the closed path by \n",
    "\n",
    "$$t: \\left[0, \\ 1\\right] \\to p(t) = (x(t), \\ y(t))$$\n",
    "\n",
    "which $p(0) = p(1)$. We suppose this path doesn't intersect it self, except at the extremities.\n",
    "\n",
    "To our example, a circular path is given by:\n",
    "\n",
    "$$x(t) = \\dfrac{1}{2} d \\cdot \\cos 2\\pi t$$\n",
    "$$y(t) = \\dfrac{1}{2} d \\cdot \\sin 2\\pi t$$\n",
    "\n",
    "Then we discretize this path in $n$ distinct points. The easiest way to do it is by uniform distribuition:\n",
    "\n",
    "$$t_i = \\dfrac{i}{n} \\ \\ \\ \\ i = 0, \\ 1, \\ \\cdots, \\ n$$\n",
    "$$p_i = \\left(x_i , \\ y_i\\right)$$\n",
    "$$x_i = \\dfrac{1}{2}d \\cdot \\cos \\dfrac{2 i \\pi}{n}$$\n",
    "$$y_I = \\dfrac{1}{2}d \\cdot \\sin \\dfrac{2 i \\pi}{n}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "diameter, npts = 50, 64\n",
    "theta = 2*np.pi*np.linspace(0, 1, npts+1)\n",
    "x = 0.5*diameter*np.cos(theta)\n",
    "y = 0.5*diameter*np.sin(theta)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing values\n",
    "\n",
    "##### Perimeter\n",
    "\n",
    "$$L = \\oint_{\\partial D} ds = \\int_{0}^{1} \\sqrt{x'(t)^2 + y'(t)^2} \\ dt \\approx \\sum_{i=0}^{n} \\sqrt{\\left(x_{i+1}-x_{i}\\right)^2 + \\left(y_{i+1}-y_{i}\\right)^2}$$\n",
    "\n",
    "##### Cross-Sectional Area\n",
    "\n",
    "By setting $P=0$ and $Q=x$ we get\n",
    "\n",
    "$$A = \\int_{D} 1 \\ dx \\ dy = \\oint_{\\partial D} x \\ dy = \\int_{0}^{1} x(t) \\cdot y'(t) \\ dt$$\n",
    "\n",
    "Since we discretized it already by $n$ lines, at the interval $\\left[t_i, \\ t_{i+1}\\right]$ we get\n",
    "\n",
    "$$u = \\dfrac{t-t_i}{t_{i+1}-t_{i}}$$\n",
    "$$\\bar{x}(u) = \\left(1-u\\right) \\cdot x_{i} + u \\cdot x_{i+1}$$\n",
    "$$\\bar{y}(u) = \\left(1-u\\right) \\cdot y_{i} + u \\cdot y_{i+1}$$\n",
    "\n",
    "$$A \\approx \\sum_{i=0}^{n} \\int_{t_i}^{t_{i+1}} x(t) \\cdot y'(t) \\ dt = \\sum_{i=0}^{n} \\int_{0}^{1} \\bar{x}(u) \\cdot \\bar{y}'(u) \\ du$$\n",
    "$$\\boxed{A \\approx \\sum_{i=0}^{n} \\dfrac{1}{2}\\left(x_{i}+x_{i+1}\\right) \\cdot \\left(y_{i+1}-y_i\\right)}$$\n",
    "\n",
    "or by choosing $P=-y$ and $Q=0$:\n",
    "\n",
    "$$\\boxed{A \\approx \\sum_{i=0}^{n} \\dfrac{1}{2}\\left(x_{i+1}-x_i\\right) \\cdot \\left(y_i+y_{i+1}\\right)}$$\n",
    "\n",
    "##### First Moments of Area\n",
    "\n",
    "Doing the same as before, but ommiting many computations:\n",
    "\n",
    "* With $Q=0$ and $P=-\\dfrac{1}{2}y^2$\n",
    "$$Q_x = \\int_{D} \\ y \\ dx \\ dy = \\int_{\\partial D} -\\dfrac{1}{2}y^2 \\ dx$$\n",
    "$$\\boxed{Q_x\\approx -\\sum_{i=0}^{n} \\dfrac{1}{6}\\left(x_{i+1}-x_{i}\\right)\\cdot \\left(y_{i}^2+y_i y_{i+1} + y_{i+1}^2\\right)} $$\n",
    "\n",
    "* With $Q=\\dfrac{1}{2}x^2$ and $P=0$\n",
    "$$Q_y = \\int_{D} \\ x \\ dx \\ dy = \\int_{\\partial D} -\\dfrac{1}{2}x^2 \\ dy$$\n",
    "$$\\boxed{Q_{y}\\approx \\sum_{i=0}^{n} \\dfrac{1}{6} \\left(x_{i}^2+x_i x_{i+1} + x_{i+1}^2\\right)\\cdot \\left(y_{i+1}-y_{i}\\right)} $$\n",
    "\n",
    "##### Second Moment of Area\n",
    "\n",
    "* With $Q=0$ and $P=-\\dfrac{1}{3}y^3$\n",
    "$$I_{xx} = \\int_{D} \\ y^2 \\ dx \\ dy = \\int_{\\partial D} -\\dfrac{1}{3}y^3 \\ dx $$\n",
    "$$\\boxed{I_{xx}\\approx -\\sum_{i=0}^{n} \\dfrac{1}{12} \\left(x_{i+1}-x_{i}\\right)\\cdot \\left(y_{i}^3+y_i^2 y_{i+1}+y_i y_{i+1}^2 + y_{i+1}^3\\right)}$$\n",
    "\n",
    "* With $Q=\\dfrac{1}{4}x^2y$ and $P=-\\dfrac{1}{4}xy^2$\n",
    "$$I_{xy} = \\int_{D} \\ xy \\ dx \\ dy = \\dfrac{1}{4}\\int_{\\partial D} -xy^2 \\ dx + x^2 y \\ dy $$\n",
    "$$\\boxed{I_{xy}\\approx \\sum_{i=0}^{n} \\dfrac{1}{24} \\left(x_i \\cdot y_{i+1}-x_{i+1} \\cdot y_{i+1}\\right)\\cdot \\left(2x_iy_i + x_{i}y_{i+1}+x_{i+1}y_{i}+2x_{i+1}y_{i+1}\\right)}$$\n",
    "\n",
    "* With $Q=\\dfrac{1}{3}x^3$ and $P=0$\n",
    "$$I_{yy} = \\int_{D} \\ x^2 \\ dx \\ dy = \\int_{\\partial D} -\\dfrac{1}{3}x^3 \\ dy$$\n",
    "$$\\boxed{I_{yy} \\approx \\sum_{i=0}^{n} \\dfrac{1}{12}\\left(x_{i}^3+x_i^2 x_{i+1}+x_i x_{i+1}^2 + x_{i+1}^3\\right)\\cdot \\left(y_{i+1}-y_{i}\\right)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "perimeter: 1.570e+02\n",
      "Cross-Sectional Area:\n",
      "    A: 1.960e+03\n",
      "First moment of area:\n",
      "    Qx: -7.350e-14\n",
      "    Qx: -1.819e-12\n",
      "Centroids:\n",
      "    xc: -9.279e-16\n",
      "    yc: -3.749e-17\n",
      "Second moment of inertia:\n",
      "    Global:\n",
      "        xx: 3.058e+05\n",
      "        xy: -1.268e-11\n",
      "        yy: 3.058e+05\n",
      "    Center:\n",
      "        xx: 3.058e+05\n",
      "        xy: -1.268e-11\n",
      "        yy: 3.058e+05\n",
      "Radii of Gyration:\n",
      "    x: 1.249e+01\n",
      "    y: 1.249e+01\n",
      "Elastic Section Moduli\n",
      "    Zxx+: 1.223e+04\n",
      "    Zxx-: 1.223e+04\n",
      "    Zyy+: 1.223e+04\n",
      "    Zyy-: 1.223e+04\n",
      "Principal axes:\n",
      "    Global:\n",
      "        I11: 3.058e+05\n",
      "        I22: 3.058e+05\n",
      "        angle: 0.0 degrees\n",
      "    Center:\n",
      "        I11: 3.058e+05\n",
      "        I22: 3.058e+05\n",
      "        angle: 0.0 degrees\n"
     ]
    }
   ],
   "source": [
    "perimeter = 0\n",
    "area = 0\n",
    "Q = {\"x\": 0, \"y\": 0}\n",
    "Ig = {\"xx\": 0, \"xy\": 0, \"yy\": 0}\n",
    "Ic = {\"xx\": 0, \"xy\": 0, \"yy\": 0}\n",
    "center = {\"x\": 0, \"y\": 0}\n",
    "radius = {\"x\": 0, \"y\": 0}\n",
    "Z = {\"xx+\": 0, \"xx-\": 0, \"yy+\": 0, \"yy-\": 0}\n",
    "\n",
    "for i in range(npts):\n",
    "    x0, x1 = x[i], x[i+1]\n",
    "    y0, y1 = y[i], y[i+1]\n",
    "    dx, dy = x1-x0, y1-y0\n",
    "    perimeter += np.sqrt(dx**2 + dy**2)\n",
    "    area += (x0+x1)*dy/2\n",
    "    # area += dx*(y[i]+y[i+1])/2 # The same\n",
    "    Q[\"x\"] -= dx*(y0**2+y0*y1+y1**2)/6\n",
    "    Q[\"y\"] += (x0**2+x0*x1+x1**2)*dy/6\n",
    "    Ig[\"xx\"] -= dx*(y0**3+y0**2*y1+y0*y1**2+y1**3)/12\n",
    "    Ig[\"xy\"] += (x0*y1 - x1*y0)*(2*x0*y0 + x1*y0 + x0*y1 + 2*x1*y1)/24\n",
    "    Ig[\"yy\"] += (x0**3+x0**2*x1+x0*x1**2+x1**3)*dy/12\n",
    "center[\"x\"] = Q[\"y\"]/area\n",
    "center[\"y\"] = Q[\"x\"]/area\n",
    "Ic[\"xx\"] = Ig[\"xx\"]-area*center[\"y\"]**2 \n",
    "Ic[\"yy\"] = Ig[\"yy\"]-area*center[\"x\"]**2 \n",
    "Ic[\"xy\"] = Ig[\"xy\"]-area*center[\"x\"]*center[\"y\"]\n",
    "radius[\"x\"] = np.sqrt(Ig[\"xx\"]/area)\n",
    "radius[\"y\"] = np.sqrt(Ig[\"yy\"]/area)\n",
    "Z[\"xx+\"] = Ic[\"xx\"]/(max(y) - center[\"y\"])\n",
    "Z[\"xx-\"] = Ic[\"xx\"]/(center[\"y\"]-min(y))\n",
    "Z[\"yy+\"] = Ic[\"yy\"]/(max(x) - center[\"x\"])\n",
    "Z[\"yy-\"] = Ic[\"yy\"]/(center[\"x\"]-min(x))\n",
    "matrix_global = np.array([[Ig[\"xx\"], Ig[\"xy\"]],\n",
    "                          [Ig[\"xy\"], Ig[\"yy\"]]])\n",
    "matrix_center = np.array([[Ic[\"xx\"], Ic[\"xy\"]],\n",
    "                          [Ic[\"xy\"], Ic[\"yy\"]]])\n",
    "eigs_global, eigvecg = np.linalg.eigh(matrix_global)\n",
    "eigs_center, eigvecc = np.linalg.eigh(matrix_center)\n",
    "angle_phi_global = 180*np.arccos(np.inner((1, 0), eigvecg[0]))/np.pi-90\n",
    "angle_phi_center = 180*np.arccos(np.inner((1, 0), eigvecg[0]))/np.pi-90\n",
    "\n",
    "print(f\"perimeter: {perimeter:.3e}\")\n",
    "print(f\"Cross-Sectional Area:\")\n",
    "print(f\"    A: {area:.3e}\")\n",
    "print(\"First moment of area:\")\n",
    "print(f\"    Qx: {Q['x']:.3e}\")\n",
    "print(f\"    Qx: {Q['y']:.3e}\")\n",
    "print(f\"Centroids:\")\n",
    "print(f\"    xc: {center['x']:.3e}\")\n",
    "print(f\"    yc: {center['y']:.3e}\")\n",
    "print(f\"Second moment of inertia:\")\n",
    "print(f\"    Global:\")\n",
    "print(f\"        xx: {Ig['xx']:.3e}\")\n",
    "print(f\"        xy: {Ig['xy']:.3e}\")\n",
    "print(f\"        yy: {Ig['yy']:.3e}\")\n",
    "print(f\"    Center:\")\n",
    "print(f\"        xx: {Ic['xx']:.3e}\")\n",
    "print(f\"        xy: {Ic['xy']:.3e}\")\n",
    "print(f\"        yy: {Ic['yy']:.3e}\")\n",
    "print(f\"Radii of Gyration:\")\n",
    "print(f\"    x: {radius['x']:.3e}\")\n",
    "print(f\"    y: {radius['y']:.3e}\")\n",
    "print(f\"Elastic Section Moduli\")\n",
    "print(f\"    Zxx+: {Z['xx+']:.3e}\")\n",
    "print(f\"    Zxx-: {Z['xx-']:.3e}\")\n",
    "print(f\"    Zyy+: {Z['yy+']:.3e}\")\n",
    "print(f\"    Zyy-: {Z['yy-']:.3e}\")\n",
    "print(f\"Principal axes:\")\n",
    "print(f\"    Global:\")\n",
    "print(f\"        I11: {eigs_global[0]:.3e}\")\n",
    "print(f\"        I22: {eigs_global[1]:.3e}\")\n",
    "print(f\"        angle: {angle_phi_global} degrees\")\n",
    "print(f\"    Center:\")\n",
    "print(f\"        I11: {eigs_center[0]:.3e}\")\n",
    "print(f\"        I22: {eigs_center[1]:.3e}\")\n",
    "print(f\"        angle: {angle_phi_center} degrees\")"
   ]
  }
 ],
 "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.10.11"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}