基于MLPG無網(wǎng)格法的土石壩非穩(wěn)態(tài)飽和滲流場分析

鄭保敬 鄒申威 雷 未 葉 永

(三峽大學(xué) 水利與環(huán)境學(xué)院，湖北 宜昌 443002)

滲流自由面作為滲流邊界的一部分，它的位置事先是未知的，并且隨著上游水位的變化而變化，因此該類問題屬于動邊界問題.解決此類問題的數(shù)值方法一般采用迭代法.滲流自由面的確定往往依賴于網(wǎng)格劃分的有限元法，其中網(wǎng)格按動態(tài)性質(zhì)分為固定網(wǎng)格和移動網(wǎng)格兩大類[1].固定網(wǎng)格在系統(tǒng)迭代處理過程中，網(wǎng)格形狀始終保持不變，根據(jù)研究內(nèi)容的差異性，又衍生出不同固定網(wǎng)格法，如耦合單元法、虛位單元法、初位流量法和剩余流量法等.而移動網(wǎng)格法在滲流迭代計算過程中，通過不間斷動態(tài)調(diào)整網(wǎng)格位置來自適應(yīng)滲流自由面位置，直至滿足最小二乘數(shù)的控制誤差.該系統(tǒng)求解方法同樣依賴移動網(wǎng)格的頻數(shù)，每一次模型迭代都將重新劃分網(wǎng)格，造成計算過程復(fù)雜冗長的問題.在解決動邊界問題時，無網(wǎng)格法可擺脫網(wǎng)格的限制，只需通過一組離散節(jié)點來構(gòu)造形函數(shù)，具有處理方便、布置靈活的特點，可有效擺脫有限元對網(wǎng)格過分依賴的局限性[2-4].

無網(wǎng)格法的研究始于上世紀(jì)七十年代光滑粒子流體動力學(xué)(SPH)方法的提出.隨后許多學(xué)者對無網(wǎng)格法進行了研究，發(fā)展至今已經(jīng)形成了很多種方法，包括再生核粒子法(RKPM)，無單元Galerkin法(EFG)，有限點法(FPM)，復(fù)變量無網(wǎng)格法以及基于邊界積分方程的無網(wǎng)格法[5]等.1998年Atluri等[6]首次提出了無網(wǎng)格局部Petrov-Galerkin法(MLPG)，該方法在構(gòu)造模型形函數(shù)和數(shù)值積分上均不需要網(wǎng)格作為求解基礎(chǔ)，是真正意義上的無網(wǎng)格法.

近幾年，有部分學(xué)者將無網(wǎng)格法運用至飽和土體的穩(wěn)態(tài)滲流問題上.如李樹忱等[7]對模型重心采用基于拉格朗日原理的配點插值無網(wǎng)格法，來模擬穩(wěn)態(tài)下的滲流自由面問題;陳銳等[8]通過比較自適應(yīng)節(jié)點與靜態(tài)均勻節(jié)點在無網(wǎng)格法下的數(shù)值解和解析解誤差，從而提出一種節(jié)點可自適應(yīng)化無網(wǎng)格法;唐雄等[9]通過無網(wǎng)格法中的物質(zhì)點法(MPM)探究了大變形巖土工程問題和多相流問題，發(fā)現(xiàn)模型計算所得的數(shù)值解和理論解吻合度很好，且多相物質(zhì)配點法能夠較為準(zhǔn)確表述非飽和場的孔隙介質(zhì)滲流特征和飽和場的孔隙介質(zhì)水力特征;楊悅[10]提出將邊界積分方程與最小二乘法相結(jié)合來求解定解問題的微分方程，通過編制Matlab程序求解變系數(shù)滲流問題，該方法不僅節(jié)點布置靈活，還能有效避免因奇異矩陣對模型精度的影響.目前無網(wǎng)格法在土石壩非穩(wěn)態(tài)滲流問題的研究相對較少，本文采用基于滑動Kriging插值的無網(wǎng)格局部Petrov-Galerkin法[11]求解水位變化情況下土石壩的穩(wěn)態(tài)和非穩(wěn)態(tài)滲流問題.該方法利用滑動Kriging插值法構(gòu)造近似函數(shù)，使算法具備Kronecker函數(shù)特性，并直接賦予本質(zhì)邊界條件，將Heaviside分段函數(shù)看作權(quán)函數(shù)，使得轉(zhuǎn)化的積分方程只剩邊界積分部分，從而減少計算量，最后采用數(shù)值算例對該方法的有效性進行驗證.

1 二維土石壩非穩(wěn)定飽和滲流問題

考慮一均質(zhì)土石壩，如圖1所示在滲流飽和區(qū)域Ω內(nèi)任意一點A的水頭h為

圖1 土石壩非穩(wěn)定滲流模型

式中:y為A點垂直坐標(biāo)分量;p為孔隙水壓力;γw為水的重度.

符合達西定律的二維非穩(wěn)定飽和滲流問題的控制方程可以表示為

不規(guī)范處方主要包括：①處方超量，如普通處方超過7天量，急診處方超過3天量，省特約處方超過7天量，無特殊治療需要處方超過1個月用量等；②臨床診斷書寫不全、病人信息不全；③醫(yī)生輸入錯誤；④醫(yī)師處方未簽名（蓋章）或與留樣不一致；⑤處方上藥品品種數(shù)量超過規(guī)定（5種）；⑥處方修改未簽名蓋章，也未注明修改日期；⑦處方用箋不規(guī)范，如門診手寫處方用的是住院處方或自費藥房處方。

邊界條件:

初始條件:

式中:kx和ky分別為沿x和y方向的滲透系數(shù);Ss=ρg(α+nβ)為單位貯存率;ρ為水的密度;α為土料的體積壓縮系數(shù);n為孔隙率;β為水的壓縮系數(shù);μ為給水度，表示為飽和土體含水層自由面其單位截面積吸收或釋放最大水量與土體總體積比率;θ為自由面外法線方向與豎直方向的夾角.

2 二維非穩(wěn)定飽和滲流問題的MLPG無網(wǎng)格法

2.1 滑動Kriging插值形函數(shù)

Kriging插值方法采用無偏估計和線性回歸模型，類似于移動最小二乘法(MLS)，該加權(quán)為近似可連續(xù)的移動，所以稱為滑動Kriging插值法.定義在問題域Ω及其邊界Γ上的場函數(shù)為h(x)，將系統(tǒng)求解域離散為N個場節(jié)點分布，且空間內(nèi)任意點x的局部支持域又含n個場節(jié)點，若已知節(jié)點的場函數(shù)值為h(x1)，h(x2)，…，h(xn)，則點x處的場函數(shù)值可以通過已知節(jié)點的場函數(shù)值來逼近

式中:Φ(x)=[φ1(x)，φ2(x)，…，φn(x)]為通過滑動Kriging插值法構(gòu)造的形函數(shù)，表達式為

式(10)中S1和S2分別為

式(12)與(13)中PT和R分別為

式(14)～(16)中:PT(x)為多項式基函數(shù);R[R(xi，xj)]為對角線為1的對稱相關(guān)矩陣;R(xi，xj)為節(jié)點xi與xj之間的相關(guān)函數(shù)，通常取Gaussian型函數(shù).

式中:θ＞0，為模型參數(shù);rij為節(jié)點xi與xj之間的距離.

2.2 非穩(wěn)態(tài)飽和滲流問題的MLPG法

采用N個場節(jié)點將整場滲流域Ω與邊界域Γ進行離散化處理，局部支持域內(nèi)的任意點I選擇加權(quán)余量法近似求解，式(2)的積分弱形式為

式中:wI為與節(jié)點I對應(yīng)的權(quán)函數(shù);Ωq為與節(jié)點I對應(yīng)的局部子域.

當(dāng)考慮求解域Γ與局部積分域邊界Γq重疊時，Γq(Γq=Γqi+Γqt+Γqu)由以下3部分構(gòu)成:Γq與Γ不重疊，記為Γqi;Γq與自然邊界重疊區(qū)域，記為Γqt;Γq與本質(zhì)邊界重疊區(qū)域，記為Γqu.對式(18)引入自然邊界條件，并令Ss=0可得

在非穩(wěn)態(tài)滲流問題中，系統(tǒng)內(nèi)場節(jié)點的物理函數(shù)隨空間域和時間域的改變而變化，其關(guān)系可近似為

式中:φk為形函數(shù);hk為給定的節(jié)點處的函數(shù)值.

將式(20)代入式(19)中得到節(jié)點I的離散系統(tǒng)方程為

式(21)的矩陣形式為

式中:

將所有編號的各節(jié)點按順序依次代入離散控制方程進行集合，從而得到總體控制方程式為

對于時間域的離散，本文采用向后差分格式

代入總體系統(tǒng)方程為

3 數(shù)值算例

3.1 穩(wěn)定滲流

例1:考慮如圖2所示的矩形均質(zhì)土石壩，壩高為8 m，寬為4 m，上游水位6 m，下游水位1 m，土體滲透系數(shù)為k=0.016 m/d.假定第一類邊界條件:HΓ1=6.0 m，HΓ2=1.0 m，Γ5為潛在溢出面邊界.第二類邊界條件:Γ3為不透水邊界，Γ4無流量補給.

為了驗證初始自由面對迭代計算結(jié)果的影響，分別假設(shè)3條初始自由面如圖2所示，各初始自由面的求解域如圖3所示.

圖2 矩形土壩滲流模型(單位:m)

圖3 場節(jié)點布置示意圖(單位:m)

在每個區(qū)域內(nèi)均布置255(17×15)個節(jié)點，計算結(jié)果的絕對誤差取0.001 m.表1給出了不同數(shù)值方法所計算的自由面上的水頭值，從表1中可以看出本文方法計算的結(jié)果和實驗誤差最小，而且本文方法對初始自由面的位置不敏感，通過迭代都能很快收斂到最終自由面的值.

表1 自由面上典型參考點在不同計算方法下水頭值比較 (單位:m)

3.2 非穩(wěn)定滲流

例2:文獻[13]中介紹的均質(zhì)土壩實驗?zāi)Ｐ?，將砂壩盛放在一個能調(diào)節(jié)兩端水位的水槽中，對不同速度水位上升和水位下降工況下的實驗結(jié)果與數(shù)值計算結(jié)果進行比較:下游水位穩(wěn)定在10 cm，上游水位由10 cm 在不同時間段(30、60、120、240、600、4 800 s)均勻上升至30 cm變化情況，計算結(jié)果與實驗結(jié)果如圖4所示;上游水位穩(wěn)定在30 cm，下游水位由30 cm在不同時間段(30、60、120、240、600、2 400 s)均勻下降至10 cm變化情況，計算結(jié)果與實驗結(jié)果如圖5所示.數(shù)值計算條件與模型試驗相同，初始時刻求解區(qū)域內(nèi)均勻布置5 151(101×51)個節(jié)點，時間步長取0.5 s.

圖4 水位上升時數(shù)值與實驗結(jié)果比較示意圖

圖5 水位下降時數(shù)值與實驗結(jié)果比較示意圖

由圖4～5可知，在不同速度水位上升和水位下降工況下，本文方法計算的滲流自由面位置與實驗結(jié)果均較為吻合，表明對非穩(wěn)定滲流問題采用MLPG法是有效可行的，且滿足精度要求.

例3:圖6為質(zhì)梯形壩示意圖，壩高17 m，頂寬6 m，底寬91 m，假設(shè)上游初始水位H1=15 m，并以勻速下降至H′1=5 m，下游水位為H2=2 m.梯形壩的貯水率為Ss=0，飽和滲透系數(shù)為k=0.861 m/d，給水度為μ=0.006.梯形壩求解區(qū)域內(nèi)均勻布置2 852(92×31)個節(jié)點，時間步長取0.01 d.圖7為上游水位在不同高度時按穩(wěn)態(tài)滲流問題計算得到的自由面位置，圖8和圖9分別為上游水位下降速度為0.1 m/d和1 m/d時的計算結(jié)果.

圖6 均質(zhì)梯形壩示意圖

圖7 不同水位下壩體穩(wěn)定滲流自由面示意圖

圖8 v=0.1 m/d自由面變化示意圖

圖9 v=1 m/d自由面變化示意圖

從圖中可以看出，當(dāng)上游水位下降速度很緩慢時瞬態(tài)結(jié)果和穩(wěn)態(tài)結(jié)果非常接近;而當(dāng)上游水位下降速度較大時自由面出現(xiàn)滯后現(xiàn)象，且后期水位在上游坡面的出滲點高于初始水位，從而造成涌水倒流現(xiàn)象，這種現(xiàn)象也符合工程實際情況.

4 結(jié) 論

本文采用MLPG無網(wǎng)格法求解均質(zhì)土石壩的非穩(wěn)態(tài)飽和滲流問題，該方法使用滿足狄里克萊函數(shù)性質(zhì)的滑動Kriging插值法構(gòu)造形函數(shù)，對模型可直接施加本質(zhì)邊界條件;Heaviside分段函數(shù)作為權(quán)函數(shù)時，系統(tǒng)內(nèi)剛度矩陣只需計算邊界積分，不需計算區(qū)域積分，大幅提高了計算效率.數(shù)值算例結(jié)果表明本文方法計算的滲流自由面最終位置與初始自由面的設(shè)置無關(guān)，只在迭代次數(shù)上有差異，因此本文方法具有良好的穩(wěn)定性，適用于求解動邊界問題.