基于流固耦合特性的非飽和膨脹土變形仿真計(jì)算

范臻輝，張春順，肖宏彬

(1. 中南大學(xué) 土木建筑學(xué)院，湖南 長(zhǎng)沙，410075；2. 悉尼大學(xué) 巖土工程研究中心，澳大利亞 悉尼，2006；3. 中南林業(yè)科技大學(xué) 土木工程與力學(xué)學(xué)院，湖南 長(zhǎng)沙，410004)

對(duì)于大多數(shù)工程來(lái)說(shuō)，將土當(dāng)作飽和土是一種合理的簡(jiǎn)化，但是，對(duì)于某些特殊區(qū)域或特殊性質(zhì)的土，這種簡(jiǎn)化將造成理論研究結(jié)果與實(shí)際結(jié)果存在誤差。如在膨脹土地基基礎(chǔ)的設(shè)計(jì)中，若單純按照膨脹土的現(xiàn)有強(qiáng)度進(jìn)行設(shè)計(jì)，則有可能將強(qiáng)度參數(shù)估計(jì)過(guò)高，導(dǎo)致工程不安全[1]；若按其最低強(qiáng)度進(jìn)行設(shè)計(jì)，又將造成浪費(fèi)[2]。因此，需合理地提出膨脹土在不同狀態(tài)下的強(qiáng)度參數(shù)，這是工程的客觀需要[3-4]。此外，膨脹土等非飽和土的變形性能也隨飽和度的變化而變化[5-6]。這些問(wèn)題都是飽和土力學(xué)難以解決的。按非飽和狀態(tài)研究膨脹土的工程力學(xué)性質(zhì)是土力學(xué)發(fā)展的趨勢(shì)[7]。在非飽和膨脹土的本構(gòu)模型研究方面，Gens和Alonso提出了一個(gè)非飽和膨脹土彈塑性概念模型(G-A模型)[8]。該模型以Alonso的一般非飽和土彈塑性模型為基礎(chǔ)，從分析膨脹的微觀機(jī)理出發(fā)，反映了膨脹土的濕脹干縮變形特性，但將土體變形分為微觀和宏觀2個(gè)結(jié)構(gòu)層次進(jìn)行計(jì)算，所建模型非常復(fù)雜。許多學(xué)者先后對(duì)G-A模型進(jìn)行了簡(jiǎn)化和改進(jìn)，取得較好的研究成果[9-12]，但這些研究側(cè)重于非飽和膨脹土的應(yīng)力-應(yīng)變本構(gòu)關(guān)系。膨脹土吸水膨脹、失水收縮等問(wèn)題需考慮水、氣兩相流體流動(dòng)和固相變形之間的相互作用，因此，研究非飽和膨脹土流-固耦合問(wèn)題具有重大理論和實(shí)際意義。在此，本文作者按照沈珠江[13]提出的思路，并借用Alonso提出的BBM模型推導(dǎo)過(guò)程建立非飽和膨脹土的彈塑性本構(gòu)模型的矩陣表達(dá)式，并選用適用的土-水特征曲線(xiàn)方程和非飽和土滲透系數(shù)方程，構(gòu)建完整的非飽和膨脹土滲流-變形耦合分析模型，以便更加準(zhǔn)確地分析膨脹土的膨脹收縮變形，進(jìn)一步評(píng)價(jià)膨脹土工程特性。

1 非飽和膨脹土本構(gòu)模型

1.1 應(yīng)力狀態(tài)變量的選擇與有效應(yīng)力原理

Bishop較早提出了非飽和土的有效應(yīng)力公式。假定作用于非飽和土的外應(yīng)力是由非飽和土骨架、孔隙水和孔隙氣共同承擔(dān)。在外應(yīng)力σ作用下，根據(jù)土粒間靜力平衡原理可得Bishop有效應(yīng)力公式：

式中：σ′為土體有效應(yīng)力；ua為孔隙氣壓力；uw為孔隙水壓力；χ為有效應(yīng)力參數(shù)， χ = Aw/A，其值介于0～1.0之間；A為剪切面面積；Aw為孔隙水面積。

從式(1)可以看出：要確定σ′，必須首先知道χ。χ一直是人們對(duì)非飽和土問(wèn)題研究的重點(diǎn)，幾十年來(lái)，人們提出的方法和得到χ的形式多種多樣[14]。在某些情況下，用飽和度Sr描述有效應(yīng)力參數(shù)χ更方便。本文采用單變量理論來(lái)描述非飽和土體的應(yīng)力狀態(tài)，并令)(rSφχ=，則式(1)可改寫(xiě)為：

式中：wauus -= ，為基質(zhì)吸力；Sr為土體飽和度。

1.2 土-水特征曲線(xiàn)

土-水特征曲線(xiàn)又稱(chēng)水分特征曲線(xiàn)，是表示非飽和土的基質(zhì)吸力與含水質(zhì)量分?jǐn)?shù)、含水體積分?jǐn)?shù)、飽和度或有效飽和度之間的關(guān)系曲線(xiàn)。對(duì)于非飽和土，土-水特征曲線(xiàn)的數(shù)學(xué)模型并不是唯一的。土的類(lèi)型不同，所得出的數(shù)學(xué)模型也有所不同，且大部分用于描述土-水特征曲線(xiàn)的數(shù)學(xué)模型都是根據(jù)經(jīng)驗(yàn)、土體結(jié)構(gòu)特征和曲線(xiàn)的形狀而建立的。戚國(guó)慶等[15]采取基質(zhì)吸力的冪函數(shù)多項(xiàng)式形式，建立了非飽和土-水特征曲線(xiàn)的通用數(shù)學(xué)模型：

應(yīng)用式(3)對(duì)非飽和膨脹土體基質(zhì)吸力與含水量或飽和度的關(guān)系進(jìn)行擬合得到相關(guān)系數(shù)，簡(jiǎn)單易行，一般由三次項(xiàng)就能達(dá)到很高精度。

1.3 非飽和土的滲透系數(shù)

對(duì)于飽和土，不論在穩(wěn)定或非穩(wěn)定滲流條件下，都認(rèn)為土的滲透系數(shù)是保持不變的常量。而非飽和土的情況較復(fù)雜，非飽和土的滲透系數(shù)除與土體的種類(lèi)、土體孔隙狀態(tài)、流體的性質(zhì)有關(guān)外，還與土體飽和程度或者土的基質(zhì)吸力密切相關(guān)。非飽和土體的滲透系數(shù)往往小于飽和土體的滲透系數(shù)，滲透系數(shù)不是常數(shù)，是隨飽和度和水的體積分?jǐn)?shù)的變化而變化的。飽和度和水的體積分?jǐn)?shù)的變化往往引起基質(zhì)吸力的變化，因此，非飽和土的滲透系數(shù)常常表征為飽和度、水的體積分?jǐn)?shù)或基質(zhì)吸力的函數(shù)。

Gardner提出非飽和土滲透系數(shù)的計(jì)算公式[16]：

式中：a和n為試驗(yàn)參數(shù)；kw為非飽和土的滲透系數(shù)；ks為這種土在飽和時(shí)的滲透系數(shù)；wρ為水的密度。式(4)物理概念比較明確，是比較常用的表達(dá)式。

1.4 非飽和膨脹土本構(gòu)模型的建立

沈珠江[13]認(rèn)為當(dāng)有效應(yīng)力原理適用時(shí)，任何一個(gè)飽和土本構(gòu)模型都可方便地推廣到非飽和土上，只需把飽和土的有效應(yīng)力公式σ′=σ-uw換成即可。其中：

按照沈珠江的思路，將修正劍橋模型中的有效應(yīng)力進(jìn)行代換，采用相關(guān)聯(lián)流動(dòng)法則，可以得到：

式中：f為屈服條件；g為塑性勢(shì)函數(shù)；p′為有效平均主應(yīng)力，；M為其飽和土的極限線(xiàn)斜率；p0為土體屈服強(qiáng)度，隨基質(zhì)吸力變化而變化。

式中： pc為參考應(yīng)力； p*為飽和土的先期固結(jié)壓力；κ為回彈指數(shù)；λ(0)為飽和情況下對(duì)于應(yīng)力的壓縮指數(shù)；λ為非飽和土相對(duì)于壓應(yīng)力p的壓縮系數(shù)，隨基質(zhì)吸力的變化而變化。

式中：r為常數(shù)， r = λ (s →∞)/λ(0)；β為控制λ隨s增長(zhǎng)速率的參數(shù)。

1.5 非飽和膨脹土彈塑性本構(gòu)模型的矩陣表達(dá)式

在一般情況下，假定非飽和膨脹土中排氣通道暢通。由于排氣速度遠(yuǎn)大于排水速度，因此，可以忽略排氣過(guò)程，即假定各處氣壓力均等于大氣壓。當(dāng)孔隙氣壓力保持不變時(shí)，可采用負(fù)孔隙水壓力代替基質(zhì)吸力[17]，將式(2)中的基質(zhì)吸力采用負(fù)孔隙水壓力代替，并采用矩陣形式可表示為：

由式(6)，(7)和(9)可知：屈服函數(shù)和塑性勢(shì)函數(shù)是有效應(yīng)力σ、孔隙水壓力wu (基質(zhì)吸力)、硬化參數(shù)*p的函數(shù)。

根據(jù)流動(dòng)法則，各應(yīng)力狀態(tài)點(diǎn)的塑性應(yīng)變?cè)隽糠较蚺c通過(guò)該點(diǎn)的塑性勢(shì)面相垂直，且塑性應(yīng)變?yōu)椋?/p>

式中：dλ為非負(fù)的標(biāo)量乘子；dεp為塑性應(yīng)變?cè)隽烤仃?；為塑性體應(yīng)變?cè)隽烤仃嚒?/p>

由應(yīng)力應(yīng)變的增量關(guān)系可知：

式中：eD為彈性剛度矩陣；edε為彈性應(yīng)變?cè)隽浚沪興為總應(yīng)變?cè)隽俊?/p>

對(duì)屈服面上的應(yīng)力狀態(tài)，根據(jù)式(12)，有：

組合式(9)，(14)，(15)和(17)得：

注意到epW 為列向量，epR 為行向量，將孔隙水壓力(或吸力)增量作為1個(gè)附加的應(yīng)變?cè)隽?，可得?/p>

則本構(gòu)方程式(18)和(19)可以轉(zhuǎn)換為：

2 非飽和膨脹土流-固耦合分析

2.1 基本假設(shè)

分析非飽和膨脹土流-固耦合有限元時(shí)，采用如下假設(shè)：(1) 土體為均質(zhì)各向同性材料；(2) 土骨架只有微小應(yīng)變；(3) 土顆粒和孔隙水不可壓縮；(4) 土顆粒和孔隙水變形不受溫度的影響；(5) 孔隙水和氣各自連通，孔隙水運(yùn)動(dòng)服從Darcy定律，而孔隙氣壓力假定為0 Pa。

2.2 控制方程

2.2.1 平衡方程

土體三維的平衡方程可以寫(xiě)成如下的緊湊形式：

應(yīng)用Green-Gauss定理和Galerkin加權(quán)殘值法，由式(22)可得：

式中： Ve為單元體積； Se為單元表面積；，為外表面張力向量； B 和 N 分別uu為應(yīng)變-位移形函數(shù)和位移形函數(shù)；U為節(jié)點(diǎn)位移向量；變量上帶一點(diǎn)“·”表示該變量對(duì)時(shí)間的導(dǎo)數(shù)。由式(11)和(18)可得：

其中：wU 為節(jié)點(diǎn)孔隙水壓力(或基質(zhì)吸力)向量；wN為孔隙水壓力(或基質(zhì)吸力)的形函數(shù)。

將單元矩陣組集起來(lái)形成總體矩陣：

2.2.2 連續(xù)方程

根據(jù)孔隙流體質(zhì)量平衡，可以導(dǎo)出流體的連續(xù)方程：

式中：wρ為孔隙流體的密度；v孔隙流體的流速向量。因?yàn)?/p>

故可以將式(25)轉(zhuǎn)換為：

對(duì)式(26)，在空間離散 uw，并應(yīng)用 Green-Gauss定理和Galerkin加權(quán)殘值法，與平衡方程類(lèi)似，組集成整體矩陣為：

由式(24)和(27)構(gòu)成了非飽和膨脹土流-固耦合分析的控制方程：見(jiàn)圖2。

圖1 計(jì)算模型網(wǎng)格劃分Fig.1 Mesh of model

2.3 流-固耦合有限元數(shù)值解法

對(duì)式(28a)和(28b)，采用交替解法求解。其基本思路是：首先將問(wèn)題在時(shí)域上進(jìn)行離散；求解其中一組方程組(如平衡方程)，將求得的結(jié)果代入另一組方程組中(如連續(xù)方程)求解，再將所得結(jié)果回代到前一組方程組中；若2次求得的值相差較大，不能滿(mǎn)足誤差要求，則重復(fù)以上步驟，直到滿(mǎn)足要求為止；然后，按上述思路繼續(xù)求解下一時(shí)步。

3 算例分析

據(jù)文獻(xiàn)[18]進(jìn)行有荷膨脹試驗(yàn)。該土樣浸水前高為2.00 cm，直徑為6.18cm，初始含水量w=20.77%，初始干密度ρd=1.76 g/cm3，垂直壓力q為25 kPa。采用軸對(duì)稱(chēng)模型，計(jì)算參數(shù)如下：A0=1.019 42，A1=0.067 996 MPa-1，A2=-13.471 64 MPa-2，A3= -23.375 64 MPa-3，a=10-5，n=3.0，ks=10-9m·s-1，M=0.8，λ(0)=0.25，κ=0.05，pc=1 MPa，p*=0.1 MPa，r=0.75，β=12 MPa-1，E=5 MPa，μ=0.3。計(jì)算模型網(wǎng)格劃分見(jiàn)圖1。

整個(gè)計(jì)算分為3個(gè)階段：第1階段，計(jì)算初始應(yīng)力分布；第2階段，計(jì)算垂直荷載；第3階段，計(jì)算浸水膨脹率。其中，第3階段浸水膨脹率的計(jì)算結(jié)果見(jiàn)圖2。

圖2 膨脹率隨時(shí)間的變化曲線(xiàn)Fig.2 Relationship between percent swell and time

從圖2可以看出：采用本文程序所得的計(jì)算值與實(shí)測(cè)值較吻合，尤其在浸水變形的前期階段，吻合程度更大?？梢?jiàn)：采用本文所編制的程序能較好地模擬膨脹土膨脹變形。

圖3所示為各時(shí)刻膨脹土樣的飽和度等值線(xiàn)圖，圖4所示為各時(shí)刻膨脹土樣的孔隙水壓力等值線(xiàn)圖。

從圖3和圖4可以看出：由于土樣是從底部浸水，土樣底部的飽和度比上部的大，同樣，土樣底部的孔隙水壓力比上部的大，隨著時(shí)間的增加，水慢慢滲入，土樣飽和度逐漸增加，孔隙水壓力也逐漸增大，由負(fù)孔隙水壓力(非飽和狀態(tài))變?yōu)檎紫端畨毫?飽和狀態(tài))；到第3 741 min時(shí)，土樣已處于飽和狀態(tài)，土體的孔隙水壓力也全部為正，這時(shí)不會(huì)再有水分進(jìn)入土體，土體處于穩(wěn)定狀態(tài)，其膨脹變形也趨于穩(wěn)定，這一點(diǎn)從圖2也可以得到驗(yàn)證。

圖3 各時(shí)刻點(diǎn)飽和度等值線(xiàn)Fig.3 Contour map of saturation degree axial

圖4 各時(shí)刻點(diǎn)孔隙水壓力等值線(xiàn)Fig.4 Contour maps of pore water pressure axial

4 結(jié)論

(1) 選擇基于有效應(yīng)力原理的單變量理論作為分析依據(jù)，推導(dǎo)出非飽和膨脹土的彈塑性本構(gòu)模型的表達(dá)式，并建立了矩陣表達(dá)式。在該矩陣表達(dá)式中，孔隙水壓力(或吸力)增量作為 1個(gè)附加的應(yīng)變?cè)隽?，可以通過(guò)設(shè)置邊界孔隙水壓(或基質(zhì)吸力)的變化來(lái)模擬膨脹土體的吸水和失水。

(2) 在選用適用的土-水特征曲線(xiàn)方程和非飽和土滲透系數(shù)方程的基礎(chǔ)上，根據(jù)流-固耦合力學(xué)的理論與方法，將膨脹土吸水膨脹或失水收縮過(guò)程視為一個(gè)動(dòng)態(tài)耦合作用過(guò)程，建立了滲流-變形耦合模型。該模型可以計(jì)算飽和和非飽和土體的滲流與固結(jié)，考慮了流-固耦合作用引起的孔隙比和滲透系數(shù)變化(包括有效應(yīng)力的影響)。室內(nèi)有荷膨脹試驗(yàn)?zāi)M結(jié)果驗(yàn)證了模型的可靠性和適用性。

[1] 董隴軍, 李夕兵, 宮鳳強(qiáng). 膨脹土脹縮等級(jí)分類(lèi)的未確知均值聚類(lèi)方法及應(yīng)用[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2008,39(5): 1075-1080.

DONG Long-jun, LI Xi-bing, GONG Feng-qiang. Unascertained average clustering method for classification of grade of shrink and expansion for expansive soils and its application[J]. Journal of Central South University: Science and Technology, 2008,39(5): 1075-1080.

[2] FAN Zhen-hui, WANG Yong-he, XIAO Hong-bin. Analytical method of load-transfer of single pile under expansive soil swelling[J]. Journal of Central South University of Technology,2007, 14(4): 575-579.

[3] XIAO Hong-bin, ZHANG Chun-shun, HE Jie, et al. Expansive soil-structure interaction and its sensitive analysis[J]. Journal of Central South University of Technology, 2007, 14(3): 425-430.

[4] 詹良通, 吳宏偉. 吸力對(duì)非飽和膨脹土抗剪強(qiáng)度及剪脹特性的影響[J]. 巖土工程學(xué)報(bào), 2007, 29(1): 82-87.

[5] ZHAN Liang-tong, WU Hong-wei. Effect of suction on shear strength and dilatancy of an unsaturated expansive clay[J].Chinese Journal of Geotechnical Engineering, 2007, 29(1):82-87.

[6] 韓華強(qiáng), 陳生水, 鄭澄鋒. 非飽和膨脹土強(qiáng)度及變形特性試驗(yàn)研究[J]. 巖土工程學(xué)報(bào), 2008, 30(12): 1872-1876.

HAN Hua-qiang, CHEN Sheng-shui, ZHENG Cheng-feng.Experimental study on strength and deformation of unsaturated expansive soils[J]. Chinese Journal of Geotechnical Engineering,2008, 30(12): 1872-1876.

[7] 范臻輝, 王永和, 肖宏彬. 南寧膨脹土膨脹變形規(guī)律的試驗(yàn)與應(yīng)用研究[J]. 鐵道學(xué)報(bào), 2007, 29(1): 86-89.

FAN Zhen-hui, WANG Yong-he, XIAO Hong-bin. Experimental and applied research on swelling laws of nanning expansive soil[J]. Journal of the China Railway Society, 2007, 29(1):86-89.

[8] 殷宗澤. 土力學(xué)學(xué)科發(fā)展的現(xiàn)狀與展望[J]. 河海大學(xué)學(xué)報(bào),1999, 27(1): 1-5.

YIN Zong-ze. Status and prospects of soil mechanics development[J]. Journal of Hohai University, 1999, 27(1): 1-5.

[9] Gens A, Alonso E E. A framework for the behavior of unsaturated expansive clays[J]. Ganadian Geotechnique Journal,1992, 29: 1013-1032.

[10] 盧再華, 王權(quán)民, 陳正漢. 非飽和膨脹土本構(gòu)模型的試驗(yàn)研究及分析[J]. 地下空間, 2001, 21(5): 379-385.

LU Zai-hua, WANG quan-min, CHEN Zheng-han. Research on the constitutive model of unsaturated expansive soil[J].Underground Space, 2001, 21(5): 379-385.

[11] 盧再華, 陳正漢. 非飽和原狀膨脹土的彈塑性損傷本構(gòu)模型研究[J]. 巖土工程學(xué)報(bào), 2003, 25(4): 422-426.

LU Zai-hua, CHEN Zheng-han. An elastoplastic damage constitutive model of unsaturated undisturbed expansive soil[J].Chinese Journal of Geotechnical Engineering, 2003, 25(4):422-426.

[12] 姚志華, 陳正漢, 黃雪峰, 等. 結(jié)構(gòu)損傷對(duì)膨脹土屈服特性的影響[J]. 巖石力學(xué)與工程學(xué)報(bào), 2010, 29(7): 1503-1512.

YAO Zhi-hua, CHEN Zheng-han, HUANG Xue-feng, et al.Influence of structural damage on yielding characteristics of expansive soils[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(7): 1503-1512.

[13] 沈珠江. 非飽和土力學(xué)實(shí)用化之路探索[J]. 巖土工程學(xué)報(bào),2006, 28(2): 256-259.

SHEN Zhu-jiang. Exploitation of practical use of unsaturated soil mechanics[J]. Chinese Journal of Geotechnical Engineering,2006, 28(2): 256-259.

[14] 邢義川, 謝定義, 李振. 非飽和土的有效應(yīng)力參數(shù)研究[J]. 水利學(xué)報(bào), 2000, 31(12): 77-81.

XING Yi-chuan, XIE ding-yi, LI Zhen. Study on effective stress parameter of unsaturated soi[J].Journal of Hydraulic Engineering,2000, 31(12): 77-81.

[15] 戚國(guó)慶, 黃潤(rùn)秋. 土水特征曲線(xiàn)的通用數(shù)學(xué)模型研究[J]. 工程地質(zhì)學(xué)報(bào), 2004, 12(2): 182-186.

QI Guo-qing, HUANG Run-qiu. An universal mathematical model of soil-water characteristic curve[J]. Journal of Engineering Geology, 2004, 12(2): 182-186.

[16] Desai C S, Christian J T. Numerical Methods in Geotechnical Engineering[M]. New York: McGraw Hill Book Co, 1977:307-329.

[17] Sheng D, Sloan S W, Gens A, Smith D W. Finite element formulation and algorithms for unsaturated soils, Part I:Theory[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2003, 27: 745-765.

[18] 范臻輝. 膨脹土地基脹縮特性及樁-土相互作用研究[D]. 長(zhǎng)沙: 中南大學(xué)土木建筑學(xué)院, 2007: 26-43.

FAN Zhen-hui. Research on swelling-shrinking characteristic and pile-soil interaction of expansive soil foundation[D].Changsha: Central South University. School of Civil Engineering and Architecture, 2007: 26-43.