MATLAB MANUAL MSc course Fluvial systems (GEO4-4436): Marion Tissier Maarten Kleinhans
MSc course Fluvial systems (GEO4-4436): MATLAB MANUAL Marion Tissier Maarten Kleinhans August 30, 2014 Contents 1 Scripts, matrices and syntax 1.1 Getting started . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Creating you own work-directory . . . . . . . . . . . . . . . . . . . . 1.1.3 Introductory exercises . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Script files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Manipulating arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Operations on arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Matlab built-in functions for arrays . . . . . . . . . . . . . . . . . . . 1.4 Visualising data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Final assignment : Analysis of the Rhine flow discharge measured at Lobith 2 Tests, loops and functions 2.1 Relational and logical operators . . . . . . . . . . . . 5 5 5 6 6 10 13 13 14 17 19 23 29 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 33 1 . . . . . . . . . . . – . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . part I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 . . . . . . 37 37 39 44 47 50 3 Numerical modelling 3.1 A very simple morphodynamic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Flow computation: backwater equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Morphodynamic model with flow computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 55 58 60 4 Creative Matlab exercise 65 2.3 2.4 Flow control . . . . . . . . . . . . . . . 2.2.1 if-statements . . . . . . . . . . 2.2.2 for-loops . . . . . . . . . . . . . 2.2.3 while-loops . . . . . . . . . . . User-defined functions . . . . . . . . . . Final assignment : Analysis of the Rhine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . flow discharge measured at Lobith 2 . . . . . – . . . . . . . . . . . . . . . . . . . . part II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Begin at the beginning • Each chapter of this manual corresponds to one of the 4-hour MATLAB practicals. Below you can find a short description of their content. First practical How to use the main features of Matlab (manipulation of arrays of data and use of script files). First data analysis and visualisation. Second practical How to write more complex programs using tests, loops. How to create your own functions and use them in a script file to analyse data. Third practical Introduction to the modelling of physical processes. During this practical you will build a simple morphodynamic model. Fourth practical Creative assignment, where you should be able to answer a research question using Matlab. • During the practicals, you will have to work in couples. Change partner every practical. • The Matlab assignments should be handed in before Friday 4:00 pm in the week that the computer practical took place. The different files (Matlab code and answers to the questions) should be combined in a zip-file and sent to [email protected]. The answers to the questions should be written as comments in the script files or (if specified) in a separate word document. 3 The evaluation of the codes will be based on three main criteria : - does the program run? - does it do the right thing? - can we read/understand it easily? (= is it commented enough?) • You will often get error messages while running your Matlab programs. Do not panic and carefully read the displayed messages. Solving them without help is part of the learning process! The devil is in the details: most errors are in syntax. To understand what is going wrong with the larger code and data matrices, it may help to play with very small and simple matrices first. 4 Chapter 1 Scripts, matrices and syntax 1.1 Getting started To start MATLAB, navigate into the program list and click on ”Matlab.exe”. 1.1.1 Overview Your Matlab window should roughly look as in Figure 1.1. If it does not, click on ”View” in the menu bar and select ”Desktop layout”, ”Default”. Four main divisions appear in the window : - ”Current folder” : displays the content of the current work-directory, - ”Workspace” : lists the variables assigned by the user, - ”Command history” : shows the last commands used, 5 - ”Command window” : where the commands are typed and the outputs displayed. In the following, we will mainly work in the ”Command window”, highlighted in Figure 1.1. 1.1.2 Creating you own work-directory The current directory is displayed on the top part of the Matlab window (see Fig. 1.1). We advise you to start every new practical with creating a new work-directory. To do that, click the ”browse for folder” button (see Fig. 1.1) right to the ”Current directory” drop-down list on the Matlab button-bar. You can then browse through the directories and create a new folder at the desired location using the ”New directory” button. For instance, create a first folder called Matlab_practicals. In this folder create a folder called Practical1 in which your will be working on today. You can also use this button to move to an already existing directory. 1.1.3 Introductory exercises We will perform in this section some introductory exercises where Matlab is basically used as an ”advanced” calculator. • Type in the command window a = 32.8 and then press ”enter”. The following lines appear: >> a = 32.8 a = 32.8000 With this command, we defined a new variable a and assigned the value 32.8 to it. • Now let’s type: >> b = 3*a^2 and then press ”enter”. We defined a new variable, b, equal to 3 times the value of a squared. Its value is displayed in the command window. The basic operations and their symbols are presented in Table 1.1. 6 Figure 1.1: Organisation of the Matlab window 7 Operation Addition Substraction Multiplication Division Exponentiation Symbol + − ∗ / ˆ(e.g. 2ˆ3 = 8) Table 1.1: Basic operations and associated symbols • Type now: >> whos and then press ”enter”. The following lines are displayed in the Matlab command window: Name Size Bytes Class a 1x1 8 double b 1x1 8 double This command displays the different variables which have been previously defined (a and b), and gives their main characteristics (see presentation of arrays, §1.3, for more explanations). These variables are now listed in the ”Workspace” sub-window, and can be visualised by double-clicking on them. • Use the command clear, and then the command whos again. What is displayed? Any interpretation? NB: To get more information about the function clear, type help clear and press ”enter” in the command window. • Type now: >> c=a*30/200 and then press ”enter”. We are defining a new variable c, as function of a, which has been deleted previously (clear). The command cannot be executed since a is not defined anymore. An error message, in red, appears in the command window. 8 Function sin(), cos(), tan() asin(), acos(), atan() abs() sqrt() round() ceil() (or floor()) log(), log10() exp() Description sine, cosine and tangent (arguments in radians) inverse sine, cosine and tangent absolute value root square (equivalent to ˆ(0.5)) round towards nearest integer round towards plus (or minus) infinity natural logarithm, logarithm base 10 exponential Table 1.2: List of some mathematical functions • Type: >> a =-4, c=a*30/200 and then press ”enter”. Two commands are written in the same line, separated by a comma. They are executed successively by Matlab. a is first defined and then used in the calculation of c. • Type now the same command, but replace the comma by a semicolon. What is the difference? And if you add a semicolon at the end of the line? • Type now: >> c_bis =abs(c) In this command, we applied the Matlab function abs() to the variable c. The result of this operation is a second variable, c_bis, which is the absolute value of c. c is called input argument of the function abs. c_bis is the output argument. Numerous elementary mathematical functions are already implemented in Matlab. The table 1.2 presents a selection of the most commonly-used elementary functions. Type help elfun in the command window for an extensive list of these functions. 9 Remarks - To obtain help in Matlab: help function name displays information about the function directly in the command window; doc function name opens a separate window where a more detailed help is displayed. The explanations usually include a few examples, which are worth studying carefully; lookfor keyword looks for all the Matlab functions related to the keyword; - The commands you’ve recently entered can be retrieved by pressing the up-arrow key; They are also displayed in the ”Command history” window. - Variable names must start by a letter, and can be followed by character chain, including letters, numbers and any symbols except those used as mathematical operators by Matlab (+,-,*,...) or punctuation signs. Variable names can contain up to 31 characters; When no variable name is specified, a temporary variable name, ans, is generated by default. - Try to use the name of a function to name one of your variables. 1.2 Script files A script, also called program, is a succession of Matlab commands written in a separate text file. Using script files allows you to organize and save series of commands. These scripts can be edited/extended/corrected anytime you want, and executed several times. Writing the commands into a script is necessary when you perform advanced data analysis, since it allows you to check/improve your methodology, but also to be sure to perform the exact same data processing for several datasets for instance. The scripts are user-created text files which can be written using Matlab editor. Type: >> edit 10 The Matlab editor window opens. In this window, you can write series of commands and save them. The resulting script is saved as *.m file. The ”new script” button (see Fig. 1.1) can also be used to open the editor. It is important to add comment lines to your scripts to make them easier to understand. This is particularly useful when scripts are re-used long after they have been created. In a script file, anything written to the right of a percentage symbol is considered as a comment and will be ignored during the execution. In the Matlab editor, the comments appear in a different color (green by default). It is recommended to always give a similar structure to the script files. An example of structure for a ”classic” script file is given below: % Header Describe the script in a few words (objective, date of creation, author) % Initialization Clear your workspace using the command clear Close all the figure windows (close all) Define the main variables, load data files used in the script % Calculations Execute the calculations % Output and visualisation Commands to display the results, plot graphics, save variables A simple example of script can be find below: 11 % Example of a very simple Matlab script % M. Tissier, 8-8-2012 % Initialization clear displacement = 12; time = 10; % displacement (m) % duration (s) % Calculation v = displacement/time; % migration speed (m/s) % Output disp([’The migration speed is ’ num2str(v) ’m/s’]) Once the script has been saved, you can run it by typing its name (without the ”.m” extension) in the ”Command window” and then press ”enter”. When a script file is executed, Matlab execute the commands in the order they are written. EXERCISE 1 The aim is to calculate flow velocity and water depth in a river from known variables and parameters. (Find out through a web search what the difference is between ”variable” and ”parameter”.) Use the law of Chezy (see Kleinhans (2005) paper for equations). The flow resistance parameter for this river is C = 44 m0.5 /s, the bankful discharge Qbf = 2500 m3 /s, the channel width W = 500 m and the channel gradient S = 1.6 × 10−4 m/m. The commands will be written/organised in a small script file called Exercise1.m, saved in your current directory. Do not forget to comment the file, and give meaningful names to the variables. HINT: 1.6 × 10−4 can be written in Matlab as 1.6e-4. 12 1.3 1.3.1 Arrays Definition The array is the fundamental form that Matlab uses to store and manipulate data. This way of storing data allows for efficient calculations and analysis over entire datasets. An array is basically a list of elements (numbers or characters, see below) organised in rows and/or columns. Different types of arrays can be defined: - The simplest array, or scalar, contains only one value (1 row, 1 column). - The vector is a one-dimensional array which is basically a row or column of values. - Two-dimensional arrays are also called matrices, and are characterized by their amount of rows and columns. - Multi-dimensional arrays (with 3 or more dimensions) can also be defined. These arrays can contain different type of data, for instance: - Integers : whole numbers (-10 -4 1 0 3 12) - Floating numbers : fractional numbers (3.271 -2345.17...) - Characters : text ASCII character (a b z / ] ?) In the two first cases, they are called double arrays, and in the third, character array. In this section, we will learn how to create arrays in Matlab, extract the relevant information from one of them, and perform calculations with these arrays. 13 1.3.2 Manipulating arrays Creating vectors and matrices • Type: >> M = [1 2 5 10; 3 4 -2 7] and then press ”enter”. A 2-by-4 matrix, called M, has been created and is now displayed in the command window. When defining an array other than a scalar, brackets [ ] should be used. The elements belonging to a same row are separated by a ”space” (or a comma), while semicolons are used to skip to the next line of the array. • Use the whos command, and check the information concerning the new variable M. The size of the array is displayed. It gives 2 × 4, which means that M has 2 rows and 4 columns. This information can also be retrived using the function size(): >> size(M) which returns a 1 × 2 array containing the number of rows and columns of M. • Type now >> u=[2; 4; 6; 3] What kind of array do you obtain? Type now [n_r n_c]=size(u); and check the values of n_r and n_c. What do they contain? • Vectors can also be created in the following way. Type >> k = [1:1:6], l = [1:2:6], m=[12:-3:0] and analyse the outputs. Here vectors are defined in the following form: [begin:interval:end], where interval can be any positive or negative number. Remark: When interval is equal to 1, it does not need to be specified. 1:6 returns the same vector as 1:1:6. • Similar commands can be applied to create matrices. Type: >> Z = [1:9; 3:11] and then press ”enter”. 14 • Type now >> Z = [1:8; 3:11] Why does it produce an error? Concatenation It can be useful to combine several related arrays in one unique array, for instance when two datasets are related and should be treated together. This can be done thanks to the operation of concatenation (=joining the arrays one after the other). • Type for instance: >> a = [1:3] >> b = [4:6] >> c = [a b] The space is used in this case to concatenate the arrays horizontally (alternatively, you can use the function horzcat(), e.g. c = horzcat(a,b)). • Type now >> d=[a ; b] In this case the arrays have been concatenated vertically (see also vertcat()). • Concatenation can be used to add an extra column (or line) to a matrix. Type >> e = [10;20] >> f = [d e] and analyse the outputs. • Vertical (horizontal) concatenations can only be performed if the arrays have the same number of columns (rows). Type: >> [f; a] Do you understand the error message? 15 Manipulating arrays by their subscripts Elements of arrays are usually identified by their position (row and column numbers), also called subscripts. The following table shows how the different elements of a 3 × 2 matrix can be identified using their subscripts: (1,1) (2,1) (3,1) (1,2) (2,2) (3,2) • Type for instance >> M(2,4) It returns the element located at the 2nd line, 4th column. Type now >> M(3,4) Do you understand the error message? • This notation can be used to alter a value of an array. Type for instance: >> M(1,2)=0 • Several elements can be selected using a vector instead of a scalar subscript. Type for instance >> M(2,[1 3]) The first and third elements of the second row of the matrix M are displayed. • To create a vector containing the entire row, type for instance: >> vectM = M(2,:) Similarly M(:,2) would refer to the second column of M. • Replace the elements of the first column of M by 10. Manipulating arrays by their indices 16 • For a vector, it is easier to define the position of an element using only one number, called its index, which is its position in the row or the column. For instance, vectM(2) is equivalent to vectM(1,2). Several elements of vectM can also be selected using a vector containing the relevant indices. Execute for instance: >> vectM([1:3]) • Similar conventions can be used for matrices. The elements are then numbered column by column, starting from the top left element and progressing downward. For a matrix of dimension 3 × 2 (note that the index of a given element will depend on the size of the matrix), we have: (1) (2) (3) (4) (5) (6) Type the following commands and compare the output: >> M(5) >> M(1,3) • How would you call the element M(1,4) using indices? 1.3.3 Operations on arrays Operations between a scalar and an array • Multiplication by a scalar: All the elements of the array are multiplied by the scalar. Type for instance : >> 2*M and check the output. Remark: 2*M is equivalent to M*2. 17 • Addition (subtraction) with a scalar: Type >> M-2 The scalar is added (subtracted) from each element of M. Operations between two matrices • Addition (subtraction) of two matrices. These operations can only be applied to matrices of identical size. The resulting matrix is obtained by adding (subtracting) their corresponding elements. Example: >> M-f • Multiplication of two matrices. This operation is executed according to the rules of linear algebra. Let’s define two matrices A and B such as: 1 3 1 0 2 A= and B = 4 1 . 5 10 7 2 2 The element (i, j) of the matrix C = A ∗ B is defined as : C(i, j) = n X A(i, k) ∗ B(k, j), k=1 where n is the number of columns of A (=number if rows of B). C is therefore equal to: 1∗1+0∗4+2∗2 1∗3+0∗1+2∗2 5 7 C= = 5 ∗ 1 + 10 ∗ 4 + 7 ∗ 2 5 ∗ 3 + 10 ∗ 1 + 7 ∗ 2 59 39 18 Operation between two matrices Addition Matlab expression A+B Multiplication A*B Element-by-element multiplication Element-by-element division Condition size of A = size of B number of columns of A = number of rows of B A.*B size of A = size of B A./B size of A = size of B Table 1.3: Summary of the main operations between two matrices A and B. The requirements in term of matrix size are also indicated. Element-by-element operations • To perform a multiplication element-by-element, it is necessary to use the operator ”.*” instead of ”*”. Type for instance >> M.*f and examine the result. Type now the same command but without the dot. Do you understand the error message? (see conditions in Table 1.3). Similarly, element-by-element division is obtained using the operator ”./”, and element-by-element exponentiation using ”.^”. • These operations can be combined. Type for instance: >> MM = vectM.^2 + vectM 1.3.4 Matlab built-in functions for arrays Elementary mathematical functions The functions defined in Table 1.2 can also be applied to arrays. In this case, the function is applied to each matrix element individually, and returns an array of the same size as the input 19 argument. Compare the outputs of the two following commands: >> M/3 >> round(M/3) Transposition Enter the following command >> v = transpose(u) and compare v to u. Now look at transpose(M) and compare it with M. Do you understand how the transposition function works? transpose() changes the rows of the matrix into columns, and vice versa. As a consequence, a matrix of size n × m will turn into a matrix of size m × n containing the same elements organised in a different way. This operation can also be performed using the symbol ’ just after the name of the array. u’ is therefore equivalent to transpose(u). Creation of arrays The functions zeros() and ones() can be used to create a matrix. Type for instance: >> C = ones(3,5) Do you understand the meaning of the arguments of this function? As expected, the function zeros() is similar to ones() but returns a matrix filled with zeros. • With the help of the functions introduced previously , create a matrix of size 20 × 30 containing only 5s. Data analysis The main Matlab functions available to perform data analysis are summarized in Table 1.4. Some examples of the use of these functions are given below. • Type the following commands: >> mean(vectM), median(vectM), sum(vectM) Each of these functions, applied to the vector vectM, returns a scalar, which is equal to (respectively) the averaged value of the vector , its median value and to the sum of its elements. • Type: >> vmax =max(vectM) 20 vmax is the largest value of the vector vectM. In some practical cases, it is important to know when the maximum value has been reached, that is to say what is the position of the largest value in the original vector. In order to get this information, the function max (or min) can be used with two output arguments. Type >> [vmax ind]=max(vectM) The second argument (ind) is the index of the maximum value in the vector. The function min works in a similar way. • Type now: >> s = sort(vectM) Compare s and vectM. • These functions can also be applied to matrices. In this case, they handle each column separately. The output is therefore a vector of length the amount of columns. Type for instance: >> mean(M) and analyse the result. When M is a matrix, mean(M) is a row vector containing the mean value of each column. • It is generally possible to specify the direction along which the computations have to be performed as an additional input argument. Type for instance: >> mean(M,1) >> mean(M,2) and compare the results. Do you understand the meaning of the second input argument? The functions median, min, max and sum work in a similar way. You can check the Matlab help to find out how to choose the direction of computation for these functions. • An alternative approach is to apply the functions after transposition of the matrix. Type for instance: >> mean(M’) and compare with mean(M,2). 21 max mean median min std var prctile sum sort Maximum value Average or mean value Median value Smallest value Standard deviation Variance Percentile values (the percent of interest should be specified) Sum of the elements Arrange the elements in ascending order Table 1.4: Some important Matlab functions for data analysis EXERCISE 2 Make a new file called Exercise2.m in your work-directory, where the commands needed to solve the following assignment will be written. 1. Create the following matrix: 16 5 A= 9 4 3 10 6 15 2 13 11 8 7 12 14 1 2. Define a vector A1 containing the first row of A. Define a scalar, sumA1, equal to the sum of its elements. 3. Create now a vector sumRows P4 where each element contains the sum of the elements of one of the row of A. For instance, sumRows(1) = k=1 A(k, 1). Can you do it in only one line of command? 4. Create a vector containing the sum of the column of A. What do you notice? 22 1.4 Visualising data How to make a simple 2D-plot : example • Define the following vectors: x = -pi:.1:pi; y = sin(x); y2 = cos(x); • Type figure,plot(x,y) The first command, figure, opens a new figure window, called Figure 1 by default since no other figure has been opened previously. The plot command is then executed in the newly created figure. • Type now plot(x,y2) The new plot has replaced the previous one in the figure window. To display the second graph on top of the first one, the command hold on should to used. Type: figure,plot(x,y) hold on,plot(x,y2,’r’) A new figure, called Figure 2 by default, is created, and contains the two plots on top of each other. For the second plot, the line color is specified using the additional argument ’r’, with r standing for red. • Title, labels on the x- and y-axis and legends can now be defined. Type the following commands one by one and analyse their effects on the figure: title(’Example plot’); xlabel(’Radians’); ylabel(’Function Value’); legend(’sine’,’cosine’); 23 • To specify the limits of the x- and y-axis, type: for instance axis([-pi pi -2 2]) Do you understand the meaning of these 4 numbers? • The figure command can also be used with an input argument which is a positive integer, as for instance figure(1). The number corresponds to the number of the corresponding figure window. If the Figure 1 does not exist yet, this command will create it. If this figure already exists, it becomes the active figure window. This means that new plot commands will be executed in this figure. Type for instance: figure(1), hold on, plot(x,-y,’.’) Because of the argument ’.’, the vectors are plotted as a cloud of points. • plot can also be used without specifying the x argument. Type: figure,plot(y) The elements are now plotted as a function of their position in the vector: y(1) at x = 1, y(2) at x = 2, etc. • Multiple plots can also be displayed in one figure. This can be done using the command subplot(a,b,c) where a, b and c are integers, and c is less than or equal to the product a*b. This command divides the figure in a*b sub-figures, organized in a rows and b columns. To obtain 2 figures on top of each other, the following lines should be used: figure subplot(2,1,1);plot(x,y) subplot(2,1,2);plot(x,y2,’g’) • It is possible to display on the figure the coordinates of a given data point using the ”data cursor” button presented in Figure 1.2. Select this button and click on a data point. • You can now close all the figures using the command: close all The function plot is the basic Matlab function for 2D plots. Numerous options are available to customize your plots. Some examples are given below: 24 - plot(x,y,’r’): red solid line - plot(x,y,’--r’): red dashed line - plot(x,y,’--r’,’linewidth’,2): thick red dashed line - plot(x,y,’*g’): the points are marked with green asterisks (no line) - plot(x,y,’+-r’): solid red line with points marked with a plus sign. The tables below show some of the colors and markers available. We advise you to look at the Matlab help (type doc plot for instance) for a more extensive description of the plotting options. You can also navigate in the help window to learn more about the numerous ways to visualise data in Matlab. red green r g plus sign asterisk Color specifiers blue b black k cyan c yellow y + * Marker specifiers square s triangles v, < or > magenta white circle cross m w o x EXERCISE 3 This exercise is based on the dataset contained in the file MPMtransportdata.xls. Meyer-Peter and Mueller (1948) (abbreviated MPM) graphically report gravel transport data and derive their famous empirical bedload transport predictor from this dataset. Wong and Parker (2006) recovered the original data and re-analysed it to find that the original fit by MPM was wrong. Here you will plot the data in nondimensional form and calculate a predicted transport. The Kleinhans 2005 paper contains all the equations necessary to solve this exercise. Remark: We assume that the density of the water is 1000 kg/m3 . 25 Figure 1.2: Figure window Remarks : • It is recommended to test parts of your script regularly. An easy way to do it is to copy-paste these parts into the command window and then press ”enter”. Alternatively, you can select the part of the text you want to run and then press ”F9”. • Always terminate the command by a semicolon if you do not want the result to be displayed in your command window. This is important when working with large datasets. 1. Make a new script file called Exercise3.m in which the series of commands necessary to solve this assignment will be written. Do not forget to clear your workspace and close the figures at the beginning of your script. 2. Download the file MPMtransportdata.xls and save it in your current work-directory. You can now load the data using the following line: 26 dataMPM = xlsread(’MPMtransportdata.xls’); Thanks to the function xlsread(), Excel sheets can be read and loaded in your workspace. The output of the function is an array called dataMPM which contains the same information as the Excel sheet, but without the description of the data (only numbers). Check its number of rows and columns and compare a few values with the initial Excel sheet. 3. Define vectors containing the channel width W, the depth h, the slope S, the median grain size D50, the specific gravity of the sediment s and the sediment transport rate qs. 4. From these vectors, calculate the total shear stress. The result will be a new vector, called tau for instance. 5. Calculate the Shields parameter using the shear stress previously defined. 6. Calculate the Einstein parameter. 7. Plot the sediment transport rate (y-axis) as a function of the shear stress (x-axis). Plot the data in log-scales for both the x- and y- axis. To do that, the function plot() can be replaced by loglog(). loglog is basically the same as the plot function, but with logarithmic axis, and can be used with the same input arguments (see Matlab help for more details). Give an appropriate title to the plot, as well as names to the x- and y-axis. HINT: To write Greek letters in the text of a Matlab figure (axis labels, title, legend, etc.), use the following characters: α β \alpha \beta τ Φ θ \tau \Phi \epsilon \theta φ λ \phi \lambda To write subscripts or superscripts, use the symbols _ and ^. For example, to display qs (m2 /s), you should write q_s (m^2/s), and to display D50 (2 characters as subscript), use D_{50}. 27 8. Plot in a new figure the data in a nondimensional form, i.e. the Einstein parameter (y-axis) as a function of the Shields parameter (x-axis). Give an appropriate title to the plot, as well as names to axis. Compare with the previous figure : does the adimensionalization seem efficient? 9. Plot the MPM prediction as a additional red line in the second figure. To plot this line, first define two vectors x and y such as x contains a range of values reachable by the Shields parameter (θ) and y = 8(x − θcr )1.5 (see eq. (30) in Kleinhans 2005). We use θcr = 0.047. 10. In the same figure, plot a vertical line corresponding to x = θcr . 11. To further analyse the data, we will now remake the figures of the questions 7 and 8 using the function scatter(). scatter is a function which allows you to plot the data as circles of different sizes and colors depending on their characteristics. Check the Matlab help to understand how to use the function scatter() and what are its input arguments. The lines below are an example of use of the function: figure; scatter(tau,qs,10,s,’filled’); set(gca,’xscale’,’log’,’yscale’,’log’); colorbar; Execute these lines one-by-one and analyse the outputs. How many colors are displayed in the figure and why? do you understand why the size of the circles does not vary? 12. Create a new figure, which will be divided in two subplots. - In the top subplot of the figure, plot sediment transport (y-axis) as a function of shear stress (x-axis), with colors depending on the grain size D∗ (Bonnefille number). - On the bottom part, plot the same figure but in its nondimensional form. Add colorbars to both plots. For a better distribution of the colors, it is recommended to use the logarithm of the Bonnefille number as the ”color vector” argument of the function scatter. 13. Look at the data plotted in non dimensional form. One set of points seems to follow a different trend than the others. To which range of D∗ does it correspond? 28 1.5 Final assignment : Analysis of the Rhine flow discharge measured at Lobith – part I In this exercise, we will perform some analysis on a dataset containing the flow discharge of the Rhine measured at Lobith (where the Rhine enters the Netherlands) for the period 1901-2000 AD. The dataset is given in the text-file Qlobith19012000.txt as an array. Each element of this array is the discharge in m3 /s for 1 day, and each column corresponds to one year (first column = 1901, last column = 2000). Preliminary analysis 1. Create a new script file called Lobith1.m in which the assignments of the exercise below are executed. Start the script by clearing your workspace and closing the windows. 2. Load the new dataset. As the data are stored as a text (ASCII) file, the following command should be used: Q=load(’Qlobith19012000.txt’); Check the size of the array Q. Do you understand how it is organised? 3. The data corresponding to the 29th of February are included for each year, leading to a number of lines of 366. Display the 20 first elements of the row corresponding to the 29th of February. What do you notice? Remark: NaN stands for ”not a number”. It is often used to represent missing values in datasets. 4. To obtain an overview of the dataset, use the following command: figure,plot(Q) In this command, the function plot is applied to a matrix instead of a vector. plot handles each column (and therefore each year of data) separately. As a result, each line appearing in the figure corresponds to the evolution of the discharge for a given year. Give a title to the figure and names to the x- and y-axis. 5. Zoom on the figure around the 29th of February (you can use the ”zoom” buttons presented in Figure 1.2 or redefine the x-axis limits using xlim()). How are the NaN-values handled by the function plot? 29 6. Create a vector called year containing the years when the data were collected, as well as a vector days containing the day numbers (1=first of January). Analysis of one year of data 7. Create a new vector called Lobith1916 containing the flow measured at Lobith in 1916. 8. Determine the mean and median discharges for this year. Compute also the maximum and minimum values of the discharge. When were these values reached? 9. Compute the mean discharge for the following year (1917). What do you obtain? To understand better what is happening, check the results of the following operations: >> NaN + 10 >> 20 * NaN This explains why we cannot compute directly the mean value for a year without 29th of February. 10. To ignore the NaN values in the computations, alternative functions can be used, such as nanmean, nansum, etc. Use the Matlab help for a description of these functions and compute the mean discharge in 1917. Statistics over the entire data set 11. Build a vector of size 1 × 366, where each element is the discharge for a given day averaged over 100 years. 12. Plot this vector as a function of the time. When are the highest mean discharges occurring (which month?)? Plot in the same figure the maximum et minimum values for each day (use different type of lines to differentiate the curves). Add names to the axis and a legend to the graph. 13. Create a vector containing the annual flood discharges for the 100 years of data. 14. Compute the mean annual flood discharge over the entire period. 30 15. Compute the median value of the annual flood data as well as the 75th, 90th, 95th percentiles (see http://en.wikipedia.org/wik for definition). 16. We will now use the annual flood discharge data to build a flood-frequency curve, which is a graph showing the relationship between flood magnitude and their recurrence interval for a specified site. The recurrence interval (T r) of a flood is a statistical measure of how often a flood of a given magnitude (QT r ) is likely to be equalled or exceeded. For instance, the ”fifty-year flood” (Q50 ) is one which will, on the average, be equalled or exceeded once in any fifty-year period. Follow the following procedure to compute the recurrence interval: - Create a vector QT r containing the annual flood discharges sorted in ascending order. - Define a vector R of same length as QT r containing the ”ranks” of the discharges. The largest discharge of QT r should have a rank equal to 1, and the smallest a rank equal to N, where N is the number of years for which we have flood data. - The recurrence interval T r (in years) can then be computed using the formula: T r = (N + 1)/R. You can now plot QT r versus T r : this is the flood-frequency curve. Use logarithmic axis for a better display of the data. Add a title and labels to the axis. 17. Read in the figure what is the magnitude of the fifty-year flood (use the data cursor button indicated Fig. 1.2 to help you). What about the 100-year flood? (the answers should be written as a comment in your script file). Remark Numerous Matlab functions are available to study the distribution of data and the occurrence of a given event. Histograms can be built using the function hist, the distribution of the data can be compared with different probability distribution with the help of probplot, etc. 31 PRODUCTS TO HAND IN Exercises 1 to 3 and final assignment (Section 1.5). Please note that only the relevant information should be displayed while running the scripts ⇒ use semicolons to avoid the display of endless lists of numbers! 32 Chapter 2 Tests, loops and functions 2.1 Relational and logical operators Relational operators A relational operator (see Table 2.1) compares two numbers (or arrays) to determine if a comparison statement (e.g., a < b) is true or false. • If the comparison statement is false, it returns a zero. Type for instance: >> 3<2 If the statement is true, such as >> 3<6 it returns 1. • Such comparisons are generally performed on arrays and can be very useful to analyse datasets. Define an array A such as: >> A = [1:2:10;10:-2:1] 33 and then type >> E = (A<5) This command returns an array E which has the same size as A, but contains only ones and zeros. Each element of A is compared to 5, and, depending on the outcome of each comparison, the elements of E are filled with 1 or 0. • Use the function sum to determine how many elements of A are smaller than 5. This method can be very useful to characterize large datasets where it is not possible to count one-by-one the elements satisfying a certain condition. • Relational operators can also be used to compare two arrays. Type: >> B = [3 7 1 8 3; 1 45 0 2 9] >> F = (A<B) and analyse the result. The comparison between A and B is done element-by-element. The results are stored in a third array, F, which has the same size as A and B. Each element of F is equal to 1 if the comparison is true at this location, and equal to 0 otherwise. • These operations can be performed using any operator described in Table 2.1. Type for instance: >> G=(A==B) Do you understand the output? • Arrays such as E, F and G, resulting from relational and logical operations, are called ”boolean” or ”logical” arrays. They can be used to select elements from another array. Type: >> a = A(E) A new vector, a, is defined, and contains the elements of A corresponding to the locations where the ”logical” matrix E was equal to one. It therefore contains the elements of A smaller than 5. The matrix E can be used to extract elements from any other matrix, as long as all the matrices have the same size. Type for instance: >> b = B(E) and interpret the result. 34 Description Less than More than Less than or equal to More than or equal to Equal to Not equal to Symbol < > <= >= == ∼= Table 2.1: Description of the different relational operators and their symbols in Matlab Of note, it would have been possible to write directly (without defining separately the logical array E): >> a = A(A<5) >> b = B(A<5) Logical operators 2.2). Comparison statements can be combined using the logical operators AND and OR (see Table • Type for instance: >> H=(B>2)&(B<10) H contains ones only if the elements of B are larger than 2 and smaller than 10. The two comparison statements should be true at the same time to return a 1. • When the operator OR is used, the output is true if either one, or both, comparison statements are true. Replace the & in the previous command by the OR symbol |. Did you expect the output? • A third logical operator, NOT, can be used. Applied to a logical array, it returns its opposite. Type for instance: ∼ H, and compare with H. 35 Description OR AND NOT Symbol | & ∼ Table 2.2: Logical operators Use of the function find The function find looks for the non-zeros elements of an array, and returns their indices. If no elements are found, it returns an empty matrix. • Type the following commands: >> u=[32 6 4 9 -6 14] >> v=(u>7) >> ind = find(v) find returns a vector, ind, containing the indices of the non-zeros values in v, that is to say the indices of the elements of u larger than 7. Again, the two last lines could have been written in one time as: >> ind = find(u>7) • Type now: >> ind2 = find(A<5) We obtain the indices of the elements of A smaller than 5 (see definitions of indices for matrices page 16). Compute now A(ind2) and compare it with A(E). Do you understand the result? EXERCISE 1 Open a new script file where the following assignments will be solved. The matrix A should be defined at the beginning of the script. 1. Compute the average value of the elements of A larger than 2 and smaller than 6. 2. Replace the elements of A smaller than 3 by 0. 36 3. Replace the elements of A larger or equal to 9 by 2 times their value. 2.2 Flow control In a simple program, the commands are executed one after the other, independently from the output of the previous one. However, it is sometimes necessary to include some sort of decision making in the program. For instance, we can decide to compute the sediment transport using different predictors depending on the sediment characteristics and flow conditions. The program should therefore be able to choose the appropriate calculation and skip the others. It can also be needed to repeat a sequence of commands several times, until a certain condition is reached for instance. The order in which the commands or group of commands are executed, also called the flow of the program, can be controlled in several ways in Matlab. In this section, the three basic flow control methods will be presented. 2.2.1 if-statements Purpose if-statements are used when a decision must be taken depending on the outcome of a comparison statement (also called conditional expression). If the comparison is true, then the group of commands are executed, and if it is not, it is skipped. Format The most simple if-statement is the if-end construction, as in the following example: x=2; if x<5 disp(’x is less than 5’) end • Execute this small program and interpret the result. • Replace x=2 by x=5 and re-run the program. 37 When the program has to choose between two alternatives, a if-else-end-constructions can be used. Read carefully the piece of code presented below (please note that in this example Theta is a scalar): if Theta<Theta_cr %no sediment transport Phi = 0; else %sediment transport computation Phi = 8*(Theta - Theta_cr)^1.5; end Based on the outcome of the comparison statement Theta<Theta_cr, either the commands between the if and else or the commands between else and end are executed. This small program is therefore simply expressing the fact that if the Shields parameter Theta is below a critical value, there is no sediment transport, and that otherwise the sediment transport is calculated using the MPM predictor. The general structure of the if-statement is presented below: 38 if expression1 group of commands evaluated if expression1 is true elseif expression2 group of commands evaluated if expression2 is true elseif expression3 group of commands evaluated if expression3 is true . . . else group of commands evaluated if none of the previous expressions is true end Only the commands associated with the first true expression are evaluated, the others are skipped. The conditional expressions can be defined using the operators presented in §2.1. 2.2.2 for-loops Purpose for-loops are used to repeat a command or a group of commands for a fixed (predetermined) number of times. Format The general structure of a for-loop is presented below: for index = start:increment:end commands to be executed in the for-loop end 39 Examples • Let’s start with a simple example to understand how these loops are working. Run the following lines and interpret the results: for j=1:8 disp(’Hello’); end If necessary, check the Matlab help to learn more about the disp function. • Loops can be for instance used to define arrays. Run now the following script and examine carefully its results: clear for j=1:9 s(j) = j^2; end s The loop is executed for each element of j, which starts at one and increases by 1 until it reaches 9. At each step, the square j is computed and stored in jth element of the s. In other words, this small script computes the squares of the numbers from 1 to 9. • Note that this for-loop could be replaced by the following lines: clear j=[1:9]; s=j.^2 These two ways of solving the problem are equivalent. However, to reduce the computational time it is recommended to avoid loops when it is possible. • Loops can also be nested, as in the following example, where we define a new matrix B and fill it in element-byelement: 40 clear B=zeros(2,3); for i=1:2 for j=1:3 B(i,j) = i*j end end In this case, every time i increases by one, the nested loop is executed 4 times (from j = 1 to j = 4). Run this script to check in which order the matrix is filled in. Do you understand why? Remark In the previous example, we first defined B as a 2 × 3 matrix of zeros, and then filled it in within the loop. This initialising step is not compulsory in Matlab (comment the second line and re-run the script), but it makes your program more efficient. Indeed, it is a way to pre-allocate memory: Matlab ”knows” from the beginning the size of the matrix you are building up, and does not need to re-size it at each iteration. The program runs faster (important when working with large arrays!). • for-loops can be combined with if-statements, as in the following example clear s = [2 1 -4 23 0]; for i=1:5 if s(i)>=2 s(i) = 3*s(i); end end s Run this script and interpret the output. In this case, the combined use of the for-loop and if-statement could have been replaced by the use of relational 41 operators. Do you understand how? Some remarks about for-loops - The first time through the loop, the variable is assigned the first value in the vector [start:increment:end]. For the second time through, MATLAB automatically assigns to the variable the second value in the vector, etc. - The number of times a loop is executed equals the number of values in the vector. - After the loop is finished executing, the loop variable still exists in memory and its value is the last value used inside the loop. EXERCISE 2 This exercise deals with the computation of the evolution of Carbon 14 in a sample due to radioactive decay. Each thousand of years a fraction α of C14 is lost. This means that, calling C0 the initial mass of C14, the mass remaining in the sample will be: - C1 = C0 (1 − α) in 1 thousand years, - C2 = C1 (1 − α) in 2 thousand years, - C3 = C2 (1 − α) in 3 thousand years, etc. We want to calculate the quantity of C14 still in the sample after 25000 years. Write a script which creates a vector C such as its first element C(1) is the initial amount of U238, C(2) the amount after 1 thousand years, etc. The last element of the vector should contain the amount of C14 after 25 thousand years. You can use the following script lines as a basis and complete them. Data: C0 = 0.5 kg; α = 0.117; 42 % Initialization ??? ??? % Calculation of the radioactive decay for i=2:26 ??? % Computation of the i^th element of the vector C end % Definition of a time vector in thousand years t = (0:25); % Visualisation figure; plot(t,C/C(1)); ylabel(’C(t)/C(0)’); xlabel(’t (10^3 years)’); EXERCISE 3 Let’s come back to exercise 3 of the last chapter (page 25) concerning the analysis of the dataset from Meyer-Peter and Mueller. We observed that part of the data was following a different trend than the others after non-dimensionalisation. A close look at it shows that these data points correspond to experiments performed using a very fine sediment (D50 = 0.3402 mm). In this case, bedforms probably developed, and the use of a grain-related shear stress instead of the total shear stress is more appropriate and might improve the prediction. 1. Copy the data file MPMtransportdata.xls in your current work directory and open a new script file. 2. Create a new vector tau2 containing the total shear stress for data such as D50 ≥ 1 mm and the grain-related shear stress otherwise (D50 < 1 mm). You can use parts of the script written during the first tutorial. Remarks: - this question can be solved using if-statements and loops or relational/logical operators only. Choose the 43 most convenient way for you! - in Matlab, the function calculating the base-10 logarithm is log10(). log() is the natural logarithm. 3. Calculate the associated Shields parameters. 4. Plot the data (using the function scatter) as the Einstein parameter as a function of the newly computed Shields parameter. Compare with the figure created during the previous tutorial. 5. Additional question for those who found the exercise easy. Now also remove wall friction and explain how you did it. See the Van Rijn fluid mechanics book or the paper by Wong and Parker 2006, J. of Hydraulic Engineering, doi 10.1061/(ASCE)0733-9429(2006)132:11(1159). 2.2.3 while-loops Purpose while-loops are used to repeat a command, or a group of commands, an indefinite number of times until the condition specified by while is no longer satisfied. Format The general structure for a while-loop is the following: while expression group of commands repeated as long as expression is true end When the program reaches the first line of the block, the conditional expression is checked. If it is true, Matlab executes the group of commands. Then Matlab jumps back to the while command and check again the expression: if it is still true, the commands are executed again. The looping process will stop when the expression becomes false. 44 Examples • Execute the following script and interpret the output: clear x=1; while x<14 x=x+4 end How many times are the commands executed and why? • Execute now: x=1; while x<14 x=x-4 end Do you understand what is happening? Remark: To terminate a computation, press Ctrl+C. This command can be very useful in case of programming errors leading to never-ending loops. Some remarks on the use of while-loops - The variables used in expression must be defined BEFORE the ’while’ loop, so you can use it to enter the loop - At least one variable used in expression must change INSIDE the while loop, or you will never exit the loop. - After the loop is finished executing, the loop variable(s) still exist in memory and its value is the last value the variable had in the while loop. 45 EXERCISE 4 In this exercise, we aim at characterizing the flow in typical experimental conditions. The channel is 0.3 m wide and has a slope of 0.006 m/m. The Nikuradse roughness length is 0.0013 m. We try to answer to the following questions : how deep will the water be when the water flow is 1 × 10−3 m3 /s? what will then be the flow velocity? This problem will be solved using an iteration process. We start with a first guess of the water depth, and correct it at each loop to make it more accurate. The looping process will stop when the desired accuracy has been reached. • Create a new script where this problem will be solved. Some help concerning the structure of the program can be found below: Initialisation - Make a first inaccurate guess to define h (should be such as 0 < h < hmax , with hmax = 3 cm) - Define also htemp to allow for the first loop to start, e.g. htemp = 0.9 ∗ h. As long as the absolute value of (h−htemp ) is larger than 1×10−6 - New guess for h : the ”real” depth should be contained between h and htemp . We try the value (h+htemp )/2 (other options are possible!). - Computation of u - Computation of htemp End of the loop A good estimation of h has been found. Calculate the associate flow velocity. Some additional comments Each iteration starts with defining a new, presumably more accurate, value of the flow depth h. We then compute the flow velocity u associated to h using the Chezy equation. This flow velocity is then used to compute a second estimation of the water depth called htemp , based on the mass conservation 46 equation. h and htemp are then compared. They should be nearly equal if our last estimation of h was good enough. • To check your code, run the script with different guesses for the initial water depth. If the while loop is correctly coded, it should always return the same final value for the depth h. 2.3 User-defined functions You have already encountered functions that are built in to MATLAB (see Tables 1.2 and 1.4 for instance). As you start to write your own programs in MATLAB you will need to create your own functions that take one or more input arguments, operate on them, and return the result (output). Function files are simply text files with a *.m extension, and can be written, as script files, using Matlab editor. They are organized in the following way: - The first line of a function m-file must always be of the following form: function [ output variables ] = function name( input variables ) This line always starts with the word function. Following that, the output variables are enclosed in square brackets [ ]. Then, the characters function_name define the name which will be used to call the function. The input variables are finally specified into brackets. NB: The function function_name should be saved in a file called function_name.m. - The definition line of a function file must be followed by a block of comments, containing at least the following fields: a short description of the function aim as well as definitions of the input and output variables (see example below). - The body of the function can then start. It is made of a succession of commands which ultimately defines the output variables. 47 • Let’s start with an example of a simple function. Copy the following lines in a new m-file and save it in your current directory as cylinder_charact.m: function [S,V] = cylinder_charact(R,h) % [S,V] = cylinder_charact(R,h); % Computation of the volume and surface area of a cylinder % input R radius (m) % h cylinder height (m) % output % S surface area (m2) % V volume (m3) a_lateral = 2*pi*R*h; a_circles = 2*pi*R^2; S = a_lateral + a_circles; V = pi*R^2*h; In this example, a function called cylinder_charact has been defined. The definition line is followed by a block of comment lines, which are meant to provide information about the function and its arguments. These lines are important since they allow the function file to be referenced by the Matlab help. • Type now in the command window >> help cylinder_charact The block of comments following the definition line is displayed. • Clear your workspace and type the following commands: >> [S1,V1] = cylinder_charact(2,3); The function cylinder_charact() is called with the input arguments 2 and 3. It returns the surface area and volume of a cylinder of radius 2 m and height 3 m. These values are assigned to the variables S1 and V1. • Type whos in the command window. The intermediate variables (=variables not defined as outputs), such as 48 a_lateral, a_circles only exist within a function and do not appear in the workspace. The use of functions is therefore important since it avoids to overload the workspace with variables which are not going to be used later. • These functions can be used in a script file, as any built in function in Matlab. Modify the function file so as it allows for non-scalar input arguments. You should then be able to run the following script: % Computation of the volume and surface area % of cylinders of increasing radius % M. Tissier, 3-09-12 % Initialisation clear h=3; % constant height of the cylinders (m) r=(1:6); % increasing radius of the cylinders (m) % Computation of the volumes S_tot and areas V_tot % of the cylinders [S_tot,V_tot] = cylinder_charact(r,h); figure; plot(r,S_tot,’o-’); xlabel(’radius (m)’);ylabel(’surface area (m^2)’) Some (additional) reasons to use functions in your scripts - Script files in which all the calculations are entered in one unique file can become too long and difficult to understand. By using functions files, you can organize better your program, and make it easier to read and understand. - A function can be called several times within a script with different input arguments: you do not need to write similar calculations over and over again... 49 EXERCISE 5 We now continue from exercise 3 of Chapter 1. Here, the objective is to compare the measured sediment transport to the predictions given by the MPM and EH formulations, using several user-defined functions. 1. Create a function calculating the (total) Shields parameter. The input parameters should be the channel slope, the densities of the water and the density of the sediment, the water depth and the median grain diameter. 2. Create two functions which calculate the non-dimensional sediment transport according to MPM and Engelund and Hansen (EH) relations. The input arguments for MPM should be the Shields parameter and its critical value. Choose the adequate arguments for EH. HINT: For the MPM predictor, do not forget to account for the fact that sediment is only transported if the Shields parameter is higher than a critical value (θcr = 0.047). 3. Create a function converting the non-dimensional transport into the dimensional one (in m2 /s). 4. Write a script which will call the functions defined above and ultimately plot the measured sediment transport against the predicted one for both MPM and EH equations. Add a line corresponding to x = y on the figure. Label the x- and y-axis and add a legend. 2.4 Final assignment : Analysis of the Rhine flow discharge measured at Lobith – part II Morphologically relevant statistics of the discharge dataset 1. Calculate the maximal and minimal daily discharge for the whole dataset. 2. Calculate the percentage of days of year 1916 where Q exceeds the mean annual flood discharge calculated in Section 1.5. 50 3. Create a function which takes as an input argument a vector and returns its number of non-NaN elements (HINT: check Matlab help for isnan and length). 4. Compute the same percentage for each year of data and store the results in a vector. Remark: For some years, there is one day less than for others. Use the function defined in question 3 to take it into account while calculating the percentages. 5. Plot the percentage of days where Q exceeds the threshold value as a function of the year (give a title and names to the axis). 6. Compute finally the same percentage but for the whole data set (over the 100 years of data). Computation of the discharge in the channel during floods 7. We consider here that the bankfull discharge, Qbankf ull , can be approximated by two times the mean value of Q over the entire dataset. Calculate Qbankf ull . 8. Calculate the number of days where the bankfull discharge has been exceeded over the 100 years of data. When the total discharge Q gets larger than the bankfull discharge Qbf , the river floods. In this case, the discharge really flowing in the channel and contributing to the channel bed shear stress is less than the total discharge Q: Qchannel = Qbankf ull + Qoverbchannel (see Figure 2.1 for definitions). The following question aims at calculating Qchannel , which is the relevant discharge for sediment transport calculations. 9. Define a matrix Qchannel of the same size as Q. The elements of this matrix should be equal to those of Q if Q < Qbf , and equal to Qbankf ull + Qoverbchannel otherwise. To check your calculations, compare Qchannel and Q. Does it match your expectations? HINT: Qoverbchannel can be expressed as a function Qoverbank , W and Wf p . To find this relation, relate the cross-sectional area Aoverbchannel to Aoverbank . As the depth-averaged flow velocity is considered to be constant 51 Figure 2.1: Skematic view of the profile of the Rhine at Lobith during a flood event. Wf p is the total width of the flood plain, W the channel width defined earlier. The total discharge Q in the river is equal to Qbankf ull + Qoverbank . The part of the overbank discharge flowing over the channel is called Qoverbchannel . 52 over the cross-section, you can then obtain a relation between the discharges Qoverbchannel and Qoverbank . Data: W = 550 m; Wf p = 3000 m; Sediment transport calculations Additional data: Channel slope S = 1.5 × 10−4 m/m; Median diameter of the sediment D50 = 2 × 10−3 m; Density of the sediment ρs = 2650 kg/m3 ; Water density ρw = 1000 kg/m3 ; 10. Calculate the water level from the discharge using the QH (stage discharge) relation given below: H = (384.73 ln(Q) − 1949.3)/100. H is the elevation in meters above mean sea level and Q the total discharge. 11. Calculate the water depth h knowing that the average bed elevation is 3 meters above sea level. 12. Apply the Engelund Hansen (1967) relation to compute the sediment transport rate associated to the channel discharge Qchannel . Use the functions defined earlier to calculate the sediment transport. 13. Calculate the sediment discharge Qs in m3 /s. Plot Qs as a function of Qchannel for the whole dataset. Choose the best representation (logarithmic axis or not). 14. Calculate the average volume (m3 ) transported through Lobith per year. To do that, you should first compute the total volume of sediment transported through Lobith for each year of the dataset, and then average these values. 15. Determine what value of flow discharge transports this same average sediment volume. This is a representative discharge for the morphodynamics and could be used in a numerical model for long-term morphological prediction instead of the full hydrograph, which is computationally more expensive. 16. Compare this representative discharge to the mean flow discharge at Lobith (calculated over the entire period covered with the dataset). 53 PRODUCTS TO HAND IN Exercises 1 to 5 and final assignment. 54 Chapter 3 Numerical modelling The answers to the questions of this practical should be organised in a short word document, in which you should include the relevant Matlab figures. Remark: To save a figure, use the Matlab function print or click on ”File” and then ”Save as” in the figure window (see Figure 1.2). 3.1 A very simple morphodynamic model We consider in this section a one-dimensional river with an initial constant slope S0 = 1/10. We want to compute the evolution of the river bed in space and time using very basic physical laws. Similar approaches are used in cellular automats of rivers and in so-called reduced-complexity models (try a search in Google Scholar to see debates and controversies). In this model, we do not consider the influence of the water motion on sediment transport: the sediment transport rate qt is assumed to depend linearly on the local bed slope S. More precisely, we assume that qt = SedP rop S, with SedP rop a proportionality coefficient. 55 To solve this problem numerically, it needs to be discretized in space and in time. The reach of the river under consideration, of length L, is divided into Nx sub-reaches of length dx = L/Nx . These sub-reaches are separated by Nx + 1 nodes, of coordinates xi = (i − 1)dx, with i integer varying between 1 and Nx + 1. To each of these nodes are associated a bed elevation and a sediment transport rate called ηi and qti (i is the node number). The bed elevation ηin at time tn is then used to compute the bed elevation ηin+1 at time tn+1 = tn + dt (concept of time discretization). This should be done in two main steps (see also Fig. 3.1): - the sediment transport rate qt is first computed as a function of the local slope at each node (proportionality coefficient given in the Matlab script). HINT : the local slope, defined as minus the gradient of the bed elevation ∂η (S = − ∂x ), can be computed using the built in function gradient (see Matlab help). - to compute the bed evolution, the sediment conservation equation (also called Exner equation) is then used. It can be written as: ∂qt ∂η =− , ∂t ∂x that is to say, in a discretized form: n ηin+1 − ηin ∂qt = , dt ∂x i and therefore ηin+1 = ηin − ∂qt ∂x n dt. i The gradient of sediment transport should be calculated and then used to update the bed elevation. • Download the script file very_simplified_model.m. It contains the structure of the model, and some commands to visualise the results. • Add in the initialisation part of the script the definition of the bed elevation vector eta at t=0. The initial bed elevation decreases linearly with a slope S0 . The bed elevation at x = L (downstream boundary) should be zero. Data: L = 10 m and S0 = 1/10. 56 Figure 3.1: Organization of our very simple morphodynamic model (no computation of the flow) • Complete the time-loop using the explanations given above. • Run the model and analyse the results. • Let’s now modify the upstream boundary condition. We set the sediment transport rate qt to zero at the upstream boundary (node 1, x = 0). This means that there is no incoming sediment flux in our model. Run the code and interpret the results. • Create a vector which will contain the time-evolution of the bed elevation at the first node. This vector should be filled element-by-element during the time loop. Plot this vector as a function of time. What kind of evolution do you obtain? 57 3.2 Flow computation: backwater equations To model the river in a more realistic way, it is necessary to compute the flow in the channel at each iteration, and then use it to compute the sediment transport and the bed evolution (see Fig. 3.3). In this section, we present the backwater equations, which will be used in Section 3.3 to compute the flow in a more complex morphodynamic model. More precisely, the flow in the channel will be computed using the backwater formulation for steady, gradually varied flow under the constraint of constant water discharge per unit width qw and an imposed water elevation ξd at the downstream boundary. Let H denote depth, ξ the water surface elevation, η the bed elevation (ξ = H + η, see also Fig. 3.2), kc the composite roughness height and S the bed slope. The backwater equation writes: dH S − Cf Fr 2 (3.1) = dx 1 − Fr 2 where S is the local bed slope and Fr the Froude number defined as Fr 2 = coefficient Cf is assumed to follow a Manning-Strickler resistance relation: Cf −1/2 = αr H kc U2 gH with U = qw H . The bed friction 1/6 (3.2) with αr = 8.1. How to solve this equation depth at the ith node): The backwater equation is first re-written in the discrete form (Hi being the water Hi+1 − Hi = Fback (H), dx with Fback a function defined as Fback (H) = S − Cf Fr 2 1 − Fr 2 58 (3.3) that is to say, using the definitions of Cf and Fr , Fback (H) = ηi+1 −ηi dx − 1 αr 2 1− H kc −1/3 2 qw gH 3 2 qw gH 3 . (3.4) Rearranging equation (3.3), we get: Hi = Hi+1 − Fback (H)dx. (3.5) In practice, the resolution is performed in two steps in the code using a predictor-corrector scheme: - a first estimation of the water depth Hpred i is computed: Hpred i = Hi+1 − Fback (Hi+1 ) dx - Hpred i is then used to obtain a more precise prediction of the water depth (corrector step): Hi = Hi+1 − 1 (Fback (Hpred i ) + Fback (Hi+1 )) . 2 The water depth Hi at the node i can be therefore computed using the water depth Hi+1 one node downstream. The equations are solved starting from the downstream boundary (node Nx + 1), where the water depth H(Nx + 1) is known: H(Nx + 1) = ξ(Nx + 1) − η(Nx + 1) = ξd − η(Nx + 1), with ξd the water surface elevation imposed by the user. • Download the Matlab program called backwater_equation.m as well as the function F_back. 59 (3.6) (3.7) • Identify the different steps in the computation of the water depth. Why is the loop starting at Nx + 1 and not at 1 as usual? • Run the script and analyse the output. • The normal water depth Hn and critical depth Hc are computed at the beginning of the script (for explanations about these depths, look at the part relative to the backwater equations in the powerpoint presentation Chapter 5 of Gary Parker’s e-book). Add on the figure two lines corresponding to H = Hc and H = Hn. Modify the legend. • Run the code for different values of ξd . Different behaviours can be observed depending on the value of Hd = ξd +η(N x+1). Make a representative figure for each type of observed behaviour (to add to your word document) and comment on your findings. Discuss in particular what is happening for values less than the critical depth. • Calculate now a second vector H2, which contains the values of water depth obtained without using a corrector step. Plot the results with and without using a corrector on the figure and compare them. This comparison should be performed for a fine spatial grid (the one given initially for instance) and a very large one. Do you understand why predictor-corrector schemes are used? 3.3 Morphodynamic model with flow computation The model studied in this section is based on the same formulation as the spreadsheet RTe-bookAdDegBW.xls from Gary Parker e-book. It allows to compute the variations of the bed level η in space and time for different upstream and downstream conditions (variables in red in Fig. 3.2). This resolution is performed in 3 steps, summarized in Fig. 3.3: - the first part computes the water depth and flow velocity, based on the backwater equations presented in Section 3.2, - the sediment transport is then computed using the flow characteristics determined in the previous step, 60 - finally, Exner equations are used to calculate the bed elevation from the sediment transport fluxes. • The folder full model contains the different m-files needed for this assignment. • A detailed description of the model is given in chapter 20 of the e-book (a word file, not the ppt).Download this document and use it to identify the different parts of the script model.m. Look through the different functions called in this script and try to understand what is going on. • Run the script. Why doesn’t the bed change? • Implement the sediment transport block in the model. Change the definition of qt and remove the command qt_u=0;, which sets the sediment transport at the upstream boundary to zero. Calculate the shear stress from the flow velocity and friction; not from the depth slope product. • Run the script again, analyse the output. Where is sediment deposited and why? • What is happening when you raise the sea water level while keeping the sediment input constant? Why? • What is happening if you decrease the sediment input? Why? • In the Rhine, both happened at the same time. What happens then? • Play with the model! 61 REMARK Your model might produce incorrect results for some values of grid size and time step. More precisely, if the hump on the bed migrates more than one grid step dx during one time step dt, the model becomes unstable. This gives the following stability condition (also known as the Courant-Friedrichs-Lewy (CFL) condition): c dt ≤ 1, dx where c is the characteristic celerity of a disturbance in bed level. In practice, it means that if you increase significantly dt (to decrease your computation time for instance) you may have to increase dx to keep the model stable. However, keep in mind that increasing dx decreases the spacial resolution of your model. 62 Figure 3.2: Initial condition and variable definitions 63 Figure 3.3: Organization of the morphodynamic model introduced in Section 3.3 64 Chapter 4 Creative Matlab exercise The last (creative) exercise gives you some freedom of subject and method and allows you to distinguish yourself (earn higher grades). Furthermore, you can perhaps use the resulting code in the delta project later in this course. The idea is that you conceive your own researchable (quantitative) question, do the research using matlab, and present it in a powerpoint to the group. Steps of the empirical research cycle, that you will go through and that you will include in the five minute powerpoint presentation, are: 1. a researchable question 2. hypothesis (based on literature, provide references) 3. method 4. results 65 5. discussion (in view of hypothesis) 6. conclusions 7. new research questions (at least one). Before you start on the research propose your research question to the supervisors early Monday afternoon (22 September), or send your research question by email to [email protected] on Monday evening 22 September the latest. Example research questions have been mentioned during class. The creative exercise will be presented using powerpoint in max 5 minutes per person (see schedule in study guide for date, time and location). Bring your powerpoint file on a USB stick. You may do the exercise on your own or in student pairs. If you work in pairs there should be two (related) research questions and investigations, and two powerpoints of 5 minutes each because each should present his/her own work. Indicate in the code (as comments) the contribution of each partner. The individual creative exercise is worth 20% and the other three exercises are worth 5% each of the final grade. Grades will be individual and will be based on the quality of all steps in the empirical cycle plus the presentation quality plus, very importantly, creativity! Grades for having all steps done well and the presentation itself can be troll (2), not yet good enough (4), just good enough (6), excellent (8) and outstanding (10) and the final grade is the numeric average of these, where steps are whimsically weighed by staff. The Matlab code is graded on the same scale separately on three criteria: does it work, does it do what it is supposed to do, and can we, experienced yet mortal and impatient lecturers, understand it. Have fun! 66
© Copyright 2024