基于水動(dòng)力模擬的湖庫水質(zhì)改善案例應(yīng)用

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

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

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

1 水動(dòng)力模型

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

1.1 水動(dòng)力模型的基本控制方程

湖庫內(nèi)水體的運(yùn)動(dòng)受質(zhì)量守恒定律、動(dòng)量守恒定律和能量守恒定律等物理學(xué)基本守恒定律的支配。水動(dòng)力模型的基本控制方程就是這些守恒定律在水體中的數(shù)學(xué)描述[4]。

1.1.1 連續(xù)性方程

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

(1)

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

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

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

1.1.2 Navie-Stokes方程

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

(2)

(3)

(4)

其中：t——時(shí)間，s；

ρ——流體的質(zhì)量密度，kg/m3；

p——壓力，Pa；

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

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

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

1.1.3 能量守恒方程

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

(5)

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

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

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

1.2 水動(dòng)力模型的數(shù)值模擬方法

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

1.2.1 空間離散方法

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

1.2.2 時(shí)間離散方法

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

1.2.3 湍流模型

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

1.3 水動(dòng)力模型的邊界條件

在數(shù)值模擬中，各種物理場對計(jì)算區(qū)域的作用和影響需要通過添加邊界條件來實(shí)現(xiàn)。因此，要建立有效的水動(dòng)力數(shù)學(xué)模型，邊界條件的準(zhǔn)確選取至關(guān)重要。在數(shù)學(xué)模型中，第一類邊界條件是指給出邊界上待求變量的分布；第二類邊界條件是指給出邊界上待求變量的梯度值；第三類邊界條件是指給出待求變量與梯度值間的函數(shù)關(guān)系[9]。

2 水動(dòng)力模型的應(yīng)用

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

2.1 MIKE 3 FM模型的基本原理

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

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

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

在X方向追趕時(shí)，求解連續(xù)性方程和X方向動(dòng)量方程，P從n-1/6時(shí)刻計(jì)算到n+1/2時(shí)刻，U從n時(shí)刻計(jì)算到n+1時(shí)刻。V和W采用兩個(gè)時(shí)間層的已知值，V采用n-2/3和n+1/3時(shí)刻的值，而W采用n-1/3和n+2/3時(shí)刻的值。

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

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

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

2.2 MIKE 3 FM模型的建立

通過輸入數(shù)據(jù)建立模型，輸入數(shù)據(jù)的多寡取決于工程精度要求和所需描述的物理現(xiàn)象本身，通常分為以下幾個(gè)部分：計(jì)算域和相關(guān)時(shí)間參數(shù)，包括網(wǎng)格地形及時(shí)間設(shè)置；校準(zhǔn)要素，包括底床阻力、渦黏系數(shù)和風(fēng)摩擦阻力系數(shù)；初始條件，如水面高程；邊界條件，包括開邊界條件和閉邊界條件；其他驅(qū)動(dòng)力，包括風(fēng)速風(fēng)向、源匯項(xiàng)和波浪輻射應(yīng)力等。

3 工程實(shí)施案例

3.1 工程概況

曹妃甸蓄水池為人造蓄水池，是當(dāng)?shù)毓┧膽?yīng)急備用水源地。蓄水池的平面圖如圖3所示，設(shè)計(jì)蓄水深度為8 m，蓄水能力為94萬m3，池體為混凝土材質(zhì)的硬質(zhì)界面。蓄水池進(jìn)水口位于東北角距離最近頂點(diǎn)146 m處，出水口位于蓄水池東南角距離最近頂點(diǎn)47 m處。

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

為探究蓄水池進(jìn)出水口位置設(shè)計(jì)對出水水質(zhì)的影響和進(jìn)水所含污染物在蓄水池的轉(zhuǎn)移分布情況，對曹妃甸蓄水池進(jìn)行三維水動(dòng)力模擬。

3.2 曹妃甸蓄水池水動(dòng)力模型的建立

蓄水池深度分布如圖4所示，其中，模型水域邊界基于設(shè)計(jì)圖紙，深度基于現(xiàn)狀實(shí)測數(shù)據(jù)。將蓄水池水平向剖分成26 456個(gè)網(wǎng)格，得到蓄水池水動(dòng)力模型水平向的計(jì)算網(wǎng)格(圖5)。

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

圖5 三維水動(dòng)力模型水平向的計(jì)算網(wǎng)格Fig.5 Grid Diagram of 3D Hydrodynamic Model in Horizontal Direction

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

表1 模型設(shè)置Tab.1 Model Setup

3.3 曹妃甸蓄水池水質(zhì)提升模擬工況設(shè)計(jì)

由于蓄水池水體生態(tài)系統(tǒng)脆弱，流動(dòng)性較差，且外源污染輸入，導(dǎo)致水體易出現(xiàn)藻類暴發(fā)，春夏季節(jié)水體pH顯著上升，夏季底泥釋放，易造成底層水質(zhì)受到有機(jī)物、錳、氨氮的復(fù)合污染。在控制外源污染輸入的情況下，改善水流流態(tài)，可以改善水質(zhì)狀況。改善水流流態(tài)主要包括流速控制、水深控制以及換水周期控制等。

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

現(xiàn)狀條件下，進(jìn)出水口的位置及優(yōu)化設(shè)計(jì)的進(jìn)出水口位置如圖6所示。

圖6 (a)現(xiàn)狀下進(jìn)出水口位置；(b)設(shè)計(jì)進(jìn)出水口位置Fig.6 (a) Existing Inlet and Outlet Location; (b) Design Inlet and Outlet Location

3.4 結(jié)果與討論

3.4.1 現(xiàn)狀方案下蓄水池水體水動(dòng)力特征及換水周期

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

表2 模型計(jì)算工況Tab.2 Working Conditions of Model Calculation

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

由圖10可知，在無風(fēng)工況下，蓄水池同時(shí)或錯(cuò)時(shí)進(jìn)出水對水體流態(tài)及換水周期的影響非常有限，換水周期呈現(xiàn)從進(jìn)水口至出水口的直線距離線性增加的趨勢，遠(yuǎn)端水體因水流滯流換水周期超過365 d。

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

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

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

圖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可知，雖然進(jìn)、出水位置的調(diào)整改善了整個(gè)蓄水池的流態(tài)，降低了水池的滯流區(qū)面積，但蓄水池日常補(bǔ)水流量相對于水池存蓄水量來說較小，蓄水池的整體換水周期由補(bǔ)水流量限制，大部分水域換水周期仍超過365 d。

由圖14可知，在實(shí)測風(fēng)速風(fēng)向條件下，蓄水池水體在各個(gè)風(fēng)場的驅(qū)動(dòng)下混合非常充分，蓄水池各水域基本處在230～250 d。由圖14可知，將蓄水池的進(jìn)出口位置調(diào)至對角，可減小部分水域的換水周期。

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

3.4.3 實(shí)際風(fēng)場條件下現(xiàn)狀與對比方案的換水周期計(jì)算

本節(jié)選用蓄水池附近測站2017年9月—2018年8月的實(shí)測風(fēng)速、風(fēng)向數(shù)據(jù)，對現(xiàn)狀方案與對比方案在實(shí)測風(fēng)場下的換水周期進(jìn)行校核計(jì)算。蓄水池附近測站的風(fēng)速、風(fēng)向的年過程線如圖15所示。

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

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

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

3.4.4 現(xiàn)狀方案下進(jìn)水污染物在蓄水池中的轉(zhuǎn)移分布

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

圖17 工況7出現(xiàn)高進(jìn)水濃度的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 結(jié)論與建議

通過對水動(dòng)力模型的介紹以及水動(dòng)力模擬在曹妃甸蓄水池水質(zhì)提升工程中的應(yīng)用，可以得出以下結(jié)論。

(1)水動(dòng)力模型的理論體系日趨完善，水動(dòng)力模擬的實(shí)際應(yīng)用得到越來越多人的認(rèn)可。水動(dòng)力模擬將在水利工程的規(guī)劃、設(shè)計(jì)、施工、管理等各個(gè)環(huán)節(jié)中發(fā)揮著巨大的作用。

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

(3)利用水動(dòng)力模型對曹妃甸蓄水池進(jìn)水污染物進(jìn)行水動(dòng)力模擬，可以得到進(jìn)水污染物在進(jìn)入蓄水池后的轉(zhuǎn)移分布情況。這為出現(xiàn)突發(fā)情況時(shí)，采取相應(yīng)的應(yīng)急措施提供了依據(jù)。