鄧岳保，謝康和

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

1 引 言

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

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

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

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

2.1 互補(bǔ)問(wèn)題及算法

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

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

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

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

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

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

壓縮曲線分段擬合情況如圖1（分為OA和AB兩段所示），線性分段的分界點(diǎn)為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′時(shí)，λ=0，對(duì)應(yīng)e-σ′關(guān)系曲線第1段；當(dāng) σ′＞ σA′ 時(shí)，λ＞ 0，對(duì)應(yīng)e-σ′關(guān)系曲線第2段。至此，控制變量λ的引入，使得e-σ′關(guān)系曲線表達(dá)式統(tǒng)一于一個(gè)表達(dá)式。

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

式中： α= av2av1。

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

若引入v=-f，可得，

據(jù)前文分析有，

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

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

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

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

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

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

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

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

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

dt時(shí)間內(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)過(guò)程中利用了一維固結(jié)過(guò)程中任一點(diǎn)豎向總應(yīng)力σ不隨時(shí)間變化的條件。

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

初始條件：t =0時(shí)，e =e0，σ′ =0；

3.2 有限差分

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

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

對(duì)式(16)進(jìn)行有限差分，

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

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

由式(18)可得：

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

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

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

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

3.3 求解

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

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

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

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

求解過(guò)程中存在收斂性問(wèn)題，具體計(jì)算時(shí)可通過(guò)改變?chǔ)禄蛘{(diào)整Δt以達(dá)到收斂目的。

4 驗(yàn) 證

4.1 參數(shù)取值

計(jì)算過(guò)程中，參數(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í)間分段數(shù)n可取任意正整數(shù)。時(shí)間和深度跨度分別為Δt和Δz ，且Δz=H/m =0.1 m。

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

4.2 與迭代法對(duì)比

用一般的迭代法與本文互補(bǔ)模型法進(jìn)行對(duì)比。一般迭代法的計(jì)算步驟如下：

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

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

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

4.3 與太沙基理論對(duì)比

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

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

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

初始條件：t =0時(shí)，u =u0；

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

5 計(jì)算與分析

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

5.1 不同時(shí)刻線性解與非線性解對(duì)比

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

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

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

5.2 非線性計(jì)算影響因素分析

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

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

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

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

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

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

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

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

（1）計(jì)算效率

普通迭代法計(jì)算過(guò)程中，判斷語(yǔ)句需要遍歷不同時(shí)刻不同深度的各個(gè)節(jié)點(diǎn)，即運(yùn)行mn次判斷。互補(bǔ)模型法迭代次數(shù)則有所減小，對(duì)照?qǐng)D9進(jìn)行說(shuō)明。圖9所示為本文方法得到的不同時(shí)刻不同深度控制變量λ解答。

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

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

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

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

6 結(jié) 論

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

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

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

[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]問(wèn)延煦, 施建勇. 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]常林越, 王金昌, 朱向榮. 多級(jí)線性荷載下飽和軟黏土一維大應(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é)理論及其試驗(yàn)研究[J]. 應(yīng)用力學(xué)學(xué)報(bào), 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]鄂建, 陳剛, 孫愛(ài)榮. 考慮低速非Darcy滲流的飽和黏性土一維固結(jié)分析[J]. 巖土工程學(xué)報(bào), 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]郭小明. 彈塑性接觸問(wèn)題的研究及其應(yīng)用[D]. 南京:東南大學(xué), 2002.

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

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

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