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

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

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

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

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

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

1 數(shù)值模型及求解方法

1.1 二維水動力模型

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

(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為謝才系數(shù)，Cf=gn2/h1/3，其中n為曼寧系數(shù)。

1.2 多組分污染物輸移及衰減反應(yīng)模型

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

(2)

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

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

R=k1ρ(BOD)

(3)

R=k1ρ(BOD)-k2D

(4)

其中

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

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

1.3 數(shù)值方法

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

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

2 模型驗證

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

2.1 解析解

式(3)與(4)中反應(yīng)項R可看作水質(zhì)變量的質(zhì)量濃度對時間的全導(dǎo)數(shù)，DO與BOD的反應(yīng)項形式為

(5)

(6)

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

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

(7)

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

(8)

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

2.2 單組分衰減模型

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

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

2.3 多組分衰減反應(yīng)模型

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

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

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

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

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

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

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

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

(a) t=20 s

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

(a) 無衰減，t=20 s

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

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

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

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

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

(a) t=300 s

(a) 無衰減，t=300 s

5 結(jié) 語

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

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