Rayleigh Ritz method

rayleigh-ritz-method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this post, we will look at a way of approximating deflection function of beams using Rayleigh Ritz method. I had been recently watching [this (click here for youtube link) lecture series by Clayton Pettit](https://www.youtube.com/watch?v=jBNi-iSKXhw&list=PLDkz31TxAOcs1HNWLFbuwmAWwFyMQnZVp&index=3) which had been a very useful introduction to finite element methods. This is one of the assignment problems provided in that series. \n",
    "\n",
    "Regarding the choice of language here, it is Julia but also python (through PyCall, but for sympy that's taken care by the `SymPy` package). The reason for this acrobatics is that python had issues converting singularity functions which tend to rise when analytically solving the system using Singularity method (if that sounds like gibberish, I would suggest \"Machine Design, an Integrated Approach by Norton\" for a quick brush up). I did initially intend to use Julia but the symbolics library doesn't have all the calculus functionalities yet. Now, let's get started!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SymPy\n",
    "using Plots\n",
    "using LaTeXStrings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@vars x1 x2;\n",
    "@vars a0 a1 a2 a3 a4 a5 a6;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following cantilever beam.\n",
    "\n",
    "![cantilever](./problem.png)\n",
    "\n",
    "The value for various parameters in the problem are defined as following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = 20e6 # Young's Modulus, KPa\n",
    "L = 6 # Length, m\n",
    "q = -45 # Uniformly Distributed Load, KNm\n",
    "P = -100 # Point Load, KN\n",
    "d = 4 # Position of UDL, m\n",
    "b = 0.3 # Width of beam, m\n",
    "h = 0.5 # Height of the beam,  m\n",
    "Ig = b*h^3/12 # Moment of Inertia, m^4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider an Euler-Bernoulli beam (with the assumptions). The deflection of the beam can be obtained by solving the loading equation,\n",
    "\n",
    "$$\\displaystyle{q= EI\\left( \\frac{d^2y}{dx_1^2} \\right) }$$\n",
    "\n",
    "with the boundary conditions, \n",
    "\n",
    "$$y(x_1=0) = 0$$\n",
    "$$\\frac{dy}{dx}(x_1=0) = 0$$\n",
    "\n",
    "Further, there are the implicit boundary conditions, $\\displaystyle{V(0^-)=0, M(0^-)=0, V(l^+)=0, M(l^+)=0 }$ which are utilised to solve the equation analytically.\n",
    "\n",
    "If the above equation is solved, the shear force $V$ and bending moment $M$ can be obtained through,\n",
    "$$\\displaystyle{M = EI\\left( \\frac{d ^2y}{d x_1^2}  \\right) }$$\n",
    "$$\\displaystyle{V = EI\\left( \\frac{d ^3y}{d x_1^3}  \\right) }$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving the above equation becomes tedious as the complexity of the loading increases. Hence, one might often resort to approximations. Rayleigh-Ritz method is one such method of approximating the deflection equation. This can be broken down into the following steps.\n",
    "\n",
    "1. Approximate the function to be solved for as a polynomial\n",
    "2. Apply boundary conditions\n",
    "3. Find the potential energy with this equation and minimize it by taking variations with respect to the parameters. \n",
    "4. Solve the arising equations to find the constants. Substitute them back to get the approximate function.\n",
    "\n",
    "Let's start with the approximation of the deflection as a polynomial. Note that, the loading equation is a 4th order equation which means that our polynomial should have a minimum degree of 4. Here, we will take a polynomial of order 6 and see how it turns out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_approx = a0 + a1*x1 + a2*x1^2 + a3*x1^3 + \n",
    "           a4*x1^4 + a5*x1^5 + a6*x1^6 # polynomial of degree 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to apply the boundary conditions on the $y_\\text{approx}$. This is done below and the following values for the first two coefficients are obtained.\n",
    "$$\\begin{aligned}\n",
    "a_0 &= 0 \\\\\n",
    " a_1 &= 0\n",
    " \\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BC1 = y_approx(x1=>0)\n",
    "BC2 = diff(y_approx, x1)(x1=>0)\n",
    "BC1, BC2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BC_sol = solve([BC1,BC2], [a0,a1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substituting the boundary conditions, $y_\\text{approx}$ is reduced to a polynomial with 5 coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_approx = y_approx(a0=>BC_sol[a0], a1=>BC_sol[a1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to find the potential energy which is given as,\n",
    "$$\\text{Potential Energy = Internal Strain Energy - External Work}$$\n",
    "$$\n",
    " \\begin{aligned}\n",
    "\t \\prod(y) &= U(y) - W(y) \\\\\n",
    "\t U &= \\int_0^L \\frac{EI}{2}\\left( \\frac{d ^2y_{\\text{approx}}}{d x_1^2}  \\right) ^2 dx_1^2 \\\\\n",
    "\t W &= \\int_0^L q y_{\\text{approx}}dx_1 + \\sum_{i=1}^{n} P_i y_{\\text{approx}, i} + \\sum_{j=1}^{m} M_j \\theta_{\\text{approx},j}\n",
    " \\end{aligned}\n",
    " $$\n",
    "\n",
    " In the above equation, the external work is given as,\n",
    "\n",
    " $$\\text{Work = W.D by UDL + W.D by Point Load + W.D by concentrated moments}$$\n",
    "\n",
    " where, the concentrated moment doesn't exist in our case. Hence, ignoring that, the rest is implemented as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Internal Strain Energy\n",
    "U = integrate((E*Ig/2 * diff(y_approx, x1, 2)^2), (x1, 0, L))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# External work done\n",
    "W1 = integrate(q*y_approx, (x1, 0, L))\n",
    "W2 = P*y_approx.subs(x1, d)\n",
    "W = W1 + W2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "∏ = U - W"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, $\\prod(y(x_1))$ is a functional which means that it takes a function $y_{\\text{approx}}(x_1)$ as input and gives a scalar as output. Our aim is to find the $y_{\\text{approx}}$ which gives the minimum value of the potential energy. In order to find that, we differentiate $\\prod$ with respect to each of the remaining coefficients ($\\displaystyle{a_2, a_3, a_4, a_5, a_6}$). We solve the resulting system of equations algebraically to get the corresponding values. Note that, if we had not applied the boundary conditions initially, the coefficients $\\displaystyle{a_0, a_1}$ can also be solved this way but that won't gaurantee that the resulting equation will satisfy the boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq1 = diff(∏, a2)\n",
    "eq2 = diff(∏, a3)\n",
    "eq3 = diff(∏, a4)\n",
    "eq4 = diff(∏, a5)\n",
    "eq5 = diff(∏, a6)\n",
    "eq_sol = solve([eq1, eq2, eq3, eq4, eq5], [a2, a3, a4, a5, a6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substituting these values back in the equation, we obtain the approximated equation for the deflection. Further, we substitute this equation in the equations for $M$ and $V$ to get the $\\displaystyle{M_{\\text{approx}}}$ and $\\displaystyle{V_{\\text{approx}}}$ respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_approx = y_approx(a2=>eq_sol[a2], a3=>eq_sol[a3], a4=>eq_sol[a4], a5=>eq_sol[a5], a6=>eq_sol[a6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_approx = E*Ig*diff(y_approx, x1, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V_approx = E*Ig*diff(y_approx, x1, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to check the quality of approximation, we will compare it against the analytical solution which is given by,\n",
    "$$\\displaystyle{y_{\\text{act}}= \\frac{1}{EI}\\left( \\frac{1}{2}M_1\\left<x_1-0 \\right>^2 + \\frac{1}{6}R_1\\left<x_1-0 \\right>^3 + \\frac{1}{24}qx_1^4 + \\frac{1}{6}P\\left<x_1-4 \\right>^3 \\right) }$$\n",
    "\n",
    "where, $<\\cdot>$ indicate the singularity functions, $R_1$ and $M_1$ are the reaction load and moment on the fixed end. The values are substituted and the $y_{\\text{act}}$ is derived below.\n",
    "\n",
    "The corresponding $M_{\\text{act}}$ and $V_{\\text{act}}$ are also calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R1 = -q*L - P\n",
    "M1 = P*d + q*L^2/2\n",
    "\n",
    "y_act = 1/(E*Ig)*(M1*x1^2/2 + R1*x1^3/6 + q*x1^4/24 + P*sympy.SingularityFunction(x1, 4, 3)/6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M_act = E*Ig*diff(y_act, x1, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V_act = E*Ig*diff(y_act, x1, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the plot below, we are plotting the actual and approximated y values. It can be seen that approximation conforms very well to the original function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "plot(x->y_approx.evalf(subs=Dict(x1=>x)), 0, L, label=\"y approximate\", lw=2,\n",
    "     lc=\"black\", title=\"Deflection\", xlabel=L\"x_1\", ylabel=\"y in m\", \n",
    "     size=(600, 200), ylim=(-0.2, 0.2), axes_style=:origin)\n",
    "plot!(x->y_act.evalf(subs=Dict(x1=>x)), 0, L, label=\"y actual\", lw=3,\n",
    "      ls=:dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the bending diagram, the approximation still seems to hold good though there seems to be a something peculiar at $x_1=4$ but not enough to judge without careful examination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(x->M_approx.evalf(subs=Dict(x1=>x)), 0, L, label=\"M approximate\", lw=2,\n",
    "     lc=\"black\", title=\"Bending Moment\", xlabel=L\"x_1\", ylabel=\"M in KNm\", \n",
    "     size=(600, 200), ylim=(-1500, 1000), axes_style=:origin)\n",
    "plot!(x->M_act.evalf(subs=Dict(x1=>x)), 0, L, label=\"M actual\", lw=3,\n",
    "      ls=:dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, while drawing the Shear force diagram, the caveats of approximation is seen. The area under the curve remains the same (aka the bending moment) but the actual curve traced is different. The analytical solution has the discontinuity at $x_1=4$ mark which is smoothed out. Also, the maximum shear force predicted is slightly lower than the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(x->V_approx.evalf(subs=Dict(x1=>x)), 0, L, label=\"V approximate\", lw=2,\n",
    "     lc=\"black\", title=\"Shear Force\", xlabel=L\"x_1\", ylabel=\"V in KN\", \n",
    "     size=(600, 200), ylim=(-100,500))\n",
    "plot!(x->V_act.evalf(subs=Dict(x1=>x)), 0, L, label=\"V actual\", lw=3,\n",
    "      ls=:dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we shall take a look at the bending and shear stresses. For an Euler-Bernoulli beam, only $\\sigma_{11}$ and $\\sigma_{12}$ exists which are given by,\n",
    "\n",
    "$$\\displaystyle{\\sigma_{11} = \\frac{Mx_2}{I}}$$\n",
    "\n",
    "$$\\displaystyle{\\sigma_{12}= \\frac{VQ}{Ib}}$$\n",
    "\n",
    "where $\\displaystyle{Q = b \\left( \\frac{t}{2}-x_2 \\right) \\left( \\frac{t}{4}+\\frac{x_2}{2} \\right) }$ for a rectangular cross sectioned beam."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "σ_11_act = -M_act*x2/Ig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "σ_11_approx = -M_approx*x2/Ig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = b*(h/2 - x2)*(h/4 + x2/2)\n",
    "σ_12_act = V_act*Q/(Ig*b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "σ_12_approx = V_approx*Q/(Ig*b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following cells, the $\\sigma_{11}$ and $\\sigma_{12}$ actual and approximate values are plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1=0:0.01:L\n",
    "X2=-h/2:0.01:h/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = contourf(X1, X2, (X1, X2)->σ_11_act.evalf(subs=Dict(x1=>X1, x2=>X2)),\n",
    "         size=(600, 200), title=\"Bending Stress, σ_11 approximate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p2 = contourf(X1, X2, (X1, X2)->σ_11_approx.evalf(subs=Dict(x1=>X1, x2=>X2)),\n",
    "         size=(600, 200), title=\"Bending Stress, σ_11 actual\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3 = contourf(X1, X2, (X1, X2)->σ_12_act.evalf(subs=Dict(x1=>X1, x2=>X2)),\n",
    "size=(600, 200), title=\"Shear Stress, σ_12 approximate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p4 = contourf(X1, X2, (X1, X2)->σ_12_act.evalf(subs=Dict(x1=>X1, x2=>X2)),\n",
    "         size=(600, 200), title=\"Shear Stress, σ_12 actual\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l =@layout[grid(2,1) a{0.05w}]\n",
    "plot(p1, p2, layout=grid(2,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(p3, p4, layout=grid(2,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the comparison above, it can be seen that the stress distributions for both $y_\\text{act}$ and $y_\\text{approx}$ remain the same. Hence, this can be considered as a viable approximation for the deflection function."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.1",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}