作者：朱永生

版权申明：本文版权归ITASCA中国公司所有，引用或转载请注明出处。

 

在土木、水力水电等工程领域，作为一种流体介质，水文条件可以使得岩土体力学特性复杂化，并通常作为一种不利因素影响工程对象的稳定性性质。水在以参与物理化学作用的方式改造岩土体物质结构、弱化强度条件的同时，还可能与固体结构发生非常复杂的耦合作用，综合引致岩土工程稳定性条件发生改变、特别是会给工程本身及其影响区带来潜在安全问题。因此，水文地质条件勘探以及其对工程稳定性条件影响评价是设计分析需涉及的一项基本且极为重要的工作内容。目前而言，试验测试结合理论模型与计算分析是目前解决该类型问题采用的较为普遍的主要解决手段。

 

数值模拟技术显然是用于开展水作用问题研究的重要计算分析方法之一。到目前为止，ITASCA系列软件以岩土体为研究对象已形成包含等效连续介质力学、非连续介质力学方法在内的两套水力学模拟技术，前一方法以FLAC/FLAC3D系列软件为代表，而UDEC/3DEC、PFC系列软件则侧重用于描述裂隙岩体的水力学特性与行为。需要提及的是，由于PFC软件同时兼有表征岩体结构及模拟破裂的技术特点，因此特别适用于开展裂隙-水耦合作用致次生破裂及扩展力学机制研究，如石油行业较为普遍的利用PFC帮助深化认识压裂机理、从而指导优化现场压裂方案。不过，PFC在裂隙-水耦合作用的力学描述及程序实现这一关键环节具有过于宽泛的开放性，对应用者的专业理论和程序开发能力提出较为苛刻的要求。

 

以ITASCA流体模拟技术体系为背景，本系列推文计划针对其中的等效连续介质流体力学分析技术开展系统性介绍，重点涉及原理、计算模式、应力-渗流耦合作用分析、耦合分析模式选择及分析效率讨论等环节。本期推文侧重介绍如下内容：

1.        水作用问题的力学描述方法

2.        基于Biot理论的流固耦合力学行为描述

3.        ITASCA系列软件等效连续介质流体分析技术功能概述

4.        计算模式与水力模型属性定义

 

1.     水作用问题的力学描述方法

与解析方法相比，除涉及更为复杂的计算数学理论外，数值模拟方法还需要侧重对水作用问题力学原理实现合理的模型描述，特别体现在：

l      流体流动模型

l      应力-渗流耦合作用机制的描述

 

1.       流体流动模型

岩土体作为水介质的赋存条件，水在其中的流动行为首先显然与岩土体的细观结构特征息息相关。

 

由于土体细观结构由颗粒组成、一般视为均匀的孔隙介质，基于等效连续介质假定的水流动模型最早被引入用于描述土体的宏观水力学特性和行为，如达西定律采用渗透系数来描述土体在给定方向的渗透特性，水介质在其中的流动速度与渗透系数和水力梯度线性相关。

 

对于岩体而言，其中的岩块和裂隙是形成岩体的两个基本要素，岩块由于内部结构致密通常具有极低的透水性能，因此裂隙成为水介质流动的主要通道。由此可见，与土体不同显著不同的是，岩体水力学特性的力学描述主要在于反映其中的裂隙作用。自建立平行板流动立方定律以后，岩体结构描述方法如基于结构面样本数据统计分析的节理网络模拟等技术推动了离散网络介质流体分析模型的发展，该类型模型直接模拟岩体中的裂隙结构、并忽略岩块的渗透性。Barenblatt于1960年首次提出适用于裂隙岩体的双重介质模型，除裂隙作为渗透通道外，该模型同时视岩块作为多孔介质也可以起到导水作用，实现水介质在裂隙与岩块直接的物质交换。

 

不过，离散网络/双重介质流体分析模型要求的岩体结构描述参数难以通过常规地质勘探手段直接获取，在巨量模型数据管理和裂接触运算效率等方面目前也客观存在无法逾越的瓶颈问题，因此限制了方法在大尺度岩体工程应用层面对相关问题的解决能力。考虑到裂隙分布一般具有优势性和成组性，其决定的岩体水力学特性在宏观尺度往往近似表现为均质各向同性或正交各向异性，此时采用等效连续介质模型开展水作用问题分析仍旧具有良好的适用性。

 

为研究某水电站坝址区边坡运行期蓄水作用及其边坡稳定性，对建立图1所示FLAC3D等效连续介质渗流-耦合作用分析模型，边坡被断层、层面等构造面划分为具有不同渗透特性的岩性地质单元。特别地，各类构造面均被处理为有厚度的数值模型单元，同时考虑其平面方向与法向渗透特性的现实差异，水介质在构造面内的流动行为采用正交各向异性水力模型描述。总体而言，等效连续介质模型迄今仍然是研究岩土体水力作用的主流方法。

 

图1 边坡地质分区及构造面模拟

 

2.       应力-渗流耦合作用机制的描述

由于岩土体中的固体部分以及水介质本身均具有变形能力，在外界条件变化扰动（力学扰动或流动边界变化）过程中，固体介质和水介质可能会发生不同程度的流-固耦合作用（应力-渗流耦合作用）。如低渗透特性软弱地基中基础施工（力学扰动）会导致土体骨架产生压缩变形，而孔隙中的水介质则同时由于未能及时排水与孔隙率降低受压缩作用形成超孔压，这部分超孔压又会反过来影响土骨架变形，周而复始表现出复杂的流-固耦合作用过程。因此，对耦合作用的力学描述是水作用问题模拟分析需重点解决的另一技术环节。

 

以有效应力原理作为背景，下一小节将进一步通过介绍孔隙介质流体质量连续性方程的分项构成来说明ITASCA软件等效连续介质流-固耦合作用模拟的理论依据即Biot原理。

 

2.       基于Boit理论的流-固耦合力学行为描述

以土体作为研究对象，太沙基（Terzaghi）和伦杜立克(Rendulic)通过试验和理论分析认为以土体为代表的孔隙固体介质中的有效应力为(符号均用张量表示)：

(1)

式中：为土体中的有效应力；P为孔隙压力，在FLAC3D程序中默认该变量压为正；为液相（水）介质在单位面积固体介质截面内的有效面积系数，一般按不利原则考虑为=1；m为单位向量，。

 

考虑一单位体积土体单元，假定土体中的总应力分量之球应力张量分量及孔压增量分别为、，且认为两者同时发生，则由此引致土元中的固态相瞬态体积变化由以下三部分构成：

1)   土元骨架轮廓体变增量。该部分体变由有效球应力张量引致：

(2)

2)  孔压变化导致的土颗粒体变增量。单位体积土元中颗粒体积总量为(1-)(为孔隙率)，则由孔压变化导致的土颗粒体积变化可写为：

(3)

3)   有效应力导致的土颗粒体变增量。土骨架在细观角度由土颗粒组成，颗粒间有效接触力变化可导致颗粒发生体变。考虑单位土元断面，设断面上有效应力改变量为，忽略颗粒间电荷力，则断面内颗粒间有效接触应力改变量可以写为：

认为颗粒间有效接触应力在土元骨架内均匀分布，则由于该项应力改变量导致的土颗粒体变为：

(4)

式中，为土体排水体积模量。

 

此外，单位土元中液相由于孔压变化导致的体变增量可以写成：

(5)

式中，为流体的体积模量。

 

设土元中的达西(Darcy)流量为，且进一步考虑体积源项及其流体密度变化行为，将式(2) ~ 式(5)各项叠加得到土元中流体（液相）的质量连续性方程：

(6)

式中：为流体密度。

 

对式(6)做变换可以写成：

(7)

式中：为Biot模量，；为Biot系数，。

 

对于饱和-非饱和土体，假定孔隙气和大气相通，且孔隙流体和孔隙气体变增量与各自的饱和度成比例，则式(2)~(5)需要考虑饱和度影响。考虑饱和度变化，得到饱和-非饱和孔隙流体的质量连续性方程：

(8)

式中：为流体饱和度。

 

若进一步考虑温度参与耦合作用，则得到流体在温度-渗流-应力(THM)三相耦合作用条件下的质量连续性方程：

(9)

式中：、分别为温度及其土体的不排水温度系数。

ITASCA等效连续介质多相耦合分析技术所采用的流体质量连续性方程服从定义：

 

 

 

对比式(9)、式(10)可见，以FLAC3D程序为代表的ITASCA等效连续介质多相耦合分析技术未考虑流体密度在作用过程中可以存在的变化行为。

3.         ITASCA系列软件等效连续介质流体分析技术功能概述

综合以上讨论可知，ITASCA系列软件等效连续介质模型为达西流相关渗流、应力-渗流耦合作用问题提供解决方案，模型假定固体材料为连续介质，且模型中的流体流动行为以及流体-固体间耦合作用行为分别服从经典达西定律和Biot固结理论。

 

以FLAC3D为例，ITASCA系列软件等效连续介质模型应力-渗流耦合分析技术具有一般性特点和功能：

1) 流体在固体介质中的流动行为服从达西流动定律，并同时满足质量平衡及其相容方程，流动介质为流体（不太适合用于气体流动分析），且认为流动介质可压缩；

2) 提供三种流动模型，即各向同性、各向异性及其流体开挖模型，其中开挖模型常用于模拟模型中的不透水材料。不同数值单元即ZONE可灵活选择不同的流动模型；

3) 多样化流体边界条件，如孔压、流量、溢出和不透水边界等；特别地，流量边界（流出或流入）可以以点源或体源的方式放置于模型内部，该边界条件常见于基坑井点降水模拟；

4) 针对流体饱和、饱和-非饱和赋存状态的不同，运算器分别提供显式和隐式两种算法加以区别求解。其中，隐式算法仅适用于饱和流体；

5)   在流体计算过程中，固体力学分析及其温度分析等模式均可同步参与进行耦合作用：

Ø   流体-固体耦合作用时互为影响服从Biot定律，其中固体变形涵盖骨架及其土颗粒变形，土颗粒变形可采用Biot系数α加以描述，但假定材料孔隙率和渗透性在耦合作用过程中为常量；

Ø   在温度-流体-固体三场复杂耦合分析中，通过不排水温度系数β来反映温度对流体和固体骨架的热膨胀效应。不过，分析技术并未反映温度作用效应可以导致固体渗透特性（如渗透率）以及流体性质（如密度）的变化。现实中存在对这些复杂行为的描述性需求，可考虑采用FISH语言进行自定义开发。

 

4.          计算模式与水力模型属性定义

以FLAC3D V7.0程序为例，本小节简要介绍ITASCA系列软件等效连续介质力学方法水作用模拟支持的可选计算模式，以及相应的模型水力学属性定义。

 

1.       计算模式

由水力学原理可知，流体压力可分解为静水压力和动水压力两个构成分量。在流体与固体两相介质中，这两部分作用力分量综合对固体介质产生浮力与渗透力作用。由于岩土体中水介质流速要慢的多，为简化分析难度，很多规范规程针对常规工程计算往往忽略对动水压力的考虑。

 

视是否考虑动水压力作用，FLAC3D程序水作用模拟区分为无渗流与渗流两种计算模式。无渗流计算模式仅用于生成对应于指定水位面条件的静水压力场，无需特别的命令来触发；渗流模式作为程序内置功能模块之一，必须通过执行命令model config fluid后方可调用。在渗流分析模式下，必须对模型中的单元赋予渗流模型即流体流动本构模型，程序提供各向同性、正交各向异性及其空流体三种本构模型供选择，分别对应于命令zone fluid cmodel assign isotropic/anisotropic/null。

 

2.       模型属性定义

由于计算模式不同，模型属性定义显然也存在明显差别。表1比较了两种计算模式所需定义模型属性以及命令格式的不同。

 

表1 模型属性定义

属性	关键字	依附对象	计算模式及命令格式
			无渗流模式	渗流模式
流体密度ρf (1)	视计算模式不同	单元	zone water density	zone initialize fluid-density
固体密度ρ (2)	density	单元	zone initialize
渗透率k (3)	permeability等	单元	无需定义	zone fluid property
孔隙率n	porosity (默认为0.5)	单元		zone fluid property
流体体积模量Kf (4)	fluid-modulus	节点		zone gridpoint initialize
比奥系数α(4)	biot (默认为1.0)	节点		zone fluid property
比奥模量M (4)	biot-modulus	节点		zone gridpoint initialize
饱和度s	saturation (默认为1.0)	节点		zone gridpoint initialize
流体抗拉强度σt (5)	fluid-tension (默认为-1e-15)	节点		zone gridpoint initialize

 

(1)  无渗流计算模式需搭配水位面人工定义功能（见命令zone water）使用；

(2)  在无渗流计算模式中，固体密度为岩土体天然密度，而渗流计算模式中的固体密度指岩土体干密度；

(3)  FLAC3D用于描述岩土体渗透特性的参数为渗透率k，其与传统渗透系数kh的关系为k=kh/，为流体重度；

(4)  当比奥系数即假定固体颗粒不可压缩时，流体介质变形特性采用体积模量Kf来定义；否则（，固体颗粒可压缩），采用比奥模量M描述流体介质变形特性；

(5) 在细颗粒土中，由于毛细管力、电化学力的作用，孔隙水在非饱和条件下可以维持一定水平的张力。在FLAC3D程序中，流体抗拉强度σt用于帮助识别孔隙介质饱和-非饱和状态。具体而言，当孔隙流体中张力超过σt时进入非饱和状态，渗透率此时遵循内置非饱和水土特征方程。不过，程序中流体张力缘由土骨架膨胀所致，内在机制本质上不同于上述毛细管力、电化学力等实际非饱和受力原理。或者说，就非饱和流体力学行为描述而言，FLAC3D采用的背景理论是基于一定假设的近似，并不追求对其精确模拟。

 

参考文献：

[1]      FLAC3D 7.0 Manual. Itasca Consulting Group, Inc.