作者：朱永生

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

 

本文案例代码可通过百度网盘下载：

链接：https://pan.baidu.com/s/1g572Sfq1aj22kapbi43mCw

提取码：gfj4

 

本文为前期推文《ITASCA等效连续介质流体分析技术(一)：方法原理》的续篇。

 

在前期推文中，以FLAC3D等效连续介质流体分析功能为例，初步介绍了模拟岩土体中流体作用的两种分析模式即渗流模式和无渗流模式的格子基本特点和一般性区别，特别是两种分析模式对流体属性定义要求的不同。

 

本期推文将进一步介绍两种分析模式在FLAC3D V7.0程序中的使用方法，意在体现不同方法实现在具体环节的差别。除流体属性定义外，这些环节还主要包括边界条件的设定，应力初始化在程序中的背景意义以及在开展具体工程分析时的处理方式等。总体而言，本期推文讨论内容涉及如下几个问题：

1       渗流分析模式下的边界条件

2       不同分析模式的实现案例 

3       讨论

 

一、      渗流分析模式下的边界条件

无渗流分析模式不涉及到对流体流动的模拟，模型中的孔压分布状态完全由人为指定，因此该分析模式也不涉及流体相关边界条件的设定。

 

此处讨论为满足考虑流体流动模拟的渗流分析模式所要求的边界条件的设定。包括FLAC、DEC在内ITASCA连续介质流体分析技术为渗流模式分析提供孔压边界和源项边界两大类边界条件。应特别留意到，不同边界条件设定的适用模型元素可存在差别，如孔压边界的设定对象仅为单元的节点，而源项边界则较为灵活。

 

1.       孔压边界

如前所述，孔压边界仅用于设定单元节点处的孔压条件。基于方便用户控制角度考虑，程序为满足这一要求主要提供了两种定义方法（命令）：

l        方法1：zone gridpoint fix pore-pressure value grad gx gy gz range …

l         方法2：zone face apply pore-pressure value grad gx gy gz range …

 

其中，方法1的操作对象直接为节点，而方法2是以节点所在的面作为媒介实现对其孔压条件的间接定义；gx、gy、gz分别为孔压沿模型坐标系各方向的分布梯度。以上方法均实现对指定操作范围内的节点处孔压水平实现量化定义：

        孔压 =  value + x * gx + y * gy + z * gz

式中，x、y、z为节点坐标。

 

方法2和方法1的主要区别在于，前者还可以帮助建立孔压水平与运算时间之间的关联，具体见手册中对fish、table、vary等关键字的描述。

 

与孔压边界定义对应地，程序同样提供移除孔压边界设定的命令：

l        方法1：zone gridpoint free pore-pressure range …

l        方法2：zone face apply-remove pore-pressure range …

 

孔压边界的另一层含义用于反映节点的透水性质，对于经孔压设定操作的节点而言，必须允许流体在该节点处流入或流出以保证其孔压维持为定值，即为透水边界；对于未经孔压设定的节点，系统默认节点部位与外界不发生流量交换，即为不透水边界，因此孔压可随平衡条件而变化。总体而言，程序默认模型中所有边界节点为不透水边界，一旦经孔压设定操作后，则转变为可透水边界。

 

2.   源项边界

源项边界与流量直接相关，用于设定流体在指定模型元素部位处的流入或流出速率，模型元素可以是节点、单元或者模型网格的边界等。具体设定命令主要有：

l         对于节点：zone gridpoint fix well v1 range …

l         对于单元：zone apply well v2 range …

l         对于模型网格边界：

i.        zone face apply discharge v3 range …

ii.       zone face apply leakage v4 v5 range …

其中，v1表示在指定节点处的流体流入/流出速率，单位为体积/每单位时间；v2为流体在指定单元处的流入/流出速率，单位为单位体积单元的流体流入或流出体积/每单位时间，可见单元源项总量是v2与单元体积之乘积；v3表示流体在模型边界面（即单元表面）部位的流入或流出速度（方向与边界面垂直），单位为长度/每单位时间；v4、v5分别为渗漏孔压和渗漏系数，综合用于确定流体在模型面边界（即单元表面）部位的渗漏速率（方向与边界垂直）：

(1)

式中，为渗漏速率，单位为长度/每单位时间；p模型面边界（即单元表面）部位的孔压；v5的单位为长度^3/(力*时间)。

 

此外，v1、v2、v3为正时表示流体流入，反之则流体流出。

 

二、   不同分析模式的实现案例

前述内容大致介绍了不同分析模式对流体属性的定义要求，以及渗流模式条件下的边界条件的设定方法，为模拟分析奠定基本概念。

 

本小节引入简单案例说明利用不同分析模式开展流体分析的方法框架和基本流程。

 

1.       问题及模型说明

材料条件：考虑位于刚性地基上一半无限分布厚度为10m的均质土层，简单视其为弹性，体积模量K与剪切模量G分别为100MPa和30MPa，干重度和孔隙率n分别为1800kg/m3和0.5，渗透率k取1.e-12m2/Pa/s。分析中取流体密度与重力加速度g分别为10m/s2。

 

分析条件：假定土层初始处于干燥状态，考虑其侧压力系数为0.5714，施加重力模拟获得平衡状态；之后考虑水位上升至土层顶部，分析因水位变化引致的土层变形特征。

 

在上述条件下，土层顶部竖向变形的理论解为：

(2)

式中，为侧限变形模量，；为土层厚度。将模型条件带入上式可得到因水位变化导致的土层顶部竖向变形为1.786×10−3m。

 

利用FLAC3D V7.0程序，选择不同分析模式模拟这一问题，在展示方法实现过程的同时，验证不同分析模式的合理性。图1为分析所用模型，模型在x、y、z（竖直方向）轴向的长度分别为2m、1m和10m，相应的网格数为2、1与10个；图2为项目中数据文件的构成说明。

 

图1 分析模型

 

图2 数据文件

 

2.       无渗流模式 — 数据文件SoilLayerHeave-NoFluid.dat

以下给出无渗流模式分析代码，并对操作环节做突出标识。

 

model new

model large-strain off

model title 'Raising the water table - not in configure fluid'

fish automatic-create off

call 'fishFunctions'

[setup]

zone create brick point 0 (0,0,0) point 1 (2,0,0) point 2 (0,1,0) point 3 (0,0,10) size 2 1 10 ratio 1 1 1

zone cmodel assign elastic

zone property bulk [m_bu] sh [m_sh]

; --- column is dry ---

zone property density [m_d]    ; m_d为土层干密度1800kg/m3

; --- boundary conditions ---

zone face apply velocity-z 0 range position-z 0

zone face apply velocity-x 0 range union position-x 0 position-x 2

zone face apply velocity-y 0 range union position-y 0 position-y 1

; --- gravity ---

model gravity 0 0 [mgrav]

; --- histories ---

history interval 10

zone history displacement-z position (0,0,10)

fish history ana_dis

; --- initial equilibrium ---

zone initialize stress-zz -1.8e5 grad 0 0 1.8e4

zone initialize stress-xx -1.029e5 grad 0 0 1.029e4

zone initialize stress-yy -1.029e5 grad 0 0 1.029e4

model solve ratio 1e-6

model save 'ini'

; --- raise water level ---

; (note: when not in configure FLUID, water density is assigned using:)

zone water density [w_d]  ; w_d为流体密度1000kg/m3

; (we can do it this way ...)

zone gridpoint initialize pore-pressure 1e5 grad 0 0 -1e4

; (or this way ...)

;water table origin 0 0 _H normal 0 0 -1

; (total stress adjustement)

zone initialize stress-xx -1e5 add grad 0 0 1e4

zone initialize stress-yy -1e5 add grad 0 0 1e4

zone initialize stress-zz -1e5 add grad 0 0 1e4

; --- use wet density below water table ---

zone property density [m_rho]  ; m_rho为土层饱和密度2300kg/m3

; --- static equilibrium ---

model solve ratio 1e-6

model save 'ncfl'

 

其中，主要参数：

l     m_bu、m_sh分别土层的体积模量与剪切模量；

l    m_d、w_d分别为土层干密度1800kg/m3和流体密度1000kg/m3；

l    m_rho为土层饱和密度=2300kg/m3，按 式    求得。其中，为饱和度，在本例中取1。

 

参考以上代码，可总结主要步骤：

1)  在初始自重平衡过程中，采用zone initialize stress命令对模型应力作初始化操作，模型底部竖向应力为       = 1800 * 10 * 10 = -1.8e5Pa，在竖向按梯度1.8e4Pa变化分布；考虑侧压力系数为0.5714，得到模型底部水平向应力       =-1.029e5Pa，竖向分布梯度为0.5714 * 1.8e4Pa = 1.029e4；进而施加重力，并迭代至平衡、得到水位变化前初始状态；

2)  利用zone gridpoint initialize pore-pressure或water table命令定义对应于水位上升至土层顶部的孔压分布；同时采用命令zone property density将土层密度改变为饱和密度；

3)  执行迭代计算，得到因水位变化引致的土层变形分布特征，分析得到土层顶部竖向变形为0.00178546m，与理论解一致。

 

3.       渗流模式A — 数据文件SoilLayerHeave-withFluidA.dat

以下给出无渗流模式分析代码，并对操作环节做突出标识。

 

model new

model title 'Raising the water table - configure fluid'

model large-strain off

fish automatic-create off

model configure fluid

call 'fishFunctions'

[setup]

zone create brick point 0 (0,0,0) point 1 (2,0,0) point 2 (0,1,0) point 3 (0,0,10) size 2 1 10 ratio 1 1 1

zone cmodel assign elastic

zone fluid cmodel assign isotropic

zone property bulk [m_bu] sh [m_sh]

; --- column is dry ---

; (initialize sat at 0)

zone gridpoint initialize saturation 0

zone property density [m_d]

; --- boundary conditions ---

zone face apply velocity-z 0 range position-z 0

zone face apply velocity-x 0 range union position-x 0 position-x 2

zone face apply velocity-y 0 range union position-y 0 position-y 1

; --- gravity ---

model gravity 0 0 [mgrav]

; --- histories ---

history interval 10

zone history displacement-z position (0,0,10)

fish history ana_dis

; --- initial equilibrium ---

zone initialize stress-zz -1.8e5 grad 0 0 1.8e4

zone initialize stress-xx -1.029e5 grad 0 0 1.029e4

zone initialize stress-yy -1.029e5 grad 0 0 1.029e4

;

model fluid active off

model mechanical active on

zone gridpoint initialize fluid-modulus 0

model solve

model save 'ini'

; --- raise water level ---

; (initialize sat at 1 below the water level)

zone gridpoint initialize saturation 1

; (note: in configure FLUID, water density is assigned using:)

zone initialize fluid-density [w_d]

; (we can do it this way ...)

zone gridpoint initialize pore-pressure 1e5 grad 0 0 -1e4

; (or this way ...)

; (note: we can use water table command in configure fluid,

; this is different from FLAC)

;water table origin 0 0 _H normal 0 0 -1

; (total stress adjustement)

zone initialize stress-xx -1e5 add grad 0 0 1e4

zone initialize stress-yy -1e5 add grad 0 0 1e4

zone initialize stress-zz -1e5 add grad 0 0 1e4

; --- Note: no need to specify wet density below water table ---

; --- static equilibrium ---

model fluid active off mech on

zone gridpoint initialize fluid-modulus 0

model solve

model save 'cflA'

 

依据上述代码，可知该渗流分析模式主要环节：

1)   通过模块配置命令model configure fluid调用流体计算功能模块；

2)  土层在自重条件下的初始应力建立过程与无渗流分析模式基本一致。注意，在自重平衡迭代时，需通过命令model fluid active off关闭流动计算，同时设定流体体积模量fluid-modulus为0，以避免因土层受自重压缩作用而形成超孔压；

3)  利用zone gridpoint initialize pore-pressure或water table命令定义对应于水位上升至土层顶部的孔压分布条件；

4)  继续关闭流体计算，仅执行力学迭代得到土层因水位变化导致的变形分布特征。分析得到土层顶部竖向变形为0.00178365m，与理论解一致。

 

可见，较之无渗流模式，在渗流模式中，我们不需要人为调整模型中岩土材料的密度。特别例如饱和密度，由程序依据流体密度、孔隙率和饱和度按式 自动计算得到。代码明确了对流体密度和饱和度的设定，但未见对孔隙率的定义，这是因为程序默认材料的孔隙率为0.5，与本次分析的土层条件一致。

 

此外，由于分析条件简单，本例中孔压分布条件由理论解直接利用zone gridpoint initialize pore-pressure或water table命令定义，或者说，水位抬升后孔压分布条件并未实际利用渗流计算得到，因此是在简单分析条件下的一种简化处理方法。若工程地质与水文地质条件复杂，需开展渗流计算以确定对应于给定条件的孔压分布状态。

 

4.   渗流模式B — 数据文件SoilLayerHeave-withFluidB.dat

与渗流模式A的区别在于，本例采用渗流计算得到水位抬升在土层中形成的孔压分布，同时考虑因此导致的土层变形。以下给出无渗流模式分析代码，并对操作环节做突出标识。

 

model new

model title 'Raising the water table - configure fluid'

model large-strain off

fish automatic-create off

model configure fluid

call 'fishFunctions'

[setup]

zone create brick point 0 (0,0,0) point 1 (2,0,0) point 2 (0,1,0) point 3 (0,0,10) size 2 2 10 ratio 1 1 1

zone cmodel assign elastic

zone fluid cmodel assign isotropic

zone property bulk [m_bu] sh [m_sh]

; --- column is dry ---

; (initialize sat at 0)

zone gridpoint initialize saturation 0

zone property density [m_d]

; --- boundary conditions ---

zone face apply velocity-z 0 range position-z 0

zone face apply velocity-x 0 range union position-x 0 position-x 2

zone face apply velocity-y 0 range union position-y 0 position-y 1

; --- gravity ---

model gravity 0 0 [mgrav]

; --- histories ---

history interval 10

zone history displacement-z position (0,0,10)

fish history ana_dis

; --- initial equilibrium ---

zone initialize stress-zz -1.8e5 grad 0 0 1.8e4

zone initialize stress-xx -1.029e5 grad 0 0 1.029e4

zone initialize stress-yy -1.029e5 grad 0 0 1.029e4

;

model fluid active off

model mechanical active on

zone gridpoint initialize fluid-modulus 0

model solve

model save 'ini'

; --- raise water level ---

; (initialize sat at 1 below the water level)

 

zone fluid cmodel assign isotropic

zone fluid prop perm 1.e-12 poros 0.5

zone gridpoint initialize saturation 1

; (note: in configure FLUID, water density is assigned using:)

zone initialize fluid-density [w_d]

; (we can do it this way ...)

zone gridpoint fix pore-pressure 1e5 grad 0 0 -1e4 range position-x -0.1 0.1

zone gridpoint fix pore-pressure 1e5 grad 0 0 -1e4 range position-x 1.9 2.1

;--------------------

model fluid active on mech on

zone gridpoint initialize fluid-modulus 2.e8

model solve

model save 'cflB'

 

可见，较之渗流模式A，本例明确指定了土层的渗透率，同时利用zone gridpoint fix pore-pressure完整设定了与水位抬升条件一致的渗流边界条件，继而开展应力-渗流完全耦合分析，计算得到土层顶部竖向变形为0.00179108m，与理论解基本一致。图3为渗流分析孔压边界条件以及最终孔压分布状态。

 

图3 孔压边界条件及最终孔压分布

 

三、  讨  论

上一小节通过一简单案例分析，介绍了利用不同分析模式（或方法）开展流体作用模拟的实现过程。对比不同方法的特点，在某些具体环节或许还存在进一步讨论的价值。

 

1.     问题的力学原理

以土层底部作为考察对象，分析其在水位抬升前后的应力变化特征，特别是决定变形行为的有效应力状态。在水位抬升前，有效应力完全由其干密度决定，可知土层底部有效应力1800*10*10=1.8e5Pa，当水位抬升至土层顶部时，底部有效应力取决于土层浮密度，即(2300 – 1000) * 10 *10=1.3e5Pa。可见，水位抬升将导致土层中有效应力降低，土层因此也会产生沿铅直向上的抬升变形。

 

2.     应力初始化及对模型响应行为的影响

类似zone initialize stress-xx等命令可实现对模型单元的应力条件进行初始化操作。需说明的是，由FLAC3D程序采用增量迭代原理决定的，程序无法给出初始化操作定义的那部分应力可以导致的模型变形量。或者说，增量计算只能获得模型总体受力与初始化应力水平两者间的不平衡力形成的变形。

 

一个极端情形是，如果初始化应力水平与模型所承受的荷载（外荷载与体力）能够平衡，则在模型中形成的变形接近为0。在上一小节案例中，在模型自重平衡阶段均采用如下应力初始化代码对模型应力进行先期赋值：

zone initialize stress-zz -1.8e5 grad 0 0 1.8e4

zone initialize stress-xx -1.029e5 grad 0 0 1.029e4

zone initialize stress-yy -1.029e5 grad 0 0 1.029e4

由于该初始化应力水平能够平衡模型所承受的荷载（仅自重），因此程序在自重模拟阶段可迅速收敛，且模型中的变形接近为0。

 

以上原理可简单理解为岩土体变形响应与应力路径具有相关性。

 

3.    孔压初始化与总应力之间的关系

ITASCA系列软件均依据有效应力原理进行固体介质变形受力分析，单元总应力由有效应力与孔压叠加形成。

 

需特别强调的是，类似zone gridpoint initialize pore-pressure或water table等孔压初始化命令仅起到调整单元孔压的作用，而不会改变单元的总应力水平。以上一节案例渗流模式A为例，在模拟水位抬升作用时，采用如下命令对模型进行孔压初始化赋值：

zone gridpoint initialize pore-pressure 1e5 grad 0 0 -1e4

该操作并不会改变总应力水平，其作用主要在于调整了总应力中孔压、有效应力的占比，如图4所示。

 

图4 孔压初始化前后模型中的应力构成变化

 

进一步考察图4，孔压初始化结果导致土层有效应力水平降低，以土层底部为参考，有效应力自1.8e5Pa降低为0.8e5Pa。或者说，在开展水位抬升引致土层变形计算前，土层底部的实际有效应力为0.8e5Pa，这显然不符合实际条件。具体而言，在水位抬升过程中，土层底部有效应力实际自1.8e5Pa降低为1.3e5Pa。依据前述增量计算方法中岩土体变形与应力路径的相关性讨论可知，有必要调整土层有效应力与实际一致，这就是在上一节无渗流模式和渗流模式分析中，除孔压初始化外，还同时采用如下命令调整模型应力场的内在原因：

zone initialize stress-xx -1e5 add grad 0 0 1e4

zone initialize stress-yy -1e5 add grad 0 0 1e4

zone initialize stress-zz -1e5 add grad 0 0 1e4

 

图5 有效应力调整对计算结果的影响

 

图5比较有效应力调整与否对水位抬升引致土层变形结果的影响。可见，在不通过人为方式干预有效应力（及不调整有效应力）的前提下，计算获得的竖向变形量为0.00178412m，且变形性质为沉降变形，这一结果显然是错误的。

 

另需提及，是否需要采用人为方式干预有效应力是有前提的。首先，本分析在于关注水位抬升引致土层变形，且无渗流模式和渗流模式A直接对水位抬升最终孔压分布进行了初始化，未真实模拟孔压演变以及相应的有效应力的动态变化过程，此为需进行人为干预的本质原因；而渗流模式B则较为真实模拟了模型中应力的变化过程，因此不需要进一步采用认为干预；另外，在无渗流模式和渗流模式A中，若只关心获得水位抬升后模型中最终应力状态，由于模型荷载条件（此处仅为重力）和孔压分布条件唯一，最终两种模式计算得到的有效应力也是一致的。

 

小  结

l  ITASCA等效连续介质流体分析技术提供的分析模式较为灵活，视问题性质选择适合的分析模式；

l 无论是无渗流模式还是渗流模式，zone gridpoint initialize pore-pressure或water table命令均可用于对模型中的孔压分布进行初始化操作；

l 在工程分析中，若关注水文地质条件变化引致岩土体变形响应分析，需依据分析模式的不同考虑是否对模型中应力状态进行必要的人为干预。总体而言，在孔压分布为初始化定义的情形下，由于该方式未真实动态地模拟模型中有效应力的变化过程，鉴于变形响应与应力路径具有相关性，因此需在变形计算前人为干预调整模型中应力分布。

 

参考文献：

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