{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# 8.2 Regularized linear softening law for stable calculation in finite element code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sympy as sp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Let us assume that the fracture energy $G_\\mathrm{f}$ is known and  we want to improve the above model to deliver stable, mesh-independent results. Knowing that the energy dissipation is only happening in the one softening element, we can require that it alwas dissipates the same amount of energy. Recall that the fracture energy dissipated by a unit area of a stress free crack is evaluated as \n",
    "\\begin{align}\n",
    "G_\\mathrm{f} = \\int_0^{\\infty} f(w) \\, \\mathrm{d}w\n",
    "\\end{align}\n",
    "where $w$ represents the crack opening.\n",
    "In the studied case of linear softening with $f$ defined as\n",
    "\\begin{align}\n",
    "f(w) = f_\\mathrm{t}\\left( 1 - \\frac{1}{w_\\mathrm{f}} w \\right)\n",
    "\\end{align}\n",
    "and the corresponding fracture energy\n",
    "\\begin{align}\n",
    "G_\\mathrm{f} = \\int_0^{w_\\mathrm{f}} f(w) \\, \\mathrm{d}w = \\frac{1}{2} f_\\mathrm{t} w_\\mathrm{f}\n",
    "\\end{align}\n",
    "This kind of model is working correctly and reproduces exactly the amount of fracture energy needed to produce the unit area of the stress free crack. \n",
    "\n",
    "However, in a finite element model, the crack is not represented as a discrete line but is represented by a strain within the softening element of the size $L_s = L / N$. Thus, the softening is not related to crack opening displacement (COD) but to a strain in the softening zone using a softening function $\\phi(\\varepsilon_\\mathrm{s})$. This model was used in the example above and delivered mesh-dependent results with varying amount of fracture energy.\n",
    "\n",
    "To regularize the finite element model let us express the crack opening displacement as a product of the softening strain $\\varepsilon_\\mathrm{s}$ and the size of the softening zone $L_\\mathrm{s}$:\n",
    "\\begin{align}\n",
    "w = \\varepsilon_\\mathrm{s} L_\\mathrm{s}\n",
    "\\end{align}\n",
    "Now, the energy dissipated within the softening zone can be obtained as an integral over the history of the softening strain as\n",
    "\\begin{align}\n",
    "G_\\mathrm{f} = L_\\mathrm{s} \\int_0^{\\infty} \\phi(\\varepsilon_\\mathrm{s}) \\, \\mathrm{d}\\varepsilon_\\mathrm{s}\n",
    "\\end{align}\n",
    "Comming back to the model with linear softening, the integral of the strain based softening function is expressed as\n",
    "\\begin{align}\n",
    "\\int_0^{\\infty} \\phi(\\varepsilon_\\mathrm{s})  \\, \\mathrm{d} \\varepsilon_\\mathrm{s} = \\frac{1}{2} \\varepsilon_\\mathrm{f} f_\\mathrm{t}\n",
    "\\end{align}\n",
    "so that\n",
    "\\begin{align}\n",
    "G_\\mathrm{f} = \\frac{1}{2} L_\\mathrm{s} \\varepsilon_\\mathrm{f}\n",
    "f_\\mathrm{t}\n",
    "\\implies\n",
    "\\varepsilon_\\mathrm{f} = \\frac{2}{L_\\mathrm{s}} \\frac{G_\\mathrm{f}}{ f_\\mathrm{t}}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "L = 10.0\n",
    "E = 20000.0\n",
    "f_t = 2.4\n",
    "G_f = 0.0125 \n",
    "n_E_list = [5,10,1000]\n",
    "# run a loop over the different discretizations\n",
    "for n in n_E_list:  # n: number of element\n",
    "    L_s = L / n\n",
    "    eps_f = 2 * G_f / f_t / L_s \n",
    "    eps = np.array([0.0, f_t / E, eps_f / n])\n",
    "    sig = np.array([0.0, f_t, 0.0])\n",
    "    g_f = eps[-1] * sig[1] / 2\n",
    "    print('Is the energy constant?', g_f)\n",
    "    plt.plot(eps, sig, label='n=%i' % n)\n",
    "    plt.legend(loc=1)\n",
    "\n",
    "plt.xlabel('strain')\n",
    "plt.ylabel('stress')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Consider a tensile test with the dimensions of a cross section $100 \\times 100$ mm and length of $1000$ mm. The measured elongation of the beam during the test was as follows\n",
    "Assume the stiffness of $E = 28000$ MPa and strength $f_t = 3$ MPa. The length of the cohesive zone at which the crack localized was measured using the digital image correlation and set to $L_s = 0.05$ mm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Task"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Given the fracture energy, linear softening function, strength, E-modulus, bar length, manually calculate the force-displacement response of a tensile test. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": []
  }
 ],
 "metadata": {
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}