BHYWI04_1D_AD_Transport.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1D Advektion-Dispersion\n",
    "Analytische Lösung\n",
    "\n",
    "Basierend auf einem XLS-Worksheet von Prof. Rudolf Liedl \\\n",
    "letzte Änderung 2021 12 08 durch Thomas.Reimann@tu-dresden.de"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aktueller Stand der der Programmierung:\n",
    "\n",
    "-\n",
    "- Text ergänzen (TR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import special\n",
    "import numpy as np\n",
    "from ipywidgets import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Daten und Parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "315e5c23240e4609a7e9adc34c828d25",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatText(value=15.0, description='LENGTH', step=1.0), IntText(value=1800, description='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.transport(l, t1, ci, c0, m, Q, n, a)>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#FUNCTIONS FOR COMPUTATION; ADS = ADVECTION, DISPERSION AND SORPTION - EVENTUALLY SET RETARDATION TO 1 FOR NO SORPTION\n",
    "\n",
    "def IC(PE,r_time):\n",
    "    IC1 = np.sqrt(0.25*PE/r_time)*(1-r_time)\n",
    "    if (IC1>0):\n",
    "        IC2 = 1-(1-special.erfc(abs(IC1)))\n",
    "    else:\n",
    "        IC2 = 1+(1-special.erfc(abs(IC1)))\n",
    "    IC3 = np.sqrt(0.25*PE/r_time)*(1+r_time)\n",
    "    if (IC3>0):\n",
    "        IC4 = 1-(1-special.erfc(abs(IC3)))\n",
    "    else:\n",
    "        IC4 = 1+(1-special.erfc(abs(IC3)))\n",
    "    if IC4 == 0:\n",
    "        IC5 = IC2\n",
    "    else:\n",
    "        IC5 = IC2+np.exp(PE)*IC4\n",
    "    IC  = 1-0.5*IC5\n",
    "    return IC\n",
    "\n",
    "def BC(PE,r_time, r_dur):\n",
    "    # BCx ist positiver puls\n",
    "    BC1 = np.sqrt(0.25*PE/r_time)*(1-r_time)\n",
    "    if (BC1>0):\n",
    "        BC2 = 1-(1-special.erfc(abs(BC1)))\n",
    "    else:\n",
    "        BC2 = 1+(1-special.erfc(abs(BC1)))\n",
    "    BC3 = np.sqrt(0.25*PE/r_time)*(1+r_time)\n",
    "    BC4 = special.erfc(BC3)\n",
    "    if BC4 == 0:\n",
    "        BC5 = BC2\n",
    "    else:\n",
    "        BC5 = BC2 + np.exp(PE) * BC4\n",
    "    \n",
    "    # BCCx ist negativer puls\n",
    "    if r_time > r_dur:\n",
    "        BCC1 = np.sqrt(0.25 * PE / (r_time - r_dur)) * (1 - (r_time - r_dur))\n",
    "        if BCC1 > 0:\n",
    "            BCC2 = 1 - (1 - special.erfc(abs(BCC1)))\n",
    "        else:\n",
    "            BCC2 = 1 + (1 - special.erfc(abs(BCC1)))\n",
    "            \n",
    "        BCC3 = np.sqrt(0.25 * PE / (r_time - r_dur)) * (1 + (r_time - r_dur))\n",
    "        BCC4 = special.erfc(BCC3)\n",
    "        if BCC4 == 0:\n",
    "            BCC5 = BCC2\n",
    "        else:\n",
    "            BCC5 = BCC2 + np.exp(PE) * BCC4\n",
    "            \n",
    "    if r_time <= r_dur:\n",
    "        BC = 0.5 * BC5\n",
    "    else:\n",
    "        BC = 0.5 * (BC5 - BCC5)\n",
    "    return BC\n",
    "    \n",
    "\n",
    "def transport(l,t1,ci,c0,m,Q,n,a):\n",
    "    # l : länge der säule, float\n",
    "    # t1 : versuchsende, float\n",
    "    # ci : initiale konzentration, float\n",
    "    # c0 : eingabekonzentration, float\n",
    "    # m : eingabemasse, float\n",
    "    # Q : durchfluss, float\n",
    "    # n : porosität, float\n",
    "    # a : dispersivität (a > 0), float\n",
    "    \n",
    "    # Data for plotting\n",
    "    t0 = 10       #Startzeit\n",
    "    dt = 10      #Zeitdiskretisierung\n",
    "    r  = 2    #Radius der Säule\n",
    "    \n",
    "    #Berechnung Zwischenergebnisse\n",
    "    A =     np.pi*r**2\n",
    "    q =     Q/A\n",
    "    v =     q/n\n",
    "    D =     a*v\n",
    "    PE =    l/a\n",
    "    dur =   m/Q/(c0-ci)\n",
    "    tPV =   l/v\n",
    "    r_dur = dur/tPV\n",
    "    r_dt =  dt/tPV\n",
    "\n",
    "    #Festlegung Zeitbereich\n",
    "    t = np.arange(t0, t1, dt)\n",
    "\n",
    "    #Berechnung Konzentration - Klammerterme\n",
    "    #Set fraction of distance\n",
    "    r_time = []\n",
    "    time   = []\n",
    "    conc   = []\n",
    "\n",
    "    #compute concentration  \n",
    "    for t in range(t0, t1, dt):      \n",
    "        # ADVECTION-DISPERSION\n",
    "        r_time = t/tPV\n",
    "        c = ci*IC(PE,r_time)+c0*BC(PE,r_time, r_dur)\n",
    "        \n",
    "        conc.append(c)\n",
    "        time.append(t)\n",
    "        \n",
    "    # measurements\n",
    "    t_obs = [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000]\n",
    "    c_obs = [1e-3, 5e-2, 8.5e-2, 9.7e-2, 9.9e-2, 9e-2, 5e-2, 1.5e-2, 2e-3, 5e-4]\n",
    "        \n",
    "# PLOT HERE        \n",
    "    fig = plt.figure(figsize=(9,6))\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(time,conc)\n",
    "    ax.scatter(t_obs, c_obs, marker=\"x\", c=\"red\", zorder=10)\n",
    "    plt.ylim(0, ci+c0)\n",
    "    plt.xlim(0,t1)\n",
    "    ax.set(xlabel='time', ylabel='conc.',\n",
    "       title='A-D transport')\n",
    "    ax.grid()\n",
    "    plt.legend(('1'), loc=4)\n",
    "    \n",
    "# COMMENTED THIS OUT - OPTIONAL SAVE IF INTENDED    \n",
    "   # fig.savefig(\"test.png\")\n",
    "    plt.show()\n",
    "    print(A, q, v, D, PE, dur, tPV, r_dur, r_dt)\n",
    "#THE FOLLOWING ALLOWS TO DEFINE THE RANGE OF SLIDERS DEPENDING ON OTHER SLIDER SETTINGS (HERE DISP AS FUNCTION OF LENGTH)\n",
    "    \n",
    "a_widget = FloatText(value=1.5,min=0.005, max=10,step=0.1,description='DISP 1', readout_format='.3f' )\n",
    "l_widget  = FloatText(value=15,min=1, max=100,step=1,description='LENGTH', readout_format='.2f' )\n",
    "\n",
    "def update_a_range(*args):\n",
    "    a_widget.min = 0.00015 * l_widget.value\n",
    "    a_widget.max = l_widget.value\n",
    "    a_widget.step = 0.01*l_widget.value\n",
    "    \n",
    "# l_widget.observe(update_D_range, 'value')\n",
    "\n",
    "interact(transport,\n",
    "         #DEFINE THE RANGE SUCH THAT NOT TOO MANY COMPUTATIONS ARE NECESSARY WHEN MOVING THE SLIDER!\n",
    "         l  = l_widget,\n",
    "         t1 = widgets.IntText(value=1800, min = 60, max = 86400, step = 60, description = 'MAX TIME'),\n",
    "         ci = widgets.FloatText(value=0,min=0,max=10,step=0.1, description='C initial',readout_format='.2f'),\n",
    "         c0 = widgets.FloatText(value=0.1,min=0, max=5,step=0.1,description='C0', readout_format='.2f' ),\n",
    "         m  = widgets.FloatText(value=10,min=0, max=1000,step=1,description='m', readout_format='.2f' ),\n",
    "         Q  = widgets.FloatText(value=0.2, min=0.00001, max=0.001, step=0.0001, description='Discharge Q',readout_format='.5f'),\n",
    "         n  = widgets.FloatText(value=0.2,min=0.02, max=0.6, step=0.01, description='PORO 1',readout_format='.2f'),       \n",
    "         a  = a_widget,\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}