dynamicForrestal.ipynb

{
 "metadata": {
  "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.7.9"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python379jvsc74a57bd0fd9ba24f3a295a5eb35c3d547b8775e76ddd3547188cc1cd81b2309f957eca06",
   "display_name": "Python 3.7.9 64-bit ('base': conda)"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sigmarr,sigmaqq,r,u,t,Y,rhot = sympy.symbols(r\"sigma_rr,sigma_\\theta\\theta,r,u,t,Y,rho_t\")\n",
    "a,c,E = sympy.symbols(\"a,c,E\")\n",
    "a = sympy.Function(\"a\")(t)\n",
    "\n",
    "adot = sympy.symbols(r\"\\dot{a}\")\n",
    "addot = sympy.symbols(r\"\\ddot{a}\")\n",
    "adot\n",
    "sigmarrf = sympy.Function(\"sigma_rr\")(r)"
   ]
  },
  {
   "source": [
    "前回力の釣り合いの式は以下に示すように、右辺が0となっていました。\n",
    "\n",
    "$$ \\frac{d\\sigma_{rr}}{dr}+\\frac{2}{r}\\left(\\sigma_{rr}-\\sigma_{\\theta\\theta}\\right) = 0$$\n",
    "\n",
    "上式は力の釣り合いの式なわけですが、標的の各部分についての運動方程式とも言えます。\n",
    "力が釣り合っていて、右辺が0ということは、侵徹中にも関わらず標的の各部分は静止をしていることを意味しており、\n",
    "その意味で、前回の結果は静的cavity expansion analysisとも呼ばれます。\n",
    "\n",
    "では、力が釣り合わず、標的の各部分が運動する状況を考えると、その$\\sigma_{rr}$はどうなるでしょうか？というのが今回の話題です。\n",
    "まず、力の釣り合いの式を、見慣れた運動方程式の形にすると、\n",
    "$$ \\frac{d\\sigma_{rr}}{dr}+\\frac{2}{r}\\left(\\sigma_{rr}-\\sigma_{\\theta\\theta}\\right)  =- \\rho_t \\left( \\frac{\\partial v}{\\partial t} +v\\frac{\\partial v}{\\partial r} \\right)$$\n",
    "となります。ここで、右辺が妙な形になっているのは、オイラー座標系での時間微分になっているためです。$v$は各部分の粒子速度であり、$r, u$と同じ向きで外側が正です。また、$v$は$u$の時間微分ですので、同じような手続きを取ると、\n",
    "\n",
    "$$ v = \\left( \\frac{\\partial u}{\\partial t} +v\\frac{\\partial u}{\\partial r} \\right) $$\n",
    "$$ \\frac{\\partial u}{\\partial t} = v \\left(1 - \\frac{\\partial u}{\\partial r} \\right)$$\n",
    "\n",
    "として$u$と$v$の関係が得られます。\n",
    "\n",
    "今、前回の記事から$u$の形は以下の通り\n",
    "$$u=r\\left[1-\\left(1-\\frac{a^3}{r^3}\\right)^{1/3}\\right]$$\n",
    "として与えられているので、$a$を\n",
    "$$a = a(t)$$\n",
    "なる関数だと思って微分して\n",
    "$$ \\frac{\\partial u}{\\partial t} = v \\left(1 - \\frac{\\partial u}{\\partial r} \\right)$$\n",
    "に代入すれば、"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "\n",
    "u = r*(1-(1-(a/r)**sympy.Rational(3,1))**(sympy.Rational(1,3)))\n",
    "dudt = sympy.diff(u,t)\n",
    "dudr = sympy.diff(u,r)\n",
    "v = (dudt/(1-dudr)).simplify()\n"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": 2,
   "outputs": []
  },
  {
   "source": [
    "$$v = \\frac{\\dot{a} a^{2}{\\left(t \\right)}}{r^{2}}$$\n",
    "を得ます。ここで、$\\dot{a}$は$a$の時間微分です。\n",
    "この結果を用いて、前回と同じくまずは降伏領域($a\\geq r \\geq c$)の$\\sigma_{rr}$を求めましょう。\n",
    "本項第2式に$\\sigma_{rr}-\\sigma_{\\theta\\theta} = Y$と$v$を代入して、"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "v = (dudt/(1-dudr)).simplify()\n",
    "dvdr = sympy.diff(v,r)\n",
    "dvdt = sympy.diff(v,t)\n",
    "rh = sympy.diff(sigmarrf,r)\n",
    "lh = (-2*Y/r-rhot*(dvdt+v*dvdr))\n",
    "\n",
    "sympy.Eq(rh,lh.subs([[sympy.diff(sympy.diff(a,t),t),addot]]).subs(sympy.diff(a,t), adot))\n"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": 3,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "Eq(Derivative(sigma_rr(r), r), -2*Y/r - rho_t*(\\ddot{a}*a(t)**2/r**2 + 2*\\dot{a}**2*a(t)/r**2 - 2*\\dot{a}**2*a(t)**4/r**5))"
      ],
      "text/latex": "$\\displaystyle \\frac{d}{d r} \\sigma_{rr}{\\left(r \\right)} = - \\frac{2 Y}{r} - \\rho_{t} \\left(\\frac{\\ddot{a} a^{2}{\\left(t \\right)}}{r^{2}} + \\frac{2 \\dot{a}^{2} a{\\left(t \\right)}}{r^{2}} - \\frac{2 \\dot{a}^{2} a^{4}{\\left(t \\right)}}{r^{5}}\\right)$"
     },
     "metadata": {},
     "execution_count": 3
    }
   ]
  },
  {
   "source": [
    "という形で得られます。ここで$\\ddot{a}$は$a$の時間に対する2階微分です。\n",
    "前回の記事と同じように積分してやれば、\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "source": [
    "Irh = sympy.integrate(sympy.diff(sigmarrf,r),(r,a,c)) \n",
    "Ilh = sympy.integrate(lh,(r,a,c)).simplify()\n",
    "Ilh.simplify()\n",
    "sigmarr_plast = sympy.Function(\"sigma_rr\")(c)-Ilh\n",
    "sympy.Eq(sympy.Function(\"sigma_rr\")(a), sigmarr_plast.simplify())\n"
   ],
   "cell_type": "code",
   "metadata": {},
   "execution_count": 4,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "Eq(sigma_rr(a(t)), 2*Y*log(c) - 2*Y*log(a(t)) + rho_t*a(t)*Derivative(a(t), (t, 2)) + 3*rho_t*Derivative(a(t), t)**2/2 + sigma_rr(c) - rho_t*a(t)**2*Derivative(a(t), (t, 2))/c - 2*rho_t*a(t)*Derivative(a(t), t)**2/c + rho_t*a(t)**4*Derivative(a(t), t)**2/(2*c**4))"
      ],
      "text/latex": "$\\displaystyle \\sigma_{rr}{\\left(a{\\left(t \\right)} \\right)} = 2 Y \\log{\\left(c \\right)} - 2 Y \\log{\\left(a{\\left(t \\right)} \\right)} + \\rho_{t} a{\\left(t \\right)} \\frac{d^{2}}{d t^{2}} a{\\left(t \\right)} + \\frac{3 \\rho_{t} \\left(\\frac{d}{d t} a{\\left(t \\right)}\\right)^{2}}{2} + \\sigma_{rr}{\\left(c \\right)} - \\frac{\\rho_{t} a^{2}{\\left(t \\right)} \\frac{d^{2}}{d t^{2}} a{\\left(t \\right)}}{c} - \\frac{2 \\rho_{t} a{\\left(t \\right)} \\left(\\frac{d}{d t} a{\\left(t \\right)}\\right)^{2}}{c} + \\frac{\\rho_{t} a^{4}{\\left(t \\right)} \\left(\\frac{d}{d t} a{\\left(t \\right)}\\right)^{2}}{2 c^{4}}$"
     },
     "metadata": {},
     "execution_count": 4
    }
   ]
  },
  {
   "source": [
    "のように、$c$での$\\sigma_{rr}$を得ます。とはいえ、上の形では見にくいので、整理すると、\n",
    "\n",
    "$$\\sigma_{rr}(a) = \\sigma_{rr}(c)-2Y\\ln{\\frac{a}{c}} +\\rho_t\\left[ -(\\ddot{a}a^2+2a\\dot{a}^2)\\left(\\frac{1}{c} - \\frac{1}{a}\\right) + \\frac{a^4\\dot{a}^2}{2}\\left(\\frac{1}{c^4}-\\frac{1}{a^4}\\right) \\right]$$\n",
    "\n",
    "となります。それでは、$\\sigma_{rr}(c)$を求めるために、前記事と同じく弾性変形領域($c\\leq r$)を解きましょう。\n",
    "やることは大して変わらなくて、結局、\n",
    "$$ \\frac{d\\sigma_{rr}}{dr}+\\frac{2}{r}\\left(\\sigma_{rr}-\\sigma_{\\theta\\theta}\\right)  =- \\rho_t \\left( \\frac{\\partial v}{\\partial t} +v\\frac{\\partial v}{\\partial r} \\right)$$\n",
    "に\n",
    "$$\\sigma_{rr}-\\sigma_{\\theta\\theta} = \\frac{4Ea^3}{3r^4}$$\n",
    "を代入して、"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "Eq(sigma_rr(a), 4*E*a**3/(9*c**3) - 2*Y*log(a) + 2*Y*log(c) + rho_t*(\\ddot{a}*a + 3*\\dot{a}**2/2))"
      ],
      "text/latex": "$\\displaystyle \\sigma_{rr}{\\left(a \\right)} = \\frac{4 E a^{3}}{9 c^{3}} - 2 Y \\log{\\left(a \\right)} + 2 Y \\log{\\left(c \\right)} + \\rho_{t} \\left(\\ddot{a} a + \\frac{3 \\dot{a}^{2}}{2}\\right)$"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "source": [
    "v = (dudt/(1-dudr)).simplify()\n",
    "dvdr = sympy.diff(v,r)\n",
    "dvdt = sympy.diff(v,t)\n",
    "lh = (-4*E*a**3/r**4/3-rhot*(dvdt+v*dvdr))\n",
    "rh = - sympy.Function(\"sigma_rr\")(c)\n",
    "sigmarr_elast = sympy.integrate(lh, (r,c))\n",
    "sigmarr_a = sigmarr_plast.subs(sympy.Function(\"sigma_rr\")(c), sigmarr_elast).simplify().collect(rhot)\n",
    "sympy.Eq(sympy.Function(\"sigma_rr\")(a), sigmarr_a)\n",
    "sympy.Eq(sympy.Function(\"sigma_rr\")(a), sigmarr_a).subs(sympy.diff(sympy.diff(a,t),t), addot).subs(sympy.diff(a,t), adot).subs(a, sympy.symbols(\"a\"))"
   ]
  },
  {
   "source": [
    "を得ます。整理すると、\n",
    "$$  \\sigma_{rr}(a) = \\frac{4Ea^3}{9c^3} -2Y\\ln{a/c}+ \\rho_t\\left(\\ddot{a}a+\\frac{3\\dot{a}^2}{2}\\right)$$\n",
    "となります。前回示したように、\n",
    "$$ \\left(\\frac{a}{c}\\right)^3 =  \\frac{3Y}{2E}$$\n",
    "ですから、\n",
    "\n",
    "$$ \\sigma_{rr}(a) =\\frac{2Y}{3}\\left[ 1+\\ln{\\frac{2E}{3Y}}\\right] + \\rho_t\\left(\\ddot{a}a+\\frac{3\\dot{a}^2}{2}\\right) $$\n",
    "\n",
    "となります。また、$\\ddot{a}$は極めて小さいこと、$\\dot{a}$、すなわちCavityの拡大速度は侵徹速度$V_z$のCavityの$r$方向成分なので$\\dot{a}=V_z\\cos{\\theta}$となることを用いて\n",
    "\n",
    "$$ \\sigma_{rr}(a) =\\frac{2Y}{3}\\left[ 1+\\ln{\\frac{2E}{3Y}}\\right] + \\rho_t\\left(\\frac{3(V_z\\cos{\\theta})^2}{2}\\right) $$\n",
    "\n",
    "として、動的なcavity expansion analysisでの$\\sigma_{rr}$を得ます。"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}