變異Henry問題檢驗(yàn)網(wǎng)格Peclet數(shù)在數(shù)值計(jì)算中的重要性

馬 倩，李海龍

(1.生物地質(zhì)與環(huán)境地質(zhì)國家重點(diǎn)實(shí)驗(yàn)室/中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083;2.地下水循環(huán)與環(huán)境演化教育部重點(diǎn)實(shí)驗(yàn)室/中國地質(zhì)大學(xué)(北京)，北京 100083)

有限元法作為求解地下水流和溶質(zhì)運(yùn)移對流-彌散方程的常用數(shù)值方法，在水流模擬中有很多優(yōu)點(diǎn)。它有固定的網(wǎng)格，通常滿足質(zhì)量守恒定律;可以精確、高效地處理以彌散為主的問題;便于編寫程序并得以實(shí)現(xiàn)等［1］。但是，對于許多野外常見的以對流為主的問題，求解易引起顯著的數(shù)值振蕩。為此，許多學(xué)者進(jìn)行了大量的計(jì)算實(shí)踐和理論分析指出［2］:對于以對流占優(yōu)的溶質(zhì)運(yùn)移問題，有限差分法和有限元法計(jì)算時(shí)皆不可避免地出現(xiàn)數(shù)值振蕩。因此，掌握消除地下水流和溶質(zhì)運(yùn)移問題中濃度振蕩的方法和理論依據(jù)是很有意義的。

Henry問題研究的是承壓含水層中咸淡水靜力平衡問題，在地下水流和溶質(zhì)運(yùn)移領(lǐng)域具有極高的代表性［3］。其常用的概念模型為:承壓含水層均質(zhì)各向同性，長2m高1m;上下頂板、底板為隔水邊界;內(nèi)陸邊界為恒定的淡水流量邊界;海岸邊界為定水頭定濃度的海水邊界，海水水位1m［4］。它是一個(gè)理想的概念模型，為提高其實(shí)際適用性及科研價(jià)值(如提高模型參數(shù)敏感性)，很多學(xué)者在Henry問題的基礎(chǔ)上，研究水文地質(zhì)條件改變后的海水入侵問題，即變異Henry問題。Frind［5］研究的海水入侵問題，水動(dòng)力彌散系數(shù)與原Henry問題不同;Huyakorn等［6］將右邊界分為定濃度和零濃度通量邊界條件兩部分;Simpon和Clement［3］發(fā)現(xiàn)減小Henry問題中左邊界的內(nèi)陸淡水補(bǔ)給可以增加密度效應(yīng);Abarca等［7］綜合考慮含水層滲透系數(shù)的各向異性及右邊界水流方向?qū)舛鹊挠绊懀l(fā)現(xiàn)變異后的Henry問題對密度變化更敏感。本研究對內(nèi)陸和海岸邊界條件改變后的Henry問題進(jìn)行模擬求解，具體為:左邊界為定水頭邊界、右邊界水位為2m，發(fā)現(xiàn)該變異Henry問題對網(wǎng)格剖分要求更高，對網(wǎng)格Peclet數(shù)更為敏感(詳見2.3.3)，為消除濃度振蕩提供了實(shí)例，同時(shí)，驗(yàn)證和說明了網(wǎng)格Peclet數(shù)在數(shù)值計(jì)算中的重要性。

MARUN是專門用于地下水流和溶質(zhì)運(yùn)移模擬的數(shù)值程序，采用Galerkin二維有限元法求解對流-彌散方程，其可靠性已經(jīng)得到各方面的驗(yàn)證［8～10］。Li等［11］利用該程序?qū)|(zhì)海灘在海潮作用下的地下水鹽動(dòng)態(tài)進(jìn)行了有效的定量模擬，Li等［12］、Guo 等［13］、Xia等［14］成功定量擬合了阿拉斯加不同海灣地下水位、鹽分和示蹤劑濃度的野外實(shí)測值，同時(shí)對海灘地下水動(dòng)態(tài)水化學(xué)特征進(jìn)行了定量刻畫。

本研究以有限元法數(shù)值計(jì)算為例，采用MARUN模擬程序?qū)ψ儺怘enry問題進(jìn)行模擬求解，得到了不同網(wǎng)格剖分及水動(dòng)力彌散系數(shù)特選節(jié)點(diǎn)處的濃度穿透曲線，發(fā)現(xiàn)了濃度振蕩現(xiàn)象的存在。同時(shí)，分析找到了濃度振蕩的原因及合理的消除方法。模擬結(jié)果及分析表明網(wǎng)格Peclet數(shù)能夠有效地判定給定的網(wǎng)格剖分是否會(huì)引起濃度振蕩，對有限元法數(shù)值計(jì)算的網(wǎng)格剖分具有指導(dǎo)意義。

1 模型介紹及網(wǎng)格剖分

1.1 模型介紹

變異Henry問題的概念模型具體為:承壓含水層均質(zhì)各向同性，長2m高1m;上下頂板、底板為隔水邊界;內(nèi)陸邊界為定水頭邊界，淡水水頭值為2.0632m;海岸邊界為定水頭定濃度的海水邊界，海水水位為2m。原Henry問題的左邊界流量為6.6×10-5m/s，本文提出的變異Henry問題所對應(yīng)的左邊界流量為1.2×10-4m/s，流速大約是原Henry問題的兩倍。模擬區(qū)域及相關(guān)模型參數(shù)取值見圖1。

圖1 變異Henry問題的模擬區(qū)域及參數(shù)取值Fig．1 Model domain and parameters used in the modified Henry problem

Henry問題使用的水動(dòng)力彌散系數(shù)為D1=1.8857×10-5m2/s，本研究考慮不同水動(dòng)力彌散系數(shù)對模擬結(jié)果的影響，故加入水動(dòng)力彌散系數(shù)D2=4.0×10-5m2/s。

1.2 網(wǎng)格剖分

本研究選用MARUN模擬程序采用Galerkin有限單元法計(jì)算模擬區(qū)域內(nèi)各節(jié)點(diǎn)的濃度。剖分區(qū)域?yàn)殚L2m高1m的二維規(guī)則矩形，故采用直角三角形單元剖分，長度方向上的單元數(shù)是高度方向上的2倍，保證單元在整個(gè)區(qū)域上均勻分布;且所有的三角形單元為“等腰直角三角形”，不存在鈍角三角形單元，是理想的三角形網(wǎng)格。研究共選取兩種網(wǎng)格剖分:［Ⅰ］20行、40列;［Ⅱ］40行、80列。網(wǎng)格剖分Ⅰ共861個(gè)節(jié)點(diǎn)，1 600個(gè)單元，網(wǎng)格間距0.05m;網(wǎng)格剖分Ⅱ共3 321個(gè)節(jié)點(diǎn)，6 400個(gè)單元，網(wǎng)格間距0.025m。

2 結(jié)果與討論

2.1 模擬求解

使用MARUN程序?qū)σ韵氯N情況的變異Henry問題進(jìn)行模擬求解，得到各情況模擬區(qū)域濃度的時(shí)空變化情況。各節(jié)點(diǎn)的初始濃度為0mg/L，初始水頭由左右邊界水頭線性插值得到。設(shè)置初始時(shí)間為0s，時(shí)間步長以1.1倍的因子增大或縮小，初始時(shí)間步長為10s，最大為 1 200s。

Case1:采用網(wǎng)格剖分Ⅰ，水動(dòng)力彌散系數(shù)為D1=1.8857×10-5m2/s;

Case2:采用網(wǎng)格剖分Ⅱ，水動(dòng)力彌散系數(shù)為D1=1.8857×10-5m2/s;

Case3:采用網(wǎng)格剖分Ⅰ，水動(dòng)力彌散系數(shù)為D2=4.0 ×10-5m2/s。

Case1與Case2僅網(wǎng)格剖分不同，模型參數(shù)、初始條件及邊界條件均相同。Case1與Case3僅水動(dòng)力彌散系數(shù)不同，其他模型參數(shù)、網(wǎng)格剖分、初始條件及邊界條件均相同。

2.2 結(jié)果

淡水從內(nèi)陸邊界流入，海水從海岸邊界驅(qū)入含水層，與排泄的淡水混合。內(nèi)陸和海岸邊界條件始終不變，所以入侵海水與淡水流場最終可以達(dá)到動(dòng)態(tài)平衡狀態(tài)，即得到承壓含水層中鹽水楔形體與淡水流場間的穩(wěn)態(tài)解。本研究模擬區(qū)域內(nèi)節(jié)點(diǎn)初始濃度均為0mg/L，故各節(jié)點(diǎn)的濃度應(yīng)先由0mg/L逐漸增加，達(dá)到穩(wěn)定狀態(tài)后不再隨時(shí)間變化［7］。

有限元法在處理以對流為主的地下水流和溶質(zhì)運(yùn)移問題時(shí)，網(wǎng)格剖分往往引起數(shù)值計(jì)算的數(shù)值振蕩，當(dāng)出現(xiàn)很陡的濃度鋒面時(shí)，也就是以對流為主時(shí)，振蕩尤為明顯［1］。為更清楚地研究和說明濃度振蕩的普遍存在，選取流速較小的節(jié)點(diǎn)A(x=1.5m，z=0.05m)處的濃度穿透曲線進(jìn)行研究。

2.2.1 不同網(wǎng)格剖分

如圖2所示，Case1、Case2為不同網(wǎng)格剖分節(jié)點(diǎn)A濃度穿透曲線的對比，兩種情況的節(jié)點(diǎn)A濃度由初始時(shí)刻的0mg/L近直線上升，至2h左右不再呈現(xiàn)明顯上升趨勢(圖2a);Case1節(jié)點(diǎn)濃度呈現(xiàn)隨時(shí)間在某值附近振蕩，并未穩(wěn)定于某一濃度值;Case2節(jié)點(diǎn)濃度最終穩(wěn)定于0.4mg/L(圖2b)。

如前所述，入侵海水與淡水流場達(dá)到動(dòng)態(tài)平衡后，濃度不再隨時(shí)間變化。而Case1節(jié)點(diǎn)濃度并未穩(wěn)定，說明使用網(wǎng)格剖分Ⅰ有限元法數(shù)值計(jì)算出現(xiàn)了數(shù)值振蕩，表現(xiàn)為節(jié)點(diǎn)處的濃度振蕩;Case2節(jié)點(diǎn)濃度達(dá)到穩(wěn)定，說明使用網(wǎng)格剖分Ⅱ有限元法數(shù)值計(jì)算未出現(xiàn)數(shù)值振蕩，表現(xiàn)為節(jié)點(diǎn)處濃度最終穩(wěn)定于某一濃度值，而不隨時(shí)間振蕩。Case1節(jié)點(diǎn)A濃度數(shù)值解明顯偏離Case2濃度數(shù)值解，表明使用網(wǎng)格剖分Ⅰ有限元數(shù)值計(jì)算存在數(shù)值誤差。即使用網(wǎng)格剖分Ⅰ有限元數(shù)值計(jì)算時(shí)既存在數(shù)值誤差，又出現(xiàn)了濃度振蕩。

圖2 三種情況的節(jié)點(diǎn)A處濃度穿透曲線Fig．2 Breakthrough curve of point A at(a)initial time，and(b)whole time of three cases

2.2.2 不同水動(dòng)力彌散系數(shù)

Case1、Case3為不同水動(dòng)力彌散系數(shù)節(jié)點(diǎn)A濃度穿透曲線的對比，兩種情況的節(jié)點(diǎn)A濃度由初始時(shí)刻的0mg/L近直線上升，至2h左右不再呈現(xiàn)明顯上升趨勢(圖2a);Case1濃度呈現(xiàn)隨時(shí)間在某值附近振蕩，并未穩(wěn)定于某一濃度值;Case3濃度最終穩(wěn)定于0.38mg/L(圖2b)。說明水動(dòng)力彌散系數(shù)為 D1=1.8857×10-5m2/s時(shí)，有限元法數(shù)值計(jì)算出現(xiàn)了數(shù)值振蕩，表現(xiàn)為節(jié)點(diǎn)處的濃度振蕩;水動(dòng)力彌散系數(shù)增大為D2=4.0×10-5m2/s時(shí)，有限元法數(shù)值計(jì)算得到了鹽水楔形體與淡水流場的穩(wěn)態(tài)解，數(shù)值解未出現(xiàn)數(shù)值振蕩，表現(xiàn)為節(jié)點(diǎn)處濃度最終穩(wěn)定于某一濃度值。另外，水動(dòng)力彌散系數(shù)為D2時(shí)，若對網(wǎng)格加密后(用網(wǎng)格Ⅱ)進(jìn)行模擬計(jì)算，則粗細(xì)網(wǎng)格剖分兩種情況下的節(jié)點(diǎn)處濃度穿透曲線幾乎完全重合，網(wǎng)格加密不再影響計(jì)算的精度。說明此時(shí)的網(wǎng)格剖分Ⅰ已滿足模擬計(jì)算精度要求，得到了精確的穩(wěn)態(tài)數(shù)值解。

2.3 討論

2.3.1 濃度振蕩消除方法

Case1、Case2僅網(wǎng)格剖分不同，數(shù)值計(jì)算得出的特選節(jié)點(diǎn)處濃度穿透曲線特征明顯不同。Case1節(jié)點(diǎn)A濃度振蕩，并未穩(wěn)定于某一濃度值;Case2節(jié)點(diǎn)A濃度最終達(dá)到穩(wěn)定。即使用網(wǎng)格剖分Ⅰ時(shí)數(shù)值計(jì)算出現(xiàn)濃度振蕩，網(wǎng)格剖分Ⅱ在網(wǎng)格剖分Ⅰ的基礎(chǔ)上將網(wǎng)格單元整體加密一倍，消除了濃度振蕩。

Case1、Case3僅水動(dòng)力彌散系數(shù)不同，特選節(jié)點(diǎn)處濃度穿透曲線特征明顯不同。Case1節(jié)點(diǎn)濃度振蕩，Case3節(jié)點(diǎn)濃度最終達(dá)到穩(wěn)定。即水動(dòng)力彌散系數(shù)為D1時(shí)數(shù)值計(jì)算出現(xiàn)濃度振蕩，而水動(dòng)力彌散系數(shù)D2在D1的基礎(chǔ)上增加約一倍，消除了濃度振蕩。

故有限元法求解地下水流和溶質(zhì)運(yùn)移問題時(shí)，若出現(xiàn)濃度數(shù)值解在某值附近振蕩的現(xiàn)象，可以通過加密網(wǎng)格將其消除，同時(shí)須考慮水動(dòng)力彌散系數(shù)設(shè)置是否合理，通過增加水動(dòng)力彌散系數(shù)也可消除濃度振蕩。另外，Case1與Case2采用不同網(wǎng)格處理同一問題，說明同一問題對網(wǎng)格精度的要求不同;Case1與Case3均采用網(wǎng)格剖分Ⅰ，僅模型參數(shù)不同，說明同一網(wǎng)格對不同模型參數(shù)的有效性不同。已有研究結(jié)果表明［3～7］，Henry問題在采用網(wǎng)格剖分Ⅰ時(shí)能得到精確的數(shù)值解，并不出現(xiàn)濃度振蕩，說明即使是研究區(qū)域相同，不同的邊界條件對網(wǎng)格精度的要求不同。網(wǎng)格剖分要具體問題具體分析，不能單純認(rèn)為邊界條件相同，網(wǎng)格就適用于不同模型參數(shù)取值。

2.3.2 理論依據(jù)

網(wǎng)格Peclet數(shù)(Pe)是衡量對流項(xiàng)在溶質(zhì)運(yùn)移問題中主導(dǎo)程度的無量綱參數(shù):

式中:v——達(dá)西速度;

D——水動(dòng)力彌散系數(shù);

Δx——網(wǎng)格間距。

對于純對流問題，D→0，Pe→∞。隨著物理彌散變得明顯，Peclet數(shù)則變小。同時(shí)，Peclet數(shù)還依賴于網(wǎng)格間距Δx;Δx減小時(shí)，Peclet數(shù)也減小。即網(wǎng)格剖分越細(xì)，Peclet數(shù)也越小，這意味著采用較小的網(wǎng)格間距可以減小數(shù)值振蕩。若Δx對應(yīng)的Pe＜2，振蕩就消除了［1］。

三種情況下對應(yīng)的最大網(wǎng)格Peclet數(shù)如下:

Case1:采用網(wǎng)格剖分Ⅰ，Δx=0.05m;D=1.8857×10-5m2/s。速度v最大值為0.0015m/s，故模擬區(qū)域內(nèi)最大網(wǎng)格Peclet數(shù)為Pe1=3.97。

Case2:采用網(wǎng)格剖分Ⅱ，Δx=0.025m;D=1.8857×10-5m2/s。速度v最大值為0.0015m/s，故模擬區(qū)域內(nèi)最大網(wǎng)格Peclet數(shù)為Pe2=1.98。

Case3:采用網(wǎng)格剖分Ⅰ，Δx=0.05m;D=4.0×10-5m2/s。速度v最大值為0.0015m/s，故模擬區(qū)域內(nèi)最大網(wǎng)格Peclet數(shù)為Pe3=1.88。

變異Henry問題的內(nèi)陸和海岸邊界條件不變，模擬區(qū)域內(nèi)各點(diǎn)的濃度最終應(yīng)達(dá)到動(dòng)態(tài)平衡狀態(tài)。而Case1濃度數(shù)值解出現(xiàn)振蕩現(xiàn)象，與真實(shí)的物理意義不符，此時(shí)Pe1＞2;Case2濃度數(shù)值解最終穩(wěn)定于某一濃度值，達(dá)到了動(dòng)態(tài)平衡狀態(tài)，符合Henry問題真實(shí)的物理意義，此時(shí)Pe2＜2;Case3濃度數(shù)值解最終穩(wěn)定于某一濃度值，符合真實(shí)的物理意義，此時(shí)Pe3＜2。即不滿足Pe＜2時(shí)，使用給定網(wǎng)格剖分進(jìn)行數(shù)值計(jì)算時(shí)，數(shù)值解會(huì)出現(xiàn)濃度振蕩(Case1);滿足Pe＜2時(shí)，使用給定的網(wǎng)格剖分進(jìn)行數(shù)值計(jì)算時(shí)，數(shù)值解不出現(xiàn)濃度振蕩(Case2、Case3)，說明Peclet數(shù)能夠有效地判定給定的網(wǎng)格剖分是否會(huì)引起濃度振蕩。使用有限元法進(jìn)行數(shù)值計(jì)算時(shí)，必須考慮Peclet數(shù)，使其盡量小于2。

2.3.3 變異Henry問題與Henry問題的對比

Henry問題與變異Henry問題的區(qū)別是，Henry問題的內(nèi)陸邊界為定流量邊界，流量為6.6×10-5m/s，變異Henry問題為定水頭邊界，水頭值為2.0632m;Henry問題的海水水位為1m，變異Henry問題的海水水位為2m;變異Henry問題左邊界平均流速為1.2×10-4m/s，大約是原Henry問題左邊界流速的兩倍，故變異Henry問題的對流作用更強(qiáng)，為滿足Peclet數(shù)小于2，要求網(wǎng)格剖分更密，即對網(wǎng)格剖分精度要求更高，更適于檢驗(yàn)網(wǎng)格剖分是否會(huì)引起濃度振蕩。海水入侵范圍較Henry問題小，50%等氯線與底板交于1.55m處(Henry問題50%等氯線與底板交于1.40m處［3～7］)。變異Henry問題的求解為消除地下水流和溶質(zhì)運(yùn)移問題數(shù)值求解中的數(shù)值振蕩提供了實(shí)例。

地下水?dāng)?shù)值模擬軟件常用的求解方法是有限元法和有限差分法，二者在數(shù)值計(jì)算時(shí)均可能出現(xiàn)數(shù)值振蕩。若出現(xiàn)數(shù)值解在某值附近振蕩的現(xiàn)象，可以通過加密網(wǎng)格將其消除，同時(shí)需要考慮模型參數(shù)設(shè)置是否合理，增加水動(dòng)力彌散系數(shù)也可消除。進(jìn)行網(wǎng)格剖分時(shí)，即使研究區(qū)域相同，不同的邊界條件、不同的水動(dòng)力彌散系數(shù)對網(wǎng)格精度的要求不同，同一網(wǎng)格對不同模型參數(shù)的有效性也不同，要具體問題具體分析，不可一概而論。

3 結(jié)論

(1)濃度振蕩可以通過加密網(wǎng)格、增加水動(dòng)力彌散系數(shù)消除。

(2)即使研究區(qū)域相同，不同的邊界條件、不同的水動(dòng)力彌散系數(shù)對網(wǎng)格精度的要求不同，同一網(wǎng)格對不同模型參數(shù)的有效性也不同。

(3)網(wǎng)格Peclet數(shù)能夠有效地判定給定的網(wǎng)格剖分是否會(huì)引起濃度振蕩，對網(wǎng)格剖分具有指導(dǎo)意義。

(4)在地下水流和溶質(zhì)運(yùn)移有限元數(shù)值計(jì)算中，要綜合考慮模型參數(shù)的合理性，分析地下水流場，對流作用較強(qiáng)的區(qū)域要加密剖分，滿足Peclet數(shù)小于2，這樣既可以消除濃度振蕩，又保證計(jì)算速度。

［1］ 鄭春苗，Gordon D Bennett.地下水污染物遷移模擬［M］.2版.北京:高等教育出版社，2009:109-118.［ZHENG C M，Gordon D Bennett.Applied Contaminant TransportModeling［M］.2nd ed.Beijing:Higher Education Press，2009:109-118.(in Chinese)］

［2］ 馮紹元.解地下水中溶質(zhì)運(yùn)移問題的特征有限元法［J］.武漢水利電力學(xué)院學(xué)報(bào)，1991，24(4):410-417.［Feng S Y.Method of Characteristic Finite Element for Solving Solute Transport in an Aquifer［J］.Journal of Wuhan institute of water conservancy and electric power，1991，24(4):410-417.(in Chinese)］

［3］ Simpson M J，T P Clement.Improving the worthiness of the Henry problem as a benchmark for densitydependent groundwaterflow models［J］. Water Resources Research，2004，40:W01504.

［4］ Henry H R.Effects of dispersion on salt encroachment in coastal aquifers［J］.US Geological Survey，1964，Water-Supply Paper 1613-C.

［5］ Frind E O.Simulation of long-term transient densitydependent transport in groundwater［J］.Advances in Water Resources，1982，5:73-88.

［6］ Huyakorn P，Andersen P F，Mercer J，et al.Saltwater intrusion in aquifers:development and testing of a three-dimensional finite element model［J］.Water Resources Research，1987，23(2):293-312.

［7］ Abarca E，Carrera J，Vila X S，et al.Anisotropic dispersvie Henry problem［J］.Advances in Water Resources，2006，30:913-926.

［8］ Boufadel M C，M T Suidan，A D Venosa.A numerical model for density-and-viscosity-dependent flows in two-dimensionalvariably-saturated porous media［J］.Journal of Contaminant Hydrology，1999a，36:1-20.

［9］ Boufadel M C，M T Suidan，A D Venosa.Numerical modeling of water flow below dry salt lakes:effect of capillarity and viscosity［J］.Journal of Hydrology，1999b，221:55-74.

［10］ Boufadel M C.A mechanistic study of nonlinear solute transport in a groundwater-surface water system under steady state and transient hydraulic conditions［J］.Water Resources Research，2000，36:2549-2565.

［11］ Li Hailong，Boufadel M C，Weaver J W.Tideinduced seawater-groundwater circulation in shallow beach aquifers［J］.Journal of Hydrology，2008，352:211-224.

［12］ Li Hailong，Boufadel M C.Long-term persistence of oil from the Exxon Valdez spill in two-layer beaches［J］.Nature Geoscience，2010，3(2):96-99.

［13］ Guo Qiana， Li Hailong， Boufadel M C， et al.Hydrodynamics in a gravel beach and its impact on the Exxon Valdez oil［J］. JournalofGeophysical Research，2010，115:C12077.

［14］ Xia Yiqiang，Li Hailong，Boufadel M C，et al.Hydrodynamic factors affecting the persistence of the Exxon Valdez oil in a shallow bedrock beach［J］.Water Resources Research，2010，46:W10528.