非飽和巖土材料的固流轉化數(shù)值模擬

呂　謙，李兆華

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

0　引　言

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

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

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

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

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

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

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

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

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

1.1.1畢肖普有效應力

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

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

(1)

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

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

(2)

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

1.1.2PLASOL彈塑性模型

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

(3)

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

(4)

(5)

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

(6)

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

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

(8)

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

1.1.3水分特征曲線

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

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

(9)

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

(11)

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

圖2　非飽和巖土介質水分特征曲線圖

1.2　固流轉化準則

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

(12)

式中σ′和ε分別為有效應力張量和應變張量。

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

(13)

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

1.3　賓漢黏性模型

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

如果

(14)

否則

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

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

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

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

圖4　FEMLIP方法計算過程

2.1　控制方程

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

(16)

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

2.2　增量本構關系

(17)

(18)

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

(19)

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

(20)

(21)

2.3　FEMLIP模型公式

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

(22)

(23)

3　FEMLIP模擬和分析

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

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

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

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

表1　彈塑性參數(shù)

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

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

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

3.2　固流轉化的模擬和分析

如果固流轉化產生，且達到破壞時土柱由彈塑性固體形態(tài)轉為黏性流體形態(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下 點A、點B、點C的基質吸力變化曲線

圖7　在飽和滲透系數(shù)1×10-3 m/s和1×10-4 m/s下 點A、點B、點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)破壞前階段基質吸力云圖及(e)～(f)破壞后階段土體流動形態(tài)

4　結　論

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

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

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

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

[2]熊威,甘忠,許旭東,等.基于粘彈塑性模型的時效成形仿真[J].塑性工程學報,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期)《俄羅斯大陸架地質調查及礦產資源開發(fā)利用現(xiàn)狀》，作者：王淑玲，王銘晗，邵明娟，吳西順，張煒，賈凌霄；工作單位：中國地質圖書館·中國地質調查局地學文獻中心；頁碼：28～34，82。該文由中國地質調查項目“地學情報綜合研究與產品研發(fā)”資助(編號：121201015000150002)。

特此說明。

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