二維水流數學模型在船廠工程防洪影響評價中的應用

邱居華，李 陽，朱士江，王海銀

(1.湖北欣瑞工程咨詢有限公司，湖北 宜昌 443000；2.湖北昊源建設工程有限公司，湖北 宜昌 443000；3.三峽大學，湖北 宜昌 443000)

1 工程概況

本船廠工程位于三峽大壩和葛洲壩兩大水利樞紐之間，長江右岸(葛洲壩庫區南岸)南沱村，距離三峽大壩下游約15.0 km。工程河段水位、流量條件直接受水利樞紐調度影響，有實測數據。

本船廠工程生產場地約33 335 m2，船臺8個(其中有4個為混凝土地面，另外4個為沙土地面)，舾裝碼頭1座，占用長江水域岸線長度約300 m，廠區分為加工區、制作區、倉儲區、人員辦公區、混凝土地面船舶修造區、沙土地面船舶修造區等；生產、測量設備及檢測設備近400臺(套)，主要生產活動為修船和造船。

2 二維水流數學模型

近年來，二維水流數學模型被運用于各種涉河建筑物的防洪影響評價計算，為工程建設提供技術支撐[1]。石少山[2]使用二維水流數學模型對親水廣場進行防洪評價計算，王麗君[3]使用二維水流數學模型對輸電線路工程進行防洪評價計算，王鑫等[4-6]使用二維水流數學模型對跨河大橋工程進行防洪評價計算，盛建興[7]、曹開興[8]等使用二維水流數學模型對碼頭工程進行防洪評價計算。

對于復雜天然河道水流運動的數值模擬，為了克服計算區域邊界復雜的困難，同時也為了提高數值計算成果的精度，多采用基于曲線網格的坐標變換方法，使計算網格貼合曲折的河道邊界。其中正交曲線變換和一般(非正交)曲線變換方法是兩種最常用的方法。與正交曲線變換相比，一般曲線變換不受計算網格必須嚴格正交的限制，網格生成也較靈活。因此本次采用一般曲線坐標變換方法實現計算邊界與物理邊界的準確耦合。

2.1 模型計算控制方程

考慮側向入匯流影響，笛卡爾坐標系下平面二維水流數學模型的控制方程如下：

水流連續方程如式(1)：

h=n+d

(1)

x方向水流運動方程如式(2)：

(2)

式中：f為柯氏力系數；η為河底高程，m；Pa為當地大氣壓，MPa；g為重力加速度，m/s2；ρ為水的密度，kg/m3；ρ0為參照水密度，kg/m3；τsx、τbx分別為風場摩擦力在x、y上的分量，MPa；sxx、sxy分別為輻射應力分量，MPa；Txx、Txy為水平黏滯力應力分量，MPa；uS為源項x方向水流流速，m/s。

y方向水流運動方程如式(3)：

(3)

式中：τsy、τby分別為底床摩擦力在x、y上的分量；syx、syy分別為輻射應力分量；Tyy為水平黏滯力應力分量；vS為源項y方向水流流速，m/s。

2.2 模型方程離散和求解

MIKE21模型的數值解法主要有空間離散方程組成，計算區域的空間離散是用有限體積法(Finite Volume Method)，其求解方式是將計算區域劃分成不重復的控制單元，單元理論上可以是任意形狀的多邊形，但在這里一般分解成三角形和四邊形單元。

淺水方程組的通用形式如式(4)：

(4)

式中：U為守恒型物理向量；F(U)為通量向量，S(U)為源項S的微分函數表達式。

在笛卡爾坐標系中，二維淺水方程組如式(5)：

(5)

式中：Fx、Fy分別為F(通量向量)在x和y方向的分量；I、V分別為無黏性通量和黏性通量。各個參數的計算如式(6)～式(11)：

(6)

(7)

(8)

(9)

(10)

(11)

式中：A為過水面積，m2。

對上述淺水方程的的第i個單元進行積分，選用Gauss原理重寫后見式(12)：

(12)

式中：Ai為單元Ωi的面積；Γi為單元的邊界；ds、dΩ分別為沿著Γi、Ai邊界的積分變量；n為糙率。

式(12)運用了單點求積法來計算面積的積分，該求積點位于單元的質點，同時運用中點求積法來計算邊界積分，如式(13)：

(13)

式中：Ui為第i個單元的U的平均值，并位于單元中心；Si為第i個單元的S的平均值，并位于單元中心；NS為單元的邊界數ΔΓj為第j個單元的長度。

3 數學模型驗證

模型驗證計算的目的在于檢驗數學模型計算方法的可行性，根據實測資料確定模型中相關參數的取值，并檢驗其精度。本次主要根據水面線和垂線平均流速分布資料對模型參數進行調試。

3.1 模型計算范圍及邊界處理

(1)工程河段地形圖。采用2016年11月實測工程河段1/2000河道地形圖。

(2)水文資料。采用2014年9月3日(Q=45 100 m3/s)一組測驗水面線及斷面垂線平均流速分布資料，共驗證4個水位點水位和4個斷面垂線流速分布情況，驗證水位點及流速分布斷面位置見圖1～圖2。

圖1 數學模型計算區域地形概況

圖2 各斷面流速分布驗證(Q=45 100 m3/s)

(3)邊界條件處理。平面二維水流模型中，邊界條件通常包括河道進出口邊界、岸邊界及動邊界等。本模型中對上述邊界分別做以下處理：①進口邊界。根據已知進口全斷面流量，給定入流單寬流量沿斷面的橫向分布。②出口邊界。給定出口斷面的水位。③地形邊界。取河道地形圖外邊線為計算的岸邊界，給定其一個比較高的高程令其流速為0 m/s，以保證計算區域的封閉。④動邊界。

本模型采用“凍結”法進行動邊界處理，即根據水位結點處河底高程來判斷該網格單元是否露出水面。若不露出，糙率取正常值；反之，更改單元的糙率(n取1010量級)。為不影響水流控制方程的求解，在露出水面的結點處需給定一個薄水層，一般給定其厚度為0.5 cm。

3.2 模型驗證計算成果

模型計算水位成果與實測成果比較見表1。模型計算水位與實測水位基本吻合，其最大相差為0.2 cm。

表1 水位驗證計算成果 m

根據實測水文資料綜合調試，得到該河段糙率變化范圍為主槽n=0.030～0.036，灘地糙率n=0.035～0.046。

圖1～圖2給出了工程河段4個流速斷面垂線流速分布驗證成果，計算與實測流速大小和沿河寬分布情況基本一致，兩者誤差一般<5.0%。圖3為根據驗證條件計算得到的河段流場。模型計算得到的流場變化平順，洪枯水條件下水流運動特征明顯，灘槽水流運動區分合理，水流運動形態與河道地形變化情況符合較好。

綜合分析以上驗證計算成果，所采用的數學模型能較好地模擬工程河段水流運動特性，驗證計算成果與實測成果吻合較好，表明數學模型的計算方法正確，模型中相關參數取值合理，可以用于船廠工程影響河道水位及流速的計算分析。

4 案例應用與分析計算

4.1 水文條件

為評價船廠工程修建后對工程河段行洪的影響，計算時主要針對長江汛期高水位條件，考慮100 a一遇設計洪水56 700 m3/s流量工況，還考慮了汛期維持下游通航的三峽控泄流量45 000 m3/s、汛期常遇流量30 000 m3/s三種工況，具體計算條件見表2。

表2 各工況計算條件

4.2 工程概化

根據船廠工程設計方案，船廠工程主要阻水部分主要為：車間平臺、船體放置平臺、實體坡道、架空下河坡道、坡道前沿的躉船等。為反映工程對河道水流運動的影響，在計算過程中通過局部糙率調整的方法進行工程概化。具體為：

(1)對于實體坡道，主要是通過修改局部地形來反映工程的影響。

(2)對于架空坡道按工程量計算各水位下的阻水面積，對于躉船按照船體尺寸及吃水計算阻水面積，根據以上兩項計算成果確定工程后的河道過水面積值，通過公式將阻水轉換為局部糙率值。

(3)局部糙率計算方法：局部阻力系數ζ計算如式(14)：

ζ=0.5(1-A工程后/A工程前)

(14)

式中：A工程后、A工程前分別為工程后、工程前的過水面積，m2。

為便于計算，將局部阻力系數轉化為糙率的形式如式(15)：

(15)

式中：n阻為糙率；H為吃水深度，m。

最后得到施工區域所在網格的局部綜合糙率系數n為如式(16)：

(16)

式中：n工程前為工程前的糙率。

4.3 計算成果分析

4.3.1 工程后斷面縮窄率

為了了解船廠工程區域內工程前后橫斷面變化及斷面過水面積縮窄情況，根據船廠工程的平面布置，選擇船廠車間平臺及船體放置平臺作為代表斷面進行分析研究。以2016年11月地形圖為基礎，計算了船廠工程后的斷面縮窄率見表3。

表3 船廠面積縮窄率計算成果

計算表明，工程后在各設計洪水位下，工程區域河道的橫斷面過水面積均略有減小，隨著水位不同面積縮窄率略有不同，船廠的面積縮窄率均<1.000%。遇100 a一遇的洪水時，船廠最大面積縮窄率約為0.285%；遇流量30 000 m3/s、45 000 m3/s洪水時最大面積縮窄率變化不大，變化大小幅度分別為0.055%、0.151%。滿足《長江流域和瀾滄江以西(含瀾滄江)區域河湖管理范圍內建設項目工程建設方案洪水影響審查技術標準》(T/CTESGS 02—2022長江技術經濟學會，2023.01.01)的要求。

4.3.2 工程對水位影響分析

船廠工程建設后，河道過流斷面被束窄，河槽水位隨之出現一定調整，模型計算表明：上游不同來流條件下，船廠工程引起的河道水位變化分布規律大致相似，均呈現出上游局部水位抬高，下游局部水位降低；且水位變化影響范圍和影響程度隨流量增加有所增加，但總體變化幅度均不大。工程前后水位變化見表4、圖4～圖5。

表4 工程前后水位變化計算成果

圖4 評價區域水位變化等值線

圖5 船廠區域水位變化等值線

4.3.3 工程對流速影響分析

工程修建以后引起的流速變化主要集中在船廠區上下游、前沿局部區域。具體影響表現為船廠上游局部區域受阻水作用，流速表現為下降，船廠區域及其下游局部由于水流繞流紊動擴散，流速也相應下降，船廠前沿水流受工程擠壓，局部流速有所增加。同時計算表明，船廠對下游流速的影響大于上游，流速減小幅度和影響范圍大于流速增加幅度及影響范圍，船廠對流速的影響隨著流量增大而增大，但總體影響較小。船廠工程前后流速變化見表5、圖6～圖7。

表5 工程前后流速變化計算成果

圖6 評價區域流速變化等值線

圖7 船廠區域流速變化等值線

4.3.4 綜合評價分析

本工程建設后，河道過流斷面略被束窄，河槽水位出現小調整。模型計算表明：本工程引起的水位變化情況為船廠工程上游局部水位抬高，下游局部水位降低，汛期遭遇100 a一遇洪水時，局部最大壅水值為0.43 cm。船廠前沿流速略有增加，船廠區域流速略有減小，汛期遭遇100 a一遇洪水時，流速最大增加值為0.06 m/s。船廠工程對水位、流速的影響量和影響范圍隨流量增加而增大。

在工程影響范圍內，本工程興建對河道過水面積縮窄率、壅水高度、流速變化的影響較小，不會對河勢及河道行洪產生明顯不利的影響。

5 結 語

(1)結合水文站實測數據成果對本文計算成果進行了驗證，驗證計算成果與實測成果吻合較好，表明本文所采用的數學模型能較好地模擬工程河段水流運動特性，數學模型的計算方法正確，模型中相關參數取值合理，適用于船廠工程影響河道水位及流速的計算分析。

(2)通過二維水流數學模型計算分析看出，工程前后河道水動力特性和流場的變化不大，水流流態及水流流速分布無明顯變化，計算成果可靠，為建設項目對防洪影響的計算與分析提供了可靠的方法，為建設項目安全運行提供了保障。