Nonlinear_elasticity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonlinear elasticity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['font.size'] = 8\n",
    "figsize =(8,4)\n",
    "mpl.rcParams['figure.figsize'] = figsize\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import riemann_tools\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets\n",
    "from clawpack import riemann\n",
    "from exact_solvers import nonlinear_elasticity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we investigate a nonlinear model of elastic strain in heterogeneous materials.  This system is equivalent to the $p$-system of gas dynamics, although the stress-strain relation we will use here is very different from the pressure-density relation typically used in gas dynamics.  The equations we consider are:  \n",
    "\n",
    "\\begin{align}\n",
    "\\epsilon_t(x,t) - u_x(x,t) & = 0 \\\\\n",
    "(\\rho(x)u(x,t))_t - \\sigma(\\epsilon(x,t),x)_x & = 0.\n",
    "\\end{align}  \n",
    "\n",
    "Here $\\epsilon$ is the strain, $u$ is the velocity, $\\rho$ is the material density, $\\sigma$ is the stress,\n",
    "and ${\\mathcal M}=\\rho u$ is the momentum. \n",
    "The first equation is a kinematic relation, while the second represents Newton's second law.  This is a nonlinear \n",
    "conservation law with spatially varying flux function, in which  \n",
    "\n",
    "\\begin{align}\n",
    "q & = \\begin{bmatrix} \\epsilon \\\\ \\rho u \\end{bmatrix}, & f(q,x) & = \\begin{bmatrix} -{\\mathcal M}/\\rho(x) \\\\ -\\sigma(\\epsilon,x) \\end{bmatrix}.\n",
    "\\end{align}  \n",
    "\n",
    "If the stress-strain relationship is linear -- i.e. if $\\sigma(\\epsilon,x)=K(x)\\epsilon$ -- then this system is equivalent to the acoustics equations that we have\n",
    "studied previously.  Here we consider instead a quadratic stress response:\n",
    "\n",
    "\\begin{align}\n",
    "\\sigma(\\epsilon,x) = K_1(x) \\epsilon + K_2(x) \\epsilon^2.\n",
    "\\end{align}\n",
    "\n",
    "We assume that the spatially-varying functions $\\rho, K_1, K_2$ are piecewise constant, taking values\n",
    "$(\\rho_l, K_{1l}, K_{2l})$ for $x<0$ and values $(\\rho_r, K_{1r}, K_{2r})$ for $x>0$.  This system has been investigated numerically in <cite data-cite=\"leveque2002\"><a href=\"riemann.html#leveque2002\">(LeVeque, 2002)</a></cite>, <cite data-cite=\"leveque2003\"><a href=\"riemann.html#leveque2003\">(LeVeque & Yong, 2003)</a></cite>, and <cite data-cite=\"2012_ketchesonleveque_periodic\"><a href=\"riemann.html#2012_ketchesonleveque_periodic\">(Ketcheson & LeVeque, 2012)</a></cite>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if we take $\\rho=1$, $\\sigma=-p$, and $\\epsilon=v$, this system is equivalent to the p-system of Lagrangian gas dynamics\n",
    "\\begin{align*}\n",
    "v_t - u_x & = 0 \\\\\n",
    "u_t - p(v)_x & = 0,\n",
    "\\end{align*}\n",
    "in which $p$ represents pressure and $v$ represents specific volume."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hyperbolic structure\n",
    "The flux jacobian is\n",
    "\n",
    "\\begin{align}\n",
    "f'(q) = \\begin{bmatrix} 0 & -1/\\rho(x) \\\\ -\\sigma_\\epsilon(\\epsilon,x) & 0 \\end{bmatrix},\n",
    "\\end{align}  \n",
    "\n",
    "with eigenvalues (characteristic speeds)\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda^\\pm(x) = \\pm \\sqrt{\\frac{\\sigma_\\epsilon(\\epsilon,x)}{\\rho(x)}} = \\pm c(\\epsilon, x).\n",
    "\\end{align}\n",
    "\n",
    "Here for the stress-strain relation we have chosen, $\\sigma_\\epsilon = K_1(x) + 2 K_2(x)\\epsilon$.\n",
    "\n",
    "It's also convenient to define the nonlinear impedance $Z(\\epsilon, x) = \\rho(x) c(\\epsilon,x)$.  Then the eigenvectors of the flux Jacobian are\n",
    "\n",
    "\\begin{align}\n",
    "R & = \\begin{bmatrix} 1 & 1 \\\\ Z(\\epsilon,x) & -Z(\\epsilon,x) \\end{bmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "Both characteristic fields are genuinely nonlinear.  Furthermore, since the characteristic speeds each have a definite sign, this system does not admit transonic rarefactions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure of the Riemann solution\n",
    "Based on the eigenstructure of the flux jacobian above, the Riemann solution will always include a left-going and a right-going wave, each of which may be a shock or rarefaction (since both fields are genuinely nonlinear).  We will see -- similarly to our analysis in [the chapter on variable-speed-limit traffic](Traffic_variable_speed.ipynb) that the jump in $\\rho$ and $K$ at $x=0$ induces a stationary wave there.  See the figure below for the overall structure of the Riemann solution.\n",
    "\n",
    "![Structure of the Riemann solution for variable-coefficient nonlinear elasticity](./figures/nonlinear_elasticity_riemann.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hugoniot loci\n",
    "From the Rankine-Hugoniot jump conditions for the system we obtain the 1-Hugoniot locus for a state $(\\epsilon^*_l, u^*_l)$ connected by a 1-shock to a state $(\\epsilon_l, u_l)$:\n",
    "\n",
    "\\begin{align}\n",
    "u^*_l & = u_l - \\left( \\frac{\\left(\\sigma_l(\\epsilon^*_l)-\\sigma_l(\\epsilon_l)\\right)(\\epsilon^*_l-\\epsilon_l)}{\\rho_l}  \\right)^{1/2}\n",
    "\\end{align}\n",
    "\n",
    "Here $\\sigma_l(\\epsilon)$ is shorthand for the stress relation in the left medium.\n",
    "Similarly, a state $(\\epsilon^*_r,u^*_r)$ that is connected by a 2-shock to a state $(\\epsilon_r, u_r)$ must satisfy\n",
    "\n",
    "\\begin{align}\n",
    "u^*_r & = u_r - \\left( \\frac{\\left(\\sigma_r(\\epsilon^*_r)-\\sigma_r(\\epsilon_r)\\right)(\\epsilon^*_r-\\epsilon_r)}{\\rho_r}  \\right)^{1/2}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integral curves\n",
    "The integral curves can be found by writing $\\tilde{q}'(\\xi) = r^{1,2}(\\tilde{q}(\\xi))$ and integrating.  This leads to\n",
    "\\begin{align}\n",
    "u^*_l  & = u_l + \\frac{1}{3 K_{2l} \\sqrt{\\rho_l}} \\left( \\sigma_{l,\\epsilon}(\\epsilon^*_l)^{3/2} - \\sigma_{l,\\epsilon}(\\epsilon)^{3/2} \\right) \\label{NE:integral-curve-1} \\\\\n",
    "u^*_r  & = u_r - \\frac{1}{3 K_{2r} \\sqrt{\\rho_r}} \\left( \\sigma_{r,\\epsilon}(\\epsilon^*_r)^{3/2} - \\sigma_{r,\\epsilon}(\\epsilon)^{3/2} \\right)\\label{NE:integral-curve-2}\n",
    "\\end{align}\n",
    "Here $\\sigma_{l,\\epsilon}$ is the derivative of the stress function w.r.t $\\epsilon$ in the left medium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The entropy condition\n",
    "For a 1-shock, we need that $\\lambda^-(\\epsilon_l,x<0) > \\lambda^-(\\epsilon^*_l,x<0)$, which is equivalent to the condition $\\epsilon^*_l>\\epsilon_l$.  Similarly, a 2-shock is entropy-satisfying if $\\epsilon^*_r > \\epsilon_r$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interface conditions\n",
    "Because the flux depends explicitly on $x$, we do not necessarily have continuity of $q$ at $x=0$; i.e. in general $q^*_l \\ne q^*_r$.  Instead, the flux must be continuous: $f(q^*_l)=f(q^*_r)$.  For the present system, this means that the stress and velocity must be continuous:\n",
    "\n",
    "\\begin{align*}\n",
    "    u^*_l & = u^*_r \\\\\n",
    "    \\sigma(\\epsilon^*_l, K_{1l}, K_{2l}) & = \\sigma(\\epsilon^*_r, K_{1r}, K_{2r}).\n",
    "\\end{align*}\n",
    "\n",
    "This makes sense from a physical point of view: if the velocity were not continuous, the material would fracture (or overlap itself).  Continuity of the stress is required by Newton's laws."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure of rarefaction waves\n",
    "For this system, the structure of a centered rarefaction wave can be determined very simply.  Since the characteristic velocity must be equal to $\\xi = x/t$ at each point along the wave, we have $\\xi = \\pm\\sqrt{\\sigma_\\epsilon/\\rho}$, or\n",
    "\\begin{align}\n",
    "\\xi^2 = \\frac{K_1 + 2K_2\\epsilon}{\\rho}\n",
    "\\end{align}\n",
    "which leads to $\\epsilon = (\\rho\\xi^2 - K_1)/(2K_2)$.  Once the value of $\\epsilon$ is known, $u$ can be determined using the integral curves \\eqref{NE:integral-curve-1} or \\eqref{NE:integral-curve-2}."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution of the Riemann problem\n",
    "Below we show the solution of the Riemann problem.  *To view the code that computes this exact solution, uncomment and execute the next cell.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load exact_solvers/nonlinear_elasticity.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "def dsigma(eps, K1, K2):\n",
    "    \"Derivative of stress w.r.t. strain.\"\n",
    "    return K1 + 2*K2*eps\n",
    "\n",
    "def lambda1(q, xi, aux):\n",
    "    eps = q[0]\n",
    "    rho, K1, K2 = aux\n",
    "    return -np.sqrt(dsigma(eps, K1, K2)/rho)\n",
    "\n",
    "def lambda2(q, xi, aux):\n",
    "    return -lambda1(q,xi,aux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_plot_function(q_l, q_r, aux_l, aux_r):\n",
    "    states, speeds, reval, wave_types = \\\n",
    "        nonlinear_elasticity.exact_riemann_solution(q_l,q_r,aux_l,aux_r)\n",
    "        \n",
    "    def plot_function(t,which_char):\n",
    "        ax = riemann_tools.plot_riemann(states,speeds,reval,wave_types,\n",
    "                                        t=t,t_pointer=0,\n",
    "                                        extra_axes=True,\n",
    "                                        variable_names=['Strain','Momentum'])\n",
    "        if which_char == 1:\n",
    "            riemann_tools.plot_characteristics(reval,lambda1,(aux_l,aux_r),ax[0])\n",
    "        elif which_char == 2:\n",
    "            riemann_tools.plot_characteristics(reval,lambda2,(aux_l,aux_r),ax[0])\n",
    "        nonlinear_elasticity.phase_plane_plot(q_l, q_r, aux_l, aux_r, ax[3])\n",
    "        plt.show()\n",
    "    return plot_function        \n",
    "        \n",
    "def plot_riemann_nonlinear_elasticity(rho_l,rho_r,v_l,v_r):\n",
    "    plot_function = make_plot_function(rho_l,rho_r,v_l,v_r)\n",
    "    interact(plot_function, t=widgets.FloatSlider(value=0.,min=0,max=1.,step=0.1),\n",
    "             which_char=widgets.Dropdown(options=[None,1,2],\n",
    "                                         description='Show characteristics'));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aux_l = np.array((1., 5., 1.))\n",
    "aux_r = np.array((1., 2., 1.))\n",
    "q_l = np.array([2.1, 0.])\n",
    "q_r = np.array([0.0, 0.])\n",
    "\n",
    "plot_riemann_nonlinear_elasticity(q_l, q_r, aux_l, aux_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximate solution of the Riemann problem using $f$-waves\n",
    "The exact solver above requires a nonlinear iterative rootfinder and is relatively expensive.  A very cheap approximate Riemann solver for this system was developed in <cite data-cite=\"leveque2002\"><a href=\"riemann.html#leveque2002\">(LeVeque, 2002)</a></cite> using the $f$-wave approach.  One simply approximates both nonlinear waves as shocks, with speeds equal to the characteristic speeds of the left and right states:\n",
    "\n",
    "\\begin{align}\n",
    "s^1 & = - \\sqrt{\\frac{\\sigma_{\\epsilon,l}(\\epsilon_l)}{\\rho_l}} \\\\\n",
    "s^2 & = + \\sqrt{\\frac{\\sigma_{\\epsilon,r}(\\epsilon_r)}{\\rho_r}}.\n",
    "\\end{align}\n",
    "\n",
    "Meanwhile, the waves are obtained by decomposing the jump in the flux $f(q_r,x>0) - f(q_l,x<0)$ in terms of the eigenvectors of the flux jacobian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = riemann.nonlinear_elasticity_1D_py.nonlinear_elasticity_1D\n",
    "\n",
    "problem_data = {'stress_relation' : 'quadratic'}\n",
    "fw_states, fw_speeds, fw_reval = \\\n",
    "            riemann_tools.riemann_solution(solver,q_l,q_r,aux_l,aux_r,\n",
    "                                           problem_data=problem_data,\n",
    "                                           verbose=False,\n",
    "                                           stationary_wave=True,\n",
    "                                           fwave=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_function = \\\n",
    "    riemann_tools.make_plot_function(fw_states,fw_speeds, fw_reval,\n",
    "                                     layout='vertical', \n",
    "                                     variable_names=('Strain','Momentum'))\n",
    "interact(plot_function, t=widgets.FloatSlider(value=0.4,min=0,max=.9,step=.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of exact and approximate solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex_states, ex_speeds, ex_reval, wave_types = \\\n",
    "            nonlinear_elasticity.exact_riemann_solution(q_l,q_r,aux_l,aux_r)\n",
    "\n",
    "varnames = nonlinear_elasticity.conserved_variables\n",
    "plot_function = riemann_tools.make_plot_function([ex_states,fw_states],\n",
    "                                                 [ex_speeds,fw_speeds],\n",
    "                                                 [ex_reval,fw_reval],\n",
    "                                                 [wave_types,['contact']*3],\n",
    "                                                 ['Exact','$f$-wave'],\n",
    "                                                 layout='vertical',\n",
    "                                                 variable_names=varnames)\n",
    "interact(plot_function, t=widgets.FloatSlider(value=0.4,min=0, max=0.9, step=0.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, there are significant discrepancies between the approximate solution and the exact one.  But in a full numerical discretization, the left- and right-going waves are averaged over neighboring cells at the next step, and the approximate solver yields an effective result quite close to that of the exact solver."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}