互補算法在一維非線性固結求解中的應用

鄧岳保，謝康和

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

1 引 言

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

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

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

2 互補條件構造

2.1 互補問題及算法

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

上式等價的分量形式表達為

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

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

2.2 e-σ'關系曲線及互補條件

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

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

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

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

若令，

則有，

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

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

式中： α= av2av1。

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

若引入v=-f，可得，

據前文分析有，

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

3 一維非線性固結求解

3.1 固結微分方程推導

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

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

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

（1）連續性條件

dt時間內微元體內水量Q的變化為

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

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

（2）達西定律

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

（3）有效應力原理

根據有效應力原理有：

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

式（16）在推導過程中利用了一維固結過程中任一點豎向總應力σ不隨時間變化的條件。

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

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

3.2 有限差分

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

圖3 差分計算網格Fig.3 Gridding of finite difference method

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

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

相應地求解條件為

由式(18)可得：

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

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

圖4 不同時刻e-σ' 關系狀態Fig.4 Correlation of e-σ' at different phases

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

3.3 求解

根據式(8)，差分網格上的任意一點有，

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

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

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

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

4 驗 證

4.1 參數取值

計算過程中，參數取值情況如下：

地基參數：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；

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

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

4.2 與迭代法對比

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

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

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

對比計算結果發現，當兩種方法的參數取值及網格劃分均相同時，相同時刻的超靜孔隙水壓力結果完全一致。圖5所示為時間段數n =50時兩種方法得到的超孔隙水壓力結果對比。

4.3 與太沙基理論對比

根據α的物理意義可知：當α=1時，分段線性壓縮曲線變為一條直線。下文將這種特殊情況時的非線性計算結果與太沙基單向固結差分法計算結果（簡稱線性解）進行對比。

太沙基單向固結理論微分方程為

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

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

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

5 計算與分析

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

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

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

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

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

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

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

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

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

（2）eA對計算結果的影響

參數 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 互補模型算法優越性說明

計算過程中發現，相比于普通迭代法，互補模型算法優越性有兩點。

（1）計算效率

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

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

由圖可知：當T =10Δt時，計算程序只迭代了4次；當T =50Δt時，程序迭代了8次；隨著時間的推移，迭代次數向m靠近。總和不同時刻迭代次數可知，互補模型算法程序總的迭代次數比mn小。這種減小幅度在網格劃分密集時非常顯著。

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

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

6 結 論

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

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

（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一維固結理論研究綜述[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]. 巖土力學, 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]梁仕華, 齊添, 謝康和, 等. 超固結土的一維固結理論及其試驗研究[J]. 應用力學學報, 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滲流的飽和黏性土一維固結分析[J]. 巖土工程學報, 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]郭小明. 彈塑性接觸問題的研究及其應用[D]. 南京:東南大學, 2002.

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

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

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