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
}