基于非線性損傷理論的改進(jìn)CVISC模型及其在FLAC3D中實(shí)現(xiàn)

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

(1.中國長江三峽集團(tuán)有限公司博士后工作站，北京 100038；2.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083；3.湘潭大學(xué)土木工程與力學(xué)學(xué)院/巖土力學(xué)與工程安全湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 湘潭 411105；4.甘肅省地質(zhì)環(huán)境監(jiān)測(cè)院，甘肅 蘭州 730050)

滑坡是斜坡巖土經(jīng)歷一定時(shí)間累進(jìn)性變形破壞的結(jié)果。因此滑坡巖土變形破壞具有不同程度的流變特征，基于流變模型研究滑坡形成過程、預(yù)測(cè)滑坡活動(dòng)趨勢(shì)一直是國內(nèi)外滑坡研究的熱點(diǎn)和難點(diǎn)問題之一[1-4]。當(dāng)滑坡以低速緩慢形式滑移或復(fù)活時(shí)，滑坡巖土變形破壞的流變特征尤為顯著[1,3,5-6]。

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

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

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

1 滑帶土非線性損傷流變本構(gòu)模型的建立

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

大量研究證實(shí)，當(dāng)應(yīng)力水平高于長期強(qiáng)度時(shí)，巖土體將發(fā)生加速蠕變破壞[20-22]。加速蠕變破壞的實(shí)質(zhì)是巖土內(nèi)部損傷量變到質(zhì)變的外在表現(xiàn)。這一過程可用非線性損傷黏塑性元件進(jìn)行描述。依照損傷力學(xué)理論，巖土破壞的有效應(yīng)力可定義為：

(1)

(2)

式中：D——巖土內(nèi)部與黏性變形相關(guān)的損傷變量；

σ——應(yīng)力；

t——流變時(shí)間；

n——與應(yīng)變速率有關(guān)的常數(shù)，可通過擬合試驗(yàn)數(shù)據(jù)來確定；

t*——巖土進(jìn)入非線性加速流變的起始時(shí)刻。

t

將含有損傷變量的黏性變形用非線性損傷牛頓體進(jìn)行刻畫。非線性損傷牛頓體與圣維南體并聯(lián)，形成一個(gè)非線性損傷黏塑性元件。將這一元件與CVISC流變模型串聯(lián)，可得到能反映蠕變?nèi)A段、改進(jìn)的CVISC非線性損傷流變模型。當(dāng)應(yīng)力小于長期強(qiáng)度時(shí)，非線性損傷黏塑性元件不起作用，模型退化為CVISC模型，可描述衰減和穩(wěn)定蠕變兩個(gè)階段；當(dāng)應(yīng)力大于長期強(qiáng)度時(shí)，非線性損傷黏塑性元件則反映加速蠕變階段應(yīng)變隨時(shí)間的變化關(guān)系。

一維應(yīng)力狀態(tài)下，改進(jìn)的CVISC非線性損傷流變模型如圖1所示。

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

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

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

(3)

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

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

ε，εP——應(yīng)變和塑性應(yīng)變；

σ∞——長期強(qiáng)度；

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

其余符號(hào)同前。

當(dāng)加載應(yīng)力σ>σ∞，t>t*時(shí)，模型各部分及損傷變量均起作用，則一維蠕變方程為：

(4)

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

要實(shí)現(xiàn)將改進(jìn)的CVISC非線性損傷流變模型應(yīng)用于FLAC3D，需將改進(jìn)的CVISC一維模型擴(kuò)展成三維模型的差分形式。遵循Perzyna[8]提出的類比原理，可進(jìn)行模型擴(kuò)展。

σ<σ∞時(shí)，非線性損傷黏塑性元件不起作用，模型退化為CVISC模型，其三維差分本構(gòu)方程已由相關(guān)文獻(xiàn)給出[25]。在此，僅討論σ>σ∞時(shí)非線性損傷黏塑性元件發(fā)揮作用的情形。此時(shí)，三維狀態(tài)下總應(yīng)變偏張量可寫為：

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

(5)

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

Maxwell體三維狀態(tài)下的偏應(yīng)力-應(yīng)變關(guān)系為：

(6)

式中：Sij——應(yīng)力偏張量；

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

類似地，可得Kelvin體的偏應(yīng)力-應(yīng)變關(guān)系為：

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

(7)

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

對(duì)于塑性元件：

(8)

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

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

δij——Kronecker符號(hào)；

g——服從Mohr-Coulomb屈服準(zhǔn)則的塑性勢(shì)函數(shù)；

λ*——僅在塑性流狀態(tài)為非零參數(shù)，其值由塑性屈服條件所確定。

對(duì)于非線性損傷黏塑性體部分，當(dāng)應(yīng)力大于長期強(qiáng)度時(shí)，應(yīng)變由非線性黏壺承擔(dān)，有：

(9)

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

對(duì)模型整體有：

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

(10)

K——體積模量；

evol——體積應(yīng)變率偏量。

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

(11)

(12)

(13)

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

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

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

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

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

(14)

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

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

模型中塑性流動(dòng)法則采用不相關(guān)聯(lián)的M-C流動(dòng)法則，當(dāng)屈服函數(shù)f<0時(shí)，需根據(jù)塑性應(yīng)變?cè)隽扛聭?yīng)力。

(15)

式中：

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

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

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

2 模型應(yīng)用與驗(yàn)證

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

2.1 泄流坡滑坡概況

泄流坡滑坡位于甘肅省舟曲縣白龍江下游約5 km處，發(fā)育在秦嶺東西向構(gòu)造帶的光蓋山—迭山蠕滑斷裂帶內(nèi)[26]?；履蟼?cè)邊界直接受蠕滑型活動(dòng)斷層控制，斷層走滑速率和擠壓速率分別為1.4 mm/a、3.7 mm/a[27]?；缕矫嫔铣书L舌狀，縱長約2.6 km，平均寬度約550 m，滑坡體積約7 150×104m3?；w物質(zhì)主要為黃土狀土以及灰?guī)r、炭質(zhì)千枚巖強(qiáng)風(fēng)化碎石土，滑帶物質(zhì)為炭質(zhì)板巖、千枚巖泥化后的黏性土(圖2)。泄流坡滑坡活動(dòng)歷史近百年，1961年9月舟曲小型地震后和1981年4月9日暴雨后的2次劇烈活動(dòng)，均堵斷白龍江。該滑坡活動(dòng)在時(shí)間上具有長期低速滑移、伴隨間歇性強(qiáng)烈活動(dòng)的特點(diǎn)，空間上具有分級(jí)分塊特征[28]。

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

2.2 滑帶土流變特征

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

2.3 改進(jìn)CVISC模型的計(jì)算參數(shù)獲取

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

圖3 正應(yīng)力400 kPa時(shí)不同角礫含量的泄流坡滑帶土蠕變曲線擬合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 基于改進(jìn)前后CVISC流變模型的滑坡活動(dòng)過程模擬

1960年2月3日舟曲5.25級(jí)地震促使滑坡變形加劇，1961年9月泄流坡滑坡中后部發(fā)生大規(guī)?？焖倩顒?dòng)[31]。為驗(yàn)證改進(jìn)CVISC非線性損傷流變本構(gòu)模型的有效性，基于改進(jìn)前后的CVISC流變模型，對(duì)泄流坡滑坡進(jìn)行地震工況下的三維數(shù)值模擬。采用擬靜力法模擬地震工況，舟曲地處地震烈度Ⅷ度區(qū)，取對(duì)應(yīng)的水平峰值加速度0.25g。

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

圖4 泄流坡滑坡三維數(shù)值模型及計(jì)算監(jiān)測(cè)點(diǎn)Fig.4 3D numerical model of the Xieliupo landslide and the location of monitoring points

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

據(jù)經(jīng)驗(yàn)值和滑帶土的流變參數(shù)初步確定。運(yùn)行模型后根據(jù)實(shí)際監(jiān)測(cè)數(shù)據(jù)對(duì)初始流變參數(shù)調(diào)參，最終確定符合實(shí)際條件的參數(shù)如表1所示。

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

表1 泄流坡滑坡流變計(jì)算參數(shù)

圖5 地震工況下計(jì)算監(jiān)測(cè)點(diǎn)速率-時(shí)步關(guān)系曲線Fig.5 Velocity curves of the monitoring points under the seismic condition

3 結(jié)論

(1)基于巖土蠕變破壞非線性特質(zhì)和內(nèi)部破壞不斷積累特征，采用非線性損傷力學(xué)理論建立流變本構(gòu)模型，較傳統(tǒng)流變模型對(duì)巖土蠕變實(shí)質(zhì)的刻畫更為合理。通過引入損傷變量，將含有損傷變量的牛頓體與圣維南體并聯(lián)，可以實(shí)現(xiàn)對(duì)巖土非線性損傷流變特性的刻畫。

(2)FLAC3D內(nèi)置的CVISC線性流變模型不能模擬巖土加速蠕變，串聯(lián)非線性損傷黏塑性元件后，改進(jìn)的CVISC模型能夠模擬應(yīng)力大于長期強(qiáng)度時(shí)的巖土加速蠕變。借助FLAC3D的開放接口，可以實(shí)現(xiàn)二次開發(fā)。

(3)甘肅泄流坡滑坡滑帶土在殘余狀態(tài)下的蠕變特征具有典型的非線性損傷流變特性。通過對(duì)滑帶殘余狀態(tài)下蠕變曲線的擬合，可獲取改進(jìn)CVISC模型的計(jì)算參數(shù)。基于改進(jìn)CVISC模型的模擬結(jié)果與滑坡實(shí)際基本一致，證實(shí)改進(jìn)CVISC模型具有較好的有效性。