ứ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ở Ks. Nguyễn viết thanh
Bộ môn Thuỷ lực - Thuỷ văn
Khoa Công trình
Trờng Đại học Giao thông Vận tải
Tóm tắt: ứng dụng phần mềm để giải quyết các bi toán kỹ thuật l một xu hớng tất yếu
hiện nay. Matlab l một trong những phần mềm có khả năng ứng dụng cao v rất tiện ích trong
các trờng đại học trên thế giới hiện nay. Bi báo ny giới thiệu một số ứng dụng của Matlab
tính toán đờng mặt nớc dòng chảy trong kênh hở.
Summary: Using software to solve the technical problems is a current tendency.
The Matlab program is one of some software which are usually applied by many universities
in the world. This article presents some applications of Matlab in the open channel flow
problems.

CT 2
I. Đặt vấn đề
Matlab là tổ hợp một cấu trúc chơng trình với các hàm toán học đã xác định trớc. ứng
dụng của Matlab rất đa dạng ta có thể xây dựng các chơng trình khác nhau để giải các bài
toán kỹ thuật, cụ thể đối với bài toán tính và vẽ đờng mặt nớc trong kênh hở ta có thể viết
chơng trình theo sơ đồ khối:

Tính chiều sâu
dọc theo dòng
chảy
Nhập
dữ liệu

Q
2
1
0
3
2
H
=
3
4
H
22
2
3
4
H
22
0
R
Vn
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
2
3

=


=


=
(2)
trong đó: S
0
- độ 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 =
;
R
H
- 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 - lu lợng (m
3
/s).

⎢
⎢
⎢
⎣
⎡
−
−
=Δ
(3)
S¬ ®å khèi:
KÕt thóc
B¾t ®Çu
y
Δ
yyy
i1i
Δ
+
=
+
i
i
Hi
i
22
i
2
i
2
i

i
1ii
i
+++
+
=
+
=
+
=

y
AR
Qn
S
A81,9
BQ
1
L
2
i
3
4
Hi
22
0
3
i
i
2

Δ+Δ=
∑

VÏ quan hÖ L = f
(y)
In kÕt qu¶

CT 2


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 +=
;
f
S - độ 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 .
i
L

Bắt đầu
y

yyy
i1i

R,
2
BB
B,
2
AA
A
1HiHi
Hi
1ii
i
1ii
i
+++
+
=
+
=
+
=











3
4
Hi
2
i
22
0
2
i
2
i
2
1i
2
1i
1i
RA
Qn
S
gA2
Q
y
gA2
Q
y
L

1i
i
1

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
i
1i
++
ì
=
+


i
i
Hi
i
22
i
2
i
2
i
P
A
R
180/2rP
)2sin(
2
r
360
2r
A
)ry(r2B
=
ìì=


=
ì=

i
3



Vẽ quan hệ L = f(y)
In kết quả
Kết thúc
CT 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');

1.6
1.8
2
chieu dai kenh
chieu cao chat long
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

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

(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

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
Ă