Matlab 在圓形斷面特征水深計算中的應用

李 蕊

（楊凌職業(yè)技術學院，陜西 楊凌 712100）

1 引言

水力計算中經(jīng)常會遇到圓形斷面正常水深、臨界水深和收縮水深的求解問題，其實質(zhì)是求解含參變量的非線性方程或超越方程，此類方程沒有求根公式，不能直接計算，傳統(tǒng)的計算方法主要有：①試算法，即利用試探法與二分法進行試算，計算繁雜，工作量大；②圖表法，即利用已經(jīng)制好的大量圖表求解，由于此類圖表種類多，使用麻煩，而且精度欠佳。

近十多年來，國內(nèi)外專家學者針對圓形斷面特征水深的計算問題，進行了大量研究，他們的研究成果主要集中在兩個方面：①引入無量綱水深，通過對特征水深方程的數(shù)學變換，得到無量綱水深的近似直接計算公式[1～9]；②采用逐次逼近原理進行迭代計算，通過選取合適的迭代初值和迭代方程計算特征水深[10～12]；另外，也有部分專家學者采用其它方法計算特征水深，如殷彥平等[13]將圓形斷面正常水深問題轉(zhuǎn)化為非線性約束優(yōu)化問題，應用混合模式搜索法求解水深；張新燕等[14]利用非線性函數(shù)模型，通過Marquardt 法建立了正常水深的直接計算公式。另外，這些研究成果大部分都集中在圓形斷面的臨界水深和正常水深的計算中，對于圓形斷面收縮水深的計算研究較少。

非線性代數(shù)方程的求解大部分都可以通過數(shù)學軟件來實現(xiàn)，Matlab 以其強大的編程及計算功能而被廣泛地應用于求解非線性代數(shù)方程中。

Matlab 是集數(shù)值分析與計算、微積分與矩陣運算、工程與科學繪圖、數(shù)字圖像處理、數(shù)字信號處理、語言編程于一體的一款工程軟件。Matlab 操作簡單，易于掌握。文中采用Matlab 中查找函數(shù)零點的命令fzero 及語言編程，對圓形斷面臨界水深、正常水深和縮水深進行編程計算，其程序簡潔明了，易于操作，而且效率和精度都非常高。

2 用Matlab 求解圓形斷面特征水深

2.1 圓形斷面臨界水深的計算

2.1.1 圓形斷面臨界水深的求解公式

水力學中臨界流的基本方程為：

如圖1 所示，圓形斷面的水力要素分別為：

圖1 圓形過水斷面

過水斷面面積：

水面寬度：

臨界水深：

式中：Q 為過水流量，m3/s；Ac為臨界流對應的過水斷面面積，m2；Bc為水面寬，m；g 為重力加速度，通常取9.81 m3/s；a 為流速分布不均勻系數(shù)，通常取1.0；θ 為臨界水深對應的圓心角，rad；d 為圓形斷面直徑，m。

將式（2）、式（3）代入式（1）得：

將上式變形得：

由此可見：式（6）為關于θ 的含參數(shù)的超越方程，理論上無解析解。因此可以利用Matlab 編程求出θ，然后代入式（4）可求出臨界水深hc。

2.1.2 工程實例

以文獻[3]為例，某圓形斷面的引水式電站輸水隧洞，洞徑d=15.0 m，求設計流量Q=1 m3/s 時的臨界水深。

利用Matlab 求解圓形斷面臨界水深程序如下：

在Matlab 語言中：alpha 表示α，theta 表示θ.

>>syms alpha Q g d theta

>>alpha=1.0；

>>Q=1；

>>g=9.81；

>>d=15；

>>f=@ (theta)(512*alpha*Q^2*sin (theta/2)./(g*d^5)).^(1/3)-theta+sin(theta)；

>>theta=fzero(f,[0.000001 2*pi])

theta=0.5440

>>hc=1/2*d*(1-cos(theta/2))

hc=0.2757

即臨界水深hc=0.2757。

用孫建公式、王正中公式和趙延風公式分別計算本例，結(jié)果見表1。

表1 臨界水深不同計算公式誤差比較

2.2 圓形斷面正常水深的計算

2.2.1 圓形斷面正常水深的求解公式

水力學中圓形斷面均勻流方程為：

圓形斷面的水力要素為：

①過水斷面面積，根據(jù)式（2）計算。

②濕周：

③正常水深：

式中：n 為粗糟系數(shù)；Q 為過水流量，m3/s；i 為底坡；Ac為發(fā)生均勻流時的過水面積，m2；Bc為水面寬度，m；p 為濕周；hc為均勻流水深，m；d 為圓形斷面直徑，m；θ 發(fā)生均勻流時的圓心角，rad。

將式（2）、式（9）代入式（8）中得：

將上式變形得：

由此可見，式（12）為關于θ 的含參數(shù)的超越方程，理論上無解析解。因此可以利用Matlab 編程求出θ，然后代入式（10）可求出正常水深hk。

2.2.2 工程實例

以文獻[7]為例，某圓形斷面的引水式電站輸水隧洞，已知斷面底坡，i=0.001 粗糟系數(shù)n=0.015，洞徑d=15.0 m，求設計流量Q=840 m3/s 時的正常水深。利用Matlab 求解圓形斷面正常水深程序如下：

在Matlab 語言中：theta 表示θ，i 表示虛數(shù)單位，所以i 用I 表示，以示區(qū)別。

>>syms n Q I d theta

>>n=0.015;

>>Q=840;

>>I=0.001;

>>d=15;

>>f=@ (theta)(2^2.6* (n*Q/sqrt(I))^0.6*theta^0.4./d^1.6)-theta+

sin(theta);

>>theta=fzero(f,[0.000001 2*pi])

theta=4.2640

>>hk=1/2*d*(1-cos(theta/2))

hk=11.4915

即正常水深hk=11.4915。

用文獻[1]和文獻[2]中的公式分別計算本例，結(jié)果見表2。

表2 正常水深不同計算公式誤差比較

2.3 圓形斷面收縮水深的計算

2.3.1 圓形斷面收縮水深的求解公式

水力學中圓形斷面收縮流的方程為：

圓形斷面的水力要素為：

①過水斷面面積，計算見式（2）。

②收縮水深：

式中：E0為上游斷面總水頭，m；Q 為過水流量，m3/s；hs為收縮水深，m；d 為圓形斷面直徑，m；φ為流速系數(shù)；θ 發(fā)生收縮流時的圓心角，rad。

將式（2）、式（14）代入式（13）中得：

將上式變形得：

求出θ 后，代入式（14）可求出收縮水深。

由此可見，式（16）為關于θ 的含參數(shù)的超越方程，理論上無解析解。因此可以利用Matlab 編程求出θ，然后代入式（14）可求出收縮水深.

2.3.2 工程實例

以文獻[9]為例，已知壩(閘)前斷面總水頭E0=12 m，圓形斷面直徑d=15.0 m，流速系數(shù)φ=0.95,求設計流量Q=500 m3/s 時的收縮水深。

利用Matlab 求解圓形斷面臨界水深程序如下：

在Matlab 語言中：phi 表示φ，theta 表示θ。

>>syms phi Q g d theta E0

>>g=9.81;

>>phi=0.95;

>>Q=200;

>>d=15;

>>E0=12;

>>f=@ (theta) (32*Q^2./ ((E0-d* (1-cos (theta/2))./2)*g*phi^2*d^4))^(1/2)-theta+sin(theta);

>>theta=fzero(f,[0.000001 4.39])

theta=1.5367

>>hs=0.5*d*(1-cos(theta/2))

hs=2.1071

即收縮水深hs=2.1071。用趙延風公式計算本例，結(jié)果見表3。

表3 收縮水深不同計算公式誤差比較

3 結(jié)論

從表1、表2 和表3 的誤差比較可以看出，應用Matlab 數(shù)學軟件求解的圓形斷面的正常水深、臨界水深和收縮水深，不僅程序簡單明了，而且計算精度高，方法更容易掌握。Matlab 作為一種強大的工程軟件，必將會廣泛應用在水力計算和水利設計中。