2.2 流固耦合理论

流固耦合理论最早用于研究土体的固结问题。1925年由太沙基提出，并由比奥加以发展和完善。载荷作用下的土体固结沉降变形不是瞬间完成的，而是按照一定的变化速率逐渐下沉。固结沉降是土体对荷载变化的逐渐适应。对多孔介质而言，在外荷载作用下，土体中将产生超静孔隙水压力，引起土体孔隙中的水力梯度发生改变，进而引起土体孔隙中的水发生流动。外荷载卸载后，由于孔隙弹性的作用，因外加荷载而变小的孔隙，在一定程度上将逐步恢复变大。在孔隙介质周围存在流体补给源情况下，逐步恢复变大的孔隙将再次饱和；否则，孔隙介质中的孔隙将保持非饱和状态。对大多数土体而言，孔隙非线性特征十分明显，其固结变形一般都是不可逆的。近年来，流固耦合理论在地基固结沉降变形、深埋高压水工隧洞稳定性、边坡稳定性分析、注水采油、水库诱发地震和核废料地质存储库研究中得到了越来越多的重视，其应用前景也越来越广阔。

在渗流和外荷载作用下，孔隙介质体内存在两种响应：一是固体介质部分产生变形，并伴随应力重分布的产生；二是孔隙介质中流体压力产生改变，引起流体扩散运动。这两种响应相互影响，直至研究区域渗流系统以及内力系统重新达到平衡。因此，流固耦合问题分析的实质就是要分析应力和渗流各自变化对饱和连续孔隙介质系统整体的影响。连续多孔介质流固耦合问题求解的控制方程包括平衡方程、连续方程和状态方程。

2.2.1 平衡方程

取多孔介质六面微元体进行研究，则该六面体上作用有应力σ_ij和应变ε_ij。在静止条件下，作用在六面体各个面上的应力及六面体体内的体力f_i将组成一个平衡体系，其平衡方程为

连续孔隙介质在应力作用产生相应的变形，描述变形关系的几何方程为

式中：u_i为位移量。

弹性介质应力和变形之间满足的本构关系为

式中:σ′_ij为有效应力;G为剪切模量;λ为拉梅常数。

根据有效应力原理，可压缩孔隙连续介质中孔隙压力和固体骨架上的有效应力之间满足下式

式中:α为比奥系数； p为孔隙内的流体压力。

联立式（2.35）和式（2.36），得

将式（2.37）代入式（2.33），可得弹性孔隙介质的平衡方程

由几何方程消去应变ε_ij，得到位移表达的平衡微分方程

以位移分量u_x、u_y、u_z表示的平衡方程为

式（2.40）所示的平衡方程中包含3个未知位移量和1个孔隙压力p共4个未知量。

2.2.2 连续方程

对于饱和多孔介质，存在两个位移矢量即固体骨架的位移矢量u_s和流体质点的位移矢量u_f，相应地有两个速度ν_s和ν_f。下标s和f分别对应于固体和流体的量。流体相对于固体的速度ν_r为

根据裘布依-福希海默关系式，孔隙介质中流体的渗流速度φν_r（φ为孔隙度）。于是，达西方程可改写成

式中:κ为渗透率张量;μ_f为黏滞系数。

在不考虑源（汇）情况下，流体连续方程可写为

将式 （2.43）展开，略去二阶小量ν_s·▽，可得以速度表达的流体连续方程

类似地，可以写出固体的连续方程

将式 （2.45）展开，略去二阶小量ν_s·▽，并将各项乘以ρf/_ρs，得

考虑到孔隙介质骨架体积应变为ε_v=u_i,i，则有

式（2.45）和式（2.47）相加可得饱和多孔连续介质体整体连续性方程

或

2.2.3 状态方程

孔隙介质在内、外应力的作用下，其密实度状态发生一定程度的改变。描述孔隙介质状态的物理量有固体骨架的孔隙度、流体的密度，其状态方程分别为

式中:c_f、c_φ分别为流体和固体孔隙的压缩系数。

随着固体介质孔隙度的变化，固体密度也会发生相应的改变。孔隙介质固体密度与孔压和应力关系为

式（2.52）中右端第二项表示孔隙流体压力所引起的固体密度的变化，第三项表示有效应力引起的固体密度的变化。将ρ_s对t求导，可得

2.2.4 控制方程

将达西方程式（2.42）和状态变化方程式（2.53）代入整体连续方程式（2.49）得孔隙介质渗流场方程为

式中：K为饱和多孔介质的整体体积模量；κ为渗透率张量。

式（2.54）与平衡微分方程（2.40）共同构成孔隙介质流固耦合控制方程。

控制方程组式（2.55）中共有4个未知量，其中3个为位移未知量（u_i,i=1,2,3），1个为孔隙压力未知量（p）。

2.2.5 耦合本构方程

对于等温多孔连续介质，基于孔隙介质线弹性假设，孔隙介质总的应变增量等于应力引起的应变增量和孔隙压力引起的应变增量的代数和，即

式中:σ′_ij0、σ′_kk0分别为初始应力和初始体积应力;p₀为初始孔隙压力；G为剪切模量;K_m为孔隙介质基质体积模量;ν为泊松比;α为比奥系数。

对式（2.56）进行逆变换得可压缩孔隙连续介质力学本构方程

式中：K₀为参考状态下的孔隙连续介质体积模量。

式（2.57）与式（2.37）本质相同，形式不同，式（2.57）为增量形式的耦合本构方程。

对于等温多孔连续介质，荷载及渗流边界条件改变下，孔隙压力增量等于孔隙介质中流体质量变化引起的孔隙压力增量和孔隙介质骨架变形引起的孔隙压力增量的代数和，即

式中：M为比奥模量；m为孔隙中的流体质量增量；ρ_fl0为参考状态下的流体密度；其余符号意义同前。

式（2.57）与式（2.58）共同组成饱和多孔连续介质的增量型耦合本构方程组。

2.2.6 定解条件

控制方程组式（2.57）和式（2.58）中包含求解位移分量的3个方程和1个求解孔隙压力的方程。依据弹性理论，符合控制方程组式（2.57）和式（2.58）的解答不唯一。如果需要获得唯一解，必须满足一定的定解条件。对于简单的孔隙介质流固耦合问题，可以得到解析解。对于大多数孔隙介质而言，在渗流和应力环境共同作用下，其解析解难以获得，在这种情况下，数值方法提供了一条有效的求解途径。

控制方程组式（2.57）和式（2.58）的定解条件，包括初始条件和边界条件，边界条件包括力学边界和渗流边界。

初始条件：研究域初始状态孔隙压力p₀和有效应力σ′₀；流体和固体初始密度ρ_f0和ρ_s0。

力学边界条件包括已知应力边界或已知位移边界，其表达式如下：

已知位移边界： u｜_s=u₀(x,t)

已知应力边界： σ｜_s=σ₀(x,t)

渗流边界包括已知流量边界和已知压力边界，其数学表达如下:

已知流量边界：Q|_s=Q₀(x,t)

已知压力边界：p|_s=p₀(x,t)

2.2.7 控制方程求解

渗流场和应力场全耦合的过程是一个动态的过程，其求解需借助有限元等数值方法完成。为实现耦合过程中渗流场和应力场之间的相互耦合，可采用间接耦合法和直接耦合法。直接耦合法是在求解流固耦合控制方程式（2.57）和式（2.58）过程中将未知变量位移和孔隙压力纳入同一组方程中，在求解位移方程中直接增加渗流引起的节点不平衡力增量，而在求解节点孔隙压力方程中增加节点位移改变引起的节点流量增量。

1.空间域离散

考虑到动态耦合问题具有时间分段计算、荷载逐级施加以及自由面迭代的非线性等特性，需建立增量形式的有限元分析格式。假定在某级荷载增量下，产生的位移增量为Δu，增量形式的位移计算有限元方程控制方程组为

其中

式中：Δp为未知孔隙水压力增量列阵；Δp^* 为已知节点孔隙水压力增量；G′为单元位移自由度选择矩阵；B′为应变矩阵：L={1,1,1,0,0,0}^T；N={N₁,N₂,…,N_i,…,N₈}；N′为节点位移形函数矩阵；，［Δf］^e分别为单元已知边界面力和体积力增量列阵，一般情况下有［Δf］^e={0,0,-Δγ}^T；ΔF_n为已知节点集中力列阵。

耦合分析的渗流控制方程组为

其中

式中：u、p分别为未知节点位移和孔隙水压力列阵；为第一类渗流边界节点已知水压力列阵；q^e为单元已知边界法向流量；G为单元孔隙水压力自由度选择矩阵；k^e为单元渗透系数张量，且k^e=［k_ij］，是应力张量的函数。

需要指出的是，K″、Q₀ 两项是与自由面相关的，对于承压流而言，此两项不存在；但对于无压渗流，对自由面边界单元需要计算K″，此项对内部单元而言为0。

2.时间域离散

连续性方程进行空间离散后所得式（2.60）包含了节点位移、孔隙水压力对时间的微分项，因此，需要将式（2.60）在时间域内进一步离散。不妨取动态计算时间步长为Δt，引入时间因子α，式（2.60）可改写为

式中:为的转置矩阵。

为了使求解迭代过程稳定收敛，通常取0.5≤ξ≤1.0。假定位移、孔隙水压力在时间段t~t+Δt内呈线性变化，则有下列各式成立：

代入式（2.61）并整理可得

考虑到一般情况下，边界补给流量Q在时间段t～t+Δt内变化不大，式（2.62）可简化为

令

于是，式（2.63）可改写为

式（2.64）即为渗流连续性方程在时间域和空间域内离散后得到的增量形式有限元方程。

3.动态全耦合有限元方程组及迭代求解

联合式（2.59）和式（2.64），可获得饱和孔隙介质渗流场与应力场弹性动态耦合定解问题的增量形式的有限元方程组：

由于为应力张量的函数，对于某一级荷载下的每一个计算时步Δt而言，当研究域内存在自由面位置，且事先不确定时，中包含的K″是随自由面位置的变动而改变的，因此，动态耦合有限元方程组（2.65）是一个强烈的非线性方程组，在每一个计算时步Δt内都需要迭代求解。

动态耦合有限元方程组（2.65）的直接迭代格式如下：

式中上标k、k+1分别为各分项在第k、k+1次迭代结束时相应的值。式（2.66）即为动态分析中实际求解的有限元代数方程组。求解过程在时间段t～t+Δt内迭代稳定之后，即可近似确定t+Δt时刻的自由面位置，并获得该时刻的位移、孔隙水压力增量；依据该时刻的应力和自由面位置重新计算系数矩阵K，以迭加后的节点位移值、孔隙水压力值作为初始条件，开始下一时步的计算，直到本级荷载终了时刻为止。

该时步计算完成后，施加荷载增量，进行荷载增量步循环，直到加载结束。需要指出的是，对承压流而言，K″是不存在的，此时式（2.66）可简化为

式中:为中扣除初流量项之后的等效流量增量列阵。

由于式（2.67）考虑了渗透张量与应力的耦合关系，即使采用弹性本构模型，式（2.67）也是一个非线性方程组，仍需在每一个计算时步内迭代求解。当然与式（2.66）相比，式（2.67）消除了自由面边界非线性的影响，非线性仅来源于随应力的变化，非线性程度减弱了，因此数值分析的计算量也减小了，收敛性也可以保证。