土体流固耦合理论研究进展

陆培毅,韩亚飞,王成华

(天津大学 建筑工程学院，天津 300350)

0 引 言

土体流固耦合是研究土体渗流场和应力场之间的相互作用，涉及到土力学，流体力学和两学科的交叉。太沙基最早提出有效应力原理并基于该原理提出了一维固结理论。此后，比奥基于弹性理论，平衡条件，连续条件等提出了真三维固结理论[1]，以上两种理论是土力学中最具有代表性的流固耦合理论。后来，诸多学者分别从耦合形式，有限元等方面对以上两种理论进行了改进和应用，形成了现在的流固耦合理论。假设条件和耦合角度的不同，得出耦合理论的适用范围和精度也不尽相同。因此有必要弄清各种耦合理论的前提条件、内在联系、区别和优缺点。所以笔者对现有的耦合理论从计算域耦合、基本未知数耦合、参数耦合和本构关系耦合4个方面进行分类，目的在深入了解每种理论并总结优缺点，以便在理论指导实际工程时做到扬长避短。

1 计算域耦合

1.1 概 述

计算域之间的耦合是指固体域和流体域之间的相互作用，按作用机理可以分为两种，一种是固体域和流体域部分或完全重叠在一起。该种问题的求解需要分别建立固体变形和流体压力的方程，然后根据两者之间的变形协调或者两者之间的相互作用力建立耦合关系。Terzaghi一维固结理论就是基于这种方法建立的。另一种是固体域和流体域完全分开，既耦合作用仅发生在两相交界面处。这种问题的求解需要建立起两相耦合面上力的平衡和位移协调的关系。问延煦等[2]提出的附加质量法以及水体位移元法就是基于这种方法建立的。

1.2 理论发展

太沙基一维固结理论是在一系列假设上建立起来的，伦杜利克把太沙基的一维固结理论推广到二维和三维，后来的研究者对太沙基一维固结理论进行改进、推广，使其更符合工程实际。改进一般基于土体的各向异性、非达西渗流和土的弹塑性本构等。问延煦等[2]已发表了太沙基一维固结理论的研究综述，且内容详实，此处不再赘述。

一维固结理论的不足已经被人们熟知，例如假设土体弹性，渗透系数不变等，但由于其原理清晰，计算简便，一直被应用至今。要深刻理解其应用假设和前提，以免造成工程事故。

罗晓辉[3]对渗流场进行了稳定渗流与非稳定渗流有限元分析，将渗流力的作用等效为节点荷载施加到划分好的土体上进行分析，然后将应力状态的变化带入渗透系数与应力状态的关系式，从而实现应力场和渗流场的耦合。

使用有限元方法，将基坑降水引起的渗流力等效成节点荷载施加到土体单元网格上来模拟渗透力对土体有效应力的影响，反之，将土中有效应力的变化与渗透系数建立起联系，实现了流固耦合分析，并以此来分析大面积降水引起的基坑稳定问题。

将渗流力简化为体力施加到土体有限元节点方法的物理意义较为明确且符合实际，将其与土体应变-渗透系数之间的关系联合应用就可以完整的模拟流固耦合现象，有较高的应用价值。较为困难的是根据不同的土质建立起有效应力和渗透系数之间的关系，有待进一步研究。

H. M. WESTERGAARD[4]提出了动水压力附加质量法公式，他将坝面视为直立且刚性的，大大化简了水与坝体之间的流固耦合作用。该理论存在假定，只是坝体和水的一种理想化模型，但反应了动水压力的一些本质特征且有计算简便的优点，适合对精度要求不高的小型工程进行初步估算。

(1)

式中：Pw(h)为水深h处的动水压力;ah为水平设计地震加速度代表值;ρ为水体质量密度;H0为库水深。

CLOUGH推广了附加质量法公式，使其适用于任意形状的坝体，并提出了适用于有限元分析的广义质量表达式：

(2)

式中：Mai为附加质量；λi是坝面点的法向矢量；Ai为该点在坝面上的隶属面积。

P. CHAKRABARTI等[5-6]分别在1972年和1973年对附加质量模型做出了修正，改进后的模型可以反应坝面的弹性，这样就更接近坝水流固耦合实际。

朱骁健等[7]在动水附加质量法公式的基础上，考虑库水的辐射阻尼耗能，采用库水阻尼元件法将库水惯性质量和阻尼元件共同作用在大坝迎水面上。阻尼元件的阻尼系数根据能量在边界的反射与透射得到，由库水密度、水中纵波速度、单个阻尼元件的影响面积和边界输入、输出能量的比值决定。将结果与实际对比，考虑水的可压缩性更符合工程实际。

动水压力法是建立在诸多假设上得出的，例如采用刚性坝面假设，忽略了坝体的弹性变形和库水的状态对动力响应的影响，这就导致其不是真正的流固耦合，且一般认为动水附加质量法夸大了动水压力作用，与实际情况存在较大差距。

E.L.WILSON等[8]首次使用水体位移元法来求解水和坝体之间的流固耦合问题，与动水压力法相比，该种方法将库水和坝体都以结点位移为主变量，使其统一在一个系统中求解坝体与库水的动力耦合问题。

陈和群等[9]通过限制水体的旋转变量实现无旋性，用水体位移元法，建立坝体与库水动力耦合的有限元法分析模型。通过分析二滩双曲拱坝与库水的耦合作用后发现水体位移元法能较好地模拟坝与库水动力实际的相互作用。

水体位移元法同时将水的节点位移作为主变量，因此考虑了水体的可压缩性，但目前研究较少,可以进行深入研究。

在计算域之间耦合的各种计算方法大致可以分为满足变形协调条件的太沙基一维固结理论、将渗流力简化为体力施加到土体有限元节点的方法、动水压力法和水体位移元法。其中，渗流力简化为节点力的方法可以结合下文中的基本未知数耦合理论进一步完善。水体位移法也应进一步研究。

2 基本未知数耦合

2.1 概 述

基本未知数之间的耦合指事先假定若干个未知量，然后根据耦合现象必须满足的物理力学条件，边界条件等列出方程，求解出未知量，从而达到耦合的目的。此种方法以比奥固结理论为代表，后人在其理论基础上进行完善。SAVAGE和BRADDOCK将三维比奥固结理论应用于对各向同性的孔隙弹性介质的分析。ZIENKIEWIEZ和SHIOMI在比奥固结理论的基础上，考虑了多孔介质的几何、边界及材料的非线性。总之，对比奥固结理论的完善主要是从有限元、本构模型、非达西渗流三方面开展的。

2.2 理论发展

比奥固结理论[1]先假定4个未知量(3个方向的位移和超静孔隙水压力)，然后通过弹性理论，线弹性应力应变关系，水流连续方程解出未知量，从而得到应力场和渗流场之间的耦合关系。

R. S. SANDHU等[10]最早使用有限元法来求解比奥固结方程。

在国内，沈珠江应用变分原理首先把比奥固结理论的有限单元法应用于固结分析。随后，邓岳保等[11]分别采用不同的方法推导了比奥固结理论有限元方程。

比奥固结理论假定土体为线弹性体，平扬等[12]将土体本构关系中的弹性矩阵替换为弹塑性矩阵，这就将该理论推广到弹塑性分析领域，并通过水土特征曲线求出自由面的渗透系数，改进了比奥固结理论中不考虑非饱和效应的假设。

王媛[13]以比奥理论为基拙 , 假定土体饱和且土骨架满足线弹性各向同性体。使用伽辽金法对比奥固结方程进行了求解,同时考虑了正交各向异性的达西渗流。

对于比奥固结理论的改进，无论是将土体视为各向同性还是各向异性，其分析原理都是研究外荷载下孔隙水压力与有效应力(或总应力)及相应的变形之间的关系。因而，这种分析方法其实与渗透力作用无关，所探讨的渗流场实质是孔隙水压力分布场，并不是渗流作用下所形成的渗透力场。

固结理论的非达西渗流多集中在一维固结理论的研究中，而对于比奥固结的非达西渗流研究较少。太沙基在1925年指出，达西定律不适用于塑性大的黏性土。刘慈群[14]在PASCAL的基础上，得出了上述问题的近似解。刘忠玉等[15]引入可以同时考虑低速渗流曲线段和较高速渗流直线段的非达西渗流方程，重新推导饱和黏土一维固结方程，并采用有限体积法对该方程进行数值求解。并探讨非达西渗流参数对固结过程的影响。计算结果表明，非达西渗流延缓了饱和黏土中孔隙水压力的消散速度，故使得地基的固结速度比太沙基一维固结理论值要慢，最后讨论了Terzaghi一维固结理论的适用范围。李传勋等[16]研究了基于指数形式、非牛顿指数的一维固结理论。

S. HANSBO[17]最早对非达西渗流进行了研究，他首次推导了考虑非达西渗流的竖井地基固结解析解。S. HANSBO的非达西渗流模型表达式：

(3)

式中：i0为起始水力梯度；iL为门槛水力梯度。

T. C. ING等[18]基于虚功原理推导了HANSBO非达西渗流模型的轴对称比奥固结有限元方程，并分析了非达西渗流对固结计算的影响。除了得出非达西渗流对固结的延缓作用这一普遍结论，还得出，只有在iL≥40且m≥1.5时，非达西渗流才对固结有明显的影响。另外，应力历史对固结速率重要影响。

在比奥固结理论基础上的非达西渗流研究中，非达西渗流能更好的描述黏性土和小荷载情况下的渗流，但由于土体的成层性导致的水平向和竖向渗透系数的差异是非达西渗流不能描述的。

3 参数耦合

3.1 概述

流固耦合是一个动态过程，在实际的渗流过程中, 由于孔隙流体压力的变化, 一方面要引起多孔介质骨架有效应力变化, 由此导致多孔介质渗透率和孔隙率的变化;另一方面, 这些变化又反过来影响孔隙流体的流动和压力的分布。参数之间的耦合就是建立起土体应力或应变与渗透系数或孔隙比之间的关系，达到耦合的目的。其中渗透系数则是应力场、渗流场相互耦合的“桥梁”，也是实现真正的渗流耦合分析的关键。对于假定渗透系数为常数的渗流耦合分析，只能体现渗流场对应力场的影响，不能体现应力场对渗流场的影响。

3.2 理论发展

Kozeny-Craman公式建立了土体渗透系数与孔隙率的关系式：

(4)

式中：K0，n0分别为初始渗透系数和孔隙率；K，n分别为耦合分析渗透系数和孔隙率。

骆祖江等[19]在耦合计算分析时，将n=(n0+εv)/(1+εv)带入式(4)得到渗透系数和体应变的关系。

李培超等[21]基于考虑孔隙率和孔隙流体压力的多孔介质有效应力原理，推导出孔隙率和渗透系数的动态变化关系，然后建立应力场和渗流场方程，得到比较完善的流固耦合数学模型。

在渗透系数和孔隙比的关系中，渗透系数是随着孔隙率变化的，因此用统一的渗透系数来表示整个渗流场的渗流特性是不合理的。

渗透系数与孔隙比关系式汇总如表1：

表1 渗透系数与孔隙比关系汇总Table 1 Summary of relationships between void ratios and permeability coefficients

注：其中Deff为土平均有效粒径；ek与当前应力状态相关的空隙比；CF为土体的粘粒含量；Ac为土体活性指数；Ck为土体渗透系数相关的指标，一般Ck=0.5e0。

冉启全等[26]研究了油藏数值模拟中的流固耦合，其中通过Kozeny-Craman方程，推导出等温渗流过程中渗透率与体积应变的关系，建立起油藏开采过程中，油藏压力变化和岩体变形之间的固耦合关系。

(5)

式中：K0，K意义同上；φ0为孔隙度；εv为体积应变。

陈晓平等[27]给出了非均质土坝渗流场和应力场耦合的数学模型，其方法是将渗流场和应力场分开求解，分别写出与应力场有关的渗流场方程和与渗流场有关的应力场方程，并采用渗透系数和应变关系的经验公式：

(6)

李筱艳等[28]根据抽水试验及沉降观测资料，建立了土体渗透系数与有效应力增量的非线性耦合模型，通过线性回归得到了它们的关系式：

K=K0exp(-λ×Δσ)

(7)

式中：λ为试验参数，反应了土体中渗透系数随有效应力变化的幅度。

马少坤等[29]采用ABAQUS中修正剑桥模型，推导出K0固结状态下土体的孔隙比随深度变化的关系，并建立渗透系数随深度线性变化的关系，然后导入修正剑桥模型进行降水开挖流固耦合分析，并与不考虑降水的开挖分析作对比，从而得到考虑孔隙比和渗透系数随深度变化时基坑降水开挖流固耦合作用下的围护结构的变形及土体变形特征。

在渗透系数和体应变的关系中，对不同的应力路径，用统一的关系式表示渗流系数与体应变的关系也是不符合实际的。因为在不同应力路径条件下，即使各主应变不同，体应变也可能会相同，但是不同的主应变对渗透系数的影响不同。例如在竖向加载时，主要表现为ε1的变化，此时对竖向渗透系数影响较大；而在侧向卸载时，主要表现为ε3的变化，此时对水平向渗透系数影响较大，但两种情况下的体应变可能是一样的。

4 本构关系耦合

4.1 概 述

土体液化和管涌现象均涉及流固耦合。土体一般在动力荷载作用下发生液化，研究液化机理必须研究土体中孔隙水压力的变化以及土体和流体之间的相互作用；管涌现象多发生在颗粒级配不良的粉土中，土体在渗流作用下，逐渐被侵蚀转化为液化细颗粒，液化细颗粒随水流失导致孔隙率增大，从而增强了土体的渗透性；渗透性的增强使得渗流速度增大，渗流速度的增大进一步加剧了土体的侵蚀。即侵蚀与渗流之间存在着耦合效应，两者相互促进，相互影响，联系两者的纽带就是孔隙率。在研究该问题时需要提出具体的本构模型来反应耦合关系。

4.2 管涌多相耦合理论的发展

大量试验和现场资料表明，管涌的发生是一个由骨架相、液相和液化细颗粒相之间的相互作用的过程。管涌产生机理如图1，管涌多相耦合示意图如图2[30]。多相对管涌的研究经历了不考虑流固耦合到考虑流固耦合再到考虑侵蚀本构关系耦合3个阶段，使对管涌的产生机理有了更深刻的认识。

周恩全[31]研究了饱和砂土液化后的流体本构模型，及将液化后的土体当作流体来研究，是一种较新颖的求解思路，他通过实验发现液化后的砂土中剪应力和孔压比之间有较好的线性关系。提出以下公式：

(8)

式中：τ为土体剪应力；pα为标准大气压；A和B为试验参数。

式(8)可以反映土体发生液化时，强度和空隙水压力的关系,但是不能反映变形和空隙水压力之间的相互影响。

I. VARDOULAKIS等[32-33]选取孔隙率为耦合参数，建立了包含渗流侵蚀本构方程、混合渗流平衡方程和力平衡方程在内的数学模型，采用伽辽金有限元法模拟了油井出砂问题。J. BELL最早提出了渗流侵蚀本构方程：

(9)

式中：m为任意时刻土骨架相转化为可动细颗粒相的速率；mcr为颗粒质量侵蚀速率；mdep为颗粒沉积速率；c为可动细颗粒浓度；qi为水相和可动颗粒相混合物的流量；φ为土体孔隙率；ρs为土骨架相密度。

管涌多相耦合多采用的有限差分和有限元模型或者不能体现管涌渐进性破坏特点或者不能考虑管涌发展过程中孔隙水压力的变化，即不能考虑管涌发展过程的渗流侵蚀耦合效应。

渗流侵蚀的耦合机理将土体分为土骨架相、液相和液化颗粒相，基于多孔介质动力学理论建立三相质量守恒微分方程，再引入可以考虑土骨架的侵蚀导致液化颗粒相质量增加的渗流侵蚀本构方程，最后得出了一维三相渗流侵蚀耦合管涌数学模型。

L.SIBILLE等[34]使用离散元模拟土颗粒，格子玻尔兹曼方法模拟水流，然后在两种方法之间建立耦合关系，从微观角度研究了管涌现象中土颗粒与土骨架分离和土颗粒随水流迁移两个阶段。分析时，土颗粒之间采用摩擦和黏结接触本构，水流和土颗粒之间的作用，用水力剪力和来反映。

图1 管涌产生机理Fig. 1 The mechanism of piping

图2 管涌的多相耦合机理Fig. 2 The multiphase coupling mechanism of piping

注：其中u为位移；n为孔隙率；k为渗透系数；p为孔隙水压力；q为渗透速度；c为细颗粒浓度；E为弹性模量；v为泊松比；φ为内摩擦角。

渗流侵蚀本构方程描述了液化细颗粒的流失量与渗流速度、管涌土体密度、孔隙率、液化细颗粒含量等物理量的关系，而现有管涌数学模型中常用的渗流侵蚀本构方程来源与石油工程中常用的砂岩的侵蚀方程，它与土木工程中常见的砂土，粉质土等的侵蚀性质不同，将其用于管涌数值模拟会产生误差。

5 结论与展望

5.1 结 论

对土体流固耦合的4种分类方法，揭示每种理论的本质，有助于加深对流固耦合问题的认识，同时也指出了流固耦合存在的问题和发展方向。

1)计算域之间耦合的理论发展较早，计算简单，但这是建立在过度简化基础上得出的，所以这种耦合精度较差。该种方法应进一步研究更复杂的耦合条件，例如不同介质的接触条件。

2)基本未知数之间耦合的理论较为完善，耦合的计算模型与实际符合的较好，缺点是计算繁琐，结合有限单元法可以在工程中得到充分的利用。

3)参数之间耦合能直观的反应耦合现象，但是参数之间的耦合多是根据某一种土的实验数据推导出来的，没有普遍性，研究有待加深。

4)本构关系耦合是研究管涌、土体液化的理论基础，但由于该理论的建立是在砂的基础上，因此使用到粉土，黏性土中还需要进一步改进和修正。

5.2 展 望

从流固耦合现有理论中存在的问题中可以看出现有耦合理论还有很多不足，可做如下改进：

1)对计算域之间耦合的理论需要进一步改进，使流体域和固体域之间既满足力的平衡条件也要满足变形协调条件。

2)可以将参数之间的耦合应用到其他3种耦合形式中，例如Terzaghi一维固结理论只是满足土体变形和液相排出体积的协调关系，而引入应力和渗透系数的关系后就可以考虑土体变形的非线性和渗流的非线性，使其更符合工程实际。

3)对本构关系耦合的研究较少，仅局限于对管涌、液化等问题的研究。以后，可以将本构关系的流固耦合推广到普通的固结问题。