鄭保敬 鄒申威 雷 未 葉 永

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

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

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

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

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

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

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

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

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

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

邊界條件:

初始條件:

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

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

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

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

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

式(10)中S1和S2分別為

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

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

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

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

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

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

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

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

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

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

式(21)的矩陣形式為

式中:

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

對(duì)于時(shí)間域的離散，本文采用向后差分格式

代入總體系統(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無(wú)流量補(bǔ)給.

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

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

圖3 場(chǎng)節(jié)點(diǎn)布置示意圖(單位:m)

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

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

3.2 非穩(wěn)定滲流

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

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

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

由圖4～5可知，在不同速度水位上升和水位下降工況下，本文方法計(jì)算的滲流自由面位置與實(shí)驗(yàn)結(jié)果均較為吻合，表明對(duì)非穩(wěn)定滲流問(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)個(gè)節(jié)點(diǎn)，時(shí)間步長(zhǎng)取0.01 d.圖7為上游水位在不同高度時(shí)按穩(wěn)態(tài)滲流問(wèn)題計(jì)算得到的自由面位置，圖8和圖9分別為上游水位下降速度為0.1 m/d和1 m/d時(shí)的計(jì)算結(jié)果.

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

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

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

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

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

4 結(jié) 論

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