integrated X0304

0c425eab · Rostislav Chudoba · dd52efea · 0c425eab · 0c425eab · 0c425eab
Commit 0c425eab authored 1 year ago by Rostislav Chudoba
--- a/exercises/X0403-Plastic_model_by_GPT.ipynb
+++ b/exercises/X0403-Plastic_model_by_GPT.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a1ba93b5-c0c3-4ca1-9adf-fac77131a5d0",
+   "metadata": {},
+   "source": [
+    "# Can Chat GPT already generate a material model?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "54719971-159c-47ed-bff2-932261fee1be",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3c9d9841-af00-4218-a6a5-548d2a94213b",
+   "metadata": {},
+   "source": [
+    "Ask Chat GPT to generate an executable code with an example. How to formulate the question as precisely as possible?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9f39ac75-b548-426f-a5d7-d767f3249ded",
+   "metadata": {},
+   "source": [
+    "\"I want to implement a user material for non-linear finite element simulation using the predictor corrector concept. Many finite element tools allow for an externally implemented method to compute the stress for a prescribed strain increment and the set of state variables. Can you generate such a method in Python assuming a onedimensional stress-strain relation and ideal plasticity? Then, please provide a sample code with a prescribed sequence of strain increments that starts from zero and  monotonically increases and another one starting at zero and monotonically decreases. Plot the stress-strain history using matplotlib.\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cc2609b4-14c3-4e67-888f-678f4b535d56",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "94b97583-0942-4fef-8c02-e934ccb422e6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "def compute_stress(strain_increment, state_variables):\n",
+    "    # Material properties\n",
+    "    E = 100.0  # Young's modulus\n",
+    "    yield_stress = 10.0  # Yield stress\n",
+    "\n",
+    "    # Unpack state variables\n",
+    "    stress_prev, plastic_strain_prev = state_variables\n",
+    "\n",
+    "    # Compute trial stress\n",
+    "    stress_trial = stress_prev + E * strain_increment\n",
+    "\n",
+    "    # Check for plastic deformation\n",
+    "    if stress_trial > yield_stress:\n",
+    "        plastic_strain = plastic_strain_prev + (stress_trial - yield_stress) / E\n",
+    "        stress = yield_stress\n",
+    "    else:\n",
+    "        plastic_strain = plastic_strain_prev\n",
+    "        stress = stress_trial\n",
+    "\n",
+    "    return stress, plastic_strain\n",
+    "\n",
+    "# Prescribed sequence of strain increments\n",
+    "strain_increments = [0.1, 0.2, 0.15, 0.25]\n",
+    "\n",
+    "# Initial state variables\n",
+    "stress_initial = 0.0\n",
+    "plastic_strain_initial = 0.0\n",
+    "state_variables = (stress_initial, plastic_strain_initial)\n",
+    "\n",
+    "# Lists to store stress and strain values\n",
+    "stress_history = [stress_initial]\n",
+    "strain_history = [0.0]\n",
+    "\n",
+    "# Perform the simulation\n",
+    "for strain_increment in strain_increments:\n",
+    "    # Compute stress and update state variables\n",
+    "    stress, plastic_strain = compute_stress(strain_increment, state_variables)\n",
+    "    state_variables = (stress, plastic_strain)\n",
+    "\n",
+    "    # Update stress and strain history\n",
+    "    stress_history.append(stress)\n",
+    "    strain_history.append(strain_history[-1] + strain_increment)\n",
+    "\n",
+    "# Plot stress-strain history\n",
+    "plt.plot(strain_history, stress_history)\n",
+    "plt.xlabel('Strain')\n",
+    "plt.ylabel('Stress')\n",
+    "plt.title('Stress-Strain History')\n",
+    "plt.grid(True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "33260f76-9733-4efc-a6d0-76d00785cfe4",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "bmcs_env2",
+   "language": "python",
+   "name": "bmcs_env2"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
+%% Cell type:markdown id:a1ba93b5-c0c3-4ca1-9adf-fac77131a5d0 tags:
+
+# Can Chat GPT already generate a material model?
+
+%% Cell type:code id:54719971-159c-47ed-bff2-932261fee1be tags:
+
+``` python
+%matplotlib widget
+```
+
+%% Cell type:markdown id:3c9d9841-af00-4218-a6a5-548d2a94213b tags:
+
+Ask Chat GPT to generate an executable code with an example. How to formulate the question as precisely as possible?
+
+%% Cell type:markdown id:9f39ac75-b548-426f-a5d7-d767f3249ded tags:
+
+"I want to implement a user material for non-linear finite element simulation using the predictor corrector concept. Many finite element tools allow for an externally implemented method to compute the stress for a prescribed strain increment and the set of state variables. Can you generate such a method in Python assuming a onedimensional stress-strain relation and ideal plasticity? Then, please provide a sample code with a prescribed sequence of strain increments that starts from zero and  monotonically increases and another one starting at zero and monotonically decreases. Plot the stress-strain history using matplotlib."
+
+%% Cell type:code id:cc2609b4-14c3-4e67-888f-678f4b535d56 tags:
+
+``` python
+%matplotlib widget
+```
+
+%% Cell type:code id:94b97583-0942-4fef-8c02-e934ccb422e6 tags:
+
+``` python
+import matplotlib.pyplot as plt
+
+def compute_stress(strain_increment, state_variables):
+    # Material properties
+    E = 100.0  # Young's modulus
+    yield_stress = 10.0  # Yield stress
+
+    # Unpack state variables
+    stress_prev, plastic_strain_prev = state_variables
+
+    # Compute trial stress
+    stress_trial = stress_prev + E * strain_increment
+
+    # Check for plastic deformation
+    if stress_trial > yield_stress:
+        plastic_strain = plastic_strain_prev + (stress_trial - yield_stress) / E
+        stress = yield_stress
+    else:
+        plastic_strain = plastic_strain_prev
+        stress = stress_trial
+
+    return stress, plastic_strain
+
+# Prescribed sequence of strain increments
+strain_increments = [0.1, 0.2, 0.15, 0.25]
+
+# Initial state variables
+stress_initial = 0.0
+plastic_strain_initial = 0.0
+state_variables = (stress_initial, plastic_strain_initial)
+
+# Lists to store stress and strain values
+stress_history = [stress_initial]
+strain_history = [0.0]
+
+# Perform the simulation
+for strain_increment in strain_increments:
+    # Compute stress and update state variables
+    stress, plastic_strain = compute_stress(strain_increment, state_variables)
+    state_variables = (stress, plastic_strain)
+
+    # Update stress and strain history
+    stress_history.append(stress)
+    strain_history.append(strain_history[-1] + strain_increment)
+
+# Plot stress-strain history
+plt.plot(strain_history, stress_history)
+plt.xlabel('Strain')
+plt.ylabel('Stress')
+plt.title('Stress-Strain History')
+plt.grid(True)
+```
+
+%% Cell type:code id:33260f76-9733-4efc-a6d0-76d00785cfe4 tags:
+
+``` python
+```
--- a/exercises/X0501 - Constant bond-slip law using damage model.ipynb
+++ b/exercises/X0501 - Constant bond-slip law using damage model.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a1ba93b5-c0c3-4ca1-9adf-fac77131a5d0",
+   "metadata": {},
+   "source": [
+    "# What damage function would render a constant bond behavior?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "54719971-159c-47ed-bff2-932261fee1be",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7436d011-2494-4db6-bc46-bda5ea94d8f1",
+   "metadata": {},
+   "source": [
+    "Consider bond behavior"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "32eb7807-4b42-4f0d-af0b-e456c1ee8241",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pylab as plt\n",
+    "import numpy as np\n",
+    "import sympy as sp\n",
+    "sp.init_printing()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1af3886e-964c-467b-aa49-1a6f5c7d90ae",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tau_bar, s, E, omega = sp.symbols(r'\\bar{\\tau}, s, E, \\omega', nonnegative=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12ad00cd-d4a3-4b6c-b65b-99b71749b713",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tau_s = sp.Piecewise(\n",
+    "    (tau_bar, s>tau_bar / E),\n",
+    "    (E * s, True),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "523c7a0f-ac25-4678-b5ac-8c9d4ba7d0ba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "constitutive_law = sp.Eq(tau_s, (1 - omega) * E * s)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b1437b9f-7a12-48c2-823b-b8dea0a5e428",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "omega_s = sp.solve(constitutive_law, omega)[0]\n",
+    "omega_s"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e225eca6-3a92-497c-b60e-9cbc110a41e3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "get_omega_s = sp.lambdify((tau_bar, E, s), omega_s)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a5330769-b6a4-4b17-a581-7a0e25e76149",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "s_range = np.linspace(1e-10, 0.01, 100)\n",
+    "omega_s = get_omega_s(8, 28000, s_range)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3b6bc2dc-48d3-4258-b067-dffd07064cac",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget\n",
+    "plt.plot(s_range, omega_s)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "31dc8d77-ad7a-47a8-9e92-30addff2c582",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "def compute_stress_tangent(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic):\n",
+    "    stress = E * (strain - plastic_strain)\n",
+    "    tangent = E + H_kinematic + H_isotropic\n",
+    "    return stress, tangent\n",
+    "\n",
+    "def calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic):\n",
+    "    stress, tangent = compute_stress_tangent(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic)\n",
+    "    if stress > sigma_yield:\n",
+    "        plastic_strain += (stress - sigma_yield) / tangent\n",
+    "    return plastic_strain\n",
+    "\n",
+    "# Material parameters\n",
+    "E = 200.0  # Young's modulus\n",
+    "sigma_yield = 100.0  # Yield stress\n",
+    "H_kinematic = 10.0  # Kinematic hardening parameter\n",
+    "H_isotropic = 5.0  # Isotropic hardening parameter\n",
+    "\n",
+    "# Strain increments\n",
+    "n_steps = 10\n",
+    "strain_increments = np.linspace(0.0, 1.5, num=n_steps)\n",
+    "\n",
+    "# Initialize arrays to store stress and plastic strain history\n",
+    "stress_history = np.zeros_like(strain_increments)\n",
+    "plastic_strain = 0.0\n",
+    "\n",
+    "# Loop over strain increments\n",
+    "for i, strain_increment in enumerate(strain_increments):\n",
+    "    strain = strain_increment\n",
+    "    stress, _ = compute_stress_tangent(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic)\n",
+    "    stress_history[i] = stress\n",
+    "    plastic_strain = calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic)\n",
+    "\n",
+    "print(\"Strain increments:\", strain_increments)\n",
+    "print(\"Stress history:\", stress_history)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ecd9da8b-9dc1-4744-8ed9-7f186469a517",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Material parameters\n",
+    "E = 200.0  # Young's modulus\n",
+    "sigma_yield = 100.0  # Yield stress\n",
+    "H = 10.0  # Hardening parameter\n",
+    "\n",
+    "# Sample usage\n",
+    "strain = 0.6  # Assuming a strain increment of 0.1\n",
+    "plastic_strain = 0.0\n",
+    "\n",
+    "stress = calculate_stress(strain, plastic_strain, E, sigma_yield, H)\n",
+    "plastic_strain = calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H)\n",
+    "\n",
+    "print(f\"Strain: {strain}\")\n",
+    "print(f\"Stress: {stress}\")\n",
+    "print(f\"Plastic Strain: {plastic_strain}\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "439faf64-45ef-401e-95a2-2fccc440829b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget\n",
+    "eps_range = np.linspace(0, 1, 20)\n",
+    "for eps in eps_range:\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "828c3df2-aabe-4727-a5dd-583830ad4871",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "def compute_stress(strain, plastic_strain, E, sigma_yield, H):\n",
+    "    stress = E * (strain - plastic_strain)\n",
+    "    return stress\n",
+    "\n",
+    "def calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H):\n",
+    "    stress = compute_stress(strain, plastic_strain, E, sigma_yield, H)\n",
+    "    if stress > sigma_yield:\n",
+    "        plastic_strain += (stress - sigma_yield) / E\n",
+    "    return plastic_strain\n",
+    "\n",
+    "# Material parameters\n",
+    "E = 200.0  # Young's modulus\n",
+    "sigma_yield = 100.0  # Yield stress\n",
+    "H = 10.0  # Isotropic hardening modulus\n",
+    "\n",
+    "# Strain increments\n",
+    "n_steps = 10\n",
+    "strain_increments = np.linspace(0.0, 1.5, num=n_steps)\n",
+    "\n",
+    "# Initialize arrays to store stress and plastic strain history\n",
+    "stress_history = np.zeros_like(strain_increments)\n",
+    "plastic_strain = 0.0\n",
+    "\n",
+    "# Loop over strain increments\n",
+    "for i, strain_increment in enumerate(strain_increments):\n",
+    "    strain = strain_increment\n",
+    "    stress = compute_stress(strain, plastic_strain, E, sigma_yield, H)\n",
+    "    stress_history[i] = stress\n",
+    "    plastic_strain = calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H)\n",
+    "\n",
+    "print(\"Strain increments:\", strain_increments)\n",
+    "print(\"Stress history:\", stress_history)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6bbbbf43-2267-4d8f-85ed-967abd1fce63",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Sample usage\n",
+    "strain = np.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]])  # Assuming a strain increment of 0.1 in each component\n",
+    "plastic_strain = np.zeros((3, 3))\n",
+    "damage = 0.0\n",
+    "\n",
+    "material_model = MaterialModel(E=1.0, nu=0.3, sigma_yield=1.0, H_kinematic=0.2, H_isotropic=0.1, damage_coeff=0.1)\n",
+    "calculate_stress_tangent_damage(strain, plastic_strain, damage, material_model)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9c8cfbf2-99f6-452c-a0f9-45b888cc160c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plastic_strain = np.zeros((3, 3))\n",
+    "damage = 0.0\n",
+    "# Sample usage\n",
+    "strain = np.zeros((3,3), np.float_)\n",
+    "eps_range = np.linspace(0, 1.1, 20)\n",
+    "material_model = MaterialModel(E=1.0, nu=0.3, sigma_yield=1.0, H_kinematic=0., H_isotropic=0., damage_coeff=0.1)\n",
+    "sig_list = []\n",
+    "for eps in eps_range:\n",
+    "    strain[0, 0] = eps\n",
+    "    stress, tangent, plastic_strain, damage = material_model.compute_stress_tangent_and_damage_increment(strain, plastic_strain, damage)\n",
+    "    sig_list.append(stress[0,0])\n",
+    "sig_range = np.array(sig_list)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e4b368fa-d6b9-4f7d-b2fd-157300acff5e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plastic_strain"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "60b3d4f2-e1cf-4f7c-bdbc-e90e6880b22f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget\n",
+    "plt.plot(eps_range, sig_range)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9f39ac75-b548-426f-a5d7-d767f3249ded",
+   "metadata": {},
+   "source": [
+    "I want to implement a user material for non-linear finite element simulation using the predictor corrector concept. Many finite element tools allow for an externally implemented method to compute the stress for a prescribed strain increment and the set of state variables. Can you generate such a method in Python assuming a onedimensional stress-strain relation and ideal plasticity? Then, please provide a sample code with a prescribed sequence of strain increments that starts from zero and  monotonically increases and another one starting at zero and monotonically decreases. Plot the stress-strain history using matplotlib."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cc2609b4-14c3-4e67-888f-678f4b535d56",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib widget"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "94b97583-0942-4fef-8c02-e934ccb422e6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "def compute_stress(strain_increment, state_variables):\n",
+    "    # Material properties\n",
+    "    E = 100.0  # Young's modulus\n",
+    "    yield_stress = 10.0  # Yield stress\n",
+    "\n",
+    "    # Unpack state variables\n",
+    "    stress_prev, plastic_strain_prev = state_variables\n",
+    "\n",
+    "    # Compute trial stress\n",
+    "    stress_trial = stress_prev + E * strain_increment\n",
+    "\n",
+    "    # Check for plastic deformation\n",
+    "    if stress_trial > yield_stress:\n",
+    "        plastic_strain = plastic_strain_prev + (stress_trial - yield_stress) / E\n",
+    "        stress = yield_stress\n",
+    "    else:\n",
+    "        plastic_strain = plastic_strain_prev\n",
+    "        stress = stress_trial\n",
+    "\n",
+    "    return stress, plastic_strain\n",
+    "\n",
+    "# Prescribed sequence of strain increments\n",
+    "strain_increments = [0.1, 0.2, 0.15, 0.25, -0.1, -0.2]\n",
+    "\n",
+    "# Initial state variables\n",
+    "stress_initial = 0.0\n",
+    "plastic_strain_initial = 0.0\n",
+    "state_variables = (stress_initial, plastic_strain_initial)\n",
+    "\n",
+    "# Lists to store stress and strain values\n",
+    "stress_history = [stress_initial]\n",
+    "strain_history = [0.0]\n",
+    "\n",
+    "# Perform the simulation\n",
+    "for strain_increment in strain_increments:\n",
+    "    # Compute stress and update state variables\n",
+    "    stress, plastic_strain = compute_stress(strain_increment, state_variables)\n",
+    "    state_variables = (stress, plastic_strain)\n",
+    "\n",
+    "    # Update stress and strain history\n",
+    "    stress_history.append(stress)\n",
+    "    strain_history.append(strain_history[-1] + strain_increment)\n",
+    "\n",
+    "# Plot stress-strain history\n",
+    "plt.plot(strain_history, stress_history)\n",
+    "plt.xlabel('Strain')\n",
+    "plt.ylabel('Stress')\n",
+    "plt.title('Stress-Strain History')\n",
+    "plt.grid(True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cc4039db-aedb-47cb-98fd-faa8be6cd63f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "increments_rev"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "33260f76-9733-4efc-a6d0-76d00785cfe4",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "bmcs_env2",
+   "language": "python",
+   "name": "bmcs_env2"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
+%% Cell type:markdown id:a1ba93b5-c0c3-4ca1-9adf-fac77131a5d0 tags:
+
+# What damage function would render a constant bond behavior?
+
+%% Cell type:code id:54719971-159c-47ed-bff2-932261fee1be tags:
+
+``` python
+%matplotlib widget
+```
+
+%% Cell type:markdown id:7436d011-2494-4db6-bc46-bda5ea94d8f1 tags:
+
+Consider bond behavior
+
+%% Cell type:code id:32eb7807-4b42-4f0d-af0b-e456c1ee8241 tags:
+
+``` python
+import matplotlib.pylab as plt
+import numpy as np
+import sympy as sp
+sp.init_printing()
+```
+
+%% Cell type:code id:1af3886e-964c-467b-aa49-1a6f5c7d90ae tags:
+
+``` python
+tau_bar, s, E, omega = sp.symbols(r'\bar{\tau}, s, E, \omega', nonnegative=True)
+```
+
+%% Cell type:code id:12ad00cd-d4a3-4b6c-b65b-99b71749b713 tags:
+
+``` python
+tau_s = sp.Piecewise(
+    (tau_bar, s>tau_bar / E),
+    (E * s, True),
+)
+```
+
+%% Cell type:code id:523c7a0f-ac25-4678-b5ac-8c9d4ba7d0ba tags:
+
+``` python
+constitutive_law = sp.Eq(tau_s, (1 - omega) * E * s)
+```
+
+%% Cell type:code id:b1437b9f-7a12-48c2-823b-b8dea0a5e428 tags:
+
+``` python
+omega_s = sp.solve(constitutive_law, omega)[0]
+omega_s
+```
+
+%% Cell type:code id:e225eca6-3a92-497c-b60e-9cbc110a41e3 tags:
+
+``` python
+get_omega_s = sp.lambdify((tau_bar, E, s), omega_s)
+```
+
+%% Cell type:code id:a5330769-b6a4-4b17-a581-7a0e25e76149 tags:
+
+``` python
+s_range = np.linspace(1e-10, 0.01, 100)
+omega_s = get_omega_s(8, 28000, s_range)
+```
+
+%% Cell type:code id:3b6bc2dc-48d3-4258-b067-dffd07064cac tags:
+
+``` python
+%matplotlib widget
+plt.plot(s_range, omega_s)
+```
+
+%% Cell type:code id:31dc8d77-ad7a-47a8-9e92-30addff2c582 tags:
+
+``` python
+import numpy as np
+
+def compute_stress_tangent(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic):
+    stress = E * (strain - plastic_strain)
+    tangent = E + H_kinematic + H_isotropic
+    return stress, tangent
+
+def calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic):
+    stress, tangent = compute_stress_tangent(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic)
+    if stress > sigma_yield:
+        plastic_strain += (stress - sigma_yield) / tangent
+    return plastic_strain
+
+# Material parameters
+E = 200.0  # Young's modulus
+sigma_yield = 100.0  # Yield stress
+H_kinematic = 10.0  # Kinematic hardening parameter
+H_isotropic = 5.0  # Isotropic hardening parameter
+
+# Strain increments
+n_steps = 10
+strain_increments = np.linspace(0.0, 1.5, num=n_steps)
+
+# Initialize arrays to store stress and plastic strain history
+stress_history = np.zeros_like(strain_increments)
+plastic_strain = 0.0
+
+# Loop over strain increments
+for i, strain_increment in enumerate(strain_increments):
+    strain = strain_increment
+    stress, _ = compute_stress_tangent(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic)
+    stress_history[i] = stress
+    plastic_strain = calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H_kinematic, H_isotropic)
+
+print("Strain increments:", strain_increments)
+print("Stress history:", stress_history)
+```
+
+%% Cell type:code id:ecd9da8b-9dc1-4744-8ed9-7f186469a517 tags:
+
+``` python
+# Material parameters
+E = 200.0  # Young's modulus
+sigma_yield = 100.0  # Yield stress
+H = 10.0  # Hardening parameter
+
+# Sample usage
+strain = 0.6  # Assuming a strain increment of 0.1
+plastic_strain = 0.0
+
+stress = calculate_stress(strain, plastic_strain, E, sigma_yield, H)
+plastic_strain = calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H)
+
+print(f"Strain: {strain}")
+print(f"Stress: {stress}")
+print(f"Plastic Strain: {plastic_strain}")
+```
+
+%% Cell type:code id:439faf64-45ef-401e-95a2-2fccc440829b tags:
+
+``` python
+%matplotlib widget
+eps_range = np.linspace(0, 1, 20)
+for eps in eps_range:
+
+```
+
+%% Cell type:code id:828c3df2-aabe-4727-a5dd-583830ad4871 tags:
+
+``` python
+import numpy as np
+
+def compute_stress(strain, plastic_strain, E, sigma_yield, H):
+    stress = E * (strain - plastic_strain)
+    return stress
+
+def calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H):
+    stress = compute_stress(strain, plastic_strain, E, sigma_yield, H)
+    if stress > sigma_yield:
+        plastic_strain += (stress - sigma_yield) / E
+    return plastic_strain
+
+# Material parameters
+E = 200.0  # Young's modulus
+sigma_yield = 100.0  # Yield stress
+H = 10.0  # Isotropic hardening modulus
+
+# Strain increments
+n_steps = 10
+strain_increments = np.linspace(0.0, 1.5, num=n_steps)
+
+# Initialize arrays to store stress and plastic strain history
+stress_history = np.zeros_like(strain_increments)
+plastic_strain = 0.0
+
+# Loop over strain increments
+for i, strain_increment in enumerate(strain_increments):
+    strain = strain_increment
+    stress = compute_stress(strain, plastic_strain, E, sigma_yield, H)
+    stress_history[i] = stress
+    plastic_strain = calculate_plastic_strain(strain, plastic_strain, E, sigma_yield, H)
+
+print("Strain increments:", strain_increments)
+print("Stress history:", stress_history)
+```
+
+%% Cell type:code id:6bbbbf43-2267-4d8f-85ed-967abd1fce63 tags:
+
+``` python
+# Sample usage
+strain = np.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]])  # Assuming a strain increment of 0.1 in each component
+plastic_strain = np.zeros((3, 3))
+damage = 0.0
+
+material_model = MaterialModel(E=1.0, nu=0.3, sigma_yield=1.0, H_kinematic=0.2, H_isotropic=0.1, damage_coeff=0.1)
+calculate_stress_tangent_damage(strain, plastic_strain, damage, material_model)
+```
+
+%% Cell type:code id:9c8cfbf2-99f6-452c-a0f9-45b888cc160c tags:
+
+``` python
+plastic_strain = np.zeros((3, 3))
+damage = 0.0
+# Sample usage
+strain = np.zeros((3,3), np.float_)
+eps_range = np.linspace(0, 1.1, 20)
+material_model = MaterialModel(E=1.0, nu=0.3, sigma_yield=1.0, H_kinematic=0., H_isotropic=0., damage_coeff=0.1)
+sig_list = []
+for eps in eps_range:
+    strain[0, 0] = eps
+    stress, tangent, plastic_strain, damage = material_model.compute_stress_tangent_and_damage_increment(strain, plastic_strain, damage)
+    sig_list.append(stress[0,0])
+sig_range = np.array(sig_list)
+```
+
+%% Cell type:code id:e4b368fa-d6b9-4f7d-b2fd-157300acff5e tags:
+
+``` python
+plastic_strain
+```
+
+%% Cell type:code id:60b3d4f2-e1cf-4f7c-bdbc-e90e6880b22f tags:
+
+``` python
+%matplotlib widget
+plt.plot(eps_range, sig_range)
+```
+
+%% Cell type:markdown id:9f39ac75-b548-426f-a5d7-d767f3249ded tags:
+
+I want to implement a user material for non-linear finite element simulation using the predictor corrector concept. Many finite element tools allow for an externally implemented method to compute the stress for a prescribed strain increment and the set of state variables. Can you generate such a method in Python assuming a onedimensional stress-strain relation and ideal plasticity? Then, please provide a sample code with a prescribed sequence of strain increments that starts from zero and  monotonically increases and another one starting at zero and monotonically decreases. Plot the stress-strain history using matplotlib.
+
+%% Cell type:code id:cc2609b4-14c3-4e67-888f-678f4b535d56 tags:
+
+``` python
+%matplotlib widget
+```
+
+%% Cell type:code id:94b97583-0942-4fef-8c02-e934ccb422e6 tags:
+
+``` python
+import matplotlib.pyplot as plt
+
+def compute_stress(strain_increment, state_variables):
+    # Material properties
+    E = 100.0  # Young's modulus
+    yield_stress = 10.0  # Yield stress
+
+    # Unpack state variables
+    stress_prev, plastic_strain_prev = state_variables
+
+    # Compute trial stress
+    stress_trial = stress_prev + E * strain_increment
+
+    # Check for plastic deformation
+    if stress_trial > yield_stress:
+        plastic_strain = plastic_strain_prev + (stress_trial - yield_stress) / E
+        stress = yield_stress
+    else:
+        plastic_strain = plastic_strain_prev
+        stress = stress_trial
+
+    return stress, plastic_strain
+
+# Prescribed sequence of strain increments
+strain_increments = [0.1, 0.2, 0.15, 0.25, -0.1, -0.2]
+
+# Initial state variables
+stress_initial = 0.0
+plastic_strain_initial = 0.0
+state_variables = (stress_initial, plastic_strain_initial)
+
+# Lists to store stress and strain values
+stress_history = [stress_initial]
+strain_history = [0.0]
+
+# Perform the simulation
+for strain_increment in strain_increments:
+    # Compute stress and update state variables
+    stress, plastic_strain = compute_stress(strain_increment, state_variables)
+    state_variables = (stress, plastic_strain)
+
+    # Update stress and strain history
+    stress_history.append(stress)
+    strain_history.append(strain_history[-1] + strain_increment)
+
+# Plot stress-strain history
+plt.plot(strain_history, stress_history)
+plt.xlabel('Strain')
+plt.ylabel('Stress')
+plt.title('Stress-Strain History')
+plt.grid(True)
+```
+
+%% Cell type:code id:cc4039db-aedb-47cb-98fd-faa8be6cd63f tags:
+
+``` python
+increments_rev
+```
+
+%% Cell type:code id:33260f76-9733-4efc-a6d0-76d00785cfe4 tags:
+
+``` python
+```
--- a/tour4_plastic_bond/4_1_PO_multilinear_unloading.ipynb
+++ b/tour4_plastic_bond/4_1_PO_multilinear_unloading.ipynb
@@ -133,8 +133,8 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "po_trc.loading_scenario.profile = 'cyclic-nonsym-incr'\n",
-    "po_trc.loading_scenario.interact()"
+    "po_trc.loading_scenario.profile = 'cyclic-nonsym-const'\n",
+    "po_trc.loading_scenario.profile_.interact()"
   ]
  },
  {

 %% Cell type:markdown id: tags:

 <a id="top"></a>
 # **4.1: Loading, unloading and reloading**

 %% Cell type:markdown id: tags:

 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=551829) part 1 -->

 %% Cell type:code id: tags:

 ``` python
 from IPython.display import YouTubeVideo
 YouTubeVideo('SCuoGg7_8GE')
 ```

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left"> <img src="../icons/start_flag.png" alt="Previous trip" width="40" height="40">
    &nbsp; &nbsp; <b>Starting point</b> </div>

 %% Cell type:markdown id: tags:

 Once we have seen that there is a relation between the shape of the constitutive law and the  stress redistribution process within a bond zone which has an immediate consequence on the shape of the pullout curve, we extend our horizon to the case of **non-monotonic loading**.

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left"> <img src="../icons/destination.png" alt="Previous trip" width="40" height="40">
    &nbsp; &nbsp; <b>Where are we heading</b> </div>

 %% Cell type:markdown id: tags:

 To motivate a more general development of the model for the bond-slip behavior let us consider the case of non-monotonic loading. What happens in the material structure of the bond if the load is reduced and then it grows again?

 %% Cell type:markdown id: tags:

 # **Motivating examples**

 %% Cell type:markdown id: tags:

 Consider again the test double-sided pullout test introduced in [notebook 3.2](../tour3_nonlinear_bond/3_2_anchorage_length.ipynb#trc_pullout_study). Within this test series studying the nonlinear bond behavior of carbon fabrics, a loading scenario with introducing several **unloading and reloading steps** has been included for the specimen with the total length of 300 and 400 mm. The obtained measured response looked as follows
 ![image](../fig/test_unloading.png)

 %% Cell type:markdown id: tags:

 **Question:** Can the  pullout models from tours 2 and 3 be used to describe such kind of behavior? What constitutive law  can reproduce this kind of behavior?

 **Answer:** The bond-slip laws presented so far only considered monotonically increasing load and cannot reproduce the real material behavior upon unloading.
 The constant-bond slip model in Tour 2 and the multilinear bond-slip model exemplified in Tour 3 did not consider any change of behavior upon unloading.
 To document this statement and to show how to introduce unloading and reloading into the interactive numerical models  that we used in Tour 3, an example has been prepared showing how to introduce a more complex type of loading into the model.

 %% Cell type:markdown id: tags:

 <a id="trc_study_monotonic"></a>
 ### **Example 1:** TRC specimen with unloading

 Let us reuse the geometrical, material and algorithmic parameters already specified in the
 [case study on carbon fabric bond](../tour3_nonlinear_bond/3_2_anchorage_length.ipynb#case_study_1)

 %% Cell type:code id: tags:

 ``` python
 %matplotlib widget
 import matplotlib.pyplot as plt
 import numpy as np
 from bmcs_cross_section.pullout import PullOutModel1D
 po_trc = PullOutModel1D(n_e_x=100, w_max=3) # mm
 po_trc.geometry.L_x=150 # [mm]
 po_trc.time_line.step = 0.02
 po_trc.cross_section.trait_set(A_m=1543, A_f=16.7, P_b=10)
 po_trc.material_model='multilinear'
 po_trc.material_model_.trait_set(
    E_m=28000, E_f=170000,
    s_data = '0, 0.1, 0.5, 1, 2, 3,   4.5, 6',
    tau_data = '0, 5.4, 4.9, 5, 6, 7, 8.3, 9'
 );
 ```

 %% Cell type:markdown id: tags:

 In addition, let us now change the `loading_scenario` to `cyclic` load.
 Note that this polymorphic attribute changes the 'shadowed' object
 which generates the time function for the scaling of the load. The
 shadowed object representing the cyclic loading scenario can then
 be accessed as an attribute `loading_scenario_`. It can be rendered
 using the `interact` method, like any other model component in the
 model tree

 %% Cell type:code id: tags:

 ``` python
-po_trc.loading_scenario.profile = 'cyclic-nonsym-incr'
-po_trc.loading_scenario.interact()
+po_trc.loading_scenario.profile = 'cyclic-nonsym-const'
+po_trc.loading_scenario.profile_.interact()
 ```

 %% Cell type:markdown id: tags:

 The parameters of the loading scenario are now visible in the rendered window. They can be
 assigned through the `trait_set` method.

 %% Cell type:code id: tags:

 ``` python
 po_trc.loading_scenario.profile_.trait_set(number_of_cycles=1,
                              unloading_ratio=0.0,
                               numbe_of_increments=200);
 ```

 %% Cell type:markdown id: tags:

 The model can now be executed and rendered.

 %% Cell type:code id: tags:

 ``` python
 po_trc.reset()
 po_trc.run()
 po_trc.interact()
 ```

 %% Cell type:markdown id: tags:

 Apparently, there is no significant difference directly visible. However, if you browse through the history, it will be obvious that unloading is running along the same path as loading, which contradicts to the behavior observed in experiments.

 %% Cell type:markdown id: tags:

 ### **Example 2:** CFRP with unloading

 Even more evident demonstration of the non-physicality of the applied model can be provided by reusing the calibrated CFRP model. Let us apply the same loading scenario again

 %% Cell type:code id: tags:

 ``` python
 A_f = 16.67 # [mm^2]
 A_m = 1540.0 # [mm^2]
 p_b = 100.0 #
 E_f = 170000 # [MPa]
 E_m = 28000 # [MPa]
 pm = PullOutModel1D()
 pm.material_model = 'multilinear'
 pm.material_model_.s_data = "0, 0.1, 0.4, 4"
 pm.material_model_.tau_data = "0, 8, 0, 0"
 pm.sim.tline.step = 0.01 # 100 time increments
 pm.cross_section.trait_set(A_f=A_f, P_b=p_b, A_m=A_m)
 pm.geometry.L_x = 300 # length of the specimen [mm]
 pm.w_max = 1.8 # maximum control displacement [mm]
 pm.n_e_x = 100 # number of finite elements
 pm.loading_scenario.profile = 'cyclic-nonsym-const'
 pm.loading_scenario.profile_.trait_set(number_of_cycles=2,
                              unloading_ratio=0.0,
                               numbe_of_increments=200)
 pm.run()
 pm.interact()
 ```

 %% Cell type:markdown id: tags:

 Browsing through the history will reveal that the debonded zone gets recovered upon unloading as the process zone travels back through the bond zone and the pullout tests gets completely healed once $w = 0$. This motivates the question: **How to introduce irreversible changes in the material structure** so that we reflect the physics behind the scenes in a correct way.

 %% Cell type:markdown id: tags:

 <a id="plasticity_and_damage"></a>
 # **Physical explanation of the bond behavior**

 %% Cell type:markdown id: tags:

 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=551829) part 2 -->

 %% Cell type:code id: tags:

 ``` python
 YouTubeVideo('btj1SgUjBeA')
 ```

 %% Cell type:markdown id: tags:

 Regarding a small segment of the bond interface with a nearly constant shear stress $\tau$ and constant
 slip $s$ let us try to describe the correspondence between the micro- and meso-scopic mechanisms that actually govern
 bond-slip relation $\tau(s)$. To classify the elementary mechanisms behind the observed bond behavior, let us can idealize
 the material structure of the bond using two types of bindings:

 - as a **series of elastic springs** that can break once they achieve their strength, and
 - as a **series of asperities** representing an unevenness or roughness of the surface area.
   Sliding of these surfaces is accompanied with a stress required to cross over the asperities. During this process the asperities can deform elastically or get abraded.

 With this simple image of the material interface we can try to associate the structure of the bond layer with the basic types of inelastic material behavior, namely, damage and plasticity. In most cases both types of material behavior occur. The question is, how to distinguish and quantify them.

 %% Cell type:markdown id: tags:

 **Damage:** Consider an infinitesimal material point of the interface exhibiting softening response upon monotonic loading. Upon unloading, this material point would unload to origin. This behavior can be reproduced by an idealization of the material point as a series of elastic springs. Assuming further a scatter of spring strength, the failure of one spring upon increasing slip displacement $s$ leads to the reduction of the bond stiffness and overall stress level $\tau$. Since all the remaining springs are elastic, they would all unload to the origin.

 %% Cell type:markdown id: tags:

 ![damage](../fig/damage.png)

 %% Cell type:markdown id: tags:

 **Plasticity:**
 On the other hand, if a bond stress is transferred by asperities, we can distinguish two phases of response. Before the stress necessary to skip over an asperity has been reached, the slip is zero. At the onset of sliding over the asperities, the slip can increase at constant level of stress. Once the stress level is reduced, the slip does not change. Or putting it differently, the stress needed to induce sliding over the surface remains constant. Once the sliding direction is changed, also the sign of the shear stress must change. The described type of behavior is referred to as ideally plastic. This kind of material response is shown as the curve (2)
 ![plasticity](../fig/plasticity.png)

 If asperities can deform before reaching the onset of sliding, an elastic deformation can be observed as indicated
 by the curve (3).

 %% Cell type:markdown id: tags:

 **Interactions:** Usually, both damage and plasticity encounter simultaneously in the material behavior. Both introduce energy dissipation. While damage is connected with the reduction of the material stiffness, plasticity leads to permanent deformation after unloading of the structure. To understand the effect of damage and plasticity at the level of a material point, we will construct simple models that can be used within the BMCS interactive sheets to understand the two types of material behavior in detail.

 %% Cell type:markdown id: tags:

 Corresponding to this classification, we will address the two categories of behavior in a more
 detailed way:

 1. starting with interactive studies for non-monotonic loading scenarios,
 2. providing a mathematical framework for the description of plasticity and damage exemplified for the bond-behavior at a material point level
 3. studying the effect of plasticity and damage at the structural level using pullout tests.

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left;width:45%;display:inline-table;"> <img src="../icons/previous.png" alt="Previous trip" width="50" height="50">
    &nbsp; <a href="../tour3_nonlinear_bond/3_2_anchorage_length.ipynb#top">3.2 Anchorage length</a>
 </div><div style="background-color:lightgray;text-align:center;width:10%;display:inline-table;"> <a href="#top"><img src="../icons/compass.png" alt="Compass" width="50" height="50"></a></div><div style="background-color:lightgray;text-align:right;width:45%;display:inline-table;">
    <a href="4_2_BS_EP_SH_I_A.ipynb#top">4.2 Basic concept of plasticity</a>&nbsp; <img src="../icons/next.png" alt="Previous trip" width="50" height="50"> </div>

 %% Cell type:code id: tags:

 ``` python
 ```

--- a/tour4_plastic_bond/4_2_BS_EP_SH_I_A.ipynb
+++ b/tour4_plastic_bond/4_2_BS_EP_SH_I_A.ipynb
@@ -199,6 +199,14 @@
    "The latter issue means, that we can directly quantify the plastic slip $s_\\mathrm{pl}$ and keep it as a memory of material response. This concept of state variable is essential for any inelastic material model."
   ]
  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## How to implement this? \n",
+    "Can Chat GPT do it? [Here is an example.](../exercises/X0403-Plastic_model_by_GPT.ipynb)\n"
+   ]
+  },
  {
   "cell_type": "markdown",
   "metadata": {},
 %% Cell type:markdown id: tags:
  
 <a id="top"></a>
 # **4.2 Basic concept of plasticity**
  
 %% Cell type:markdown id: tags:
  
 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=615712) part 1 -->
  
 %% Cell type:code id: tags:
  
 ``` python
 from IPython.display import YouTubeVideo
 YouTubeVideo('OLMa5bYffSA')
 ```
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/start_flag.png" alt="Previous trip" width="40" height="40">
    &nbsp; &nbsp; <b>Starting point</b> </div>
  
 %% Cell type:markdown id: tags:
  
 With the basic distinction of damage and plasticity we need a mathematical instrument that can correctly describe both of these cases. We start with the latter one.
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/destination.png" alt="Previous trip" width="40" height="40">
    &nbsp; &nbsp; <b>Where are we heading</b> </div>
  
 %% Cell type:markdown id: tags:
  
 In fact, we will revisit the constant-bond slip interface addressed using analytical models in **Tour 2**. In contrast to the analytical models formulated there, we will now include general loading scenarios. The development will allow us to include also hardening and softening in the framework of plasticity theory. The one-dimensional interface provides an ideal framework to explain the applied assumptions and the way how to transform them into mathematical expressions that can be used to simulate the material behavior.
  
 %% Cell type:markdown id: tags:
  
 # **Plastic interactive exercise**
  
 Before starting with a theoretical description, let us gain some intuition by observing an elementary example of plastic behavior at the level of a single material point. For this purpose, we import an explorer which contains a plastic bond-slip model with four parameters:
  
 | Parameter name | Symbol | Description |
 | - | - | - |
 | E | $E$ | material stiffness |
 | tau_bar | $\bar{\tau}$ | yield stress |
 | K | $K$ | isotropic hardening modulus |
 | gamma | $\gamma$ | kinematic hardening modulus |
  
 %% Cell type:code id: tags:
  
 ``` python
 %matplotlib widget
 from plastic_app.bs_model_explorer import BSModelExplorer
 bs = BSModelExplorer(name='plasticity explorer')
 bs.bs_model.trait_set(E=100, K=0, gamma=0, tau_bar=20)
 bs.n_steps=200
 bs.interact()
 ```
  
 %% Cell type:markdown id: tags:
  
 Let us now render the model `bs`. The load is imposed by setting the parameter `s_1` and running the calculation without resetting the control bar. To see directly the effect, we impose a loading and unloading step before the model is rendered.
 <a id="plastic_explorer"></a>
  
 %% Cell type:code id: tags:
  
 ``` python
 bs.s_1 = 3 # load to 3 mm
 bs.run()
 bs.s_1 = 0 # unload to 0
 bs.run()
 bs.s_1 = 3 # load to 3 mm
 bs.run()
 bs.s_1 = 0 # unload to 0
 bs.run()
 bs.s_1 = 3 # load to 3 mm
 bs.run()
 bs.s_1 = 0 # unload to 0
 bs.run()
 bs.interact() # render
 ```
  
 %% Cell type:markdown id: tags:
  
 **Explore!** Find out the meaning of the parameters $K$ and $\gamma$ by investigating the model.
  
 %% Cell type:markdown id: tags:
  
 ## **Observations**
 * By setting the parameters $K$ and $\gamma$ equal to zero, yielding is horizontal, and recovers the case of **ideal plasticity**, or the already known stick-slip bond behavior discussed in **Tour 2**
  
 * The hardening moduli $K$ and $\gamma$  both **change the size of the elastic domain during the yielding**.
  
 * While **isotropic hardening** controlled by $K$ enlarges the whole elastic domain upon yielding, the **kinematic hardening** $\gamma$ keeps the size of the elastic domain constant. Instead, it induces the shift of the whole elastic domain in the direction of yielding.
  
 %% Cell type:markdown id: tags:
  
 # **Let us take a pencil and describe it**
  
 %% Cell type:markdown id: tags:
  
 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=615712) part 2 -->
  
 %% Cell type:code id: tags:
  
 ``` python
 YouTubeVideo('BlsIRMT7SEo')
 ```
  
 %% Cell type:markdown id: tags:
  
 ![image.png](attachment:8d669bba-4aed-447d-9dd3-e793dafe1313.png)
  
 %% Cell type:markdown id: tags:
  
 ![image.png](attachment:0dbede00-6617-4020-b5f1-a27f01ad82b5.png)
  
 %% Cell type:markdown id: tags:
  
 **Summary:** The above derivation concludes that upon yielding, the rate of stress must be zero. Even though simple, the describe set of equations allow us to decide
 - whether or not we are in a reversible state, and
 - if in irreversible regime, directly quantify the amount of changes in the material state.
  
 The latter issue means, that we can directly quantify the plastic slip $s_\mathrm{pl}$ and keep it as a memory of material response. This concept of state variable is essential for any inelastic material model.
  
 %% Cell type:markdown id: tags:
  
+## How to implement this?
+Can Chat GPT do it? [Here is an example.](../exercises/X0403-Plastic_model_by_GPT.ipynb)
+
+%% Cell type:markdown id: tags:
+
 # **Plastic model with isotropic and kinematic hardening**
  
 %% Cell type:markdown id: tags:
  
 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=615712) part 3 -->
  
 %% Cell type:code id: tags:
  
 ``` python
 YouTubeVideo('QOALP75E4uU')
 ```
  
 %% Cell type:markdown id: tags:
  
 Let us now express the manual derivation of the ideally plastic behavior using the algebra system. This step will help us to apply the same logic to more general assumptions including changes of elastic domain during the yielding process. In particular, we want to use the same skeleton of derivation to include the effects of hardening, both isotropic and kinematic.
  
 %% Cell type:markdown id: tags:
  
 The required packages include again plotting `matplotlib`, symbolic algebra `sympy`, and numerical manipulation of arrays `numpy`.
  
 %% Cell type:code id: tags:
  
 ``` python
 import matplotlib.pyplot as plt
 import sympy as sp
 import numpy as np
 sp.init_printing()
 ```
  
 %% Cell type:markdown id: tags:
  
 ## **Elastic, ideally plastic behavior**
 Assuming a constant yielding stress $\tau_\mathrm{Y}$ and no hardening the possible paths along which the stress strain states can develop are depicted in Figure 1.
 ![image.png](attachment:image.png)
  
 %% Cell type:markdown id: tags:
  
 &nbsp;<font color='blue'>
 **Naming conventions:**
 - Variables with trailing underscore (e.g. `f_tau_`) denote `sympy` expressions.
 - Variables denoting `sympy` symbols (e.g. `tau`) have no underscore at the end and have a name which is close to the mathematical symbol (e.g. $\tau$)
 - Mathematical symbols defined as string in `sp.symbols(r'\tau')` use `latex` syntax to introduce Greek symbols, super and subindexes. This makes the pretty printing of expressions possible.
 - In an implemented algorithm at the end of the notebook, the Python variables containing the numerical values of the material parameters $E_b$, $\tau_\mathrm{Y}$, etc. are denoted with a leading underscore `_E_b` and `_tau_Y` to avoid name collisions within the notebook
 </font>
  
 %% Cell type:markdown id: tags:
  
 ### Yield condition
 <a id="yield_function"></a>
 Elastic domain is defined by the inequality equation with the stress variable $\tau$ and yielding stress $\tau_Y$
  
 ![image.png](attachment:image.png)
 \begin{align}
 f := | \tau| - \tau_Y \le 0
 \label{eq:f_perfect_plasticity}
 \end{align}
 which we will denote `f_tau_`
  
 %% Cell type:code id: tags:
  
 ``` python
 tau = sp.symbols(r'\tau')
 tau_Y = sp.symbols(r'tau_Y', real=True, nonnegative=True)
 f_tau_ = sp.sqrt( tau*tau ) - tau_Y
 f_tau_
 ```
  
 %% Cell type:markdown id: tags:
  
 ### Elastic behavior
 <a id="elastic_law"></a>
 Expression evaluating the stress from the elastic slip is denoted as `tau_`
  
 ![image.png](attachment:image.png)
 \begin{align}
 \tau = E_\mathrm{b} (s - s_\mathrm{pl})
 \label{eq:elastic_behavior}
 \end{align}
  
 %% Cell type:code id: tags:
  
 ``` python
 s, s_pl = sp.symbols('s, s_pl') # total slip and plastic slip
 E_b = sp.symbols('E_b', positive=True) # E modulus
 tau_ =  E_b * (s - s_pl) # stress depends on the elastic part of slip
 tau_
 ```
  
 %% Cell type:markdown id: tags:
  
 As a next step, to describe the state with a memory, generally reflecting the relation of the current material state and its yield condition we need to define the evolution equation and satisfy the consistency condition.
  
 %% Cell type:markdown id: tags:
  
 ### Evolution equation
 Let us be more precise in defining the goal of the derivation: Inelastic behavior is path dependent. That means: <font color="brown"> **during yielding, the value of stress does not depend only on the value of slip but also on the history of loading in a material point.**</font> Thus, we cannot get an explicit relation between the stress and slip, we only can obtain a relation between their trends, or more precisely, rates. From now on, we do not consider the values of the variables but the rate of their change.
  
 This will allow us to formulate the needed criteria
 * <font color="green">to find out **in which direction does the yielding process proceed** from a current state, and</font>
 * <font color="green">to account for **the loading history that a material point** experienced in the past.</font>
  
 Our goal is to establish the relation between stress and strain (or shear and slip) in a rate form,<br> i.e.
 <font color="red"> $\dot{\tau}$ and $\dot{s}$</font>, which is a short form of
 <font color="red"> $\displaystyle{\frac{\partial \tau}{\partial t}}$ and $\displaystyle{\frac{\partial s}{\partial t}}$, </font> respectively.
  
 %% Cell type:markdown id: tags:
  
 Regarding Figure 2, we can directly answer the question, that the sign of the plastic slip rate $\dot{s}_\mathrm{pl}$ upon  positive or negative yielding, i.e. when $\tau = \tau_\mathrm{Y}$ or $\tau = -\tau_\mathrm{Y}$, respectively, will be linked with the sign of the current stress $\mathrm{sign}({\tau})$
  
 ![image.png](attachment:image.png)
  
 %% Cell type:markdown id: tags:
  
 **Compass for the inelastic domain:** Let us postulate, that the amount of yielding can be mathematically controlled by a nonnegative variable $\lambda$. Yielding proceeds in the direction normal to the yield surface in stress space, i.e. $\partial f/ \partial \tau$
  
 ![image.png](attachment:image.png)
 \begin{align}
 \dot{s}_\mathrm{pl}
 \label{eq:s_p_rage} = \lambda \frac{\partial f}{\partial \tau}
 \end{align}
  
 %% Cell type:code id: tags:
  
 ``` python
 lambda_ = sp.symbols(r'\lambda', nonnegative=True)
 dot_s_pl_ = lambda_ * f_tau_.diff(tau)
 dot_s_pl_
 ```
  
 %% Cell type:markdown id: tags:
  
 But how far do we go along the direction given by the gradient $\partial f / \partial \tau$? The amount of yielding
 is now controlled by a new unknown variable $\lambda$?
  
 <a id="kuhn_tucker"></a>
 The idea that helps to sort this out and to mathematically distinguish the loading steps into elastic and inelastic ones is provided by the **Kuhn-Tucker condition**
  
 ![image.png](attachment:image.png)
  
 \begin{align}
  \lambda f = 0, \; \lambda > 0,\; f \le 0
  \label{eq:kuhn_tucker}
 \end{align}
  
 which acts as a switch between either growing $\lambda > 0$ (yielding) or negative $f < 0$ (elastic loading, or unloading) or both $\lambda = 0$ and $f=0$ (neutral loading). Still, it does not provide the additional condition to get $\lambda$ resolved in case of yielding. The last needed equation can be obtained by realizing that even the yielding condition $f$ is a subject of change during yielding so that we need to constrain its evolution $\dot{f}$.
  
 %% Cell type:markdown id: tags:
  
 ### Consistency condition
 If the material is yielding, then $f = 0$ and $\lambda \le 0$. Moreover, <font color="green">the value of $f$ must **remain zero** during the whole yielding process</font>. This means that the rate of the yield function, i.e.
 \begin{align}
 \frac{ \mathrm{d} f}{\mathrm{d} t} = \dot{f} = 0
  \label{eq:consistency}
 \end{align}
 must remain zero as well. This is the sought condition that can be used to resolve for $\lambda$ and, thus, the evolution of the state variable $s_\mathrm{pl}$. **That's all** regarding the criteria and system becomes solvable. What follows is application of chain derivatives and substitutions.
  
 To construct the rate of the yield condition let us recall that it depends on stress, which in turn depends on the control slip and on the plastic slip
 \begin{align}
 f:= f( \tau (s, s^\mathrm{pl}) )
 \nonumber
 \end{align}
 Thus, to obtain $\dot{f}$ we  can apply the chain rule
 \begin{align}
 \dot{f} = \frac{\mathrm{d} f}{\mathrm{d} t} &=
 \frac{\partial f}{\partial \tau} \frac{\mathrm{d} \tau}{\mathrm{d} t}
 = \frac{\partial f}{\partial \tau}
 \left(
 \frac{\partial \tau}{\partial s} \frac{\mathrm{d} s}{\mathrm{d} t} +
 \frac{\partial \tau}{\partial s^\mathrm{pl}} \frac{\mathrm{d} s^\mathrm{pl}}{\mathrm{d} t}
 \right)
 = \frac{\partial f}{\partial \tau} E_\mathrm{b}\left( \dot{s} - \dot{s}^\mathrm{pl} \right)
  \label{eq:f_chain_derivatives} \\
 &= \mathrm{sign}(\tau) E_\mathrm{b} \dot{s} - E_\mathrm{b} \lambda = 0
 \end{align}
  
 %% Cell type:markdown id: tags:
  
 ![image.png](attachment:image.png)
  
 %% Cell type:markdown id: tags:
  
 To realize this derivation in `sympy` let us transform the elastic relation from absolute values into the rate form
 <a id="elastic_rate"></a>
 \begin{align}
 \dot{\tau} = E_\mathrm{b}(\dot{s} - \dot{s}_\mathrm{pl})
 \label{eq:elastic_rate}
 \end{align}
  
 %% Cell type:code id: tags:
  
 ``` python
 dot_tau = sp.Symbol(r'\dot{\tau}') # will be used later just to display the derived expressions next to this symbol
 dot_s, dot_s_pl = sp.symbols(r'\dot{s}, \dot{s}^\mathrm{pl}') # control slip and plastic slip rates
 dot_tau_ = E_b * (dot_s - dot_s_pl)
 dot_tau_
 ```
  
 %% Cell type:markdown id: tags:
  
 To obtain the rate of the yield function
 \begin{align}
 \dfrac{\partial f}{\partial \tau} \dot{\tau}
 \end{align}
 substitute differentiate the expression `f_tau_` with respect to `tau` using the `.diff` method and multiply by
 the expression `dot_tau_` with `dot_s_pl` substituted by the expression `dot_s_pl_` to obtain the result derived above manually
  
 %% Cell type:code id: tags:
  
 ``` python
 dot_f_ = f_tau_.diff(tau) * dot_tau_.subs(dot_s_pl, dot_s_pl_)
 sp.simplify(dot_f_)
 ```
  
 %% Cell type:markdown id: tags:
  
 Now, by setting this expression to zero and resolving for $\lambda$ we again obtain
  
 %% Cell type:code id: tags:
  
 ``` python
 lambda_solved = sp.solve( dot_f_, lambda_)[0]
 lambda_solved
 ```
  
 %% Cell type:markdown id: tags:
  
 which is equivalent to
  
 %% Cell type:markdown id: tags:
  
 ![image.png](attachment:image.png)
  
 %% Cell type:markdown id: tags:
  
 ### Solution
  
 **As already said**, the conclusion from the shown derivation is rather trivial - for perfectly plastic model without hardening or softening, the plastic slip grows along with the control slip with the same rate. Thus, the resolved evolution equation for plastic slip reduces to the equality
  
 %% Cell type:markdown id: tags:
  
 ![image.png](attachment:image.png)
  
 %% Cell type:code id: tags:
  
 ``` python
 {dot_s_pl : dot_s_pl_.subs(lambda_, lambda_solved) }
 ```
  
 %% Cell type:markdown id: tags:
  
 **Resolved evolution equation:** As we will show in the next extension of the plastic model accounting for isotropic hardening, the derivation of a evolution equation is the actual goal in the formulation of an inelastic model. The numerical simulation of the response is done by integrating the evolution equations along the time axis.
  
 %% Cell type:markdown id: tags:
  
 **However**, in this continuous representation it is possible to relate the rate of stress and strain (slip) explicitly rendering the result that, upon yielding, the stress rate is zero, i.e. $\tau = \tau_\mathrm{Y}$.
  
 %% Cell type:markdown id: tags:
  
 ![image.png](attachment:image.png)
  
 %% Cell type:code id: tags:
  
 ``` python
 {dot_tau: sp.simplify(dot_tau_.subs(dot_s_pl, dot_s_pl_).subs(lambda_, lambda_solved))}
 ```
  
 %% Cell type:markdown id: tags:
  
 ## **Hardening - elastic domain changes during yielding**
  
 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=615712) part 4 -->
  
 %% Cell type:code id: tags:
  
 ``` python
 YouTubeVideo('Nu86DjtFKQk')
 ```
  
 %% Cell type:markdown id: tags:
  
 The derivation scheme presented above can be used for more general assumptions describing the inelastic range of the material behavior. In this section, we will include the next dissipative mechanism, namely isotropic hardening. In simple terms we now allow the growth of the elastic domain with a variable $Z$.
  
 ![image.png](attachment:image.png)
 ### Yield condition
 Correspondingly, the yield condition takes the form
 \begin{align}
 f = | \tau | - (\tau_Y + Z)
 \label{eq:f_iso_hardening}
 \end{align}
  
 %% Cell type:markdown id: tags:
  
 Using `sympy`, we introduce this variable and define the yield function as `f_tau_Z_`.
  
 **Side remark:** <font color="gray"> To document that the derivation is systematic and can incorporate further assumption, we prepare also the variable $X$ for inclusion of kinematic hardening. In the accompanying video, we will show the procedure of extending the model with the kinematic hardening behavior. The model can be extended by uncommenting the lines referring to kinematic hardening, i.e.
 \begin{align}
 [f = | \tau - X | - (\tau_Y + Z)]
 \label{eq:f_iso_hardening}
 \end{align}</font>
  
 %% Cell type:code id: tags:
  
 ``` python
 Z = sp.symbols('Z')
 # X = sp.symbols('X') # prepared for kinematic hardening
 f_tau_Z_ = sp.sqrt( (tau)**2 ) - (tau_Y + Z) # extend the definition of elastic domain with to kinematic hardening
 f_tau_Z_
 ```
  
 %% Cell type:markdown id: tags:
  
 ### Hardening behavior
 The hardening force $Z$ is an additional stress measure. It must be accompanied with a corresponding kinematic variable, a hardening slip and a material parameter - hardening modulus $K$. Let us assume linear hardening relation
 <a id="hardening_law"></a>
 \begin{align}
 Z = K z \\
 \color{gray}{ [X = \gamma \alpha] }
 \label{eq:isotropic_hardening}
 \end{align}
 This relation will be needed to solve for the consistency condition in the rate form. Thus, in `sympy` we shall introduce directly $\dot{Z} = K \dot{z}$
  
 %% Cell type:code id: tags:
  
 ``` python
 K = sp.symbols(r'K', positive=True )
 z = sp.symbols(r'z')
 Z_ = K * z
 dot_z, dot_Z = sp.symbols(r'\dot{z}, \dot{Z}')
 dot_Z_ = K * dot_z
 # gamma = sp.symbols(r'\gamma', positive=True)
 # alpha = sp.symbols(r'alpha')
 # X_ = gamma * alpha
 # dot_alpha, dot_X = sp.symbols(r'\dot{\alpha}, \dot{X}')
 # dot_X_ = gamma * dot_alpha
 {dot_Z: dot_Z_, dot_tau: dot_tau_} # , dot_X: dot_X_}
 ```
  
 %% Cell type:markdown id: tags:
  
 ### Evolution equations
 As in case of ideal plasticity, the yielding process is controlled by the plastic multiplier $\lambda$ which must be non-negative. Following the same arguments as above, the  hardening slip $z$ is assumed to evolve in a normal direction with respect to the yield surface, i.e.
 \begin{align}
 \frac{\mathrm{d} f }{ \mathrm{d} Z} = -1
 \end{align}
 To let the hardening slip $\dot{z}$ grow in a positive direction we introduce the evolution equation as
 \begin{align}
 \dot{z} &= - \lambda  \frac{\mathrm{d} f }{ \mathrm{d} Z} = \lambda.
 \end{align}
 The other state variable $s_\mathrm{pl}$ is obtained using the same arguments as
 \begin{align}
 \dot{s}_\mathrm{pl} &= \lambda \frac{\partial f}{\partial \tau} = \lambda \, \mathrm{sign}(\tau)
 \end{align}
  
 %% Cell type:markdown id: tags:
  
 **Side remark:** <font color='blue'>
 The choice of the sign in the evolution equation might seem somewhat arbitrary but it is not. If a hardening stress variable, $Z$ in this case, appears in the yield condition with a negative sign we need to take the minus sign in the evolution equation to let the hardening slip grow positive consistently with the growth of $\lambda$.
 </font>
 <a id="evolution_equations"></a>
  
 %% Cell type:code id: tags:
  
 ``` python
 dot_s_pl_ = lambda_ * f_tau_Z_.diff(tau)
 dot_z_ = - lambda_ * f_tau_Z_.diff(Z)
 #dot_alpha_ = - lambda_ * f_tau_Z_.diff(X)
 # let us show the result as a dictionary to see associate the obtained expression to the coresponding state variable
 {dot_s_pl: dot_s_pl_, dot_z: dot_z_} #, dot_alpha: sp.simplify(dot_alpha_) }
 ```
  
 %% Cell type:markdown id: tags:
  
 ### Consistency condition
 The chain rule applied to express $\dot{f}$ needed for the consistency condition renders an expression which constraints the admissible evolution of both $\dot{s}_\mathrm{pl}$ and $\dot{z}$
 \begin{align}
 \dot{f} =
     \frac{\partial f}{\partial \tau} \dot{\tau}
 +
    \frac{\partial f}{\partial Z} \dot{Z}
    \label{eq:f_isotropic}
 \end{align}
  
 %% Cell type:code id: tags:
  
 ``` python
 dot_f = sp.symbols(r'\dot{f}')
 dot_f_tau_Z_ = f_tau_Z_.diff(tau) * dot_tau_ + f_tau_Z_.diff(Z) * dot_Z_ # + + f_tau_Z_.diff(X) * dot_X_
 dot_f_tau_Z_ = sp.simplify(dot_f_tau_Z_)
 {dot_f : dot_f_tau_Z_}
 ```
  
 %% Cell type:markdown id: tags:
  
 where the rates $\dot{s}_\mathrm{pl}$ and $\dot{z}$ can be substituted using the [evolution equations](#evolution_equations) which after simplification yields
  
 %% Cell type:code id: tags:
  
 ``` python
 dot_f_lambda_ = dot_f_tau_Z_.subs(dot_s_pl, dot_s_pl_).subs(dot_z, dot_z_) # .subs(dot_alpha, dot_alpha_)
 {dot_f : sp.simplify(dot_f_lambda_)}
 ```
  
 %% Cell type:markdown id: tags:
  
 The only unkown in this expression is now the plastic multiplier $\lambda$. By setting $\dot{f}$, i.e. `dot_f_lambda_`,  to zero we can finally express $\lambda$ as
  
 %% Cell type:code id: tags:
  
 ``` python
 lambda_solved = sp.solve(dot_f_lambda_, lambda_)[0] # the solution is returned as a list with one entry, therefore [0]
 sp.simplify(lambda_solved)
 ```
  
 %% Cell type:markdown id: tags:
  
 ### Solution
  
 The solved $\lambda$ can be substituted into the [evolution equations](#evolution_equations) to obtain the direct relation between the rate of state variables and the change of the control displacement $\dot{s}$
 \begin{align}
 \dot{s}_\mathrm{pl} &= \dfrac{E_\mathrm{b} \dot{s}}{ E_\mathrm{b} + K }
 \end{align}
 \begin{align}
 \dot{z} = \dfrac{E_\mathrm{b} \dot{s}}{ E_\mathrm{b} + K } \mathrm{sign}(\tau)
 \end{align}
  
 %% Cell type:code id: tags:
  
 ``` python
 dot_s_pl_solved = dot_s_pl_.subs(lambda_, lambda_solved)
 dot_z_solved = dot_z_.subs(lambda_, lambda_solved)
 # dot_alpha_solved = dot_alpha_.subs(lambda_, lambda_solved)
 {dot_s_pl: sp.simplify(dot_s_pl_solved), dot_z: sp.simplify(dot_z_solved)}
 ```
  
 %% Cell type:markdown id: tags:
  
 The rate form of the elastic constitutive relation can be obtained by substituting the plastic slip evolution to render
 \begin{align}
 \dot{\tau} = \dfrac{E_\mathrm{b} K \dot{s}}{ E_\mathrm{b} + K }
 \end{align}
  
 %% Cell type:code id: tags:
  
 ``` python
 dot_tau_solved = dot_tau_.subs(dot_s_pl, dot_s_pl_solved)
 {dot_tau: sp.simplify(dot_tau_solved)}
 ```
  
 %% Cell type:markdown id: tags:
  
 ![image.png](attachment:image.png)
  
 ### **Summary:**
 - The yielding process including the growth of the elastic domain is now uniquely described by
 the evolution of two state variables $\dot{s}_\mathrm{pl}$ and $\dot{z}$ depending on the rate of prescribed slip $\dot{s}$.
 - We have derived a direct relation between the rate of stress and the rate of slip using differentials.</br>
 <font color="darkblue">To simulate the behavior for a particular loading history we have to integrate these differential rate equations over the whole time line.</font>
 - In a way, we have just mathematically rediscovered the wheel:</br> <font color="red"> It is impossible to predict something about the future without accounting for what happened in the past.</a>
 - Thus, the question remains:</br>
 <font color="green">**How to transform the differential continuous evolution equations into an algorithm that can simulate the material response to variable loading?**</font>
  
 %% Cell type:markdown id: tags:
  
 # **Numerical iterative solution**
  
 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=615712) part 5 -->
  
 %% Cell type:code id: tags:
  
 ``` python
 YouTubeVideo('hfP8hcqWfbg')
 ```
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/evaluate.png" alt="Evaluate" width="40" height="40">
    &nbsp; &nbsp; <b>How to get numbers?</b> </div>
  
 %% Cell type:markdown id: tags:
  
 To move through an inelastic space of a material, let us now consider a discrete instance of time $t_n$ with the history represented by known values of $s_{n}$ and $s^{\mathrm{pl}}_{n}$ and $z_n$.
  
 %% Cell type:markdown id: tags:
  
 ## Discrete integration of evolution equations
  
 %% Cell type:markdown id: tags:
  
 ### Relaxing the yield condition
 In a continuous case we used the consistency condition to explicitly glue the state onto the yield surface
 \begin{align}
 \dot{f}(\tau(s, s_\mathrm{pl}(\lambda), z(\lambda)) &= 0.
 \end{align}
 Thus, it was impossible to reach an inadmissible state beyond the yield locus. In discrete case, we have to relax this requirement and allow intermediate steps into the inadmissible region during iterations.
  
 %% Cell type:markdown id: tags:
  
 To introduce discrete time stepping procedure, let us prescribe an increment of total control slip $\Delta s$ to achieve the state at $t_{n+1}$ as
 \begin{align}
 s_{n+1} = s_n + \Delta s.
 \end{align}
 Since the state variables $s^\mathrm{pl}_{n+1}, z_{n+1}$ are unknown, we start the search for an admissible state by evaluating the yield function with the values known from the previous step, which will result in a step into an inadmissible range
 \begin{align}
 f(s_{n+1}, s^{\mathrm{pl}}_n, z_n) > 0
 \end{align}
 ![image.png](attachment:image.png)
  
 %% Cell type:markdown id: tags:
  
 Indeed, by taking the state variables $s^\mathrm{pl}_n, z_n$ from the last step as a first estimate, the yield function $f(s_{n+1}, s^{\mathrm{pl}}_n, z_n)$ can deliver positive values in this trial iteration step.
  
 **The "trial" step beyond the admissible domain $f \le 0$ must be accompanied with a "return mapping" algorithm that iteratively returns back to an admissible state located on the yield surface.**
  
 %% Cell type:markdown id: tags:
  
 ### Return mapping using discrete consistency condition
  
 %% Cell type:markdown id: tags:
  
 Given an inadmissible trial state $k$ violating the yield condition, i.e. $f_k > 0$, let us introduce a linearized approximation of its change along the plastic multiplier $\lambda$ around the iteration state $k$.
 \begin{align}
 f_{k+1} &= f_{k} + \left. \frac{\partial f}{\partial \lambda} \right|_k \Delta \lambda
 \end{align}
 In this form, we can reach an admissible state of yielding recovering the condition $f_{n+1} = 0$ by iterating over $k$ to identify the value of the plastic multiplier $\lambda$.
 Note that in initial iteration $k = 0$ the state from previous step is reused, i.e. $f(s_{n+1}, s_n^\mathrm{pl}, z_n)$.
  
 %% Cell type:markdown id: tags:
  
 In the linearized form, we can transform the yield condition to a recurrent formula
 \begin{align}
 \left. \frac{\mathrm{d} f}{\mathrm{d} \lambda}\right|_k \Delta \lambda &= -f_k,
 \hspace{1cm} f_k \rightarrow 0 \; \;\mathrm{for}\;\; k = 1\ldots\infty
 \end{align}
 This is another application of the [Newton method](../extras/newton_method.ipynb) for iterative solution of a nonlinear equation visited in **Tour 2**.
  
 %% Cell type:markdown id: tags:
  
 In the visual representation of the return mapping algorithm, shown in the figure,
 we can recognize the process of approaching the zero level of the yield function $f$.
  
 ![image.png](attachment:image.png)
  
 Note that the predictor, i.e. the slope of the yield condition, is negative and $\Delta \lambda > 0$. In every step, the plastic multiplier is updated:
 \begin{align}
 \lambda_{k+1} &= \lambda_k + \Delta \lambda, \, \lambda_0 = 0 \nonumber \\ \nonumber
 \end{align}
 As we will see in the simple implementation below, the admissible state is reached in just one step in case of the considered linear isotropic hardening.
 In more general plasticity models and two- or three-dimensional stress states, an iterative return mapping can require several steps and the level of admissibility expressed by the yield function change as well during the return mapping.
  
 %% Cell type:markdown id: tags:
  
 **We have the iterative scheme, let us derive its components:** Two more questions must be addressed to define a general numerical algorithm for plasticity:
 <font color="brown">
 * **Update of state variables $s^\mathrm{pl}_{k+1}$ and $z_{k+1}$ in each iteration**
 * **Expression of the predictor $\mathrm{d} f / \mathrm{d} \lambda$ in terms of the state variables**
 </font>
  
 %% Cell type:markdown id: tags:
  
 ### State update using evolution equations
 In every iteration step the state variables $s_\mathrm{pl}$ and $z$ must be updated using the discrete [evolution equations](#evolution_equations), i.e.
  
 <a id="discrete_evolution_equations"></a>
 \begin{align}
 s^\mathrm{pl}_{k+1} &= s^\mathrm{pl}_n + \lambda_{k+1}
 \left. \frac{\partial f}{\partial \tau} \right|_k \nonumber \\
 z_{k+1} &= z_n - \lambda_{k+1} \left. \frac{\partial f}{\partial Z} \right|_k
 \label{eq:discrete_evolution}
 \end{align}
  
 %% Cell type:markdown id: tags:
  
 ### Predictor operator
 Recalling that $f(\tau(s,s^{\mathrm{pl}}(\lambda), z(\lambda)))$ the chain rule delivers the expression
  
 \begin{align}
 \left. \frac{\partial f}{\partial \lambda} \right|_k
 &=
 \left.\frac{\partial f}{\partial \tau} \right|_k
 \left.\frac{\partial \tau}{\partial s^{\mathrm{pl}}} \right|_k
 \left.\frac{\partial s^{\mathrm{pl}}} {\partial \lambda} \right|_k
 +
 \left.\frac{\partial f}{\partial Z} \right|_k
 \left.\frac{\partial Z}{\partial z} \right|_k
 \left.\frac{\partial z}{\partial \lambda} \right|_k
 \end{align}
  
 after substituting the derivatives of the constitutive laws for [elastic behavior](#elastic_law) and [isotropic hardening](#hardening_law) w.r.t. $s^{\mathrm{pl}}$ and $z$, respectively and, of the [discrete evolution equations](#discrete_evolution_equations)  w.r.t. $\lambda$ we obtain
  
 \begin{align}
 \left. \frac{\partial f}{\partial \lambda} \right|_k
 &= -
 \left. \frac{\partial f}{\partial \tau} \right|_k
 E_\mathrm{b}
 %\left. \frac{\partial \tau}{\partial s^{\mathrm{pl}}} \right|_k
 \left. \frac{\partial f}{\partial \tau} \right|_k -
 \left. \frac{\partial f}{\partial Z} \right|_k
 K
 %\left. \frac{\partial Z}{\partial z} \right|_k
 \left. \frac{\partial f}{\partial Z} \right|_k
 \end{align}
 which can be finally expressed as
 \begin{align}
 \left. \frac{\partial f}{\partial \lambda} \right|_k
 &= -
 \left. \frac{\partial f}{\partial \tau} \right|_k^2 E_\mathrm{b} -
 \left. \frac{\partial f}{\partial Z} \right|_k^2 K = -(E_\mathrm{b} + K)
 \end{align}
  
 %% Cell type:markdown id: tags:
  
 ### Time stepping algorithm
 Substituting the obtained predictor back into the recurrent formula we obtain the solution for $\Delta \lambda$
  
 \begin{align}
 f_k + \left. \frac{\partial f}{\partial \lambda} \right|_k \Delta \lambda =
 f_k - (E_\mathrm{b} + K) \Delta \lambda = 0
 \implies
 \Delta \lambda = \frac{f_k}{E_\mathrm{b}+K}
 \end{align}
  
 Apparently, the derivative of $f$ with respect to $\lambda$ is constant in the present model. This means that the solution can be found in a single iteration step. This gives us the chance to derive an explicit analytical formulas for return mapping in a time step $s_{n+1} = s_n + \Delta s$ using the state variables $s^\mathrm{pl}_n, z_n$ from the last time step as follows:
 <div style="background-color:#dfffef;text-align:left">
    <a id='return_mapping'></a>
 Given $s_{n+1}, s^{\mathrm{pl}}_n, z_n$
 \begin{align}
 \tau_{k} &= E_b(s_{n+1} - s^{\mathrm{pl}}_n) \nonumber \\
 Z_k &= K z_n \\
 f_k &= | \tau_k | - Z_k - \tau_{\mathrm{Y}} \nonumber
 \end{align}
 if $f_k > 0$
 \begin{align}
 \Delta \lambda &= \frac{f_k}{E_\mathrm{b} + K} \\
 s^\mathrm{pl}_{n+1} &= s^\mathrm{pl}_{n} + \Delta \lambda \; \mathrm{sign}(\tau_k)
 \nonumber \\
 z_{n+1} &= z_{n} + \Delta \lambda \nonumber \\
 \tau_{n+1} &= E_\mathrm{b} (s_{n+1} - s^\mathrm{pl}_{n+1})
 \end{align}</font></div>
  
 %% Cell type:markdown id: tags:
  
 <a id="trilinear_material_model"></a>
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/work.png" alt="Coding intermezzo" width="40" height="40">
    &nbsp; &nbsp; <b>Coding intermezzo:</b> we have an algorithm - let us implement it </div>
  
 %% Cell type:markdown id: tags:
  
 ## Example implementation
  
 %% Cell type:markdown id: tags:
  
 ### Cyclic loading history
  
 Consider a loading prescribed in terms of the slip function $s(t)$ in the form
 \begin{align}
 s(t) = s_\max \, \theta(t)
 \end{align}
 where $s_\max$ represents a prescribed amplitude and $\theta(t)$ a time history expressed by a $\sin$ function as
 \begin{align}
 \theta(t) = \sin(2 n_\mathrm{cycles} \, \pi \, t)
 \end{align}
  
 %% Cell type:markdown id: tags:
  
 The next cell delivers an array with 10 cycles, i.e. $n_\mathrm{cycles} = 10$, with an amplitude of $s_\max = 3$ with 500 data points, i.e. `n_steps` = 500.
  
 %% Cell type:code id: tags:
  
 ``` python
 import numpy as np
 n_cycles, s_max, n_steps = 10, 3, 500 # load history parameters
 t_arr = np.linspace(0,1,n_steps) # time range t in (0,1)
 theta = np.sin(2*n_cycles * np.pi * t_arr) # load history with unloading
 s_n1_arr = s_max * theta # load history
 s_n1_arr[:10] # show just the 10 first entries in the array
 ```
  
 %% Cell type:code id: tags:
  
 ``` python
 fig, ax = plt.subplots(1,1, figsize=(10,3))
 ax.plot(t_arr, s_n1_arr)
 ```
  
 %% Cell type:markdown id: tags:
  
 With this loading history at hand and assuming the material parameters. To distinguish the numerical parameters from the symbols above, a prefix underscore is used in the numerical code.
 | parameter description | mathematical symbol | variable name | value | unit |
 | - | - | - | - | - |
 | Elastic stiffness | $E_\mathrm{b}$ | `_E_b` | 10 | MPa |
 | Hardening modulus | $K$ | `_K` | 0.2 | MPa |
 | Yield stress | $\tau_Y$ | `_tau_Y` | 1 | MPa |
  
 %% Cell type:markdown id: tags:
  
 ### Example of the time integration loop
  
 The following loop is evaluating the admissible state for all entries from the load history `s_n1_arr`. First, the trial state is evaluated using the values of `s_pl_k` and `z_k` from the previous step. If the yield condition is violated, the [return mapping](#return_mapping) formulas are evaluated. The resulting stress is appended to the recorded stress history in `tau_list`
 <a id="implementation_return_mapping"></a>
  
 %% Cell type:code id: tags:
  
 ``` python
 _E_b, _K, _tau_Y = 1, 0.2, 1 # material parameters
 tau_list = [] # list to record the stresses
 s_pl_k, z_k = 0, 0 # initialization of trial states
 for s_n1 in s_n1_arr:
    ### START: Material model implementation as a user subroutine in a finite element code
    tau_k = _E_b * (s_n1 - s_pl_k) # elastic trial step
    Z_k = _K * z_k # isotropic hardening
    f_k = np.abs( tau_k ) - (_tau_Y + Z_k)
    if f_k > 0: # inelastic step - return mapping
        delta_lambda_k = f_k / (_E_b + _K)
        s_pl_k += delta_lambda_k * np.sign(tau_k)
        z_k += delta_lambda_k # to save lines n=n+1 is shortend to k
        tau_k = _E_b * (s_n1 - s_pl_k)
    ### END: Material model implementation as a user subroutine in a finite element code
    tau_list.append(tau_k) # record the calculated stress
 ```
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/view.png" alt="Evaluate" width="40" height="40">
    &nbsp; &nbsp; <b>Plot the loading history and stress-slip history</b> </div>
  
 %% Cell type:code id: tags:
  
 ``` python
 %matplotlib widget
 fig, (ax_t, ax_tau) = plt.subplots(1,2,figsize=(10,5))
 fig.canvas.header_visible = False
 ax_t.plot(t_arr, s_n1_arr,color='black'); ax_t.set_ylabel(r'$s = s_\max \theta(t)$'), ax_t.set_xlabel('$t$')
 ax_tau.plot(s_n1_arr, tau_list, color='green'); ax_tau.set_ylabel(r'$\tau$'), ax_tau.set_xlabel('$s$');
 ax_t.set_title('loading history'); ax_tau.set_title('stress-slip');
 ```
  
 %% Cell type:markdown id: tags:
  
 ### Material model implemented
  
 %% Cell type:markdown id: tags:
  
 **Almost a new material model:** The code within the loop represents an implementation of a material model within a non-linear finite-element code accepting the imposed value of strain (slip) $s_{n+1}$ denoted as `s_n1` in the code and returning the corresponding admissible value of stress (shear) `tau`.
  
 In fact, one more step is missing. To embed the new material model into a finite-element calculation, we would also need to return the material stiffness relevant for the time $t_{n+1}$. Recall, that in the notebook [3.1 - nonlinear bond](../tour3_nonlinear_bond/3_1_nonlinear_bond.ipynb#trilinear_material_model) we already explained that an efficient time stepping requires not only the value of the residuum, i.e. a spatial integral of stresses, but also its derivatives, i.e. the derivative of stress w.r.t. strains or slip in our case. Thus, to complete the task of providing a plastic model with hardening, we need express the derivative
 \begin{align}
 E_{\mathrm{b},n+1} &=
 \dfrac{\mathrm{d} \tau_{n+1}}{\mathrm{d} s_{n+1}}
 \end{align}
 This expression is called algorithmic stiffness and can be obtained by resolving the chain derivatives in the [solved return mapping algorithm](#return_mapping)
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/enhancement.png" alt="Question" width="50" height="50">
    &nbsp; &nbsp; <b>DIY extension</b> </div>
  
 %% Cell type:markdown id: tags:
  
 **Task:** To show that the derivation procedure is general, let us now extend the with kinematic hardening $X = \gamma \alpha$. The following video shows how to do it.
  
 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=615712) part 6 -->
  
 %% Cell type:code id: tags:
  
 ``` python
 YouTubeVideo('WROT3aN-BR0')
 ```
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/question_puzzle.png" alt="Question" width="50" height="50">
    &nbsp; &nbsp; <b>Stimulating questions</b> </div>
  
 %% Cell type:markdown id: tags:
  
 1. Use the model with both isotropic and kinematic hardening at the [top of this notebook](#plastic_explorer) to find and extend the answers to questions in exercise [X0401]
 2. What behavior is obtained for one loading cycle with $K = 0.1, \gamma = -0.1$ and for $K = -0.1, \gamma = 0.1$ - for monotonic and for cyclic loading?
 3. Use the model to find an explanation why do the interchanged signs for kinematic and isotropic hardening lead to a completely different cyclic response.
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/exercise.png" alt="Run" width="50" height="50">
    &nbsp; &nbsp; <a href="../exercises/X0401 - Bond-slip model classification (unloading and reloading).pdf"><b>Exercise X0401:</b></a> <b>Bond-slip model classification (unloading and reloading)</b>
 </div>
  
 %% Cell type:code id: tags:
  
 ``` python
 YouTubeVideo('1M0SmlF2bFI')
 ```
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left"> <img src="../icons/exercise.png" alt="Run" width="40" height="40">
    &nbsp; &nbsp; <a href="../exercises/X0402 - Nonlinear bond behavior modeled using plasticity.pdf"><b>Exercise X0402:</b></a> <b>Nonlinear bond behavior modeled using plasticity</b>
 </div>
  
 %% Cell type:markdown id: tags:
  
 <div style="background-color:lightgray;text-align:left;width:45%;display:inline-table;"> <img src="../icons/previous.png" alt="Previous trip" width="50" height="50">
    &nbsp; <a href="4_1_PO_multilinear_unloading.ipynb#top">4.1: Loading, unloading and reloading</a>
 </div><div style="background-color:lightgray;text-align:center;width:10%;display:inline-table;"> <a href="#top"><img src="../icons/compass.png" alt="Compass" width="50" height="50"></a></div><div style="background-color:lightgray;text-align:right;width:45%;display:inline-table;">
    <a href="4_3_PO_trc_cfrp_cyclic.ipynb#top">4.3 Pullout and cyclic loading</a>&nbsp; <img src="../icons/next.png" alt="Previous trip" width="50" height="50"> </div>
  
 %% Cell type:code id: tags:
  
 ``` python
 ```
--- a/tour4_plastic_bond/4_3_PO_trc_cfrp_cyclic.ipynb
+++ b/tour4_plastic_bond/4_3_PO_trc_cfrp_cyclic.ipynb
@@ -199,7 +199,7 @@
   "source": [
    "## Discussion of the study\n",
    "\n",
-    " - The predicted response reproduces the experimental response reasonably well. For the **short bond length** used in the test, the slip and shear profiles are nearly uniform along the bond zone.  - **Task:** Change the bond length to 200 mm and compare the response with the pull-out curve depicted in the [test results](#trc_test_results)\n",
+    " - The predicted response reproduces the experimental results reasonably well. For the **short bond length** used in the test, the slip and shear profiles are nearly uniform along the bond zone.  - **Task:** Change the bond length to 200 mm and compare the response with the pull-out curve depicted in the [test results](#trc_test_results)\n",
    " - The main focus of the study at hand was the question, if the **unloading stiffness** corresponds to the initial stiffness, similarly to the test response. Since this is the case, the conclusion can be made, that the bond behavior in the tested TRC specimens is primarily governed by plasticity.\n",
    " - Since the model uses **linear hardening**, the pull-out response has a linear inelastic branch up to the point of unloading. The non-linear shape of the experimentally obtained curve cannot be reproduced in a better way. **Non-linear hardening** function instead of a single hardening modulus $\\gamma$ can be theoretically included into the framework of plasticity. However, such enhancement requires a non-trivial systematic calibration procedure which would be able to uniquely associate the inelastic mechanisms governing the response to the hardening process, which is a non-trivial task. Such effort would only be justified for a material combination with large amount of practical applications. "
   ]

 %% Cell type:markdown id: tags:

 <a id="top"></a>

 %% Cell type:markdown id: tags:

 # **4.3: Cyclic pullout of textile fabrics and CFRP sheets**

 %% Cell type:markdown id: tags:

 <!-- [![title](../fig/bmcs_video.png)](https://moodle.rwth-aachen.de/mod/page/view.php?id=619022) -->

 %% Cell type:code id: tags:

 ``` python
 from IPython.display import YouTubeVideo
 YouTubeVideo('ItmQQ9yQDjI')
 ```

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left"> <img src="../icons/start_flag.png" alt="Previous trip" width="40" height="40">
    &nbsp; &nbsp; <b>Starting point</b> </div>

 %% Cell type:markdown id: tags:

 Once we have seen that there is a relation between the shape of the constitutive law and the  stress redistribution process within a bond zone which has an immediate consequence on the shape of the pullout curve, we extend our horizon to the case of **non-monotonic loading**.

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left"> <img src="../icons/destination.png" alt="Previous trip" width="40" height="40">
    &nbsp; &nbsp; <b>Where are we heading</b> </div>

 %% Cell type:markdown id: tags:


 To motivate a more general development of the model for the bond-slip behavior let us consider the case of non-monotonic loading. What happens in the material structure of the bond if the load is reduced and than it grows again?

 %% Cell type:markdown id: tags:

 <a id="trc_test"></a>
 # **TRC pullout test modeled using plastic material model**

 %% Cell type:markdown id: tags:

 ## Test results

 %% Cell type:markdown id: tags:

 <a id="trc_test_results"></a>
 Consider again the test double-sided pullout test that was used to motivate the development of plastic material model in [notebook 4.1](4_1_PO_multilinear_unloading.ipynb#trc_study_monotonic) which demonstrated that no unloading branch could be reproduced. We reprint the test results for convenience here with the goal to examine the ability of the plastic model to reproduce the unloading stiffness of the specimen with the total length of 300 mm, i.e. $L_\mathrm{b} = 150$ mm.

 ![image](../fig/test_unloading.png)

 %% Cell type:markdown id: tags:

 ## Pullout model setup

 Let us first construct the model with the same geometry as specified in [notebook 3.2](../tour3_nonlinear_bond/3_2_anchorage_length.ipynb#trc_parameters). Note again that the maximum control displacement $w_\mathrm{max}$ is
 the half of the crack opening displacement measured in the test. Therefore `w_max` has been set to 3 mm.

 %% Cell type:code id: tags:

 ``` python
 %matplotlib widget
 import matplotlib.pyplot as plt
 import numpy as np
 from bmcs_cross_section.pullout import PullOutModel1D
 po_trc = PullOutModel1D(n_e_x=100, w_max=3) # mm
 po_trc.geometry.L_x=150 # [mm]
 po_trc.time_line.step = 0.02
 po_trc.cross_section.trait_set(A_m=1543, A_f=16.7, P_b=10);
 ```

 %% Cell type:markdown id: tags:

 ## Material model with kinematic hardening

 To apply the elastic-plastic model in this pullout model, set the the `elasto-plasticity` option for the `material_model` attribute. Then, the plastic model with the parameters $\tau_Y = \bar{\tau}, K, \gamma$ introduced in the [notebook 3.2](../tour3_nonlinear_bond/3_2_anchorage_length.ipynb#trc_pullout_study) can be set to the attribute `material_model_` using the `trait_set` method as follows

 | parameter | name | value | unit |
 | - | - | -: | - |
 | bond stiffness | `E_b` | 6.4 | MPa |
 | elastic limit stress | `tau_bar` | 4 | MPa |
 | isotropic hardening modulus | `K` | 0 | MPa |
 | kinematic hardening modulus | `gamma` | 1.0 | MPa |

 %% Cell type:markdown id: tags:

 Using these parameters, we aim to reproduce the same shape of the [bond slip law that was used previously in notebook 3.2](../tour3_nonlinear_bond/3_2_anchorage_length.ipynb#trc_parameters), represented as a tri-linear function for the simulation of the monotonically increasing load.

 %% Cell type:code id: tags:

 ``` python
 po_trc.material_model='elasto-plasticity'
 po_trc.material_model_.trait_set(
    E_m=28000, E_f=170000, E_b=6.4, K = 0, gamma=1.0, tau_bar=4
 );
 po_trc.material_model_.interact()
 ```

 %% Cell type:markdown id: tags:

 ## Load scenario

 Finally, following the same mimic, we change the loading scenario from `monotonic` to `cyclic` and can set the parameters of the cyclic loading through the attribute with the trailing underscore `loading_scenario_`.
 To verify the loading scenario we directly render the model component using the `interact` method.

 %% Cell type:code id: tags:

 ``` python
 po_trc.loading_scenario.profile = 'cyclic-nonsym-const'
 po_trc.loading_scenario.profile_.trait_set(number_of_cycles=2,
                              unloading_ratio=0.8,
                            number_of_increments=200,
                              );
 po_trc.loading_scenario.interact()
 ```

 %% Cell type:markdown id: tags:

 ## Run the simulation
 Apparently, we apply directly the full control displacement, then reduce it to 80% and reload again to the maximum. Our intention is not to simulate the whole history displayed in the [test results](trc_test_results) but to compare the unloading stiffness reproduced by the model.
 The model is now ready for execution so that let us run it and render the user interface to inspect the results.

 %% Cell type:code id: tags:

 ``` python
 po_trc.reset()
 po_trc.run()
 po_trc.interact()
 ```

 %% Cell type:markdown id: tags:

 ## Discussion of the study

- - The predicted response reproduces the experimental response reasonably well. For the **short bond length** used in the test, the slip and shear profiles are nearly uniform along the bond zone.  - **Task:** Change the bond length to 200 mm and compare the response with the pull-out curve depicted in the [test results](#trc_test_results)
+ - The predicted response reproduces the experimental results reasonably well. For the **short bond length** used in the test, the slip and shear profiles are nearly uniform along the bond zone.  - **Task:** Change the bond length to 200 mm and compare the response with the pull-out curve depicted in the [test results](#trc_test_results)
 - The main focus of the study at hand was the question, if the **unloading stiffness** corresponds to the initial stiffness, similarly to the test response. Since this is the case, the conclusion can be made, that the bond behavior in the tested TRC specimens is primarily governed by plasticity.
 - Since the model uses **linear hardening**, the pull-out response has a linear inelastic branch up to the point of unloading. The non-linear shape of the experimentally obtained curve cannot be reproduced in a better way. **Non-linear hardening** function instead of a single hardening modulus $\gamma$ can be theoretically included into the framework of plasticity. However, such enhancement requires a non-trivial systematic calibration procedure which would be able to uniquely associate the inelastic mechanisms governing the response to the hardening process, which is a non-trivial task. Such effort would only be justified for a material combination with large amount of practical applications.

 %% Cell type:markdown id: tags:

 <a id="cfrp_test"></a>
 # **CFRP sheet test with softening plastic material**

 It is interesting to examine the ability of the plastic model with negative hardening, i.e. softening, to simulate the cyclic response of the CFRP sheets.

 ## Pullout model setup
 Let us construct the model with these parameters and check the shape of the bond-slip relation

 %% Cell type:code id: tags:

 ``` python
 %matplotlib widget
 import matplotlib.pyplot as plt
 import numpy as np
 from bmcs_cross_section.pullout import PullOutModel1D

 A_f = 16.67 # [mm^2]
 A_m = 1540.0 # [mm^2]
 p_b = 100.0 #
 E_f = 170000 # [MPa]
 E_m = 28000 # [MPa]
 pm = PullOutModel1D()
 pm.sim.tline.step = 0.01 # 100 time increments
 pm.cross_section.trait_set(A_f=A_f, P_b=p_b, A_m=A_m)
 pm.geometry.L_x = 300 # length of the specimen [mm]
 pm.w_max = 2.8 # maximum control displacement [mm]
 pm.n_e_x = 100 # number of finite elements
 ```

 %% Cell type:markdown id: tags:

 ## Material model with isotropic softening

 The bond-slip law identified during the [**Tour 3:** (CFRP sheet test)](../tour3_nonlinear_bond/3_1_nonlinear_bond.ipynb#cfrp_bond_slip) can be reproduced by setting the parameters of
 the plastic material model with isotropic hardening to:
 | parameter | name | value | unit |
 | - | - | -: | - |
 | bond stiffness | `E_b` | 80 | MPa |
 | elastic limit stress | `tau_bar` | 8 | MPa |
 | isotropic hardening modulus | `K` | -20 | MPa |
 | kinematic hardening modulus | `gamma` | 0 | MPa |

 %% Cell type:code id: tags:

 ``` python
 pm.material_model='elasto-plasticity'
 pm.material_model_.trait_set(
    E_m=28000, E_f=170000, E_b=80, K = -20, gamma=0, tau_bar=8
 )
 pm.material_model_.interact()
 ```

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left"> <img src="../icons/remark.png" alt="Previous trip" width="40" height="40">
    &nbsp; &nbsp; <b>What is the difference compared to a trilinear bond-slip law used before?</b> </div>

 %% Cell type:markdown id: tags:

 Apparently, the shape is identical to the [trilinear bond slip law](../tour3_nonlinear_bond/3_1_nonlinear_bond.ipynb#cfrp_trilinear_bond) law applied in the first analysis of this test setup. However, now unloading in inelastic regime would be governed by initial stiffness.
 It is important to remark that the model has been augmented with an additional feature that prevents softening into negative stresses. To explain this minor but essential extension recal the shape of the yield condition for plastic model with isotropic hardening
 \begin{align}
 f = | \tau | - (\tau_Y - Z)
 \end{align}
 where the hardening stress is given as $Z = K z$. With the negative value of hardening modulus $K$ and implicitly positive value of the hardening state variable $z$ we can notice that the size of the elastic domain $(\tau_Y - Z)$ shrinks with increasing $z$. Without any treatment, the elastic domain could reach a negative size, which is nonphysical. Thus, the criterion of the positive size of the elastic domain is used to define a third branch in the plastic model.

 **In fact,** with this necessary extension, we introduce **damage into a plastic model**, because we explicitly set the material stiffness at a predefined value of deformation, namely at the value of slip inducing the loss of elastic domain. As we will lear in Tour 5, relating the loss of stiffness to the measure of strain is the general concept behind damage models. :

 In this context, where we introduced an argument of physical admissibility, it is worth reminding that inelastic behavior implicitly induces **dissipation of energy**. A sound formulation of material models is only possible in the framework of thermodynamic laws that can guarantee an admissible behavior for any loading history. The evaluation of energy dissipation and its importance for objectivity and robustness of material models will be illuminated in **Tour 6**.

 With the material model at hand, let us explore the response of the CFRP sheets in a cyclic loading scenario with symmetric, increasing amplitudes

 %% Cell type:markdown id: tags:

 ## Load scenario

 To see the effect of propagating softening plasticity upon unloading, let us apply four loading cycles oscillating applying the pull-out and push-in load in each cycle. Note that with this,
 a compressive stress is transferred by the FRP sheet upon push-in loading, which is unrealistic.
 Such response might be expected from a test on CFRP bars. Still the purpose of the study is
 to show the consequences of the applied assumptions on a broad scale of loading configuration and check the plausibility of the model response for a wider range of loading scenarios.

 %% Cell type:code id: tags:

 ``` python
 pm.loading_scenario.profile = 'cyclic-sym-incr'
 pm.loading_scenario.profile_.trait_set(number_of_cycles=2,
                              unloading_ratio=0.5,
                            number_of_increments=200,
                              );
 pm.reset()
 pm.run()
 pm.interact()
 ```

 %% Cell type:markdown id: tags:

 ## Discussion of the study

 - By checking the `history` showing the field variables along the bond zone, we can observe the alternating tension and compression profiles in the reinforcement and in the matrix.
 - The shear transfer is again taking place only on a short length due to the softening nature of the bond behavior.
 - The yielding of the bond within the process zone occurs both under pullout and under push-in load as apparent from the existence of the hysteretic loop.
 - In contrast to the [study of TRC test](#trc_results), the unloading stiffness does not correspond to the initial stiffness. This question is discussed in a more detail below.

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left"> <img src="../icons/question_puzzle.png" alt="Question" width="50" height="50">
    &nbsp; &nbsp; <b>Stimulating question</b> </div>

 %% Cell type:markdown id: tags:

 # **Qualitative comparison of the two studies**

 %% Cell type:markdown id: tags:

 ## **How to interpret unloading stiffness at structural level?**

 In the classification of plasticity and damage models explained in  [4.1](4_1_PO_multilinear_unloading.ipynb#damage_and_plasticity) we have learned that in a single material point, plastic behavior unloads with initial stiffness and induces a persistent deformation $s_\mathrm{pl}$ at zero level of load. On the other hand, unloading of material point governed by damage behavior would recover all deformation and reach zero deformation at zero stress. **Can we assume such behavior of both class of models at the structural level?**

 Comparing the slopes of the first unloading branch obtained for the [TRC test](#trc_test) and of the [CFRP test](cfrp_test), we can notice that in the former case, the unloading stiffness is indeed equal to the initial stiffness. In this case, governed by the hardening bond-slip law, the shear stress increases and decreases upon loading and unloading in all material points simultaneously. On the other hand, in the case of [CFRP test](cfrp_test) governed by softening, the debonded zone reaches zero stress level and does not contribute to stress transfer any more. As a result, the effective length of the bond zone becomes shorter, so that the structural stiffness upon unloading must be smaller than the initial stiffness.

 **The conclusion** is, that an application of a plastic model with softening,
 is only valid up to the point at which the elastic domain vanished. From that point on,
 a full damage must be assumed. Therefore, the kind of model used in the simulation of the CFRP test can be regarded as a first version of a combined damage-plasticity model.

 %% Cell type:markdown id: tags:

 ## **When to apply isotropic and when kinematic hardening?**

 In the two examples above, both kinematic and isotropic hardening has been used. In the first yielding step and upon the first unloading, there is no difference between the two parameters. They will only show a different behavior upon cyclic loading scenario. A correct distinction is only possible based on a test delivering several unloading and reloading cycles. Correct identification of the parameters is a demanding task which calls for a systematic and automated support in terms of calibration algorithms and validation procedures.

 %% Cell type:markdown id: tags:

 <div style="background-color:lightgray;text-align:left;width:45%;display:inline-table;"> <img src="../icons/previous.png" alt="Previous trip" width="50" height="50">
    &nbsp; <a href="../tour4_plastic_bond/4_2_BS_EP_SH_I_A.ipynb#top">4.2 Basic concept of plasticity</a>
 </div><div style="background-color:lightgray;text-align:center;width:10%;display:inline-table;"> <a href="#top"><img src="../icons/compass.png" alt="Compass" width="50" height="50"></a></div><div style="background-color:lightgray;text-align:right;width:45%;display:inline-table;">
    <a href="../tour5_damage_bond/5_1_Introspect_Damage_Evolution_Damage_initiation.ipynb#top">5.1 Pullout simulation using damage model</a>&nbsp; <img src="../icons/next.png" alt="Previous trip" width="50" height="50"> </div>

 %% Cell type:code id: tags:

 ``` python
 ```