BestTutorial2016Sediment DynamicsPARTA-Jo hohai university pdf

HOHAI UNIVERSITY REPORT ON SEDIMENT PROCESS MARCH 18, 2016 SUBMITTED TO: ASSOC. PROF. XI LI SUMITTED BY: ISHWAR JOSHI M2015076 School of Water Conservancy and Hydropower Engineering 1 Tutorial 1/2: The tutorial contains two parts: 1. To determine the settling velocity as a function of concentration in the combined flocculation and hindered settling range by fitting an empirical relation to experimental data on C(z,t) from a settling column. 2. To use settling velocity formula obtained in the first part to simulate the time-evolution of the concentration profile during settling. Theory: The empirical relation for settling velocity in the region of hindered settling and flocculation is given by WS = (2+2) Where, a, b, n, m are sediment dependent empirical constants. Fig: Settling Velocity variation with suspension concentration 2 Solution in Matlab: Settling Velocity and Concentration Profile Simulations Number of elevation for settling column concentration sampling, nele, is: nele = 8 Number of times for concentration sampling, ntime, is: ntime = 7 Input file name is: setm2.txt Estimated coefficient, m, is: m = 2 Estimated coefficient, n, is: n = 1.4 Estimate suspension concentration for maximum settling velocity, C2 (kg/m3), is: C2 = 6 Calculate coefficient, b, based on equation (4.18) page 4-28 is: b = 8.176622 Estimate maximum settling velocity, Ws2 (m/s), is: Ws2 = 0.02 Calculate coefficient, a, based on equation (4.17) page 4-26 is: a = 0.043003 Do you want to replot for y-axis range? Y/N Y: Y Minimum settling velocity for velocity-concentration plot: y_min = 1e-6 Maximum settling velocity for velocity-concentration plot: y_max = 1e-2 3 Do you want to estimation of the four parameters for settling velocity formula? Y/N Y: Y Estimated coefficient, a, is: a = 0.15 Estimated coefficient, b, is: b = 4 Estimated coefficient, a, is: m = 2 Estimated coefficient, n, is: n = 1.5 Limiting free settling concentration, C1 (kg/m3), is: C1 = 0.2 Do you want to simulate concentration profiles? Y/N Y: Y Time-step for calculation (seconds) is: time = 10 Fig: Initial Profile 4 Fig: Plot of Settling Velocity versus Sediment Concentration Fig: Time Evolution of Concentration Profile 5 The Source code in Matlab is shown in Appendix A. Tutorial 3: Runge Kutta Method: Theory: Runge Kutta is an explicit solution to sedimentation concentration in a basin with one entrance. The explicit 4th-order Runge-Kutta method is given by Ci+1 = Ci + (tF(t, Ci ) + 2k2 + 2 k3 + k4)/6 Where, i =1,2, .refers to the time step, k1 = tF(t, Ci ) , k2 =t F(t+t/2, Ci +k1/2 ) , k3 = t F(t+t/2, Ci +k2/2 ) , k4 = F(t+t/2, Ci +k3) and t is the finite time difference. Solving Runge Kutta Method using Matlab: Solution in Matlab: LGKT h=5 x0=1 y0=2 input is? X=2 Y=3 LGKT h=1 x0=1 y0=1 input is? X=1 Y=2 result is=2.000,0.3782234,0.4000000 6 The Source code in Matlab is shown in Appendix A. Tutorial 4: Numerical Simulation of Suspended Sediment Concentration and Depth of shoaling in a basin with one entrance. Theory: The Sediment concentration in basin with one entrance is given by the equation = F(t, C) = (0 )0 + The above equation can be solved by using explicit 4th-order Runge-Kutta method is given by Ci+1 = Ci + (tF(t, Ci ) + 2k2 + 2 k3 + k4)/6 Where, i =1,2, .refers to the time step, k1 = tF(t, Ci ) , k2 =t F(t+t/2, Ci +k1/2 ) , k3 = t F(t+t/2, Ci +k2/2 ) , k4 = F(t+t/2, Ci +k3) and t is the finite time difference. The settled sediment mass per unit area in the basin , is given by Sm(t) = 0 In Finite form Sm(t) = =1 The Shoaling thickness St(t) = () Solution in Mat Lab: Running the matlab file “sdt0.m”: Calculation of Sediment Rate in a Basin with Single Entrance Screen input Sediment concentration outside basin, ci (kg/m3) = ci = 0.1 Four parameters for the settling velocity formula, a, b, m, n a = 0.15 b = 4.0 m = 2.0 n = 1.5 Limit free settling concentration, c1 (kg/m3) = c1 = 0.2 Tidal-mean water depth in basin, h (m) = h = 5.0 Tidal amplitude, a0 (m) = a0 = 1.0 Tidal period, tp (hr) = tp = 12.5 7 Total simulation time, tt (day) = tt = 5 Output time step, ot (hr) = ot = 0.5 Calculation time step, dt (s) = dt = 60 Bed dry density, rhob (kg/m3) = rhob = 300 Do you want to replot for the figure with different of ranges of concentration and time? Y/N Y: Y Maximum time for the plot (day) = max_x1 = 5 Maximum suspended concentration for the plot = max_y1 = 0.05 This is the Runge - Kutta Order Four Method. Input the function F(t,y) in terms of t and y For example: y-t2+1 Input left and right endpoints on separate lines. A = 0 B = 5 Input the initial condition 0.012 Input a positive integer for the number of subintervals 200 Choice of output method 1. Output to screen 2. Output to text file Please enter 1 or 2 2 Output file output.txt created successfully 8 Do you want to replot for the figure with different of ranges of concentration and time? Y/N Y: Y Maximum time for the plot (day) = max_x2 = 5 Maximum shoaling thickness for the plot = max_y2 = 5 9 The Source code in Matlab is shown in Appendix A. Modelling Framework done in Surfer 10 APPENDIX A Tutorial 1/2: clc clear all format short disp(Settling Velocity and Concentration Profile Simulations) disp(Number of elevation for settling column concentration sampling, nele, is:) nele = input(nele = ); disp(Number of times for concentration sampling, ntime, is:) ntime = input(ntime = ); disp(Input file name is: setm2.txt) % Read file setm2.txt fid = fopen(setm2.txt,r); filetemp = fscanf(fid,%f,8 9); filedata = filetemp.; fclose(fid); %* First estimation for settling velocity formula * disp(Estimated coefficient, m, is:) m = input(m = ); disp(Estimated coefficient, n, is:) n = input(n = ); disp(Estimate suspension concentration for maximum settling velocity, C2 (kg/m3), is:) C2 = input(C2 = ); disp(Calculate coefficient, b, based on equation (4.18) page 4-28 is:) b=C2*(2*m/n-1)0.5); str = sprintf(b = %f,b); disp(str) disp(Estimate maximum settling velocity, Ws2 (m/s), is:) Ws2 = input(Ws2 = ); disp(Calculate coefficient, a, based on equation (4.17) page 4-26 is:) temp_b=(b(n-m)*(2*m/n-1)(m-0.5*n)/(2*m/n-1)0.5); a=(Ws2*1/temp_b); str = sprintf(a = %f,a); disp(str) %Replotting option for y-axis range - Insert Yes or No. If Yes: reply = input(Do you want to replot for y-axis range? Y/N Y: , s); if isempty(reply)|strcmp(reply,Y) reply = Y; disp(Minimum settling velocity for velocity-concentration plot:) ymin = input(y_min = ); disp(Maximum settling velocity for velocity-concentration plot:) ymax = input(y_max = ); for i=2:9 11 for j=2:8 cc(i,j)=filedata(i,j); ww(i,j)=(a*(cc(i,j).n)/(cc(i,j).2+b.2).m; end end loglog(cc,ww,o) xlabel(Sediment Concentration (g/l); ylabel(Settling Velocity (m/s); ylim(ymin ymax) grid off end for i=2:9 for j=2:8 cc(i,j)=filedata(i,j); ww(i,j)=(a*(cc(i,j).n)/(cc(i,j).2+b.2).m; end end loglog(cc,ww,o) title(Fig.1 Experimental data,fontsize,10,fontweight,bold); xlabel(Sediment Concentration (g/l); ylabel(Settling Velocity (m/s); grid off %Option for further estimation of the four parmeters for settling velocity formula - Insert Yes or No. If Yes: reply = input(Do you want to estimation of the four parameters for settling velocity formula? Y/N Y: , s); if isempty(reply)|strcmp(reply,Y) reply = Y; clear a b m n; disp(Estimated coefficient, a, is:) a = input(a = ); disp(Estimated coefficient, b, is:) b = input(b = ); disp(Estimated coefficient, a, is:) m = input(m = ); disp(Estimated coefficient, n, is:) n = input(n = ); disp(Limiting free settling concentration, C1 (kg/m3), is:) C1 = input(C1 = ); c0=C1-0.1; for i=1:1:1000 c(i)=c0+0.1*i; Ws(i)=(a*(c(i).n)/(c(i).2+b2).m); end loglog(c,Ws) title(Fig.2 Calculated of settling velocity and experimentaldata,fontsize,10,fontweight,bold); 12 xlabel(Sediment Concentration (g/l); ylabel(Settling Velocity (m/s); xlim(10-1 102) ylim(10-6 10-2) grid off hold on for i=2:9 for j=2:8 cc(i,j)=filedata(i,j); ww(i,j)=(a*(cc(i,j).n)/(cc(i,j).2+b.2).m; end end loglog(cc,ww,o) hold off end %Option to simulate concentration profiles - Insert Yes or No. If Yes: reply = input(Do you want to simulate concentration profiles? Y/N Y: , s); if isempty(reply)|strcmp(reply,Y) reply = Y; disp(Time-step for calculation (seconds) is:) time = input(time = ); aa = importdata(setm2.txt); bb = aa(2:9,1); cx = aa(2:9,2:8); x = logspace(-2,3); semilogx(cx,bb,-); hold on semilogx(cx,bb,o); hold off axis(1e-2,1e3,0,2); title(Fig.3 Time-evolution of the concentration profile,fontsize,10,fontweight,bold); xlabel(Sediment Concentration (g/l); ylabel(Elevation (m); end Tutorial 3: function = LGKT(h,x0,y0,X,Y) format long h=input(h=); x0=input(x0=); y0=input(y0=); disp(input is?); X=input(X=);Y=input(Y=); n=round(Y-X)/h); i=1;x1=0;k1=0;k2=0;k3=0;k4=0; 13 for i=1:1:n x1=x0+h; k1=-x0*y02;%k1=y0-2*x0/y0;% k2=(-(x0+h/2)*(y0+h/2*k1)2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);% y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);% x0=x1;y0=y1; y=2/(1+x02);%y=sqrt(1+2*x0);% fprintf(result is=%.3f,%.7f,%.7fn,x1,y1,y); end Tutorial 4: clc clear all format short disp(Calculation of Sediment Rate in a Basin with Single Entrance) %Option for input data - select screen input of file input disp(Screen input) disp(Sediment concentration outside basin, ci (kg/m3) =) ci = input(ci = ); disp(Four parameters for the settling velocity formula, a, b, m, n) aa = input(a = ); bb = input(b = ); m = input(m = ); n = input(n = ); disp(Limit free settling concentration, c1 (kg/m3) =;) c1 = input(c1 = ); disp(Tidal-mean water depth in basin, h (m) =) hh = input(h = ); 15 disp(Tidal amplitude, a0 (m) =;) a0 = input(a0 = ); disp(Tidal period, tp (hr) =) tp = input(tp = ); disp(Total simulation time, tt (day) =) tt = input(tt = ); disp(Output time step, ot (hr) =) ot = input(ot = ); disp(Calculation time step, dt (s) =) dt = input(dt = ); disp(Bed dry density, rhob (kg/m3) =) rhob = input(rhob = ); %Option for replotting the figure with different of ranges of concentration and time - Insert Yes or No. If Yes: 14 reply = input(Do you want to replot for the figure with different of ranges of concentration and time? Y/N Y: , s); if isempty(reply)|strcmp(reply,Y) reply = Y; disp(Maximum time for the plot (day) =) max_x1 = input(max_x1 = ); disp(Maximum suspended concentration for the plot =) max_y1 = input(max_y1 = ); end % RUNGE-KUTTA (ORDER 4) % % TO APPROXIMATE THE SOLUTION TO THE INITIAL VALUE PROBLEM: % Y = F(T,Y), A=T= B fprintf(1,Left endpoint must be less than right endpointn); else OK = TRUE; end; end; fprintf(1,Input the initial conditionn); ALPHA = input( ); OK = FALSE; while OK = FALSE fprintf(1,Input a positive integer for the number of subintervalsn); N = input( ); if N = 0 15 fprintf(1,Number must be a positive integern); else OK = TRUE; end; end; if OK = TRUE fprintf(1,Choice of output method:n); fprintf(1,1. Output to screenn); fprintf(1,2. Output to text filen); fprintf(1,Please enter 1 or 2n); FLAG = input( ); if FLAG = 2 NAME = output.txt; OUP = fopen(NAME,wt); else OUP = 1; end; % STEP 1 H = (B*24-A)/N; T = A; W = ALPHA; fprintf(OUP, %5.3f %11.7fn, T, W); % STEP 2 for I = 1:N % STEP 3 % use K1, K2, K3, K4 for K(1), K(2), K(3), K(4) RESP. K1 = H*F(T,W); K2 = H*F(T+H/2.0, W+K1/2.0); K3 = H*F(T+H/2.0, W+K2/2.0); K4 = H*F(T+H,W+K3); % STEP 4 % compute W(I) W = W+(K1+2.0*(K2+K3)+K4)/6.0; % compute T(I) T = A+I*