楊曉宇，董建華*，吳曉磊

(1.蘭州理工大學(xué) 土木工程學(xué)院，甘肅 蘭州 730050；2.蘭州理工大學(xué) 西部土木工程防災(zāi)減災(zāi)教育部工程研究中心，甘肅 蘭州 730050)

電滲透系數(shù)與土顆粒大小無關(guān)，因此利用電滲法對軟土進(jìn)行處理時能夠得到不錯的效果。土體的電滲現(xiàn)象最早由俄國科學(xué)家Reuss發(fā)現(xiàn)，但在相當(dāng)長的時間內(nèi)一直未得到工程應(yīng)用，直到20世紀(jì)30年代，Casagrande[1]首次將電滲排水法應(yīng)用于鐵路路基的加固并取得了良好的效果，這一新的方法立刻引起了工程界廣泛的關(guān)注；Bjerrum等[2]給出了電滲排水加固法中電學(xué)參數(shù)的設(shè)計方法，極大地指導(dǎo)了后續(xù)電滲法的應(yīng)用研究。

基礎(chǔ)性的試驗研究[3-4]是促進(jìn)人們對電滲法認(rèn)知的重要手段。電勢差是電滲流發(fā)生的主要驅(qū)動力，不同的電壓會導(dǎo)致不同強(qiáng)度的電勢分布[5-9]，使得土體中孔隙水的定向流動產(chǎn)生差異，產(chǎn)生不同的電滲加固效果。平行排布、梅花、矩形是最為常見的排布形式。王柳江等[10]試驗研究指出，電滲對象中陽極的數(shù)量越多，加固后土體強(qiáng)度越高。李一雯等[11]的研究結(jié)果表明，平行排布產(chǎn)生的裂縫較為規(guī)律，平行錯位形式下排水效果最好，而梅花形排布情況下土體裂縫開展情況較為嚴(yán)重。陶燕麗[12]、張雷[13]等試驗研究指出陽極對于材料的敏感程度要高于陰極。鄭凌逶等[14]對電滲排水加固法的研究與應(yīng)用進(jìn)行了全面的總結(jié)，為后續(xù)的研究提供了方向性指導(dǎo)。眾多電滲試驗研究的成果，為后續(xù)電滲固結(jié)理論研究提供了必要支持。

Esrig[15]最先提出了1維模型的電滲固結(jié)理論，水力滲流與電滲流滿足線性疊加原理[16]；后續(xù)學(xué)者在此基礎(chǔ)上進(jìn)一步拓展了電滲固結(jié)理論的研究內(nèi)容[17-19]。李瑛等[20]建立了軸對稱形式下的電滲固結(jié)理論；王軍[21]、王柳江[22]等在1維電滲固結(jié)模型中引入了堆載作用的影響；吳輝[23]、胡黎明[24]等在電滲固結(jié)計算中考慮了電導(dǎo)率及土體參數(shù)的非線性，并進(jìn)行了數(shù)值求解。也有學(xué)者開展了多物理場耦合作用下的電滲固結(jié)研究并進(jìn)行了數(shù)值計算[25-27]。蘇金強(qiáng)等[28]針對平行排布形式情況，基于“分塊思路”建立了2維電滲固結(jié)理論。電滲固結(jié)理論研究為電滲法工程應(yīng)用提供了必要的理論與技術(shù)支持。

然而，電滲試驗研究[5-12]表明電滲過程存在有效電勢衰減現(xiàn)象，一方面，有效電勢的衰減使得電滲后期效率低下，電滲法能耗較高在一定程度上限制了電滲法的應(yīng)用；另一方面，電極間距相當(dāng)時，電勢場分布的2維效應(yīng)顯著，采用1維固結(jié)理論計算會高估電滲加固效果，而已有的電滲固結(jié)理論沒有考慮有效電勢衰減及電勢場的2維特性帶來的影響。針對上述問題，本文在前人研究的基礎(chǔ)上，依據(jù)平行、矩形排布形式的組成規(guī)律，提出了點狀電極單元的概念，并將其作為2維電滲固結(jié)理論的研究對象；針對研究對象中多電極的有效電勢衰減情況，依據(jù)黎曼求和方式的思路，將有效電勢衰減轉(zhuǎn)化為離散、多態(tài)的階梯曲線，并結(jié)合電滲與有效電勢衰減在時間上的統(tǒng)一性，采取了初始條件“繼承”的策略，將2維電滲固結(jié)全時間域上的求解轉(zhuǎn)化為對單一黎曼狀態(tài)段中間變量方程的反復(fù)迭代；基于Galerkin法給出了方程的有限元形式，并采用Рython語言進(jìn)行了數(shù)值實現(xiàn)；為了驗證該算法及確定有效電勢變化規(guī)律，進(jìn)行了電滲試驗研究，利用實測電勢變化進(jìn)行了多態(tài)下的計算，該方法能夠體現(xiàn)出有效電勢變化對電滲孔壓的影響，且方法簡潔、思路巧妙。

1 點狀電極單元

電滲排水固結(jié)法中，平行排布是最為常用的一種電極排布形式，當(dāng)異性電極間距遠(yuǎn)大于同性電極間距時，通常會將其簡化為1維模型，Esrig理論便是在這一基礎(chǔ)上建立的。電滲機(jī)理指出土中孔隙水會沿著電勢梯度下降的方向流動，即從陽極流向陰極，使得陽極區(qū)的土體加固效果較好。王柳江等[10]對這一問題進(jìn)行了試驗研究，結(jié)果表明陽極數(shù)量越多則土體加固效果越好?；谶@種考慮，將平行排布形式中相應(yīng)的陰極替換為陽極，可以得到矩形排布形式，即1根陰極周圍環(huán)繞8根陽極，如圖1所示。

圖1 矩形電極排布形式Fig.1 Rectangular electrode arrangement

矩形排布的分布特性決定了該形式無法被簡化為1維模型。無論是對平行形式還是對矩形形式的研究，依據(jù)排布形式的組成情況與重復(fù)規(guī)律，都可以歸結(jié)為對最小電極單元的研究。每一個最小電極單元由4根電極組成，矩形排布形式下最小電極單元是非對稱的，含有3根陽極與1根陰極，如圖2所示，其中，b為單元寬度，h為單元高度。

圖2 非對稱點狀電極單元Fig.2 Asymmetric punctiform electrode unit

最小電極單元是一個完整的2維電滲單元，能夠體現(xiàn)出電勢場與滲流場的2維空間描述。最小電極組成單元中，電極自身的大小遠(yuǎn)遠(yuǎn)小于單元尺寸，所以可以將電極當(dāng)作點狀電極，這樣單元的電勢分布可以按照靜電場的邊值問題進(jìn)行求解，電勢場分布計算與電偶極子類似。

點狀電極單元依據(jù)電極的極性可以分為對稱和非對稱兩種形式，單元中邊界分為同性邊界與異性邊界。電極極性的非對稱特性使得單元電勢場分布與均勻電場存在較大差異，此時不應(yīng)簡化為1維模型。

2 有效電勢衰減及多態(tài)處理

電滲過程中化學(xué)變化產(chǎn)生的氣體及陽極區(qū)土體孔隙水的定向流動，使得陽極與土體之間接觸電阻增大，造成施加在土體上的有效電勢發(fā)生衰減，這一現(xiàn)象已在諸多電滲試驗研究[5-12]中得到了證實。

2.1 有效電勢衰減下電滲固結(jié)計算的關(guān)鍵問題

事實上，有效電勢衰減的現(xiàn)象早已被試驗驗證，但是電滲固結(jié)理論中一直缺乏對該現(xiàn)象的考量，根本原因是在電滲固結(jié)控制方程中引入中間變量時，如果電勢分布與時間有關(guān)，那么中間變量中電勢項對于時間的導(dǎo)數(shù)不為0，使得變換后的方程仍然具有兩個場函數(shù)，無法進(jìn)行求解。Esrig理論中電滲固結(jié)的控制方程如下[15]：

式中：u(x,t) 為 孔隙水壓力，Рa；V(x)為電勢分布函數(shù)，V；kex為x方向的電滲透系數(shù)， m2/(V·s)；khx為x方向的水力滲透系數(shù)， m/s ； γw為水的重度， N/m3；mv為土體體積壓縮系數(shù)， Pa-1。

式（1）中，電勢分布函數(shù)僅是關(guān)于空間坐標(biāo)的函數(shù)，即表明電勢分布與時間無關(guān)，為了進(jìn)行方程的求解，引入中間變量 ζ(x,t)，表示如下：

引入中間變量后，控制方程（1）等價于式（3）:

式中，Cv為土體固結(jié)系數(shù)，Cv=khx/mvγw， m2/s。

變換后得到的中間變量方程式（3）是最簡形式的拋物線型偏微分方程，此時該方程是可以進(jìn)行解析求解的，在得到中間變量場的分布函數(shù) ζ(x,t)之后，利用式（2）可得到電滲的孔壓分布。

上述是Esrig理論對電滲控制方程的求解思路與方法，可以看到控制方程能夠進(jìn)行求解的關(guān)鍵在于中間變量 ζ(x,t) 中 電勢分布項與時間無關(guān)，這樣ζ(x,t)對時間t求導(dǎo)時，電勢項可以當(dāng)作常量，其對時間的導(dǎo)數(shù)為0，從而可以得到最簡形式的拋物線型偏微分方程。然而，有效電勢的衰減使得電勢分布函數(shù)與時間相關(guān)，這一客觀特性與中間變量函數(shù)引入時的要求相矛盾，這一矛盾關(guān)系制約著有效電勢衰減下電滲固結(jié)方程的求解，是求解計算的難點也是關(guān)鍵點。此外，如果直接在中間變量表達(dá)中引入電勢的時間效應(yīng)，會導(dǎo)致數(shù)學(xué)邏輯的嚴(yán)密性不足。上述分析表明，有效電勢衰減影響下的電滲固結(jié)計算的關(guān)鍵性問題在于如何處理有效電勢衰減的數(shù)學(xué)描述與中間變量函數(shù)引入時的限制條件之間的矛盾關(guān)系。

2.2 有效電勢衰減的離散多態(tài)方法

在2維模型中，有效電勢是電勢場邊值問題的邊界條件，有效電勢Ve(t)隨時間變化，引起電勢場分布的變化，已有的研究成果中對于該問題的解答是完全迎合電勢衰減的描述，而弱化中間變量引入過程中的數(shù)學(xué)要求。實際上，完全可以從一個相反的角度來考慮這個問題，即通過某種有效電勢衰減的描述，未完全滿足中間變量引入時的要求。

在微積分理論中，黎曼積分是一個十分重要、普遍的概念，它給出了定積分的精準(zhǔn)定義，而左側(cè)黎曼求和形式是定積分計算的基本方法和思路。對于任何一個有效電勢變化函數(shù)Ve(t)，該函數(shù)是時間的函數(shù)；當(dāng)采用左側(cè)黎曼求和形式時，如圖3所示，在積分區(qū)段內(nèi)進(jìn)行子狀態(tài)段的劃分，每一個狀態(tài)段內(nèi)函數(shù)值均采用區(qū)間段左端點處的函數(shù)值，即圖3中的Ve1、Ve2、Ve3，經(jīng)過黎曼求和形式的劃分，能夠保證在子區(qū)間內(nèi)函數(shù)值保持為常數(shù)，不會隨自變量發(fā)生變化，使得在引入中間變量時對時間的導(dǎo)數(shù)為0，從而可以獲得標(biāo)準(zhǔn)的拋物型偏微分方程。

圖3 左側(cè)黎曼求和形式Fig.3 Left Riemann summation form

黎曼求和的處理方式使得任意一個狀態(tài)段內(nèi)的函數(shù)值是一個常數(shù)，同時在整個積分域上又能逼近原函數(shù)的變化。這一特性十分適合用于描述有效電勢的衰減，因為在子狀態(tài)段內(nèi)有效電勢是一個常數(shù)，那么中間變量引入時的條件自然可以得到滿足，僅從這一點來看，在單一的子狀態(tài)段內(nèi)，該處理方式能夠解決這一核心矛盾關(guān)系。

2.3 電勢場分布

電滲法中施加的是直流電，因此電勢場滿足靜電場理論，將x-y坐標(biāo)系原點取在單元左下角的陰極處，單元模型如圖2所示。對有效電勢進(jìn)行離散、多態(tài)處理后，記第k黎曼狀態(tài)段內(nèi)場狀態(tài)變量電勢分布函數(shù)為Vk(x,y)，其滿足Laplace方程，如式（4）所示：

有效電勢衰減函數(shù)的離散、多態(tài)處理，是在時間域上的分割，而有效電勢是電勢場的邊界條件，電勢分布是一個邊值問題，不同的黎曼段內(nèi)的有效電勢不同，導(dǎo)致各個狀態(tài)段的電勢場分布存在一定的差異。邊界條件滿足Dirichlet條件，見式（5）：

對式（4）采用Galerkin法進(jìn)行離散，電勢場的空間域記作 Ω ， 邊界為 ?Ω，單元類型采用三節(jié)點三角形單元，則單元內(nèi)任一點處的場函數(shù)的近似值滿足式（6）：

式中：N為形函數(shù)矩陣，N=[N1N2N3]；vk=[v1v2v3]T，其中，v1、v2、v3為單元節(jié)點處的電勢值。

方程（4）對應(yīng)的等效積分強(qiáng)的形式表示如下：

利用分部積分進(jìn)行處理，離散時考慮邊界條件為完全齊次，則分部積分的降階項在邊界 ?Ω處為0，邊界條件的影響在方程求解時利用“乘大數(shù)法”進(jìn)行考慮?？刂品匠痰牡刃Хe分弱形式表示如下：

將式（6）代入式（8），得到電勢場控制方程的有限元形式表示如下：

通過式（9）可以獲得每一個單元內(nèi)的剛度矩陣；然后，通過集成得到總體的剛度矩陣，邊界條件采用“乘大數(shù)法”引入。有效電勢的衰減對于電勢場的影響體現(xiàn)在邊界條件的變化，而離散多態(tài)的處理方法能夠為黎曼狀態(tài)段內(nèi)電勢分布計算提供修正的邊界條件。隨著黎曼狀態(tài)段數(shù)的增加，離散形式的描述會不斷逼近原始的有效電勢衰減函數(shù)，與此相應(yīng)地需要付出更多次的穩(wěn)態(tài)求解。

取單元大小為b=h=0.4 m，陽極處電勢值為48 V，依據(jù)前述有限元形式（8）進(jìn)行求解，可以得到非對稱電極單元的電勢場分布，如圖4所示。

圖4中箭頭為各點處的電勢梯度的方向（場強(qiáng)方向），非對稱情況下電勢的分布與均勻電場差異很大，同時，與軸對稱模型的電勢分布也存在明顯區(qū)別。

對稱情況下，單元中存在2個陽極、2個陰極，電勢場分布如圖5所示，該情況下電勢梯度在單元中部區(qū)域，梯度方向與均勻電場差異較大。無論是對稱還是非對稱形式，電極單元的電勢場分布具有明顯的2維特性，電滲中電勢梯度決定了電滲流的大小及流動方向，因而采用2維電滲固結(jié)模型是較為合理的。

圖5 對稱電極單元的電勢梯度方向場Fig.5 Potential gradient direction field of symmetric electrode element

3 2維電滲固結(jié)的多態(tài)繼承方法

利用第2.2節(jié)黎曼求和的思路，對有效電勢衰減函數(shù)進(jìn)行離散多態(tài)處理后，可以得到多個黎曼狀態(tài)段；在任意一個黎曼狀態(tài)段內(nèi)，有效電勢Ve是一個常量，每一個狀態(tài)段都對應(yīng)一個電勢分布狀態(tài)。此時，在狀態(tài)段內(nèi)電滲固結(jié)方程引入中間變量進(jìn)行變換是完全成立的，但對于整個時間域上的求解，還需要考慮各個黎曼狀態(tài)段內(nèi)初始條件在時間維度上的先后順序。

3.1 多態(tài)繼承的基本思路

有效電勢衰減是一個相對于時間而言的概念，因此離散多態(tài)的處理也只是在整個時間域上進(jìn)行有限狀態(tài)段的劃分；電滲固結(jié)問題是一個初邊值問題，從時間維度上來看，電滲孔壓的非穩(wěn)態(tài)與有效電勢的衰減是同一個過程，在時間維度上具有統(tǒng)一性，對于有效電勢離散多態(tài)的分割實際上也是對電滲固結(jié)非穩(wěn)態(tài)的一個分割。

基于上述認(rèn)識，可以認(rèn)為多態(tài)處理后得到多個黎曼狀態(tài)段，第k狀態(tài)段必然是發(fā)生在第k+1個狀態(tài)段之前的。如果令每一個狀態(tài)段的總時間長度均為dT，那么在狀態(tài)段k內(nèi)，場狀態(tài)量在經(jīng)歷了 dT時間的變化之后，會得到該狀態(tài)段末時刻的場函數(shù)值ζk(x,y,dT) ；而對于下一狀態(tài)段 （k+1）來說，它的發(fā)生是在狀態(tài)段k的基礎(chǔ)之上的，自然 ζk(x,y,dT)就是k+1 狀態(tài)段內(nèi)非穩(wěn)態(tài)問題的初始條件 ζk+1(x,y,0)。若從第1個狀態(tài)段開始進(jìn)行各個狀態(tài)段的遍歷迭代計算，并記錄每一個黎曼狀態(tài)段的非穩(wěn)態(tài)結(jié)果，最終可以獲得整個電滲過程在全時間域上的解答。為了簡化計算，可以取各個黎曼狀態(tài)段的時間長度均為 dT，又因為中間變量場的邊界條件全程保持不變，這樣將對整個時間域上的求解轉(zhuǎn)化為了對單一黎曼狀態(tài)段上的反復(fù)迭代求解。

3.2 方法的有限元形式

由于該方法需要進(jìn)行反復(fù)的迭代計算，因此，采用有限元法進(jìn)行數(shù)值實現(xiàn)。點狀電極單元中2維電滲固結(jié)控制方程的建立依然承認(rèn)Esrig理論中的基本假設(shè)[15]，此外加入以下假設(shè)：

1）不考慮單元外部的孔隙水補(bǔ)給；

2）單元中滲流場以電滲流為主，非對稱單元中，異性邊界上孔隙水的流動主要為電滲流，方向為電勢梯度方向；

3）電極尺寸遠(yuǎn)小于單元幾何尺寸，將其理想化為“質(zhì)點”；

4）土體性質(zhì)為各向同性介質(zhì)，即kex=key=ke，khx=khy=kh。

水力滲流滿足達(dá)西定律，電滲流與水力滲流滿足線性疊加原理[16]。經(jīng)過多態(tài)處理后的，使得方程的求解域被限制在了子狀態(tài)段內(nèi)，子狀態(tài)段內(nèi)有效電勢為常量，在狀態(tài)段的時間范圍內(nèi)電勢場是與時間無關(guān)的，本狀態(tài)段內(nèi)的電勢分布記作Vk(x,y)。依據(jù)連續(xù)性原理，得到第k黎曼狀態(tài)段內(nèi)單元的2維電滲固結(jié)控制方程式（10）：

式中，下標(biāo)k表示物理量屬于第k狀態(tài)段。

由于在第k狀態(tài)段內(nèi)有效電勢值是一個常量，則該狀態(tài)段內(nèi)的電勢分布函數(shù)Vk(x,y)與時間無關(guān)，因此引入中間變量 ζk(x,y,t)是完全成立的，見式（11）：

由此得到控制方程（10）等效變換后的形式，見式（12）：

網(wǎng)格劃分及單元類型均與前述的電勢場分布計算保持一致，單元的位移模式見式（13）：

該問題為非穩(wěn)態(tài)問題，采用Galerkin法進(jìn)行空間項離散，將近似解 ζˉk(x,y,t)代入式（12）后，得到等效積分強(qiáng)形式，見式（14）：

采用分部積分進(jìn)行處理，在齊次邊界條件下，分部積分產(chǎn)生的降階積分項在 ?Ω處為0，得到中間變量方程對應(yīng)的等效積分弱形式，見式（15）：

將單元位移模式式（13）代入式（15）后，得到式（16）：

這里采用的是三節(jié)點三角形單元，所以單元形函數(shù)為關(guān)于x、y的一次項，可以直接積分得出依據(jù)形函數(shù)的具體形式得到B；當(dāng)積分點放在單元的重心處時，可以得到S的具體表述式如式（17）所示：

式中，A為有限單元的面積。

上述計算已經(jīng)完成了空間坐標(biāo)的半離散。為了對時間項進(jìn)行離散。采用有限差分形式進(jìn)行時間項的離散，考慮到迭代的穩(wěn)定性，這里采用后向差分形式，對式（16）進(jìn)行后向差分處理，黎曼狀態(tài)段內(nèi)時間步數(shù)為n，因為每一個狀態(tài)段的時間均為 dT，得到差分的時間步長為 Δt=dT/n，最終中間變量控制方程的有限元形式見式（18）：

式中， ζk(i-1) 表 示i-1 時 刻的場函數(shù)值向量， ζk(i)為i時刻的場函數(shù)值向量， Δt為差分的時間步長。

式（18）給出了單一狀態(tài)段內(nèi)中間變量控制方程的有限元形式，在邊界條件及初始條件下，可以完成求解。第k狀態(tài)段內(nèi)的初始條件需要繼承第k-1黎曼狀態(tài)段末時刻的值，即為式（19）：

式中： ζk(i=0) 為第k狀態(tài)段內(nèi)初始時刻時各節(jié)點的中間變量向量，從0開始計數(shù)； ζk-1(i=n) 表示第k-1狀態(tài)段末時刻的各節(jié)點處中間變量向量，從0開始計數(shù)，該狀態(tài)段內(nèi)時間步數(shù)為n， 末時刻i=n。

依據(jù)第3.2節(jié)的假設(shè)1）～3），保持陰極開孔進(jìn)行排水，而其他邊界不排水，得到非對稱單元中中間變量方程的邊界條件，見式（20）：

通過上述有限元形式及（初始）邊界條件，遍歷每一個狀態(tài)段可以獲得全時間域上的有效電勢衰減下的中間變量函數(shù)分布，然后通過式（13）求解出電滲孔壓值。

4 試驗驗證與算例分析

4.1 點狀電極單元模型試驗

為探究點狀非對稱單元的電滲固結(jié)特性，同時驗證提出算法的可靠性，進(jìn)行了非對稱點狀電極單元的電滲固結(jié)試驗研究，模型尺寸為 40 cm×40 cm，電滲電源電壓為48 V，陽極電極采用鐵電極，陰極開孔為排水邊界。設(shè)置了16個電勢測點（包含4個電極），位置及編號如圖6所示。圖6中：a0、a3、a12、a15分別為電極界面電勢探針，布置在距離電極表面約1 cm處，監(jiān)測有效電勢；孔壓計采用動態(tài)微型負(fù)孔隙水壓力計，型號為SCYG318，探頭直徑8 mm，量程-100～100 kРa，其距離上邊界、右側(cè)邊界均為6.65 cm，利用RS485模塊可讀取孔壓傳感器信號，孔壓采樣間隔60 s。試驗土體為重塑的粉質(zhì)黏土，初始土樣平均含水率約為29.4%，土樣粒徑分布如圖7所示，試驗裝置如圖8所示。

圖6 電滲模型及測點分布Fig.6 Electroosmotic model and distribution of measuring points

圖7 土樣粒徑分布Fig.7 Soil particle size distribution

圖8 非對稱點狀電極單元電滲試驗裝置Fig.8 Asymmetric punctiform electrode unit electroosmosis test

4.2 試驗實測結(jié)果與分析

模型通電約51 h后，孔壓計位置處土體裂縫的開展使得孔壓讀數(shù)為0，隨即停止了試驗。電滲過程中電流初期呈上升趨勢，這是由土體的含鹽量較高所導(dǎo)致的[8]。電滲結(jié)束時電流下降了約45%，試驗過程中總電源電壓保持48 V不變，則總模型的電阻在增大，如圖9所示。由圖9可知：初始時刻模型總電阻約為35.87 Ω，51 h電滲后模型總電阻增長到64.86 Ω，并且模型的電阻在后期增大較快，導(dǎo)致電滲后期電滲效率較低；排水速度隨著電滲的進(jìn)行逐步降低，在48 h后，1 h內(nèi)的排水量小于2.0 g。

圖9 電滲排水量、電流、電阻隨時間的變化曲線Fig.9 Curves of electroosmotic discharge, current and resistance with time

由于電滲排水加固法具有加固不均勻的特性，試驗后在模型陽極-陰極對角線上進(jìn)行了含水率測試，結(jié)果見表1。電滲后，陽極附近土體含水率下降幅度較大，加固效果較好，測點15附近含水率下降9.85%；陰極由于是排水邊界，在電勢梯度作用下，模型中的孔隙水從陽極向陰極發(fā)生定向流動，使得陰極附近土體的含水率下降幅度較小，陰極處土體含水率下降約2.98%，排水加固效果不明顯。

表1 土體含水率Tab.1 Soil moisture contents

電滲過程中通過測量金屬19個探針位置處的電勢以確定點狀電極單元的實際電勢分布；其中，a0、a3、a12、a15處的電勢值用于計算土體與電極之間因接觸所損失的電壓，從而得到真正施加在土體上用于電滲固結(jié)的有效電勢。圖10為界面電壓隨時間變化曲線。由圖10可知：陰極處的界面電壓變化幅度較小，總體呈下降趨勢，中后期基本保持不變，這主要是由于陰極作為排水邊界，土體中的水在陰極匯聚排出，使得陰極附近的土體含水率較高，土體電導(dǎo)率下降較小，從初始時刻的8.3 V下降到3.2 V。界面損失電壓的降低，對于電滲效率具有積極作用。距離陰極較近的2個陽極a3與a12，變化趨勢較為一致；界面電壓隨著電滲的進(jìn)行，增長幅度較大，a3從3.7 V增加到21.7 V，增加幅度為486.5%；a12從6.6 V變化到23.5 V，增長幅度為256.1%。距離陰極較遠(yuǎn)陽極的界面電壓a15總體呈下降趨勢，但下降幅度較小，從初始時刻的13.2 V下降到6.4 V，下降幅度約為51.5%。

圖10 界面電壓-時間變化曲線Fig.10 Curves of interface voltage with time change

將陰極與土體的接觸界面確定為零勢面，即可得到施加在土體上用于電滲作用的有效電勢（電壓），3個陽極處的有效電勢變化曲線如圖11所示。

圖11 有效電勢-時間變化曲線Fig.11 Curves of effective potential with time change

由圖11可知：電滲初始時刻，真正作用在土體上的電壓分別為陽極a3為36.0 V，陽極a12為33.1 V，陽極a15為26.5 V；電滲結(jié)束后，陽極a3為23.1 V，陽極a12為21.3 V，陽極a15為38.4 V。有效電壓的變化受到電極與土體接觸界面的影響（圖10）。土體在電場作用下失水，而受到滲透距離的影響，使得距離陰極（排水邊界）較近的陽極區(qū)域土體含水率下降幅度較大，導(dǎo)致該區(qū)域的電阻增長較大，界面損失電壓較高，有效電壓隨著電滲的進(jìn)行發(fā)生衰減；陽極a15附近排水路徑較遠(yuǎn)，土體電阻增加幅度較小，低于其他兩個陽極處，而電源總電壓保持不變，所以分擔(dān)的有效電壓是在增大的。點狀電極單元具有顯著的電勢2維效應(yīng)，并且電滲過程中土體受到物理、化學(xué)反應(yīng)的影響，使得有效電勢不斷發(fā)生變化，從而導(dǎo)致電勢分布時刻在變化，電滲的機(jī)理指出電滲固結(jié)過程中的負(fù)孔壓幅值受到電勢值的影響，如果忽略有效電勢的變化，那么電滲固結(jié)計算后的負(fù)孔壓值將被高估?？傮w上，有效電勢近似呈線性變化，通過離散的數(shù)據(jù)電可以擬合出試驗時間內(nèi)有效電勢的連續(xù)函數(shù)，進(jìn)而采用多態(tài)方式進(jìn)行處理。擬合得到的有效電勢變化連續(xù)函數(shù)為：

圖12為初始時刻實測電極單元電勢分布。圖13為電滲51 h后電極單元電勢分布。電滲的進(jìn)行伴隨著土體電學(xué)特性的變化，有效電勢發(fā)生變化，導(dǎo)致了點狀電極單元的電勢場分布發(fā)生了重分布現(xiàn)象。觀察初始與電滲結(jié)束的電勢分布，可以發(fā)現(xiàn)：陰極區(qū)附近土體電勢隨著電滲進(jìn)行逐步降低，而陽極附近電勢是在增加，兩者的殘差矩陣如圖14所示。2維單元的電勢變化主要沿著最遠(yuǎn)端陽極向陰極發(fā)生變化。

圖12 初始時刻實測電極單元電勢分布Fig.12 Initial measured electrode element potential distribution

圖13 電滲51 h后單元電勢分布Fig.13 Potential distribution after electroosmosis for 51 hours

圖14 電勢殘差矩陣分布Fig.14 Distribution of residual matrix of electric potential

4.3 算法驗證

為了驗證本文提出算法的可靠性及實用性，依據(jù)試驗中提取到的有效電勢變化曲線，依據(jù)實測的電勢變化曲線進(jìn)行電勢場及電滲固結(jié)孔壓的計算。土體電滲參數(shù)依據(jù)文獻(xiàn)[29]的研究成果，土體電滲參數(shù)見表2。電滲中有效電勢按照圖11中曲線確定；多態(tài)分段數(shù)量m確定為1、10和20段，當(dāng)狀態(tài)段數(shù)為m=1時，就是不考慮電勢變化下的計算結(jié)果，為了保證具有對比性，3種計算方案中均保證相同的時間步長，總時間取51 h，總時間步數(shù)為n=10 000，計算方案詳細(xì)參數(shù)見表3。單元尺寸為b=h=0.4 m，網(wǎng)格單元采用三節(jié)點的Delaunay三角形單元，模型初始孔壓值依據(jù)孔壓計實測數(shù)據(jù)進(jìn)行確定，初始孔壓值為1.465 kРa。

表2 土體電滲參數(shù)Tab.2 Soil electroosmotic parameters

表3 計算方案Tab.3 Computing schemes

采用Рython語言自行編程開發(fā)了РyEcFem庫，并進(jìn)行電勢多態(tài)的分段，以及多態(tài)下電勢穩(wěn)態(tài)分布和2維電滲固結(jié)的多態(tài)繼承的計算，本系統(tǒng)具有完整的后處理功能。狀態(tài)段數(shù)為10時的電極有效電勢的多態(tài)分割曲線如圖15所示。

圖15 有效電勢多態(tài)分割Fig.15 Effective potential polymorphism segmentation

利用實測的有效電勢數(shù)據(jù)，通過擬合構(gòu)造得到近似的連續(xù)函數(shù)，如圖11所示，具體表達(dá)式見式（21）。非對稱單元的3電極處的有效電勢變化接近一次函數(shù)，在獲得擬合函數(shù)后，利用多態(tài)方法進(jìn)行處理，具體如圖15所示。多態(tài)方法的優(yōu)勢在于，對函數(shù)形式?jīng)]有具體要求，在獲得擬合函數(shù)的表達(dá)式基礎(chǔ)上，均可以對其進(jìn)行多態(tài)的分割處理。

在電極有效電勢進(jìn)行分割的基礎(chǔ)上，進(jìn)行每一個狀態(tài)段的電勢分布計算，由于篇幅有限，這里僅給出狀態(tài)段數(shù)為10初始狀態(tài)與末時刻（51 h）的電勢殘差分布，如圖16所示。與電勢殘差矩陣相比，兩者變化具有相似的規(guī)律：較遠(yuǎn)處的陽極附近電勢值在增加，其余兩個陽極處電勢值下降，電勢分布下降區(qū)域主要分布在單元中部。這一結(jié)果反映出多態(tài)繼承算法下的電勢場能夠較好地符合實際情況。

圖16 有效電勢殘差矩陣分布（狀態(tài)段數(shù)為10）Fig.16 Distribution of residual matrix of electric potential (status segment number is 10)

采用微型孔壓傳感器進(jìn)行了電滲孔壓的監(jiān)測，受室內(nèi)模型尺寸的影響，隨著排水的進(jìn)行，土體裂縫在持續(xù)的開展，在49.6 h時孔壓計附近處土體裂縫貫通，使得孔壓計讀數(shù)回彈至0；此外，需要注意的是，電滲過程中的電化學(xué)反應(yīng)會對傳感器進(jìn)行腐蝕，陰極附近對于金屬的腐蝕強(qiáng)烈，腐蝕會使陰極處的孔壓計讀數(shù)在約10 h左右出現(xiàn)異常。

孔隙水壓力（簡稱為孔壓）實測值如圖17所示，在49.6 h達(dá)到-21.53 kРa，陽極15處有效電勢隨著電滲的進(jìn)行增加，如圖11所示。因此，對于理論計算來說，如果不考慮這一種變化，那么得到的孔壓幅值是偏小的，考慮電勢變化得到的孔壓更符合實際情況。電滲后多態(tài)繼承算法在分段數(shù)為10時的孔壓值達(dá)到-29.01 kРa，而不考慮衰減的結(jié)果為-14.11 kРa。

圖17 實測模型孔隙水壓力曲線Fig.17 Measured model pore water pressure curves

4.4 多態(tài)方法的特性

由于第4.3節(jié)的計算是依據(jù)試驗而進(jìn)行的，時間較短，土體電滲固結(jié)過程不夠充分，為了探究多態(tài)繼承算法數(shù)值結(jié)果的特性，設(shè)置電滲時間為1×108s，有效電勢變化仍然按照圖11來確定，其他參數(shù)保持不變，進(jìn)行狀態(tài)數(shù)段數(shù)為1、10、20下的計算并對比。圖18為狀態(tài)段數(shù)為20時的電滲過程中第1 000步、5 000步和8 000步的負(fù)孔壓分布云圖。

圖18 有效電勢衰減下孔壓分布云圖（狀態(tài)段數(shù)為20）Fig.18 Pore pressure distribution under effective potential attenuation (status segment number is 20)

對于點狀電極單元而言，由于電極的非對稱性，使得電勢分布及電滲孔壓分布與1維理論存在較大的區(qū)別；距離排水邊界越近，滲流路徑越短，土體孔隙水排出越容易，從而失水導(dǎo)致土體電導(dǎo)率發(fā)生變化，進(jìn)而導(dǎo)致界面接觸電壓改變，較遠(yuǎn)處陽極受到滲流路徑的影響，有效電勢在增大，使得該位置處的土體負(fù)孔壓值最大，加固效果最好，其他兩個陽極處次之。圖19為不考慮電勢衰減下的電滲固結(jié)孔壓分布。

圖19 不考慮有效電勢衰減下孔壓分布云圖Fig.19 Pore pressure distribution without considering effective potential attenuation

對比圖19與圖18的結(jié)果可知，不考慮有效電勢的影響得到的負(fù)孔壓幅值要高于衰減條件下的結(jié)果，說明出多態(tài)繼承算法能夠體現(xiàn)出電滲過程中有效電勢變化對電滲固結(jié)的影響。

圖20為不同位置處孔隙水壓力變化曲線。由圖20可知：距離陽極較近的點，電滲產(chǎn)生的負(fù)孔壓值較大，符合經(jīng)典電滲固結(jié)理論中電勢值決定負(fù)孔壓大小的原理；位于單元中部的點距離陽極較遠(yuǎn)，電勢值較低，因此負(fù)孔壓值要小于陽極處的?？紤]有效電勢變化下的計算結(jié)果負(fù)孔壓值要低于有效電勢不變的結(jié)果，即經(jīng)典理論中如果不考慮電勢變化，會高估電滲產(chǎn)生的負(fù)孔壓值。隨著分段數(shù)目的增多，由于繼承所帶來的突變現(xiàn)象逐漸減弱，曲線趨于平滑，高電勢區(qū)域繼承帶來的突變越明顯。因此，本文提出的算法能夠反映出有效電勢對于電滲的影響，并且計算逼近具有穩(wěn)定性。

圖20 不同位置處孔壓變化曲線Fig.20 Pore water pressure curves at different locations

電滲固結(jié)與水力排水固結(jié)不同，電滲排水過程會產(chǎn)生負(fù)孔壓，因此，孔壓的消散往往指的是從正的超孔壓發(fā)展為負(fù)孔壓，描述這個過程的固結(jié)度定義與普通水力排水固結(jié)中的固結(jié)度定義不同。固結(jié)度本質(zhì)反映的是實時孔壓的變化與孔壓變化最大幅度之間的比值，由于有效電勢的衰減變化，負(fù)孔壓并非在最后時刻達(dá)到極值，使得電滲中土體固結(jié)度的定義需要被修正，見式（22）：

式中，U為平均固結(jié)度，u0初 始時刻的孔壓值，u(x,y,t)max為電滲孔壓幅值的最大值。

式（22）分母表示總的孔壓變化范圍，即初始孔壓到負(fù)孔壓幅值最大值之間的距離；分子項反映實時孔壓變化的程度。通過這樣的定義可以反映出電勢變化對固結(jié)度的影響。點狀電極單元整體平均固結(jié)度及平均孔壓曲線如圖21所示。由圖21可知：在電滲后期固結(jié)度曲線會出現(xiàn)回彈現(xiàn)象，說明在之前的電勢分布下已經(jīng)達(dá)到負(fù)孔壓極值，由于有效電勢的衰減使得后續(xù)難以維持這一種較高的幅值狀態(tài)，會逐步降低回彈，土中的單元平均孔壓曲線反映出這一現(xiàn)象，最終固結(jié)度必定小于1。m=10時最后時刻的固結(jié)度為0.980，m=20時最終固結(jié)度衰減到0.979；而不考慮衰減時，最終時刻固結(jié)度為1，曲線單調(diào)增長。

圖21 單元平均固結(jié)度及平均孔壓曲線Fig.21 Unit average consolidation degree and average pore water pressure curve

圖22為單元平均孔壓與經(jīng)典結(jié)果（分段數(shù)為1時的計算結(jié)果）的偏差曲線。m=10與m=20這兩種情況下差別較小，與經(jīng)典結(jié)果相比，主要在電滲后期殘差較大，分別達(dá)到9.160 kPa與9.516 kPa，意味不考慮有效電勢變化的結(jié)果比多態(tài)結(jié)果孔壓高出9.160 kPa與9.516 kPa，如圖22（a）所示。多態(tài)算法可以有效地體現(xiàn)出有效電勢對電滲中孔壓的影響，針對有效電勢變化下的電滲固結(jié)求解是具有合理性的。狀態(tài)段數(shù)的增加導(dǎo)致逼近誤差的降低，為了體現(xiàn)狀態(tài)段數(shù)目對計算的影響，進(jìn)行了m=10與m=20下孔壓偏差研究，如圖22（b）所示。偏差表現(xiàn)出初始條件繼承帶來的突變現(xiàn)象，隨著時間呈鋸齒態(tài)增加，從10段增大到20段，計算結(jié)果的偏差最大值為0.353 kPa；曲線鋒與谷之間的差值隨著時間不斷降低，從初期的0.12 kPa降低到0.07 kPa。因此，狀態(tài)段數(shù)沒有必要設(shè)定太多，在曲線已經(jīng)較為平滑的情況下，再增加狀態(tài)數(shù)目對提升計算精度貢獻(xiàn)并不高。

圖22 單元平均孔壓與經(jīng)典結(jié)果的偏差曲線Fig.22 Deviation curve between unit average pore pressure and classical results

5 結(jié) 論

針對同性電極間距較大排布形式下的電滲固結(jié)計算問題，將實際模型抽象為點狀電極單元，有利于該問題的理論研究，并針對有效電勢衰減條件下的電滲固結(jié)計算，提出了多態(tài)繼承的算法，建立了點狀電極單元的2維電滲固結(jié)的控制方程，利用Galerkin法進(jìn)行了方程的離散，給出了算法完整的有限元格式，基于Python語言進(jìn)行了數(shù)值實現(xiàn)；同時，進(jìn)行了非對稱形式下點狀電極單元的電滲試驗，并與所提出算法的數(shù)值結(jié)果進(jìn)行了對比分析。主要結(jié)果如下：

1）異性電極間距與同性電極間距相差不大時，單元的電勢分布具有明顯的2維效應(yīng)，與1維電滲固結(jié)理論有較大區(qū)別；點狀電極單元能夠描述不同極性匹配下的電滲模型，有利于電滲固結(jié)理論的研究。

2）利用黎曼求和思路進(jìn)行有效電勢變化曲線的處理，使得在某一狀態(tài)段內(nèi)滿足中間變量引入的條件限制，并利用初始條件繼承的方式，將全時間域內(nèi)的求解轉(zhuǎn)化為子狀態(tài)段內(nèi)的迭代計算。

3）非對稱點狀電極單元的電滲試驗表明，單元整體上有效電勢是在降低的，而滲流路徑的長短會對有效電勢的衰減變化產(chǎn)生影響。距離陰極（排水邊界）較遠(yuǎn)的陽極處，有效電勢隨時間增長；而其余兩個陽極距離陰極距離相同且較近，有效電勢處于衰減狀態(tài)。電滲過程中單元內(nèi)電勢分布隨有效電勢的變化而變化，主要衰減區(qū)域處于單元兩個陽極對角線附近，陰極電勢衰減幅度較小。

4）多態(tài)繼承算法下的電滲固結(jié)計算結(jié)果能夠反映出有效電勢變化對電滲中負(fù)孔壓分布的影響。初始條件繼承時造成的孔壓突變幅度受到狀態(tài)段數(shù)的影響，越接近陽極，幅度越大。當(dāng)孔壓曲線較為平滑后，繼續(xù)增大狀態(tài)段數(shù)，數(shù)值解精度提升有限，不利于計算效率的提升。