李 蕊

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

1 引言

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

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

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

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

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

2.1 圓形斷面臨界水深的計(jì)算

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

水力學(xué)中臨界流的基本方程為：

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

圖1 圓形過水?dāng)嗝?/p>

過水?dāng)嗝婷娣e：

水面寬度：

臨界水深：

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

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

將上式變形得：

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

2.1.2 工程實(shí)例

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

利用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。

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

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

2.2 圓形斷面正常水深的計(jì)算

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

水力學(xué)中圓形斷面均勻流方程為：

圓形斷面的水力要素為：

①過水?dāng)嗝婷娣e，根據(jù)式（2）計(jì)算。

②濕周：

③正常水深：

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

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

將上式變形得：

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

2.2.2 工程實(shí)例

以文獻(xiàn)[7]為例，某圓形斷面的引水式電站輸水隧洞，已知斷面底坡，i=0.001 粗糟系數(shù)n=0.015，洞徑d=15.0 m，求設(shè)計(jì)流量Q=840 m3/s 時(shí)的正常水深。利用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。

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

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

2.3 圓形斷面收縮水深的計(jì)算

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

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

圓形斷面的水力要素為：

①過水?dāng)嗝婷娣e，計(jì)算見式（2）。

②收縮水深：

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

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

將上式變形得：

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

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

2.3.2 工程實(shí)例

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

利用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。用趙延風(fēng)公式計(jì)算本例，結(jié)果見表3。

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

3 結(jié)論

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