基于SPH-SWE方法對(duì)河流保守型污染物擴(kuò)散的模擬

田麗蓉，余金玲，吳玉帥，顧聲龍

(青海大學(xué)水利電力學(xué)院，青海 西寧 810016)

光滑粒子流體動(dòng)力學(xué)(Smoothed Particle Hydrodynamics，SPH)是一種基于拉格朗日描述下的無(wú)網(wǎng)格方法，將連續(xù)介質(zhì)離散成為許多粒子，通過(guò)追蹤離散粒子的運(yùn)動(dòng)軌跡來(lái)描述連續(xù)介質(zhì)的運(yùn)動(dòng)過(guò)程。SPH方法最初由Lucy等[1]和Gingold等[2]在研究天體物理問(wèn)題中提出的一種方法。它已經(jīng)成功應(yīng)用到水動(dòng)力學(xué)領(lǐng)域，而且是解決水力學(xué)領(lǐng)域諸多問(wèn)題較成熟的數(shù)值模擬工具。SPH方法在處理大變形，自由表面流[3]和多相流[4]等問(wèn)題上具有明顯的數(shù)值優(yōu)勢(shì)。近年來(lái)，基于淺水波方程的SPH方法也得到了很大的發(fā)展[5]。Cleary等[6]將SPH方法應(yīng)用于熱傳導(dǎo)問(wèn)題，后續(xù)研究大多是基于此文獻(xiàn)提出的SPH擴(kuò)散模型來(lái)模擬，研究不同環(huán)境下污染物的擴(kuò)散。Zhu等[7]研究了溶質(zhì)在多孔介質(zhì)中的傳輸。饒登宇等[8]建立描述孔隙水流動(dòng)的SPH水動(dòng)力模型和描述溶質(zhì)分子擴(kuò)散的擴(kuò)散模型，得出彌散度與流速變異系數(shù)、迂曲度、迂曲路徑差以及不均勻系數(shù)大致呈正相關(guān)，與孔隙率呈負(fù)相關(guān)。沈嘉淵[9]則使用SPH方法和一階精度的歐拉向前方法模擬污染物擴(kuò)散問(wèn)題，通過(guò)不同時(shí)間積分方法的對(duì)比分析，得出精度最高的時(shí)間積分法。然而，在淺水環(huán)境下模擬污染物擴(kuò)散的研究較少，采用SPH-SWE(Shallow Water Equation，SWE)方法研究污染物擴(kuò)散，利用無(wú)網(wǎng)格性質(zhì)求解淺水中污染物的輸運(yùn)過(guò)程比網(wǎng)格法求解SWE更具有優(yōu)勢(shì)[10]。由于污染物輸運(yùn)方程的二階性質(zhì)，傳統(tǒng)的網(wǎng)格方法在模擬過(guò)程中存在網(wǎng)格畸形等問(wèn)題。無(wú)網(wǎng)格方法中SPH法的發(fā)展迅速，在流體力學(xué)領(lǐng)域具有顯著數(shù)值優(yōu)勢(shì)。它不依賴網(wǎng)格，在計(jì)算中以質(zhì)點(diǎn)描繪整個(gè)模擬區(qū)域，因此在模擬污染物擴(kuò)散時(shí)具有獨(dú)特的優(yōu)越性。本研究使用開(kāi)源代碼SWE-SPHysics(http：//www.sphysics.org)求解水流運(yùn)動(dòng)方程，并對(duì)其二次開(kāi)發(fā)，近似離散求解污染物的輸運(yùn)方程，設(shè)置一維/二維均勻流在不同擴(kuò)散系數(shù)下的污染物案例，并進(jìn)行模擬，在驗(yàn)證結(jié)果合理的前提下，對(duì)實(shí)際地形污染物擴(kuò)散的過(guò)程進(jìn)行預(yù)測(cè)模擬。

1 模擬方法

1.1 控制方程

忽略地球自轉(zhuǎn)引起的科氏力和流體粘度，拉格朗日描述下的SWE方程形式如下：

(1)

(2)

(3)

式中：ρ為密度，v為速度，g為重力加速度(g=9.81 m/s2)，b為河床粒子海拔，c為污染物濃度，D為污染物的濃度擴(kuò)散系數(shù),Sf=nv|v|/dw4/3，其中n為曼寧糙率系數(shù)。

1.2 水深求解

采用變密度法和變光滑長(zhǎng)度法求解水深，利用牛頓迭代法求解流體粒子密度和水深。流體粒子密度ρl：

ρl=ρw/d

(4)

(5)

(6)

式中：ρw表示流體體積密度(ρw=1 000 kg/m3)，m為流體粒子質(zhì)量，下標(biāo)“i”和“j”分別為i粒子和j粒子，W為核函數(shù)，d為水深，ρ0和h分別為初始密度和光滑長(zhǎng)度，“dm”為維度，下標(biāo)“0”表示初始狀態(tài)。

1.3 速度求解

淺水環(huán)境下，動(dòng)量方程的SPH離散格式如下：

(7)

式中：ti=Ti/mi，bi為i河床粒子梯度，ki=(b)。

(8)

(9)

修正矩陣計(jì)算公式：

(10)

1.4 濃度求解

淺水環(huán)境下，污染物輸運(yùn)方程的SPH離散格式如下：

(11)

式中：D為污染物的濃度擴(kuò)散系數(shù)，D也可是偏導(dǎo)數(shù)Dxx、Dxy、Dyx和Dyy，變擴(kuò)散系數(shù)公式：

DL=6.09×u*×h

(12)

DT=0.6×u*×h

(13)

(14)

(15)

(16)

2 模擬工況設(shè)計(jì)

2.1 一維均勻流模擬

設(shè)置長(zhǎng)10 000 m的矩形水槽，d=1 m，粒子間距△x=1 m，共10 000個(gè)流體粒子。模擬D為0、5、50 m2/s時(shí)污染物的擴(kuò)散過(guò)程。出入流邊界，vin/out=2 m/s，din/out=1 m，濃度為0 kg/m3，下標(biāo)“in/out”表示入流/出流，初始流體充滿計(jì)算域，v0=2 m/s，d0=1 m，c0設(shè)置如下：

(17)

一維分布源的解析解：

(18)

式中：c0為初始濃度，x1為濃度投放位置。

為驗(yàn)證SPH-SWE擴(kuò)散模型在不同擴(kuò)散系數(shù)的工況計(jì)算結(jié)果是否合理，設(shè)計(jì)模擬一維均勻流連續(xù)分布源。主要參數(shù)包括：其污染源類型都是分布源，在2 500～3 500 m處投放污染物濃度為1 kg/m3，擴(kuò)散系數(shù)分別取0、5、50 m2/s的3種工況。

2.2 經(jīng)典潰壩案例

設(shè)置經(jīng)典干、濕潰壩案例，不考慮擴(kuò)散影響。干、濕河床的計(jì)算區(qū)域長(zhǎng)均為2 000 m，△x=10 m，上游d0和c0分別為10 m和0.7 kg/m3，下游d0和c0分別為5 m和0.5 kg/m3，下游無(wú)水，v0=0 m/s，模擬瞬時(shí)潰壩后污染物濃度擴(kuò)散過(guò)程，與解析解[11]對(duì)比驗(yàn)證。

(19)

2.3 二維均勻流下的污染物擴(kuò)散模擬

設(shè)置二維均勻流的驗(yàn)證工況，計(jì)算域1 200 m×400 m，△x=1 m，初始狀態(tài)時(shí)，流體粒子充滿計(jì)算域，坡度為0.001，初始vi和d分別為2.9 m/s和5 m，vin、din和cin分別為2.9 m/s、5 m和0 kg/m3，分別模擬D為0，10 000 m2/s時(shí)污染物濃度的分布情況。

2.4 實(shí)際地形模擬

選取位于青海省互助縣南門(mén)峽河上南門(mén)峽水庫(kù)作為研究區(qū)域，庫(kù)區(qū)海拔高程在2 700 m以上，占地面積約為218 km2，地形圖如圖1所示。鄭仙佩[12]對(duì)南門(mén)峽水庫(kù)區(qū)域進(jìn)行了潰壩模擬研究，取△x為10、15、20 m進(jìn)行模擬，并對(duì)模擬結(jié)果進(jìn)行了收斂性分析，得出的3種結(jié)果基本一致，驗(yàn)證了SPH-SWE模型的收斂性和結(jié)果的可靠性。

圖1 南門(mén)峽水庫(kù)Fig.1 Nanmenxia reservoir

模擬污染物擴(kuò)散系數(shù)取0 m2/s和變擴(kuò)散系數(shù)時(shí)，南門(mén)峽水庫(kù)瞬間潰壩后污染物的擴(kuò)散過(guò)程。計(jì)算區(qū)域?yàn)? 380 m×10 800 m，比降0.012，糙率0.02。河床△x=30 m，共89 167個(gè)河床粒子，流體△x=20 m，共2 664個(gè)流體粒子，污染物坐標(biāo)位置(X：2 450 m，Y：9 255 m)，污染物投放半徑為100 m，濃度為1 kg/m3，其它區(qū)域co為0 kg/m3，模擬時(shí)長(zhǎng)600 s。

3 結(jié)果分析

3.1 一維均勻流模擬結(jié)果分析

圖2為一維瞬時(shí)點(diǎn)源的解析解和模擬值對(duì)比圖。圖2a為當(dāng)t=2 000 s時(shí)，計(jì)算域內(nèi)流體粒子的水深-位移圖和速度-位移圖，水深和速度分別為常數(shù)1 m和2 m/s。圖2b是Dx=0 m2/s，t為0、1 300、2 600 s時(shí)，解析解與模擬結(jié)果吻合。圖2c和圖2d分別是Dx取5 m2/s和50 m2/s，t為0、1 300、2 600 s，模擬值與解析解一致，上述3種工況均采用SPH-SWE擴(kuò)散模型模擬污染物的擴(kuò)散過(guò)程，其結(jié)果合理，驗(yàn)證了擴(kuò)散系數(shù)大小對(duì)模擬結(jié)果無(wú)影響。

圖2 一維瞬時(shí)點(diǎn)源的解析解和模擬值對(duì)比圖Fig.2 Comparison of analytical and simulation results with 1D instantaneous point source

3.2 經(jīng)典潰壩模擬結(jié)果分析

圖3為t=50 s時(shí)干河床潰壩后水深和濃度分布圖，圖3a為干河床潰壩計(jì)算域內(nèi)的水深變化趨勢(shì)圖，圖3b為污染物濃度為常數(shù)0.7 kg/m3。圖4為t=50 s時(shí)濕河床潰壩后水深和濃度分布圖，采用MATLAB求得解析解，將解析解與模擬值進(jìn)行對(duì)比。圖4a中，水深的數(shù)值結(jié)果接近解析解，在不連續(xù)處未出現(xiàn)數(shù)值震蕩，證明了水流計(jì)算結(jié)果的合理性。圖4b為污染物濃度的模擬值和解析解一致，驗(yàn)證了SPH-SWE模型具有較準(zhǔn)確的計(jì)算結(jié)果。

圖3 干河床潰壩后水深和濃度分布圖Fig.3 Distribution of water depth and concentration after the dam break of dry bed

圖4 濕河床潰壩后水深和濃度分布圖Fig.4 Distribution of water depth and concentration after the dam break of wet bed

3.3 二維均勻流模擬結(jié)果分析

圖5為二維均勻流中D取0 m2/s時(shí)，僅考慮移流作用下污染物隨時(shí)間變化的濃度分布圖。由圖5可知，當(dāng)t為0、50、120、200、300 s時(shí)，污染物擴(kuò)散到達(dá)最右端的距離是均勻水流速度與時(shí)間的乘積，分別對(duì)應(yīng)200、345、548、786、1 076 m，驗(yàn)證了移流擴(kuò)散過(guò)程中污染物濃度隨時(shí)間擴(kuò)散的模擬結(jié)果具有可靠性。

圖5 移流作用下二維均勻流中污染物濃度分布圖Fig.5 Distribution of pillutant concentration in two-dimensional uniform flow under the effects of advection

圖6為D取10 000 m2/s時(shí)考慮移流和擴(kuò)散作用下污染物隨時(shí)間變化的濃度分布圖。由圖6可知，當(dāng)t為0、50、120、200、300 s時(shí)的污染物濃度分布，與圖5污染物擴(kuò)散源的中心位置保持一致，由于擴(kuò)散效應(yīng)導(dǎo)致污染物的擴(kuò)散范圍比圖5更大。根據(jù)模擬結(jié)果的合理性，驗(yàn)證了SPH-SWE擴(kuò)散模型在模擬二維水流仍然具有很好的穩(wěn)定性，而且不需考慮移流擴(kuò)散項(xiàng)的求解，只需對(duì)其余部分離散求解，結(jié)果表明SPH-SWE方法在污染物擴(kuò)散模擬方面具有很好的優(yōu)勢(shì)。

圖6 移流擴(kuò)散作用下二維均勻流中污染物濃度分布圖Fig.6 Distribution of pillutant concentration in two-dimensional uniform flow under the effects of advection and diffusion

3.4 實(shí)際地形預(yù)測(cè)結(jié)果分析

圖7分別是當(dāng)t為600 s時(shí)，南門(mén)峽水庫(kù)潰壩后在D取0 m2/s和變擴(kuò)散系數(shù)條件下，污染物擴(kuò)散的濃度分布圖。鄭仙佩[12]對(duì)南門(mén)峽水庫(kù)潰壩的水流運(yùn)動(dòng)過(guò)程進(jìn)行了細(xì)致的討論，且驗(yàn)證粒子間距取20 m的方案是可行的，這里不做重復(fù)討論，著重考慮潰壩時(shí)污染物擴(kuò)散濃度分布情況。

圖7 南門(mén)峽水庫(kù)潰壩后濃度分布圖Fig.7 Distribution of dam break concentration in Nanmenxia reservoir

圖7a為當(dāng)D=0 m2/s，時(shí)間取600 s時(shí)污染物在實(shí)際潰壩后的濃度分布圖。由圖7a可知，污染物的遷移過(guò)程主要與水流流速有著很大的關(guān)系，污染物主要集中在水流的中間位置。但僅考慮污染物的移流過(guò)程，濃度會(huì)出現(xiàn)局部不連續(xù)，實(shí)際水流不會(huì)出現(xiàn)此情況。由于此模擬工況的前提假設(shè)是不考慮擴(kuò)散效應(yīng)，污染物在輸運(yùn)過(guò)程中只與流體粒子的運(yùn)動(dòng)軌跡有關(guān)。圖7b為在變擴(kuò)散系數(shù)下時(shí)間取600 s時(shí)，南門(mén)峽水庫(kù)潰壩后污染物濃度的分布圖，由圖7b可知，在沿河道水流間位置污染物濃度達(dá)到峰值0.075 6 kg/m3。其次，濃度值依次向上、下游遞減，不會(huì)出現(xiàn)圖7a中濃度不連續(xù)的情況，濃度未出現(xiàn)數(shù)值奇異點(diǎn)，符合實(shí)際河流中污染物擴(kuò)散規(guī)律。圖7b污染物的變擴(kuò)散系數(shù)與地形起伏、水流速度和水深有關(guān)，接近實(shí)際寬淺河流中污染物擴(kuò)散的過(guò)程，驗(yàn)證了SPH-SWE擴(kuò)散模型在復(fù)雜水流環(huán)境中具有很好的穩(wěn)定性，得出SPH-SWE方法在污染物輸運(yùn)過(guò)程的預(yù)測(cè)具有很好的研究?jī)r(jià)值，可較好預(yù)測(cè)可溶保守型污染物在河流中的濃度分布情況。

4 討論與結(jié)論

污染物在水體中的遷移擴(kuò)散規(guī)律是水環(huán)境分析的基礎(chǔ)內(nèi)容，對(duì)研究河流中污染物擴(kuò)散規(guī)律具有重要意義。劉圣勇[13]采用有限差分法對(duì)一維河流污染物擴(kuò)散進(jìn)行模擬。許媛媛[14]采用有限體積法對(duì)二維淺水方程和物質(zhì)輸運(yùn)方程進(jìn)行求解，該模型采用高精度的迎風(fēng)重構(gòu)技術(shù)處理網(wǎng)格界面的污染物輸運(yùn)，雖然計(jì)算精度有所提高，但是計(jì)算效率較低和求解格式較復(fù)雜。本文采用SPH-SWE法對(duì)河流中可溶保守型污染物擴(kuò)散進(jìn)行模擬，由于SPH方法具有拉格朗日的特性，在模擬污染物對(duì)流擴(kuò)散過(guò)程具有優(yōu)勢(shì)。通過(guò)設(shè)計(jì)一維和二維均勻流工況，模擬保守型溶質(zhì)污染在不同擴(kuò)散系數(shù)下的遷移擴(kuò)散過(guò)程，數(shù)值模擬結(jié)果和解析解一致；同時(shí)選取經(jīng)典潰壩案例進(jìn)行驗(yàn)證，模擬結(jié)果與解析解高度吻合，表明了SPH-SWE法對(duì)河流中污染物的遷移擴(kuò)散模擬研究具有很大潛力；通過(guò)預(yù)測(cè)南門(mén)峽水庫(kù)潰壩后污染物濃度分布情況，得到在沿河道水流中間位置污染物濃度達(dá)到峰值0.075 6 kg/m3，可預(yù)測(cè)污染物在實(shí)際地形中的濃度變化趨勢(shì)和規(guī)律。本研究主要有以下幾點(diǎn)結(jié)論：

(1)SPH-SWE擴(kuò)散模型模擬一維、二維均勻流和經(jīng)典潰壩工況，避免了污染物濃度奇異值的出現(xiàn)，驗(yàn)證了該模型具有很好的穩(wěn)定性。

(2)驗(yàn)證了出入流邊界條件下，SPH-SWE擴(kuò)散模型可以很好地模擬可溶保守型污染物的遷移擴(kuò)散過(guò)程。

(3)實(shí)際地形預(yù)測(cè)污染物模擬結(jié)果合理，表明了SPH-SWE法可用于實(shí)際地形中河流污染物擴(kuò)散規(guī)律的研究，可較好預(yù)測(cè)可溶保守型污染物在河流中的濃度分布情況。