Ứng dụng chương trình Matlab trong việc giải một số bài toán tính đường mặt nước dòng chảy trong kênh hở

AR
Qn
S ==⇒ (1) CT 2 
Từ ph−ơng trình vi phân cơ bản của dòng chảy ổn định, thay đổi dần không áp: 
3
2
23
4
H
22
0
2
3
4
H
22
0
2
0
gA
BQ
1
AR
Qn
S
gA
BV
1
R
Vn
S
Fr1
SS
dx
dy
−
−
=
−
−
=
−
−= (2) 
trong đó: S0 - độ dốc đáy kênh; 
n - hệ số nhám Manning; 
V - vận tốc trung bình mặt cắt, 
A
Q
V = ; 
RH - bán kính thuỷ lực; 
B - chiều rộng mặt thoáng; 
A - diện tích mặt cắt −ớt; 
Fr - số Frút, 
gy
V
Fr = ; 
Q - l−u l−ợng (m3/s). 
 Sai phân ta có: y
AR
Qn
S
gA
BQ
1
L
23
4
H
22
0
3
2
Δ
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
−
−
=Δ (3) 
Sơ đồ khối: 
Kết thúc
Bắt đầu 
yΔ
yyy i1i Δ+=+
i
i
Hi
i
22
i
2
i
2
i
P
A
R
180/2rP
)2sin(
2
r
360
2r
A
)ry(r2B
=
πìθì=
θ−θπ=
−−ì=
2
RR
R,
2
BB
B,
2
AA
A 1HiHiHi
1ii
i
1ii
i
+++ +=+=+= 
y
AR
Qn
S
A81,9
BQ
1
L
2
i
3
4
Hi
22
0
3
i
i
2
1i Δì
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
−
ì
−
=Δ + 
1i
i
1
i1i LLL ++ Δ+Δ= ∑ 
Vẽ quan hệ L = f(y)
In kết quả 
 CT 2
 • Ph−ơng pháp 2 
Từ ph−ơng trình vi phân f0 SSdx
de −= chuyển thành ph−ơng trình sai phân f0 SSL
e −=Δ
Δ
. 
Dẫn đến: 
f0 SS
e
L −
Δ=Δ 
trong đó: e - năng l−ợng đơn vị mặt cắt, 
2
i
2
ii
gA2
Q
ye += ; fS - độ dốc thuỷ lực trung bình. 
Sơ đồ khối: cấu trúc ch−ơng trình Matlab về cơ bản không thay đổi ta chỉ thay đổi công thức 
tính độ tăng chiều dài . iLΔ
Bắt đầu 
yΔ
yyy i1i Δ+=+
i
i
Hi
22
i
i
2
i
2
i
P
A
R;)2sin(
2
r
360
2r
A
180/2rP;)ry(r2B
=θ−θπ=
πìθì=−−ì=
2
RR
R,
2
BB
B,
2
AA
A 1HiHiHi
1ii
i
1ii
i
+++ +=+=+= 
⎥⎥
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎢⎢
⎣
⎡
−
⎟⎟⎠
⎞
⎜⎜⎝
⎛ +−⎟⎟⎠
⎞
⎜⎜⎝
⎛ +
=Δ +
+
+
3
4
Hi
2
i
22
0
2
i
2
i2
1i
2
1i
1i
RA
Qn
S
gA2
Q
y
gA2
Q
y
L 
1i
i
1
i1i LLL ++ Δ+Δ= ∑ 
Vẽ quan hệ L= f(y); In kết quả 
Kết thúc 
CT 2 
 • Ph−ơng pháp 3: Ph−ơng pháp tính lặp (TRAP) 
[ ])1i(F)i(Fl
2
1
yy 1ii1i ++Δ+= ++ 
suy ra 
)1i(F)i(F
y2
l i1i ++
Δì=Δ + 
với 
i3
i
2
fi0
B
gA
Q
1
SS
)i(F
−
−= 
Ch−ơng trình Matlab chỉ thay đổi công thức tính iLΔ : 
Bắt đầu
yΔ
yyy i1i Δ+=+
i
i
Hi
i
22
i
2
i
2
i
P
A
R
180/2rP
)2sin(
2
r
360
2r
A
)ry(r2B
=
πìθì=
θ−θπ=
−−ì=
i3
i
2
fi0
B
gA
Q
1
SS
)i(F
−
−= 
)1i(F)i(F
y2
L i1i ++
Δì=Δ + 
1i
i
1
i1i LLL ++ Δ+Δ= ∑ 
Vẽ quan hệ L = f(y) 
In kết quả 
Kết thúc 
 CT 2
 Ví dụ: 
Một dòng chảy có l−u l−ợng 10 m3/s dọc theo một cống tròn có bán kính r = 1 m. Nếu tại 
mặt cắt 1, chiều sâu dòng chảy là 0,8 m. Tính chiều dài kênh từ mặt cắt 1 đến mặt cắt 2 có độ 
sâu 1,2 m về phía hạ l−u? Kênh làm bằng bêtông (n = 0,014) có độ dốc dọc không đổi S0 = 
0,003. 
• Ch−ơng trình Matlab 1 
% ----------------Tinh duong mat nuoc trong kenh ho------------------- 
clear all;close all; 
c=pi./180; 
delta_y=input('do tang chieu sau(m)\n'); 
r=input('ban kinh mat cat kenh(m)\n'); 
n=input('do nham kenh\n'); 
q=input('luu luong kenh m^3/s\n'); 
so=input('do doc day kenh\n'); 
y=[0.80: delta_y: 1.20]; 
%chieu sau dong chay thay doi tu 0.8 den 1.2m 
for i=1:length(y); 
 al_pha=(180./pi)*acos((y(i)-r)/r); 
 the_ta=180-al_pha; 
 b(i)=2*sqrt(r.^2-(y(i)-r).^2); 
 %chieu rong dinh mat cat kenh 
 A(i)=r.^2*the_ta*c-r.^2/2.*sin(2*the_ta*c); 
 %dien tich mat cat ngang 
 P(i)=r.*2.*the_ta*c; 
 %chu vi uot cua mat cat kenh 
 rh(i)=A(i)./P(i); 
 %ban kinh thuy luc 
CT 2 end 
delta_l(1)=0; 
l(1)=0; 
%cac gia tri ban dau trong day "delta_l" va "l" bang 0 
for i=1:length(y)-1; 
 A_av(i)=(A(i)+A(i+1))./2; 
 b_av(i)=(b(i)+b(i+1))./2; 
 rh_av(i)=(rh(i)+rh(i+1))./2; 
 %gia tri trung binh dt,be rong dinh,ban kinh thuy luc 
 delta_l(i+1)=((1-(q.^2.*b_av(i))./(9.81.*A_av(i).^3))./(so-
(n).^2.*(q.^2./(rh_av(i).^(4./3).*A_av(i).^2)))).*delta_y; 
 %xac dinh chieu dai theo y va delta_y dua tren cac gia tri trung binh 
 l(i+1)=sum(delta_l(1:i))+delta_l(i+1); 
 %chieu dai kenh 
end 
plot(l,y);%Ve quan he l va y 
grid; 
xlabel('chieu dai kenh'); 
ylabel('chieu cao chat long'); 
title('dong chay trong kenh ho'); 
axis([0 100 2 4]); 
fprintf('\n\n\nsection y, m delta_y, m delta_l, m l, m\n'); 
for i=1:length(y); 
 fprintf('%1.0f %2.2f %1.2f %4.2f %4.2f\n', i, y(i), 
delta_y, delta_l(i), l(i)); 
end 
%---------------------Ket thuc---------------------------- 
 • Nhập dữ liệu: 
do tang chieu sau (m): 0.1 
ban kinh mat cat kenh (m): 1 
do nham kenh: 0.014 
luu luong kenh m^3/s: 10 
do doc day kenh: 0.003 
Kết quả: 
 Section y(m) delta_y(m) delta_l(m) l(m) 
 1 0.80 0.10 0.00 0.00 
 2 0.90 0.10 27.09 27.09 
 3 1.00 0.10 25.79 52.88 
 4 1.10 0.10 23.98 76.86 
 5 1.20 0.10 21.51 98.37 
0 10 20 30 40 50 60 70 80 90 100
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
chieu dai kenh
ch
ie
u 
ca
o 
ch
at
 lo
ng
dong chay trong kenh ho
CT 2
• Ch−ơng trình Mablab 2: 
% ----------------Tinh duong mat nuoc trong kenh ho------------------- 
clear all;close all; 
c=pi./180; 
delta_y=input('do tang chieu sau(m)\n'); 
r=input('ban kinh mat cat kenh(m)\n'); 
n=input('do nham kenh\n'); 
 q=input('luu luong kenh m^3/s\n'); 
so=input('do doc day kenh\n'); 
y=[0.80: delta_y: 1.20]; 
%chieu sau dong chay thay doi tu 0.8 den 1.2m 
for i=1:length(y); 
 al_pha=(180./pi)*acos((y(i)-r)/r); 
 the_ta=180-al_pha; 
 b(i)=2*sqrt(r.^2-(y(i)-r).^2); 
 %chieu rong dinh mat cat kenh 
 A(i)=r.^2*the_ta*c-r.^2/2.*sin(2*the_ta*c); 
 %dien tich mat cat ngang 
 P(i)=r.*2.*the_ta*c; 
 %chu vi uot cua mat cat kenh 
 rh(i)=A(i)./P(i); 
 %ban kinh thuy luc 
end 
delta_l(1)=0; 
l(1)=0; 
%cac gia tri ban dau trong day "delta_l" va "l" bang 0 
for i=1:length(y)-1; 
 A_av(i)=(A(i)+A(i+1))./2; 
 b_av(i)=(b(i)+b(i+1))./2; 
 rh_av(i)=(rh(i)+rh(i+1))./2; 
 %gia tri trung binh dt,be rong dinh,ban kinh thuy luc 
 y(i+1)=y(i)+delta_y; 
 delta_l(i+1)=(((y(i+1)+q.^2./(2*9.81.*A(i+1).^2))-
(y(i)+q.^2./(2*9.81*A(i).^2)))/(so-
(n.^2*q.^2./(A_av(i).^2*rh_av(i).^(4./3))))); 
 %xac dinh chieu dai theo y va delta_y dua tren cac gia tri trung 
binh CT 2 
 l(i+1)=sum(delta_l(1:i))+delta_l(i+1); 
 %chieu dai kenh 
end 
plot(l,y);%Ve quan he l va y 
grid; 
xlabel('chieu dai kenh'); 
ylabel('chieu cao chat long'); 
title('dong chay trong kenh ho'); 
axis([0 100 0 2]); 
fprintf('\n\n\nsection y, m delta_y, m delta_l, m l, m\n'); 
for i=1:length(y); 
 fprintf('%1.0f %2.2f %1.2f %4.2f %4.2f\n', i, 
y(i), delta_y, delta_l(i), l(i)); 
end 
%---------------------Ket thuc---------------------------- 
• Nhập dữ liệu: 
do tang chieu sau (m): 0.1 
ban kinh mat cat kenh (m): 1 
do nham kenh: 0.014 
luu luong kenh m^3/s: 10 
 do doc day kenh: 0.003 
Kết quả: 
 Section y(m) delta_y(m) delta_l(m) l(m) 
1 0.80 0.10 0.00 0.00 
2 0.90 0.10 27.48 27.48 
3 1.00 0.10 26.10 53.58 
4 1.10 0.10 24.23 77.81 
5 1.20 0.10 21.72 99.53 
0 10 20 30 40 50 60 70 80 90 100
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
chieu dai kenh
ch
ie
u 
ca
o 
ch
at
 lo
ng
dong chay trong kenh ho
CT 2
• Ch−ơng trình Matlab 3: 
% ----------------Tinh duong mat nuoc trong kenh ho------------------- 
clear all;close all; 
c=pi./180; 
delta_y=input('do tang chieu sau(m)\n'); 
r=input('ban kinh mat cat kenh(m)\n'); 
n=input('do nham kenh\n'); 
q=input('luu luong kenh m^3/s\n'); 
so=input('do doc day kenh\n'); 
y=[0.80: delta_y: 1.20]; 
%chieu sau dong chay thay doi tu 0.8 den 1.2m 
for i=1:length(y); 
 al_pha=(180./pi)*acos((y(i)-r)/r); 
 the_ta=180-al_pha; 
 b(i)=2*sqrt(r.^2-(y(i)-r).^2); 
 %chieu rong dinh mat cat kenh 
 A(i)=r.^2*the_ta*c-r.^2/2.*sin(2*the_ta*c); 
 %dien tich mat cat ngang 
 P(i)=r.*2.*the_ta*c; 
 %chu vi uot cua mat cat kenh 
 rh(i)=A(i)./P(i); 
 %ban kinh thuy luc 
 F(i)=(so-(n).^2.*(q.^2./(rh(i).^(4./3).*A(i).^2)))./((1-
(q.^2.*b(i))./(9.81.*A(i).^3))); 
end 
delta_l(1)=0; 
l(1)=0; 
%cac gia tri ban dau trong day "delta_l" va "l" bang 0 
for i=1:length(y)-1; 
 delta_l(i+1)=2*delta_y./(F(i)+F(i+1)); 
 %xac dinh chieu dai theo y va delta_y 
 l(i+1)=sum(delta_l(1:i))+delta_l(i+1); 
 %chieu dai kenh 
end 
plot(l,y);%Ve quan he l va y 
grid; 
xlabel('chieu dai kenh'); 
ylabel('chieu cao chat long'); 
title('dong chay trong kenh ho'); 
axis([0 100 0 2]); 
fprintf('\n\n\nsection y, m delta_y, m delta_l, m l, m\n'); 
for i=1:length(y); 
 fprintf('%1.0f %2.2f %1.2f %4.2f %4.2f\n', i, 
y(i), delta_y, delta_l(i), l(i)); 
CT 2 
end 
%---------------------Ket thuc---------------------------- 
Nhập dữ liệu: 
do tang chieu sau (m): 0.1 
ban kinh mat cat kenh (m): 1 
do nham kenh: 0.014 
luu luong kenh m^3/s: 10 
do doc day kenh: 0.003 
Section y (m) delta_y (m) delta_l (m) l (m) 
 1 0.80 0.10 0.00 0.00 
 2 0.90 0.10 27.11 27.11 
 3 1.00 0.10 25.79 52.90 
 4 1.10 0.10 23.93 76.83 
 5 1.20 0.10 21.38 98.21 
 0 10 20 30 40 50 60 70 80 90 100
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
chieu dai kenh
ch
ie
u 
ca
o 
ch
at
 lo
ng
dong chay trong kenh ho
III. Kết luận CT 2
Ph−ơng pháp 1 và 3 cho kết quả khá giống nhau; ph−ơng pháp 1 cho kết quả l = 98,37 m, 
ph−ơng pháp 3 cho kết quả l = 98,21 m. 
Nh− vậy, cùng một cấu trúc chung của ch−ơng trình Matlab tính toán đ−ờng mặt n−ớc có 
thể đ−ợc cải biến đ−a vào nhiều ph−ơng pháp giải khác nhau để so sánh, đánh giá và chọn 1 
ph−ơng pháp thích hợp. 
Ví dụ trên ta thấy đ−ợc sự lợi ích của ch−ơng trình máy tính, nó có thể giúp ta giải quyết 
hàng loạt bài toán t−ơng tự bằng cách thay đổi dữ kiện đầu vào và thay đổi công thức tính phù 
hợp với bài toán, −u điểm của việc giải bài toán trên máy tính là cho kết quả khá chi tiết và hình 
ảnh trực quan sinh động qua đó giúp cho ng−ời học tiếp cận dần với ph−ơng pháp học mới. 
Tài liệu tham khảo 
[1] Irving H. Shames. Mechanics of Fluid. NXB Mc Graw – Hill, 2003. 
[2] Young W. Hwon & Hyochoong Bang. The Finite Element Method Using Matlab, 2000. 
[3] Phùng Văn Kh−ơng, Trần Đình Nghiên (chủ biên), Bùi Thị Vinh. Thuỷ lực. NXB Giao thông 
Vận tải, 1994Ă