基于水動力模擬的湖庫水質改善案例應用

王澤民，王巖波，沈昌明，于 蕾，李春明，郎華偉

(1. 唐山市曹妃甸供水有限責任公司, 河北唐山 063200；2. 上海同濟環境工程科技有限公司, 上海 200092)

據統計，我國1.0 km2以上的自然湖泊共2 693個，總面積為81 414.6 km2，約占全國國土面積的0.9%[1]。新中國成立以來，我國共修建各類水庫8.5萬余座，總庫容達6 000多億m3[2]。如今，湖庫等水體水質變差，富營養化問題突出，而人們對水環境的要求越來越高，越來越多的水環境治理工程提上日程。進行水環境整治投資巨大、影響深遠，為避免決策失誤，對湖庫綜合整治的各項措施進行深入分析、充分論證十分必要。應用水動力學和水質模型，對水體的水動力學規律和水質情況進行分析預測，再根據模擬預測的結果對湖庫進行有效的治理和管理，這對水環境質量提高以及生態文明建設具有十分重要的意義。

1 水動力模型

水動力學模型源于19世紀圣維南的理論，通過研究提出了圣維南方程組，奠定了非恒定流的理論基礎。20世紀中葉，計算機誕生后，水動力模型迎來了大發展。20世紀50—60年代，建立了許多一維數學模型，主要研究水流運動規律；20世紀70年代，二維水動力學模型得到充分的發展與研究應用；20世紀80年代至今，三維水動力學模型得到迅速發展，且很多模型已經成功地運用于水體[3]。

1.1 水動力模型的基本控制方程

湖庫內水體的運動受質量守恒定律、動量守恒定律和能量守恒定律等物理學基本守恒定律的支配。水動力模型的基本控制方程就是這些守恒定律在水體中的數學描述[4]。

1.1.1 連續性方程

連續性方程是質量守恒定律在流體力學中的具體表述形式，對流體采用連續介質模型，速度和密度都是空間坐標及時間的連續可微函數。不可壓縮流體的連續性微分方程如式(1)。

(1)

其中：vx——水流速度在X方向上的分量，m/s；

vy——水流速度在Y方向上的分量，m/s；

vz——水流速度在Z方向上的分量，m/s。

1.1.2 Navie-Stokes方程

通過動量守恒定律可以得到不可壓縮流體的Navie-Stokes方程(簡稱N-S方程)，其單位質量流體的運動微分方程如式(2)～式(4)。

(2)

(3)

(4)

其中：t——時間，s；

ρ——流體的質量密度，kg/m3；

p——壓力，Pa；

fx——外部體積力在X方向上的分量，N/kg；

fy——外部體積力在Y方向上的分量，N/kg；

fz——外部體積力在Z方向上的分量，N/kg。

1.1.3 能量守恒方程

對于不可壓縮流體，能量指的是機械能。理想、不可壓縮流體在重力場中做穩定流動，沿流線或者無旋流場中作流束運動時，單位重量流體的位能、壓力能和動能之和是常數，即機械能是守恒的，且它們之間可以相互轉換。表示流體能量關系的伯努利方程的微分形式如式(5)。

(5)

其中：v——流速，m/s；

g——重力加速度，g=9.81 m/s2;

z——流體的位壓頭，m。

1.2 水動力模型的數值模擬方法

將描述水體流動的控制方程網格化、離散化，并在生成的網格內通過積分等形式求得精度較高的近似解[5]。目前，主要的方法為空間離散方法[6]、時間離散方法[7]和湍流模型[8]。

1.2.1 空間離散方法

空間離散方法包括有限差分法、有限單元法和有限體積法等。有限差分法是直接將微分轉化為代數問題進行近似求解；有限單元法則是利用加權余量的方法將描述各類物理場的泊松方程轉化為求解特定泛函數的極值；有限體積法的核心控制方程是流量、動量等物理量基于積分形式的守恒方程。

1.2.2 時間離散方法

時間離散方法又稱時間積分，主要為顯式、隱式和半隱式3種格式。顯式格式的特點是在每個待求解的方程中只存在一個未知數，即可以通過直接計算來求得方程的解；隱式格式要求同時求解在同一時刻所有網格上的未知量，需要方程組聯立等方法進行求解；半隱式格式是指在同一方程中，對激發快波的項用隱式表示，對描述慢波的項用顯式表示差分格式。

1.2.3 湍流模型

在自然情況下，水體的流動形式大多為湍流。湍流是不確定性和確定性的統一體，根據湍流的性質，建立附加條件，使方程組封閉，可構成湍流模型。k方程模型是指在時均流動的偏微分方程組之外，增加一個和湍動流速尺度有關的偏微分方程，如增加變量單位質量流體的湍動動能來全面反映湍流運動等[7]。k-ε方程模型在增加湍動能k外又增加了一個確定湍動特征長度L的變量，這樣可構成湍動能量k和湍動特征長度L的的偏微分方程，或者湍動能量k和湍動能量耗損率ε的偏微分方程[7]。其中，k-ε模型在模擬計算中具有良好的穩定性和精度，廣泛應用于模擬計算的軟件中。

1.3 水動力模型的邊界條件

在數值模擬中，各種物理場對計算區域的作用和影響需要通過添加邊界條件來實現。因此，要建立有效的水動力數學模型，邊界條件的準確選取至關重要。在數學模型中，第一類邊界條件是指給出邊界上待求變量的分布；第二類邊界條件是指給出邊界上待求變量的梯度值；第三類邊界條件是指給出待求變量與梯度值間的函數關系[9]。

2 水動力模型的應用

水動力模型是水環境研究的重要工具，隨著對水體污染物研究的不斷深入以及數學手段在環境領域研究中應用，水動力模型得到了長足的發展。目前，較著名的二維模型為丹麥的MIKE 21、荷蘭的Delf3tD、美國的RMA模型等；較著名的三維模型為丹麥的MIKE 3、美國的POM模型、歐盟的COHERENS模型等。其中，MIKE 3 FM作為一款可以解決帶自由表面的三維流動問題的通用模型，可以勝任與內陸湖泊、河流及景觀水體相關的模型研究工作。

2.1 MIKE 3 FM模型的基本原理

MIKE 3 FM模型的數學基礎建立在包括了紊流影響、密度變化、鹽度和溫度平衡的雷諾平均化的N-S方程之上，采用交替方向隱式迭代法對質量及動量守恒方程進行積分，并采用雙精度掃描法對其產生的數學矩陣進行求解[10]。各差分項在X、Y和Z方向上的交錯網格布置如圖1所示，方程組的時間中心差分法如圖2所示。

圖1 X、Y、Z方向上差分網格[10]Fig.1 Difference Grid in X, Y, Z Directions[10]

圖2 時間中心差分方法[10]Fig.2 Time Center Difference Methods[10]

在X方向追趕時，求解連續性方程和X方向動量方程，P從n-1/6時刻計算到n+1/2時刻，U從n時刻計算到n+1時刻。V和W采用兩個時間層的已知值，V采用n-2/3和n+1/3時刻的值，而W采用n-1/3和n+2/3時刻的值。

在Y方向追趕時，求解連續性方程和Y方向動量方程，P從n+1/6時刻計算到n+5/6時刻，V從n+1/3時刻計算到n+4/3時刻。U采用剛從X方向追趕時獲得的n和n+1時刻的值，而W采用n-1/3和n+2/3時刻的值。

最終在Z方向追趕時，求解連續性方程和Z方向動量方程，P從n+3/6時刻計算到n+7/6時刻，W從n+2/3時刻計算到n+5/3時刻。U采用剛從X方向追趕時獲得的n和n+1時刻的值，而V采用剛從Y方向追趕時獲得的n+1/3和n+4/3時刻的值。

綜合上述3個方向的追趕過程，可保證時間中心差分位于n+1/2時刻。在X方向上的追趕過程減少了Y和Z方向維數，因此，稱為向下追趕過程，而Y和Z方向上的追趕過程則稱為向上追趕過程。采用該數值方法可保證工程應用中沒有矢量、動量和能量的失真，差分精度達到二階。

2.2 MIKE 3 FM模型的建立

通過輸入數據建立模型，輸入數據的多寡取決于工程精度要求和所需描述的物理現象本身，通常分為以下幾個部分：計算域和相關時間參數，包括網格地形及時間設置；校準要素，包括底床阻力、渦黏系數和風摩擦阻力系數；初始條件，如水面高程；邊界條件，包括開邊界條件和閉邊界條件；其他驅動力，包括風速風向、源匯項和波浪輻射應力等。

3 工程實施案例

3.1 工程概況

曹妃甸蓄水池為人造蓄水池，是當地供水的應急備用水源地。蓄水池的平面圖如圖3所示，設計蓄水深度為8 m，蓄水能力為94萬m3，池體為混凝土材質的硬質界面。蓄水池進水口位于東北角距離最近頂點146 m處，出水口位于蓄水池東南角距離最近頂點47 m處。

圖3 曹妃甸蓄水池平面圖Fig.3 Plan of Caofeidian Reservoir

為探究蓄水池進出水口位置設計對出水水質的影響和進水所含污染物在蓄水池的轉移分布情況，對曹妃甸蓄水池進行三維水動力模擬。

3.2 曹妃甸蓄水池水動力模型的建立

蓄水池深度分布如圖4所示，其中，模型水域邊界基于設計圖紙，深度基于現狀實測數據。將蓄水池水平向剖分成26 456個網格，得到蓄水池水動力模型水平向的計算網格(圖5)。

圖4 蓄水池深度分布Fig.4 Reservoir Depth Distribution

圖5 三維水動力模型水平向的計算網格Fig.5 Grid Diagram of 3D Hydrodynamic Model in Horizontal Direction

模型參數、初始條件和邊界條件設置如表1所示。

表1 模型設置Tab.1 Model Setup

3.3 曹妃甸蓄水池水質提升模擬工況設計

由于蓄水池水體生態系統脆弱，流動性較差，且外源污染輸入，導致水體易出現藻類暴發，春夏季節水體pH顯著上升，夏季底泥釋放，易造成底層水質受到有機物、錳、氨氮的復合污染。在控制外源污染輸入的情況下，改善水流流態，可以改善水質狀況。改善水流流態主要包括流速控制、水深控制以及換水周期控制等。

本研究通過在定常風條件下，對現狀及對比方案下的不同進出水工況進行水動力模擬，進而分析在實際風場條件下，現狀以及對比方案的進出水口位置設計對水流速度和換水周期的影響。同時，考慮外源污染輸入時，污染物的轉移分布情況，設計工況條件具體如表2所示。

現狀條件下，進出水口的位置及優化設計的進出水口位置如圖6所示。

圖6 (a)現狀下進出水口位置；(b)設計進出水口位置Fig.6 (a) Existing Inlet and Outlet Location; (b) Design Inlet and Outlet Location

3.4 結果與討論

3.4.1 現狀方案下蓄水池水體水動力特征及換水周期

現狀方案計算工況根據進出水方式的不同分為兩組，為表2中工況1和工況2。各工況的流速分布、矢量圖及對應的換水周期如圖7～圖10所示。

表2 模型計算工況Tab.2 Working Conditions of Model Calculation

由圖7～圖9可知，在無風工況下，由于進水流量相對于水池存蓄水量來說較小(按存蓄水量/補水流量計算得到的換水周期為183 d)，通過進出水驅動的流速量級在0.01 cm/s(取、排水口附近流速為0.03～0.05 m/s)，表、中、底層流態相近，流速大小從底層至表層逐步增大。同時，進出水與錯時進出水工況下，離進、出水口較遠的水域均存在較大面積的滯流區。

由圖10可知，在無風工況下，蓄水池同時或錯時進出水對水體流態及換水周期的影響非常有限，換水周期呈現從進水口至出水口的直線距離線性增加的趨勢，遠端水體因水流滯流換水周期超過365 d。

3.4.2 對比方案下蓄水池水體水動力特征及換水周期分析

對比方案計算工況根據進出水方式的不同分為兩組，為表2中工況3和工況4。各工況的流速分布、矢量圖及對應的換水周期如圖11～圖14所示。

對比圖11～圖13與圖7～圖9，對比方案將進、出水位置調整至對角，蓄水池整體流態明顯改善，進、出水口間的流場分布得更為均勻，有效地降低了水池的滯流區面積。

圖7 表層流速分布 (a)工況1；(b)工況2Fig.7 Diagram of Surface Layer Velocity Distribution (a) Condition 1; (b) Condition 2

圖8 中層流速分布 (a)工況1；(b)工況2Fig.8 Diagram of Middle Layer Velocity Distribution (a) Condition 1; (b) Condition 2

圖9 底層流速分布 (a)工況1；(b)工況2Fig.9 Diagram of Bottom Layer Velocity Distribution (a) Condition 1; (b) Condition 2

圖10 換水周期分布 (a)工況1；(b)工況2Fig.10 Diagram of Water Change Period Distribution (a) Condition 1; (b) Condition 2

圖11 表層流速分布 (a)工況3；(b)工況4Fig.11 Diagram of Surface Layer Velocity Distribution (a) Condition 3; (b) Condition 4

圖12 中層流速分布 (a)工況3；(b)工況4Fig.12 Diagram of Middle Layer Velocity Distribution (a) Condition 3; (b) Condition 4

圖13 底層流速分布 (a)工況3；(b)工況4Fig.13 Diagram of Bottom Layer Velocity Distribution (a) Condition 3; (b) Condition 4

對比圖14與圖11可知，雖然進、出水位置的調整改善了整個蓄水池的流態，降低了水池的滯流區面積，但蓄水池日常補水流量相對于水池存蓄水量來說較小，蓄水池的整體換水周期由補水流量限制，大部分水域換水周期仍超過365 d。

由圖14可知，在實測風速風向條件下，蓄水池水體在各個風場的驅動下混合非常充分，蓄水池各水域基本處在230～250 d。由圖14可知，將蓄水池的進出口位置調至對角，可減小部分水域的換水周期。

圖14 換水周期分布 (a)工況3；(b)工況4Fig.14 Diagram of Water Change Period (a) Condition 3; (b) Condition 4

3.4.3 實際風場條件下現狀與對比方案的換水周期計算

本節選用蓄水池附近測站2017年9月—2018年8月的實測風速、風向數據，對現狀方案與對比方案在實測風場下的換水周期進行校核計算。蓄水池附近測站的風速、風向的年過程線如圖15所示。

圖15 蓄水池附近測站的風速、風向的年過程線Fig.15 Annual Process Line of Wind Speed and Direction at a Measuring Station near the Reservoir

校核計算工況基于2017年9月—2018年8月實測風速過程，根據不同進、出水位置分為兩組工況，為表2中工況5和工況6?，F狀與對比方案在實際風場條件下計算得到的換水周期如圖16所示。

圖16 換水周期分布 (a)工況5；(b)工況6Fig.16 Diagram of Water Change Period (a) Condition 5; (b) Condition 6

3.4.4 現狀方案下進水污染物在蓄水池中的轉移分布

工況7模擬在現狀進出水口位置條件下，選擇歷史高進水濃度(CODCr為30 mg/L)、盛行風(南風4 m/s)、蓄水池同時進出水5 00 m3/h時，模擬計算高濃度水在蓄水池的遷移過程。圖17為出現高進水濃度1、2、3、4、7、14、21、28 d后的COD分布。模擬結果顯示，進水污染物隨著水流運動，不斷與蓄水稀釋混合，在出現高濃度補水14 d后輸運至整個蓄水池并造成污染。

圖17 工況7出現高進水濃度的CODCr分布 (a)1 d后；(b)2 d后；(c)3 d后；(d)4 d后； (e)7 d后；(f)14 d；(g)21 d后；(h)28 d后Fig.17 Diagram of CODCr Distribution with High Concentration Inflow in Condition 7 (a) after 1 Day; (b) after 2 Days; (c) after 3 Days; (d) after 4 Days; (e) after 7 Days; (f) after 14 Days; (g) after 21 Days; (h) after 28 Days

4 結論與建議

通過對水動力模型的介紹以及水動力模擬在曹妃甸蓄水池水質提升工程中的應用，可以得出以下結論。

(1)水動力模型的理論體系日趨完善，水動力模擬的實際應用得到越來越多人的認可。水動力模擬將在水利工程的規劃、設計、施工、管理等各個環節中發揮著巨大的作用。

(2)利用水動力模型對曹妃甸蓄水池進出水口更改前后進行水動力模擬，可以得到更改進出水口，即將進出水口調至對角，可以促進水體在蓄水池內的流動，減少了部分水域的換水周期。這為更改進出水口這一實際工程提供了有效的理論支撐。

(3)利用水動力模型對曹妃甸蓄水池進水污染物進行水動力模擬，可以得到進水污染物在進入蓄水池后的轉移分布情況。這為出現突發情況時，采取相應的應急措施提供了依據。