西安地區(qū)地裂縫帶Q3原狀黃土流變特性試驗研究

鄧亞虹 ，李 麗 ，慕煥東 ，王 鵬，李飛霞

（1.長安大學(xué) 地質(zhì)工程系，陜西 西安 710054；2.長安大學(xué) 國土資源部巖土工程開放研究實驗室，陜西 西安 710054；3.長安大學(xué) 西部礦產(chǎn)資源與地質(zhì)工程教育部重點實驗室，陜西 西安 710054；4.甘肅省交通科學(xué)研究院有限公司，甘肅 蘭州 730050；5.甘肅省地礦局第三地質(zhì)礦產(chǎn)勘查院，甘肅 蘭州 730050）

1 引 言

地裂縫是一種危害嚴(yán)重的地質(zhì)災(zāi)害，在世界上許多國家普遍存在，尤以我國和美國最為嚴(yán)重。自20 世紀(jì)中、后期以來，我國的汾渭盆地已有56 個縣市、170 處發(fā)現(xiàn)地裂縫，總計約400 條，已成為我國乃至全世界地裂縫災(zāi)害最為發(fā)育、災(zāi)害最為嚴(yán)重的地區(qū)之一。而西安是汾渭盆地中地裂縫發(fā)育最具典型代表的地區(qū)，也是遭受地裂縫災(zāi)害最為嚴(yán)重的城市[1-3]。目前西安地裂縫造成的經(jīng)濟(jì)損失已超過50 億，嚴(yán)重制約著城市的規(guī)劃和發(fā)展。

目前，對地裂縫成因機(jī)制的研究已趨于成熟。美國是開展地裂縫研究最早的國家，構(gòu)造成因觀提出者Leonard[4]認(rèn)為，地裂縫及地面破裂與附近的地震活動有關(guān)。Lofgren[5-6]用豎向和橫向滲透應(yīng)力解釋了地面變形和地裂縫的形成。綜合成因的提出者Holzer[7]認(rèn)為，構(gòu)造條件和開采地下水是影響地裂縫形成和活動的兩個主要因素。我國對地裂縫的研究主要是以張家明的構(gòu)造成因說[2]、易學(xué)發(fā)的水成說[8]以及劉國昌、彭建兵等的綜合成因說[9-10]為代表。土流變特性的研究也早已開始，在Тердаги 著的《土力學(xué)原理》[11]一書中所列黏土的長期試驗數(shù)據(jù)表明，黏土表現(xiàn)出明顯的彈性后效特性。土的蠕變本構(gòu)是以室內(nèi)流變試驗為基礎(chǔ)，得到土體應(yīng)力-應(yīng)變隨時間的變化規(guī)律，可以定量表征土體的流變特性。對黃土蠕變本構(gòu)的研究是土流變力學(xué)理論中最基本也是最重要的組成部分，其選取原則是在能反映黃土流變規(guī)律的前提下盡量選擇最簡單的模型。最常用的模型形式就是元件模型，既直觀、簡單，又能全面反映流變介質(zhì)的各種流變特性。如夏才初等[12-14]提出了統(tǒng)一流變力學(xué)模型，涵蓋了介質(zhì)的所有流變性態(tài)。王艷婷[15]用五元件廣義的Kelvin 模型描述Q2黃土的黏彈性特性。李良權(quán)等[16]提出了一維、三維非線性黏彈塑性蠕變模型。但與地裂縫的流變特性相關(guān)的研究則幾乎處于空白狀態(tài)，考慮到汾渭盆地地裂縫多發(fā)育于黃土介質(zhì)中，探索隱伏地裂縫的破裂擴(kuò)展模式就必須考慮到黃土介質(zhì)的力學(xué)特性，這就需要通過室內(nèi)試驗手段來實現(xiàn)。已有大型物理模型試驗和數(shù)值模擬結(jié)果均顯示隱伏地裂縫的破裂擴(kuò)展具有反傾特征，與實際地裂縫剖面結(jié)構(gòu)揭示的傾向一致性不符[17-18]。而地裂縫又是一種緩變型地質(zhì)災(zāi)害，考慮其發(fā)育介質(zhì)的流變特性，通過室內(nèi)蠕變試驗方法來探索隱伏地裂縫的破裂擴(kuò)展模式及其力學(xué)機(jī)制應(yīng)該是更為合理的途徑。

本文以西安地區(qū)長安地裂縫（fc1）帶Q3原狀黃土為研究對象，采用分級加、卸載方式進(jìn)行室內(nèi)三軸流變試驗，得到了不同圍壓下地裂縫帶Q3原狀黃土的蠕變曲線及應(yīng)變速率與時間關(guān)系曲線。并基于統(tǒng)一流變力學(xué)模型理論，分析了不同圍壓、不同應(yīng)力水平下的蠕變曲線特征，建立了能夠全面描述各流變性態(tài)的三維蠕變本構(gòu)模型，對模型及參數(shù)進(jìn)行辨識，并建立了相應(yīng)的三維蠕變本構(gòu)方程。

2 試驗方案

試驗所用Q3原狀黃土取自西安地區(qū)長安地裂縫（fc1）金滹沱探槽東壁裂縫上盤，取樣點距主裂縫水平距離6.7 m，取土深度5.5 m，其基本物理力學(xué)性質(zhì)指標(biāo)見表1。

表1 基本物理力學(xué)性質(zhì)指標(biāo)Table 1 Basic physical and mechanical properties

試驗在CSS-2901TS 三軸流變試驗機(jī)上進(jìn)行，土樣為直徑為39.1 mm，高為80 mm 的圓柱形試樣[19]。采用同一含水率（按天然含水率進(jìn)行統(tǒng)一配置），固結(jié)壓力σ3分別取100、200、300 kPa。試驗控制條件為固結(jié)不排水，采用分級循環(huán)加、卸載方式，試驗加載速率為0.001 kPa/s，卸載速率為0.005 kPa/s。按照常規(guī)三軸壓縮試驗獲取的抗剪強(qiáng)度值，將應(yīng)力分為5～8 級進(jìn)行加載。對于每一級荷載設(shè)定相同的加、卸載穩(wěn)定標(biāo)準(zhǔn)，即24 h 內(nèi)加載蠕變量或卸載回彈蠕變量小于0.005 mm 或以恒定變形速率發(fā)展時即可認(rèn)為穩(wěn)定，然后卸載至圍壓狀態(tài)，觀測24 h 內(nèi)無滯后恢復(fù)時，再施加下一級荷載，依此類推逐級進(jìn)行。當(dāng)蠕變曲線出現(xiàn)等速蠕變階段后，勢必會出現(xiàn)加速蠕變階段，此時試樣的應(yīng)變速率急劇增大，試樣在短時間內(nèi)即達(dá)到破壞[20]。

3 試驗結(jié)果及分析

3.1 蠕變曲線

按照“坐標(biāo)平移法”對試驗數(shù)據(jù)進(jìn)行處理，得到了地裂縫帶Q3原狀黃土在100、200、300 kPa 圍壓下的蠕變加卸載曲線，如圖1 所示。

圖1 不同圍壓下試樣的蠕變加、卸載試驗曲線Fig.1 Creep loading and unloading curves of samples under different confining pressures

對比3 個圍壓下的蠕變曲線可知，同一圍壓下，應(yīng)力越大，試樣達(dá)到穩(wěn)定的時間越長，反之越短；不同圍壓的相同應(yīng)力下，試樣的蠕變變形量基本一致，且到達(dá)穩(wěn)定所需的時間也基本一致；圍壓越大，試樣達(dá)到破壞時所需應(yīng)力越大；蠕變曲線一旦出現(xiàn)等速蠕變階段，則必將出現(xiàn)加速蠕變階段，最終導(dǎo)致試樣破壞。相應(yīng)的試樣破壞形態(tài)如圖2 所示。

從圖2 可以看出，100、200 kPa 圍壓下試樣破壞時鼓脹型特征不明顯，300 kPa 圍壓下試樣破壞形態(tài)呈明顯鼓脹型；隨著圍壓的增大，試樣的變形量逐漸增大。

3.2 應(yīng)變速率與時間關(guān)系曲線

為了更好地分析蠕變特性，從相應(yīng)的蠕變曲線中選擇數(shù)條典型曲線進(jìn)行分析，繪制相應(yīng)的應(yīng)變速率與時間關(guān)系曲線，如圖3 所示。

圖2 不同圍壓下試樣的破壞形態(tài)Fig.2 Failure modes of samples under different confining pressures

圖3 不同圍壓時低應(yīng)力及高應(yīng)力下應(yīng)變速率與時間關(guān)系曲線Fig.3 Curves of strain rate vs.time under low and high stresses and different confining pressures

圖3(a)、3(b)、3(c)分別為100、200、300 kPa圍壓下低應(yīng)力及較高應(yīng)力（左圖）和高應(yīng)力（右圖）水平下的應(yīng)變速率與時間關(guān)系曲線。從圖可以看出，在低應(yīng)力及較高應(yīng)力水平下，應(yīng)變速率隨時間逐漸減小，最終趨于0，符合減速蠕變階段特征，且隨著應(yīng)力水平的增大，相應(yīng)的應(yīng)變速率也越大；在高應(yīng)力水平下，應(yīng)變速率隨時間先減小，后恒定，最后急劇增大，符合蠕變曲線經(jīng)歷減速蠕變、等速蠕變和加速蠕變3 個階段的過程；加速蠕變段最終應(yīng)變速率比初始應(yīng)變速率大，且隨著破壞應(yīng)力的增大，相應(yīng)的應(yīng)變速率也越大，破壞所需時間越短。

對比3 個圍壓下的應(yīng)變速率與時間關(guān)系曲線可知，無論圍壓大小，試樣在破壞之前的各級應(yīng)變速率均是由大到小變化，最終為0；一旦出現(xiàn)應(yīng)變速率持續(xù)相等的情況，最終必將出現(xiàn)應(yīng)變速率急劇變大直至試樣破壞；試樣最終破壞時的應(yīng)變速率總是大于加載時的初始應(yīng)變速率；隨著圍壓的增大，試樣出現(xiàn)加速蠕變即達(dá)到破壞所需應(yīng)力越大，應(yīng)變速率也越大。

4 蠕變本構(gòu)模型的建立

4.1 模型辨識

力學(xué)模型理論是根據(jù)不同情況將材料抽象成彈簧、黏壺和滑塊等元件組成的復(fù)雜體系，元件之間不同的組合代表材料的不同蠕變特性?；玖W(xué)模型是由這3 個元件的串并聯(lián)組成的，共形成7 個模型[21]。將其中4 個與時間有關(guān)的模型即黏彈性、黏彈塑性、黏性和黏塑性稱為基本流變力學(xué)模型[12]。為了研究與時間有關(guān)的變形形態(tài)，考慮材料一般均具有瞬時彈性變形，從而可以將彈簧與基本流變力學(xué)模型進(jìn)行串聯(lián)組合，共形成15 個流變力學(xué)模型。統(tǒng)一流變力學(xué)模型（如圖4 所示）是將彈簧與4 個基本流變力學(xué)模型同時串聯(lián)組成的模型，它是包含所有流變形態(tài)的最復(fù)雜最完整的一個模型[12-14]。

圖4 統(tǒng)一流變力學(xué)模型[14]Fig.4 Unified rheological model[14]

基于統(tǒng)一流變力學(xué)模型理論，結(jié)合此次試驗蠕變曲線對模型進(jìn)行辨識：

（1）在加載瞬間均有瞬時變形產(chǎn)生，說明具有瞬時彈性變形性態(tài)，即本構(gòu)模型中應(yīng)串聯(lián)有單獨的彈簧元件。

（2）在很低應(yīng)力下，蠕變曲線為減速蠕變類型，且能夠完全回彈，說明具有黏彈性形態(tài)。

（3）在低應(yīng)力及較高應(yīng)力水平下，蠕變試驗曲線均為減速蠕變，說明不具有完全黏性形態(tài)，即模型中沒有單獨的黏壺存在。

（4）在高應(yīng)力水平下，蠕變試驗曲線出現(xiàn)等速蠕變階段，且最終出現(xiàn)加速蠕變階段，說明低應(yīng)力下該元件不起作用，當(dāng)應(yīng)力達(dá)到一定值時才起作用，為黏塑性體。

（5）分析較高應(yīng)力下加卸載曲線，卸載后蠕變變形有且僅有部分回彈，則必有黏彈塑性體。

于是辨識出的本構(gòu)模型與統(tǒng)一流變力學(xué)模型相比僅少了一個黏性體，是由西原體和村山體串聯(lián)而成，命名為“改進(jìn)西原模型”，如圖5 所示。

圖5 改進(jìn)西原模型（一維）Fig.5 Improved Nishihara model(one-dimensional)

根據(jù)模型的特點，屈服應(yīng)力的大小不同，蠕變曲線表現(xiàn)出兩種形式，如圖6 所示。從圖可以看出，在低應(yīng)力水平下，無論屈服應(yīng)力大小，蠕變曲線均僅有減速蠕變階段，且卸載后變形完全回彈；在較高應(yīng)力水平下，圖6(a)的蠕變曲線具有減速蠕變和等速蠕變階段，卸載后變形部分回彈；圖6(b)的蠕變曲線僅具有減速蠕變階段，卸載后變形部分回彈；在高應(yīng)力水平下，無論屈服應(yīng)力大小，蠕變曲線均具有減速蠕變和等速蠕變階段。

圖6 改進(jìn)西原模型蠕變曲線特征Fig.6 Creep characteristics of improved Nishihara model

從前述試驗結(jié)果可以看出，蠕變曲線在低應(yīng)力及較高應(yīng)力水平下均表現(xiàn)為減速蠕變階段，符合圖6(b)情況，相應(yīng)的一維本構(gòu)方程為

式中：ε(t)為蠕變量；t為施加應(yīng)力的時間；σ0為施加的應(yīng)力；E0、E1、E2分別為彈性體、黏彈性體及黏彈塑性體中彈簧的彈性模量；η0、η1、η2分別為黏彈性體、黏塑性體及黏彈塑性體中黏壺的黏滯系數(shù)；σS1、σS2分別為黏塑性體和黏彈塑性體的屈服應(yīng)力。

考慮實際工程問題均為復(fù)雜的三維問題，進(jìn)行蠕變試驗研究及數(shù)值模擬時也是三維狀態(tài)，因此，需要將一維本構(gòu)方程推廣為三維形式。

采用德魯克-普拉格（D-P）（廣義米塞斯）屈服準(zhǔn)則，可以得到平面應(yīng)變狀態(tài)下的三維蠕變本構(gòu)模型[22-23]，如圖7 所示，相應(yīng)的本構(gòu)方程見式（2）。

圖7 改進(jìn)西元模型（三維）Fig.7 Improved Nishihara model(three-dimensional)

式中：σ1、σ2為施加的應(yīng)力；K為彈性體積模量；G0為彈性剪切模量；G1、G2分別為黏彈性及黏彈塑性剪切模量；H1、H2、H3分別為黏彈性、黏塑性和黏彈塑性體積模量；F為屈服函數(shù)，F(xiàn)=，其中，I1、J2分別為第一應(yīng)力不變量、第二應(yīng)力偏量不變量，α為與黏聚力c 有關(guān)的系數(shù)，k為與內(nèi)摩擦角φ 有關(guān)的系數(shù)，分別為

4.2 參數(shù)辨識

根據(jù)蠕變試驗數(shù)據(jù)，分別對一維和三維本構(gòu)方程各參數(shù)進(jìn)行一一辨識。

4.2.1 一維本構(gòu)方程參數(shù)辨識方法[13]

通過蠕變曲線上的瞬時變形 ε0，可以求得不同應(yīng)力水平下對應(yīng)的彈性模量 E0，即

通過低應(yīng)力水平即σ0≤σS1、σS2時，蠕變曲線達(dá)到穩(wěn)定時的最終蠕變量 εc(∞)確定黏彈性參數(shù) E1和η1，即按式（1）中低應(yīng)力對應(yīng)的本構(gòu)方程，當(dāng)，可得

在蠕變曲線上任取一點(εc(t),t)，可得

通過高應(yīng)力水平即σ0≥σS1、σS2時，蠕變曲線直線段斜率k1來確定黏塑性參數(shù)η2，即

通過較高應(yīng)力水平即σS2≤σ0≤σS1時，對應(yīng)的蠕變和回彈曲線及已辨識的 E1和η1求得黏彈塑性參數(shù)。按式（1）中較高應(yīng)力對應(yīng)的本構(gòu)方程，可得：

在蠕變加載狀態(tài)下，蠕變量 εc(t)與加載歷時t的關(guān)系為

在蠕變卸載狀態(tài)下，可恢復(fù)的蠕變量 εce(t)與卸載歷時t 的關(guān)系為

式（7）減式（8），可得卸載狀態(tài)下不可恢復(fù)的蠕變量 εcir(t)與卸載歷時t 的關(guān)系為

當(dāng) t→∞時，

于是可得

將式（13）代入式（11），在蠕變曲線上任取一點(εcir(t),t)代入，得

4.2.2 三維本構(gòu)方程參數(shù)辨識方法[22-23]

根據(jù)彈性剪切模量 G0、彈性體積模量K、彈性模量E 及泊松比μ 之間的關(guān)系可得，

當(dāng) t=0時，式（2）中高應(yīng)力對應(yīng)的本構(gòu)方程變?yōu)?/p>

當(dāng) t→∞時，式（2）中高應(yīng)力對應(yīng)的本構(gòu)方程變?yōu)?/p>

將式（17）、（18）代入式（2）中高應(yīng)力對應(yīng)的本構(gòu)方程中，可得

移項后得

對上式兩邊取對數(shù)：

根據(jù)高應(yīng)力下蠕變曲線直線段斜率k2，即式（2）中高應(yīng)力對應(yīng)的本構(gòu)方程與時間呈線性關(guān)系的黏塑性體的系數(shù)確定H2：

在蠕變加載狀態(tài)下，蠕變量 εc(t)與加載歷時t的關(guān)系為

在蠕變卸載狀態(tài)下，可恢復(fù)的蠕變量 εce(t)與卸載歷時t 的關(guān)系為

二者相減可得卸載狀態(tài)下不可恢復(fù)的蠕變量εcir(t)與卸載歷時t 的關(guān)系：

式中：F1、F2分別為加載時、卸載后的屈服函數(shù)。

當(dāng) t→∞時，式（25）、（27）分別化為

可得

將式（29）代入式（27），在蠕變曲線上任取一點(εcir(t),t)，得

依據(jù)上述步驟，選取適合的蠕變試驗曲線便可辨識出三維本構(gòu)方程的所有參數(shù)，見表2。

從表2 可以看出，隨著圍壓增大，試樣達(dá)到破壞所需的應(yīng)力越大；隨著應(yīng)力增大，除黏彈性參數(shù)不變外，其他參數(shù)均呈現(xiàn)先增大后減小的現(xiàn)象；相同應(yīng)力條件下，圍壓越大，相應(yīng)的參數(shù)越小。

為驗證參數(shù)的合理性，將由上述參數(shù)所得本構(gòu)方程繪制的曲線與相對應(yīng)的試驗數(shù)據(jù)點進(jìn)行對比分析，如圖8 所示。

從圖中可以看出，模型曲線與試驗數(shù)據(jù)點擬合的較好，說明通過辨識后選取的模型和確定的模型參數(shù)都較好。限于篇幅，在此僅給出了模型參數(shù)的辨識方法。

圖8 不同圍壓下模型曲線與試驗數(shù)據(jù)點對比圖Fig.8 Comparison between simulated curves and experimental data points under different confining pressures

5 結(jié) 論

（1）在天然含水率狀態(tài)下，西安地區(qū)地裂縫帶Q3原狀黃土的蠕變曲線表現(xiàn)為：加載瞬間均有瞬時變形產(chǎn)生，在低應(yīng)力水平下，蠕變曲線表現(xiàn)為減速蠕變階段，卸載時變形完全回彈。

（2）在較高應(yīng)力水平下，蠕變曲線表現(xiàn)為減速蠕變階段，卸載時變形部分回彈。

（3）在高應(yīng)力水平下，蠕變曲線表現(xiàn)為等速蠕變階段，且最終出現(xiàn)加速蠕變階段直至試樣破壞。

（4）通過模型曲線與試驗數(shù)據(jù)點的對比可知，所得本構(gòu)模型、參數(shù)均合理，能夠較好地描述西安地區(qū)地裂縫帶Q3原狀黃土的流變特性。

[1]李永善.西安地裂縫[M].北京:地震出版社,1986.

[2]張家明.西安地裂縫研究[M].西安:西北大學(xué)出版社,1990.

[3]李永善,等.西安地裂縫及渭河盆地活動斷層研究[M].北京:地震出版社,1992.

[4]LEONARD R J.An earth fissure in southern Arizona[J].Journal of Geology,1929,37(8):765－774.

[5]LOFGREN B E.Hydraulic stresses cause ground movement and fissures,Picacho,Arizona[J].Geological Society American Abstracts Programs,1978,10(3):113.

[6]LOFGREN B E.Earth cracks caused by horizontal stresses[J].Earth Observation System,1984,65:882－883.

[7]HOLZER T L,GABRYSCH R K.Effect of water-level recoveries on fault creep,Houston,Texas[J].Ground Water,1987,25(4):392－397.

[8]易學(xué)發(fā).西安市地面沉降及地裂縫成因的討論[J].地震,1984,12(6):50－54.YI Xue-fa.A discussion on the ground subsidence and the genesis of ground fissure in Xi’an city[J].Earthquake,1984,12(6):50－54.

[9]劉國昌.西安的地裂縫[J].西安地質(zhì)學(xué)院學(xué)報,1986,8(4):9－22.LIU Guo-chang.The ground fissures in Xi’an,Shaanxi[J].Journal of Xi 'an Institute of Geology,1986,8(4):9－22.

[10]彭建兵,等.渭河盆地活動斷裂與地質(zhì)災(zāi)害[M].西安:西北大學(xué)出版社,1992.

[11]SCOTT R F.土力學(xué)原理[M].張孟喜譯.武漢:華中科技大學(xué)出版社,2007.

[12]夏才初,王曉東,許崇幫,等.用統(tǒng)一流變力學(xué)模型理論辨識流變模型的方法和實例[J].巖石力學(xué)與工程學(xué)報,2008,27(8):1594－1600.XIA Cai-chu,WANG Xiao-dong,XU Chong-bang,et al.Method to identify rheological models by unified rheological model theory and case study[J].Chinese Journal of Rock Mechanics and Engineering,2008,27(8):1594－1600.

[13]夏才初,許崇幫,王曉東,等.統(tǒng)一流變力學(xué)模型參數(shù)的確定方法[J].巖石力學(xué)與工程學(xué)報,2009,28(2):425－432.XIA Cai-chu,XU Chong-bang,WANG Xiao-dong,et al.Method for parameters determination with unified rheological mechanical model[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(2):425－432.

[14]夏才初,孫鈞.統(tǒng)一流變力學(xué)模型的若干問題研究[C]//第六次全國巖石力學(xué)與工程學(xué)術(shù)大會論文集.武漢:中國學(xué)術(shù)期刊電子出版社,2000.

[15]王艷婷.黃土流變特性試驗分析及本構(gòu)模型的研究[D].西安:長安大學(xué),2006.

[16]李良權(quán),徐衛(wèi)亞,王偉.基于西原模型的非線性黏彈塑性流變模型[J].力學(xué)學(xué)報,2009,41(5):671－680.LI Liang-quan,XU Wei-ya,WANG Wei.A nonlinear viscoelastoplastic rheological model based on Nishihara’s model[J].Chinese Journal of Theoretical and Applied Mechanics,2009,41(5):671－680.

[17]彭建兵,陳立偉,黃強(qiáng)兵,等.地裂縫破裂擴(kuò)展的大型物理模擬試驗研究[J].地球物理學(xué)報,2008,51(6):1826－1834.PENG Jian-bing,CHEN Li-wei,HUANG Qiang-bing,et al.Large-scale physical simulative experiment on ground fissure expansion mechanism[J].Chinese Journal of Geophysics,2008,51(6):1826－1834.

[18]徐繼山,彭建兵,馬學(xué)軍,等.邢臺市隆堯地裂縫發(fā)育特征及成因分析[J].工程地質(zhì)學(xué)報,2012,20(2):160－169.XU Ji-shan,PENG Jian-bing,MA Xue-jun,et al.Characteristic and mechanism analysis of ground fissuresin Longyao,Xingtai[J].Journal of Engineering Geology,2012,20(2):160－169.

[19]南京水利科學(xué)研究院.SL237-1999 土工試驗規(guī)程[S].北京:中國水利水電出版社,1999.

[20]袁聚云,徐超,趙春風(fēng),等.土工試驗與原位測試[M].上海:同濟(jì)大學(xué)出版社,2004.

[21]周德培.流變力學(xué)原理及其在巖土工程中的應(yīng)用[M].成都:西南交通大學(xué)出版社,1995.

[22]王芝銀,李云鵬.巖體流變理論及其數(shù)值模擬[M].科學(xué)出版社,2008.

[23]孫鈞.巖土材料流變及其工程應(yīng)用[M].中國建筑工業(yè)出版社,1999.