{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "
\n", "\n", "# Exploratory Computing with Python\n", "*Developed by Mark Bakker*\n", "## Water Management Notebook 2: Manning's equation and emptying reservoir" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Exercise 1: Manning's equation\n", "The average velocity of water flowing in a canal may be estimated with Manning's equation. Manning's equation for the average velocity $v$ is\n", "\n", "$v=\\frac{1}{n}R_h^{2/3}S^{1/2}$\n", "\n", "where $n$ is Manning's coefficient, $R_h$ is the hydraulic radius, and $S$ is the slope of the bottom of the canal. \n", "A detailed description of Manning's equation can be founde here.\n", "\n", "\n", "We will consider canals with a trapezoidal cross-section as shown below. The slope of the banks is 1:1.\n", "\n", "\n", "\n", "Our goal in this exercise is to compute the water height $h$ for a given design discharge $Q_d$ given the width $w$, the roughness coefficient $n$, and the bottom slope $S$. This exercise consists of 6 steps." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 1*. Write a function that computes the hydraulic radius. The function should take two arguments: the height $h$ of the water above the bottom of the canal, and the width $w$ of the canal at the bottom (see Figure). The function should return the hydraulic radius. To test your function, the hydraulic radius for $h=2$ m, and $w=4$ m, is $R_h=1.2426$ m. (Note: in engineering practice it is, of course, ridiculous to give an answer with so many digits. But for testing code it is useful.)" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 2*. Write a function that computes the average velocity in the canal. The function should take four arguments: the height $h$, the width $w$, the roughness coefficient $n$, and the bottom slope $S$. Inside the function, call the function for the hydraulic radius that you wrote for step 1. To test your function, the average velocity for $h=2$ m, $w=4$ m, $n=0.025$, $S=0.005$ is $v=3.2692$ m/s. Remember that `2/3=0`, while `2.0/3=0.666...`. " ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 3*. Write a function that computes the total discharge in the canal. The function should take four arguments: the height $h$, the width $w$, the roughness coefficient $n$, and the bottom slope $S$. Inside the function, call the function for the average velocity. To test your function, the total discharge for $h=2$ m, $w=4$ m, $n=0.025$, $S=0.005$ is $Q=39.23$ m$^3$/s." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 4*. Make a graph of total discharge $Q$ vs. the height of the water $h$, where $h$ varies from 1 to 4 m. Draw three lines for three different width $w=4$ m, $w=6$ m, $w=8$ m. Otherwise, use $n=0.025$, and $S=0.005$. Add a legend, labels on the axes, etc. " ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 5*. When you know the design discharge $Q_d$, you can find the corresponding water height $h$ from the graph you just produced. This is not very convenient though. Especially since the graph is for specific values of the roughness coefficient $n$ and the slope $S$. We want to find the water height that corresponds to a given $Q$. We need one additional function, let's call it `Qdiff` that computes the difference between the discharge computed with the discharge function you wrote in Step 3 and the design discharge you want. Inside the function, call the function for the total discharge that you wrote in Step 3. The input arguments of `Qdiff` are the height $h$, the width $w$, the roughness coefficient $n$, the bottom slope $S$, and the design discharge $Q_d$. `Qdiff` should return the difference between the design discharge and the discharge you compute with the specified values of $h$, $w$, $n$, and $S$. Test your `Qdiff` function. Use `fsolve` (see previous Notebook) to find the water height $h$ corresponding to the design discharge and the parameters of the canal $w$, $n$, $S$. Show that your code produces the correct result. " ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 6*. Compute the water height for $Q_d$ = 20, 40, 60, 80, 100 m$^3$/s, for $w=6$ m, $n=0.02$, $S=0.004$ and print them to the screen. \n", "Make a graph for $h$ vs $Q_d$." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Polynomial fit\n", "The function `polyfit` may be used to fit a polynomial through a set of data. Given a set of $x$ and $y$ values, `polyfit` computes the parameters $p$ of a polynomial of degree $N$ that goes through the data\n", "\n", "$f(x) = p_0 x^N + p_1 x^{N-1} + ... + p_{N-1}x + p_N$ \n", "\n", "Once you have computed the parameters $p_n$, you can evaluate the polynomial using the `polyval` function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = array([1,3,4,6,8,9])\n", "y = array([1,4,2,6,5,7])\n", "p = polyfit(x,y,5)\n", "print p\n", "x2 = linspace(0,10,100)\n", "y2 = polyval(p,x2)\n", "plot(x,y,'bo')\n", "plot(x2,y2)\n", "ylim(-1,8)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "When the degree of the polynomial is equal to the number of data points minus one, the polynomial goes exactly through the data points as shown above where there are 6 data points and a polynomial of degree 5 (hence, 6 $p$ values) is fitted through the data. Higher order polynomials often wiggle through the data and may give unreasonable values outside the range of the data. A smoother polynomial is obtained when the degree of the polynomial is smaller than the number of data points minus 1. The best fit is obtained by computing the parameters such that the sum of squared differences between the fitted line and the data is minimized. For example, let's fit a polynomial of order 3 through the data of the previous example" ] }, { "cell_type": "code", "collapsed": false, "input": [ "x = array([1,3,4,6,8,9])\n", "y = array([1,4,2,6,5,7])\n", "p = polyfit(x,y,3)\n", "print 'the values of p are:',p\n", "x2 = linspace(0,10,100)\n", "y2 = polyval(p,x2)\n", "plot(x,y,'bo')\n", "plot(x2,y2)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "###Exercise 2. Emptying a reservoir \n", "Consider a reservoir filled with water that is emptied by opening a gate at the bottom of the reservoir. \n", "The water level in the reservoir is measured with a pressure transducer on the bottom of the reservoir while it is being emptied. \n", "The rate of change of the water level $h$ in the reservoir is theoretically related to the water level in the reservoir as\n", "\n", "$\\frac{\\text{d}h}{\\text{d}t} = \\alpha \\sqrt{h-h_0}$\n", "\n", "where $\\alpha$ and $h_0$ are constants. Rearrangement of the differential equation gives\n", "\n", "$\\frac{1}{\\sqrt{h-h_0}}\\text{d}h = \\alpha \\text{d}t$\n", "\n", "and integration gives\n", "\n", "$\\sqrt{h-h_0} = 2\\alpha t + \\beta$\n", "\n", "where $\\beta$ is a constant of integration. Rearrangment gives\n", "\n", "$ h = 4\\alpha^2 t^2 + 4\\beta\\alpha t + \\beta^2 + h_0$\n", "\n", "Hence, $h$ vs. $t$ is a parabola.\n", "\n", "*Step 1*. The raw data of the experiment is provided in the file `emptying_reservoir.csv`. A `csv` file is an ascii file that you can open in a standard editor (csv stands for comma separated values). Load the column labeled as 'Pressure' from the data file. Note the units of pressure. Pressure measurements include both air pressure and the pressure caused by the water column. Note: if you want to specify one colomn to read with the `usecols` keyword, put square brackets around the column number to make it a list. This may be a bit confusing, but `usecols` expects a sequence. `(3)` is not a sequence, it is just the number 3, while `[3]` is a sequence as it is the data type list. Alternatively, the somewhat funny syntax `(3,)` is a tuple and would work. If you want to read more than one column, for example columns 2 and 3, you can type either `(2,3)` (which is a tuple of two values) or `[2,3]` (which is a list of two values)." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 2* The pressure transducer starts recording pressures before the start of the experiment (when the reservoir is openened and starts to empty) and is removed after a bit more than an hour. The pressure is recorded for some time after the pressure transducer is removed from the reservoir. Plot the pressure and determine at what index the experiment starts and at what index the experiment ends. Zooming is, unfortunately, not yet supported for inline graphs in IPython notebooks (v.1.0). Use the `xlim` and `ylim` functions to zoom in on the beginning and end of the experiment to determine the index of the beginning and end of the experiment. " ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 3*. Create a new array that only contains the pressure data during the experiment. Make a graph for pressure against time. Set time equal to zero at the beginning of the experiment and determine the time between measurements from the data file." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 4*. Fit a parabola through the data using polyfit. Plot in one graph the data of the experiment and the fitted polynomial. Do you get a nice fit?" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Step 5*. Write a Python function for the change of water level with time $\\text{d}h/\\text{d}t$ based on the parabola of Step 4. Plot $\\text{d}h/\\text{d}t$ vs. $t$ based on the parabola. Next, compute the derivative numerically using the measured water levels during the experiment. Plot the numerical derivative based on the data on the same graph. Do you still get a good match between the model and the data?" ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Answers to the exercises" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Answers to Exercise 1" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 1\n", "def Rh(h,w):\n", " A = (w * h) + (h * h) # m^2, is the cross sectional area of the flow\n", " P = w + 2 * h * sqrt(2.) # m, is the wetted perimeter of the channel\n", " rv = A / P\n", " return rv\n", "Rh(2,4)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 2\n", "def Vavg(h,w,n,S):\n", " rv = Rh(h,w)**(2.0/3) * sqrt(S) / n\n", " return rv\n", "Vavg(2,4,0.025,0.005)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 3\n", "def Qtot(h,w,n,S):\n", " A = (w+h) * h\n", " return A * Vavg(h,w,n,S)\n", "Qtot(2,4,0.025,0.005)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 4\n", "h = linspace(1,4,100)\n", "Q = Qtot(h,4,0.025,0.005)\n", "plot(h,Q,label='w=4')\n", "Q = Qtot(h,6,0.025,0.005)\n", "plot(h,Q,label='w=6')\n", "Q = Qtot(h,8,0.025,0.005)\n", "plot(h,Q,label='w=8')\n", "legend(loc='best')\n", "xlabel('h (m)')\n", "ylabel('Q (m$^3$/s)')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 5\n", "def Qdiff(h,w,n,S,Qd):\n", " Qt = Qtot(h,w,n,S)\n", " return Qd - Qt\n", "from scipy.optimize import fsolve\n", "h50 = fsolve(Qdiff,3,args=(4,0.025,0.005,50))\n", "print 'design discharge is 50'\n", "print 'computed water height is ',h50\n", "print 'discharge computed with Qtot is ',Qtot(h50,4,0.025,0.005)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 6\n", "w = 6.0\n", "n = 0.02\n", "S = 0.004\n", "h = zeros(5)\n", "Qd = arange(20,101,20)\n", "for i in range(5):\n", " h[i] = fsolve(Qdiff,1,args=(w,n,S,Qd[i]))\n", "print 'Design discharge is ',Qd\n", "print 'Corresponding water height is ',h\n", "plot(Qd,h,'b-o')\n", "xlabel('Design discharge (m$^3$/d)')\n", "ylabel('Water height (m)')\n", "title('w=6 m, n=0.02, S=0.004')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 1\n", "\n", "Answers to Exercise 2" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 1\n", "head = loadtxt('emptying_reservoir.csv',skiprows=54,usecols=[1],delimiter=\";\")" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 2\n", "subplot(211)\n", "plot(head)\n", "xlim(170,190)\n", "ylim(1115,1120)\n", "subplot(212)\n", "plot(head)\n", "xlim(325,335)\n", "ylim(1040,1080)\n", "print 'Data from index 184 to 332'" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 3\n", "headnew = head[184:332]\n", "print len(headnew)\n", "time = arange(0,len(headnew)*30,30)\n", "plot(time,headnew)\n", "xlabel('time (seconds')\n", "ylabel('head (cm)')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 12 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 4\n", "a = polyfit(time,headnew,2)\n", "hfit = polyval(a,time)\n", "plot(time,headnew,'b',label='data')\n", "plot(time,hfit,'r',label='parabola')\n", "xlabel('time (seconds')\n", "ylabel('head (cm)')\n", "legend(loc='best')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "code", "collapsed": false, "input": [ "# Step 5\n", "def dhdt(t):\n", " return 2*a[0]*t + a[1]\n", "dhdtnum = (headnew[1:] - headnew[:-1]) / 30.0\n", "plot(0.5*(time[1:]+time[:-1]), dhdtnum,label='numerical 1') \n", "dhdtnum2 = (headnew[2:] - headnew[:-2]) / 60.0\n", "plot(time[1:-1],dhdtnum2,'r',label='numerical 2')\n", "dhdtpoly = dhdt(time)\n", "plot(time,dhdtpoly,'k',label='parabola')\n", "legend(loc='best')\n", "xlabel('time (seconds')\n", "ylabel('head (cm)')" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Back to Exercise 2" ] } ], "metadata": {} } ] }