Skip to content

Instantly share code, notes, and snippets.

@carlos-adir
Created June 25, 2023 09:30
Show Gist options
  • Save carlos-adir/67478fa6650ffda1c8e63b7f0cb5d695 to your computer and use it in GitHub Desktop.
Save carlos-adir/67478fa6650ffda1c8e63b7f0cb5d695 to your computer and use it in GitHub Desktop.
Compute section properties of a circle by using boundary method
Display the source blob
Display the rendered blob
Raw
{
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment