nonlinear_elasticity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "75baa715",
   "metadata": {},
   "source": [
    "# Integration of nonlinear elasticity equations\n",
    "\n",
    "Soil is often modeled as having pressure-dependent shear and bulk moduli. Integrating these equations poses some challenges in return mapping algorithms. This notebook explores a common functional form for pressure dependence of shear and bulk modulus and derives updating equations for the elastic stresses.\n",
    "\n",
    "## Functional form\n",
    "\n",
    "$K = K_{ref}\\left(\\frac{p}{p_a}\\right)^m$  \n",
    "\n",
    "$G = G_{ref}\\left(\\frac{p}{p_a}\\right)^n$  \n",
    "\n",
    "where  \n",
    "\n",
    "$K_{ref}$ is the reference bulk modulus, and $G_{ref}$ is the reference shear modulus at $p = p_a$.  \n",
    "\n",
    "Note that  \n",
    "\n",
    "$K_{ref} = G_{ref}\\frac{1-2\\nu}{2\\left(1+\\nu\\right)}$  \n",
    "\n",
    "## Rate equations\n",
    "\n",
    "$\\dot{p} = K \\dot{\\epsilon_v^e}$  \n",
    "\n",
    "$\\dot{q} = 3G \\dot{\\epsilon_q^e}$  \n",
    "\n",
    "## Integrating rate equation for $\\epsilon_v^e$ term\n",
    "\n",
    "Note that\n",
    "\n",
    "$\\int_{t_i}^{t_{i+1}} \\frac{dp}{dt} dt = p_{i+1} - p_{i}$\n",
    "\n",
    "Rate equation after substituting expression for $K$\n",
    "\n",
    "$\\frac{\\dot{p}}{p^m} = \\frac{K_{ref}}{p_a^m}\\dot{\\epsilon_v^e}$  \n",
    "\n",
    "Integrate both sides  \n",
    "\n",
    "$\\int_{p_i}^{p_{i+1}} p^{-m} dp = \\int_{\\epsilon_v,i^e}^{\\epsilon_{v,i+1}^e} \\frac{K_{ref}}{p_a^m} d\\epsilon_v^e$ \n",
    "\n",
    "$\\frac{p_{i+1}^{1-m} - p_{i+1}^{1-m}}{1-m} = \\frac{K_{ref}}{p_a^m} d\\epsilon_v^e$  \n",
    "\n",
    "Solve for $p_{i+1}$  \n",
    "\n",
    "$p_{i+1} = exp\\left[\\frac{ln\\left[p_i^{1-m} + \\frac{ \\left(1-m\\right) K_{ref}}{p_a^m} d\\epsilon_v^e\\right]}{1-m}\\right]$  \n",
    "\n",
    "note: if m = 1  \n",
    "\n",
    "$p_{i+1} = p_i exp\\left(\\frac{K_{ref}}{p_a^m} d\\epsilon_v^e\\right)$  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3db5e47d",
   "metadata": {},
   "source": [
    "## Validation using explicit integration and many time steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "a3e452e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Explicit integration using 10000 time steps: 336.894964135377\n",
      "Exact solution: 336.9198188534109\n",
      "Forward Euler method using one time step: 221.46747898043233\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# validate rate equation for epsilon_v^e term using explicit integration\n",
    "Gref = 2.0*200**2\n",
    "nu = 0.3\n",
    "Kref = Gref*(1-2*nu)/(2*(1+nu))\n",
    "m = 1.0\n",
    "pa = 101.325\n",
    "\n",
    "N = 10000\n",
    "deps_ve = 0.01\n",
    "epsilon_ve = np.linspace(0,deps_ve,N)\n",
    "po = 100.0\n",
    "p_ex = po\n",
    "for i in range(1,N):\n",
    "    depsilon_ve = epsilon_ve[i] - epsilon_ve[i-1]\n",
    "    K = Kref*(p_ex/pa)**m\n",
    "    p_ex += K*depsilon_ve\n",
    "print('Explicit integration using ' + str(N) + ' time steps: ' + str(p_ex))\n",
    "\n",
    "if(m!=1):\n",
    "    pi1 = np.exp(np.log(po**(1-m) + (1-m)*Kref/pa**m*deps_ve)/(1-m))\n",
    "else:\n",
    "    pi1 = po*np.exp(Kref/pa*deps_ve)\n",
    "\n",
    "print('Exact solution: ' + str(pi1))\n",
    "\n",
    "print('Forward Euler method using one time step: ' + str(po + Kref*(po/pa)**m*deps_ve))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2838bb02",
   "metadata": {},
   "source": [
    "# Integrating rate equations for $\\epsilon_q^e$ term\n",
    "\n",
    "$\\frac{\\dot{q}}{\\dot{p}} = \\frac{3G\\dot{\\epsilon_q^e}}{K\\dot{\\epsilon_v^e}}$  \n",
    "\n",
    "$\\frac{G}{K} = \\frac{G_{ref}}{K_{ref}}\\left(\\frac{p}{p_a}\\right)^{n-m}$  \n",
    "\n",
    "$\\dot{q} = \\frac{3G_{ref}}{K_{ref}}\\frac{\\dot{\\epsilon_q^e}}{\\dot{\\epsilon_v^e}}\\dot{p}\\left(\\frac{p}{p_a}\\right)^{n-m}$  \n",
    "\n",
    "integrate both sides  \n",
    "\n",
    "$\\int_{q_i}^{q_{i+1}}dq = \\frac{3G_{ref}}{K_{ref}}\\frac{\\dot{\\epsilon_q^e}}{\\dot{\\epsilon_v^e}}\\int_{p_i}^{p_{i+1}}\\left(\\frac{p}{p_a}\\right)^{n-m}dp$  \n",
    "\n",
    "evaluate and solve for $q_{i+1}$  \n",
    "\n",
    "$q_{i+1} = q_i + \\frac{3G_{ref}}{K_{ref}}\\frac{\\dot{\\epsilon_q^e}}{\\dot{\\epsilon_v^e}}\\frac{p_a^{m-n}}{1+n-m}\\left(p_{i+1}^{1+n-m} - p_i^{1+n-m}\\right)$  \n",
    "\n",
    "note, if m = n, we have $\\nu$ = const, and the equation reduces to:  \n",
    "\n",
    "$q_{i+1} = q_i + \\frac{3G_{ref}}{K_{ref}}\\frac{\\dot{\\epsilon_q^e}}{\\dot{\\epsilon_v^e}}\\left(p_{i+1}-p_{i}\\right)$  \n",
    "\n",
    "furthermore, in the case that $\\dot{\\epsilon_v^e} = 0$, G is constant during shearing and we have:  \n",
    "\n",
    "$q_{i+1} = q_i + 3G_{ref}\\left(\\frac{p_i}{p_a}\\right)^n\\dot{\\epsilon_q^e}$  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d824fee5",
   "metadata": {},
   "source": [
    "## Validation using explicit integration and many time steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "62a7a6b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Explicit integration using 10000 time steps: p = 336.894964135377\n",
      "Exact solution: p = 336.9198188534109\n",
      "Explicit integration using 10000 time steps: q = 6559.895027281174\n",
      "Exact solution: q = 6560.2274055912985\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# validate rate equation for epsilon_v^e term using explicit integration\n",
    "Gref = 2.0*200**2\n",
    "nu = 0.3\n",
    "Kref = Gref*(1-2*nu)/(2*(1+nu))\n",
    "m = 1.0\n",
    "n = 0.5\n",
    "pa = 101.325\n",
    "\n",
    "N = 10000\n",
    "deps_ve = 0.01\n",
    "deps_qe = 0.02\n",
    "epsilon_ve = np.linspace(0,deps_ve,N)\n",
    "epsilon_qe = np.linspace(0,deps_qe,N)\n",
    "po = 100.0\n",
    "qo = 0.0\n",
    "p_ex = po\n",
    "q_ex = qo\n",
    "for i in range(1,N):\n",
    "    depsilon_ve = epsilon_ve[i] - epsilon_ve[i-1]\n",
    "    depsilon_qe = epsilon_qe[i] - epsilon_qe[i-1]\n",
    "    K = Kref*(p_ex/pa)**m\n",
    "    G = Gref*(p_ex/pa)**n\n",
    "    p_ex += K*depsilon_ve\n",
    "    q_ex += 3*G*depsilon_qe\n",
    "print('Explicit integration using ' + str(N) + ' time steps: p = ' + str(p_ex))\n",
    "\n",
    "if(m!=1):\n",
    "    pi1 = np.exp(np.log(po**(1-m) + (1-m)*Kref/pa**m*deps_ve)/(1-m))\n",
    "else:\n",
    "    pi1 = po*np.exp(Kref/pa*deps_ve)\n",
    "print('Exact solution: p = ' + str(pi1))\n",
    "\n",
    "if(m==n):\n",
    "    qi1 = qo + 3*Gref/Kref*deps_qe/deps_ve*(pi1-po)\n",
    "elif(deps_ve ==0):\n",
    "    qi1 = qo + 3*Gref*(po/pa)**n*deps_qe\n",
    "else:\n",
    "    qi1 = qo + 3*Gref/Kref*deps_qe/deps_ve*pa**(m-n)/(1+n-m)*(pi1**(1+n-m) - po**(1+n-m))\n",
    "print('Explicit integration using ' + str(N) + ' time steps: q = ' + str(q_ex))\n",
    "print('Exact solution: q = ' + str(qi1))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}