呂　謙，李兆華

(1.中國礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院，北京 100083；2.深部巖土力學(xué)與地下工程國家重點(diǎn)實(shí)驗(yàn)室，北京 100083)

0　引　言

一般來說，非飽和巖土材料表現(xiàn)出固態(tài)的彈塑性力學(xué)特性，可以由各種現(xiàn)有的彈塑性本構(gòu)模型來描述[1]。然而，在破壞發(fā)生后，隨著動(dòng)能的突增，一些近飽和的松散巖土體也往往會(huì)表現(xiàn)出類似流體的力學(xué)特征。從微觀角度講，滑坡的發(fā)生過程可以解釋為巖土顆粒系統(tǒng)由“強(qiáng)接觸”向“弱接觸”轉(zhuǎn)變的過程，即抗剪強(qiáng)度突降而剪切應(yīng)變快速發(fā)展。為了實(shí)現(xiàn)滑坡從啟動(dòng)、擴(kuò)展到停滯全過程的模擬，許多學(xué)者針對非飽和巖土材料建立了一個(gè)考慮水力耦合作用的黏彈塑性固流轉(zhuǎn)化模型[2-3]。該模型既可以描述非飽和巖土介質(zhì)在破壞前的固體彈塑性力學(xué)特性，也可以描述破壞后的流體黏性力學(xué)特性，介于兩種狀態(tài)之間，引入了一個(gè)破壞準(zhǔn)則作為該模型固流轉(zhuǎn)化的判別準(zhǔn)則。該固流轉(zhuǎn)化統(tǒng)一模型詳細(xì)情況可參見文獻(xiàn)[3]。

現(xiàn)有的數(shù)值方法無法在同一個(gè)框架下準(zhǔn)確地處理固流兩個(gè)問題，于是一些學(xué)者將滑坡現(xiàn)象分成兩個(gè)不同的階段來進(jìn)行研究[4]。處理此類問題要求數(shù)值方法能夠準(zhǔn)確計(jì)算各種力學(xué)行為，尤其是滑動(dòng)過程中的大變形問題。為滿足上述要求，多種顆粒點(diǎn)法應(yīng)運(yùn)而生。物質(zhì)點(diǎn)有限元法(PFEM)[5]是一種使用mesh-less有限元插值函數(shù)的方法。

此外，物質(zhì)點(diǎn)法(MPM)吸取了拉格朗日方法和歐拉方法的長處，計(jì)算域被離散為一系列可自由移動(dòng)的物質(zhì)點(diǎn)，牛頓第二定律作為控制方程并用以計(jì)算固定網(wǎng)格上的速度場。這種方法較頻繁地應(yīng)用在滑坡模擬中[6]，但存在缺點(diǎn)。

確定了宏觀本構(gòu)模型后，只要流動(dòng)過程中沒有明顯的不連續(xù)現(xiàn)象，基于連續(xù)介質(zhì)力學(xué)的數(shù)值方法仍是有效的。本文選擇拉格朗日積分點(diǎn)有限元法(FEMLIP)對非飽和巖土介質(zhì)的固流轉(zhuǎn)化現(xiàn)象進(jìn)行研究，該方法由質(zhì)點(diǎn)網(wǎng)格法發(fā)展而來[7]，吸取拉格朗日和歐拉方法兩者的優(yōu)勢，既可以計(jì)算存儲(chǔ)彈塑性變量的問題又可以處理大變形問題[8-9]。

1　考慮水力耦合的非飽和巖土材料固流轉(zhuǎn)化統(tǒng)一模型

對于降雨誘發(fā)的滑坡，細(xì)粒土在破壞前表現(xiàn)出固體力學(xué)特性，而破壞后往往表現(xiàn)為流體的黏性流動(dòng)。由李兆華等[3]提出的固流轉(zhuǎn)化統(tǒng)一模型可用于這一類型的模擬。如圖1所示，這個(gè)三維模型可由一維的形式表現(xiàn)出來。圖1中“d2W”指二階功[9]。

圖1　非飽和巖土材料固流轉(zhuǎn)化三維模型的一維形式

如圖1所示，飽和或非飽和多孔介質(zhì)使用彈塑性模型來描述，水力耦合計(jì)算通過引入有效應(yīng)力和水分特征曲線來實(shí)現(xiàn)。需要指出的是，該模型使用二階功作為破壞準(zhǔn)則，可以描述分岔破壞理論中的擴(kuò)散破壞、應(yīng)變集中破壞和塑形極限破壞。

1.1　考慮水力耦合的彈塑性模型

1.1.1畢肖普有效應(yīng)力

對于非飽和多孔介質(zhì)，關(guān)于應(yīng)力變量選取的問題在國際范圍內(nèi)一直沒有定論。本文采用的固流統(tǒng)一模型使用畢肖普有效應(yīng)力來描述非飽和多孔介質(zhì)，見式(1)。

σ′=σ-uam+χ(ua-uw)m

(1)

式中：σ′、σ、ua和uw分別為有效應(yīng)力矢量、總應(yīng)力矢量、孔隙氣壓力和孔隙水壓力；χ為畢肖普參數(shù)，這里mT是一個(gè)標(biāo)量。mT=(1,1,1,0,0,0)，s=ua-uw是基質(zhì)吸力。畢肖普參數(shù)通常表達(dá)為χ=Sr，Sr是飽和度，但該表達(dá)式所反應(yīng)的是一個(gè)理想介質(zhì),而不是實(shí)際的多孔介質(zhì)。根據(jù)Alonso等的提議[10]，χ可以由aχ和nχ兩個(gè)參數(shù)確定，通過適當(dāng)?shù)囟x這兩個(gè)參數(shù)可以描述非飽和土浸潤過程中的塑形破壞現(xiàn)象。

新的畢肖普參數(shù)χ的公式，見式(2)。

(2)

式中：s為基質(zhì)吸力;Patm為大氣壓力。通過適當(dāng)定義aχ和nχ,畢肖普參數(shù)χ可以描述非飽和土的很多特性。

1.1.2PLASOL彈塑性模型

這里使用彈塑性模型PLASOL來模擬非飽和巖土材料破壞前的力學(xué)行為。該模型采用Van Eekelen準(zhǔn)則作為屈服面及塑性極限破壞準(zhǔn)則，VE屈服面接近于摩爾庫倫屈服面。該屈服面公式，見式(3)。

(3)

m=a(1+bsin(3β))n

(4)

(5)

式中：參數(shù)n用來控制屈服面形狀，根據(jù)Van Eekelen的建議將其默認(rèn)為-0.299;參數(shù)rc和re分別為三軸壓縮和拉伸試驗(yàn)簡化半徑，計(jì)算公式見式(6)。

(6)

式中φe為三軸拉伸路徑下的修正內(nèi)摩擦角。

區(qū)別于屈服面的不相關(guān)聯(lián)塑性流動(dòng)法則計(jì)算見式(8)。

(8)

式中：m′與屈服面公式中的m形式相同;φc和φe由壓縮、拉伸路徑下的剪脹角ψc和ψe替代。

1.1.3水分特征曲線

水分特征曲線(WRC)反應(yīng)了非飽和巖土介質(zhì)中基質(zhì)吸力和飽和度或含水量之間的水力關(guān)系。

如圖2所示，一個(gè)完整的水分特征曲線通常包括邊界干燥曲線、邊界浸潤曲線和掃描曲線。這里我們使用修正VanGenuchten-Mualem模型，該模型考慮了回滯效應(yīng)及孔隙率對飽和度的影響，其邊界曲線公式見(9)。

(9)

式中：下標(biāo)v為浸潤過程;d為干燥過程;Srres,Srsat和Srv分別為殘余飽和度、飽和狀態(tài)下的飽和度及當(dāng)前飽和度;av和nv為兩個(gè)參數(shù)，av定義如(10)所示;saev與n相關(guān)，定義如式(11)所示。

(11)

式中：saev0為參考孔隙率n0對應(yīng)的空氣進(jìn)入量值;λ為一個(gè)物理參數(shù)。如此，建立了一個(gè)基質(zhì)吸力及孔隙率相關(guān)的水力模型。

圖2　非飽和巖土介質(zhì)水分特征曲線圖

1.2　固流轉(zhuǎn)化準(zhǔn)則

根據(jù)分岔破壞理論，巖土材料的塑性流動(dòng)具有明顯不相關(guān)聯(lián)性。根據(jù)分岔破壞的定義，極小的擾動(dòng)荷載可能導(dǎo)致突然或不連續(xù)的反應(yīng)。這里我們選用該準(zhǔn)則作為巖土材料的破壞準(zhǔn)則和固流轉(zhuǎn)化準(zhǔn)則，其公式見式(12)。

(12)

式中σ′和ε分別為有效應(yīng)力張量和應(yīng)變張量。

二階功由有效應(yīng)力增量計(jì)算得來，以便分析土骨架的失穩(wěn)問題。如果對于所有加載方向，當(dāng)d2w>0,則材料處于穩(wěn)定狀態(tài)。否則在二階功小于等于零的加載方向上可能發(fā)生應(yīng)變集中破壞或擴(kuò)散破壞。式(12)所表述的局部二階功在各積分點(diǎn)上計(jì)算,并用以判斷材料的局部破壞情況，為了分析邊值問題，按式(13)計(jì)算整體二階功。

(13)

式中：ωi為積分點(diǎn)i的數(shù)值權(quán)重;Ji為雅可比矩陣行列式。式(13)求出的整體二階功,通常被用來判斷巖土體的整體破壞情況。

1.3　賓漢黏性模型

巖土材料在破壞后可以由單一閥值sy的賓漢黏性模型來描述。根據(jù)Duvaut和Lions及Balmforth和Craster的公式,三維簡化賓漢黏性模型表達(dá)，見式(14)。

如果

(14)

否則

如圖3所示，賓漢應(yīng)力閥值通常小于Van Eekelen彈性極限，而二階功在彈性極限范圍內(nèi)總是正值，故一般情況下二階功可認(rèn)為是巖土材料固流轉(zhuǎn)化的準(zhǔn)則。

2　基于拉格朗日積分點(diǎn)有限元的水力耦合方法

本研究使用拉格朗日積分點(diǎn)(FEMLIP)有限元法進(jìn)行數(shù)值分析。如圖4(a)、圖4(b)所示，計(jì)算域通過歐拉網(wǎng)格進(jìn)行離散化，而巖土材料采用拉格朗日物質(zhì)點(diǎn)進(jìn)行空間離散化。各個(gè)計(jì)算步結(jié)束時(shí)，物質(zhì)點(diǎn)根據(jù)插值所得的速度場更新位置以體現(xiàn)材料的變形，而網(wǎng)格本身不會(huì)發(fā)生變形，如圖4(d)所示。這樣可以處理大變形問題而避免網(wǎng)格畸變現(xiàn)象。

圖3　主應(yīng)力空間下的固流統(tǒng)一模型

圖4　FEMLIP方法計(jì)算過程

2.1　控制方程

非飽和巖土材料是一種可壓縮的三相介質(zhì)，即由巖土框架、水和空氣組成的固相、液相和氣相三相介質(zhì)。簡化后的控制方程由平衡方程和液相連續(xù)方程組成，平衡方程的表示形式見式(15)。液相連續(xù)方程由水的質(zhì)量守恒方程和廣義達(dá)西定律組成，見式(16)。

(16)

式中：S和V分別為面力f作用的邊界和單元體積;σ為柯西總應(yīng)力張量；b為重力引起的體應(yīng)力矢量；g為重力加速度；ρw為水密度；bw為孔隙水的體力矢量；uw為孔隙水壓力；k為非飽和巖土體的滲透率矩陣。

2.2　增量本構(gòu)關(guān)系

(17)

(18)

應(yīng)當(dāng)注意，為了簡化起見，假設(shè)當(dāng)Δt足夠小時(shí),χ是恒定的。含水量表示為θ=n·Sr，用修正Van Genuchten-Mualem模型建立含水量-基質(zhì)吸力關(guān)系，見式(19)。

(19)

組合式(19)和式(20)，矩陣關(guān)系如式(20)所示。

(20)

(21)

2.3　FEMLIP模型公式

利用線性和常量函數(shù)N和Nw分別插值速度和水壓力場，基于控制方程和式(20)，離散有限元方程見式(22)。

(22)

(23)

3　FEMLIP模擬和分析

本文給出了一個(gè)重力狀態(tài)下初始穩(wěn)定的非飽和土柱的數(shù)值分析。該土柱面積為6 m×6 m，初始基質(zhì)吸力為500 kPa,且均勻分布。土體用一系列物質(zhì)點(diǎn)離散化。所有邊界均可自由滑動(dòng)，土體頂部表面的降雨(持續(xù)79 h)以1.8 mm/h(5×10-7m/s)的恒定入滲水力邊界條件模擬。

系數(shù)為簡化起見，在FEMLIP模型中，邊坡的滲透系數(shù)假定為各向同性，非飽和土的滲透率k由k=kr·ks確定，ks為飽和土的滲透系數(shù)，kr為相對水力傳導(dǎo)系數(shù)。在本文中，使用Van Genuchten提出的公式確定kr。

3.1　吸力和滲透性的建模與分析

表1和表2為模擬所需要的參數(shù)，在表2中，ks為飽和土的初始滲透系數(shù)。使用式(22)控制非飽和土的滲透系數(shù)。水力梯度和非飽和土的滲透系數(shù)使得水從介質(zhì)的頂部向底部滲透。選擇邊坡表面(A)、中部(B)和底部(C)三個(gè)位置來觀察滲透效果，如圖5(a)、圖5(b)所示。需要指出，這次分析暫不考慮固流轉(zhuǎn)化。

表1　彈塑性參數(shù)

表2　液壓和水力耦合參數(shù)

圖5　不同水壓邊界條件下的非飽和土柱

圖6給出了點(diǎn)A、點(diǎn)B、點(diǎn)C在兩種不同初始滲透系數(shù)ks=1×10-3m/s和ks=1×10-4m/s基質(zhì)吸力變化情況。因?yàn)橛晁枰^長時(shí)間才能滲透至比較深的地下，所以表面的基質(zhì)吸力比底部小。土體表面的基質(zhì)吸力變化量取決于水量，然而，根據(jù)達(dá)西定律， 滲透地面的水量部分受非飽和土壤滲透率的控制。圖7為非飽和滲透系數(shù)相對于時(shí)間的變化。對于較高的飽和滲透系數(shù)ks=1×10-3m/s，使用式(22)確定初始非飽和滲透率為5×10-7m/s，近似等于入滲水量。這表明水可以完全滲透進(jìn)土體。而對于較低飽和滲透系數(shù)ks=1×10-4m/s，所計(jì)算的非飽和滲透率約等于6×10-8m/s，明顯小于水量。因此，在這種情況下，只有部分水滲入土體。

3.2　固流轉(zhuǎn)化的模擬和分析

如果固流轉(zhuǎn)化產(chǎn)生，且達(dá)到破壞時(shí)土柱由彈塑性固體形態(tài)轉(zhuǎn)為黏性流體形態(tài)。圖5(b)所示在長為5 m的E-F邊界，以入滲速度為10-6m/s(3.6 mm/h)下模擬固流發(fā)生的過程。最初和最終的有效黏聚力c′0為0 kPa、10 kPa、25 kPa和c′f為0 kPa、25 kPa、45k Pa。表3和表4列出了用于模擬的附加參數(shù)。

圖6　在飽和滲透系數(shù)1×10-3 m/s和1×10-4 m/s下 點(diǎn)A、點(diǎn)B、點(diǎn)C的基質(zhì)吸力變化曲線

圖7　在飽和滲透系數(shù)1×10-3 m/s和1×10-4 m/s下 點(diǎn)A、點(diǎn)B、點(diǎn)C的非飽和滲透系數(shù)曲線

表3　彈塑性參數(shù)

土體參數(shù)Eνφe0=φc0φef=φcfψe=ψcBpBcρηsy值50.351525100.010.02110015012

表4　水力耦合參數(shù)

圖8　有效凝聚力對破壞的影響

圖9　黏度系數(shù)對水平位移的影響

圖10　(a)～(d)破壞前階段基質(zhì)吸力云圖及(e)～(f)破壞后階段土體流動(dòng)形態(tài)

4　結(jié)　論

本研究將非飽和巖土材料固流轉(zhuǎn)化模型應(yīng)用到拉格朗日積分點(diǎn)有限元中，并提出一個(gè)新的可用于模擬高速滑坡或泥石流的FEMLIP模型。

為了驗(yàn)證該FEMLIP模型，首先對一個(gè)具有均勻分布初始基質(zhì)吸力的非飽和土柱進(jìn)行模擬，分析了不同滲透系數(shù)對基質(zhì)吸力變化的影響。之后對水力邊界條件下土體固流轉(zhuǎn)化現(xiàn)象進(jìn)行模擬，并分析了黏聚力和黏度系數(shù)對固流轉(zhuǎn)化和流動(dòng)速度的影響?；谘芯拷Y(jié)果， FEMLIP模型在描述非飽和土

體的固體和流體狀態(tài)行為方面是有效的。

[1]楊光華,姚捷,溫勇.考慮擬彈性塑性變形的土體彈塑性本構(gòu)模型[J].巖土工程學(xué)報(bào),2013(8):1496-1503.

[2]熊威,甘忠,許旭東,等.基于粘彈塑性模型的時(shí)效成形仿真[J].塑性工程學(xué)報(bào),2012(2):33-37.

[3]LI Z H,DUFOUR F,DARVE F.Hydro-elasto-plastic modeling with a solid/fluid transition[J].Computers and Geotechnics,2016,75:69-79.

[4]ALONSO E E,GENS A,DELAHYTE C.Influence of rainfall on the deformation and stability of a slope in over consolidated clays:a case study[J].Hydrogeology Journal,2003,11:174-192.

[5]ZHANG X,KRAB benhoft K,SHENG D,et al.Numerical simulation of a flow-like land slide using the particle finiteelement method[J].Computational Mechanics,2015,55(1):167-177.

[6]ABE K,SOGA K,BANDARA S.Material point method for coupled hydromechanical problems[J].Journal of Geotechnical and Geoenvironmental Engineering,2014,140(3)：1-15.

[7]MORESI L,DUFOUR F,MUHLHAUS H B.Mantle convection modeling with viscoelastic/brittle lithosphere:numerical methodology and platetectonic modeling[J].Pureand Applied Geophysics,2002,159(10):2335-2356.

[8]HARLOW F H.The particle-in-cell computing method for fluid dynamics[J].Computer Methods in Physics,1964,3:319-343.

[9]PRIME N,DUFOUR F,DARAVE F.Unified model for gromateria solid/fluid states and the transition in between[J].Journal of Engineering Mechanics,2014，140(6):682-694.

[10]ALONSO E E,PEREIRA J M,VAUNAT J,et al.Amicrostructurally based effective stress for unsaturated soils[J].Géotechnique,2010,60(12):913-925.

更正說明

我刊2016年第25卷第11期(總第231期)《俄羅斯大陸架地質(zhì)調(diào)查及礦產(chǎn)資源開發(fā)利用現(xiàn)狀》，作者：王淑玲，王銘晗，邵明娟，吳西順，張煒，賈凌霄；工作單位：中國地質(zhì)圖書館·中國地質(zhì)調(diào)查局地學(xué)文獻(xiàn)中心；頁碼：28～34，82。該文由中國地質(zhì)調(diào)查項(xiàng)目“地學(xué)情報(bào)綜合研究與產(chǎn)品研發(fā)”資助(編號：121201015000150002)。

特此說明。

中國礦業(yè)雜志社編輯部