射陽河口污染物輸運(yùn)數(shù)值模擬

魏 廣

（南京揚(yáng)子水利自動(dòng)化技術(shù)開發(fā)總公司，江蘇 南京 210012）

0 引言

射陽河是江蘇省蘇北地區(qū)排水入海的最大天然干河，隨著經(jīng)濟(jì)發(fā)展和人口的不斷增長，以及廢水排放量逐年增加，射陽河流域水環(huán)境壓力越來越大，近海水域水環(huán)境污染日益受到了人們的重視。射陽河口處設(shè)有閘門，開閘排放污染物的時(shí)間與污染物在河口近海地區(qū)的輸運(yùn)擴(kuò)散規(guī)律和濃度分布特征有密切關(guān)系。為了有效地控制污染，合理地確定水環(huán)境容量及允許排放量，科學(xué)地制定污染物排放時(shí)間和過程，需要深入進(jìn)行水動(dòng)力狀況及污染物在水環(huán)境中的輸運(yùn)遷移規(guī)律的研究，而數(shù)值模擬是解決上述問題的重要途徑。

1 二維水流水質(zhì)耦合數(shù)學(xué)模型建立及求解

1.1 控制方程

守恒型二維淺水方程[1]與對(duì)流擴(kuò)散方程耦合的矢量表達(dá)式為：

式中：q = [ h, h u, h v, h Ci]T為守恒物理量；F (q) = [ h u,h u2+ g h2/ 2, h u v, h u Ci]T為 x 向通量；G (q) = [ h v, h u v,h v2+ g h2/ 2, h v Ci]T 為 y 向通量；b (q) = [0, g h (s0x- sfx),g h (s0y- sfx), Δ(DkΔ) h Ci)) - KCig h Ci+ Si]T為源匯項(xiàng)；x，y，t 分別為空間及時(shí)間坐標(biāo)系；h 為水深；u 和 v 分別為 x 和 y 向沿水深積分平均流速分量；g 是重力加速度；s0x和 sfx分別是 x向的水底和摩阻坡度；s0y和 sfy分別是 y 向的水底和摩阻坡度；Ci為污染物沿水深積分的垂線平均濃度；本模型包括 8 個(gè)組分（COD，NH3-N，BOD，DO，溫度，鹽度，TP 及 TN）方程，可同時(shí)也可單獨(dú)模擬。其中，KCi是污染物綜合降解系數(shù)；Si是污染物源匯項(xiàng)，其內(nèi)容隨組分而變；Dk為污染物的擴(kuò)散系數(shù)；Δ為梯度算子。

1.2 方程離散

在任意形狀的單元 Ω，對(duì)式（1）進(jìn)行積分離散求得[2]

式中：A 為單元 Ω 的面積；m 為單元邊總數(shù)； Lj為第 j 單元邊的長度。

b (q) = [0, A g h (s0x- sfx)，A g h (s0y- sfy)，∑Dk(Δ hCi)nLj+ Si- AKCihCi]T，q = T(φ) · q，T (φ) 和T (φ)-1為坐標(biāo)旋轉(zhuǎn)變換及逆變換矩陣[1-2]，其中 φ 為坐標(biāo)軸旋轉(zhuǎn)角度。

式 (2) 采用顯格式時(shí)間離散，其求解歸結(jié)為如何確定法向通量，可通過解算局部一維黎曼初值問題的外法向數(shù)值通量（記為?LR）而得到[3]。

1.3 水利設(shè)施處的通量計(jì)算

使用計(jì)算域內(nèi)的水利設(shè)施作為模型的內(nèi)邊界，水利設(shè)施處的法向通量的計(jì)算不能應(yīng)用黎曼近似解，而要采用與水利設(shè)施相應(yīng)的出流公式[3]。其通量公式概況如下：

式中：Qm是水工建筑物出流；h 為水深；un，vτ分別為法向及切向流速，滿足 un= u cos φ + v sin φ 和 vτ=v cos φ - u sin φ。模型中包含的水工建筑物為閘門。本模型中采用經(jīng)驗(yàn)公式：Qg= C1g· C2g· Hg· Ng· Lg·式中：C1g為閘門系數(shù)；C2g是淹沒系數(shù)；Hg為閘門啟閉高度；Ng是總的閘門個(gè)數(shù)；Lg是每個(gè)閘門的長度；當(dāng)尾水位高于閘門溢洪道堰頂高程時(shí)，Δ Z 為閘上水位和尾水位之差，反之為閘上水位和溢洪道堰頂高程之差。

1.4 動(dòng)邊界的處理

在河口地區(qū)，隨著潮水漲落，淺灘的漫、露交替出現(xiàn)，高潮位時(shí)淹沒，低潮位時(shí)露出。在進(jìn)行該區(qū)域數(shù)值模擬時(shí)，模型必須具備捕捉這種動(dòng)態(tài)過程的能力。動(dòng)邊界是水平計(jì)算域中有水與無水區(qū)域的界線。在該二維模型中，根據(jù)有限體積法特點(diǎn)，采用干濕動(dòng)邊界法，以水深為判別標(biāo)準(zhǔn)設(shè)定單元界面的類型，并運(yùn)用相應(yīng)的方法計(jì)算跨界面的法向通量，以保證水量平衡。在每一時(shí)間步判斷單元界面2 側(cè)單元的水深 hL，hR，當(dāng) hL和 hR均為 0 時(shí)，則通量為 0；當(dāng) hL和 hR有 1 個(gè)為 0 時(shí)，則結(jié)合 2 側(cè)單元的底高程，交界面類型可能為固壁邊界、跌水、漫流等形態(tài)的邊界，根據(jù)不同類型可選擇瞬時(shí)潰壩解析解、堰流公式等方式估算通量。當(dāng) hL和 hR均不為 0 時(shí)，則按正常的方法估算通量。在實(shí)際建模時(shí)，設(shè)定 1 個(gè)水深限值，以此來判斷單元的干濕情況。若相鄰單元的水深均小于水深限值，則認(rèn)為該單元為干單元，不進(jìn)行計(jì)算。

2 模型驗(yàn)證

2.1 射陽河口原型觀測

為了考察河口二維水流水質(zhì)耦合數(shù)學(xué)模型的可靠性、單元概化處理的合理性和模型參數(shù)選擇的適當(dāng)性，采用射陽河口 2006年8月24日 6:00 ～ 8 月25日 10:00（大潮）和 2007年1月3日 17:00 ～1月4日 19:00（大潮）的潮位、流向和平均流速等實(shí)測資料對(duì)模型進(jìn)行驗(yàn)證。利用平均絕對(duì)誤差、均方差和相對(duì)均方差，分析潮位、流速和流速的計(jì)算值與實(shí)測值的擬合程度。

本次計(jì)算范圍全長約 74 km，射陽河口以北25 km，以南 49 km，河口至外海處 48 km，計(jì)算區(qū)域面積約 4070 km2。前處理工作實(shí)現(xiàn)了網(wǎng)格的剖分，對(duì)計(jì)算水域采用無結(jié)構(gòu)任意四邊形網(wǎng)格布置，以貼合黃海天然岸線邊界。對(duì)河口附近區(qū)域采用網(wǎng)格局部加密技術(shù)。網(wǎng)格總數(shù)為 5943，節(jié)點(diǎn)總數(shù)為 6076，網(wǎng)格尺度變幅范圍為 350 m×500 m～900 m×1000 m。

根據(jù)網(wǎng)格單元的大小，設(shè)定時(shí)間步長為 Δ t = 10 s。

2.2 水動(dòng)力計(jì)算結(jié)果驗(yàn)證及分析

采用上述網(wǎng)格布置及邊界條件等，得出水動(dòng)力模擬計(jì)算結(jié)果。大潮潮位、流向和垂線流速計(jì)算值與實(shí)測值誤差統(tǒng)計(jì)表如表1 所示。

表1 潮位計(jì)算值與實(shí)測值誤差統(tǒng)計(jì)表

上述驗(yàn)證結(jié)果分析表明：射陽河口二維水流水質(zhì)耦合數(shù)學(xué)模型對(duì)區(qū)域的概化處理合理，模型選用的參數(shù)基本重現(xiàn)了射陽河口區(qū)域的潮位及潮流變化過程，計(jì)算結(jié)果可信。能夠較好地模擬射陽河口的水流運(yùn)動(dòng)，可在此基礎(chǔ)上模擬污染物的輸運(yùn)擴(kuò)散。

3 射陽河口污染物輸運(yùn)數(shù)值模擬

3.1 計(jì)算條件選取

3.2 濃度場模擬

計(jì)算步長：水質(zhì)模擬時(shí)取 10 s。

邊界條件：輸入污染物 CODMn排放量 823.31 g/s。

初始值的選?。簡卧獌?nèi)初始污染物濃度取為0 mg/L。

針對(duì)射陽河口實(shí)際情況，高潮位時(shí)不開閘，低潮位時(shí)開閘，在排放量保持一定的情況下，計(jì)算模擬 2006年8月 31日落潮不同時(shí)段污染物濃度場分布，確定何時(shí)開閘排放污染物對(duì)環(huán)境影響最小。設(shè)計(jì)工況條件表如表2 所示。

表2 工況條件

通過數(shù)值模擬結(jié)果分析可以得出在整個(gè)潮流周期中：

1）工況 1 污染物變化幅度為1.0～7.0 mg/L，4.0 mg/L 的污染帶面積在開閘 3 h 時(shí)最大，為 3.0 km2。由流場分布可知，退潮時(shí)，流場是由西北至東南向，污染物在潮流的作用下向東南方移動(dòng)。漲憩過后剛開始退潮時(shí)，水流流速較小，污染物擴(kuò)散速度偏慢，隨著時(shí)間的推移，水流流速逐漸變大，污染物擴(kuò)散速度也變快。由于前 3 h 內(nèi)污染物是連續(xù)不斷排放，此后關(guān)閘，因此開閘 3 h 時(shí)污染物濃度和影響范圍大。污染物不斷向東南方移動(dòng)，但污染物濃度和影響范圍逐漸變小。

2）工況 2 污染物變化幅度為1.0～6.0 mg/L，4.0 mg/L 的污染帶面積在開閘 3 h 時(shí)最大，為 2.1 km2。開閘的 3 h 和關(guān)閘后的前 2 h 還屬退潮時(shí)期，污染物在潮流的作用下向東南方移動(dòng)。關(guān)閘的第 3 小時(shí)，污染物不再繼續(xù)向東南方移動(dòng)。由潮流特征，可知落潮 3 h 后，水流流速較大，污染物擴(kuò)散速度也隨之變大。所以相同位置對(duì)應(yīng)時(shí)刻污染物的濃度較工況 1 小，且相應(yīng)污染物影響范圍比工況 1 小。當(dāng)閉閘不排污時(shí)，由于潮流和水體自凈功能的共同作用，之前 3 h 所排污染物濃度和影響范圍逐漸變小。

3）工況 3 污染物變化幅度為 5.0～12.0 mg/L。開閘 3 h 時(shí)濃度和影響范圍最大，9.0 mg/L 的污染帶面積為 3.8 km2。此時(shí)污染物排放的時(shí)間屬于落急至落憩點(diǎn)之間這段時(shí)間。污染物在橫、縱向的擴(kuò)散都加大，且濃度值相對(duì)較大，但關(guān)閘后緊接著漲潮，流向相反，污染物在潮流作用下向西北方移動(dòng)。先前排放擴(kuò)散的污染物回蕩、疊加而使污染物濃度變大，且污染物向河口聚積，不易擴(kuò)散，污染物濃度和影響范圍縮小的速度慢，因而對(duì)污染排放不利。

4）比較各個(gè)工況的最大污染帶面積，工況 2 中最大污染帶面積小，4.0 mg/L 的最大污染帶面積為 2.1 km2；工況 1 中 4.0 mg/L 的最大污染帶面積為3.0 km2，比工況 2 情況差，但還能為環(huán)境所接受，沒超標(biāo)；工況 3 中 9.0 mg/L 的污染帶面積為 3.8 km2，超標(biāo)。且由 3 種工況比較得出，在工況 2 中，閉閘后污染物濃度和影響范圍縮小最快，因此工況 2 對(duì)環(huán)境影響最小，在實(shí)際中，可按此種情況開閘排污。這為有效控制射陽河口近岸海域污染、科學(xué)制定污染物總量控制規(guī)劃及射陽閘調(diào)度運(yùn)行具有重要的應(yīng)用價(jià)值。

4 結(jié)語

本文采用基于有限體積法及黎曼近似解的二維水流水質(zhì)模型，模擬射陽河口的水動(dòng)力場和污染物濃度場分布。將潮位、流速和流向的計(jì)算值與實(shí)測值進(jìn)行比較，驗(yàn)證結(jié)果表明所建立射陽河口二維水量水質(zhì)耦合模型較好地重現(xiàn)了射陽河口區(qū)域的潮位及潮流變化過程，說明模型中的參數(shù)取值合理，計(jì)算結(jié)果可信，為進(jìn)一步研究該海域提供一定技術(shù)支持。

[1]Zhao D H, Shen H W, Tabios Ⅲ GO, et al. Finite-Volume Two-Dimensional Unsteady-Flow Model for River Basins[J].Hydr. Engrg, ASCE, 1994, 120 (7): 863-883.

[2]趙棣華，徐葆華. 平面二維水流—水質(zhì)有限體積法及黎曼近似解模型[J]. 水科學(xué)進(jìn)展，2000, 11(4): 368-373.

[3]Zhao D H, Shen H W, Lai J S, et al. Approximate Riemann Solvers in FVM for 2D Hydraulic Shock Waves Modeling[J].Hydr. Engrg, ASCE, 1996,122 (12): 693-702.