基于流固耦合特性的非饱和膨胀土变形仿真计算

范臻辉，张春顺，肖宏彬

(1. 中南大学 土木建筑学院，湖南 长沙，410075；2. 悉尼大学 岩土工程研究中心，澳大利亚 悉尼，2006；3. 中南林业科技大学 土木工程与力学学院，湖南 长沙，410004)

对于大多数工程来说，将土当作饱和土是一种合理的简化，但是，对于某些特殊区域或特殊性质的土，这种简化将造成理论研究结果与实际结果存在误差。如在膨胀土地基基础的设计中，若单纯按照膨胀土的现有强度进行设计，则有可能将强度参数估计过高，导致工程不安全[1]；若按其最低强度进行设计，又将造成浪费[2]。因此，需合理地提出膨胀土在不同状态下的强度参数，这是工程的客观需要[3-4]。此外，膨胀土等非饱和土的变形性能也随饱和度的变化而变化[5-6]。这些问题都是饱和土力学难以解决的。按非饱和状态研究膨胀土的工程力学性质是土力学发展的趋势[7]。在非饱和膨胀土的本构模型研究方面，Gens和Alonso提出了一个非饱和膨胀土弹塑性概念模型(G-A模型)[8]。该模型以Alonso的一般非饱和土弹塑性模型为基础，从分析膨胀的微观机理出发，反映了膨胀土的湿胀干缩变形特性，但将土体变形分为微观和宏观2个结构层次进行计算，所建模型非常复杂。许多学者先后对G-A模型进行了简化和改进，取得较好的研究成果[9-12]，但这些研究侧重于非饱和膨胀土的应力-应变本构关系。膨胀土吸水膨胀、失水收缩等问题需考虑水、气两相流体流动和固相变形之间的相互作用，因此，研究非饱和膨胀土流-固耦合问题具有重大理论和实际意义。在此，本文作者按照沈珠江[13]提出的思路，并借用Alonso提出的BBM模型推导过程建立非饱和膨胀土的弹塑性本构模型的矩阵表达式，并选用适用的土-水特征曲线方程和非饱和土渗透系数方程，构建完整的非饱和膨胀土渗流-变形耦合分析模型，以便更加准确地分析膨胀土的膨胀收缩变形，进一步评价膨胀土工程特性。

1 非饱和膨胀土本构模型

1.1 应力状态变量的选择与有效应力原理

Bishop较早提出了非饱和土的有效应力公式。假定作用于非饱和土的外应力是由非饱和土骨架、孔隙水和孔隙气共同承担。在外应力σ作用下，根据土粒间静力平衡原理可得Bishop有效应力公式：

式中：σ′为土体有效应力；ua为孔隙气压力；uw为孔隙水压力；χ为有效应力参数， χ = Aw/A，其值介于0～1.0之间；A为剪切面面积；Aw为孔隙水面积。

从式(1)可以看出：要确定σ′，必须首先知道χ。χ一直是人们对非饱和土问题研究的重点，几十年来，人们提出的方法和得到χ的形式多种多样[14]。在某些情况下，用饱和度Sr描述有效应力参数χ更方便。本文采用单变量理论来描述非饱和土体的应力状态，并令)(rSφχ=，则式(1)可改写为：

式中：wauus -= ，为基质吸力；Sr为土体饱和度。

1.2 土-水特征曲线

土-水特征曲线又称水分特征曲线，是表示非饱和土的基质吸力与含水质量分数、含水体积分数、饱和度或有效饱和度之间的关系曲线。对于非饱和土，土-水特征曲线的数学模型并不是唯一的。土的类型不同，所得出的数学模型也有所不同，且大部分用于描述土-水特征曲线的数学模型都是根据经验、土体结构特征和曲线的形状而建立的。戚国庆等[15]采取基质吸力的幂函数多项式形式，建立了非饱和土-水特征曲线的通用数学模型：

应用式(3)对非饱和膨胀土体基质吸力与含水量或饱和度的关系进行拟合得到相关系数，简单易行，一般由三次项就能达到很高精度。

1.3 非饱和土的渗透系数

对于饱和土，不论在稳定或非稳定渗流条件下，都认为土的渗透系数是保持不变的常量。而非饱和土的情况较复杂，非饱和土的渗透系数除与土体的种类、土体孔隙状态、流体的性质有关外，还与土体饱和程度或者土的基质吸力密切相关。非饱和土体的渗透系数往往小于饱和土体的渗透系数，渗透系数不是常数，是随饱和度和水的体积分数的变化而变化的。饱和度和水的体积分数的变化往往引起基质吸力的变化，因此，非饱和土的渗透系数常常表征为饱和度、水的体积分数或基质吸力的函数。

Gardner提出非饱和土渗透系数的计算公式[16]：

式中：a和n为试验参数；kw为非饱和土的渗透系数；ks为这种土在饱和时的渗透系数；wρ为水的密度。式(4)物理概念比较明确，是比较常用的表达式。

1.4 非饱和膨胀土本构模型的建立

沈珠江[13]认为当有效应力原理适用时，任何一个饱和土本构模型都可方便地推广到非饱和土上，只需把饱和土的有效应力公式σ′=σ-uw换成即可。其中：

按照沈珠江的思路，将修正剑桥模型中的有效应力进行代换，采用相关联流动法则，可以得到：

式中：f为屈服条件；g为塑性势函数；p′为有效平均主应力，；M为其饱和土的极限线斜率；p0为土体屈服强度，随基质吸力变化而变化。

式中： pc为参考应力； p*为饱和土的先期固结压力；κ为回弹指数；λ(0)为饱和情况下对于应力的压缩指数；λ为非饱和土相对于压应力p的压缩系数，随基质吸力的变化而变化。

式中：r为常数， r = λ (s →∞)/λ(0)；β为控制λ随s增长速率的参数。

1.5 非饱和膨胀土弹塑性本构模型的矩阵表达式

在一般情况下，假定非饱和膨胀土中排气通道畅通。由于排气速度远大于排水速度，因此，可以忽略排气过程，即假定各处气压力均等于大气压。当孔隙气压力保持不变时，可采用负孔隙水压力代替基质吸力[17]，将式(2)中的基质吸力采用负孔隙水压力代替，并采用矩阵形式可表示为：

由式(6)，(7)和(9)可知：屈服函数和塑性势函数是有效应力σ、孔隙水压力wu (基质吸力)、硬化参数*p的函数。

根据流动法则，各应力状态点的塑性应变增量方向与通过该点的塑性势面相垂直，且塑性应变为：

式中：dλ为非负的标量乘子；dεp为塑性应变增量矩阵；为塑性体应变增量矩阵。

由应力应变的增量关系可知：

式中：eD为弹性刚度矩阵；edε为弹性应变增量；εd为总应变增量。

对屈服面上的应力状态，根据式(12)，有：

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

注意到epW 为列向量，epR 为行向量，将孔隙水压力(或吸力)增量作为1个附加的应变增量，可得：

则本构方程式(18)和(19)可以转换为：

2 非饱和膨胀土流-固耦合分析

2.1 基本假设

分析非饱和膨胀土流-固耦合有限元时，采用如下假设：(1) 土体为均质各向同性材料；(2) 土骨架只有微小应变；(3) 土颗粒和孔隙水不可压缩；(4) 土颗粒和孔隙水变形不受温度的影响；(5) 孔隙水和气各自连通，孔隙水运动服从Darcy定律，而孔隙气压力假定为0 Pa。

2.2 控制方程

2.2.1 平衡方程

土体三维的平衡方程可以写成如下的紧凑形式：

应用Green-Gauss定理和Galerkin加权残值法，由式(22)可得：

式中： Ve为单元体积； Se为单元表面积；，为外表面张力向量； B 和 N 分别uu为应变-位移形函数和位移形函数；U为节点位移向量；变量上带一点“·”表示该变量对时间的导数。由式(11)和(18)可得：

其中：wU 为节点孔隙水压力(或基质吸力)向量；wN为孔隙水压力(或基质吸力)的形函数。

将单元矩阵组集起来形成总体矩阵：

2.2.2 连续方程

根据孔隙流体质量平衡，可以导出流体的连续方程：

式中：wρ为孔隙流体的密度；v孔隙流体的流速向量。因为

故可以将式(25)转换为：

对式(26)，在空间离散 uw，并应用 Green-Gauss定理和Galerkin加权残值法，与平衡方程类似，组集成整体矩阵为：

由式(24)和(27)构成了非饱和膨胀土流-固耦合分析的控制方程：见图2。

图1 计算模型网格划分Fig.1 Mesh of model

2.3 流-固耦合有限元数值解法

对式(28a)和(28b)，采用交替解法求解。其基本思路是：首先将问题在时域上进行离散；求解其中一组方程组(如平衡方程)，将求得的结果代入另一组方程组中(如连续方程)求解，再将所得结果回代到前一组方程组中；若2次求得的值相差较大，不能满足误差要求，则重复以上步骤，直到满足要求为止；然后，按上述思路继续求解下一时步。

3 算例分析

据文献[18]进行有荷膨胀试验。该土样浸水前高为2.00 cm，直径为6.18cm，初始含水量w=20.77%，初始干密度ρd=1.76 g/cm3，垂直压力q为25 kPa。采用轴对称模型，计算参数如下：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。计算模型网格划分见图1。

整个计算分为3个阶段：第1阶段，计算初始应力分布；第2阶段，计算垂直荷载；第3阶段，计算浸水膨胀率。其中，第3阶段浸水膨胀率的计算结果见图2。

图2 膨胀率随时间的变化曲线Fig.2 Relationship between percent swell and time

从图2可以看出：采用本文程序所得的计算值与实测值较吻合，尤其在浸水变形的前期阶段，吻合程度更大。可见：采用本文所编制的程序能较好地模拟膨胀土膨胀变形。

图3所示为各时刻膨胀土样的饱和度等值线图，图4所示为各时刻膨胀土样的孔隙水压力等值线图。

从图3和图4可以看出：由于土样是从底部浸水，土样底部的饱和度比上部的大，同样，土样底部的孔隙水压力比上部的大，随着时间的增加，水慢慢渗入，土样饱和度逐渐增加，孔隙水压力也逐渐增大，由负孔隙水压力(非饱和状态)变为正孔隙水压力(饱和状态)；到第3 741 min时，土样已处于饱和状态，土体的孔隙水压力也全部为正，这时不会再有水分进入土体，土体处于稳定状态，其膨胀变形也趋于稳定，这一点从图2也可以得到验证。

图3 各时刻点饱和度等值线Fig.3 Contour map of saturation degree axial

图4 各时刻点孔隙水压力等值线Fig.4 Contour maps of pore water pressure axial

4 结论

(1) 选择基于有效应力原理的单变量理论作为分析依据，推导出非饱和膨胀土的弹塑性本构模型的表达式，并建立了矩阵表达式。在该矩阵表达式中，孔隙水压力(或吸力)增量作为 1个附加的应变增量，可以通过设置边界孔隙水压(或基质吸力)的变化来模拟膨胀土体的吸水和失水。

(2) 在选用适用的土-水特征曲线方程和非饱和土渗透系数方程的基础上，根据流-固耦合力学的理论与方法，将膨胀土吸水膨胀或失水收缩过程视为一个动态耦合作用过程，建立了渗流-变形耦合模型。该模型可以计算饱和和非饱和土体的渗流与固结，考虑了流-固耦合作用引起的孔隙比和渗透系数变化(包括有效应力的影响)。室内有荷膨胀试验模拟结果验证了模型的可靠性和适用性。

[1] 董陇军, 李夕兵, 宫凤强. 膨胀土胀缩等级分类的未确知均值聚类方法及应用[J]. 中南大学学报: 自然科学版, 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] 詹良通, 吴宏伟. 吸力对非饱和膨胀土抗剪强度及剪胀特性的影响[J]. 岩土工程学报, 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] 韩华强, 陈生水, 郑澄锋. 非饱和膨胀土强度及变形特性试验研究[J]. 岩土工程学报, 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] 范臻辉, 王永和, 肖宏彬. 南宁膨胀土膨胀变形规律的试验与应用研究[J]. 铁道学报, 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] 殷宗泽. 土力学学科发展的现状与展望[J]. 河海大学学报,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] 卢再华, 王权民, 陈正汉. 非饱和膨胀土本构模型的试验研究及分析[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] 卢再华, 陈正汉. 非饱和原状膨胀土的弹塑性损伤本构模型研究[J]. 岩土工程学报, 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] 姚志华, 陈正汉, 黄雪峰, 等. 结构损伤对膨胀土屈服特性的影响[J]. 岩石力学与工程学报, 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] 沈珠江. 非饱和土力学实用化之路探索[J]. 岩土工程学报,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] 邢义川, 谢定义, 李振. 非饱和土的有效应力参数研究[J]. 水利学报, 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] 戚国庆, 黄润秋. 土水特征曲线的通用数学模型研究[J]. 工程地质学报, 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]. 长沙: 中南大学土木建筑学院, 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.