FLAC3D渗流说明书中文版.docx

1、FLAC3D渗流说明书中文版FLAC3D流固耦合（手册翻译）1.1简介FLAC3D通过具有渗流性的实体（比如土）来模拟流体的流动。流动模型的建立可以独立于力学计算而自动完成，或者说可以与力学模型同时建立，这样就可以考虑流体与土体之间的相互作用。流固耦合的一种类型是“固结”，即:空隙水压力逐渐消散而导致土体的沉降。这个过程包括两种力学反映：一，空隙水压的改变导致有效应力的变化，这将影响到土体的力学反映（如：有效应力的减小可能导致塑性区的产生）；二，力学实体中某一区域的流动会随着空隙水压的改变而改变。该程序可以计算完全饱和情况下的流动，也可以模拟具有自由水面的流动。模拟具有自由水面的流动时，自由水

2、面以上的部分空隙水压等于0，气相将不参与计算。对于不考虑毛细水压力颗粒较粗的材料可以采用这种模拟方法。流体计算就有以下特点：1 根据各项同性和各项异性的渗流计算，相应采用两种流体运动定律。流动中的null材料用来模拟流动范围内的非渗流材料。2 不同区域可以拥有不同的流动模型（isotropic, anisotropic or null）和模型参数。3 可以事先指定流体的压力、流量、非渗流区边界条件。4 流体源可以以电源，也可以以体源的形式插入到材料中，这些源对应于流体的流入或流出，可以随着时间而变化。5 对于完全饱和流动，可以采用显式和隐式两种算法，但对于非饱和流动则只能采用显示计算。6 任何

3、力学和温度计算模型都可以与流体模型一起使用，在耦合计算中，可以考虑饱和体的压缩性和热膨胀性。7.流体与力学计算的耦合通过提供比奥系数来实现。8与温度的耦合计算可以通过提供线性热膨胀系数t和不排水温度系数（undrained thermal coefficient,可能翻译的不对）来实现。9.热-流体计算以线性理论为基础，假定材料参数为常数，不考虑对流。流体与实体的温度保持局部平衡。非线性行为可以采用FISH语言改变孔隙压力、材料特性来实现。动力水压力的产生以及循环荷载引起的液化也可以用FLAC3D模拟（3.4.4节）。FLAC3D不考虑部分饱和体中的毛细水压、粒子之间的电、化学力。但是可以根据

4、材料的局部饱和情况，孔隙率或其他相关变量利用FISH语言施加内部应力来实现上述几种作用力。与此类似，由于气体溶解引起的流体刚度变化的影响也是不能直接模拟的，但是FISH编写一个压力，时间的函数就可以用来改变局部流体的模量本章主要包括7个方面的内容：1. 1.2节和1.3节中分别介绍流体力学，流-力耦合以及流-热耦合的数学理论模型和数学表达式。2. 1.4中介绍流体力学模拟采用的计算模型和相关的命令。3. 1.5节介绍流体力学模拟采用的材料参数，各参数的单位。4. 1.6介绍流体力学模型采用的不同的边界条件，初始条件，流源和汇集。5. 1.7介绍流体分析的一般过程，该部分还包括几个例子，建议读者

5、在做自己的流体分析以前首先研究清楚我们提供的例子。6. 1.8是几个考核题，用来验证flac3d流体分析的准确性。7. 1.9总结了所有与流体计算相关的命令，以供参考。读者应该首先熟悉力学计算，然后才能将流体，力一起来考虑。流-力耦合是一个复杂的问题，需要读者有很好的洞察力来判断结果的正确性。在开始计算一个大的项目之前，有必要在一个网格比较少的例子上作试验，尝试不同的边界条件和建模策略。“浪费”在试验上的时间，必定可以通过计算时间来弥补回来。（1.2 1.3 我看就没有必要翻译了，有兴趣看看朱伯芳院士的书就可以了，理论是一样的）1.2 FLAC3D流体固体方程描述(部分翻译)大多数工程分析，包

6、括瞬时变形和孔隙水压力消散机制通过非耦合的技术来模拟。FLAC3D除了提供这些非耦合的计算模式，也提供流体固体力学耦合计算的选项。在该方法中，多孔介质的力学响应通过瞬态流来研究。虽然本章节主要考虑变形流体消散的问题，流体热固体力学耦合作用的方程完整的列出。FLAC3D中力学变形流体消散的描述在准静态Biot理论的框架下完成，而且可以应用到多孔介质中遵循Darcy定律的单相渗流的问题。不同类型的流体，包括气体和水，可以用这个模型表述。描述多孔介质中流体渗流的变量是孔隙水压力，饱和度和特定排水向量的三个分量。相关的变量遵循流体的质量守恒定律，达西定律，和流体响应孔隙水压力改变，饱和度改变，体积应变

7、改变和温度改变的本构模型。1.3数值方程(部分翻译)流体连续性方程：在FLAC3D的数值方法中，渗流区域离散为由8节点定义的砖块状的zone。孔隙压力和饱和度设定为节点变量。每个zone可以划分为5个四面体，四面体中的孔隙水压力和饱和度假定为线性变化。数值方法基于流体连续性方程的节点有限差分方程。该方程可以与导致牛顿定律节点形式的固体力学常应力方程并行计算。1.4 流固耦合的计算模式孔隙水压的计算模型和命令取决于节点是否已经设置了流动（也就是说，CONFIG fluid命令是否已经采用）1.4.1和1.4.2将分别介绍这两种情况。方便起见：下面所有的命令在文末1.9都作了汇总。1.4.1节点没

8、有设置渗流分析如果不给出CONFIG fluid命令，仍旧可以在节点上设置孔隙水压力。这种计算模式下孔隙水压力不变，但是当材料为弹塑性体时有可能不收敛。使用initial pp 命令加a gradient, 或用 WATER table 命令可以指定节点上孔隙水压力的分布。如果使用WATER table 命令，水面一下部分的孔隙静水压力的分布程序将自动进行计算。这种计算模式，必须用WATER density设置流体密度，SET gravity 设置重力加速度。当用水面定义好以后，采用PRINT water命令可以输出流体密度，水面的位置，PLOT water命令可以显示出水面来两种情况下微元体

9、的孔隙压力计算是将个节点的值平均，从而计算有效应力。这种计算模式，体力的计算不考虑流体压力：水面以上和以下材料的天然重度和饱和重度要有用户自己指定。PRINT gp pp 和PRINT zone pp可以输出节点和体的孔隙水压力。PLOT contour pp 命令用来绘制节点孔隙水压的云图。在用户手册中的3.6.3有这种计算模式的一个简单例子供读者学习。1.4.2节点设置渗流分析如果使用了CONFIG fluid命令，就可以进行流体的瞬态分析。地下水面和孔隙压力都会发生改变。在CONFIG fluid计算模式下，有效应力（静止孔隙压力的分布）和不排水计算分析可以同时考虑，此外，可以进行完全耦

10、合的计算分析，这时孔隙水压的改变会导致变形同时体积应变又会引起孔隙水压的变化。如果节点配置了渗流，微元体的孔隙水压仍旧等于其各个节点上孔隙水压的平均值，但是这种计算模式用户必须指定材料的干密度（包括水面以上和以下的部分），因为在FLAC计算体力的时候考虑了流动的影响。当设置CONFIG fluid后，zone必须赋予渗流模型。MODEL fl_ isotropic命令指定为各向同性流动模式MODEL fl _anisotropic.命令指定为各向异性流动模式MODEL fl_null命令指定为不透水区域,注意力学分析指定的null区域不能自动成为流体中的不透水区域。流体的性质可以赋予单元体和节

11、点，单元体的流体性质可以用命令赋值，各向同性体包括：渗透系数、孔隙率、Biot(比奥)系数。对于各向同性的流动，渗透系数采用关键词perm指定。各项异性流动采用关键词：，指定三个主方向的渗透系数，用关键词fdip, fdd, frot.指定三个主方向。三个主方向的渗透系数与，相对应且服从右手定则。其中fdip和fdd分别是k1 k2两方向组成的平面的倾角和倾向角，其中，倾角是以从xy-面向z-轴负方向转动为正。倾向角是倾向在xy-面同轴正向的夹角（顺时针）。是k1-k2面内倾向顺时针转到k1方向的角度。图三个渗透主方向的定义默认条件下，比奥系数等于1，孔隙率为0.5。节点的流体性质用INITI

12、AL命令付值。其性质包括：流体密度，流体体积模量，比奥模数，饱和度。节点的性质可以给定变化范围，流体密度可以采用WATER命令一次性赋值。下面表1-1是相关参数及其关键字：表1-1 相关的参数的确定性质关键字范围命令渗透系数（各向同性）permin zonesPROPERTY主方向渗透系数（各项异性）k1in zonesPROPERTY主方向渗透系数（各项异性）K2in zonesPROPERTY主方向渗透系数（各项异性）K3in zonesPROPERTY渗透主方向定义（各项异性）Fdipin zonesPROPERTY渗透主方向定义（各项异性）Fddin zonesPROPERTY渗透主方

13、向定义（各项异性）frotin zonesPROPERTY孔隙率porosin zonesPROPERTY比奥系数Biot_cin zonesPROPERTY流体模量fmodat gridpointsINITIAL比奥模量biot modat gridpointsINITIAL饱和度satat gridpointsINITIAL流体拉极限ftensat gridpointsINITIAL流体密度fdensat gridpointsINITIALdensgloballyWATER有一点需要注意：流体的压缩性可以通过两种方式定义（设置CONFIG fluid）：（1）定义比奥系数 ，比奥模量；（2

14、）定义流体模量和孔隙率。第一种情况考虑到土颗粒的压缩性，当比奥系数等于1时，颗粒不具有压缩性，第二种情况假定土颗粒不可压缩。PRINT zone property用命令输出实体单元的性质，PRINT gp命令输出节点的性质。如给定了流体密度以及水面的位置，可以用PRINT water命令输出，流体流动的性质可以采用PLOT bcontour property命令或PLOT block property命令输出云图。对于各项异性流动，流体区域的整体渗透方向量可以使用实体单元的关键字kxx, kyy, kzz, kxy,kxz, kzz 输出。节点初始孔隙压力分布，CONFIG fluid和non

15、-CONFIG fluid的做法是相同的。即：采用INITIAL pp和WATER table命令。FIX pp命令和FREE pp命令可以为选定的节点施加的接触约束。源和汇可以采用APPLY命令。1.6节会介绍流体边界和初始条件。流体计算采用的是SET fluid and SOLVE命令。不同的计算过程可以通过选择。例如：SET fluid on or off命令为采用或取消流体力学的计算。这些命令和关键字的选择取决于流动分析的耦合程度。1.7节将介绍不同耦合程度应该采用的计算过程。该部分还给出了几个例子，包括：仅作流动的分析和流固耦合的分析等。流体力学的计算可以采用好几种性是查看其结果。P

16、RINT gp pp andPRINT zone pp命令可以输出节点和实体单元的孔隙水压力。HISTORY gp pp and HISTORY zone pp命令可以查看孔隙水压的变化。瞬态计算中，可以用HISTORY fltime命令记录流动时间的变化，从而可以输出，孔隙水压随时间的变化趋势。PLOT contour pp命令输出节点孔隙压力的云图。PLOT contour saturation输出饱和度云图。PLOT fluid输出流体排放方向量。CONFIG fluid模式计算的一般信息可以用P RINT fluid命令输出。某些流体变量还可以采用FISH语言编程得到。在1.9.2中列

17、出了这些变量，与节点相关的一个变量，gp flow只能通过FISH函数来得到它的值。1.5渗流分析中涉及到的属性和单位FLAC3D渗流计算中涉及的属性为渗透系数，流体的密度，或者Biot系数和Biot模量（颗粒可压缩土体中的渗流），或者流体体积模量和孔隙率（只适用于颗粒不可压缩的土体）。1.5.1渗透系数FLAC3D中各向同性渗透系数，在文献中也称为流度系数。该系数为达西定律中压力的系数，而且与水利传导系数有关。相关公式见手册。FLAC3D用到的K与我们通常说的渗透系数（cm/sec）的换算关系：K(in SI units)=Kh(in cm/sec)*1.02*10-6如果沿着节点渗透系数有

18、所变化，计算时间步由最大的渗透系数决定。在稳定流计算问题中，最好减少渗透系数的差异以提高收敛速度。例如，对于最后的稳定状态，20:1的渗透系数差别与200:1的渗透系数差别几乎对结果没有影响。1.5.2土体质量FLAC3D中涉及的土体密度有3种：干密度，饱和密度和流体的密度。只有重力荷载需要考虑的时候才设置土体的重度。如果FLAC3D设置了CONFIG fluid，那么干密度必须设置。FLAC3D会根据已知的流体密度（f）、空隙率（n）和饱和度（s）自动计算每个单元的饱和重度。唯一一种需要设置饱和重度的情况是在计算有效应力时没有设置CONFIG fluid。WATER TABLE(INITIA

19、L PP)命令指定地下水位的位置。水面以上的单元体（ZONES）设定干密度，水面以下的单元体（ZONES）设定饱和密度。土体的密度由INITIAL density命令设置。全局流体密度可以功过WATER density命令设置，或者可以用INITIAL fdensity命令设置不同位置的不同流体密度。所有的密度是zone的变量和质量密度（Kg/m3）。1.5.3流体模量1.5.3.1 Biot系数和Biot模量Biot系数定义为当孔隙压力改变时，吸收流体的体积与排出土体的流体体积的比例。该系数可以由测定排水体积模量的排水试验来决定。该系数的变化区域为，3n/(2+n)与1之间，n是孔隙率。在土

20、体骨架不可压缩的情况下，Biot系数为1。该值为FLAC3D的默认值。（计算公式参考手册）可压缩骨架的计算模式可以由命令SET fluid biot on打开。Biot系数是zone的属性，用PROPERTY命令定义。Biot模量是节点变量，由INITIAL命令定义。1.5.3.2流体体积模量在土体骨架的压缩性可以忽略的分析中，用户可以选择系统默认的Biot系数（=1），或者定义Kf/n等于Biot模量；或者直接定义流体的体积模量Kf。流体的“可压缩性” Cf，是体积模量的倒数（Kf=1/ Cf)。例如，对于室温下的纯水，其体积模量为2Gpa。在实际土体中，孔隙水含有溶解的空气气泡，这降低了水

21、的体积模量。对于地下水问题，考虑到空气含量的不同，水的体积模量在不同的节点可能有不同的模量。当流体体积模量作为输入参数时，对于颗粒不可压缩的土体，Biot模量由程序自动计算（式1.82）。计算中，孔隙率（zone的属性）由节点平均体积在节点上计算。Biot系数在渗流区域设置为1,忽略已经赋予的数值。1.5.3.3流体模量和计算收敛速度如果考虑稳定渗流，完全饱和渗流，Biot模量对于数值收敛来说不是很重要。因为，系统的响应时间和时间步的长度与Biot模量成反比。对于包含浸润线的系统，较低的体积模量可以加速收敛到稳定状态的速度。因为，饱和度的变化涉及到时间间隔，而不是Biot模量和时间间隔的乘积。

22、在流固耦合比较强烈的系统中，比较难以估计收敛速度和计算时间。我们可以做如下评价：高的模量比低的模量导致更快的收敛速度。无论如何，从数值计算的观点，没有必要将Biot模量设置为大于20倍的土体的体积模量。1.5.3.4排水和不排水分析的流体模量在FLAC3D中，当设置了CONFIG fluid命令，而且选择了输入流体的模量（或者Biot模量），必须为土体设置排水体积模量。土体的表象体积模量（不排水体积模量）由程序计算，而且随着之间更新。不排水分析也可以不设置CONFIG fluid命令来分析。在这种情况下，直接设置土体的不排水体积模量。不排水流体体积模量计算：1.5.4孔隙率孔隙率是一个无量纲的

23、数，由孔隙的体积与土体的总体积的比值定义。程序默认孔隙率为0.5.孔隙率的取值范围为01,但是必须谨慎对待小于0.2的情况。因为流体表象刚度和空隙率是成比例的（Kf/n）。对于，孔隙率比较小的情况，表象刚度与骨架刚度相比可能相差很多，导致了计算收敛缓慢。这种情况下，可以考虑折减流体刚度（Kf）。FLAC3D中的孔隙率用于计算饱和重度和Biot模量（当流体的模量作为输入参数时）。由于更新孔隙率非常耗时，而且仅影响到瞬态流浸润线的坡度，FLAC3D不会更新孔隙率。1.5.5饱和度饱和度定义为流体所占的体积与所有孔隙体积的比值。在FLAC3D中，如果任何一点的饱和度不是精确的等于1，该点的孔隙水压力

24、为0。流体中溶解的空气或者气泡，可以通过降低流体的刚度来实现，但是必须在饱和度为1的情况下实施。虽然在部分饱和区域的初始孔隙水压力为0,但是孔隙中的流体还是要考虑重量的。这些流体在自重下流动。用户可以设置初始饱和度，但是程序为了遵守质量守恒定律，会自动更新。有一点需要指出：饱和度不是一个独立的变量，其节点上的值不能固定。1.5.7流体的抗拉极限在细粒土中，孔隙水可以承受明显的拉力（负超孔隙水压力）。在FLAC3D中，负超孔隙水压力能够发展起来，直到发生“去饱和”现象（孔隙水从孔隙中流出，土体不再饱和）。负超孔隙水压的极限用INITIAL ftens命令设置，程序默认值为1015。需要指出：负超

25、孔隙水压与毛细压力，电化学力是不一样的。后者可以在本构模型上增加有效应力表现出来。负超孔隙水压与土体是否由颗粒组成没有关系。负超孔隙水压只是由于包含水的土体体积膨胀所引起。1.6渗流边界条件、初始条件、源与汇默认的边界条件是不透水边界，所有的节点上流体可以流动的。（例如，孔隙水压力在这些节点上是自由的，表现在孔隙水可以从临近单元流进或流出网格点）这种条件可以通过作用在节点上的命令FREE PP明确指定说明，相反的条件可以通过FIX PP命令施加在节点上实现。一般情况，如果孔隙水压力固定，流体会从外部边界进入或者流出。如下所示两种边界条件：1. 孔隙压力自由（FREE PP）：这是一个不可渗透边

26、界，而且是默认选项。节点上的流体与外界没有交换。孔隙压力与饱和度根据式1.69和式1.70计算。该项计算依赖于当前饱和度的值和流体是否已经受拉破坏(降低到流体抗拉极限)。2. 固定了的孔隙水压力（FIX PP）：这是一个可渗透边界，沿着该边界流体可以流到模型区域外。只有当孔隙压力设置为零时，饱和度饱和度才会变化。否则，饱和度必须为1.孔隙水压力的设置值不能小于流体的抗拉极限。否则系统将重置这样的压力到流体抗拉极限。如上所示，有几种边界条件的组合是不可能实现的。例如，节点的孔隙水压力设定值小于先前设定的流体的抗拉极限。FLAC3D在执行任何一个step之前会纠正这个错误。孔隙水压力可以通过FIX

27、 pp固定。或者，用APPLY pp命令设置外部或者内部的边界条件。需要指出，如果不是设定外边界面上的孔隙水压力，必须用interior关键词。用APPLY命令的一个好处就是，可以用history控制孔隙水压力。单相渗流的边界条件可用APPLY命令应用到单个range或者几个range，节点，zone的表面或者zone。APPLY pwell命令将预定义的流入或者流出量应用到边界的节点上。如果加了interior关键词，该条件可应用到内部节点。APPLY DICHARGE和APPLY LEAKAGE命令可以在边界单元的表面指定流体通量或渗漏边界条件。APPLY VWEEL命令用于在指定范围流入

28、单元的流体体积流速。所有的这些命令，除了APPLY LEAKAGE，都可以通过HISTORY命令监测。固定孔隙压力的节点的作用有点类似于井。在节点上的渗流量没有现成的命令可以利用去监测。然而，FISH中的gp_flow记录了节点的不平衡流量；可以写一个简单的程序监控节点的流量。孔隙压力，孔隙率，饱和度和流体属性的初始分布可以用INITIAL命令或者PROPERTY命令定义。如果设置了重力加速度，以上变量的初始分布必须和由重力加速度和节点上的水的密度，饱和度和孔隙率引起的重力梯度一致。如果初始的孔隙压力与静水压力不匹配，以开始计算zone内就产生渗流。因此，在开始模拟的时候必须先运算几个step

29、以检验这种可能性是否存在。如果模型中包含有接触面，有效应力将会沿着接触面初始化。例如，当节点上的应力初始化是孔隙水压力的存在将会计入接触面的应力（即：在节点应力初始化时，认为接触面应力包含孔压）。water lable命令将包含沿着接触面的孔压，这是因为定义在单元节点上的孔压也在接触面节点上。有一点需要注意：如果接触面的上下两面连在一起，在没有阻力时，将发生穿越接触面的流体流动。但程序不对沿着接触面的流体流动（裂隙流）进行计算。1.7单渗流和流固耦合问题（重点看例题）FLAC3D既能进行单渗流分析，也能进行固流耦合分析。耦合分析可由FLAC3D内置力学模型完成。但要注意，渗流模型中的空单元并不

30、是力学空单元。必须用命令model fl_null给单元赋予流体空属性。对于耦合过程，FLAC3D提供了几种计算模式。其中之一是假设孔压一旦被赋予便不再改变。该方法并不要求任何额外空间存储计算过程。在1.4.1章节将讨论这部分内容的相关命令。除此之外涉及到渗流的计算模式都要求使用命令config fluid。而命令model fl_iso使所有单元中都能发生渗流。在1.4.2章节将讨论有关与config fluid相关的命令。 不同的耦合计算模式在下面讨论，越是复杂的越需要更多的计算机内存和时间。一般情况下，在能跟所模拟问题的物理过程相似的情况下，应使用尽可能简单的模式。计算模式的选择根据以下

31、几个方面确定。1.7.1时间比例 对所需模拟的渗流或耦合问题用FLAC3D估计与涉及的不同进程相关的时间比例是非常有用的。对有关研究问题的时间度量和扩散性的认识有助于估计最大网格宽度、最小区域尺寸、时步大小和计算可行性。如果不同进程的时间比例相差太大，则很可能采用一种简单的（非耦合）方法。（完全耦合分析过程会在1.7.2章节详细讨论） 时间比例可用特征时间给出。以下这些由量纲分析得出的定义，都是基于解析的连续源理论表达式。它们可用于得出FLAC3D分析的大致时间比例。力学过程特征时间注：Ku是未排水体积模量，G是剪切模量，是质量密度，Lc是特征长度（如模型的平均尺寸）流体扩散过程特征时间注：Lc是特征长度（如通过模型的流体路径平均长度），c是流体扩散率，定义为流度系数k和储水系数S的比值。流体扩散率 FLAC3D中使用了取决于控制过程的储水系数的几种形式：1. 流体存储系数2. 地下潜水相存储系数3. 弹性存储系数注：M是比奥模量，是比奥系数，K是排水体积模量，G是剪切模量，w是流体密度，g是重力加速度，Lp是特征储水长度。以上定义，有几点特性值得注意：（1）因为FLAC3D中显式的时步对应于最小区域中信息从一个节点传到下一节点所需要的时间，时步的大小可用