Flamant.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flamant problem (Vertical point loading on half space)\n",
    "\n",
    "Consider an infinite half-space in $x \\geq 0$. A point load, $P$, acts vertically at the origin. We wish to find the stress distribution in the half-space. This problem is of a special kind, and we need to approach it a bit differently. First, we note that if we consider a semi-circle around the origin, then the following \"force-balance\" must hold:\n",
    "\n",
    "\\begin{gather}\n",
    "\\int_{-\\pi/2}^{\\pi/2} \\sigma_{rr} \\cos\\theta \\; r {\\rm d} \\theta = - P\n",
    "\\end{gather}\n",
    "\n",
    "Now, note that in the RHS we have a constant value and in the LHS we have an integration with respect to $\\theta$ while in the integrand, we have apparently an $r$-dependence. Now, the only possibile way in which this $r$-dependence can disappear is if $\\sigma_{rr}$ is proportional to $\\displaystyle \\frac{1}{r}$. Thus, $\\sigma_{rr}$ must be of the form $\\displaystyle \\sigma_{rr} = \\frac{h(\\theta)}{r}$, where it is important to note that $h(\\theta)$ must be symmetric in $\\theta$ about $\\theta = 0$. We will use this form of $\\sigma_{rr}$ to motivate the form of the Airy stress function. \n",
    "\n",
    "Towards that end, we first remind ourselves of the expression of $\\sigma_{rr}$ in terms of a generic Airy stress function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from polarUtilities import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\frac{\\partial}{\\partial r} \\varphi{\\left(r,\\theta \\right)}}{r} + \\frac{\\frac{\\partial^{2}}{\\partial \\theta^{2}} \\varphi{\\left(r,\\theta \\right)}}{r^{2}}$"
      ],
      "text/plain": [
       "Derivative(varphi(r, theta), r)/r + Derivative(varphi(r, theta), (theta, 2))/r**2"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi = sym.Function('varphi')(r,theta)\n",
    "sigma_rr(phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, if our unknown Airy stress function is $\\varphi$, we must have:\n",
    "\n",
    "\\begin{gather*}\n",
    "\\frac{1}{r} \\frac{\\partial \\varphi}{\\partial r} + \\frac{1}{r^2} \\frac{\\partial^2 \\varphi}{\\partial \\theta^2} = \\frac{h(\\theta)}{r}, \n",
    "\\end{gather*}\n",
    "\n",
    "from which we obtain:\n",
    "\\begin{gather*}\n",
    "\\frac{\\partial \\varphi}{\\partial r} + \\frac{1}{r} \\frac{\\partial^2 \\varphi}{\\partial \\theta^2} = h(\\theta). \n",
    "\\end{gather*}\n",
    "\n",
    "From the above, we can deduce that $\\varphi$ must be of the form $\\varphi = r g(\\theta)$. We set it up in our code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\varphi = r g{\\left(\\theta \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = sym.Function('g')(theta)\n",
    "phi = r*g\n",
    "\n",
    "from IPython.display import Math, Latex\n",
    "display(Math(r'\\varphi = {}'.format(sym.latex(phi))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we take a look at the biharmonic of this $\\phi$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{g{\\left(\\theta \\right)} + 2 \\frac{d^{2}}{d \\theta^{2}} g{\\left(\\theta \\right)} + \\frac{d^{4}}{d \\theta^{4}} g{\\left(\\theta \\right)}}{r^{3}}$"
      ],
      "text/plain": [
       "(g(theta) + 2*Derivative(g(theta), (theta, 2)) + Derivative(g(theta), (theta, 4)))/r**3"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polarbiharmonic(phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the polar biharmonic equation becomes very easy to solve analytically. Nevertheless, we proceed with the solution through SymPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle g{\\left(\\theta \\right)} = \\left(C_{1} + C_{2} \\theta\\right) e^{- i \\theta} + \\left(C_{3} + C_{4} \\theta\\right) e^{i \\theta}$"
      ],
      "text/plain": [
       "Eq(g(theta), (C1 + C2*theta)*exp(-I*theta) + (C3 + C4*theta)*exp(I*theta))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.dsolve(polarbiharmonic(phi),g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We rewrite $g(\\theta)$ in a slightly different form as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle g(\\theta) = A \\cos{\\left(\\theta \\right)} + B \\theta \\cos{\\left(\\theta \\right)} + C \\sin{\\left(\\theta \\right)} + D \\theta \\sin{\\left(\\theta \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A, B, C, D = sym.symbols('A, B, C, D')\n",
    "gnew = A*sym.cos(theta) + B*theta*sym.cos(theta) + C*sym.sin(theta) + D*theta*sym.sin(theta)\n",
    "display(Math(r'g(\\theta) = {}'.format(sym.latex(gnew))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we find the expressions of $\\sigma_{rr}$, $\\sigma_{\\theta\\theta}$, and $\\sigma_{r\\theta}$ as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr} = - \\frac{2 B \\sin{\\left(\\theta \\right)}}{r} + \\frac{2 D \\cos{\\left(\\theta \\right)}}{r}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{\\theta\\theta} = 0$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{r\\theta} = 0$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "phi = r*gnew\n",
    "sigmarr = sigma_rr(phi)\n",
    "display(Math(r'\\sigma_{{rr}} = {}'.format(sym.latex(sigmarr))))\n",
    "\n",
    "sigmatt = sigma_tt(phi)\n",
    "display(Math(r'\\sigma_{{\\theta\\theta}} = {}'.format(sym.latex(sigmatt))))\n",
    "\n",
    "sigmart = sigma_rt(phi)\n",
    "display(Math(r'\\sigma_{{r\\theta}} = {}'.format(sym.latex(sigmart))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing the form of $\\sigma_{rr}$ just obtained with the previous general form $\\displaystyle \\sigma_{rr} = \\frac{h(\\theta)}{r}$, we find that $h(\\theta) = -2B\\sin\\theta + 2D\\cos\\theta$. But we had already noted that $h(\\theta)$ must be symmetric in $\\theta$ about $\\theta=0$. Therefore, we must have $B=0$. Consequently, the expression of $\\sigma_{rr}$ becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr} = \\frac{2 D \\cos{\\left(\\theta \\right)}}{r}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigmarr = sigmarr.subs(B,0)\n",
    "display(Math(r'\\sigma_{{rr}} = {}'.format(sym.latex(sigmarr))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unknown $D$ can be readily found by using the \"force-balance\" condition. Thus, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\pi D = - P$"
      ],
      "text/plain": [
       "Eq(pi*D, -P)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "P = sym.symbols('P')\n",
    "lhs = sym.integrate(r*sigmarr*sym.cos(theta),(theta,-sym.pi/2,sym.pi/2))\n",
    "rhs = -P\n",
    "eq = sym.Eq(lhs,rhs)\n",
    "display(eq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution for $D$ can be easily found now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle D = - \\frac{P}{\\pi}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Dsoln = -P/sym.pi\n",
    "display(Math(r'D = {}'.format(sym.latex(Dsoln))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can find the expression of $\\sigma_{rr}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sigma_{rr}^{\\rm soln} = - \\frac{2 P \\cos{\\left(\\theta \\right)}}{\\pi r}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigmarr_soln = sigmarr.subs(D,Dsoln)\n",
    "display(Math(r'\\sigma_{{rr}}^{{\\rm soln}} = {}'.format(sym.latex(sigmarr_soln))))"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.8.3"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}