Add files via upload

Stoker_solution.m

@@ -0,0 +1,167 @@
+%**************************************************************************
+%       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({"Analytical solution (Stoker, 1957) 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
+
+