基于非線性損傷理論的改進CVISC模型及其在FLAC3D中實現

蔣 樹，文寶萍，蔣秀姿，李瑞冬，趙 成

(1.中國長江三峽集團有限公司博士后工作站，北京 100038；2.中國地質大學(北京)水資源與環境學院，北京 100083；3.湘潭大學土木工程與力學學院/巖土力學與工程安全湖南省重點實驗室，湖南 湘潭 411105；4.甘肅省地質環境監測院，甘肅 蘭州 730050)

滑坡是斜坡巖土經歷一定時間累進性變形破壞的結果。因此滑坡巖土變形破壞具有不同程度的流變特征，基于流變模型研究滑坡形成過程、預測滑坡活動趨勢一直是國內外滑坡研究的熱點和難點問題之一[1-4]。當滑坡以低速緩慢形式滑移或復活時，滑坡巖土變形破壞的流變特征尤為顯著[1,3,5-6]。

巖土體流變本構模型構建有多種方式，包括基于試驗數據擬合的經驗模型[7]，基于流變力學的解析模型[8]和基于黏彈塑性元件的元件模型[9]等。其中，元件模型因其物理意義明確、簡單直觀，在國內外研究中應用最廣。元件模型研究中，將巖土流變特性看作是彈性、塑性和黏滯性疊加的結果，認為虎克彈性體、圣維南塑性體和牛頓黏性體三個基本元件能夠表征巖土的彈性、塑性和黏滯性狀，通過將三個三基本元件按不同方式進行組合，實現對巖土體各種流變行為的擬合及預測[10-11]。由于虎克彈性體和牛頓黏性體反映應力-應變速率的線性關系，所以這類元件的線性關系組合不能描述具有典型非線性特性的巖土蠕變破壞。于是，許多學者嘗試將這些元件進行非線性組合建立反映巖土蠕變破壞性狀的非線性流變本構模型。鄧榮貴等[12]提出了應力與蠕變加速度成正比的非牛頓流體粘滯阻尼元件，將其與傳統模型結合，實現對巖石加速蠕變的刻畫。曹樹剛等[13]用非牛頓黏滯體替代西原模型中與塑性體并聯的牛頓黏滯體，提出了能夠描述巖石蠕變三階段的改進西原模型。然而，巖土變形破壞過程是內部微破裂逐漸積累過程，這些模型不能體現對巖土蠕變本質的刻畫。損傷力學出現后，基于損傷力學原理構建能夠反映巖土時效變形過程中損傷積累、裂紋擴展特征的流變本構模型成為流變力學研究的熱點[14]。朱維申等[15]將周維垣等[16]針對節理巖體提出的損傷斷裂模型與西原流變模型結合，研究巖質邊坡變形破壞發展過程。肖洪天等[17]基于三峽船閘花崗巖的裂紋擴展試驗，提出了裂紋流變擴展計算公式，以此為基礎，建立了裂隙巖體損傷流變力學模型。周峙等[18]基于巴東組泥質粉砂巖的室內三軸試驗，從Mohr-Coulomb準則出發建立了粉砂巖變形破裂全過程的損傷軟化統計本構模型。然而，這些模型僅適合于裂隙巖質邊坡的流變特征分析。目前針對滑坡巖土的非線性損傷流變模型在國內外研究尚較少涉及。

FLAC3D(Fast Lagrangian Analysis of Continua in 3 Dimensions)是由美國Itasca公司開發的、基于連續介質快速拉格朗日算法的有限差分數值模擬軟件[19]，內置Mohr-Coulomb、Drucker-Prager等25種彈塑性本構模型和CVISC、Burger等8種流變模型。Itasca公司在FLAC3D中提供了用戶接口和所有本構模型的源代碼，方便用戶對軟件內置本構模型進行修改和二次開發。CVISC流變模型是由Burger流變模型與Mohr-Coulomb塑性模型組合而成的元件模型，可以描述剪切狀態下巖土流變特征，尤為適合模擬以剪切破壞為主的滑坡巖土蠕變行為。但是，該模型為基礎元件的線性組合模型，只能描述滑坡巖土的減速蠕變和等速蠕變，不能刻畫滑坡巖土的加速蠕變。若以CVISC模型為基礎，基于損傷力學理論對其進行改進，建立描述滑坡巖土蠕變三階段的非線性損傷流變模型；依照FLAC3D的代碼規則，對改進模型進行二次開發，便可實現對滑坡巖土蠕變破壞全過程的數值模擬。

滑帶是滑坡的控制單元，滑帶蠕變特性控制滑坡變形破壞特征。本文以具有長期緩慢活動、并伴隨間歇性劇烈活動特征的甘肅舟曲泄流坡滑坡為例，針對該滑坡滑帶土的流變特征，在FLAC3D內置的CVISC流變模型中引入非線性損傷黏塑性元件，構建可描述滑坡加速蠕變過程的非線性損傷流變本構模型，通過FLAC3D開放的用戶接口實現對本構模型的二次開發。通過對比改進前后CVISC模型對泄流坡滑坡的數值模擬結果，驗證模型的有效性。

1 滑帶土非線性損傷流變本構模型的建立

1.1 損傷變量與非線性黏塑性元件

大量研究證實，當應力水平高于長期強度時，巖土體將發生加速蠕變破壞[20-22]。加速蠕變破壞的實質是巖土內部損傷量變到質變的外在表現。這一過程可用非線性損傷黏塑性元件進行描述。依照損傷力學理論，巖土破壞的有效應力可定義為：

(1)

(2)

式中：D——巖土內部與黏性變形相關的損傷變量；

σ——應力；

t——流變時間；

n——與應變速率有關的常數，可通過擬合試驗數據來確定；

t*——巖土進入非線性加速流變的起始時刻。

t

將含有損傷變量的黏性變形用非線性損傷牛頓體進行刻畫。非線性損傷牛頓體與圣維南體并聯，形成一個非線性損傷黏塑性元件。將這一元件與CVISC流變模型串聯，可得到能反映蠕變三階段、改進的CVISC非線性損傷流變模型。當應力小于長期強度時，非線性損傷黏塑性元件不起作用，模型退化為CVISC模型，可描述衰減和穩定蠕變兩個階段；當應力大于長期強度時，非線性損傷黏塑性元件則反映加速蠕變階段應變隨時間的變化關系。

一維應力狀態下，改進的CVISC非線性損傷流變模型如圖1所示。

圖1 改進CVISC非線性流變模型示意圖Fig.1 Modified CVISC non-linear rheological model

其中，模型第四部分為非線性損傷黏塑性元件，其余同CVISC模型。

當加載應力σ<σ∞時，模型退化為CVISC模型；當加載應力σ≥σ∞、t≤t*時，模型中非線性損傷變量不起作用。根據疊加原理，一維蠕變方程為：

(3)

式中：EM，ηM——Maxwell彈性模量和黏滯系數；

EK，ηM——Kelvin彈性模量和黏滯系數；

ε，εP——應變和塑性應變；

σ∞——長期強度；

ηR——非線性損傷黏塑性元件黏度；

其余符號同前。

當加載應力σ>σ∞，t>t*時，模型各部分及損傷變量均起作用，則一維蠕變方程為：

(4)

1.2 三維模型的差分形式及其在FLAC3D中的嵌入

要實現將改進的CVISC非線性損傷流變模型應用于FLAC3D，需將改進的CVISC一維模型擴展成三維模型的差分形式。遵循Perzyna[8]提出的類比原理，可進行模型擴展。

σ<σ∞時，非線性損傷黏塑性元件不起作用，模型退化為CVISC模型，其三維差分本構方程已由相關文獻給出[25]。在此，僅討論σ>σ∞時非線性損傷黏塑性元件發揮作用的情形。此時，三維狀態下總應變偏張量可寫為：

eij=(eM)ij+(eK)ij+(eP)ij+(eR)ij

(5)

式中：eij——應變總偏張量，下標M、K、P、R分別代表Maxwell體、Kelvin體、Mohr-Coulomb體和非線性損傷黏塑性元件的應變偏張量。

Maxwell體三維狀態下的偏應力-應變關系為：

(6)

式中：Sij——應力偏張量；

GM，ηM——Maxwell體的剪切模量和黏滯系數。

類似地，可得Kelvin體的偏應力-應變關系為：

Sij=2ηK(eK)ij+2GK(eK)ij

(7)

式中：GK，ηK——Kelvin體的剪切模量和黏滯系數。

對于塑性元件：

(8)

式中：(eP)ij——塑性偏應變率；

(eP)vol——塑性體的體積應變率偏量；

δij——Kronecker符號；

g——服從Mohr-Coulomb屈服準則的塑性勢函數；

λ*——僅在塑性流狀態為非零參數，其值由塑性屈服條件所確定。

對于非線性損傷黏塑性體部分，當應力大于長期強度時，應變由非線性黏壺承擔，有：

(9)

式中：(eR)ij——非線性損傷黏塑性元件的應變速率偏量。

對模型整體有：

σ0=K(evol-(eP)vol)

(10)

K——體積模量；

evol——體積應變率偏量。

將式(5)～(7)、(9)、(10)寫為增量形式，聯立即可得改進CVISC模型的三維差分形式：

(11)

(12)

(13)

式中：(SN)ij，(SO)ij——新、舊應力偏張量；

Δeij，Δ(eP)ij，Δt——應變總偏張量、Morh-Coulomb體應變偏張量以及時間的增量形式；

(eK,O)ij——Kelvin體應變偏張量的老值。

同理，式 (10)球應力的差分形式為：

(σN)0=(σO)0+K(Δevol-Δ(eP)vol)

(14)

式中：(σN)0，(σO)0——應力球張量變化率的新老值；

Δevol，Δ(eP)vol——體積應變率偏量、塑性體體積應變率偏量的增量形式。

模型中塑性流動法則采用不相關聯的M-C流動法則，當屈服函數f<0時，需根據塑性應變增量更新應力。

(15)

式中：

通過式 (11)和式 (14)計算的僅考慮黏彈性變形的新應力偏量和應力值；

(SN)i，(σN)0——考慮塑性部分的新偏應力值和應力值。

將推導得到的非線性損傷流變本構模型三維差分形式、應力更新及修正公式，利用FLAC3D軟件提供的本構模型二次開發程序接口，采用C++語言在Visual studio 2005平臺上對CVISC模型源代碼進行修改。修改數據項包括：初始化材料參數和關鍵求解函數，每一時步均調用一次求解函數，通過重載函數，根據子單元狀態進行塑性判斷與修正，并計算得到新的應力值，進而求得不平衡力、節點速率和節點位移。程序文件編寫完成后，將自定義本構模型代碼編譯成動態鏈接庫文件，在FLAC3D軟件中調用該文件即可應用自定義本構模型。

2 模型應用與驗證

泄流坡滑坡是我國著名的巨型低速滑坡，其活動特征具有典型的流變特性，因此將改進CVISC非線性損傷流變模型應用于該滑坡，以驗證模型的有效性。

2.1 泄流坡滑坡概況

泄流坡滑坡位于甘肅省舟曲縣白龍江下游約5 km處，發育在秦嶺東西向構造帶的光蓋山—迭山蠕滑斷裂帶內[26]?；履蟼冗吔缰苯邮苋浠突顒訑鄬涌刂疲瑪鄬幼呋俾屎蛿D壓速率分別為1.4 mm/a、3.7 mm/a[27]。滑坡平面上呈長舌狀，縱長約2.6 km，平均寬度約550 m，滑坡體積約7 150×104m3?；w物質主要為黃土狀土以及灰巖、炭質千枚巖強風化碎石土，滑帶物質為炭質板巖、千枚巖泥化后的黏性土(圖2)。泄流坡滑坡活動歷史近百年，1961年9月舟曲小型地震后和1981年4月9日暴雨后的2次劇烈活動，均堵斷白龍江。該滑坡活動在時間上具有長期低速滑移、伴隨間歇性強烈活動的特點，空間上具有分級分塊特征[28]。

圖2 泄流坡滑坡簡化剖面圖Fig.2 Cross section of the numerical model of the Xiliupo landslide

2.2 滑帶土流變特征

該滑坡的長期活動特性，表明其滑帶強度已降至殘余狀態，因而其緩慢持續活動特征受殘余狀態下滑帶土的流變行為控制[29]。蔣秀姿[30]對泄流坡滑坡滑帶土殘余狀態下的蠕變行為進行了系統研究，發現當剪應力小于殘余強度時，剪應變經過一段時間的減速增長后趨于定值，表現為衰減蠕變特征；當剪應力大于殘余強度時，剪應變經過減速蠕變后進入加速蠕變階段，直至蠕變破壞。據此，滑帶土中微裂隙在加速蠕變過程中形成、并不斷擴展，損傷積累直至蠕變破壞。因此，該滑坡滑帶土在殘余狀態下的蠕變行為具有非線性損傷流變性質。

2.3 改進CVISC模型的計算參數獲取

盡管遵循相同本構模型的滑坡具有相似的活動模式，但是模型參數刻畫著各個滑坡間行為的差異。如前所述，滑帶蠕變行為控制低速滑坡的活動特征?；诖?，依據泄流坡滑帶在殘余狀態下的蠕變曲線，擬合該滑坡的改進CVISC模型，從而獲取模型計算參數。對于應力水平低于或超過殘余強度的改進CVISC模型，分別采用未進入和進入加速蠕變階段的曲線分別擬合。圖3為不同顆粒級配下泄流坡滑帶土試樣的典型蠕變曲線與擬合的改進CVISC模型蠕變曲線。從圖3可以看出，不同應力水平下改進的CVISC蠕變曲線與試驗曲線擬合良好，由此獲取的滑帶模型計算參數列于表1。

圖3 正應力400 kPa時不同角礫含量的泄流坡滑帶土蠕變曲線擬合Fig.3 Fitting of the creep curves of the Xieliupo slip zone soil with different gravel content under the normal stress of 400 kPa

2.4 基于改進前后CVISC流變模型的滑坡活動過程模擬

1960年2月3日舟曲5.25級地震促使滑坡變形加劇，1961年9月泄流坡滑坡中后部發生大規?？焖倩顒覽31]。為驗證改進CVISC非線性損傷流變本構模型的有效性，基于改進前后的CVISC流變模型，對泄流坡滑坡進行地震工況下的三維數值模擬。采用擬靜力法模擬地震工況，舟曲地處地震烈度Ⅷ度區，取對應的水平峰值加速度0.25g。

滑坡三維模型依據實測工程地質圖和勘探資料建立(圖4)。在斷層位置設置界面單元，斷層上盤設定速率邊界條件，模擬斷層作用。在地表設置8個計算監測點JC1-JC8(圖4)。

圖4 泄流坡滑坡三維數值模型及計算監測點Fig.4 3D numerical model of the Xieliupo landslide and the location of monitoring points

巖土體的基本物理力學參數均來源于室內試驗測定，滑帶流變力學參數通過前述改進CVISC模型與試驗曲線擬合獲得，每條曲線都可擬合出1組流變參數，得到流變參數隨角礫含量變化的關系，再將滑帶土實際角礫含量13.75%代入關系式中來確定初始流變參數。黃土狀土以及風化碎石土層的流變力學參數則依

據經驗值和滑帶土的流變參數初步確定。運行模型后根據實際監測數據對初始流變參數調參，最終確定符合實際條件的參數如表1所示。

基于改進前后CVISC模型的數值模擬結果均顯示，施加地震荷載前，各監測點活動速率恒定，滑坡整體處于穩定蠕變狀態(圖5)。施加地震荷載后，基于CVISC流變模型的模擬結果顯示，各監測點的位移速率短期增加后，很快趨于恒定，表明地震作用造成滑坡活動速率加快，但并未出現加速蠕變，這與滑坡活動歷史不符(圖5)?；诟倪M的非線性損傷CVISC模型的模擬結果顯示，滑坡中下部3個監測點JC2、JC3和JC4的位移速率增大一定幅度后趨于勻速發展，而坡腳處JC1點和中上部JC5、JC6和JC7點的位移速率經過前期緩慢增加后，呈現急劇增長趨勢，量值達施加地震荷載前的2～3倍，反映滑坡中上部和坡腳處出現了局部大規模加速蠕變破壞特征，這與滑坡曾出現大規模分塊滑移的歷史基本一致。如此，證實了基于非線性損傷理論的改進CVISC模型具有較好的有效性。

表1 泄流坡滑坡流變計算參數

圖5 地震工況下計算監測點速率-時步關系曲線Fig.5 Velocity curves of the monitoring points under the seismic condition

3 結論

(1)基于巖土蠕變破壞非線性特質和內部破壞不斷積累特征，采用非線性損傷力學理論建立流變本構模型，較傳統流變模型對巖土蠕變實質的刻畫更為合理。通過引入損傷變量，將含有損傷變量的牛頓體與圣維南體并聯，可以實現對巖土非線性損傷流變特性的刻畫。

(2)FLAC3D內置的CVISC線性流變模型不能模擬巖土加速蠕變，串聯非線性損傷黏塑性元件后，改進的CVISC模型能夠模擬應力大于長期強度時的巖土加速蠕變。借助FLAC3D的開放接口，可以實現二次開發。

(3)甘肅泄流坡滑坡滑帶土在殘余狀態下的蠕變特征具有典型的非線性損傷流變特性。通過對滑帶殘余狀態下蠕變曲線的擬合，可獲取改進CVISC模型的計算參數?；诟倪MCVISC模型的模擬結果與滑坡實際基本一致，證實改進CVISC模型具有較好的有效性。