基于自適應網格的急流條件下污染物輸運高效高精度模擬

曹 引，冶運濤，梁犁麗，趙紅莉，蔣云鐘，王 浩

（1.東華大學 環境科學與工程學院，上海 201620；2.中國水利水電科學研究院 水資源研究所，北京 100038）

1 研究背景

近年來，突發水污染事件頻發，導致局部水域水質惡化，嚴重威脅周圍及下游區域居民的用水安全，同時會造成巨大的經濟損失和嚴重的生態環境問題［1］。突發水污染事件的影響范圍和水流條件密切相關，如果在暴雨洪水和潰壩水流等急流條件下突發水污染事件，釋放的污染物將隨著急速演進的洪水快速運移和擴散，導致水質在短時間內迅速惡化，可能造成自來水廠停水和工業停產等問題，威脅下游居民的用水安全，極易引發社會恐慌［2］。

水污染事件突發后，需要了解污染物在水體中的運移和擴散規律，及時判斷突發水污染事件的危害范圍和持續時間，以助于快速評估突發水污染事件的危害并及時預警［3-4］。水動力水質模型是評估突發水污染事件危害程度的有效手段，在突發水污染事件管理中得到有效應用［5］。陶亞等［6-7］分別采用基于有限差分方法的平面二維水動力水質模型和EFDC水動力及污染物模型預測了河流突發水污染事故后下游城市的應急響應時間，模擬了不同應對措施的處理效果；饒清華等［8］基于有限元法構建了模擬閩江下游突發水污染事件中污染物遷移擴散的水動力及水質模型，定量模擬了污染物的時空分布；Dong等［9］將基于有限體積法的MIKE21模型用于突發水污染事件的風險分析，模擬了沿海河流和近岸海域化學污染物輸移擴散過程。目前研究多集中于緩流條件下突發水污染事件中污染物的輸運模擬，準確模擬暴雨洪水和潰壩洪水等急流條件下污染物的輸運規律充滿了挑戰［5］:（1）突發水污染事件的地點和時間的隨機性，需要快速獲取水污染事件發生區域的地形、初始條件和出入流邊界條件，為水質模型提供建模數據；（2）洪水演進過程中會產生激波和動態變化的干濕邊界，模型需要能夠有效捕捉激波和模擬干濕邊界；（3）實際應用中水流常流經復雜地形區域，如山谷、河流和城市區域等，模型必須能夠有效描述復雜的地形條件；（4）突發水污染事件危害評估和預警要求很高的時效性，模型必須具有很高的計算效率。水質模型構建所需的地形條件可以通過SRTM 30米DEM數據、歷史地形觀測數據、人工測量和插值等方式快速獲取；初始條件和出入流邊界可以通過上下游水文站、水質在線監測站或者便攜式水質分析儀測量獲取。水質模型構建后，求解模型控制方程獲得穩定和諧的數值解是應用水動力水質模型模擬污染物輸運過程的關鍵，相比于有限元法和有限差分方法，Godunov有限體積法由于可以有效捕捉急流條件下突發水污染事件中產生的激波間斷和污染物濃度間斷［10］，適合急流條件下污染物的輸運模擬［11］。水動力水質模型利用網格離散計算區域，結構網格和非結構網格是水動力水質模型最常用的網格類型。非結構網格對復雜地形和邊界具有較強的擬合能力，但生成過程復雜［12］；結構網格生成簡單，但對復雜地形的擬合能力較差，為了提升模型模擬精度，必須增加網格數量，這勢必會降低模型計算效率。基于結構網格的自適應網格技術可以根據水流特性和污染物濃度分布自適應調整網格大小，保證模型模擬精度的前提下可以明顯提升模型計算效率［5，13］。

本文基于自適應結構網格構建了適用于復雜地形和急流條件下突發水污染事件中污染物輸運模擬的二維水流-輸運耦合模型，模型采用能夠有效捕捉激波的Godunov有限體積法離散二維水流-輸運控制方程，利用具有時空二階精度的MUSCL-Hancock方法求解，界面通量采用HLLC近似黎曼求解器求解，為了保證干濕邊界處模型的穩定性和守恒性，采用非負重構技術［11］和局部底部高程修正方法［13］重構界面兩側黎曼變量。最后分別利用水槽試驗、物理模型和實際算例檢驗了模型模擬突發水污染事件中污染物輸移的精度和穩定性。

2 自適應網格生成技術

本文采用的自適應網格為遞歸的層次結構網格（圖1），葉網格（如網格（i，j）和（i，j+1））劃分層級lev為0，網格最大劃分層級為div_max，葉網格大小為Δx0和Δy0。對葉網格進行細化可以得到子網格，子網格大小為Δx=Δx0/2lev，Δy=Δy0/2lev。任何網格和其鄰居網格大小必須滿足2倍關系，即網格大小必須是其鄰居網格大小的2倍、1倍和1/2倍（2∶1準則），該過程需要存儲鄰居網格信息，為了減少存儲要求，鄰居網格信息采用網格之間的相對關系來確定［14］。自適應網格可以快速生成，具體生成步驟如下:

（1）生成能覆蓋計算域的葉網格（i，j），i=1，2，…，M，j=1，2，…，N，M和N分別為x和y方向上的網格數，設置葉網格劃分層級為0，即lev（i，j）=0；

（2）根據研究區邊界種子點識別參與計算的網格單元；

（3）依次計算內部網格單元（i，j，is，js）（圖1）的水位梯度（gradη）和污染物守恒濃度梯度（gradqc）（式（1）—式（4）），取兩者最大值gradΦ=max（gradη，gradqc）作為該網格劃分層級的判定梯度；

（4）如果gradΦ大于網格劃分閾值Φsub，或者網格（i，j，is，js）位于干濕邊界處，設置葉網格的劃分層級為lev（i，j）=lev（i，j）+1，如果lev（i，j）=div_max，則保持lev（i，j）=div_max；

（5）采用三角形線性插值獲取子網格中心高程，然后根據質量守恒定律和動量守恒定律確定子網格水位、流量和污染物濃度［15］；

（6）判斷并調整子網格（i，j，is，js）大小和其8個鄰居網格大小之間滿足2∶1準則［13］；

（7）如果gradΦ，is=1，…，Ms，js=1，…，Ms，Ms=2lev（i，j）均小于網格合并閾值Φcoar，則葉網格的劃分層級設置為lev（i，j）=lev（i，j）-1，如果lev（i，j）=0，則保持lev（i，j）=0，合并得到的母網格中心狀態變量為4個子網格中心狀態變量的平均值，該過程同樣保證網格（i，j，is，js）大小和其8個鄰居網格大小滿足2∶1準則。

式中:ηeast、ηwest、ηnorth、ηsouth和qceast、qcwest、qcnorth、qcsouth分別為單元（i，j，is，js）東、西、北、南4個方向的水位和污染物守恒濃度，當網格（i，j，is，js）的劃分層級和其鄰居網格的劃分層級相同時，則網格（i，j，is，js）鄰居方向的水位和污染物守恒濃度直接取鄰居網格中心的水位和污染物守恒濃度，否則通過自然鄰近插值獲取該方向的水位和污染物守恒濃度［16-17］。

圖1 自適應網格示意圖

自適應網格生成技術常采用水位、水深、流速或者污染物質量、濃度梯度等設置網格劃分閾值，閾值設置存在兩種形式:相對閾值［5］和絕對閾值［13］。相對閾值通過對所有計算網格梯度grad進行降序排列，取某個分位數Sa處的梯度作為閾值，相對閾值在計算過程中會不斷變化；相對閾值設置簡單，但當多數網格均存在較大梯度時，基于相對閾值的自適應網格技術只能識別和細化部分（計算網格總數×Sa）大梯度網格，影響模型模擬精度；當水流趨于平穩，大部分網格梯度均較小時，基于相對閾值的自適應網格技術仍會對部分網格（計算網格總數×Sa）進行細化，影響計算效率。絕對閾值采用固定梯度作為閾值，對梯度大于該閾值的網格進行細化，直至達到最大劃分層級；當多數網格均存在較大梯度時，基于絕對閾值的自適應網格技術能夠更準確地識別并細化大梯度網格，當多數網格梯度均較小時，低于絕對閾值的子網格將會被識別和合并。相比于相對閾值，絕對閾值可以更好地識別和細化大梯度網格、合并小梯度網格，提高模型模擬精度和計算效率。

3 控制方程

在滿足靜水壓力和忽略水體垂向加速條件下，利用水深平均的三維Navier-Stokes方程可以推導得到二維淺水方程［13，18］，加上描述物質輸移的對流-擴散方程［19］，二維水流-輸運方程形式如下:

式中:t為時間，x和y分別為水平橫向和縱向坐標；U為守恒變量；E為對流項通量，E=（F，G）；F、G分別代表x、y方向的對流通量；S為源項。

其中U、F、G和S表示如下:

式中:η、h和zb分別為水位、水深和河底高程，η=zb+h；u和v分別為x和y方向流速；g為重力加速度；ρ為水密度；Dx和Dy表示x和y方向的擴散系數，c為物質的垂線平均濃度，qin和cin分別表示點源的流量強度和物質垂線平均濃度；τbx和τby為x和y方向的床面摩擦應力，表示因床面摩擦導致的水流能量耗散（式（7））。

式中:Cf為床面摩擦系數，Cf=gn2/h1/3，n為Manning糙率系數，s/m1/3。在涉及干濕邊界計算的案例中，設置一個最小水深，當單元水深小于最小水深時，Cf取0。

4 數值求解

采用有限體積法離散控制方程，在任意控制體Ω上對式（5）進行積分:

通過高斯-格林公式將∫Ω??E dΩ轉化為沿控制體邊界的面積分，得到式（9）:

式中:U為單元均值；S為控制體Ω的邊界；n為垂直于面元dS的單位向量。在笛卡爾坐標系下，式（9）中的曲面積分項可以轉換為:

式中:Δx和Δy分別為笛卡爾坐標系中網格單元（i，j）的左右和上下邊長，F（i+1/2，j）、F（i-1/2，j）、F（i，j+1/2）和F（i，j-1/2）分別為單元東、西、南和北4個界面的界面通量（圖2）。

采用自適應網格時，網格單元和鄰居網格單元可能處于不同的劃分層級，當網格單元右側的鄰居網格劃分層級大于該網格單元的劃分層級時，此時通過東邊界面的通量的計算需要w1、e1、w2和e2位置的狀態變量，e1和e2位置狀態變量可以直接獲取，w1和w2位置狀態變量可以通過自然鄰近插值獲取（圖2）［16-17］。

由式（9）—式（10）可以得到由n時層到n+1時層推進的守恒有限體積公式:

為了使模型數值解具有時空二階精度，采用MUSCL-Hancock方法求解控制方程，MUSCL-Hancock方法包括預測步和校正步兩步。

圖2 自適應網格通量計算示意圖

預測步中，利用n時層單元（i，j）和4個鄰近單元中心狀態變量U重構單元（i，j）4個界面左右的狀態變量，以界面（i+1/2，j）為例，界面左側狀態變量UL（i+1/2，j）和右側狀態變量UR（i+1/2，j）重構如下:

式中:φ()r為坡度限制器，限制器根據β（1≤β≤2）取值不同分為Roe’s minmod限制器（β=1）和Roe’s superbee限制器（β=2）。本文采用Roe’s minmod限制器。

基于重構后的界面左右狀態變量，計算界面數值通量F*，然后代入式（14）計算n+1/2時層單元（i，j）中心狀態變量值

校正步中，基于預測步獲取的n+1/2時層單元（i，j）和4個鄰近單元中心狀態變量12采用非負線性重構技術［11］和局部底部高程修正方法［13］對界面左右狀態變量進行空間重構，利用HLLC近似黎曼求解器［20］計算界面通量F*n+12，然后代入式（15）計算n+1時層單元（i，j）中心狀態變量值

5 結果與討論

為了檢驗模型模擬急流條件下污染物輸移的能力，分別利用水槽試驗、物理模型和實際案例等經典算例檢驗模型模擬精度、穩定性以及模擬效率。3個算例中網格細化和粗化的絕對閾值分別設置為0.08（Φsub）和0.05（Φcoar）。所有算例中重力加速度g取9.81 m/s2，水體密度ρ取1000 kg/m3。

5.1 均勻濃度的水槽潰壩水流-輸運模擬本算例［20］常用于檢驗模型的計算精度以及動態變化干濕邊界處理的有效性。計算域為75 m×30 m，大壩位于x=16 m，忽略大壩厚度。糙率n=0.018 s/m1/3。底部高程計算公式為:

大壩上游初始水位和物質濃度分別取1.875 m和1 mg/L，流速設為0；下游水深為0，即干河床；模擬污染物輸移過程中不考慮物質降解。x和y方向的擴散系數Dx和Dy均取0.5 m2/s。假設t=0 s時大壩瞬時全潰。由于初始條件為均勻濃度水體，因此，水體的物質濃度不隨時間變化，始終為1 mg/L。

初始條件（t=0 s）及不同時刻（t=10、20、30 s）的流場、物質濃度計算結果如圖3所示，不同時刻的網格分布和網格數如圖4和圖5所示。由圖3可知，潰壩水流演進過程中水體的物質濃度始終保持為1 mg/L，模型模擬準確。由圖4可以看出，網格劃分層級隨著潰壩水流演進不斷的變化，由于地形和初始條件的對稱性，水體流速和水質狀態時刻保持對稱，因此自適應生成的網格同樣具有對稱性；由圖5可以看出t=14 s時網格最多，網格數為18 411，隨后網格數逐漸降低，最后網格數穩定在11 226附近。自適應網格技術可以自動捕捉高水位梯度和高污染物濃度梯度所在區域以及干濕邊界區域，對這些區域的網格進行細化，提升模型對這些區域水流和物質輸移的捕捉能力。此外，自適應網格技術可以提高模型模擬效率，t=300 s時基于自適應網格的模型模擬時間為1055.8 s，如果采用自適應網格中最大劃分層級的網格作為固定網格，模型模擬時間為2370.5 s，自適應網格技術可以節約55.5%的計算時間。

圖3 不同時刻潰壩水流流場和污染物濃度模擬結果

圖4 不同時刻網格分布

5.2 Toce河潰壩物理模型案例Toce河模型［21］是意大利國家電力公司ENEL按1∶100的比例尺對意大利米蘭市Toce河上游約5 km范圍內建立的物理模型，該模型常用于驗證各種潰壩數學模型的精度和穩定性。該物理模型長約50 m，寬約11 m。模型提供空間分辨率為5 cm的高程點，比較精確地描述了Toce河的真實地形。Toce河物理模型的地形、觀測點位置以及出入流邊界位置如圖6所示，Toce河中部有一個空水庫，水庫在靠近河流一側有個開口。Toce物理模型入流流量過程如圖7所示。模型初始水深為0 m，出口為自由出流邊界（圖6）；糙率取ENEL推薦的0.0162 s/m1/3；模型模擬初始網格大小為0.4 m×0.4 m，最大細分水平設置為2。假設從潰壩發生后第25秒開始，位于［7.868，5.882］處的一個點源開始釋放污染物，污染物釋放速度為1 m/s，污染物濃度為1 kg/m3，污染物的擴散系數為0.01 m2/s。

圖5 自適應網格數隨時間變化圖

圖6 Toce河物理模型地形

圖7 入流流量過程

圖8 水位實測值和模擬值

圖9 t=30s和60s模擬水深分布

圖10 t=30s和60s模擬水質分布

圖11 t=30s和60s模型網格分布

圖8為t=0～180 s模擬時間內P4、P9和P21點水位模擬值和觀測值對比圖，3個觀測點處模擬水位和觀測水位基本一致，模型具有很高的模擬精度。圖9為t=30 s和60 s的模擬水深分布，可以看出隨著t=18 s潰壩開始（圖7），潰壩洪水流速很快，快速向下游演進。圖10為t=30 s和60 s模擬污染物濃度分布，可以看出，急流條件下如果突發水

污染事故，污染物會隨著水流迅速向下游擴散。圖11和圖12為t=30 s和60 s的網格分布以及模擬過程中網格數的變化情況，潰壩開始前，計算網格為2461，t=18 s潰壩開始后，潰壩水流向下游快速演進，洪水經過區域具有較高的水位梯度，加上點源污染物釋放后的運移產生的污染物濃度梯度，導致網格數快速增加，t=78 s時達到最大值16 930個，隨后逐漸降低，穩定在15 300附近。t=180 s時基于自適應網格的模型模擬時間為365.3 s，相較于基于固定網格的模型模擬時間（435.2 s），節省了16.1%的模擬時間。由于潰壩水流經過的大部分區域均具有較高的水位或污染物濃度梯度，因此，需要相對較細的網格捕捉水流和污染物輸運特性，自適應網格可以自動捕捉具有高水位梯度和高污染物濃度梯度以及干濕邊界區域，細化該區域網格，保證該區域模型模擬精度。

圖12 t=0～180s自適應網格數隨模擬時間變化圖

5.3 Malpasset潰壩實際案例法國Malpasset大壩建于1956年，最大庫容為5500萬m3，主要用于農業灌溉和城市用水。由于1959年11月底連降暴雨，Malpasset大壩水位短時間內迅速上升，大壩最終于12月2日21時14分潰決，潰壩歷時很短，可認為是瞬時潰壩，Malpasset水庫及下游地形如圖13所示。潰壩水流迅速向下游演進，最終進入大海，海面高程為0 m；潰壩前，大壩上游水位為100 m，下游水深為0 m。假設潰壩900 s后位于（9660 m，3074 m）的位置突發污染物排放，排放的污染物濃度為1 kg/m3，排放強度為1 m/s，污染物擴散系數為0.5 m2/s。利用基于自適應網格的二維水流-輸運耦合模型模擬潰壩水流演進和污染物遷移擴散，初始網格大小為80 m×80 m，最大細分水位設為2，Manning糙率系數取0.033 s/m1/3。

圖13 Malpasset水庫及下游地形

圖14 Malpasset潰壩洪水演進過程

圖15 Malpasset潰壩水流污染物運移和擴散過程

圖16 不同時刻網格分布

模型模擬潰壩后1000 s和1800 s的水深和污染物濃度分布分別如圖14和圖15所示。潰壩后洪水向下游快速演進，t=1800 s時水流已經到達下游平原，和Hou等［22］模擬結果一致。t=900 s時水流已經經過（9660 m，3074 m），隨著污染物的釋放，污染物隨水流往下游遷移，同時由于濃度梯度的存在污染物逐漸向周圍擴散，污染物濃度沿著水流演進方向逐漸降低。圖16為1000 s和1800 s模型網格分布，自適應網格技術可以識別水位和污染物濃度梯度，自適應調整網格大小，同時可以捕捉干濕邊界變化，在干濕邊界處對網格進行細化，提高水位梯度和污染物濃度梯度較大區域以及干濕邊界處數值模擬精度。t=3000 s時基于自適應網格的模型模擬時間為860.4 s，相較于基于固定網格的模型模擬時間（1515.1 s），節省了43.2%的模擬時間。

6 結論

本文基于自適應網格技術構建了一種適用于復雜地形和急流條件下突發水污染事件中污染物輸運模擬的二維水流-輸運耦合模型，分別利用水槽試驗、物理模型以及實際算例驗證了模型精度和穩定性。驗證結果表明:（1）基于自適應網格技術的二維水流-輸運耦合模型可以捕捉高水位梯度、高物質濃度梯度和干濕邊界區域，并對這些區域網格進行細化，模擬精度高；（2）自適應網格技術在保證模型模擬精度的同時，可以提升模擬效率；（3）基于自適應網格技術的二維水流-輸運耦合模型能夠準確高效地模擬急流條件下污染物的輸運過程，適用于突發水污染事件的評估、預警和應急管理，具有推廣應用價值。未來將在模型中考慮COD、BOD、TN和TP等水質指標，研究突發水污染事件中多種目標污染物的輸運規律。