周鳳璽 ， 趙 曼， 應(yīng) 賽

(1. 蘭州理工大學(xué)土木工程學(xué)院， 甘肅 蘭州 730050； 2. 西部土木工程防災(zāi)減災(zāi)教育部工程研究中心， 甘肅 蘭州 730050； 3. 蘭州理工大學(xué)技術(shù)工程學(xué)院， 甘肅 蘭州 730050； 4. 長江師范學(xué)院土木建筑工程學(xué)院， 重慶 408100)

0 引言

非飽和土在自然界中廣泛存在， 尤其是在干旱與半干旱地區(qū). 由于受氣候條件的影響， 存在著若干種具有特殊性質(zhì)的土類， 如膨脹土、 黃土與殘積土等特殊土， 它們均具有非飽和土的基本特性. 在工程實(shí)踐中， 如在路堤及土壩工程的建設(shè)過程中， 通常涉及到膨脹土、 濕陷性土的變形分析以及邊坡穩(wěn)定與滑坡等問題. 但由于土體通常處于非飽和狀態(tài)， 孔隙水和孔隙氣共存， 使得非飽和土的力學(xué)特性與飽和土相比較顯得非常復(fù)雜. 以往大都通過室內(nèi)直剪試驗(yàn)與三軸試驗(yàn)， 來獲得非飽和土的強(qiáng)度參數(shù)以及應(yīng)力-應(yīng)變曲線等宏觀力學(xué)性質(zhì)， 或通過離散單元法， 從細(xì)觀角度分析土顆粒的運(yùn)動(dòng)規(guī)律以及孔隙水的作用力. 其中， 三軸剪切試驗(yàn)的優(yōu)點(diǎn)在于可以控制圍壓， 應(yīng)力狀態(tài)、 應(yīng)變測(cè)量方便可靠， 能夠反映試樣受力變形直到破壞的全過程， 但是， 室內(nèi)三軸剪切試驗(yàn)成本高、 制樣復(fù)雜， 技術(shù)要求高， 且耗時(shí)長. 相比之下， 三維離散元數(shù)值模擬不僅能通過控制顆粒級(jí)配、 孔隙率等生成數(shù)值試樣， 通過借助標(biāo)定細(xì)觀參數(shù)的方式考慮含水率、 剪切強(qiáng)度等因素， 還具備成本低、 效率高、 易建模等明顯的優(yōu)勢(shì).

離散單元法是從細(xì)觀尺度分析顆粒材料力學(xué)性質(zhì)以及解決非連續(xù)介質(zhì)力學(xué)問題的數(shù)值計(jì)算方法. 自從Cundall等[1-2]提出離散單元法后， 利用離散單元法對(duì)土體力學(xué)的研究已經(jīng)取得一系列成果[3-7]. Liu等[8]通過離散單元法對(duì)非飽和土在各向同性壓縮和雙軸剪切過程中的破壞機(jī)理進(jìn)行數(shù)值模擬， 確定顆粒間不同毛細(xì)黏聚力下的屈服應(yīng)力； Jiang等[9]應(yīng)用二維離散元法研究非飽和土的抗剪強(qiáng)度并提出兩個(gè)能夠考慮大吸力范圍和顆粒級(jí)配影響的抗剪強(qiáng)度函數(shù)； 蔣明鏡等[10]應(yīng)用三維離散元法， 考慮毛細(xì)水的作用， 研究吸力對(duì)密實(shí)和松散顆粒材料強(qiáng)度的影響； 在此基礎(chǔ)上， 許自立[11]通過三維離散元法建立非飽和土三軸剪切試驗(yàn)和單軸拉伸試驗(yàn)的離散元計(jì)算模型， 研究非飽和土的強(qiáng)度問題.

在非飽和土離散元數(shù)值模擬的過程中， 大多采用顆粒流軟件PFC中內(nèi)嵌的線性平行黏結(jié)模型或線性接觸黏結(jié)模型來模擬土中的基質(zhì)吸力， 這兩種模型適用于模擬黏結(jié)材料的力學(xué)行為. 但是， 在程序運(yùn)行的過程中， 需要人工不斷地調(diào)節(jié)黏結(jié)參數(shù)來近似模擬土中的吸力， 其本身不符合實(shí)際中非飽和土的基質(zhì)吸力. 對(duì)于蔣明鏡等[9-10]所提出的毛細(xì)水模型， 是將毛細(xì)水作用考慮入顆粒間作用力， 在PFC3D中編寫fish語句， 每當(dāng)有新的顆粒接觸時(shí)， 將新產(chǎn)生的接觸模型變?yōu)槊?xì)水模型. 然而， 對(duì)于這種用戶自定義的接觸模型， 需要基于第三方編譯平臺(tái)編寫準(zhǔn)確且符合實(shí)際的程序代碼， 其技術(shù)難度高， 且耗時(shí)長、 難推廣.

本文基于顆粒流理論， 借助離散元數(shù)值模擬工具PFC， 綜合考慮非飽和土的孔隙水壓力、 含水率等因素， 基于Hill接觸模型， 建立非飽和土的顆粒流模型， 并從細(xì)觀角度深入分析非飽和土的力學(xué)特性， 旨在對(duì)利用數(shù)值方法高效解決非飽和土強(qiáng)度相關(guān)問題進(jìn)行探索.

1 離散元模型

1.1 Hill材料

Hill材料被定義為顆粒集合體， 且Hill接觸模型存在于所有的顆粒-顆粒間接觸處， 所組成的Hill材料類似于不飽和的顆粒材料. 首先， Hill接觸模型可能只存在于顆粒-顆粒間的接觸處， 可將水分添加到Hill材料中， 如果水分狀態(tài)是濕的， 則顆粒間存在液橋， 如果水分狀態(tài)干燥， 則不存在液橋. 其次， Hill接觸模型提供了一個(gè)微小的、 非線性彈性(無張力)和摩擦界面， 具有壓縮表面相互作用力的行為， 并可攜帶吸力， 故采用Hill接觸模型來考慮濕顆粒材料.

Hill模型是在液橋模型的啟發(fā)下建立的， 這種模型機(jī)械地捕捉顆粒間分散的液體對(duì)顆粒系統(tǒng)宏觀響應(yīng)的影響. 只有當(dāng)水分狀態(tài)為濕時(shí)， 才會(huì)出現(xiàn)吸力Fm， 這是一種在法線方向的吸引力. 當(dāng)顆粒接觸時(shí)最大， 并且呈指數(shù)形式衰減， 直到接觸間隙達(dá)到一個(gè)臨界值2scr時(shí)， 液橋破裂， 此時(shí)水分消失吸力減小為零. 在不飽和的顆粒材料中， 考慮到與表面張力有關(guān)的吸力， 即在顆粒間接觸處保持孔隙水的作用， 從而考慮了水分效應(yīng)在Hill材料中的作用. 通過以下方法可在選定的顆粒之間加入吸力來模擬Hill材料中的水分. 每個(gè)Hill接觸都有水分狀態(tài)， 可以是濕的， 也可以是干的， 當(dāng)顆粒-顆粒間產(chǎn)生新的接觸時(shí)， 被指定為干的Hill接觸模型， 在模擬的過程中， 通過調(diào)用函數(shù)@him_makeWet(ψ,gm)， 可以將水分隨時(shí)加入到Hill材料中， 其中ψ為吸力，gm為粒間距離. 在此函數(shù)被調(diào)用后， 具有吸力的Hill接觸模型將存在于彼此的粒間距離內(nèi)的所有顆粒之間.

(a)顆粒-顆粒接觸示意圖 (b)Hill’s法向和切向方向流變接觸模型圖1 濕的Hill接觸模型的力學(xué)行為和流變特性Fig.1 Behavior and rheological components of the wet hill contact model

1.2 力-位移定律

Hill接觸模型的力-位移定律：

Fc=Fs+Fm;Mc≡0

(1)

式中：Fs為顆粒間的表面相互作用力；Fm為吸力. 表面相互作用力由赫茲力和阻尼力組成：

Fs=Fh+Fd

(2)

另外， 表面相互作用力Fs可分解為法向力和切向力：

(3)

吸力作用在法線方向：

Fm=-Fmnc

(4)

式中：Fm>0為張力(吸力)且Fm≥0.

2 非飽和土室內(nèi)三軸剪切試驗(yàn)的離散元模擬

2.1原型試樣顆粒級(jí)配

本文使用的試驗(yàn)數(shù)據(jù)來自室內(nèi)非飽和重塑黃土三軸剪切試驗(yàn)[12]， 該試驗(yàn)采用FYS30型應(yīng)力-應(yīng)變控制式非飽和三軸儀， 進(jìn)行非飽和黃土等吸力三軸剪切試驗(yàn). 試驗(yàn)中選用新疆西部伊犁地區(qū)重塑黃土， 天然含水量為4%～16.1%， 主要以粉土為主(含量高達(dá)81.5%)， 含有少量的砂粒(0.5%)和黏粒(18%)， 其基本物理性質(zhì)見表1.

表1 伊犁黃土試樣基本物理性質(zhì)

2.2 離散元模型

建立三軸剪切試驗(yàn)的顆粒流模型可分為試樣生成、 預(yù)固結(jié)和豎向加載3個(gè)階段. 以室內(nèi)非飽和土重塑黃土三軸剪切試驗(yàn)為基礎(chǔ)， 進(jìn)行離散元模擬. 模型邊界由圓柱筒和上下兩面墻體組成， 墻體尺寸與室內(nèi)試驗(yàn)中所制備的土樣大小相一致， 直徑39.1 mm， 高80 mm， 土粒相對(duì)密度為2.72. 由于試樣尺寸與顆粒粒徑相差較大， 直接生成顆粒工作量大、 計(jì)算時(shí)間過長， 且生成模型不夠直觀， 因此考慮對(duì)顆粒粒徑進(jìn)行適當(dāng)比例放大， 當(dāng)模型最短邊長度為顆粒平均粒徑的30倍以上時(shí)， 可忽略尺寸效應(yīng)[13]. 本文離散元試樣顆粒級(jí)配如圖2所示， 顆粒粒徑在0.6～6.0 mm之間， 顆?？倲?shù)在5 000～8 000之間.

通過伺服機(jī)制， 以迫使模型接近真實(shí)的物理過程， 即通過對(duì)模型邊界條件的調(diào)整， 使得顆粒體系間的接觸盡可能快地達(dá)到理想狀態(tài)， 然后在其基礎(chǔ)上進(jìn)行加載分析， 完成模型試樣的預(yù)固結(jié)過程. 在PFC3D中， 由于墻體在其接觸處是允許發(fā)生一定重疊的剛體， 故減小側(cè)面墻體的剛度來模擬橡皮膜的柔性邊界， 墻體剛度取試樣顆粒剛度的1/10[14]. 在數(shù)值模擬的過程中， 顆粒與墻體之間通過接觸發(fā)生作用， 設(shè)顆粒與墻體間的摩擦系數(shù)為零， 采用線性接觸模型， 顆粒與顆粒間采用Hill接觸模型. 對(duì)于每一個(gè)顆粒， 都能滿足運(yùn)動(dòng)定律， 由于剛性墻體沒有質(zhì)量， 因此墻體受到的不平衡力無法影響墻體的運(yùn)動(dòng)， 即墻體不符合運(yùn)動(dòng)定律， 不能直接施加外力， 只能施加速度， 在PFC三軸剪切試驗(yàn)中， 通過對(duì)墻體施加速度來達(dá)到相應(yīng)的圍壓. 非飽和土三維離散元試樣數(shù)值模型如圖3所示.

圖2 離散元試樣顆粒級(jí)配曲線Fig.2 Particle size distribution of numerical sample

圖3 非飽和土三維離散元試樣數(shù)值模型Fig.3 Numerical model of three-dimensional discrete element specimen of unsaturated soil

剪切階段中， 墻體內(nèi)部上、中、下分別設(shè)置半徑為19.55 mm的測(cè)量圓， 記錄并監(jiān)控不同時(shí)步下模型內(nèi)部應(yīng)力、 應(yīng)變、 孔隙率以及配位數(shù)等變化過程. 此外， 控制側(cè)面墻體的圍壓恒定， 對(duì)固結(jié)后的模型試樣在垂直方向上施加豎向應(yīng)力. 此時(shí)上下墻體按恒定的速率對(duì)模型試樣進(jìn)行加載， 使得試樣在周圍壓力和豎向應(yīng)力共同作用下產(chǎn)生豎向和徑向變形， 從而獲取模型試樣的應(yīng)力-應(yīng)變關(guān)系曲線. 當(dāng)豎向應(yīng)力出現(xiàn)峰值或軸向應(yīng)變達(dá)到15%時(shí)， 則認(rèn)為試樣已經(jīng)發(fā)生破壞[15].

2.3 模型參數(shù)標(biāo)定

圖4給出干密度ρd=1.39 g·cm-3， 飽和度Sr=60%的試樣， 在相同吸力不同圍壓下應(yīng)力-應(yīng)變曲線的離散元數(shù)值結(jié)果， 并與室內(nèi)三軸剪切試驗(yàn)結(jié)果進(jìn)行比較. 從圖4可看出， 數(shù)值模擬結(jié)果與室內(nèi)試驗(yàn)實(shí)測(cè)曲線較為吻合， 可模擬非飽和土的三軸剪切試驗(yàn). 離散元分析中的材料參數(shù)如表2所示.

(a) 吸力與黏聚力的關(guān)系

(b)吸力與內(nèi)摩擦角的關(guān)系

(c) Su=70 kPa， σ3=300 kPa

PFC模型中的細(xì)觀參數(shù)與宏觀參數(shù)之間具有高度非線性的關(guān)系， 且相互之間具有一定耦合性， 即使只調(diào)整一個(gè)細(xì)觀參數(shù)也可能會(huì)同時(shí)影響到多個(gè)宏觀特性. 因此， 細(xì)觀參數(shù)的標(biāo)定順序至關(guān)重要. 對(duì)于Hill接觸模型， 首先， 確定模型所處的圍壓和吸力的大??； 其次標(biāo)定楊氏模量和泊松比， 二者是通過控制顆粒法向和切向剛度的大小， 進(jìn)而影響曲線的斜率； 最后， 對(duì)摩擦系數(shù)進(jìn)行標(biāo)定， 控制曲線峰值高度. 另外， 局部阻尼系數(shù)根據(jù)工程實(shí)踐經(jīng)驗(yàn)所得取常數(shù)值0.7[16]， 其他值與室內(nèi)試驗(yàn)各參數(shù)相對(duì)應(yīng). 此時(shí)三軸剪切試驗(yàn)數(shù)值模型所表現(xiàn)出的宏觀特性， 將與相應(yīng)的室內(nèi)試驗(yàn)結(jié)果接近或基本吻合， 再根據(jù)宏細(xì)觀參數(shù)之間相互影響的一般規(guī)律， 對(duì)細(xì)觀參數(shù)做進(jìn)一步調(diào)整， 即可快速標(biāo)定出數(shù)值模型細(xì)觀參數(shù).

表2 離散元試樣細(xì)觀參數(shù)標(biāo)定結(jié)果

3 非飽和土離散元分析

3.1 應(yīng)力-應(yīng)變關(guān)系曲線

對(duì)于干密度ρd=1.39 g·cm-3， 飽和度Sr=60%的試樣， 在吸力Su=70 kPa和各級(jí)圍壓(100， 200和300 kPa)下， 按恒定的加載速率進(jìn)行豎向壓縮， 通過調(diào)節(jié)模型參數(shù)， 得到同密度、 同吸力， 不同圍壓下的應(yīng)力-應(yīng)變曲線. 在此基礎(chǔ)上， 改變吸力的大小， 根據(jù)應(yīng)力-應(yīng)變曲線的變化規(guī)律， 分別研究圍壓和吸力對(duì)非飽和土抗剪強(qiáng)度的影響.

圖5(a)所示為ρd=1.39 g·cm-3的試樣在吸力Su=70 kPa時(shí)不同圍壓下的應(yīng)力與應(yīng)變關(guān)系， 該試樣呈應(yīng)變軟化特性. 當(dāng)圍壓σ3=100 kPa時(shí)， 最大主應(yīng)力差值在100 kPa附近， 當(dāng)圍壓σ3=200 kPa時(shí)， 最大主應(yīng)力差值在250 kPa附近， 當(dāng)圍壓σ3=300 kPa時(shí)， 最大主應(yīng)力差值在400 kPa附近， 說明圍壓對(duì)試樣的影響較為明顯. 試樣達(dá)到相同軸向應(yīng)變時(shí)， 圍壓σ3=300 kPa的情況下施加的豎向應(yīng)力大于圍壓為200 kPa下施加的豎向應(yīng)力; 同樣， 圍壓σ3=200 kPa下施加的豎向應(yīng)力大于圍壓為100 kPa下的情況， 即增大圍壓有利于提高土樣的強(qiáng)度. 對(duì)比圖5(a)～(d)， 當(dāng)吸力從70 kPa逐步增加到200 kPa， 試樣的應(yīng)力-應(yīng)變曲線由應(yīng)變軟化向應(yīng)變硬化發(fā)展， 說明增大吸力可以改變應(yīng)力-應(yīng)變特性.

(a) Su=70 kPa， ρd=1.39 g·cm-3

(b) Su=100 kPa， ρd=1.39 g·cm-3

(d) Su=200 kPa， ρd=1.39 g·cm-3

3.2 強(qiáng)度特性

根據(jù)密度ρd=1.39 g·cm-3的試樣， 在吸力Su=20、 50、 70、 100、 120和150 kPa下的應(yīng)力-應(yīng)變曲線， 利用莫爾-庫倫強(qiáng)度理論， 獲得不同吸力下試樣的抗剪強(qiáng)度包線， 從而得到非飽和土模型試樣的抗剪強(qiáng)度參數(shù)， 如圖6所示. 同一密度的試樣， 隨著吸力的增加， 其黏聚力和內(nèi)摩擦角均呈非線性增加， 在高吸力段， 黏聚力增加的速度大于低吸力段， 此外， 隨著吸力的增大， 內(nèi)摩擦角增加速度放緩.

(a) 吸力與黏聚力的關(guān)系

(b) 吸力與內(nèi)摩擦角的關(guān)系

3.3 吸力與含水率的關(guān)系

以干密度ρd=1.39 g·cm-3， 飽和度Sr=60%的離散元試樣， 在吸力Su=70 kPa、 圍壓σ3=200 kPa的條件下所得的應(yīng)力-應(yīng)變曲線為基礎(chǔ)， 結(jié)合室內(nèi)試驗(yàn)中非飽和土試樣在相同密度不同飽和度下的應(yīng)力-應(yīng)變曲線， 對(duì)離散元試樣的吸力值進(jìn)行調(diào)整， 標(biāo)定飽和度為40%和20%時(shí)的應(yīng)力-應(yīng)變曲線， 從而獲得不同飽和度下離散元試樣的應(yīng)力-應(yīng)變關(guān)系曲線， 如圖7所示.

由數(shù)值模擬所得的離散元試樣在不同飽和度下的吸力值， 可繪制出Hill接觸模型下吸力與含水率的關(guān)系曲線， 如圖8所示. 在7%～21%的變化范圍內(nèi)， 含水率隨吸力的增大呈線性遞減.

圖7 不同飽和度下離散元試樣的應(yīng)力-應(yīng)變關(guān)系曲線Fig.7 Stress-strain curves of discrete element specimens under different saturation

圖8 吸力與含水率的關(guān)系曲線Fig.8 Relationship between suction and water conten

4 結(jié)語

1) 在剪切過程中其抗剪強(qiáng)度隨圍壓和吸力的增大顯著提高， 且隨著吸力的增大， 應(yīng)力-應(yīng)變曲線由應(yīng)變軟化向應(yīng)變硬化發(fā)展.

2) 黏聚力和內(nèi)摩擦角均隨吸力的增大呈非線性增加， 低吸力段黏聚力增加速度緩于高吸力段， 且隨吸力的增大內(nèi)摩擦角增大速度變緩.

3) 隨著模型試樣含水率的減小， 吸力線性增大， 應(yīng)用Hill接觸模型的顆粒流模型可以較好模擬非飽和土的力學(xué)行為.