"بسمه تعالی"
تمرین 1 – هیدرولیک پیشرفته
"ترسیم پروفیل سطح آب با روش های اویلر و گام به گام مستقیم توسط نرم افزار MATLAB"
مساله:
یک کانال ذوزنقه ای با مشخصات زیر، آ ب را با دبی 28 m^3⁄s حمل می کند. اگر کانال ها در انتها به یک افتادگی ناگهانی برسد، با استفاده از دو روش گام به گام مستقیم و روش عددی اویلر، ضمن نوشتن برنامه ی کامپیوتری، پروفیل سطح آب را رسم نمایید.
عرض کانال: b=6.1 m
شیب جانبی مقطع: z=2
شیب کف کانال: S_0=0.001
عدد مانینگ: n= 0.025
"روش گام به گام مستقیم"
این روش در کانال های منشوری با مقاطع منظم و یا نامنظم به کار می رود. و اساس آن بر معادله ی انرژی استوار است.
این روش به محاسبه ی فاصله از روی عمق می پردازد.
معادله ی انرژی بین دو مقطع (1) و (2) که به فاصله ی xΔ از هم قرار دارند:
H_1-h_f=H_2
y_1+α_1 (v_1^2)/2g+z_1-h_f=y_2+α_2 (v_2^2)/2g+z_2
y_1+α_1 (v_1^2)/2g+S_0 ∆x-S_f ∆x=y_2+α_2 (v_2^2)/2g
∆x=([y_2+(α_2 v_2^2)/2g]-[y_1+(α_1 v_1^2)/2g])/(S_0-S_f )=(E_2-E_1)/(S_0-S_f ) (1)
در رابطه (1)، S_(0 ) و S_f شیب طولی کانال و شیب خط انرژی می باشند که S_f مقدار میانگین دو شیب خط انرژی مربوطه می باشد و از رابطه ی مانینگ محاسبه می شود.
S_f=(S_f1+S_f2)/2=(S_f ) ̅
S_f=(Q^2 n^2)/(R^(4⁄3) A^2 )
در این مثال α=1 فرض شده است.
1-محاسبه ی عمق نرمال و عمق بحرانی:
به منظور تحلیل نیمرخ سطح آب، عمق نرمال و عمق بحرانی در کانال محاسبه می شوند. از آنجا که مقطع کانال ذوزنقه ای می باشد، برای محاسبه ی عمق نرمال یا از رابطه ی (2) و یا از نمودارهای مربوط به کانال ذوزنقه ای استفاده می شود. در اینجا از رابطه ی (2) استفاده شده است.
(Q^2 T)/(g A^3 )=1 (2)
T= b+2zy , A=(b+zy)y
با قرار دادن مقادیر مربوط به مساله در روابط فوق تنها مجهول مساله y_C بوده که با روش سعی و خطا قابل محاسبه است.
محاسبه ی عمق نرمال از رابطه ی مانینگ امکان پذیر است.
Q=1/n R^(2⁄3) S_0^(1⁄2) (3)
R=(b+zy)y/(b+2y√((1+z^2)))
معادله ی (1) با استفاده از روش سعی و خطا قابل محاسبه است. در اینجا هر دو معادله ی فوق با برنامه نویسی محاسبه شده اند:
% % Normal depth calculation
f=inline('((1/0.025)*((6.1+2.*y0).*y0)*sqrt(0.001).*(((6.1+2.*y0).*y0)/(6.1+2.*sqrt(5).*y0)).^(2/3))-28','y0'); % Q=(1/n)*R^(2/3)*S(1/2)
y0=fzero(f,[0.01,20]) % y0=1.8938
% % critical depth calculation
g=inline('((28^2).*(6.1+4.*yc))./(9.81.*((6.1+2.*yc).*yc).^3)-1','yc'); % ((Q^2*T)/(g*A^3))=1
yc=fzero(g,[0.01,10]) % yc=1.1322
مقادیر به دست آمده برای y_C و y_0 به ترتیب برابر با 1.1322 و 1.8938 می باشند. بنابراین:
y_0>y_C → جریان زیر بحرانی است
جریان زیر بحرانی بوده پس محاسبات از پایین دست آغاز می شود ونقطه ی شروع در همسایگی مقطع کنترل (y_C) خواهد بود. از مطالب گفته شده می توان نتیجه گرفت نوع پروفیل جریان M_2 است.
2- تعیین dy و محاسبه ی y و x در مقاطع مختلف
از آنجایی که اساس روش گام به گام مستقیم محاسبه ی فاصله از روی عمق می باشد. بنابراین با تعریف dy به عنوان گام های محاسبه و y_2=y_1+dy، می توان به ازای هر عمق، مساحت (A)، شعاع هیدرولیکی(R)، انرژی مخصوص(E)و شیب خط انرژی (S_f) را تعیین نموده و پس از مشخص نمودن مقدار ∆E و (S_f ) ̅ در حد فاصل دو عمق ∆x از رابطه ی (1) محاسبه می شود.
در نهایت با دست داشتن مختصات (X,Y) در هر نقطه پروفیل سطح آب ترسیم می شود.
متن برنامه-روش گام به گام مستقیم
clc;
clear all;
% % Normal depth calculation
f=inline('((1/0.025)*((6.1+2.*y0).*y0)*sqrt(0.001).*(((6.1+2.*y0).*y0)/(6.1+2.*sqrt(5).*y0)).^(2/3))-28','y0'); % Q=(1/n)*R^(2/3)*S(1/2)
y0=fzero(f,[0.01,20]) % y0=1.8938;
% % critical depth calculation
g=inline('((28^2).*(6.1+4.*yc))./(9.81.*((6.1+2.*yc).*yc).^3)-1','yc'); % ((Q^2*T)/(g*A^3))=1
yc=fzero(g,[0.01,10]) % yc=1.1322;
if y0<yc %if y0<yc then the flow is supercritical else the flow is subcritical
yinit=1.01*y0
else
yinit=1.01*yc
end
% % Problem information
Q=28; %Discharge (m^3/s)
n=0.025; %Manning Roughness Coefficient
g=9.81;
s0=0.001; % Bottom slope
N=100;
% yinit=1.01*yc % initial guess
h=(y0-yinit)/N % computational step
b=6.1;% Bottom Width (m)
z=2;%Side Slope
syms Y
A=b.*Y+z.*Y.^2
A1=subs(A,yinit,Y);
T=b+2*z.*Y;
T1=subs(T,yinit,Y)
P=b+2.*Y.*sqrt(1+z^2);
P1=subs(P,yinit,Y)
R=inline('A./P','A','P');
R1=R(A1,P1)
Sf=inline('0.025^2*28^2./(R.^(4/3).*A.^2)','R','A')
Sf1=Sf(R1,A1)
E=inline('Y+((28^2)./((2*9.81).*A.^2))','Y','A');
E1=E(yinit,A1);
Y=zeros(N+1,1);
X=zeros(N+1,1);
R=zeros(N+1,1);
A=zeros(N+1,1);
E=zeros(N+1,1);
Sf=zeros(N+1,1);
SSf=zeros(2,N);
MSf=zeros(N,1);
Y(1)=yinit;
X(1)=0;
A(1)=A1;
R(1)=R1;
E(1)=E1;
Sf(1)=Sf1;
hb(1)=0
hw(1)=Y(1)
for i=1:N
Y(i+1)=Y(i)+h;
A(i+1)=b.*Y(i+1)+z.*Y(i+1).^2;
T(i+1)=b+2*z.*Y(i+1);
P(i+1)=b+2.*Y(i+1).*sqrt(1+z^2);
R(i+1)=A(i+1)./P(i+1);
E(i+1)=Y(i+1)+((Q^2)./((2*9.81).*A(i+1).^2));
Sf(i+1)=n^2*Q^2./(R(i+1).^(4/3).*A(i+1).^2);
SSf(1,i)=Sf(i); SSf(2,i)=Sf(i+1);
MSf(i)=mean(SSf(:,i));
DeltaX(i)=(E(i+1)-E(i))/(s0-MSf(i));
X(i+1)=X(i)+DeltaX(i); %distance from beginning
hb(i+1)=X(i+1).*(-s0); % Bed elevation
hw(i+1)=hb(i+1)+Y(i+1); % Water surface elevation
end
figure(1);
hc=hb+yc; % the critical depth
h0=hb+y0; % the normal depth
plot (X,hw,'LineWidth',2);
hold on
plot(X,hb,'k-','LineWidth',2);
plot(X,hc,'–r');
plot(X,h0,'-.m' );
h=legend('Water surface','Canal bed','Critical depth','Normal depth');
text(-500,2.5,'M2','FontSize',12)
hold off
پروفیل سطح آب با روش گام به گام مستقیم
"روش اویلر"
در این روش، عمق جریان در فواصل طولی یکسان و یا در موضع خاصی از کانال محاسبه می شوند. بنابراین فاصله از نقطه شروع به عنوان داده ی معلوم مساله مطرح می شود.
dy/dx=f(y)
dy/dx=∆y/∆x=(y_2-y_1)/∆x
y_2=y_1+dy/dx ∆x (4)
در صورتی که dy/dx و y_1 معلوم باشند مقدار y_2 مشخص می شود. در شروع محاسبات با توجه به معلوم بودن مقدار y_1 (در اینجا y_C) ، مقدار dy/dx برابر ├ dy/dx┤|_(x=x_1 ) در نظر گرفته شده و در رابطه ی (4) قرار داده می شود. یعنی:
y_2=y_1+├ dy/dx┤|_(x=x_1 ) ∆x (5)
منظور از ├ dy/dx┤|_(x=x_1 ) محاسبه ی شیب در مختصات (x_1,y_1 ) یا در مقدار y_1 می باشد. سپس در صورتی که فرض شود تغییرات شیب در فاصله ی ∆x خطی است، مقدار dy/dx در رابطه ی (5) با میانگین گیری آن بین دو مقدار y_1 و y_2 (محاسبه شده) اصلاح و مقدار اصلاح شده در رابطه ی (4) به منظور تعیین y_2 جدید قرار داده می شود:
y_2=y_1+(├ dy/dx┤|_(x=x_1 )+├ dy/dx┤|_(x=x_2 ) ) ∆x/2 (6)
با توجه به اینکه معادله ی حاکم بر جریان متغیر تدریجی طبق رابطه ی (7) بیان می شود:
dy/dx=(S_0-S_f)/(1-〖Fr〗^2 )=f(y) (7)
لذا می توان روش ارائه شده را در حل معادله ی (7) به کار برد و در محاسبات را از مقطع کنترل (در اینجا y_C ) تا هر طول دلخواه انجام داد.
در اینجا محاسبات را برای طول 1000 متری و در فاصله های dx=1 m انجام می دهیم. سپس پروفیل رسم شده را با پروفیل های حاصل ازdx=0.01,0.1, 4 ,6 m مقایسه می کنیم.
متن برنامه-روش اویلری
clear all;
Q=28; %Discharge (m^3/s)
s0=0.001; % Bottom slope
n=0.025; %Manning Roughness Coefficient
z=2 %Side Slope
b=6.1; % Bottom Width (m)
% % Normal depth calculation
syms y0
f=inline('((1/0.025)*((6.1+2.*y0).*y0)*sqrt(0.001).*(((6.1+2.*y0).*y0)/(6.1+2.*sqrt(5).*y0)).^(2/3))-28','y0'); % Q=(1/n)*R^(2/3)*S(1/2)
y0=fzero(f,[0.01,10]) %y0= 1.8938;
% % critical depth calculation
syms yc
g=inline('((1*28^2).*(6.1+4.*yc))./(9.81*sqrt(1/(1+(0.001)^2)).*((6.1+2.*yc).*yc).^3)-1','yc'); % ((Q^2*T)/(g*A^3))=1
yc=fzero(g,[0.01,10]) %yc=1.1322;
if y0<yc %if y0<yc then the flow is supercritical else the flow is subcritical
yinit=1.01*y0
else
yinit=1.01*yc
end
for dx=[0.01 0.1 1 4 6];
dx=-dx;
L=1000;
N=abs(L/dx);
A=zeros(N,1);
T=zeros(N,1);
P=zeros(N,1);
R=zeros(N,1);
Sf=zeros(N,1);
f=zeros(N,1);
Yn=zeros(N,1);
fn=zeros(N,1);
Y=zeros(N+1,1);
hb=zeros(N+1,1);
hw=zeros(N+1,1);
X=zeros(N+1,1);
Y(1)=yinit; % initial guess
hb(1)=0;% The first bed elevation
hw(1)=Y(1);% The first water elevation
X(1)=0;
for i=1:N
A(i)=b.*Y(i)+z.*Y(i).^2; % Area m^2
T(i)=b+2*z.*Y(i); % Top Width
P(i)=b+2.*Y(i).*sqrt(1+z^2); % Wetted Perimeter
R(i)=A(i)./P(i); % Hydraulic radius
Sf(i)=n^2*Q^2./(R(i).^(4/3).*A(i).^2); % Energy slope
f(i)=(s0-Sf(i))./(1-Q^2.*T(i)./9.81./A(i).^3); % f=(s-sf/(1-fr^2)), fr=Q^2*T/g*A^3
Yn(i)=Y(i)+f(i).*dx; % first depth
A(i)=b.*Yn(i)+z.*Yn(i).^2;
T(i)=b+2*z.*Yn(i);
P(i)=b+2.*Yn(i).*sqrt(1+z^2);
R(i)=A(i)./P(i);
Sf(i)=n^2*Q^2./(R(i).^(4/3).*A(i).^2);
fn(i)=(s0-Sf(i))./(1-Q^2.*T(i)./9.81./A(i).^3);%Calculate new (s-sf/(1-fr^2)) by using Yn(i)
Y(i+1)=Y(i)+(fn(i)+f(i)).*dx/2; % modified depth
hb(i+1)=-dx.*(i).*s0; % Bed elevation
hw(i+1)=hb(i+1)+Y(i+1);% Water surface elevation
X(i+1)=dx*i; %the distance from begining point
end
figure(1);
hc=hb+yc; % the critical depth
h0=hb+y0; % the normal depth
plot (X,hw,'LineWidth',0.5);
hold on
plot(X,hb,'k-','LineWidth',2);
plot(X,hc,'–r');
plot(X,h0,'-.m' );
xlabel('distance')
ylabel('elevation')
h=legend('Water surface','Canal bed','Critical depth','Normal depth');
keyboard
end
hold off
پروفیل سطح آب با روش اویلری
در حالت dx=1
حالت dx=0.1
حالت های dx=0.01, 0.1, 1, 4,6
مقایسه ی پروفیل ها نشان می دهد که هر چقدر میزان dx کوچکتر باشد نتایج به مقدار واقعی نزدیک تر خواهد بود و با بزرگتر شدن شکل به دست آمده غیرمنطقی به نظر می رسد. پروفیل های حاصل از dx=0.01,0.1,1 انطباق خوبی دارند و با بزرگ شدن مقادیر غیر حقیقی می گردد.