耦合高效高精度水動力模型的多組分污染物輸移及衰減反應模型

欒廣學，侯精明，楊 露，李丙堯，郭敏鵬，杜穎恩，馬 鑫

(西安理工大學省部共建西北旱區生態水利國家重點實驗室，陜西 西安 710048)

隨著我國城鎮化進程快速發展，各類水環境問題相繼浮現，尤其是突發污染事故頻繁發生[1-4]，不僅對城鄉居民的飲用水安全和生態環境造成了嚴重的影響，而且也影響了地區的經濟發展[5]。如2012年三友化工污染門事件、2012年廣西鎘污染事件和2014年蘭州自來水苯含量超標事件，對水市場、水環境產生了嚴重的影響，造成大量的生命財產損失。為有效應對各類水污染事故，及時預測污染事故嚴重程度并開展預報預警工作，對污染事故進行模擬顯得尤其重要[6-10]。

在洪水演進及其伴隨污染物輸移模擬研究中，污染物的衰減反應與模型的計算精度和效率成為研究重點[11-12]。以計算機為依托構建數值模型，精確掌握水體有機污染物的輸移及衰減反應規律，是當前研究自然水體水動力過程和水環境特征的重要手段[13-15]。目前，能夠模擬水環境問題的水質模型有很多，如SWAT、SWMM和MIKE21 FM等[16]。張京等[17]應用SWAT模型對義烏江流域的徑流、泥沙與水質過程進行了模擬。姚煥玫等[18]基于SWMM模型構建了南寧市區地表徑流及非點源污染精細化雨水徑流模型。李添雨等[19]采用MIKE21模型構建了沙河水庫庫區二維水動力水質模型，對沙河水庫水量水質變化情況進行了模擬。以上模型均可對河流水庫水質變化進行模擬，但模擬時存在輸入地形網格精度低、計算效率不高等缺點[20]。李一平等[21]對地表水環境數學模型進行了系統的分析。陶亞等[22]以深圳灣為例對突發污染物事故的輸移擴散規律和影響因素進行了模擬分析。龍巖等[23]采用HEC-RAS一維水動力水質模型對南水北調中線單組分水溶性突發污染事件的輸移擴散規律進行了模擬預測研究。司鵠等[24]對三峽庫區突發污染物的輸移及衰減規律進行了模擬分析。以上研究主要為單一污染物的輸移特性，且采用的大多為國外開發的商業軟件，無法根據實際模擬對象對源代碼進行修改調整。由此可見，能夠對多組分水污染事故輸移及衰減反應過程進行模擬預測的國內自主開發的模型還不夠成熟。

國內自主開發的水動力水質模型對流域水質模擬時必須解決的關鍵問題是提高水質模擬性能和功能。基于以上模型對污染物輸移過程模擬的不足，本文全面耦合基于圖形處理器(graphics processing unit， GPU)加速技術的二維水動力模型與多組分污染物輸移及衰減反應的Streeter-Phelps模型，建立一個可應用于復雜地形河道入流與潰壩洪水演進全過程與多組分伴隨污染物輸移及衰減反應的高性能耦合模型(GPU accelerated surface water flow and associated transport, GAST模型)，以期合理高效地模擬各類突發水污染事件并進行預警和評估。

1 數值模型及求解方法

1.1 二維水動力模型

GAST模型水動力模塊控制方程為二維淺水方程[25]，主要針對具有自由液面且以平面運動為主的水流，只考慮水平方向流速，忽略垂向運動。忽略了運動黏性項、紊流黏性項、風應力和科氏力的二維非線性淺水方程守恒形式為

(1)

其中

式中：t為時間，s；q為變量矢量；h為水深，m；qx、qy分別為x、y方向的單寬流量，m2/s；F、G分別為x、y方向的通量矢量；g為重力加速度，m/s2；u、v分別為x、y方向的流速，m/s；S是源項矢量；i為入滲源項；zb為河床底面高程；Cf為謝才系數，Cf=gn2/h1/3，其中n為曼寧系數。

1.2 多組分污染物輸移及衰減反應模型

污染物輸移控制方程為二維對流擴散方程：

(2)

式中：C為污染物垂線平均質量濃度，mg/L；Dx、Dy分別為x、y方向的擴散系數；qin為點源排放的流量強度，m/s；Cin為點源的物質垂線平均質量濃度，mg/L；R為轉化項，如污染物生化反應、生物作用等。

以BOD-DO耦合模型為例，僅考慮衰減與考慮生化反應情況時轉化項R可分別表示為

R=k1ρ(BOD)

(3)

R=k1ρ(BOD)-k2D

(4)

其中

D=ρs(DO)-ρ(DO)

式中：k1為生化需氧量BOD耗氧速率，s-1；k2為溶解氧DO復氧速率，s-1；ρ(BOD)為BOD質量濃度，mg/L；D為氧虧值，mg/L；ρ(DO)為DO質量濃度，mg/L；ρs(DO)為飽和溶解氧質量濃度，mg/L。

1.3 數值方法

本文模型對水動力及污染物控制方程進行離散所采用的方法為基于Godunov格式的有限體積法，該方法保證了守恒性且能有效解決潰壩激波等非連續問題。采用可自動滿足Godunov格式的HLLC近似黎曼求解器計算單元界面上的質量通量和動量通量以解決如沖擊波等的間斷問題[26]。通過靜水重構來實現干濕邊界處全穩條件，采用格式自適應方法來保障干濕交替過程模擬的穩定性。為達到復雜地形上的全穩條件，底坡源項使用底坡通量法處理復雜地形引起的動量不守恒問題。采用TVD-MUSCL方法進行數值重構，并采用兩步Runge-Kutta法進行時間步長的推進，以保證時間積分的二階精度[26]。離散網格類型為均勻結構網格。采用水深變化和水深值共同作為判別條件，以有效解決復雜地形干濕界面處的負水深和極端高流速等非物理現象所造成的計算失穩和物質動量的不守恒等問題。為提高模型的計算效率，通過CUDA語言自主編程實現GPU并行計算，以達到高速運算的目的[27]，圖1為GPU加速計算流程圖。

圖1 GPU加速計算流程Fig.1 GPU accelerated calculation flowchart

2 模型驗證

為驗證污染物輸移及衰減模型的計算精度，本文通過建立理想地形算例對污染物衰減的數值解與解析解進行對比分析，模型處理擴散項的合理性已得到充分驗證[28]，故模型驗證忽略污染物的擴散項，僅對模型的衰減過程進行驗證。構建一個邊長為50 m、底坡為0的正方形作為理想地形，以驗證模型在靜水狀態下的污染物衰減過程。

2.1 解析解

式(3)與(4)中反應項R可看作水質變量的質量濃度對時間的全導數，DO與BOD的反應項形式為

(5)

(6)

其邊界條件為：在x=0處，ρ(BOD)=ρ0(BOD)，ρ(DO)=ρ0(DO)，D=D0。因數值模擬時間步長很短，因此對上述反應項使用有限差分法求解可得一個步長內的解析解形式為

ρ(BOD)=ρ0(BOD)-k1ρ(BOD)dt

(7)

ρ(DO)=ρ0(DO)-k1ρ0(DO)dt+k2ρ0(DO)dt

(8)

式(7)為單組分衰減模型解析解，式(7)與(8)為多組分衰減反應模型解析解。

2.2 單組分衰減模型

單組分衰減模型驗證算例四周均為固壁邊界，點源污染物設置在x=25 m、y=25 m處，質量濃度為1 mg/L、初始水深為1 m、衰減系數取0.1，計算點源污染物質量濃度的衰減變化，計算時長為50 s。圖2為點源中心位置單組分污染物質量濃度在線性衰減情況下的變化趨勢。結果表明：在線性衰減情況下，污染物質量濃度隨時間的推移逐漸降低，因其與自身質量濃度有關，因此在污染物衰減速率不變時衰減量隨污染物質量濃度的減小而減小；單組分污染物質量濃度的衰減過程數值解與解析解基本吻合，無明顯數值震蕩現象，表明模型穩定性好且精度高。

圖2 單組分污染物質量濃度衰減變化過程Fig.2 Change process of attenuation of single-component pollutant concentration

2.3 多組分衰減反應模型

多組分衰減模型驗證算例采用BOD和DO兩種污染物組分耦合模型，忽略縱向離散作用，且不考慮飽和DO質量濃度對復氧速率的作用。多組分點源污染物位置設置在x=25 m、y=25 m處，將BOD和DO初始質量濃度均設置為1 mg/L，水深為1 m，BOD耗氧系數取0.05，DO復氧系數取0.02，計算多組分點源污染物質量濃度的衰減反應變化，計算時長為50 s。

靜水流場中多組分污染物衰減反應模擬結果對比與數值解及解析解部分計算結果如圖3所示。由圖3可見，在忽略對流項與擴散項的情況下，BOD和DO相互反應的數值解與解析解在模擬時間內的任何時刻相對差值均極小。計算得出，BOD數值解與分析值的相對誤差為1.8%；DO數值解與解析解的相對誤差為1.0%。顯然，模型模擬多組分污染物的數值模擬結果與理論值相比偏差較小，在允許范圍內，具有較高的穩定性和精度。模型實現了多組分污染物衰減反應功能，BOD的質量濃度逐漸降低且衰減速率逐漸變緩；DO質量濃度在26 s之前降低且衰減速率逐漸變小，在26 s達到轉折點，隨后其質量濃度開始回升，且回升速率逐漸增加。本算例BOD質量濃度的衰減只與其自身質量濃度有關，不受DO質量濃度影響，且為線性衰減，故其質量濃度逐漸降低且總量衰減速度逐漸變緩；DO質量濃度的變化受BOD質量濃度的影響且其復氧速率與其自身濃度有關，BOD質量濃度較大時，其耗氧速率較DO的復氧速率大，故DO質量濃度降低，因BOD耗氧速率比其自身的衰減速率小，因此隨著BOD質量濃度的減小，其耗氧速率逐漸接近DO的復氧速率，并且在26 s時刻相等，DO質量濃度到達最低值，隨后BOD質量濃度繼續降低，耗氧速率越來越微弱，故復氧速率逐漸占主導地位，DO質量濃度開始回升且回升速率越來越快。

圖3 靜水流場中多組分污染物相互反應的數值解與解析解對比Fig.3 Comparison of numerical and analytical solutions of multi-component pollutant interaction in static water flow field

3 Toce河道多組分點源輸移及衰減過程模擬

對Toce河道上釋放的多組分點源污染物輸移及衰減過程進行模擬，河道DEM數據分辨率為0.05 m，網格單元共計206 640個。河道數字地形高程如圖4所示，左邊界為入流口，入流寬度設置為3.4 m，右邊界為自由出流的開邊界，其余為閉邊界。入流流量過程線如圖5所示，河道曼寧系數取0.016 2 s/m1/3，初始水深設為0 m，點源位置設置在x=7.868 m、y=5.882 m處，如圖4所示。點源同時釋放兩種污染物，用C1和C2表示，排放質量濃度為1 mg/L，強度為0.8 m/s，設置C1衰減系數為 0.3 s-1，C2無衰減，污染物隨水流向下游推進。

圖4 Toce河道數字地形高程及點源位置Fig.4 Digital terrain elevation and point source location of the Toce River

圖5 入流流量過程線Fig.5 Inflow hydrograph

采用本文模型模擬計算3 min內的水流演進及其伴隨污染物輸移及衰減過程。在t=20 s和t=50 s 時的水深模擬結果如圖6所示，污染物有衰減與無衰減輸移在t=20 s，40 s，60 s、80 s時的質量濃度分布情況如圖7所示，多組分污染物在x=11.125 m、y=5.825 m處的有衰減與無衰減過程質量濃度對比如圖8所示。

(a) t=20 s

由圖6可知水流向下游推進速度很快。圖7表明，在t=20 s時水流即流過點源位置，污染物隨著水流向下游推進，質量濃度逐漸降低。在t=40 s時，污染物隨水流向周圍低洼處輸移；在t=60 s時，部分水流進入旁側低洼區，污染物隨水流開始分流；在t=80 s時，污染物即將隨水流從河道下游邊界流出，污染物C1隨水流遷移過程中因自身的衰減，質量濃度明顯低于無衰減情況，C1到達下游出口位置時間明顯較C2延后。圖7也顯示了污染物在點源位置釋放并向下游輸移的整個過程，污染物分布與干濕界面相匹配，在復雜地形的河道內始終隨水流向下游輸移，由此表明：點源污染物排放不僅會造成嚴重的水體污染，而且對下游水體的影響很大。圖8表明，在x=11.125 m、y=5.825 m處，C1的質量濃度始終低于C2的質量濃度，由于復雜地形河道的影響，C1與C2的質量濃度差值呈不規則趨勢。

(a) 無衰減，t=20 s

圖8 多組分污染物有衰減與無衰減輸移過程質量濃度對比Fig.8 Comparison of concentrations of multi-component pollutants with and without attenuation

4 Malpasset潰壩多組分點源輸移及衰減過程模擬

本節模擬計算Malpasset河道潰壩水流運動及其伴隨污染物的輸移、衰減過程。河道數字地形高程如圖9所示，網格精度為10 m，網格單元共計 1 581 714 個，壩前水位設置為100 m，多組分點源污染物排放口設置于x=5.257 km、y=6.685 km處，如圖9所示，排放質量濃度為1 mg/L，強度為1 m/s，河道曼寧系數取0.016 2 s/m1/3，點源同時釋放兩種污染物，用C3和C4表示，點源排放質量濃度為 1 mg/L，強度為0.8 m/s，設置C3衰減系數為 0.000 1 s-1，C4無衰減。圖10為t=300 s和t=1 800 s時的水深分布，圖11為污染物無衰減與衰減輸移(t=300 s，1 500 s，3 000 s)質量濃度分布。

圖9 Malpasset河道數字地形高程及點源位置Fig.9 Digital terrain elevation and location of point source of the Malpasset River

由圖10可見，潰壩發生后，水流迅速向下游推進。由圖11可見，多組分污染物隨潰壩水流向下游推進，1 500 s時即到達下游位置。河道下游C3的質量濃度明顯小于C4的質量濃度，由此表明污染物因其自身的衰減反應，隨水流到達下游時，質量濃度削減效果明顯。本文模型在GTX1080計算機上運行0.37 h即可完成2 h水流及污染物輸移、衰減反應過程模擬。

(a) t=300 s

(a) 無衰減，t=300 s

5 結 語

本文將基于GPU加速技術的二維水動力模型與多組分污染物輸移及衰減反應模型進行耦合，構建了一個可應用于復雜地形河道入流與潰壩洪水演進全過程與多組分污染物輸移及衰減反應的高性能耦合模型。模型實現多組分污染物輸移、衰減反應模擬功能，適用于大區域復雜地形河道的水質模擬；引入GPU加速技術提高了模擬效率，計算1 581 714個均勻結構網格2 h的潰壩及污染物輸移、衰減反應模擬，僅需運行0.37 h。

通過對理想地形及兩個復雜地形河道的水質模擬，結果表明，模型計算結果與解析解基本吻合，在入流與潰壩算例中的模擬結果符合實際的物理過程，運行速度快且計算效率高，具有較好的魯棒性，能高效地模擬各類水污染事件并進行預警和評估，可為敏感水域水環境的治理和保護提供技術與數據支撐。