基于自適應(yīng)網(wǎng)格的急流條件下污染物輸運(yùn)高效高精度模擬

曹 引，冶運(yùn)濤，梁犁麗，趙紅莉，蔣云鐘，王 浩

（1.東華大學(xué) 環(huán)境科學(xué)與工程學(xué)院，上海 201620；2.中國(guó)水利水電科學(xué)研究院 水資源研究所，北京 100038）

1 研究背景

近年來(lái)，突發(fā)水污染事件頻發(fā)，導(dǎo)致局部水域水質(zhì)惡化，嚴(yán)重威脅周圍及下游區(qū)域居民的用水安全，同時(shí)會(huì)造成巨大的經(jīng)濟(jì)損失和嚴(yán)重的生態(tài)環(huán)境問(wèn)題［1］。突發(fā)水污染事件的影響范圍和水流條件密切相關(guān)，如果在暴雨洪水和潰壩水流等急流條件下突發(fā)水污染事件，釋放的污染物將隨著急速演進(jìn)的洪水快速運(yùn)移和擴(kuò)散，導(dǎo)致水質(zhì)在短時(shí)間內(nèi)迅速惡化，可能造成自來(lái)水廠停水和工業(yè)停產(chǎn)等問(wèn)題，威脅下游居民的用水安全，極易引發(fā)社會(huì)恐慌［2］。

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

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

2 自適應(yīng)網(wǎng)格生成技術(shù)

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

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

（2）根據(jù)研究區(qū)邊界種子點(diǎn)識(shí)別參與計(jì)算的網(wǎng)格單元；

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

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

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

（6）判斷并調(diào)整子網(wǎng)格（i，j，is，js）大小和其8個(gè)鄰居網(wǎng)格大小之間滿足2∶1準(zhǔn)則［13］；

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

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

圖1 自適應(yīng)網(wǎng)格示意圖

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

3 控制方程

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

式中:t為時(shí)間，x和y分別為水平橫向和縱向坐標(biāo)；U為守恒變量；E為對(duì)流項(xiàng)通量，E=（F，G）；F、G分別代表x、y方向的對(duì)流通量；S為源項(xiàng)。

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

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

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

4 數(shù)值求解

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

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

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

式中:Δx和Δy分別為笛卡爾坐標(biāo)系中網(wǎng)格單元（i，j）的左右和上下邊長(zhǎng)，F(xiàn)（i+1/2，j）、F（i-1/2，j）、F（i，j+1/2）和F（i，j-1/2）分別為單元東、西、南和北4個(gè)界面的界面通量（圖2）。

采用自適應(yīng)網(wǎng)格時(shí)，網(wǎng)格單元和鄰居網(wǎng)格單元可能處于不同的劃分層級(jí)，當(dāng)網(wǎng)格單元右側(cè)的鄰居網(wǎng)格劃分層級(jí)大于該網(wǎng)格單元的劃分層級(jí)時(shí)，此時(shí)通過(guò)東邊界面的通量的計(jì)算需要w1、e1、w2和e2位置的狀態(tài)變量，e1和e2位置狀態(tài)變量可以直接獲取，w1和w2位置狀態(tài)變量可以通過(guò)自然鄰近插值獲?。▓D2）［16-17］。

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

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

圖2 自適應(yīng)網(wǎng)格通量計(jì)算示意圖

預(yù)測(cè)步中，利用n時(shí)層單元（i，j）和4個(gè)鄰近單元中心狀態(tài)變量U重構(gòu)單元（i，j）4個(gè)界面左右的狀態(tài)變量，以界面（i+1/2，j）為例，界面左側(cè)狀態(tài)變量UL（i+1/2，j）和右側(cè)狀態(tài)變量UR（i+1/2，j）重構(gòu)如下:

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

基于重構(gòu)后的界面左右狀態(tài)變量，計(jì)算界面數(shù)值通量F*，然后代入式（14）計(jì)算n+1/2時(shí)層單元（i，j）中心狀態(tài)變量值

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

5 結(jié)果與討論

為了檢驗(yàn)?zāi)Ｐ湍M急流條件下污染物輸移的能力，分別利用水槽試驗(yàn)、物理模型和實(shí)際案例等經(jīng)典算例檢驗(yàn)?zāi)Ｐ湍M精度、穩(wěn)定性以及模擬效率。3個(gè)算例中網(wǎng)格細(xì)化和粗化的絕對(duì)閾值分別設(shè)置為0.08（Φsub）和0.05（Φcoar）。所有算例中重力加速度g取9.81 m/s2，水體密度ρ取1000 kg/m3。

5.1 均勻濃度的水槽潰壩水流-輸運(yùn)模擬本算例［20］常用于檢驗(yàn)?zāi)Ｐ偷挠?jì)算精度以及動(dòng)態(tài)變化干濕邊界處理的有效性。計(jì)算域?yàn)?5 m×30 m，大壩位于x=16 m，忽略大壩厚度。糙率n=0.018 s/m1/3。底部高程計(jì)算公式為:

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

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

圖3 不同時(shí)刻潰壩水流流場(chǎng)和污染物濃度模擬結(jié)果

圖4 不同時(shí)刻網(wǎng)格分布

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

圖5 自適應(yīng)網(wǎng)格數(shù)隨時(shí)間變化圖

圖6 Toce河物理模型地形

圖7 入流流量過(guò)程

圖8 水位實(shí)測(cè)值和模擬值

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

圖10 t=30s和60s模擬水質(zhì)分布

圖11 t=30s和60s模型網(wǎng)格分布

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

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

圖12 t=0～180s自適應(yīng)網(wǎng)格數(shù)隨模擬時(shí)間變化圖

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

圖13 Malpasset水庫(kù)及下游地形

圖14 Malpasset潰壩洪水演進(jìn)過(guò)程

圖15 Malpasset潰壩水流污染物運(yùn)移和擴(kuò)散過(guò)程

圖16 不同時(shí)刻網(wǎng)格分布

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

6 結(jié)論

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