0 引言
岩土工程施工阶段引排水改变地下水分布，令土体有效应力重新分布，有效应力场反过来又影响地下水运动，这种应力场和渗流场相互作用的现象被称为流固耦合［1］。20 世纪中期，比奥按弹性力学理论推导出方程组，由于求解困难，当时未引起重视。近些年来，相关研
究成果显著。一方面，专家学者研制出能模拟流固耦合的模型试验系统。比如: 文献［2］通过模型试验得到海底隧道施工中围岩力学行为和渗流场变化规律; 文献［3］对比了模型试验和数值计算之间的差异。另一方面，基于流固耦合的数值计算如火如荼。比如: 文献［4］将流固耦合用于双圆异性断面盾构隧道施工分析; 文献［5］研究了基坑开挖时支护变形和孔隙水压的分布。目前，研究正向着裂隙介质非均匀流、非饱和渗流以及多场耦合方向发展［6 － 7］。
虽然成果丰硕，但在流固耦合基础应用方面还存在一些问题。我国流固耦合数值模拟主要基于FLAC3D平台，软件内置4 种渗流模型。基坑、隧道开挖是一个复杂、循环往复的过程，既存在土体开挖卸荷，也存在引排水改变渗流场，难以根据用户手册选择适合的模型，很多文献选择模型的标准也不尽相同。文献［8］模拟地铁开挖时，先计算应力场，再进行流固耦合或渗流计算;文献［9］模拟基坑开挖时，先等渗流场稳定后，再进行应力计算。上述处理方法均有一定道理，到目前尚未有文献针对4 种模型受力特点进行详实论证。为了解4 种模型的受力本质，确保计算结果能指导工程，本文将以一各项同性弹性立方体为例，对各模型的受力特点进行分析。
1 流固耦合概述
1． 1 相关理论
如图1 所示，取一微元体，建立基本微分平衡方程:
式中: σ'x、σ'y、σ'z、τxy、τxz、τyz为有效应力在x、y、z 轴上的分量; u 为孔隙水压力; γ'为体积力。
图1 微元体受力示意图
根据胡可定律，可以用位移us、vs、ws 表示σ'x、σ'y、σ'z、τxy、τxz、τyz，由于式( 1) 共有4 个未知数( us、vs、ws、u) ，无法获得唯一解。假设水和土颗粒不可压缩，土微元体单位时间内体积的变化量等于水的流入量与流出量之差，根据达西渗透定律得到一连续方程。4 个未知数对应4 个方程组，理论上可解出未知数。前人已就具体推导过程做了详细的说
明［10 － 11］，本文不再展开讲解。
1． 2 计算原理
FLAC3D 计算过程简化如图2 所示，左边表示应力场计算的过程，右边表示考虑渗流场后的新增步骤。表1 为FLAC3D 用户手册渗流模型选择标准，主要是根据特征时间和扰动类型确定的。
图2 FLAC3D 建模流程
软件根据运动方程、平衡方程、本构方程、相容方程以及边界条件进行求解［1］。
1) 运动方程是由达西定律v = k•i 推导而来，渗流速度
式中: kil为绝对渗透系数; k∧ ( s) 为相对渗透系数; p为孔隙水压; pf为流体密度; gj为重力加速度。
2) 平衡方程是在小变形情况下，单位时间内微元体含水量变化值等于流入量与流出量之差:
式中: qυ为微元体流体流入量; ζ 为单元流体体积变化值。
3) 本文只针对饱和土体渗流模型，其本构方程为
式中: M 为比奥模量; α 为比奥系数; ε 为应力场引起的体积应变。
4) 应变率和速度梯度应满足相容方程:
表1 FLAC3D 渗流模式选择
1． 3 工程应用难点
FLAC3D 英文手册给出渗流特征时间的计算公式，特征时间可通过式( 6) —( 8) ［1］计算。
式中: Lc为特征长度; Vf为渗流区域体积; Af为渗流区域表面积; tc为特征时间; k 为FLAC3D 渗透系数;n 为孔隙率; Kf为流体体积模量; K 为围岩排水条件
体积模量; G 为围岩排水条件剪切模量; α 为比奥系数。
对于简单的一维渗流，容易得到特征长度和特征时间。如图3 所示，现有一厚度为20 m 的饱和土层，底部为一不透水刚体，顶部作用均布荷载，水通过顶部排除，所以特征长度取土体厚度，即20 m［12］; 但在实际工程中，地下水的运动是三维的，难以确定特征时间和渗流模型。
太沙基把土的渗透固结简化为一个装满水的容器，容器上方有一开小孔的活塞，活塞下面用一根弹簧支撑。当荷载刚作用上去时，小孔中的水来不及排出，由水和弹簧共同承担荷载。当荷载作用长时间后，多余的水分排出，弹簧承担全部的荷载［13］。本文的研究思路是将土的渗流固结过程分段，建立各段与渗流模型之间的对应关系。例如对应排水固结阶段的渗流模型，理论上不会出现超孔隙水压。
图3 一维渗流算例
2 模型设计
2． 1 模型试算
为了保证建模正确性，现进行模型试算。根据相关文献［3、14］，建立青岛胶州湾海底隧道流固耦合模型，如图4 所示。由于部分参数( 开挖时间、加固区尺寸)未在文献中找到，所以未知参数按经验取值。
图4 青岛胶州湾海底隧道模型
图5 与图6 为试算结果同原文的对比，结果存在一定的差异，主要是由于部分参数取值不同而引起的，但围岩竖向位移、内力均在同一数量级，变形收敛趋势是相同的。文献［3］还就数值计算结果同模型试验进行了对比，若不考虑裂隙水、围岩产状等因素，数值解和试验结果相当接近。
图5 青岛胶州湾海底隧道试算结果
图6 青岛胶州湾海底隧道试算结果
流固耦合计算费时，可以沿隧道中心面取一半建模，以加快计算速度，如图7 所示。需要注意隧道中心面处孔隙水压力边界条件，通常设置为不透水边界，即孔隙水压是可变化的。
图7 圆形隧道模型
2． 2 参数设计
本文的目的是为了分析4 种渗流模式的受力本质，从而与实际工程建立联系，就如同钢桁梁桥总体受力虽然复杂，但每根钢杆依然准守基本的力学公式。
本文以简单模型为例，更便于分析说明，所得结论，亦适用于复杂模型。
现设一100 m × 100 m × 10 m 立方体，水位高出地面10 m。初始平衡以后，在顶部施加一竖向均布荷载q，大小为500 kPa。在模型中心从上到下依次取点A—J 10 个监测点，监测其初始竖向总应力σ0、初始孔隙水压p0、加载后竖向总应力σ1、加载后孔隙水压p1。模型示意图见图8。
图8 模型示意图
为了防止模型屈服破坏，力学模式采用理想弹性模型，渗流模式采用各项同性渗流模型，参数取整，详见表2 和表3。
表2 模型物理力学参数
表3 模型水力学参数
2． 3 工况设计
FLAC3D 内置4 种渗流模型，各模型力学进程、流体进程、流体体积模量Kf有所差别。通常流体不考虑压缩变形，土体的变形是由于流体迁移和土颗粒骨架空隙间的变形。FLAC3D 将流体考虑成可压缩变形的，相关资料请参考FLAC3D 用户手册。现根据不同进程和流体模量设置4 种工况，见表4 。
表4 计算工况
3 计算结果及分析
3． 1 工况1 计算结果
工况1 是在初始平衡后，顶部施加荷载500 kPa，维持力学进程开启和流体进程关闭，流体体积模量Kf为0，计算至平衡。工况1 应力竖向分布曲线如图9所示。由图9 可知: 曲线p0
和曲线p1重合，说明加载后孔隙水压不变; 有效应力增大约为500 kPa，即外界施加的荷载由土颗粒骨架承担，水不分担压力。
图9 工况1 应力竖向分布曲线
3． 2 工况2 计算结果
工况2 同工况1 类似，维持力学进程开启和流体进程关闭。与工况1 不同的是，此时需要将流体体积模量Kf设置为真实值。工况2 应力竖向分布曲线见图10。由图10 可知，外界荷载一部分由土颗粒骨架传递，一部分由流体承担。
在试算时还发现流体承担的荷载大小，同流体和固体的刚度比Ｒk有关。新增1 组对比模型，将流体的体积模量降低10 倍，以分析刚度比Ｒk对结果的影响。
刚度比
式中: Kf为流体体积模量; n 为孔隙率; K 为固体体积模量; G 为固体切变模量。
图10 工况2 应力竖向分布曲线
图10 中曲线σ1( 柔性) 和曲线σ1( 刚性) 重合。通过对比发现，流体分担外荷载比例同刚度比Ｒk成正比。因此可作出一个假设: 土颗粒和流体分担外力原理类似图11。刚度比Ｒk就是2 根弹簧的弹性系数k比值，当刚度比Ｒk等于0 时，代表流体的弹簧弹性系数k 等于0，外荷载全由另一根弹簧———土颗粒骨架全部承担。工况1 正是刚度Ｒk等于0 的特殊状态。
图11 流体－ 固体刚度比示意图
3． 3 工况3 计算结果
由于工况3 是模拟孔隙水压变化引起的内力重分布，所以初始平衡之后，不再施加500 kPa 荷载，而是将顶部孔隙水压提高500 kPa。另外，工况3 和工况4需要设置时间，这里将时间取一个大值，使模型渗流运动场趋于稳定。
工况3 计算过程可视作2 个独立的步骤: 第1 步是孔隙水压变化引起渗流场重分布，有效应力场不发生变化; 第2 步是渗流场重分布后，有效应力场变化发生。工况3 应力竖向分布曲线见图12。图12 中曲线p1( M off /F on) 和曲线p1( M on /F off) 重合，说明有效引力场改变并未引起孔隙水压力重分布。结果再次证明本文3． 2 节假设模型的正确性———Ｒk = 0 时，流体不参与受力的。
图12 工况3 应力竖向分布曲线
在分析实际工程时，孔隙水压的改变通常是人为控制的。比如开挖隧道时，将开挖后的单元孔隙水压设置为0，围岩周围单元的孔隙水压会根据已知孔隙水压力边界条件重新分布。
3． 4 工况4 计算结果
工况4 中流体参与分担外荷载，并且分担的数值与时间有关。现在设置2 组对比模型，计算时间分别为6 000 s 和60 000 s，其余参数相同。工况4 应力竖向分布曲线见图13。图13 中曲线σ1( 6 000 s) 和曲线σ1( 60 000 s) 重合，随着时间的的增加，超孔隙水压消散，有效应力增大。结合图9 和图10，工况4 可以近似看作工况1 和工况2 的中间状态。当计算时间足够长时，结果近似于工况1; 当计算时间十分短时，工况4的结果近似于工况2 结果。
由此得到这样一个结论———模型D 适用于任何情况。分析透水性好的砂土排水固结或是透水性差的黏土不排水固结，都可以选用模型D，只要时间按真实情况取值即可; 但流固耦合相当费时，所以模型D 只是“理论上”适用于各种情况。根据大量的试算还发现，不仅合理的渗流模型能节约时间，采取合理的网格形状，避免单元尺寸过小，亦能减小单元数量和增大时间步time-step，从而缩短时间。
图13 工况4 应力竖向分布曲线
4 结论与建议
1) 模型A、模型C 外荷载全部由土颗粒骨架承担，模型B、模型D 流体也分担外荷载。
2) 流体分担外荷载比例与刚度比Ｒk成正比，模型D 流体分担比例还与作用时间有关。如透水砂层在外荷载作用下长期变形，宜选用模型A; 透水不良黏土层突然受外力作用，宜选用模型B; 抽水引起基础变形，宜选用模型C; 黏土层排水固结过程可采用模型D。如果将土的渗透固结过程看作是用手在压一个装满水的容器的话［13］，那么FLAC3D 中4 个模式对应关系参考图14( 笔者注: 图11 与图14的模型是在大量试算基础上，为了便于理解而假设的定性模型) 。
图14 对应关系示意图
3) FLAC3D 流固耦合正是通过调整刚度比Ｒk，打开和关闭力学－ 流体进程来实现的; 但是由于流固耦合十分费时，一般通过调整流体体积模量Kf或比奥模量M 减少计算时间。此外，手册还指出计算渗流稳定状态时，流体比奥模量要满足M ＞ αLzρωg /n［1］，以确保流体自由面稳定。相关理论还待作进一步研究。
摘自：隧道建设