examples/sand_transport.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pySand\n",
    "#### Sand transport examples - created with version 1.6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stokes model should be used for vertical or deviated sand transport (typically in wells) which estimates critical (mix) velocity for sand transport. The hydro model should be used for horisontal sand transport (typically for subsea flowlines) and estimates two critical velocities, one for transition from stationary beds to moving beds and one for the transition from moving beds to fully dispersed flow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pysand.transport import stokes, hydro\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deviated well using stokes model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function stokes in module pysand.transport:\n",
      "\n",
      "stokes(rho_m, mu_m, d_p, angle, rho_p=2650)\n",
      "    Sand lifting rate calculation model for vertical or deviated (up to ~80 deg) inclination\n",
      "    Based on Stokes law for settling velocities, but corrected for turbulent flow around sand grains\n",
      "    Mixture fluid properties and inclination at point of interest (e.g. down hole in a well)\n",
      "    :param rho_m: Mixture density [kg/m3]\n",
      "    :param mu_m: Mixture viscosity [kg/ms]\n",
      "    :param d_p: Particle diameter [mm]\n",
      "    :param angle: Inclination from vertical [deg]\n",
      "    :param rho_p: Particle density [kg/m3], default = 2650 (quartz sand)\n",
      "    :return: Particle settling velocity at point of interest [m/s]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(stokes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Input parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_m = 350  # Mix density [kg/m3]\n",
    "mu_m = 1e-4  # Mix viscosity [kg/ms]\n",
    "d_p = .4  # Particle diameter [mm]\n",
    "angle = 24  # Angle from vertical [deg]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Critical velocity for transport [m/s]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.24"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stokes(rho_m, mu_m, d_p, angle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Horisontal flowline using hydro model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function hydro in module pysand.transport:\n",
      "\n",
      "hydro(D, rho_l, mu_l, d_p, e=5e-05, rho_p=2650)\n",
      "    Equinor sand transport model for horizontal pipelines\n",
      "    Based on T. Søntvedt (1995) and R. Schulkes (2002) work in Hydro\n",
      "    :param D: Pipe diameter [m]\n",
      "    :param rho_l: Liquid density [kg/m3]\n",
      "    :param mu_l: Dynamic liquid viscosity [Pa.s]\n",
      "    :param d_p: Sand grain diameter [mm]\n",
      "    :param e: Roughness of pipe [m] Optional. Default value: 5e-5 m\n",
      "    :param rho_p: Particle density [kg/m3], default = 2650 (quartz sand)\n",
      "    :return: list of critical velocities for transitions from stationary to moving bed\n",
      "    and from moving bed to dispersed transport [m/s]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(hydro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Input parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "D, rho_l, mu_l, d_p = 0.12, 780, 1.2e-3, .7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Critical velocities [m/s] for moving dunes and dispersed flow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.78, 1.44]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hydro(D, rho_l, mu_l, d_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sensitivity study for varying inclination with Stokes model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inclination [deg] in study:  [0, 10, 20, 30, 40, 50, 60, 70, 80]\n"
     ]
    }
   ],
   "source": [
    "angles = list(range(0, 81, 10))\n",
    "print('Inclination [deg] in study: ', angles)\n",
    "velocities = []\n",
    "for angle in angles:\n",
    "    velocities.append(stokes(rho_m, mu_m, d_p, angle))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Critical velocity [m/s]')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU9b3/8deHEEKAQIBEdgSRVRHQSN231rrUtdUW7CLV1tveert7i7f9aZfbakv31vbW3WoV0FJFa6XWulTbSoIoyKbIkpCwBNnClvXz++OcwBAmyQQzOTOZ9/PxmEdmvmeZz5yZnM853+8536+5OyIikrm6RB2AiIhES4lARCTDKRGIiGQ4JQIRkQynRCAikuGUCEREMpwSQYYws/8xs7tbmP5xM/trO7yPm9mx73U9TdZ5jplteI/rGG5mu80sq51iOt3M3g7XeUV7rLMjmNk6M/tA1HFIalEiSFNmdo2ZlYQ7oo1m9hczO6O5+d39B+7+mXDZEeEOu2vM9D+4+wc7IvYouHupu/dy93oAM3vBzD7zHlb5XeDX4Tofb58oo2Vm95vZ/0YdR3sws2+b2UNRx5EulAjSkJl9Ffg58ANgADAc+A1weTPzd41XLu/J0cCyI1kwXb+PdIk7XeJMKe6uRxo9gD7AbuDqFub5NvAY8BCwC/hMWPZQOL0U8HA9u4FTgRnAyzHrOA54FtgGbAb+JyyfCvwL2AFsBH4NdItZzoFj48Q0DShpUvYVYH74PAf4cRjbZuD/gNxw2jnAhpjlxgMvhDEsAy6LmZYL/ARYD+wEXg7LRoSxdQW+D9QD+8PP/2vgDuAnTeJ7EvhynM/yDtAA7AuXzwEGA/PD7bUa+GxL30ecdV4MLAeqgHLg62F5X+ApoBLYHj4fGrPcC8D3gFfCZf8KFMRM/2S4Ld4FvgmsAz4Q5/1vAGqBmvAzPRmWrwO+ASwBqsPtNzPcBlVhzFfGrGdGuM1/HMa7FrioyfQ14bJrgY/HlL8C/Cr83lYC749Zri3b98bwc9SGn+WNqP9vU/0ReQB6tPELgwuBOqBrC/N8O/wnuILgrC+XQxPBCMKdYswyMwgTAZBHsJP/GtA9fP2+cNpJwCnhDmEEsIKYnSXNJ4Ie4T//6JiyYmBa+Pzn4T96v/D9ngRuC6edQ5gIgOxwR/A/QDfgvHC9Y8PpdxDsHIcAWcBpBDvqQz5zOM9nYmKZClQAXcLXBcBeYEAz23gdMTtU4EWCs7LuwGSCHff7m/s+4qxvI3Bm+LwvcGL4vD/wkXD75QGPAo/HLPcCwU55TPg9vwDcHk6bQLAjPCvcBj8l+O0clgjC+e8H/jfO53wdGMbBxHw1wY65C/AxYA8wKOZ3VAt8Ntz+nw+3qwE9CXbUjd/VIOC4mOXqCA4OssP17gT6Hcn2Jeb3rkfrD1UNpZ/+wFZ3r2tlvn+5++Pu3uDu+9r4HpcAm9z9J+6+392r3P1VAHdf5O7/dvc6d18H/A44u7UVuvte4AlgOoCZjQbGAfPNzAh2HF9x923uXkVQ7TUtzqpOAXoR7Oxq3P3vBEfJ082sC3Ad8CV3L3f3enf/p7tXJxDfQoIdz/vDomnAC+6+ubVlzWwYcAbwjXB7vQ7cTXA03qi176MWmGBmvd19u7u/Fsb1rrv/0d33htvl+xy+ve9z97fC9c4l2FECXAU85e4vhdvg/xGcybTVL929rDFud3/U3SvCzzIHeJsgkTZa7+53edAe8wDBDn9AOK0BON7Mct19o7vHVq9tAX7u7rXhelcBH2qn7SstUCJIP+8CBQnUg5a9h/cYRnCUeRgzG2NmT5nZJjPbRbDDLkhwvQ8TJgLgGoIj271AIcER7yIz22FmO4BnwvKmBgNl7h67Q1tPcAZQQHDEGDf2BDwAfCJ8/gngwQSXGww0JrCmMTVq7fv4CEH10Hoze9HMTgUwsx5m9jszWx9u75eA/CZXP22Keb6XIFE2xnXgfd19D8Hvp60Oid3MPmVmr8d8V8dz6G/gQDzh9wvQK3z/jwGfAzaa2Z/NbFzMcuXuHtsL5vrwM7TH9pUWKBGkn38R1G23dsliS93KttblbBkwqplpvyWovx3t7r0JqmislfU1+itBEptMkBAeDsu3EtS3H+fu+eGjj7v3irOOCmBYePTfaDhBvfpWgm3TXOyx4m2Dh4DLzWwSQTtEolcDVQD9zCwvTkwtvd/Bie7F7n45cFT4vnPDSV8DxhJUzfUmqOaBxLb5RoKkHixg1oPgjLLZMForN7OjgbsI6uH7u3s+8GaC8eDuC9z9fIKzhJXhuhoNCc8OGw0n2LZHsn3VrXIbKBGkGXffCdwC3GFmV4RHjNlmdpGZ/SjB1VQSnKIf08z0p4CBZvZlM8sxszwze184LY+gnnd3eDT3+TbEXkfQqDeLoC3g2bC8gWCH8DMzOwrAzIaY2QVxVvMqQZ30f4ef+xzgUmB2uJ57gZ+a2WAzyzKzU80sJ856Njf9/O6+gaDd4kHgj4lWMbh7GfBP4DYz625mJwDXA39IZHkz6xbex9HH3WsJtm99ODmPIEnuMLN+wK2JrDP0GHCJmZ1hZt0ILnlt6X/+sG0SR0+CnWxlGPunCc4IWmVmA8zsMjPrSdDwvJuDnxOCJPjF8Hu9miAZP32E23czMKLJAYM0QxspDbn7T4GvAt8i+IcsIzhCS+gINjxd/z7wSnh6f0qT6VXA+QQ72E0EdcDnhpO/TlCtU0Ww857TxvAfBj4APNqkneMbBI3A/w6rQP5GcCTcNPYa4DLgIoIzgN8An3L3lTHxLSXYoW8Dfkj83/kvgKvMbLuZ/TKm/AFgIolXCzWaTtAgXQH8CbjV3Z9tw/KfBNaFn/1zHKyi+jlB4+dW4N8EVWYJCevfv0CwzTcSXMXT0o159xC0U+wws7i/JXdfTnBV1r8IdrYTCa72SUQXgjOcCoLv5mzgP2OmvwqMJvis3weucvfGqqy2bt9Hw7/vmtlrCcaXsezQKjmRzGZmZxFUEY1o0g4hSWRmMwiu4mr2pkhJHp0RiITMLBv4EnC3koBkEiUCEcDMxhPcoDaIoDpGJGOoakhEJMPpjEBEJMOlXedMBQUFPmLEiKjDEBFJK4sWLdrq7vFu0ky/RDBixAhKSkqiDkNEJK2Y2frmpqlqSEQkwykRiIhkOCUCEZEMp0QgIpLhlAhERDJc2l01JCKSaR5fXM6sBauo2LGPwfm53HTBWK6YMqT1BROkRCAiksIeX1zOzfOWsq826LG7fMc+bp63FKDdkoGqhkREUtisBasOJIFG+2rrmbVgVbu9hxKBiEgKq9gRf3yk5sqPhBKBiEgKG5yf26byI6FEICKSwr5w7uFDcOdmZ3HTBYcN4HfElAhERFJYdV0wRlJhXg4GDMnP5bYPT9RVQyIimcDdmb2wjBOG9mH+jckbxVNnBCIiKeqNDTtZtbmKj508LKnvo0QgIpKi5hSXkpudxWWTBif1fZKaCMzsQjNbZWarzWxmnOlHm9lzZrbEzF4ws6HJjEdEJF3sqa5j/usVXHLCIPK6Zyf1vZKWCMwsC7gDuAiYAEw3swlNZvsx8Ht3PwH4LnBbsuIREUknTy2pYE9NPdOmJrdaCJJ7RjAVWO3ua9y9BpgNXN5kngnAc+Hz5+NMFxHJSLOLyzj2qF6cOLxv0t8rmYlgCFAW83pDWBbrDeAj4fMrgTwz65/EmEREUt6qTVUsLt3BtJOHYWZJf79kJoJ40XuT118HzjazxcDZQDlQd9iKzG4wsxIzK6msrGz/SEVEUsjs4lKys4wPn9gxzabJTAQbgNjKraFARewM7l7h7h929ynAN8OynU1X5O53unuRuxcVFhYmMWQRkWjtr63nT4vL+eBxA+nXs1uHvGcyE0ExMNrMRppZN2AaMD92BjMrMLPGGG4G7k1iPCIiKW/Bsk3s2FvL9JOHd9h7Ji0RuHsdcCOwAFgBzHX3ZWb2XTO7LJztHGCVmb0FDAC+n6x4RETSwZziMob2zeW0UR3XXJrULibc/Wng6SZlt8Q8fwx4LJkxiIiki/Xv7uGf77zL184fQ5cuyW8kbqQ7i0VEUsTckjK6GFxV1LH31ioRiIikgLr6Bh4t2cC5Y49iUJ/2G2sgEUoEIiIp4PlVlWypqk56B3PxKBGIiKSAOcWlFOblcO64ozr8vZUIREQitmnnfv6+cgtXnzSU7KyO3y0rEYiIROyxRWU0OHy0qOOrhUCJQEQkUg0NzpySMk49pj8jCnpGEoMSgYhIhP75zruUbdvXId1NN0eJQEQkQrOLS+mTm80Fxw2MLAYlAhGRiGzbU8Nfl23myilD6J6dFVkcSgQiIhH50+JyauobIq0WAiUCEZFIuDuzF5YyeVg+4wb2jjQWJQIRkQi8VrqDt7fsZloEdxI3pUQgIhKBOcWl9OiWxSWTBkcdihKBiEhHq9pfy5NvbOSySYPplZPU0QASokQgItLBnnxjI/tq6yPpYC4eJQIRkQ42p7iUsQPymDwsP+pQACUCEZEOtbxiF29s2Mm0qcMw67hRyFqiRCAi0oHmFJfSrWsXrpwyJOpQDlAiEBHpIPtr6/nT4nIuPG4g+T26RR3OAUoEIiId5Jk3N7Frf13kdxI3pUQgItJBHllYytH9e3DKyP5Rh3IIJQIRkQ6wpnI3r67dxkeLhtGlS2o0EjdSIhAR6QBzSzaQ1cW4+qShUYdyGCUCEZEkq61v4LFFGzhv3FEc1bt71OEcRolARCTJnluxha27q1Oig7l4kpoIzOxCM1tlZqvNbGac6cPN7HkzW2xmS8zs4mTGIyIShTnFpQzoncPZYwqjDiWupCUCM8sC7gAuAiYA081sQpPZvgXMdfcpwDTgN8mKR0QkChU79vHiW5V8tGgYXbNSsxImmVFNBVa7+xp3rwFmA5c3mceBxhEZ+gAVSYxHRKTDPVqygQaHjxalZrUQJDcRDAHKYl5vCMtifRv4hJltAJ4G/iveiszsBjMrMbOSysrKZMQqItLu6hucuSVlnHFsAcP69Yg6nGYlMxHEu1DWm7yeDtzv7kOBi4EHzeywmNz9TncvcveiwsLUrGMTEWnqldVbKd+xL+XuJG4qmYlgAxD76YdyeNXP9cBcAHf/F9AdKEhiTCIiHWZ2cSl9e2Rz/oQBUYfSomaHxjGz+Qksv83dZzQzrRgYbWYjgXKCxuBrmsxTCrwfuN/MxhMkAtX9iEja27q7mmeXb+ZTp44gp2tW1OG0qKUx0sYDn2lhuhFcFRSXu9eZ2Y3AAiALuNfdl5nZd4ESd58PfA24y8y+QlBtNMPdm1YfiYiknT+9Vk5tvafsvQOxWkoE33T3F1ta2My+09J0d3+aoBE4tuyWmOfLgdMTiFNEJG24O48Ul3LS0X0ZPSAv6nBa1WwbgbvPbVpmZl3MrHdL84iIZLqS9dtZU7knZcYkbk2rjcVm9rCZ9TaznsByYJWZ3ZT80ERE0tPshWX0yunKJScMijqUhCRy1dAEd98FXEFQzTMc+GRSoxIRSVM799Xy56UVXDZ5MD26tVT7njoSSQTZZpZNkAiecPdaDr8fQEREgPlvVLC/tiEtGokbJZIIfgesA3oCL5nZ0cCuZAYlIpKu5hSXMn5QbyYO6RN1KAlrNhGY2almZu7+S3cf4u4Xh5d2lgLndlyIIiLp4c3ynbxZvovpU4dhllqjkLWkpTOCa4FFZjbbzGaY2UAAD9R1THgiIuljdnEpOV27cPmkpt2qpbZmWzLc/XMAZjaOoCvp+82sD/A88AzwirvXd0iUIiIpbl9NPU8sruDiiYPo0yM76nDapNU2Andf6e4/c/cLgfOAl4GrgVeTHZyISLp4eulGqqrr0qqRuFFCnc6ZWV8zO4Gg24lNwH3uXpTUyERE0sjs4lJGFvRk6sh+UYfSZq1e5Gpm3wNmAGuAhrDYCc4OREQy3uotuylet52ZF41Lq0biRonc7fBRYFQ4ypiIiDQxt6SMrl2Mj5w4NOpQjkgiVUNvAvnJDkREJB3V1DXwx0Ub+MD4ARTm5UQdzhFJ5IzgNmCxmb0JVDcWuvtlSYtKRCRN/G3FZt7dU8PHUnwUspYkkggeAH4ILOVgG4GIiACzi8sY3Kc7Z41O32F0E0kEW939l0mPREQkzWzYvpd/vF3JF88bTVaX9GskbpRIIlhkZrcB8zm0aui1pEUlIpIG5pZsAODqovRsJG6USCKYEv49JaZMl4+KSEarb3AeLSnjzNGFDO3bI+pw3pNWE4G7q4M5EZEmXnq7ko0793PLJROiDuU9a6n30UtaWziReUREOqPZC0vp37Mb7x8/IOpQ3rOWzghmmVk50FILyA+Ap9o3JBGR1Lalaj/PrdjCdWeMpFvXhHrqSWktJYLNwE9bWf7tdoxFRCQtzHutnLoGT5vB6VvTUjfU53RgHCIiacHdmVNcxtQR/RhV2CvqcNpF+p/TiIh0oFfXbmPt1j2d5mwAlAhERNpkTnEZed27cvHEQVGH0m6UCEREErRzby1PL93IFZOHkNstK+pw2k2ricDMSszsC2bWt60rN7MLzWyVma02s5lxpv/MzF4PH2+Z2Y62voeISEd5/PVyqusaOlW1ECR2RjANGAwUhwPZX2AJjLxgZlnAHQTjHU8AppvZIXdeuPtX3H2yu08GfgXMa/MnEBHpAO7OIwtLOX5Ib44f0ifqcNpVImMWr3b3bwJjgIeBe4FSM/uOmbU0JttUYLW7rwkHtZkNXN7C/NOBRxIPXUSk4ywt38nKTVVMO3l41KG0u0THLD4B+AkwC/gjcBWwC/h7C4sNAcpiXm8Iy+Kt/2hgZHPrM7MbwiqqksrKykRCFhFpV48sLKN7dhcumzw46lDaXSJjFi8CdgD3ADPdvbEH0lfN7PSWFo1T5s3MOw14zN3r40109zuBOwGKioqaW4eISFLsqa5j/uvlfGjiYHp3z446nHaXSO+jV7v7mtgCMxvp7mvd/cMtLLcBiG1RGQpUNDPvNOALCcQiItLh/rx0I3tq6pmexqOQtSSRqqHHEixrqhgYbWYjzawbwc5+ftOZzGws0Bf4VwLrFBHpcHOKyxhV2JOTjm7zxZNpodkzAjMbBxwH9DGz2CP/3kD31lbs7nVmdiOwAMgC7nX3ZWb2XaDE3RuTwnRgtrurykdEUs5bm6tYtH4737x4PAlcMJmWWqoaGgtcAuQDl8aUVwGfTWTl7v408HSTsluavP52IusSEYnCnOIysrOMD58Y91qXTqGlTueeAJ4ws1PdXdU2IpJxquvqmffaBj44YSD9e+VEHU7StFQ19N/u/iPgGjOb3nS6u38xqZGJiETsr8s2s31vbae7k7iplqqGVoR/SzoiEBGRVDOnuIwh+bmccWxB1KEkVUtVQ0+Gfx/ouHBERFJD2ba9vLx6K189fwxdunTORuJGiXQ696yZ5ce87mtmC5IblohItOYUl9HF4KqThkYdStIlch9Bobsf6BXU3bcDRyUvJBGRaNXVN/DoojLOHlPI4PzcqMNJukQSQb2ZHehlKewXSNf8i0in9eJblWzeVc20qZ2vg7l4Euli4pvAy2b2Yvj6LOCG5IUkIhKt2cVlFPTK4bxxmVH50WoicPdnzOxE4JSw6CvuvjW5YYmIdLzHF5dz+19WsmnXfnrldOXPSzZyxZTOeyNZo0TOCABOIzgTaPRUEmIREYnM44vLuXneUvbVBp0g766u4+Z5SwE6fTJI5Kqh24EvAcvDx5fM7LZkByYi0pFmLVh1IAk02ldbz6wFqyKKqOMkckZwMTDZ3RsAzOwBYDFwczIDExHpSBU79rWpvDNJaIQygo7nGnWuwTpFRIC87vGPizPh8tFEzghuAxab2fMEo46dhc4GRKQTmVtcxq79dWSZUR/TI35udhY3XTA2wsg6RiJXDT1iZi8AJxMkgm+4+6ZkByYi0hGeXb6ZmfOWcOboAi6fNJif/e1tKnbsY3B+LjddMLbTNxRDy72PntikaEP4d7CZDXb315IXlohI8r265l1ufPg1Jg7N5/8+cRI9c7pyVVHn7mk0npbOCH7SwjQHzmvnWEREOszyil185vclDO2by30zTqZnTqJX03c+LfU+em5HBiIi0lFK393LtfctpFdOV35//fvo17Nb1CFFKpH7CHqY2bfM7M7w9WgzuyT5oYmItL8tVfv55L2vUlvfwO+vm8qQDLgqqDWJXD56H1BDcHcxBG0F/5u0iEREkmTX/lpm3FvMll3V3DvjZEYPyIs6pJSQSCIYFQ5ZWQvg7vsIrh4SEUkb+2vr+ewDJby1uYrffuJEThzeN+qQUkYirSM1ZpZL2PW0mY0CqpMalYhIO6qrb+CLjyzm1bXb+MW0yZwzNjN6FU1UIong28AzwDAz+wNwOjAjiTGJiLQbd+dbj7/JX5dv5tZLJ3D55M5/X0BbJXJD2V/NbBFBN9QGfEndUItIuvjxX1cxu7iMG889lk+fPjLqcFJSq4nAzOYDjwDz3X1P8kMSEWkf97y8ljuef4fpU4fxtQ+OiTqclJVIY/FPgDOB5Wb2qJldZWbdkxyXiMh78vjicr731HIuPG4g/3vFRMx0jUtzWk0E7v6iu/8ncAxwJ/BRYEsiKzezC81slZmtNrOZzczzUTNbbmbLzOzhtgQvIhLP86u28PVH3+CUY/rx82mTyeqiJNCShO6pDq8auhT4GHAi8EACy2QBdwDnE9x7UGxm8919ecw8owl6Mj3d3bebmZryReQ9WbR+O59/aBFjB+Zx16eK6J6dFXVIKS+RNoI5wPsIrhy6A3ihcZCaVkwFVrv7mnA9s4HLCUY5a/RZ4A533w7g7gmdaYiIxPPW5iquu7+YAb27c/+np5LXPTvqkNJCImcE9wHXuHt9q3MeaghQFvN6A0FCiTUGwMxeAbKAb7v7M01XZGY3ADcADB8+vI1hiEgmKN+xj0/ds5BuXbvw4HXvozAvJ+qQ0kYibQTPHEESgPh3H3uT112B0cA5wHTgbjPLP2wh9zvdvcjdiwoLC48gFBHpzLbtqeGT97zKnpo6fn/dVIb37xF1SGkl0aEqj8QGILZj76FARZx5nnD3WndfC6wiSAwiIgnZU13Hp+9bSPn2fdxz7cmMH9Q76pDSTjITQTEw2sxGmlk3YBowv8k8jwPnAphZAUFV0ZokxiQinUhNXQOfe2gRb1bs4o5rTmTqyH5Rh5SW2jJC2SFaG6HM3evM7EZgAUH9/73uvszMvguUuPv8cNoHzWw5UA/c5O7vtvVDiEjmaWhwvjr3df7x9lZmXXUCH5gwIOqQ0pa5N622DycEg9U3x909khHKioqKvKSkJIq3FpEU4e58e/4yHvjXem6+aBz/cfaoqENKeWa2yN2L4k3TCGUiknZ+9ffVPPCv9dxw1jFKAu0g0RvKjgcmAAe6lnD33ycrKBGR5jz07/X89Nm3+MiJQ5l54biow+kUErmh7FaCyzsnAE8DFwEvA0oEItKhnl66kf/3xJu8f9xR3P6RiXRR1xHtIpGrhq4C3g9scvdPA5MA3akhIh3qldVb+fLs1zlpeF9+fc2JZGcl86LHzJLIltwXdilRZ2a9CTqcOya5YYmIHLR0w05u+H0JIwt6cs+1J5PbTf0HtadE2ghKwrt97wIWAbuBhUmNSkQktKZyNzPuW0h+j278/vqp9Omh/oPaWyIjlP1n+PT/zOwZoLe7L0luWCIisHnXfj55T3Dc+eD1UxnQW0OhJEOrVUNmdqWZ9QFw93VAqZldkezARCSz7dxby6fuWciOvTXc/+mpHFPYK+qQOq1E2ghudfedjS/cfQdwa/JCEpFMt6+mnuseKGbt1j3c9akiJg7tE3VInVoibQTxkkVC9x+IiLRVbX0DX3j4NV4r3c4d15zIaccWRB1Sp5fIGUGJmf3UzEaZ2TFm9jOCRmMRkXbV0OB847El/H3lFr53+fFcPHFQ1CFlhEQSwX8BNcAc4FFgP/CFZAYlIpnH3fnB0yuYt7icr54/hk+ccnTUIWWMRK4a2gPEHXheRKS9/O6lNdz98lquPfVo/uu8Y6MOJ6O01A31z939y2b2JIePLIa7X5bUyEQkY8wtLuP2v6zk0kmDufXS4zBT1xEdqaUzggfDvz/uiEBEJDM9u3wzM+ct4czRBfzk6knqPygCLXVD3dggPNndfxE7zcy+BLyYzMBEpPN6fHE5sxasomLHPhwY1i+X//vESXTrqv6DopDIVr82TtmMdo5DRDLE44vLuXneUsrDJABQWVXNs8s3RxpXJmupjWA6cA0w0sxixxrOAzScpIgckVkLVrGvtv6Qsv21DcxasIorpgyJKKrM1lIbwT+BjUAB8JOY8ipAfQ2JSJvt2FtD+Y59cadVNFMuyddSG8F6YD1waseFIyKdkbvz56Ub+fb8Zc3OMzg/twMjkljNthGY2cvh3yoz2xXzqDKzXR0Xooiks4od+/jMAyXc+PBiBufn8vULxpCbfeh4ArnZWdx0wdiIIpSWzgjOCP/mdVw4ItJZ1Dc4D/5rHbMWrKLB4VsfGs+M00bQNasLQ/N7HLhqaHB+LjddMFbtAxFq8c5iM+sCLHH34zsoHhHpBFZtqmLmvCUsLt3BWWMK+f4VxzOsX48D06+YMkQ7/hTSYiJw9wYze8PMhrt7aUcFJSLpaX9tPb95fjW/ffEdeuV05Wcfm8QVk4foTuEUl0h30oOAZWa2ENjTWKguJkQk1sK125g5bwlrKvdw5ZQhfOtD4+nfKyfqsCQBiSSC7xzpys3sQuAXQBZwt7vf3mT6DGAWUB4W/drd7z7S9xORjrdzXy23/2UljywsZWjfXB64bipnjymMOixpg5ZuKDsWGODuLzYpP4uDO+5mmVkWcAdwPrABKDaz+e6+vMmsc9z9xjZHLiKRe+bNjdzyxDK27q7ms2eO5Cvnj6FHN41blW5a+sZ+DvxPnPK94bRLW1n3VGC1u68BMLPZwOVA00QgImlm86793PLEmyxYtpkJg3pz97VFnDA0P+qw5Ai1lAhGuPthdxC7e4mZjUhg3UOAspjXG4D3xZnvI+FZxlvAV9y9rOkMZnYDcAPA8OHDE3hrEUmGhgbn4YWl/PAvK6mpb+AbF47jM2eOJDtLncWls8ChKhgAABHYSURBVJYSQfcWpiVyC2C8ywSajmvwJPCIu1eb2eeAB4DzDlvI/U7gToCioqLDxkYQkeRbvWU3N89bQvG67Zw2qj8/uHIiIwp6Rh2WtIOWEkGxmX3W3e+KLTSz60lszOINwLCY10OBitgZ3D2287q7gB8msF4R6UA1dQ389oV3uOP51eR2y+JHV53A1ScN1SWhnUhLieDLwJ/M7OMc3PEXAd2AKxNYdzEw2sxGEjQuTyPozfQAMxvk7hvDl5cBK9oQu4gk2aL127l53hLe2rybS04YxK2XHkdhni4J7Wxa6mJiM3CamZ0LNN5Z/Gd3/3siK3b3OjO7EVhAcPnove6+zMy+C5S4+3zgi2Z2GVAHbEPjHIikhKr9tcxasIoH/72eQb27c++MIs4bNyDqsCRJzD29qtyLioq8pKQk6jBEOq2/Ld/M/3viTTbt2s+1p47g6xeMpVeOLglNd2a2yN2L4k3TtysiAGyp2s93nlzOn5dsZOyAPO74+ImcOLxv1GFJB1AiEMlw7s7ckjK+/+cV7K9t4Gvnj+E/zh6l8YMziBKBSAZbu3UPN89bwr/XbGPqyH7c9uGJjCrsFXVY0sGUCEQyUG19A3e+tIZfPPc2OV278IMrJzLt5GF06aJLQjOREoFIhnmjbAff+OMSVm6q4sLjBvKdy49jQO+W7h+Vzk6JQKQTe3xx+YGRwAb26c6Yo3rxj9VbKczL4XefPIkLjhsYdYiSApQIRDqpxxeXc/O8peyrrQdg4879bNy5n9NH9eO3nyyid/fsiCOUVKHLAkQ6qVkLVh1IArHWvbtPSUAOoUQg0gktLt1O+Y59cadVNFMumUtVQyKdREOD8/eVW7jzpTUsXLcN4/DufgEG5yfSebBkEiUCkTS3v7aexxeXc9c/1vBO5R6G5OdyyyUT6NEti+88ufyQ6qHc7CxuumBshNFKKlIiEElTO/bW8NC/13P/P9ezdXc1xw3uzS+mTeZDEwfRNRwopnt21oGrhgbn53LTBWO5YsqQiCOXVKNEIJJmyrbt5Z6X1zK3pIy9NfWcPaaQ/zjrGE4d1f+wMQKumDJEO35plRKBSJpYumEnv3vpHZ5eupGsLsZlk4bw2bNGMm5g76hDkzSnRCCSwhoanBffquR3L73Dv9dsIy+nK5896xg+fdpIBvbR3cDSPpQIRFJQdV09T7xewV0vreHtLbsZ1Kc737x4PNOmDiNP9wBIO1MiEEkhO/fV8vCrpdz3ylq2VFUzbmAeP/vYJC45YTDZWbrtR5JDiUAkBZTv2Me9L69l9sJS9tTUc+boAn589STOHF2gQeIl6ZQIRCK0rGInd760hqeWbMSASycN5jNnjuS4wX2iDk0yiBKBSAdzd/7x9lbufGkNL6/eSs9uWXz6tBFcd8ZI3fUrkVAiEOkgNXUNPLWkgjtfWsPKTVUM6J3DzIvGMX3qcPrkqgFYoqNEIJJkVftreWRhKfe+vI5Nu/YzdkAeP756EpdNGqxxgSUlKBGIJMnGnfu475V1PPJqKVXVdZw2qj+3fWQi54wpVAOwpBQlApH3KHYUsMH5uVzzvuG8U7mb+a9X4MDFEwdxw5nHMHGoGoAlNSkRiLwHTUcBK9+xj1kLVpGdZXzy1KO57vSRDOvXI+IoRVqmRCDSBu7OlqpqVmzcxYqNVfzyubfjjgJW0CuHWy89LoIIRdouqYnAzC4EfgFkAXe7++3NzHcV8ChwsruXJDMmkURV19WzestuVmysYsXGXazcFOz8t+2paXXZTTv3d0CEIu0jaYnAzLKAO4DzgQ1AsZnNd/flTebLA74IvJqsWERa4u5U7q4+uMMPj/bfqdxNXUMwxldO1y6MHZjH+eMHMG5QHuMH9Wb8wN5c/Mt/xB0SUvcDSDpJ5hnBVGC1u68BMLPZwOXA8ibzfQ/4EfD1JMYiAgTX8gdH+QeP8Fds3MW7MUf5g/p0Z/yg3rx//FHBDn9Qb0b073FgsJdYN10w9pA2AtAoYJJ+kpkIhgBlMa83AO+LncHMpgDD3P0pM2s2EZjZDcANAMOHD09CqNIZVYZ1+bE7/NVbDh7ld+vahbED8jhv3MEd/riBefTt2S3h92gc9EWjgEk6S2YiiHeh9IGxtM2sC/AzYEZrK3L3O4E7AYqKiuKNxy0ZoOllmo073Jq6Bt6p3H3IDn/Fxiq27q4+sOzA3t0ZPyiPcxt3+gPzGFnQM+5RfltpFDBJd8lMBBuAYTGvhwIVMa/zgOOBF8KbawYC883sMjUYS1PxLtP82tw3+NEzK6ncXU1t/cGj/DEDenHO2MLwKD+PcQN7068NR/kimSaZiaAYGG1mI4FyYBpwTeNEd98JFDS+NrMXgK8rCQhAbX0Dayr3hEf3u7j/n+uorms4ZJ56d97dU8P1ZxzD+LAB95h2OsoXySRJSwTuXmdmNwILCC4fvdfdl5nZd4ESd5+frPeW9LJtT82BHX5sXX5NfbDj75bV5cDzpmrqGph50biODFek00nqfQTu/jTwdJOyW5qZ95xkxiLRq6tvYM3WPYfs8Fds3MWWqoN1+YV5OYwf1JszRxccaMA9prAn58x6QZdpiiSJ7iyWpNjeeJS/6eAO/+0tu6kJq3eys4xjj8rjjGMP7vDHDcqjoFdO3PXpMk2R5FEikMM0d3VOPHX1DazduueQHf7KjVVs2nXwztqCXjmMH5THjNNGHGi8HVXYq01dMOsyTZHkUSKQQ8S7OufmeUsBOGds4SFVOis27eKtzQeP8rt2MY49qhenjup/YIc/flBvCvPiH+W3lS7TFEkOJQI5xA+fWXlYJ2r7auv56tzXaYi5g6N/z26MH9SbT51y9IGqnWOPattRvoikBiWCDLC3po6tVTVU7t5PZVUNW3dXU1lVzdbd1THPa6isqo7bkyZAg8PMi8YduBmrMC9Hg6uIdBIZkQjaUuedLnHtq6ln6+5qthy2Qz+4Y298vrcm/s69X89uFPTqRkGvHKYMz6egVw6PlpSxa3/dYfMOyc/lc2ePek+fV0RSU6dPBC3VeUeZDOLFNXPeEt7dU82U4X3ZWlVN5e5qtjZzBL+nmZ173x7ZFPTKoTAvh0lD8w88L+jVjYK8HArD1/16diM7zo1XE4f00dU5IhnG3NOr656ioiIvKUn85uPTb/973OvPu3YxRhb0bM/Q2mTt1j0HOj9rTX7jzr1XDgXhTj3YuQc79cJewfP+veLv3NsqVc+gROTImdkidy+KN63TnxFUxEkCAHUNzugBvTo4moPe3rK72Wn3zTiZgl45FOR1o3/PnA5vgNXVOSKZpdMngsH5uXHPCIbk5/Kbj58UQUSB5s5UhuTncu64oyKISEQyVae/1u+mC8aSm511SFkq1Hmnalwiknk6/RlBqt6RmqpxiUjm6fSNxSIi0nJjcaevGhIRkZYpEYiIZDglAhGRDKdEICKS4ZQIREQyXNpdNWRmlcD6I1y8ANjajuG0F8XVNoqr7VI1NsXVNu8lrqPdvTDehLRLBO+FmZU0d/lUlBRX2yiutkvV2BRX2yQrLlUNiYhkOCUCEZEMl2mJ4M6oA2iG4mobxdV2qRqb4mqbpMSVUW0EIiJyuEw7IxARkSaUCEREMlzGJAIzu9DMVpnZajObGWEc95rZFjN7M6asn5k9a2Zvh3/7RhDXMDN73sxWmNkyM/tSKsRmZt3NbKGZvRHG9Z2wfKSZvRrGNcfMunVkXDHxZZnZYjN7KlXiMrN1ZrbUzF43s5KwLBV+Y/lm9piZrQx/Z6dGHZeZjQ23U+Njl5l9Oeq4wti+Ev7m3zSzR8L/haT8vjIiEZhZFnAHcBEwAZhuZhMiCud+4MImZTOB59x9NPBc+Lqj1QFfc/fxwCnAF8JtFHVs1cB57j4JmAxcaGanAD8EfhbGtR24voPjavQlYEXM61SJ61x3nxxzzXnU3yPAL4Bn3H0cMIlgu0Ual7uvCrfTZOAkYC/wp6jjMrMhwBeBInc/HsgCppGs35e7d/oHcCqwIOb1zcDNEcYzAngz5vUqYFD4fBCwKgW22RPA+akUG9ADeA14H8HdlV3jfb8dGM9Qgp3EecBTgKVIXOuAgiZlkX6PQG9gLeEFKqkSV5NYPgi8kgpxAUOAMqAfwQBiTwEXJOv3lRFnBBzcqI02hGWpYoC7bwQI/0Y6aLGZjQCmAK+SArGF1S+vA1uAZ4F3gB3uXhfOEtX3+XPgv4GG8HX/FInLgb+a2SIzuyEsi/p7PAaoBO4Lq9LuNrOeKRBXrGnAI+HzSONy93Lgx0ApsBHYCSwiSb+vTEkEFqdM183GYWa9gD8CX3b3XVHHA+Du9R6cug8FpgLj483WkTGZ2SXAFndfFFscZ9Yofmenu/uJBFWhXzCzsyKIoamuwInAb919CrCHaKqn4grr2i8DHo06FoCwTeJyYCQwGOhJ8H021S6/r0xJBBuAYTGvhwIVEcUSz2YzGwQQ/t0SRRBmlk2QBP7g7vNSKTYAd98BvEDQhpFvZo1jbkfxfZ4OXGZm64DZBNVDP0+BuHD3ivDvFoL67qlE/z1uADa4+6vh68cIEkPUcTW6CHjN3TeHr6OO6wPAWnevdPdaYB5wGkn6fWVKIigGRoct7t0ITgHnRxxTrPnAteHzawnq5zuUmRlwD7DC3X+aKrGZWaGZ5YfPcwn+QVYAzwNXRRWXu9/s7kPdfQTB7+nv7v7xqOMys55mltf4nKDe+00i/h7dfRNQZmZjw6L3A8ujjivGdA5WC0H0cZUCp5hZj/B/s3F7Jef3FVXDTEc/gIuBtwjql78ZYRyPENT51RIcJV1PULf8HPB2+LdfBHGdQXCauQR4PXxcHHVswAnA4jCuN4FbwvJjgIXAaoLT+ZwIv9NzgKdSIa7w/d8IH8saf+tRf49hDJOBkvC7fBzomyJx9QDeBfrElKVCXN8BVoa/+weBnGT9vtTFhIhIhsuUqiEREWmGEoGISIZTIhARyXBKBCIiGU6JQEQkwykRSNows91HuNw5Mb2DXmZH2Pts2Hvmf8a8Hmxmjx3JuuKsu7HH0MMGJo+N/wjX/byZ7Y63bhFQIpAM4+7z3f32I1w8HziQCNy9wt2vamH+tjrX3UvacX0AuPu5BNfvi8SlRCBpJzxCfiGmb/s/hHdfYmYnm9k/LRi/YGHjXbYxy84ws1+Hz+83s1+G868xs6vC8l5m9pyZvRYepV8eLn47MCrst36WmY2wcFyJsK/4+8L5F5vZuTHvN8/Mngn7kP9Rgp/xwvCzvQx8OKa8pwVjWhSH73N5WN7DzOaa2ZKwn/pXdQYgiera+iwiKWkKcBxBXyuvAKeb2UJgDvAxdy82s97AvlbWM4jgrupxBN0KPAbsB650911mVgD828zmE3SSdrwHHeA19tLa6AsA7j7RzMYR9P45Jpw2OYy3GlhlZr9y99jecA9hZt2Buwj6L1odfqZG3yTozuK6sOuNhWb2N+DzwHZ3P8HMjie4M1wkITojkHS10N03uHsDwU5vBDAW2OjuxQDuvssPdtnbnMfdvcHdlwMDwjIDfmBmS4C/EXT1O6C5FYTOIOgGAHdfCawHGhPBc+6+0933E/QXc3Qr6xpH0OHY2x7c+v9QzLQPAjPDbrlfALoDw8P3nx2+/5sE3TiIJERnBJKuqmOe1xP8lo22d8sbu57GbqQ/DhQCJ7l7bdjDaPdW1hOvC+p479EYa2ua+xwGfMTdVx1SGFaNiRwJnRFIZ7ISGGxmJwOYWV5Ml71t0YdgrIHasK6/8Qi+CshrZpmXCBIIYZXQcIJRro7ESmCkmY0KX0+PmbYA+K+YNpEpYfnLwEfDsgnAxCN8b8lASgTSabh7DfAx4Fdm9gbBaGatHcnH8wegyIKB3z9OsGPG3d8FXrFgMPFZTZb5DZBlZksJ6vRnuHs1RyCsQroB+HPYWLw+ZvL3gGxgSdhQ/b2Y9y8Mq7O+QVA1tPNI3l8yj3ofFUkBYfVTkbtvPcLls4Bsd98fnkk8B4wJkyNm9gLw9WRcnirpT20EIqmhEnjOzK4/wp11D+B5C0aZM+DzMUngeYJ+7GvbLVrpVHRGICKS4dRGICKS4ZQIREQynBKBiEiGUyIQEclwSgQiIhnu/wMoi7AuUieG4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(angles, velocities, 'o-')\n",
    "plt.title('Critical velocity for sand transport')\n",
    "plt.xlabel('Inclination [deg]')\n",
    "plt.ylabel('Critical velocity [m/s]')"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}