(完整word版)FLAC3D流力耦合作用

1、FLAC 3D流力耦合作用 1. 1 耦合作用简介 . 1 1. 2 数学模型描述 . 2 1.2.1 规定和定义 . 2 1.2.2 流体重量平衡方程 . 3 1.2.3 流动法则 . 4 1.2.4 力学结构法则 . 4 1.2.5 边界及初始条件 . 5 1. 3 数值公式 . 5 1.3.1 空间导数的有限差分近似 . 5 1.3.2 质量平衡方程的节点公式 . 6 1.3.3 显式有限差分公式 . 8 1.3.3.1 稳定标准 . 9 1.3.4 隐式有限差分公式 . 9 1.3.4.1 收敛准则 . 11 1.3.5 力学时间步和力学稳定性 . 12 1.3.6 总应力修正 . 1

2、2 1. 4 流动耦合问题的属性和单位 . 12 1.4.1 渗透系数 . 13 1.4.2 Biot 系数：和 Biot 模数 M . 13 1.4.3 流体体积模量 . 14 1.4.4 孔隙率 . 14 1.4.5 密度 . 14 1.4.6 流体张力限 . 15 1. 5 单一流动问题和耦合流动问题 . 15 1.5.1 恒定孔压（用于有效应力计算） . 15 1.5.2 建立了孔压分配的单一流动计算 . 16 1.5.3 非流动，力学变形产生的孔隙压力 . 16 1.5.4 耦合流动和力学计算 . 17 1. 6 对于渗流分析的输入指导 . 18 1.6.1 FLAC3D 命令 .

3、18 1.6.2 FISH 变量 . 21 1.7 验证举例 . 22 1.7.1 在限制层内的不稳定地下水流动 . 22 1.7.2 单方向固结 . 25 1.7.3 穿透浅含水层限制边界的井水流动 . 29 1.1耦合作用简介 FLAC3D允许在饱和多孔材料中进行流体流动的瞬时模拟。 流动计算可以脱 离FLAC3D中的力学计算独立进行，也可以与其他力学模型进行耦合计算，以控 制流一一固耦合作用的影响，其计算具有如下特征。 1. 提供了在各向同性条件下的流体运动法则，也提供了在流动区域中的 无渗流材料的流动零模型。 2. 不同的区域可以有不同的流动模型和法则。 3. 流体压力、流量通量、松散

4、和不透水边界可以设定。 4. 流体源可以以点源和体积源插入在材料中。这些流源可以是随时间变 化的流入源和流出源。 5. 对于流体计算既提供了显式算法也提供了隐式算法模式。 6. 任何力学和热学模型都可以和流体模型进行耦合计算。在耦合计算中， 流体和固体骨架的压缩和热膨胀都是允许发生的。 7. 由于骨架变形产生的流一一力耦合效果通过 Biot 固结系数来实现。 8. 流体与热之间的耦合作用可以通过热力学线性膨胀系数 :t 和不排水 热力学系数一：来实现。 9. 流动法则是基于线性理论的，其假设恒定的材料属性并忽略对流。流 体与固体骨架的温度保持局部平衡。非线性变形则要通过 FISH 函数 设定孔

5、压和材料的属性来实现。 1.2数学模型描述 在FLAC3D中的变形传播增量的准确描述过程在 Biot半静态线性理论框架内 进行。下面详细介绍水-热-力耦合与热-孔压-力的方程。这些方程可以被被应用 于孔隙介质中的达西流动定律。各种各样的流体-包括空气和水-都可以应用本模 型。包含在描述经过孔隙介质的流动的变量是孔隙压力和比流量矢量的三个组成 部分。这些变量是通过流体运动的重量平衡方程-达西定律和描述流动对于孔隙 压力、体积应变和温度的改变的反应机制的方程联系在一起的。 孔隙压力和温度 的影响被包含在完全水热力耦合的力学结构法则增量中。 假设体应变率和温度率已知，把重量平衡方程-达西定律-代入流

6、动结构关系 中，根据孔隙压力产生一个不同的方程，此方程可以根据具体的几何特性、属性、 边界和初始条件解出。 1.2.1 规定和定义 作为一个习惯符号，符号ai表示笛卡儿坐标系中矢量a的 i 部分；Aj表示 张量A的（i, j ）部分。而f,i用来代表与人相联系的f的偏导数部分（符号f表 示变量的矢量或标量部分）。 爱因斯坦的求和方程只应用于有下标i,j,k的变量，下标的取值分别为 1,2, 3 代表空间组成。当用于矩阵方程时，下标可以取任何值。我们用国际单位制来 描述变量参数。 F面这个无量纲参数来描述瞬时流动特性，特性长度 Lc: 流体扩散率: (1-2) 这里 k 是渗透率，M是 Biot

7、 模数，是 Biot 系数，：1=KG4，K 是体 积模量，而 G 是排水弹性材料的剪切模量。 Biot 系数考虑孔隙材料土骨架的压缩性。如果:取单位 1,那么土骨架被认 为是无压状态的，而 Biot 模数 M 等于kw/，这里kw是流体的体积模量，而n是 孔隙压力。流体的扩散率成为下面的式子： (1-3) 注意，对于单一流动不考虑耦合效应的计算，流体扩散率成为下面的式子： c = kK （ 1-4） 特性时间： (1-5) 1.2.2 流体重量平衡方程 对于小变形的问题，流体的重量平衡方程表述为： 一4叩 qv （ 1-6） ct 这里qi为比流量矢量单位是m/s, qv是体积流量源强度，单

8、位是1 sec】，而 是 每单位体积孔隙材料的流体体积或流体容量。 在FLAC3D中，流体体积变量的改变与孔隙压力 p、力学体积应变：和温度 Lc Vf af (1-1) 这里Vf表示流动区域的体积, a f表示流动区域的表面积。 tc T 的改变成线性关系。流体结构法则表述为： 1 P T -:t M ：:t ：:t ：:t 这里M是 Biot 模数，单位是N/m2】，用是 Biot 系数，是不排水热力系 数，单位1 C,他们都考虑土骨架和流体的热力膨胀因素。 把（1-6）代入（1-7） 得到： -qi,i qv （ 1-8） M 盘 这里 qv = qv T （ 1-9） ct 盘 1.2

9、.3 流动法则 基本的法则是定义比流量矢量与孔隙压力之间关系的达西定律。对于均质， 各向同性土体和恒定流密度条件下，这个法则以下列形式给出： qi kb - ?fXjgj 1 （ 1-10） 这里 k 是渗透系数，单位是m4Ns，陥是流体密度，单位kg/m3，gi， i =1,3是重力加速度的三个分量，单位是m s2 。 为了供将来参考，量： （ 1-11） 被定义为重力头，而 ?f-定义为压力头。 1.2.4 力学结构法则 关于孔隙压力的体应变的影响在流体力学法则中反映出来。 依次地，孔隙压 力的改变引起力学变形的发生：孔隙固体结构方程的增量形式表述如下形式： q ；巾 r =Hij*ij“

10、 - j） （ 1-12） 这里是关联循环应力增量， 是给定的函数，I是总应变，“ 1 是热 (1-7) 应变，而 r 是克罗内克尔增量。 尤其，弹性关系可以用下式表示： 匚 ij i0 P - Po ； ij 二 2G ；ij - * 2 ；kk - ；kk （ 1-13） 这里：2 = K 2G 3，：二0和p0代表初始状态 125 边界及初始条件 用（1-10）代替（1-8）中的qi，产生不同的流动方程，假设 qv已知。初始 条件与给定的压力区一致。边界条件通常根据孔隙压力或垂直与边界的比流量矢 量来表述。在 FLAC3D中，定义了四个类型的条件：（1）给定孔隙压力；（2）给 定垂直于边

11、界的比流量矢量；（3）渗透边界；（4）不渗透边界。 在FLAC3D中，渗透边界有如下形式： (1-14) 这里qn是垂直于边界的比流量矢量的外垂直线方向的部分， h 是渗透系数，单位 m3/N -sec， p 是在边界表面的孔隙压力，而 pe是在渗透层的孔压。 默认情况下，边界条件是不透水边界。 1.3数值公式 在 FLAC3D中，重力平衡方程式（1-8）与达西定律（1-10），对于给定的q： 两式组合在一起用基于介质离散化的有限差分方法解出， 离散化区域由两层恒 定的比流量四面体覆盖层组成。数值计算是根据重力平衡方程的节点公式进行 的。这些公式即为牛顿力法的节点公式。它是通过用孔隙压力、比流

12、量矢量和孔 隙压力梯度分别代替速度矢量、应力张量和应变率张量来得到。初等常微分方程 的解系是用与时间一致的显式和隐式离散化两种模式解出。其主解描述如下。 1.3.1 空间导数的有限差分近似 根据习惯，四面体节点的指局部由 1 到 4 编号的节点。面n与编号为n的节 点相对。上标f与面f上的独立变量的值相联系。线性孔隙压力变量和恒定流体 密度被假设在四面体内部；根据应用高斯离散定理得到的节点孔压值得到压力头 梯度的表达方式： 4 （p PfXjgj ）j = - Z （p1 Pf X： gi nj1 S（） （ 1-15） 3V -1 这里n 1是垂直于面 I 外单位矢量，S 是面的表面积，V

13、是四面体的体积。 对于数值精确性，量Xi x1，这里x1与四面体顶点之一的坐标一致，在（1-15） 的压力头公式中可以用Xi来代替，而且并不影响压力头梯度的值，根据（1-15）， 有下面公式： P - ?fXigi ,j 二 4 -p 1 nj S 3V j （1-16） 这里节点量p 定义如下： p 轴=p _Pf（x： -x1 gi （1-17） 1.3.2 质量平衡方程的节点公式 质量平衡方程（1-8）可以被描述为： qi,i b =0 （1-18） 这里 屮=丄生 （1-19） M 它等于瞬时“体力”，B 用在力学节点公式。应用这种相似，节点流量, Q： m3 s, n =1,4,等于

14、四面体的比流量和体积源强度 b.，可以表述为： mn 竝 4 dt 这里 4M n 原则上， 在每个节点的重力平衡方程节点形式都是根据需要建立的, 边界的等节点流量（-Qe）与节点贡献流量（QW ）的总和为零。 式（1-21）的比流量矢量部分与根据流动法则（1-10）得到的压力头梯度相 联系。依次地，压力头梯度部分可以根据式（1-16）四面体节点孔压来表述。 为了节省计算时间，在用式（1-21）在四边形每个节点n遇到的网格区获得 的网（1-20 Qn qi nSn （1-21 （1-22 作用于 格区贡献流量Q；时，采用网格区公式。局部网格区的矩阵集合把节点值 Qn与节点压力头p n联系起来。

15、因为这些矩阵是对称的，36 个组成部分被计算；这 些计算在大变形模式下每十个计算步更新一次。根据网格区矩阵的定义，我们得 到： Qn 二 MnjP j （ 1-23） 这里 P”是对于网格区的节点孔隙压力头的局部矢量。 依次地，总节点值 Q；通过网格区贡献的累加来获得。我们可以写出每个节 点的总流量公式： Q； =CnjP j （ 1-24） 这里 C 1 是总矩阵，P 1 节点压力头的总矢量。 返回到我们前面的考虑因素，我得到： 八 Qen QW =0 （ 1-25） 这里，为了符号的简单化，符号 1 用来代表在节点n与该节点相遇的所有 网格区的贡献流量之和。（一个网格区的贡献由网格区内所有

16、的四面体流量之 和）。应用（1-20）,（ 1-9），（ 1-25），我们得到: 这里 Q；是用式（1-24）和（1-17）得到节点孔隙压力的函数，7 ；pp是体积源、 边界通量和点源贡献的已知流量。有如下形式： n Q；pp qvV Qw （ 1-27） -4 而送 Qtnm是节点热力贡献流量。定义如下: (1-28) 方程（1-26）是节点n的重力平衡方程的节点形式；右边 dpn dt Q; Q：PP “ Q；hJ (1-26) Q； +Qapp+Qt：m 是指不 平衡流量。它由两部分组成：流体不平衡流量 Q； ：pp ;不平衡热力流量 Qthm。当流体保持静止时（流体计算停止关闭），不平

17、衡流消失，孔隙压力的 改变只由力学或热变形引起。对于独立的流体计算，不平衡热力流量等于零。 在 FLAC3D中，Biot 模数是节点属性，利用（1-22）我们得到： 使每个节点包含在离散化状态。这些方程在 FLAC 3D 中的解形成一个初等常微分方程的解系， 对于给定的 血，应用显式或隐式有限 dt 差分方法解出。 1.3.3 显式有限差分公式 在显式有限差分公式中，在节点的量P - Pv的值在一个时间段 t内假设承线 性变化。原则上，在式（1-31）左边的导数应用有限差分方法解出，在时间 t 的 流体不平衡流量也被估算。由初始孔隙压力区开始，在增加时间内的节点孔隙压 力的值被更新，假设孔隙压

18、力是变化的，有下式： P；沁=P； P：t. 二P； （1-33） 这里 二 “ QT 七 Q；pp J （ 1-34） 如果时间步保持在一个有限值以下，显式差分的数值稳定性才能保持。这里 vn八 （1-28）代入（1-26）,我们写出: 这里 n Pv Vn n - Q；ppl - Z V、 - -心卩忌*10，（口詞也卩思） （1-58） 时间步的大小必须在收敛和隐式解的精确性之间进行选择。尽管 Cran k-Nicolson 方法对于所有正值的 t是稳定的，Jacobi 迭代收敛方法是无条件保证的，除非矩 阵A是完全主对角矩阵一一例如，当： n AjjAkk （ 1-59） j衣 对于

19、1 汨冬 n。根据式（1-52）定义的Anj，和由（1-60）定义的n，足够多的条 件对于足够小的 t 总是满足的。如果没有达到 Jacobi 迭代收敛条件，就会出现 错误信息。那就必须减小隐式时间步的大小或者用显式的方法。 尽管隐式方法也具有二阶精确性，在选择合适的时间步时要十分小心。 实际 上,它的值与孔隙压力的波动波长相比必须保持足够小。 有代表性的是，显式差 分通常是更早地应用于计算过程或在它的扰动阶段。 （隐式差分方法当超精确度 时可以应用显式时间步的值），然而隐式方法更适合于模拟残余。(1-55) (1-56) 135 力学时间步和力学稳定性 流体的存在将明显增加介质的力学体积模量

20、， 依次地，也将影响在力学数值 稳定性方面节点重力密度的大小。通过考虑不排水孔隙 -力条件（一：，：T ）下的 等温线，可以发现割线模量的上限值。那么流体结构法则（ 1-7）的增量形式如 下: 二 p - - ： M 二 2 （ 1-60） 采用弹性法则（1-13）的增量形式来描述在时间段内的应力-应变关系： 1 H ip 二 K （ 1-61） 3 把式（1-60）代入（1-61）得到: 1 2 H = K 叱 ML】 （1-62） 3 1.3.6 总应力修正 在FLAC3D公式中，节点孔隙压力首先被计算，网格区孔隙压力通过四面体 的体积平均值得到，在本文中，总应力修正由两部分组成：一部分是

21、对于流体的 2; 部分是针对于热-力耦合的也 0ith。可以写成如下形式： 二：j 二:C j * 二 ij （1-63） -ii - -Pnt / ij （1-64） =aM -pAT).-ij （1-65） 画线符号表示上面提到的网格区的平均值。 1.4流动耦合问题的属性和单位 在 FLAC3D中，与流动相联系的属性是对于可压缩流体介质渗透系数 k，流 体重力密度 4，Biot 系数，Biot 模数M，流体体积模量Kf和孔隙率n。热耦 合参数排水线性热膨胀系数:t和不排水热系数：。 所有的热-孔-力数量必须以同意的单位给出，系统不提供转化单位的功能。 1.4.1 渗透系数 各向同性渗透系数

22、 k （在国际单位制中的单位是 m2 / Pa - sec ）。它是达西定 律中的压力系数并与水压传导率kh m s相联系： khkgf （ 1-66） 这里 g 是重力加速度。注意，“本质渗透系数”被定义为渗透系数的和动力粘 滞系数u的产生物。下列的数值转换可以被用于得到国际单位制的 k : k 三 * （单位 cm2） 9.9 10， k 三 kh （单位 cm sec） 1.02 10 k=k （单位毫达西） 9.8 10 这里 u =1.011 103kg/m-sec， ：?w = 1000kg / m3， g = 9.81m/sec2。 这个渗透系数是用 PROPERTY 命令定义的

23、网格区属性。 1.4.2 Biot 系数口和 Biot 模数 M Biot 系数 a 定义为一个材料单元中的流体体积与当孔隙压力返回到它的初 始状态时的体积改变值的比率。它同样可以在获取材料体积模量 K 的排水实验来 获取。它的变化范围介于和 1 之间。n 是孔隙比。在不可压缩骨架的组成 2n 十 3 的特殊情况下，=1。此值是默认值。 对于一个理想的孔隙材料，Biot 系数与材料固体部分的体积模量Ks相联 系： K :=1 （ 1-67） Ks Biot 模数 M 定义为： K -K M 二 u 2 （ 1-68） 这里Ku是材料的不排水体积模量。 对于理想的孔隙材料，Biot 模数 M 与

24、流体体积模量K f相联系: Kf M - （1-69） n +（a n【1 a K f / K 这里 n 是孔隙比。因此，对于一个不可压缩的土骨架的情况下（:.=1）： (1-70) 对于可压缩土骨架的计算模式用 SET FLUID BIOT ON 。这里 Biot 系数: 是用PROPERTY 命令定义的网格区参数。Biot 模数 M 是用 INITIAL 命令定义 的网格点变量。 1.4.3 流体体积模量 对于包含不可压缩土骨架的流动计算，流体体积模量定义如下： p Kf （ 1-71） AV/V 这里中是对于体应变 V/V 的压力改变值。 流体的“压缩性系数” Cf与Kf是相互影响的，K

25、f =Cf，例如，对于室 温下的纯净水 Kf = 2 109Pa （国际单位）。在真实土体中，孔隙水可能包含一些 溶解的空气和空气泡，这充分的减少了流体的割线体积模量。 经典土力学教材提 供了这个减少值的大小。对于地下水问题，水的体积模量在不同的网格区可能不 同。FISH 函数可以根据一些原则来来改变局部模数。但是，使用者应该避免使 用不合理的假设。 当流体的体积模量Kf以输入的方式给出时，对于不可压缩骨架的 Biot 模数 M 用式（1-70）在内部进行计算。在这个计算中，网格区参数 -孔隙率在节点处 用节点体积平均值来估算。然后 Biot 系数在整个流动区域内设为 1,并与任何 给定的参数

26、值无关。 1.4.4 孔隙率 孔隙率n是取值由 0 到 1 的一个无量纲的数，定义为孔隙体积与单元总体积 的比值。默认值是 0.5。此参数是 FLAC3D用来计算介质的饱和密度，还有就是 当流体的体积模量 Kf已知的情况下用来估算 Biot 模数 M。它是一个用 PROPERTY 命令定义的网格区参数。 1.4.5 密度 在不同的环境下可以给出三种不同的密度，固体材料的干密度 6，饱和密 度二，及流体密度几 只有在计算有效应力时才应用到饱和密度（不是在 CONFIG FLUID 模式）, 用 WATER TABLE 命令（或者 INITIAL PP 命令）定义地下水 TABLE。干密度 定义含

27、水面以上的土体密度，饱和密度定义含水面以上土体密度。 对于一个流体计算（CONFIG FLUID 模式），不管什么时候，重力矢量部分 都是用 SET GRAVITY 命令初始化定义的，而干密度和流体密度也必须给出。 然后饱和密度在 FLAC3D内部用公式亠=4 来进行计算。 饱和或干密度用INITIAL DENSITY命令给出， 流体密度可以用 WATER DENSITY命令进行整体施加， 或者它允许用 INITIAL FDENSITY 命令使其随 位置的改变而改变。所有的密度都是网格区变量。 1.4.6 流体张力限 在细粘土中，孔隙水能保持足够的张力，在 FLAC3D中，默认情况下，负孔 隙

28、压力将会发展。 注意，负孔隙压力并不是所谓的张力是由于毛细作用或电化学力产生的， 后 者的力最好用结构模型内部的有效应力增量来代替。 负的流体压力与土骨架的组 成成分无关，它只与用流体填充后的体积膨胀率有关。 用户可以用 SET FLIUD PCUT ON 命令来把流体张力限设置为零，当流体 压力也回落为零时，流体内开始有气体和气泡，逻辑性组织网格点的孔压变为负 值。并协调总应力和网格区孔压。 1.5单一流动问题和耦合流动问题 FLAC3D具有单独进行流动分析和流动耦合分析的功能，耦合分析可以和任 何内置的力学模型相耦合。然而，需要注意的是，流动模型不是内部自动置零， 而力学模型确是内部自动置

29、零。在非流动区域必须用 MODEL FL-NULL 命令进 行定义。 耦合程序可以分为几个阶段。应用 WATER TABLE 命令就是其中之一，而 且在计算过程中并不需要额外的分配内存。而要用 CONFIG FLUID 来设定网格 流动计算，用 MODEL FL-ISO 命令使所有的网格区发生流动。 对于耦合分析的不同阶段将在下面详细介绍。 1.5.1 恒定孔压（用于有效应力计算） 在固定的潜水层表面以下，用 WATER TABLE 或 INITIAL PP 命令来定义初 始流体静力学孔压分配。（可供选择的是，可以利用 FISH 函数来进行流体静力 学孔压分配）在这种情况下，对于流动网格不需要

30、形成。水的密度必须提供，饱 和密度，干密度也必须给定。 对于没有应变发生的区域， 孔压的分配与初始状态保持一致， 它保持恒定并 且 不受力学变形的影响。流动没有发生。 对于依靠主应力产生孔压分配的材料对于孔压分配的影响是失效的， 这个模 拟方法对于估计边坡稳定性是十分有用的，为了计算的目的，含水（ WATER TABLE ）状况设定并给出。 1.5.2 建立了孔压分配的单一流动计算 单一的流动计算来确定在一些独立系统的流动和压力分配。 例如，估算由于 排水渠的开挖或排水井的作用导致的地下水的改变可能是必须的。 另外，对于耦 合计算来说， 可能需要进行初始孔隙压力的分配， 但是在区域边界只提供水

31、压条 件。在两种情况下，FLAC3D在没有任何力学计算的情况下只在流动状态下进行。 接下来， 力学计算可以继续做也可以不做。 然而值得注意的是， 在这里完全饱和 流动是唯一的流动模式。 尤其，施加的水压边界条件必须与随流动计算无影响的 地下水表面一致。 在命令程序的第一步是键入 CONFIG FLUID 命令，以致有更多的内存分配 给流动计算。 用 SET MECH OFF 命令停止力学计算。 然后必须在隐式和显式差 分计算之间进行选择。 默认情况下， 选择显式有限差分模式， 但是隐式差分模式 可以通过 SET FLUID IMPLIcIT ON 或 SET FLUID IMPLIcIT OF

32、F 命令开关。 在显式模式下，流动时间步自动计算，但是可以用 SET FLUID DT 命令来 选择更小的流动时间步。 但是在隐式模式下， 时间步的大小必须由用户自己来定 义。也由 SET FLUID DT 设定。当对比渗透性而言，用隐式方法更为有效。 流动模型和流动属性必须分配给所有的网格区（ ZONE），初始和边界条件 被分配给完全的流动问题。在单一流动或流 -力耦合问题中的流动区域用非零流 动模型定义。例如，流动边界条件可以用 APLLY 和 RANGE 组合命令来分配。 （记住，对于零力学区域其流动模型并不自动置零） 。 可以用 STEP 命令来执行给定的流动时间步，当一个特定的流动时

33、间达到时， 为了自动阻止流动计算，可以用 SOLVE AGE 命令来进行设定。另外， SET FLUID AGE和 SET FLUID STEP 可以用来设定流动时间限或最大运行步。然后 给出 SOLVE 命令。 如果在力学计算中需要进行孔隙压力的分配，那么假设孔隙压力保持恒定， 而SET FLUID OFF MECH ON 命令应该被给出。而 Biot 系数也应该被设定为零， 以阻止由于力学变形引起的超孔隙压力的产生。 1.5.3 非流动，力学变形产生的孔隙压力 一个耦合系统的短期反应机制在 FLAC 3D中也可以在限制流动下分析，但是允 许力学计算产生孔隙压力。如果给出 SET FLUID

34、 OFF 命令， Biot 模数或流体模 数给出一个真实值， 那孔隙压力将作为一个力学变形的结果而产生。 例如， 由于 底角荷载产生的瞬时孔隙压力可以用这种方法进行计算。 如果流体的体积模量比 固体的体积模量更大， 计算程序收敛的速度会很慢。 逐渐的减少流体的体积模量 而不影响计算结果是可以实现的。 在例.1 中，描述包含在一个盒子中的弹塑性材 料在地窖荷载的作用下空隙压力建立的过程。 盒子的左边界是线性对称的， 而沿 上表面的孔隙压力设定为零以阻止孔隙压力的产生。 默认情况下，孔隙率为 0.5, 渗透性不需要所以没有流动发生。 在加载期间有一个大数量的塑性流动发生， 法向应力用 FISH 函

35、数 RAMP 逐 渐施加，图一出示了由作用力引起的孔隙压力轮廓和矢量图。 实现在现实中一个 很短的时间段内塑性流动的发生是非常重要的。这里词“ FLOW”与“ FLOW FLUID ”具有误导性，它是在瞬时发生的。因此，不排水分析可以实现。 1.5.4 耦合流动和力学计算 在FLAC3D中全水-力耦合计算的发生具有两个方向： 孔隙压力的变化引起 力学应变从而影响应力；反过来，应力变化则影响孔隙压力。 在进行任何耦合分析之前比较时间比例是很有好处的。 因为时间步必须与从 一个网格点传到另一个网格点需要的时间一致。 经典地，与力学计算联系的时间步具有如下的形式： K 4/3G Lc 这里是固体密度

36、。 因此，流-力的时间步比率具有如下形式: _ K 4/3G Lc i kM 在大多数情况下，M 近似的等于10 pa，但是渗透系数可能不同，典型值是花 岗岩 10 J9m2 / Pa -sec，石灰岩 10 47 m2 / P - sec，砂岩 10_l5m2 / P sec，粘土 10 m2 / Pa-sec，沙子10 m2 / Pa-sec。对于岩石和土其数量级为103kg /m3。 这里 K 4/3G 大约为1010 Pa。对于岩石而言，式（1-73）流体与力学计算的时 间步的比率大约为：对于沙子Lc，粘土 106 Lc，砂岩108 Lc，石灰岩1010 Lc和花 岗岩1012 Lc。

37、如果我们排斥渗透系数比粘土大的材料，可以观察到，即使 Lc保 持10*的数量级，这个比率会保持一个很大值。事实上，当对比流体传播影响时， 力学影响可以假设为一个瞬时过程，这也是 FLAC3D中采用的方法。这里没有任 何时间是与力学时间步联系在一起的。 在大多数的模拟中，初始力学条件必须与在耦合计算开始之前必须达到的平 衡状态相一致。具有代表性地，在小的无量纲流动时间内，为了使包含不同时间 步的系统达到协调一致，对于每一个流动时间步，一定数量的力学计算时间步必 须获得。(1-72) (1-73) 对于大量无量纲的流动时间步， 如果系统接近持续流动状态，在没有力 学扰动的的介质中，一定数量的时间步

38、可以获得。一个合适的数值模拟可以在单 一流动和耦合计算之间人工控制切换。 用 SOLVE 命令的组合可以避免复杂的程序： SET MECH FORCE 将设定一 个最大不平衡力的限制，低于次此值则认为达到平衡态。SET MECH SUBSTEP n SLAVE 命令表明力学模块作为计算的从属部分，必须运行 n 个子循环（或，如 果检测到平衡后，也可低于这个值）。如果对于每一个流动时间不只需要一个力 学时间步就可以达到平衡态，那么这个流动子时间步的值翻倍，但是不能超过 m 值（默认m=100）。不管什么时候，一个力学时间步对于每一个流动群的子步的连 续性被打破，那么系统恢复到初始的设定状态。如果

39、 AGE 命令在 SOLVE 中设 定，那么计算将一直持续直到用 AGE 定义的流动时间达到为止。 在第三步，当力学模块和流动模块都运行时，那么就用到了 STEP 命令。在 这种情况下，对于每一个流动时间步都将获取一个力学时间步。 这里，流动时间 步假设很小以致于使一个力学时间步就可以达到平衡态。 为了描述一个完全的耦合分析，就上例我们继续进行模拟。在计算程序过程 中，屏幕应该可以被观察，在 SOLVE 命令执行后，有 7 个变量被更新：1当前 时间步；2对于主控制程序的总时间步；3在附属控制程序下的总时间步；4当 前激活的程序；5当前的子循环数；6当前最大不平衡力；7.总时间。 不平衡力限制

40、为1.5 103，足够的力学计算时间步以保证这个最大不平衡力 低于这个限制。然而，在本例中，用 SET MECH SUBSTEP 定义的时间步为 100。 如果真实的时间步总是等于用 SET MECH SUBSTEP 定义的值，那么肯定有些东 西错了。要么不平衡力限制或 SET MECH SUBSTEP 太低，要么就是系统不稳定 达不到平衡。解的质量依赖于不平衡力限制：一个小的不平衡力限制将得到一个 平滑的准确的反映，但是运行很慢。 1.6对于渗流分析的输入指导 1.6.1 FLAC 3D 命令 下列命令用来进行流动问题的分析，注意，在进行流动分析前必须输入 CONFIG FLUID 命令。

41、Apply keywordvalue vkeywordxra nge Apply 命令用来给模型的内外边界设定流动边界条件。 用户必须确定设 定的 keyword 的类型。三个 keyword 的类型用来分配流动边界条件。 网格点类型的 keyword: PWELL v V 是流动速率，单位是m3/sec，分配给在 RANGE 区域内的每个网格点。 本命令用来沿流动边界设定恒定的流入和流出值。用 INTERIOR 这个 keyword 作用于内部的网格点。当一个新的井被设定时，旧井将被代替。 面类型的 keyword: DISCHARGE v 作用于法向边界的特定流出矢量的流体通量部分 v,单

42、位 m/sec。 Leakage v1 v2 V1 渗透层的孔隙压力 V2 渗透系数 h,单位是m3/Nsec ZONE 类型 keyword: VWELL v 流动的体积速率 v,（每单位时间每 ZONE 体积的流体体积），对于流入 情况下 v0。 CONFIG FLUID FIX ppvrange 在点处的孔隙压力被固定。 FREE PP RANGE . 释放孔隙压力变化固定。 HISTORY vid nhxnstep=nkeyword x y z 一定量的网格点变量可以被记录， 应用下列关键词。 gp pp 点（x,y,z）孔隙压力的记录 zone pp 点（x,y,z） 点的 zone

43、 孑 L 隙压力的记录 真实时间记录 fltime 流动的真实的时间的记录 INITIAL keywordv keyword value Contour pp 孔隙压力轮廓。 Flow 具体的流量矢量。 PRINT keywordv keyword vrange gp v keyword BIOT 网格点 BIOT 模数。 Fmodulus 网格点流动模数。 PP 网格点孔隙压力 Zone v keyword Fdensity 流体密度。 PP 孔隙压力。 PROP PROPERTY keyword value vi=i1,i2 j=j1,j2 此命令为用 MODEL 命令定义的流动模型分配属

44、性。 (1) BIOT_C BIOT 系数:。 (2) Permeability 各向同性渗透系数 k (3) Porosit 孔隙率 n。 U_thc 不排水热力系数。 SET KEYWORDVKEYWORD VALUE FLUID KEYWORD v KEYWORD VALUE AGE T 用 SOLVE 命令解题时的流动计算时间限制。 BIOT ON/OFF 如果是 ON，则用到 BIOT 系数和模数M的 流动计算，如果是 OFF 则此时用到流体模量Kf和孔 隙率n，而且=1，默认是 OFF。 DT T 用 T 定义流动时间步，此时间步在隐式方法下必须定 义，默认下，FLAC3D自动计算

45、显式模式下的时间步。 此命令允许用户选择不同的时间步。 如果FLAC3D确 定用户选择的值对于数值稳定分析而言太大时，那么 时间步将被减少到一个稳定值。 Implicit on/off 流动隐式方法。 ON/OFF ON 流动计算停止，OFF 流动计算开始。 PCUT ON/OFF 默认下，在 FLAC3D中进行负孔隙压力计算， 如果是ON，那么零孔隙压力设定被停止。 RATIO VALUE 用 SOLVE 命令进行流动计算时，设定流动比率限制。 默认情况下，RATIO 被定义为在模型中所有网格点的 平均不平衡水流强度与平均水流强度之和。当计算时 低于此值时，那么计算就停止，默认的值为1.0

46、10*。 STEP VALUE 规定流动计算步。 SUBSTEP value 设定在耦合计算中的流动子步。默认值 100,此时流动计 算是藕荷计算中的附属部分。 MECHANICAL keywordv keyword value 此命令为耦合计算设定力学参数。 FORCE value 设定最大不平衡力限制。 ON /OFF 力学计算开关。 RATIO 意义同上，默认的值为1.0 10。 STEP VALUE 意义同上。 SUBSTEP valuevslave：意义同上。 SOLVE keyword value v keyword value AGE T 规定所有计算程序的最大时间限制。 CLO

47、CK T 是计算机运行时间限制。 WATER keyword value 此命令规定一个地下水位的位置和属性， 在有效应力计算中获得恒定 的流体静理学孔压分配 （不是在 CONFIG FLUID ） 。计算期间， FLAC 3D 应用在结构模型中的有效应力。 Zone 的孔隙压力计算出来并保存。孔 压并不受 Zone 体积改变的影响，也不受任何水流动的影响。注意，当 用到这个命令时，地下水表面以下饱和的介质密度和地下水表面以上的 干介质密度都要定义。应用下列关键词： DENSITY value 地下水密度。 TABLE keyword value 此命令为所有地下水位下的网格点定义孔压。 孔压

48、梯度方向沿给出的重 力矢量方向（ SET GRAVITY ）。地下水水位面可以以两个形式定义： 一个无限平面或者多边凸平面集合。 对于一个无限平面， 应用下列关键 词： ORIGIN X Y Z 平面上坐标为（X , Y, Z）的点。 NORMAL nx ny nz 用单位矢量 nx ny nz 来定义的平面的法向 矢量，并沿此方向孔压增加。 另外，可以用多边行平面定义，应用到下列关键词： face x1,y1,z1 xn,yn,zn vface 用从点 x1,y1,z1 到点 xn,yn,zn 进行定义的多边形平面，必须是 一个凸多边形， 它可以有任何数量的节点， 但是它却要分划成三 角形进

49、行保存，只有那些沿着重力方向的内凹面的网格点的孔压 才分配。并不检查面的交搭和交叉。 1.6.2 FISH 变量 在 FISH 函数中，下面的标量变量也提供以帮助进行流动分析。 Fldt 流动时间步 Fltime 流动时间 Flow_ratio 当前的流动比率 下面的网格点变量可以通过 FISH 函数进行改变： gp_pp 网格点孔压 下面的网格点变量可以通过 FISH 函数进行改变： gp_flow 在网格点的不平衡流量 下面的网格区域（ZONE）变量可以通过 FISH 函数进行改变: z_pp zone 孔压 z_qx 比流量矢量的 X 部分 z_qy 比流量矢量的 Y 部分 z_qz 比

50、流量矢量的 Z 部分 1.7验证举例 1.7.1 在限制层内的不稳定地下水流动 一个长堤坝，宽 100m并建立在一个浅层的饱和土层上。堤坝宽与土层的厚 度比较起来是非常大的，而且与土体的渗透系数 k =10J2比较起来，其渗透系数 是可以忽略的。土体的 BIOT 模数值为 M =10GPa。在均匀土层里，达到初始稳 定状态。目的 是研究地下水位上升 2m 时孔隙压力的变化情况。它要与堤坝水流 高水位边的孔隙压力 P1二H （这里水的密度 二1000 kg/m3重力加速度 g =10m/s2 ）上升保持一致。 土层内的流动被假设为单方向的，模型宽 100m。超孔隙压力 p 初始状态设 为零，而在

51、模型的一个末端突然增加到零。可压缩分析解具有如下形式： P（?t） =1 eF?（Si 血% （ 1-74） 兀心 n 这里z轴沿着堤坝宽方向， * Figure M.4 Confined flow in a soil layer ? = ， ? = ， t? = ct/L2和 c = Mk。 P1 L 在 FLAC3D模型中，土层被定义为具有 25 个 ZONE 的柱体。超孔隙压力在 z=0 表面上定义为2 104 Pa，而且在 z=100m 平面上设定为零，网格图形如下: 在 FISH 函数中运行的解析解，与时间分别为t?=0.05,0.1,0.2和 1.0。解析解 和数值孔压结果被保存在

52、 TABLE 中。 例 3 包含了用显式模式解得的 FLAC3D数据文件。例 4 包含了用隐式模式解 得的FLAC3D数据文件。图 6 展示了显式有限差分模式解与数值解的比较；图 7 展示了隐式有限差分模式解与数值解的比较。标准化超孔压力 p/pi与标准化 距离z/L的比较在两图中都有展示。TABLE2，4，6 中包含了超孔隙压力的解 JOib lile ： Uivlrac grDindhnkii I iiin cmirrd lnirr: r pKiilinelirKl Gr4lton|： 展旅： V: MKOP 4000 : MOOT Z: I Z： BT.TOT DE1：ZSI4e4aQ2

53、 htn：H Ar:2Z3i Grid Ax& L imlt Ikcca Co-rai 6-ig Grciii bi:. FsfrTE-qidE. hfriE-safiLJSA Figure M.S ILAC 3 h gh d for fl uul flow in 001 13332 7/ed Feb 12 1W Job IMcUwtpa gKMndaiM low inn ccnfnpd I刃: bnEiltr iod Table I Jitantd Livsde - lCXXje*W0c- 2HtlTTd IIVIUP 3-Wimrd Livst/le lOOOrCQO 4 Lrrs

54、lJe Livsde ioaje*coo 0 Jiivrrd LivsMe 1.( (r0e-or TXtanrd LrwsMe IOQJFCDO 1.000P4000 1.000*OCD 1.000P4000 1.00Dr000 1.000P4000 1.000*000 1.000 4000 1.000 *000 IfcRaCcfBifrnGfciu he. Mnpgis. h*TF$Ofc)USA Figure M.7 Comparison of excess pore pressures for the implicit-solution algorithm (analytical va

55、lues = lines; numerical values = crosses) 1.7.2 单方向固结 饱和土层厚度 20m ,如下图，位于一个坚硬的不可压缩的基础之上。在不 & 心 心 L 丄 山 丄 i H Figure M.8 One-dintension al consolidatioji 排水条件下有一个恒定的表面何载 Pz =105，土体基值均匀弹性。应用各向同性 的达西流动定律。初始压力由流体建立，随着时间的继续，流体通过土层表面排 出，把荷载传递到基础上。对于单方向的固结问题可以用 BIOT 固结理论解。 孔隙压力的扩散方程具有如下形式： 这里 1 a2 S - 存

56、储系数 M - K 4G 3 边界条件是：二zz二-PzHt , p=0，在 Z 二 H ( H(t)是步骤函数)处。Uz=0, =0，在 z =0 处。因为应力恒定。所以上式成为如下形式: :z . .2 :一 P _ C - 2 Z 具有如下边界条件：在 z = H 处 p=0，在 z = 0 处, -2 :z :d二； ：S (1-75) (1-75) SP :z 7. =0 由土层荷载引起的孔压初始值 P0可以由流体结构法则(P二M -zz ) 得到，同时需要考虑不排水 条件（ =0），并用单方 向的力学结构法则 （ zz = 1 - -P）根据应力的形式表示应变。并以下面形式给出:

57、a 1 p0 Pz C(1 S 孔隙压力的解具有如下形式： Po J.： sin am? mt? p =2 e =0 Pz mo am这里 (1-76) (1-77) 31 aI(2m 1) ? H - z ? , and H ?二 H2 垂直位移Uz通过考虑平衡方程F-JzzZ=0与力学结构方程-zz 立。根据；zz = -:u/；:z，我们得到上面的组合，通过考虑（1-77） 们可以得到： p0 ： cos（am?） “ Uz = 2a _ 竺 |送一 e mt 出 :-1 ；zz -p 来建 和边界条件我 （1-78） 这里Uz被定义为Uz二卫. HPz 用下面的材料属性来描述本例: 干

58、体积模里 K 5 108Pa 干剪切模量 G 8 2 10 Pa BIOT 模数 M 9 6 10 Pa BIOT 系数 a 1.0 渗透系数 k 苗0 Pa - s e c Pz 在 FLAC3D中本问题的网格模型是沿z轴向上，由 28 个尺寸为单位 1 的 ZONE组成的柱壮结构。柱底固定，并且在 X 和丫方向的横向位移都是固定的。 一个压力作用在柱顶。 U f IIK iltT-riftarrJ EEtiwIitalhvi llrtB3CDiE3H h nciiai pp IO LraiiMfe gUKSe-onzA BCm-ODI 2 F &H i ndfim ppcd krs

59、-rMv - .一* dTVE* V3L F 辭I 卯 i1l 吐归虫 4033 Figure M. 10 Comparison between 席打“(y/iffff anil n urrt erica I values of pore pressure in a( (me*duneHsignal consulidatiun test ELAC3l 2d)0 吕lep44W hfcdrlPermcliF O:H:Q Tir Feb 11 I看 atue-ou 卄国H耳* Lresle - Job Tile-: Ore- dmmaandlcorsaiickiiiCTii lilb 创 pwiw 辺l刑电 lir nl .越hliMnid4P|ilTl Hilton- 与有效应力有关的孔压传递显示在图 12 中 ItEcnCornJiiqGrocL 阮. hirEqxJfc. Irfe-iieEaEiJL Figure Mr 11 Cf/mparisan between (UHilytieal urul vmtterietti values far vertical displacetnent in a one -di WJ e H si on til cons