互補算法在一維非線性固結(jié)求解中的應(yīng)用

鄧岳保，謝康和

（浙江大學(xué) 巖土工程研究所，杭州 310058）

1 引 言

由于二向、三向固結(jié)理論在指標(biāo)測定與計算求解方面存在諸多困難，Terzaghi單向固結(jié)理論在工程中仍被廣泛應(yīng)用。該理論建立在一系列簡化和假設(shè)基礎(chǔ)之上，如假定土體是均質(zhì)飽和的理想彈性材料、土層的壓縮和水的滲流只沿著豎向發(fā)生、滲透系數(shù)和壓縮系數(shù)是常數(shù)、大面積的外荷載瞬時施加等。這些假定是對實際情況的理想化，因而理論結(jié)果與實際情況常常存在差異。對此，國內(nèi)外學(xué)者通過對Terzaghi理論基本假設(shè)進行修正，得到考慮多種影響因素的固結(jié)理論，如考慮固結(jié)性狀的非線性、固結(jié)荷重隨時間改變等。目前，這些研究仍在進一步深入[1－6]。

筆者嘗試將互補算法引入到非線性固結(jié)問題研究中?；パa算法即計算互補問題的方法，是數(shù)學(xué)規(guī)劃的一個重要分支。研究者對其廣泛關(guān)注始于 20世紀(jì)60年代中期。40多年來，互補問題己發(fā)展成為數(shù)學(xué)、經(jīng)濟、工程等多個學(xué)科非常受歡迎的工具。其在工程領(lǐng)域中的應(yīng)用包括彈塑性力學(xué)問題、接觸力學(xué)問題、斷裂力學(xué)問題、潤滑問題、最優(yōu)控制問題及交通平衡問題等[7－10]。上述研究通過挖掘互補關(guān)系，可使得待求問題模型化為互補問題，從而最終歸結(jié)為利用互補算方法求解互補問題。應(yīng)用互補算法的優(yōu)越性在于，一方面使得所求問題有自然、清晰的數(shù)學(xué)描述、另一方面可運用其豐富而實用的理論去分析和求解問題[10]。

下文以一維非線性固結(jié)為例進行推導(dǎo)。首先將固結(jié)土體的壓縮特性曲線進行分段線性擬合，然后挖掘其中的互補條件，構(gòu)造互補模型，并將其應(yīng)用于固結(jié)微分方程的差分求解。最后，通過與迭代法解答和太沙基進行對比以驗證該法的合理性。

2 互補條件構(gòu)造

2.1 互補問題及算法

互補問題的數(shù)學(xué)描述是：給定函數(shù)F(x)，求矢量x，使?jié)M足如下的方程和不等式條件，

上式等價的分量形式表達(dá)為

式(1)中，當(dāng) F(x)是線性函數(shù)，即 F(x)=Mx+q（M為n階矩陣，q為n階向量）時，上述問題為線性互補問題，否則為非線性互補問題?；パa問題的名稱來自于式(1)中的第3個方程，稱為互補條件，即互補對( xi， Fi(x))中至少有一個變量為0。

對互補問題的深入研究，促進了互補算法的發(fā)展。目前，計算互補問題的方法主要包括Lemke算法、方程組類算法和內(nèi)點法等。有關(guān)這些算法的研究可參考文獻[8, 10－11]。下面將在e-σ′關(guān)系擬合曲線中尋找互補條件并構(gòu)造互補方程。

2.2 e-σ'關(guān)系曲線及互補條件

不論是太沙基的單向固結(jié)理論，還是目前規(guī)范中的沉降計算，均取壓縮曲線上的壓縮系數(shù)為常數(shù)，相當(dāng)于用直線擬合e-σ′關(guān)系曲線。然而實際上，根據(jù)常規(guī)的側(cè)限壓縮試驗或三軸試驗，不同類型土(尤其是高壓縮性土)的 e-σ′關(guān)系曲線均反映出非線性的特點。為此，本文用分段線性來進行逼近。出于簡化推導(dǎo)過程，選擇分段數(shù)為 2（該法可方便地推廣到多分段數(shù)）。

壓縮曲線分段擬合情況如圖1（分為OA和AB兩段所示），線性分段的分界點為A(eA，σA′)。擬合曲線可描述為

式中：av1、av2分別為OA段和AB段的壓縮系數(shù)，由曲線形式可知 av1≥ av2。

圖1 土體e-σ' 曲線及其擬合曲線Fig.1 Soil e-σ' curve and its fitting curve

若令，

則有，

控制變量λ的物理意義如圖1所示。當(dāng) σ′ ≤σA′時，λ=0，對應(yīng)e-σ′關(guān)系曲線第1段；當(dāng) σ′＞ σA′ 時，λ＞ 0，對應(yīng)e-σ′關(guān)系曲線第2段。至此，控制變量λ的引入，使得e-σ′關(guān)系曲線表達(dá)式統(tǒng)一于一個表達(dá)式。

將式(5)中σ′表達(dá)式代入式(4)并令左右相減得g，則：

式中： α= av2av1。

結(jié)合式(4)、(6)和(7)分析 f和λ之間的關(guān)系：當(dāng)σ′ ≤σA′ 時，λ=0，由 e-σ′曲線可知eA＜e，故f＜ 0；當(dāng) σ′＞ σA′ 時，λ＞0，由式(4)和式(7)的關(guān)系可得 f=0。

若引入v=-f，可得，

據(jù)前文分析有，

式(9)表明，v和λ之間為互補關(guān)系。這樣，e-σ′關(guān)系曲線中的互補條件已經(jīng)獲得。

3 一維非線性固結(jié)求解

3.1 固結(jié)微分方程推導(dǎo)

在如圖2所示厚度為H的均質(zhì)飽和土層上施加無限寬廣均布荷載 p，該荷載不隨時間變化，土中附加應(yīng)力沿深度均勻分布(即面積 abcd)，頂層面為排水邊界，底層面為不透水邊界。均質(zhì)土層的初始孔隙比e0、滲透系數(shù)k不隨時間變化。超孔隙水壓力u、孔隙比e以及有效應(yīng)力σ′均隨時間變化而變化。分析中其他假設(shè)條件與太沙基假定相同，僅壓縮系數(shù)為非常數(shù)。壓縮系數(shù)在不同的應(yīng)力階段有不同的值，見式(3)。

圖2 飽和軟土一維固結(jié)分析Fig.2 One-dimensional consolidation analysis of saturated soil

考察土層頂面以下z深度的微元體dxdydz在dt時間內(nèi)的變化。

（1）連續(xù)性條件

dt時間內(nèi)微元體內(nèi)水量Q的變化為

式中：q為單位時間內(nèi)流過單位水平橫截面的水量。dt時間內(nèi)微元體內(nèi)孔隙體積Vv的變化為

dt時間內(nèi)微元體內(nèi)水量的變化等于微元體內(nèi)孔隙體積的變化，即則：

（2）達(dá)西定律

式中：i為水頭梯度；h為超靜水頭；γw為水重度。結(jié)合式(12)和(13)，有：

（3）有效應(yīng)力原理

根據(jù)有效應(yīng)力原理有：

結(jié)合式(14)、(15)可得，

式（16）在推導(dǎo)過程中利用了一維固結(jié)過程中任一點豎向總應(yīng)力σ不隨時間變化的條件。

式(16)即為本文所要用到的一維非線性固結(jié)微分方程。該微分方程求解條件包括初始條件和邊界條件。

初始條件：t =0時，e =e0，σ′ =0；

3.2 有限差分

式(16)形式上與太沙基一維固結(jié)微分方程相同，故可采用類似的有限差分法求解。地基差分計算網(wǎng)格劃分如圖3所示。圖中地基深度方向分段數(shù)為m，時間分段數(shù)為n。時間和深度跨度分別為Δt和Δz，且Δz=H/m。

圖3 差分計算網(wǎng)格Fig.3 Gridding of finite difference method

對式(16)進行有限差分，

式中：i =1,2,…, m，j =1, 2,…。處理后得，

相應(yīng)地求解條件為

由式(18)可得：

式(19)表明，在單向固結(jié)情況下，對于地基中的任意一點，只要知道該點及其上下節(jié)點處的孔隙比和有效應(yīng)力（或超靜水壓力），就可以求得經(jīng)過一個時段Δt后相同位置的孔隙比。根據(jù)太沙基單向理論假設(shè)，若e與σ′線性相關(guān)，則可由各節(jié)點處的初始孔隙比和初始有效應(yīng)力，通過反復(fù)利用式(19)，逐點推算出后續(xù)任意時刻的孔隙比和有效應(yīng)力（或超靜水壓力），獲得解答。

然而，當(dāng)e與σ′并非線性相關(guān)，而是在不同應(yīng)力階段有不同比例系數(shù)時，情況則發(fā)生變化。對照圖4進行分析。當(dāng)荷載施加后，地基中產(chǎn)生與外荷載相等的超靜孔隙水壓力，地基中有效應(yīng)力為 0。隨著時間推移，超孔隙水逐漸消散，有效應(yīng)力逐漸增大。這一過程中，由于地基上層離排水面近，孔隙水消散迅速，有效應(yīng)力增大快，孔隙比與有效應(yīng)力關(guān)系率先進入e-σ′擬合曲線下一分段，而地基下層土體的孔隙比與有效應(yīng)力關(guān)系還處于 e-σ′擬合曲線第一分段。這樣，地基中存在一個界面，界面以上土體e-σ′曲線斜率為av2，界面以下土體e-σ′曲線斜率為av1。隨著時間的推移和超靜孔隙水的消散，該界面不斷向下移動。

圖4 不同時刻e-σ' 關(guān)系狀態(tài)Fig.4 Correlation of e-σ' at different phases

上述分析表明，若e與σ′為非線性關(guān)系，則在反復(fù)利用式(19)時，需要反復(fù)判斷e與σ′關(guān)系所處的階段。下文利用e-σ′關(guān)系曲線的互補關(guān)系解決上述問題。

3.3 求解

根據(jù)式(8)，差分網(wǎng)格上的任意一點有，

結(jié)合式(19)和(20)得：

這樣，當(dāng)初始孔隙比和有效應(yīng)力確定后，利用式(19)可獲得下一時段的孔隙比，然后利用式(21)求解該時段的控制變量和有效應(yīng)力，進而獲得全部解答。下文稱該法為互補模型法，其總的求解思路如下：

①給定參數(shù) m、β，可確定Δz、Δt，完成網(wǎng)格劃分；②由初始孔隙比e0和初始有效應(yīng)力σ0′，根據(jù)式(19)確定經(jīng)過Δt后的孔隙比e1；③將e1代入式(21)得q1，然后利用互補算法可確定Δt時刻的λ1，進而得到σ1′、u1；④根據(jù)Δt時刻計算結(jié)果和式(19)確定經(jīng)過2Δt后的孔隙比e2，用上一步中相同的方法確定σ2′和u2。⑤依此類推，可得經(jīng)過nΔt后的en、σn′和un。

求解過程中存在收斂性問題，具體計算時可通過改變β或調(diào)整Δt以達(dá)到收斂目的。

4 驗 證

4.1 參數(shù)取值

計算過程中，參數(shù)取值情況如下：

地基參數(shù)：H =1 m，k =1.0 m/d，γw=10 kN/m3，e0=1.0，eA=0.8，av1=1.0，α=av2/av1=0.8；

網(wǎng)格劃分：深度方向分段數(shù) m=10，時間分段數(shù)n可取任意正整數(shù)。時間和深度跨度分別為Δt和Δz ，且Δz=H/m =0.1 m。

系數(shù)β=0.3，由此確定Δt。另外，仿效文獻[12]，取荷載參數(shù)p =1，則可得相對于外荷載的無量綱超靜孔隙水壓力u。

4.2 與迭代法對比

用一般的迭代法與本文互補模型法進行對比。一般迭代法的計算步驟如下：

(1)確定初始量，包括初始孔隙比 e0、初始有效應(yīng)力σ0′和初始超靜孔隙水壓力u0；(2)由式(18)確定經(jīng)過Δt后的孔隙比 e1，再由式(2)運用判斷語句確定有效應(yīng)力σ1′，進而得到Δt時刻的超靜孔隙水壓力 u1；(3)根據(jù)Δt時刻計算結(jié)果和式(18)確定經(jīng)過2Δt后的孔隙比e2，用上一步中相同的方法確定有效應(yīng)力σ2′和超靜孔隙水壓力u2；(4)依此類推，可得經(jīng)過nΔt后的孔隙比en、σn′和un。

圖5 超靜孔隙水壓力u對比結(jié)果Fig.5 Result contrast of u of two methods

對比計算結(jié)果發(fā)現(xiàn)，當(dāng)兩種方法的參數(shù)取值及網(wǎng)格劃分均相同時，相同時刻的超靜孔隙水壓力結(jié)果完全一致。圖5所示為時間段數(shù)n =50時兩種方法得到的超孔隙水壓力結(jié)果對比。

4.3 與太沙基理論對比

根據(jù)α的物理意義可知：當(dāng)α=1時，分段線性壓縮曲線變?yōu)橐粭l直線。下文將這種特殊情況時的非線性計算結(jié)果與太沙基單向固結(jié)差分法計算結(jié)果（簡稱線性解）進行對比。

太沙基單向固結(jié)理論微分方程為

式中：av=av1= av2。求解條件如下，

初始條件：t =0時，u =u0；

太沙基單向固結(jié)理論差分解求解過程可參考文獻[12]。對比該法和本文方法的計算結(jié)果發(fā)現(xiàn)：當(dāng)參數(shù)取值及網(wǎng)格劃分相同時，相同時刻的超靜孔隙水壓力結(jié)果亦完全一致。這樣，兩種驗證結(jié)果均表明本文方法的合理性。

5 計算與分析

參數(shù)取值同上文驗證部分。限于篇幅，下文僅給出超靜孔隙水壓力計算結(jié)果，并對其影響因素進行簡要分析。最后對本文方法的優(yōu)越性給予說明。

5.1 不同時刻線性解與非線性解對比

參數(shù)取值同上，不同時段線性解答與本文非線性計算得到超靜孔隙水壓力結(jié)果對比情況如圖6。

圖6 不同時刻超靜孔隙水壓力u情況Fig.6 Solutions of u at different phases

由圖6可知，由于考慮了壓縮系數(shù)隨有效應(yīng)力增大而減小，非線性計算得到超靜孔隙水壓力消散快于線性計算結(jié)果，且隨時間推移，兩者相差越來越明顯。

5.2 非線性計算影響因素分析

（1）非線性程度對本文計算結(jié)果的影響

參數(shù)α反映壓縮曲線的非線性程度。改變α的取值，令其分別等于1.0、0.8和0.6。3種情況下，時間分段數(shù)n =50時，不同深度處超靜孔隙水壓力分布如圖7所示。由圖可知，隨著非線性程度的增加，即壓縮曲線上后一分段的壓縮系數(shù)變小，孔隙水壓力消散增快，固結(jié)速率加快。

圖7 非線性程度對超靜孔隙水壓力u計算影響Fig.7 Effect of nonlinear level on the solution of u

（2）eA對計算結(jié)果的影響

參數(shù) eA反映壓縮曲線上進入下一線性段的門檻值。eA越接近初始孔隙比e0，表明越容易進入壓縮曲線的下一分段。改變eA的取值，令其分別等于1.0、0.8和0.6。3種情況下，n =50時，不同深度處超靜孔隙水壓力分布如圖8所示。由圖可知，隨著eA的減小，即進入壓縮曲線后一分段的臨界孔隙比變小，孔隙水壓力消散減緩。

圖8 eA對超孔隙水壓力u計算影響Fig.8 Effect of eA on the solution of u

5.3 互補模型算法優(yōu)越性說明

計算過程中發(fā)現(xiàn)，相比于普通迭代法，互補模型算法優(yōu)越性有兩點。

（1）計算效率

普通迭代法計算過程中，判斷語句需要遍歷不同時刻不同深度的各個節(jié)點，即運行mn次判斷?；パa模型法迭代次數(shù)則有所減小，對照圖9進行說明。圖9所示為本文方法得到的不同時刻不同深度控制變量λ解答。

圖9 不同時刻不同深度控制變量λ 解答Fig.9 Solution of λ at different phases

由圖可知：當(dāng)T =10Δt時，計算程序只迭代了4次；當(dāng)T =50Δt時，程序迭代了8次；隨著時間的推移，迭代次數(shù)向m靠近?？偤筒煌瑫r刻迭代次數(shù)可知，互補模型算法程序總的迭代次數(shù)比mn小。這種減小幅度在網(wǎng)格劃分密集時非常顯著。

（2）控制變量求解的意義

控制變量λ解答一方面能反映迭代次數(shù)；另一方面還能反映不同時刻、不同深度地基土體所處的壓縮狀態(tài)。即：當(dāng)λ=0為~時，地基土體處于壓縮曲線的第1分段；λ=0時，則處于第2分段。仍對照圖 9進行說明。由圖可知：當(dāng) T =10Δt時，地基0.4 m以上土體壓縮狀態(tài)處于壓縮曲線的第2分段，該部分土體進入非線性壓縮狀態(tài)；當(dāng)T =50Δt時，地基0.8 m以上土體壓縮狀態(tài)處于壓縮曲線的第2分段；而當(dāng)T =100Δt時，地基土體全都處于壓縮曲線的第2分段。

6 結(jié) 論

（1）互補模型算法與普通迭代法得到的結(jié)果一致，但在計算效率方面，由于互補模型算法是以最優(yōu)收斂方向逼近真實解答，故其計算效率高，當(dāng)其應(yīng)用于復(fù)雜問題求解時在速度上具有優(yōu)勢。

（2）在線性計算程序中嵌入互補模型算法即可進行非線性分析，程序容易實現(xiàn)。

（3）互補模型算法中通過對控制變量的求解，可判斷地基土體所處的壓縮狀況。

[1]XIE K H, QI T, DONG Y Q. Nonlinear analytical solution for one-dimensional consolidation of soft soil under cyclic loading[J]. Journal of Zhejiang University:Science A, 2005, 7(8): 1358－1364.

[2]問延煦, 施建勇. Terzaghi一維固結(jié)理論研究綜述[J].西部探礦工程, 2003, (2): 1－4 WEN Yan-xu, SHI Jian-yong. Current research of Terzaghi one-dimensional consolidation theory[J]. West China Exploration Engineering, 2003, (2): 1－4.

[3]HSU T W, LU S C. Behaviour of one-dimensional consolidation under time-dependent loading[J]. Journal of Engineering Mechanics, ASCE, 2006, 132(4): 457－462.

[4]常林越, 王金昌, 朱向榮. 多級線性荷載下飽和軟黏土一維大應(yīng)變固結(jié)解析解[J]. 巖土力學(xué), 2009, 30(8):2343－2347.CHANG Lin-yue, WANG Jin-chang, ZHU Xiang-rong.An analytical solution of 1-D finite strain consolidation of saturated soft clay under multistep linear loading[J]. Rock and Soil Mechanics, 2009, 30(8): 2343－2347.

[5]梁仕華, 齊添, 謝康和, 等. 超固結(jié)土的一維固結(jié)理論及其試驗研究[J]. 應(yīng)用力學(xué)學(xué)報, 2009, 26(2): 268－274.LIANG Shi-hua, QI Tian, XIE Kang-he, et al.Experimental verification of one-dimensional linear over-consolidation theory[J]. Chinese Journal of Applied Mechanics, 2009, 26 (2): 268－274.

[6]鄂建, 陳剛, 孫愛榮. 考慮低速非Darcy滲流的飽和黏性土一維固結(jié)分析[J]. 巖土工程學(xué)報, 2009, 31(7):1114－1118.E Jian, CHEN Gang, SUN Ai-rong. One-dimensional consolidation of saturated cohesive soil considering non-Darcy flows[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(7): 1114－1118.

[7]LEMKE C E, HOWSON J T. Equilibrium points of bimatrix games[J]. Journal of the Society for Industrial and Applied Mathematics, 1964, 12: 413－423.

[8]FERRIS M C, PANG J S. Engineering and economic applications of complementarity problems[J]. SIAM Review, 1997, 39(4): 669－713.

[9]郭小明. 彈塑性接觸問題的研究及其應(yīng)用[D]. 南京:東南大學(xué), 2002.

[10]何素艷. 互補問題算法研究及其在力學(xué)中的應(yīng)用[D].大連: 大連理工大學(xué), 2003.

[11]陳寶林. 最優(yōu)化理論與方法(第二版)[M]. 北京: 清華大學(xué)出版社, 2005.

[12]李廣信. 高等土力學(xué)[M]. 北京: 清華大學(xué)出版社,2004.