function = explicit(mx) %% Problem Number: 49 L = 45; tmax = 80; U = 0.8; D = 1.2; Xs = 4; Cs = 10; XFlux = 36; % % % % % % % m s m/s m2/s m kg/m m % Calculating dx and dt based of mx input dx = L/(mx-1); dt = ((dx)^2)/(2*10*D); % Stability coefficients Stab1 = (U*dx)/(2*D); Stab2 = (D*dt)/((dx)^2); % Stability conditions if Stab1 > 1 || Stab2 > 1/2 error('dx and dt do not satisfy the stability conditions outlined for the explicit scheme') end % Initialising values and vectors mt = nearest((tmax/dt)+1); t = linspace(0, tmax, mt); B = zeros(mx, 1); % Initial conditions B(nearest((Xs/dx)+1)) = Cs; % Plotting vectors Yvector = 0 : dx : L; v = linspace(0,Cs,200); %% Explicit Scheme % Numerical representation Ex1 = (D*dt)/(dx^2) + (U*dt)/(2*dx); Ex2 = ((-1*2*D*dt)/(dx^2)) + 1; Ex3 = (D*dt)/(dx^2) - (U*dt)/(2*dx); % Boundary conditions ExB1 = (U*dt)/dx; ExB2 = ((-1*U*dt)/dx) + 1; % Concentration matrix ExplicitConc = zeros(mx,length(t)); ExplicitConc(:,1) = B; % For loop to fill in concentrations for j = 2:size(ExplicitConc,2) for i = 1:size(ExplicitConc,1) if i == 1 ExplicitConc(i,j) = 0; elseif i == mx ExplicitConc(i,j) = ExB1*ExplicitConc(i-1,j-1) + ExB2*ExplicitConc(i,j-1); else ExplicitConc(i,j) = Ex1*ExplicitConc(i-1,j-1) + Ex2*ExplicitConc(i,j-1) + Ex3*ExplicitConc(i+1,j-1); end end end % Extracting Xflux and last rows for plotting XFluxExplicit = ExplicitConc(nearest((XFlux/dx)+1),:)*U; mxExplicit = ExplicitConc(mx,:)*U; %% Figures % Explicit Plots % 2D contour plot of explicit concentrations subplot(1,3,1) contour(t, Yvector, ExplicitConc,v) title('Figure 1 - Explicit Concentration Contour Plot') xlabel('Time (s)') ylabel('Distance Downstream (m)') % Line plot of explicit flux at Xflux distance subplot(1,3,2) plot(t, XFluxExplicit) title('Figure 2 - Explicit Concentration Flux at 36m') xlabel('Time (s)') ylabel('Flux (kg/sm^2)') % Line plot of explicit flux at final distance subplot(1,3,3) plot(t, mxExplicit) title('Figure 3 - Explicit Concentration Flux at L') xlabel('Time (s)') ylabel('Flux (kg/sm^2)') end Contents 1 Objectives 3 2 Governing equation 3 3 Method 3 3.1 3.2 3.3 3.4 3.5 Discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A \u001cnite di\u001berence representation of the governing equation . . . . . . . . . . . . 3.2.1 Explicit \u001cnite di\u001berence scheme . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Implicit \u001cnite di\u001berence scheme . . . . . . . . . . . . . . . . . . . . . . . Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Boundary conditions in the explicit scheme . . . . . . . . . . . . . . . . . 3.3.2 Boundary conditions in the implicit scheme . . . . . . . . . . . . . . . . Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A solution method for the numerical representation . . . . . . . . . . . . . . . . 3.5.1 Advancing the solution by one time step from a known concentration C n to the new unknown concentration C n+1 . . . . . . . . . . . . . . . . . . 3.5.2 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3 Computing the concentrations at all the time levels from C 0 (initial condition) to the unknown C \u001cnal . . . . . . . . . . . . . . . . . . . . . . . . . 4 4 4 5 5 5 5 6 6 6 7 8 4 Presenting the results 8 5 Some hints for undertaking the assignment 8 6 Assignment presentation 9 7 Marking criteria 9 7.1 7.2 7.3 7.4 7.5 7.6 7.7 Report \u001cgures . . . . . . . . . . . Introduction . . . . . . . . . . . . Methodology [3 marks] . . . . . . Results and discussion . . . . . . Conclusion and recommendations References . . . . . . . . . . . . . Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 . 10 . 10 . 10 . 10 . 11 . 11 1 Objectives This assignment has been designed with three speci\u001cc goals: To give you practical experience implementing most of the numerical techniques that you will be learning about this semester. To instruct you in the writing, testing and use of a reasonably complex computer program. To give you practice in applying numerical methods to practical engineering systems. The speci\u001cc task that needs to be completed requires you to write your own computer program that is capable of predicting the environmental impact of sewage discharges into a river. A conceptual diagram of the 1D spatial domain relevant to the problem is shown in Figure 1 (along with the 2 boundary conditions). A pollutant is discharged into the river at the initial C=0 U 2C =0 x2 i = mx i=1 0 i1 i i+1 L x Figure 1: Schematic of the Problem time (t = 0) at a location Xs to be speci\u001ced. The pollutant is advected by the river which has a \u001dow velocity U in the x-direction and also di\u001bused in the x direction at a rate proportional to the di\u001busivity coe\u001ecient D. The banks of the river are considered impermeable. For the purposes of this assignment, the river will be assumed to be one dimensional of length L. The river has a constant velocity U and a di\u001busivity coe\u001ecient D. 2 Governing equation This problem can be represented by the advection-di\u001busion equation for the sewage concentration C(x, t): C C 2C +U =D 2 t x x (1) This governing equation is e\u001bectively a statement of mass conservation of the e\u001fuent. The two processes represented by this equation are advection of the pollutant by the river \u001dow and spreading or dispersion of the pollutant by turbulence & velocity gradients (even though we will adopt a constant \u001dow velocity U throughout). 3 Method As we shall see during the course of this semester, every numerical model must include the following components: 3 3.1 Discretisation The region of interest (the so-called domain or solution domain) must be discretised \u0015 that is broken down into a set of representative points. Finite di\u001berences require that these points be located on line with spacings between mesh points of x. The x direction will have mx discrete points, numbered from 1 to mx. Time will also be discretised into equally spaced timesteps t. 3.2 A \u001cnite di\u001berence representation of the governing equation We will adopt \u001cnite di\u001berence expressions for the partial derivatives that are accurate to second order accuracy in space and \u001crst order in time. The \u001cnite di\u001berence form of each partial derivative can be written, using n as the time index, and i as the space index. By convention (and to help readability), the time index n is written as a superscript n and the space index i as subscript i. In the equations which follow, you will see that we have not yet speci\u001ced the time levels of the concentrations in the spatial derivatives. C C n+1 C n \u0002 + O(t) t t C Ci+1 Ci1 \u0002 + O(x2 ) x 2x 2 C Ci+1 2Ci + Ci1 \u0002 + O(x2 ) x2 x2 (2) (3) (4) In this assignment, you are required to implement two (2) versions of the \u001cnite di\u001berences equations: one called explicit, and the other one called implicit. 3.2.1 Explicit \u001cnite di\u001berence scheme In the explicit \u001cnite di\u001berence model, we choose to use the known concentrations at the old time level n in the spatial derivatives in eq (3) and (4). Equation (1) becomes: n n C n Ci1 C n 2Cin + Ci1 Cin+1 Cin + U i+1 =D i+1 t 2x x2 (5) All concentrations Cin are known at the old time level n and eq (5) is reorganised to compute Cin+1 , moving the unknown (at the new level n + 1) to the left hand side, and what is known (at time level n) to the right hand side: Cin+1 \u0002 Cin+1 = =Cin n n n n 2Cin + Ci1 Ci+1 Ci1 Ci+1 U t + Dt 2x x2 Dt U t + x2 2x \u0003 \u0002 n Ci1 + Dt U t x2 2x \u0003 \u0002 n Ci+1 \u0003 Dt + 2 2 + 1 Cin x How we can incorporate the boundary conditions is considered in Section 3.3. 4 (6) (7) 3.2.2 Implicit \u001cnite di\u001berence scheme In the implicit \u001cnite di\u001berence model, we now choose to use the unknown at the new time level n+1 in the spatial derivative in equations (3) and (4). Equation (1) becomes: n+1 n+1 C n+1 Ci1 C n+1 2Cin+1 + Ci1 Cin+1 Cin + U i+1 =D i+1 t 2x x2 (8) Cin are known at the old time level n and eq (8) is reorganised to compute n+1 Ci , moving the unknowns (at the time level n + 1) to the left hand side and what is known All concentrations (at the time level n) on the right hand side: Cin+1 \u0002 n+1 n+1 n+1 n+1 Ci1 2Cin+1 + Ci1 Ci+1 Ci+1 + U t =Cin Dt 2 2x x Dt U t 2 x 2x \u0003 \u0002 n+1 Ci1 Dt U t + 2 + x 2x \u0003 (9) \u0002 n+1 Ci+1 \u0003 Dt + 2 2 + 1 Cin+1 =Cin x (10) 3.3 Boundary conditions The domain must be bounded and at each boundary there are potential e\u001bects of the boundary on the interior concentrations in the domain. The \u001cnite di\u001berence mesh must be constructed so that all boundary points coincide with the boundary of the mesh. An equation (i.e boundary condition) must be speci\u001ced at every boundary point. 3.3.1 Boundary conditions in the explicit scheme At the western in\u001dow boundary, the in\u001dow concentration Ci at these points is 0 at all times. C1 =0 (11) At the eastern out\u001dow boundary, the boundary condition is complex. As a close approximation, we will assume the\u0003second spatial derivative of the concentration is zero normal \u0002 to the boundary 2C =0 x2 and combine this (using eq (4)) with the governing equation (eq (7)) to yield: \u0002 n+1 Cm x 3.3.2 = U t x \u0003 \u0002 n Cm x 1 \u0003 U t n + + 1 Cm x x (12) Boundary conditions in the implicit scheme At the western in\u001dow boundary, the in\u001dow concentration C1 =0 5 Ci is 0 at all times. (13) At the eastern out\u001dow boundary, the boundary condition is complex. As a close approximation, we will assume the\u0003second spatial derivative of the concentration is zero normal \u0002 2 to the boundary xC2 = 0 and combine this (using eq (4)) with the governing equation (eq (10)) to yield: \u0002 U t x \u0003 \u0002 n+1 Cm x 1 + \u0003 U t n+1 n + 1 Cm = Cm x x x (14) 3.4 Initial conditions The river is considered unpolluted at the initial time t = 0, so Ci0 =0 except (15) at one location xs where the sewage is released by accident: C (t = 0, x = xs ) =S0 (16) The discrete location of this accidental release is located at the discrete point which is the closest to (xs) 3.5 A solution method for the numerical representation The solution method to be adopted here is to write a Matlab program which assembles the full set of linear equations that incorporates an equation for each boundary and each internal point. (Note that if the in\u001dow and out\u001dow boundaries are used, you can reduce the number of equations). 3.5.1 Advancing the solution by one time step from a known concentration C n to the new unknown concentration C n+1 An explicit \u001cnite di\u001berence scheme is by de\u001cnition "explicit". When all the equations are combined (eq (7), (12) and (11)) the expression for Cin+1 depends only on all Cin , which are all known. The solution process requires the program to sweep through every node of the mesh and apply the previous relevant equations and boundary conditions. Explicit scheme Implicit scheme In this case, solution of the unknown concentrations at time level requires the solution of a set of simultaneous equations. A linear problem can always be written as a matrix-vector product of the kind: AC =B (n + 1) (17) where A is a coe\u001ecient matrix, C the vector of unknown concentrations at the new time level, and B the known RHS vector. In the problem we want to solve in this assignment, the domain is a line domain, containing mx points in the x direction. Because the sewage concentration is to be computed at each discrete 6 point i, the number of unknowns is mx . They need to be organised into a vector of unknown C which will contain mx rows. To form the product AC , the matrix A will need to have mx columns. And \u001cnally, to be able to solve the system to obtain a unique solution, we need as many equations as unknowns, so the matrix A will have mx rows and B will also have mx rows. A is consequently a square matrix of mx rows and mx columns. To \u001cnd a solution, the matrix equation de\u001cned in eq (17) is to be solved where A is the coe\u001ecient matrix as seen above, C is our (unknown) column vector of the sewage concentration values at each grid point, and B is a column vector containing the values from the previous concentrations corresponding to each grid point at the time level n. Speci\u001ccally, row i of A will contain the concentration coe\u001ecients from the corresponding equation for the point xi = ix. Almost all of the values in row i will be zero except for each column that corresponds to a point included in our \u001cnite di\u001berence equation. The process then becomes one of forming the (large) matrix in Matlab and then asking Matlab to compute the solution value at each mesh point. For the purposes of this assignment, you just compute A1 B to obtain C . 3.5.2 Example As an example, let us assume that we have U = 1, D = 0.5, L = 100, mx = 6 and t = 1. Therefore x = 20. A will be a 6 by 6 matrix. C and B will both have 6 rows. The \u001crst row in A will correspond to the point i = 1, which is on the western end. In matrixvector equation, the \u001crst line will be of the form 1 C1 = 0. All entries in the \u001crst row of A will be zero, except for column 1 which will contain the value 1. The corresponding value in row 1 of B will be zero. If the \u001crst time step is computed, from t = 0 to t = t, then the point i = 6 is the eastern end. Therefore, all entries in the last row of A will be zero except for columns 5 and 6 which will contain the values 0.05 and 1.05. The corresponding value in the 6th row of B will be 0. The point i = 3 is buried in the interior of the mesh. Therefore, all entries in the 3rd row of A will be zero except for columns 2, 3, 4, which will contain the values 0.0263, 1.0025, 0.0238 respectively. The corresponding value in the 3rd row of B will be 1. 7 3.5.3 Computing the concentrations at all the time levels from C 0 (initial condition) to the unknown C \u001cnal There are stability conditions associated with these numerical schemes on both t and x. explicit scheme: and U x 1 2D 1 Dt 2 x 2 (18) (19) Remember, errors are O(t, x2) implicit scheme: unconditionally stable, and remember, errors are proportional to (t, x2 ). A time stepping loop will need to be written in your code to solve consecutively each new set of unknown concentrations at each time level until the \u001cnal time t = t\u001cnal, using at each stage, the method detailed in the Section 3.5.1 4 Presenting the results Once your scripts have been written, you should consider presenting your results showing: 4-6 line plots of the concentration in time between t0 and t\u001cnal C.U at a position to compute the sewage concentration \u001dux per unit area (kg.s1 .m2 ) as a function of time. \u0004You can also study the mass conservation in time of the sewage by western end computing Mass(t) = eastern end C(t, x)dx and plot it. use your imagination! 5 Some hints for undertaking the assignment 1. Completing this assignment will require solid work over the next weeks. Leaving it until the last moment will result in disaster. 2. The problem can be broken down into the following components: sorting out the geometry. for the implicit model, assembling the coe\u001ecient matrix and RHS vector. inverting the coe\u001ecient matrix & multiplying it with the RHS to obtain the concentration values at each time step. plotting the initial concentration values. showing the time evolution by plotting concentration values at several other times. brie\u001dy discuss the advantages ans disadvantages of \u001cnite di\u001berence explicit and implicit models. space-time contour or 3D plots of the concentration. 8 3. The beauty of Matlab is that each step can be done interactively and checked without actually doing any programming. The program just follows the set of steps that you take to solve the problem interactively. 4. Each one of the steps above can be checked out by using dummy data to ensure that you understand what each command does. 6 Assignment presentation For the assignment presentation, you must prepare a brief engineering report including: 1. Introduction 2. Numerical methods 3. Results and Discussion 4. Conclusions 5. References Matlab includes some neat plotting packages including contouring which should give you a neat colour contour diagram of your concentration values within the solution domain. Calls to these packages should be made from your program to make the results \u001cgures in your report pop up automatically on the screen. You should suitably illustrate your report. You are also required to e-mail a copy of your m \u001cle to x.barthelemy@unsw.edu.au and this code will be used to verify that it produces the outputs and for automatic comparison with other student programs to ensure that it is your own work. 7 Marking criteria The assignments will be marked in 6 categories: \u001cgures, introduction, methodology, results, conclusions, and references. 7.1 Report \u001cgures Sketch of the problem [2 marks] Conventional plots of concentration and \u001duxes [2 marks] Innovative plots to make results particularly clear [1 mark] De\u001cnition sketches can be very simple to create, and they give the reader an instant appreciation of the problem. There is already a sketch provided in the assignment problem sheet. Marks will be awarded for the plots regardless of their accuracy. Marks will be mostly awarded to 1D plots or 2D space-time, provided they actually made the results clearer. 9 7.2 Introduction Clear speci\u001ccation of the problem [2 marks] Marks will be awarded for introductions that actually mention the speci\u001ccs of the problem, e.g. contamination problems. 7.3 Methodology [3 marks] Discretisation Boundary conditions Numerical representation Solution method Marks will be awarded to students who explained the general procedure of applying the numerical solution method, rather than simply describing matrix operations in Matlab . 7.4 Results and discussion Correct results [7 marks] Marks will be awarded for correctly plotting contours showing the distribution of concentrations. Discussion of maximum resolution achievable [2 marks] Marks will be awarded for discussing the computational limitations of their program, e.g. how \u001cne could the mesh be made before running out of memory? Discussion on the signi\u001ccance of the results [4 marks] Students will be awarded marks for identifying the consequences of releasing sewage in the environment. Comments of the results and their signi\u001ccance are more important than precise location and/or value. 7.5 Conclusion and recommendations Resolution adequacy [2 marks] Marks will be awarded for discussion on mesh resolution, and the trade-o\u001b between accuracy and computation time. Sensitivity [2 marks] Marks will be awarded for computing the problem using di\u001berent mesh resolutions and comparing the results. Other sensible remarks or recommendations [1 mark] Marks will be awarded for insightful comments describing the practical implications of the problem. 10 7.6 References References other than supplied reading material and websites [2 marks] Marks will be only awarded when students referenced at least one book or journal article. \u0015 Journals (for details on speci\u001cc research) \u0015 Textbooks (for background information) 7.7 Comments Some students will not manage to build a program successfully. Figures which are produced with another student's program will not be awarded any marks. This is a challenging assignment, but take note of the clear warning that it will take a signi\u001ccant amount of time to complete. 11