wu.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import pandas as pd\n",
    "from scipy import interpolate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "athm_pr = 100000                           # [Pa]\n",
    "dm = 3.8E-10                               # [-]\n",
    "R = 8.314462                               # [J/(mol - K)]\n",
    "M = 1.6E-2                                 # [kg/mol]\n",
    "Na = 6.0221415E23                          # [1/mol]\n",
    "b = -1\n",
    "dpdl = 0.1                                 # [MPa/m]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "properties = pd.read_csv('z,cg,n.csv')\n",
    "p = properties['P'].values*1000000\n",
    "t = properties['T'].values\n",
    "f_z = interpolate.bisplrep(p, t, properties['Z'].values, s=2)\n",
    "f_rho = interpolate.bisplrep(p, t, properties['Density(kg/m3)'].values, s=2)\n",
    "f_vis = interpolate.bisplrep(p, t, properties['Viscosity(uPa-s)'].values/1000000, s=2)\n",
    "f_cg = interpolate.interp1d(p, properties['cg(1/Mpa)'].values/1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "var_ranges = pd.read_csv('variability_ranges.csv').set_index('Variable')\n",
    "\n",
    "p = (var_ranges['Min']['Pressure']+var_ranges['Max']['Pressure'])/2  # [Pa]\n",
    "por = (var_ranges['Min']['Porosity']+var_ranges['Max']['Porosity'])/2       # [-]\n",
    "ro = (var_ranges['Min']['Radius']+var_ranges['Max']['Radius'])/2            # [m]\n",
    "tor = (var_ranges['Min']['Tortuosity']+var_ranges['Max']['Tortuosity'])/2   # [-]\n",
    "op = (var_ranges['Min']['Overburden Pressure']+var_ranges['Max']['Overburden Pressure'])/2      # [Pa]\n",
    "q = (var_ranges['Min']['q']+var_ranges['Max']['q'])/2   # [-]\n",
    "t = (var_ranges['Min']['t']+var_ranges['Max']['t'])/2   # [-]\n",
    "k = (var_ranges['Min']['k']+var_ranges['Max']['k'])/2   # [-]\n",
    "fd = (var_ranges['Min']['Fractal Dimension']+var_ranges['Max']['Fractal Dimension'])/2   # [-]\n",
    "ih = (var_ranges['Min']['Isosteric heat']+var_ranges['Max']['Isosteric heat'])/2   # [J/mol]\n",
    "lp0 = (var_ranges['Min']['lp0']+var_ranges['Max']['lp0'])/2   # [Pa]\n",
    "a0 = (var_ranges['Min']['a0']+var_ranges['Max']['a0'])/2   # [-]\n",
    "a1 = (var_ranges['Min']['a1']+var_ranges['Max']['a1'])/2   # [-]\n",
    "beta = (var_ranges['Min']['beta']+var_ranges['Max']['beta'])/2   # [-]\n",
    "T = (var_ranges['Min']['T']+var_ranges['Max']['T'])/2   # [K]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wu(p,ro,por,tor,op,q,t,k,fd,ih,lp0,a0,a1,beta,T):\n",
    "    Z = interpolate.bisplev(p,T,f_z)   # [-]\n",
    "    vis = interpolate.bisplev(p,T,f_vis)   # [Pa s]\n",
    "    cg = f_cg(p)   # [1/MPa]\n",
    "    dzdp = (interpolate.bisplev(p+10000,T,f_z)-interpolate.bisplev(p-10000,T,f_z))/20000\n",
    "\n",
    "    porosity = por*((op-p)/athm_pr)**-q         # [-]\n",
    "    radius = ro*((op-p)/athm_pr)**-t            # [m]\n",
    "    kn = (vis/(p*2*radius))*(math.pi*Z*R*T/(2*M))**0.5   # [-]\n",
    "    wv = 1/(1+kn)   # [-]\n",
    "    wk = 1-wv   # [-]\n",
    "    lp = lp0 * math.exp(-ih/(R*T))   # [Pa]\n",
    "    theta = (p/Z)/(lp+p/Z)   # [-]\n",
    "    rad = theta*dm   # [m]\n",
    "    fc1 = (porosity/tor)*(1-rad/radius)**2   # [-]\n",
    "    fc2 = fc1 * (((1-rad/radius)**-2)-1)   # [-]\n",
    "    alpha = a0*2*math.atan(a1*kn**beta)/math.pi   # [-]\n",
    "    B1 = wv*fc1*radius**2*p*M*(1+alpha*kn)*(1+4*kn/(1+kn))/(8*vis*Z*R*T)   # [s]\n",
    "    B2 = wk*fc1*2*radius*(dm/radius)**(fd-2)*((8*Z*M/(math.pi*R*T))**0.5)*p*cg/(3*Z)   # [s]\n",
    "    if k>=1:\n",
    "        gam = 0\n",
    "    else:\n",
    "        gam = 1 \n",
    "    ds0 = 8.29E-7 * T**0.5 * math.exp(-ih**0.8/(R*T))    # [m2/s]\n",
    "    ds = ds0 * ((1-theta) + 0.5*k*theta*(2-theta) + gam*0.5*k*(1-k)*(theta**2))/((1-theta+0.5*k*theta)**2)    # [m2/s]\n",
    "    Csc = 4*theta*M/(math.pi*dm**3*Na)    # [kg/m3]\n",
    "    B3 = fc2*ds*Csc/(p)   # [s]\n",
    "    D1 = (B1*R/M)/(1/(T*Z)+p*dzdp/(T*Z**2))   # [m2/s]\n",
    "    D2 = (B2*R/M)/(1/(T*Z)+p*dzdp/(T*Z**2))   # [m2/s]\n",
    "    D3 = (B3*R/M)/(1/(T*Z)+p*dzdp/(T*Z**2))   # [m2/s]\n",
    "    DT = D1+D2+D3   # [m2/s]\n",
    "\n",
    "    J = (B1+B2+B3)*dpdl*86400*365/1000\n",
    "\n",
    "    return B1, B2, B3, J, DT, kn, dzdp, Z, cg, vis,  \n"
   ]
  }
 ],
 "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.10"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "6849aba498b485528751ffc0985cbf38aa555fdf46032f476b156a6d588b32a5"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}