Interactive Thermal Stress Profiles, Cinar Turhan.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a1aa8f44",
   "metadata": {},
   "source": [
    "# Interactive Thermal Stress Profiles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a42a0e4",
   "metadata": {},
   "source": [
    "<hr style=\"height:2px;border-width:0;color:black;background-color:black\">\n",
    "\n",
    "#### Çınar Turhan, Hildebrand Department of Petroleum and Geosystems Engineering, The University of Texas at Austin\n",
    "**[LinkedIn](https://www.linkedin.com/in/cinarturhan/) - [GitHub](https://github.com/cinarturhan)**\n",
    "\n",
    "\n",
    "#### PGE 383 Advanced Geomechanics Final Project\n",
    "#### Instructor: Dr. Nicolas Espinoza \n",
    "**[GitHub](https://github.com/dnicolasespinoza/)**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "543c3b91",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from ipywidgets import interactive\n",
    "from ipywidgets import widgets                            \n",
    "from ipywidgets import Layout\n",
    "from ipywidgets import Label\n",
    "from ipywidgets import VBox, HBox\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6f86cf09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assigning Variables:\n",
    "sigmaHmax=12 # Maximum Effective Stress [MPa]\n",
    "sigmaHmin=3 # Minimum Effective Stress [MPa]\n",
    "a=0.2 # Diameter of Interest [m]\n",
    "E=10e3 # Young's Modulus [MPa]\n",
    "PR=0.2 # Poisson's Ratio \n",
    "UCS=30 # Unconfined Compressive Strength [MPa]\n",
    "Ts=2 # Tensile Strength [MPa]\n",
    "Pw=30 # Wellbore (Mud) Pressure [MPa]\n",
    "Pp=30 # Pore Pressure [MPa]\n",
    "alpha = 1e-5 # Linear Coefficient of Thermal Expansion [1/K]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "83d89902",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmaHmax = widgets.IntSlider(min=12, max = 70, value=12, step = 7, description = '$\\sigma_{Hmax}$',orientation='horizontal',style = {'description_width': 'initial'}, continuous_update=False)\n",
    "sigmaHmin = widgets.IntSlider(min=3, max = 30, value=3, step = 3, description = '$\\sigma_{Hmin}$',orientation='horizontal',style = {'description_width': 'initial'}, continuous_update=False)\n",
    "Pw = widgets.IntSlider(min=10, max = 60, value=30, step = 10, description = '$P_w$',orientation='horizontal',style = {'description_width': 'initial'}, continuous_update=False)\n",
    "Trock = widgets.IntSlider(min=100, max = 300, value=100, step = 20, description = '$T_{rock}$',orientation='horizontal',style = {'description_width': 'initial'}, continuous_update=False)\n",
    "Tfluid = widgets.IntSlider(min=50, max = 300, value=100, step = 20, description = '$T_{fluid}$',orientation='horizontal',style = {'description_width': 'initial'}, continuous_update=False)\n",
    "l = widgets.Text(value='                                                                                   Interactive Thermal Stress Profiles, Çınar Turhan',\n",
    "                 layout=Layout(width='auto', height='auto'))\n",
    "ui = widgets.HBox([sigmaHmax,sigmaHmin,Trock,Tfluid,Pw])\n",
    "ui2 = widgets.VBox([l,ui])\n",
    "widgets.AppLayout(justify_content ='center')\n",
    "\n",
    "def Interactive_Thermal_Stress(sigmaHmax,sigmaHmin,Trock,Tfluid,Pw):\n",
    "    \n",
    "    # Defining the coordinates and dataframe for assigning stress values to each coordinate for the heatmap\n",
    "    r = np.linspace(a,3.1,300);theta = np.linspace(0,359,359);r_arr=[];theta_arr=[]\n",
    "    for r in r:\n",
    "        for th in theta:\n",
    "            r_arr.append(r);theta_arr.append(th)\n",
    "    df = pd.DataFrame();df['r']=r_arr;df['theta']=theta_arr;df['x']=r_arr*np.cos(np.deg2rad(theta_arr));df['y']=r_arr*np.sin(np.deg2rad(theta_arr))\n",
    "    df['sigmarr'] = (Pw - Pp)*(a**2/df['r']**2)+(sigmaHmax+sigmaHmin)/2*(1-(a**2/df['r']**2))+(sigmaHmax-sigmaHmin)/2*(1-(4*a**2/df['r']**2)+(3*a**4/df['r']**4))*np.cos(np.deg2rad((2*df['theta'])))\n",
    "    df['sigmatt'] = -(Pw - Pp)*(a**2/df['r']**2)+(sigmaHmax+sigmaHmin)/2*(1+(a**2/df['r']**2))-(sigmaHmax-sigmaHmin)/2*(1+(3*a**4/df['r']**4))*np.cos(np.deg2rad(2*df['theta']))\n",
    "    df['sigmatt_T'] = -(Pw - Pp)*(a**2/df['r']**2)+(sigmaHmax+sigmaHmin)/2*(1+(a**2/df['r']**2))-(sigmaHmax-sigmaHmin)/2*(1+(3*a**4/df['r']**4))*np.cos(np.deg2rad(2*df['theta']))+alpha*E/(1-PR)*(Tfluid-Trock)\n",
    "    df['sigmart'] = (sigmaHmax-sigmaHmin)/2*(1+(2*a**2/df['r']**2)-(3*a**4/df['r']**4))*np.sin(np.deg2rad(2*df['theta']))\n",
    "    \n",
    "    # Defining the coordinates and dataframe for assigning stress values to each coordinate, for heatmap\n",
    "    r = np.linspace(0.1,3.1,300);theta = np.linspace(0,359,359);r_arr=[];theta_arr=[]\n",
    "    for r in r:\n",
    "        for th in theta:\n",
    "            r_arr.append(r);theta_arr.append(th)\n",
    "    df2 = pd.DataFrame();df2['r']=r_arr;df2['theta']=theta_arr;df2['x']=r_arr*np.cos(np.deg2rad(theta_arr));df2['y']=r_arr*np.sin(np.deg2rad(theta_arr))\n",
    "    df2['sigmarr'] = (Pw - Pp)*(a**2/df['r']**2)+(sigmaHmax+sigmaHmin)/2*(1-(a**2/df['r']**2))+(sigmaHmax-sigmaHmin)/2*(1-(4*a**2/df['r']**2)+(3*a**4/df['r']**4))*np.cos(np.deg2rad((2*df['theta'])))\n",
    "    df2['sigmatt'] = -(Pw - Pp)*(a**2/df['r']**2)+(sigmaHmax+sigmaHmin)/2*(1+(a**2/df['r']**2))-(sigmaHmax-sigmaHmin)/2*(1+(3*a**4/df['r']**4))*np.cos(np.deg2rad(2*df['theta']))\n",
    "    df2['sigmatt_T'] = -(Pw - Pp)*(a**2/df['r']**2)+(sigmaHmax+sigmaHmin)/2*(1+(a**2/df['r']**2))-(sigmaHmax-sigmaHmin)/2*(1+(3*a**4/df['r']**4))*np.cos(np.deg2rad(2*df['theta']))+alpha*E/(1-PR)*(Tfluid-Trock)\n",
    "    df2['sigmart'] = (sigmaHmax-sigmaHmin)/2*(1+(2*a**2/df['r']**2)-(3*a**4/df['r']**4))*np.sin(np.deg2rad(2*df['theta']))\n",
    "    \n",
    "    # Defining the dataframe for only considering r = 0.1\n",
    "    mask=df2['r']==0.1\n",
    "    df3=df2[mask];y1 = df3['sigmarr'];y2 = df3['sigmatt'];y3 = df3['sigmart'];y4 = df3['sigmatt_T']\n",
    "    \n",
    "    # Creating UCS and Ts arrays for plotting.\n",
    "    UCSarr=[]; TSarr=[]\n",
    "    for i in range(0,len(np.linspace(0.1,1,359))):\n",
    "        UCSarr.append(UCS);TSarr.append(Ts)\n",
    "        \n",
    "    # Finding the failure types and zones\n",
    "    mask=df3['sigmatt_T']>UCS;mask2=df3['sigmatt_T']<Ts;shear_failure_zone = df3[mask];tensile_failure_zone = df3[mask2]\n",
    "    x_shear = a*np.cos(np.deg2rad(shear_failure_zone['theta']));y_shear = a*np.sin(np.deg2rad(shear_failure_zone['theta']))\n",
    "    x_tensile = a*np.cos(np.deg2rad(tensile_failure_zone['theta']));y_tensile = a*np.sin(np.deg2rad(tensile_failure_zone['theta']))\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    plt.title('Stress profile along $r = 0.1$ m',fontsize=14)\n",
    "    plt.plot(df3['theta'],UCSarr,'-.',label='UCS',color='darkred');plt.plot(df3['theta'],TSarr,'--',label='Ts',color='b')\n",
    "    plt.plot(df3['theta'],y1,label='$\\sigma_{r-r}$',color='purple');plt.plot(df3['theta'],y2,label='$\\sigma_{\\\\theta-\\\\theta}$',color='darkorange')\n",
    "    plt.plot(df3['theta'],y3,label='$\\sigma_{r-\\\\theta}$',color='green');plt.plot(df3['theta'],y4,label='$\\sigma_{\\\\theta-\\\\theta}$ with $\\Delta_T$',color='red')\n",
    "    plt.grid(which='both');plt.xlabel('$\\\\theta$',fontsize=12);plt.ylabel('Stress [MPa]',fontsize=12)\n",
    "    plt.legend(loc='upper right',fontsize=12);plt.xticks([0,45,90,135,180,225,270,315,360]);plt.xlim(0,360)\n",
    "    \n",
    "    plt.subplot(122)\n",
    "    plt.title('$\\sigma_{\\\\theta-\\\\theta}$ with $\\Delta_T$',fontsize=16)\n",
    "    sns.scatterplot(x='x',y='y',hue=df['sigmatt_T'],data=df,palette='icefire',s=30,label=None);normalized = plt.Normalize(df['sigmatt_T'].min(), df['sigmatt_T'].max())\n",
    "    plt.scatter(x_shear,y_shear,s=30,color='red',marker='x',label='Shear Failure');plt.scatter(x_tensile,y_tensile,s=30,color='blue',marker='x',label='Tensile Failure')\n",
    "    colormapmappable = plt.cm.ScalarMappable(cmap=\"icefire\", norm=normalized);colormapmappable.set_array([]);plt.colorbar(colormapmappable)\n",
    "    plt.xlim(-a*2,a*2);plt.ylim(-a*2,a*2);plt.gca().set_aspect('equal', adjustable='box');plt.gca().set_facecolor(\"lightgrey\");\n",
    "    handles, labels = plt.gca().get_legend_handles_labels();plt.legend([handles[-2],handles[-1]],[labels[-2],labels[-1]])\n",
    "\n",
    "    \n",
    "    plt.subplots_adjust(left=0.5, bottom=0.5, right=2.5, top=1.5, wspace=0.3, hspace=0.3)\n",
    "\n",
    "interactive_plot = widgets.interactive_output(Interactive_Thermal_Stress, {'sigmaHmax':sigmaHmax,'sigmaHmin':sigmaHmin,'Pw':Pw,'Trock':Trock,'Tfluid':Tfluid})\n",
    "interactive_plot.clear_output(wait = True) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6093ebe0",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2b3616cdc7044f9d9c57c9f185dcd007",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                                                                   …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "27520cf4a6014d998364647e4d83a5e3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui2, interactive_plot)    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e075299",
   "metadata": {},
   "source": [
    "<hr style=\"height:2px;border-width:0;color:black;background-color:black\">\n",
    "\n",
    "### References\n",
    "\n",
    "Espinoza, D. N. (2020). 1.4 WP4: Solution of Navier's Equation for the Stress Field. Advanced Geomechanics - Class Notes.<br>\n",
    "Espinoza, D. N. (2021). 6.2 Kirsch solution for stresses around a cylindrical cavity. Introduction to Energy Geomechanics. <br>\n",
    "Espinoza, D. N. (2022). Handwritten notes Advanced Geomechanics course.<br>\n",
    "Zoback, M. D. (2007). Chapter 6 - Compressive and tensile failures in vertical wells. In Reservoir Geomechanics (pp. 167–205). Cambridge University Press.<br>\n",
    "\n",
    "Some keywords of the interactive plot are implemented from the examples in GeostatsGuy workflows. Special Thanks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d3f8c66",
   "metadata": {},
   "source": [
    "<hr style=\"height:2px;border-width:0;color:black;background-color:black\">\n",
    "\n",
    "### About Çınar Turhan\n",
    "\n",
    "*Petroleum and Natural Gas Engineer - BSc, METU 2022*<br>\n",
    "*MSc Student at Hildebrand Department of Petroleum and Geosystems Engineering, University of Texas at Austin*<br>\n",
    "*Graduate Research Assistant at RAPID Research Consortium.* <br>\n",
    "\n",
    "**Research Interests:**\n",
    "* Well Integrity\n",
    "* Well Repurposing\n",
    "* Carbon Storage\n",
    "* Geothermal Energy\n",
    "* Hole Cleaning\n",
    "* Machine Learning\n",
    "\n",
    "Feel free to contact me about anything related to this workflow, or to discuss my research. <br>\n",
    "**[LinkedIn - Cinar Turhan](https://www.linkedin.com/in/cinarturhan/)**\n",
    "<hr style=\"height:2px;border-width:0;color:black;background-color:black\">\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}