From d3de39a06b2c8b86cd41427daabd60ec0e62a39a Mon Sep 17 00:00:00 2001 From: Payam Sarkhosh <91587210+psarkhosh@users.noreply.github.com> Date: Mon, 25 Oct 2021 21:29:40 -0600 Subject: [PATCH] Delete Stoker_solution.m --- Stoker_solution.m | 167 ---------------------------------------------- 1 file changed, 167 deletions(-) delete mode 100644 Stoker_solution.m diff --git a/Stoker_solution.m b/Stoker_solution.m deleted file mode 100644 index 58d6a33..0000000 --- a/Stoker_solution.m +++ /dev/null @@ -1,167 +0,0 @@ -%************************************************************************** -% Stoker's (1957) analytical solution for ideal dam-break -% code generated by: Payam Sarkhosh -% Research assistant at Prof. Yee-Chung Jin's Lab -% University of Regina, Saskatchewan, Canada -% Fall 2021 -%************************************************************************** -clc -clear all -clf - -%******************************* inputs *********************************** -disp('*******************************************************************') -disp('******* Please enter the following inputs, and press ENTER ********') -disp('*******************************************************************') -disp(' ') -h0 = input(' initial upstream water depth: h0 (m) >> '); -hd = input(' initial downstream water depth: hd (m) >> '); -Lc = input(' channel length: Lc (m) >> '); -Lr = input(' reservoir length: Lr (m) >> '); -T = input(' total simulation time T (s) >> '); - -n=2000; %Number of spaceintervals -m=T*200; %Number of time intervals - -%************************************************************************** -dx=Lc/n; %space step -dt=T/m; %Time step -g=9.81; %Gravitational acceleration -x=zeros(n,1); %X vector -u=zeros(n,1); %Flow velocity vector -h=zeros(n,1); %Flow depth vector - -x_UBC=-Lr; -x_DBC=Lc-Lr; -x0=0; - -%************************ Plotting initial condition plot ***************** -x1=-Lr; -x2=0; -xDam2=linspace(x1,x2,n); -Dam2=h0+1e-50*xDam2; - -x1=x0+1e-10; -xDam1=linspace(x0,x1,n); -Dam1=h0*(xDam1-x0)/1e-10; - -%*************************** Mesh generation ****************************** -x(1)=x_UBC; -for i=1:n-1 - x(i+1)=x_UBC+i*dx; - if abs(x(i)-x0)<0.5*dx - i_0=i; - end -end -x(n)=Lc-Lr; -h(1)=h0; - -%****************** defenition of constant values************************** -x2_end=-1e-20; -hA=1e+20; - -Cup=(g*h0)^0.5; % wave speed at upstreamstream -Cdown=(g*hd)^0.5; % wave speed at downstream - -for k=1:m - Time=k*dt; - - %**************** Newton-Raphson iteration *************************** - f_CB=1; df_CB=1; CB=10*h0; - - while abs(f_CB/CB)>1e-10 - f_CB= CB*hd - hd*(((8*CB^2)/Cdown^2 + 1)^(1/2) - 1)*(CB/2 - Cup... - +((g*hd*(((8*CB^2)/Cdown^2 + 1)^(1/2) - 1))/2)^(1/2)) ; - df_CB= hd -hd*((2*CB*g*hd)/(Cdown^2*((8*CB^2)/Cdown^2 + 1)^(1/2)... - *((g*hd*(((8*CB^2)/Cdown^2 + 1)^(1/2) - 1))/2)^(1/2)) + 1/2)... - *(((8*CB^2)/Cdown^2 + 1)^(1/2) - 1) - (8*CB*hd*(CB/2 - Cup... - + ((g*hd*(((8*CB^2)/Cdown^2 + 1)^(1/2) - 1))/2)^(1/2)))/... - (Cdown^2*((8*CB^2)/Cdown^2 + 1)^(1/2)) ; - CB=CB-f_CB/df_CB; - end - %*************** Newton-Raphson iteration end ************************* - - hA=0.5*hd*((1+8*CB^2.0/Cdown^2.0)^0.5-1); - if hd==0 - CB=0; hA=0; - end - - X2_start=(2*(g*h0)^0.5-3*(g*hA)^0.5)*Time; - X2_end=CB*Time; - uA=2*Cup-2*(g*hA)^0.5; - - for i=2:n - - h(i)=(2*(g*h0)^0.5-x(i)/Time)^2.0/9/g; - u(i)=2*(x(i)/Time+(g*h0)^0.5)/3; - - h(1)=h(2); - u(1)=u(2); - - %****************************************************************** - if h(i)>=h0 - i_A=i; - h(i)=h0; - u(i)=0; - end - if hA==0 && h(i)>h(i-1) - h(i)=0; - u(i)=0; - end - if hA>0 - if x(i)<=X2_end && h(i)<=hA - i_B=i; - h(i)=hA; - u(i)=uA; - elseif x(i)>X2_end - h(i)=hd; - u(i)=0; - end - end - end - - if (rem(Time,1/m)==0) - subplot(2,1,1) - plot(xDam2,Dam2,'--k','LineWidth',1), hold on - plot(xDam1,Dam1,'--k','LineWidth',1) - plot(x,h,'b','LineWidth',3) - xlim([x_UBC x_DBC]) - ylim([0 1.1*h0]) - - y_label=ylabel('water depth (m)'); - set(y_label,'position',get(y_label,'position')-[0.2 0 0]); - set(gca,'FontSize',14) - - Time=Time+0.0001; - title({"Stoker's (1957) solution for ideal dam-break problem" - ['t = ',num2str(Time,'%.2f'),' s']},'FontSize',15) - Time=Time-0.0001; - hold off - - subplot(2,1,2) - brown = [0.5, 0, 0]; - plot(x,u,'Color',brown,'LineWidth',3) - xlim([x_UBC x_DBC]) - ylim([0 (g*h0)^0.5*2.2]) - x_label=xlabel('x (m)'); - set(x_label,'position',get(x_label,'position')+[0.15 -0.1 0]); - y_label=ylabel('velocity (m/s)'); - set(y_label,'position',get(y_label,'position')-[0.2 0 0]); - set(gca,'FontSize',14) - - fig=figure(1); - - if(Time == T) - disp(' ') - disp(['******* Outputs at t = ',num2str(Time),' s **********']) - T2 = table(x,h,u); - format short - disp(T2); - disp(['******* Outputs at t = ',num2str(Time),' s **********']) - end - format long - end - hold off -end - -