Fluid-Mechanical Interaction.doc

1 流体－固体耦合与单相渗流 1.1介绍 FLAC3D模拟了流体流过可渗透的介质，例如土体。渗流模型可以独立于通常FLAC3D的固体力学计算，而只考虑渗透；或者为了描述流体和固体的耦合特性，与固体模型并行计算。固结是一类流固耦合的现象，在固结过程中孔隙随压力逐渐消散，从而导致了固体的位移。这种行为包含了两种力学效应。其一，孔隙水压力的改变导致了有效应力的改变，有效应力的改变影响了固体的力学性能，例如有效应力的降低可能引发塑性屈服；其二，土体中的流体对孔隙体积的变化产生反作用，表现为孔隙水压力的变化。 本程序不仅可以解决完全饱和土体中的渗流，也可以分析有浸润线定义的饱和与非饱和区的渗流计算。该条件下，浸润面以上的土体的孔隙水压力为零，气体的压力考虑成负的。这种方法用于颗粒比较粗的毛细现象可以忽略的土体。 渗流分析中有如下的特征 1. 对应于渗流各向同性和各向异性材料采用不同的定律。渗流区域中的不可渗透的区域用流体的null材料定义。 2. 不同的zone可以赋予不同的渗流模型和属性。 3. 流体压力，涌入量，渗漏量和不可渗透边界都可以定义。 4. 土体中可以加入抽水井，考虑成点源或者体积源。 5. 计算完全饱和土体中的渗流问题，可以采用显式差分法或者隐式差分法；而非饱和渗流问题只能采用显式差分法。 6. 渗流模型可以和固体力学模型和传热模型耦合。 7. 流体和固体的耦合程度依赖于土体颗粒（骨架）的压缩程度，用Biot系数表示颗粒的可压缩程度。 由于循环荷载引起的动水压力和液化问题也可以用FLAC3D模拟。FLAC3D不考虑毛细现象，土体颗粒间的电化学作用力。然而，可以根据土体的饱和度，孔隙率，或者其他的变量，通过编写一段FISH语言来考虑这种力。类似的，由于液体中溶解了空气而引起的液体刚度变化，也不能显式的模拟，而通过FISH将液体刚度表示为压力，时间和其他变量的函数。 这以章节可以分为七个主要部分 1. 数学模型描述和相应的数值方法（单相渗流和流固耦合计算）。 2. 计算模式和渗流分析相关的命令。 3. 渗流分析所需要的材料属性的单位。 4. 不同边界条件，初始条件，和流体的源和汇的的描述。 5. 求解单相渗流问题和流固耦合问题的推荐方法。 用户在尝试解决流固耦合问题之间，强烈建议先熟悉一下FLAC3D求解固体力学问题的步骤。流固耦合的力学行为通常非常复杂，需要用户对计算结果是否正确具有远见和判断。在开始计算一个大的项目之前，有必要在一个网格比较少的例子上作试验，尝试不同的边界条件和建模策略。“浪费”在试验上的时间，必定可以通过计算时间来弥补回来。 1.2 FLAC3D流体－固体方程描述 大多数工程分析，包括瞬时变形和孔隙水压力消散机制通过非耦合的技术来模拟。FLAC3D除了提供这些非耦合的计算模式，也提供流体固体力学耦合计算的选项。在该方法中，多孔介质的力学响应通过瞬态流来研究。虽然本章节主要考虑变形－流体消散的问题，流体－热－固体力学耦合作用的方程完整的列出。 FLAC3D中力学变形－流体消散的描述在准静态Biot理论的框架下完成，而且可以应用到多孔介质中遵循Darcy定律的单相渗流的问题。不同类型的流体，包括气体和水，可以用这个模型表述。 描述多孔介质中流体渗流的变量是孔隙水压力，饱和度和特定排水向量的三个分量。相关的变量遵循流体的质量守恒定律，达西定律，和流体响应孔隙水压力改变，饱和度改变，体积应变改变和温度改变的本构模型。 1.3数值方程部分翻译 流体连续性方程 在FLAC3D的数值方法中，渗流区域离散为由8节点定义的砖块状的zone。孔隙压力和饱和度设定为节点变量。每个zone可以划分为5个四面体，四面体中的孔隙水压力和饱和度假定为线性变化。 数值方法基于流体连续性方程的节点有限差分方程。该方程可以与导致牛顿定律节点形式的固体力学常应力方程并行计算。 1.4流固耦合的计算模式 孔隙水压力分析所需的计算模式和命令依靠于节点是否设置了渗流分析。 1.4.1节点没有设置渗流分析 如果没有设置CONFIG fluid命令，也可以在节点上设置孔隙水压力。在该计算模式中，孔隙水压力不变，但是由有效应力决定的具有弹塑性力学性能的材料有可能破坏。 可以用INITIAL pp命令或者用WATER table命令，设置节点上的孔隙水压力分布。 如果设置了WATER table命令，在浸润线以下的静水压力分布由程序直接计算。在该情况下，流体密度和重力加速度必须设置。流体密度和浸润线的位置可以用PRINT water命令显示。 Zone的孔隙压力通过平均节点值来确定。孔隙水压力用于生成有本构模型中的效应力。在该计算模式下，计算体积力时，流体（质量）部分不能自动考虑，湿重度和干重度必须由用户手动设置。 1.4.2节点设置渗流分析 如果设置了CONFIG fluid命令，可以进行瞬态渗流分析。孔隙水压力的改变随着浸润线的改变而改变。有效应力计算和不排水计算可以在CONFIG fluid模式下计算。而且，也可以进行完全流固耦合计算。在该过程中，孔隙水压力的改变产生了变形，体积应变导致了孔隙水压力的改变。 如果节点配置了渗流，zone的孔隙压力可以通过平均节点的孔隙水压力值获得。在该模式下，土体干密度必须由用户手动设置（浸润线上下都设置干密度）。FLAC3D在该模式下，自动将流体质量计入体积力部分。 当设置CONFIG fluid后，zone必须赋予渗流模型。各向同性流用MODEL fl_isotropic命令来定义，各向异性渗流由MODEL fl_anisotropic定义，不可渗透材料由MODEL fl_null定义。有一点需要注意开挖掉的单元不会自动设置为不可渗透单元。 流体的属性不是赋值给zone就是赋值给grid。Zone中的流体属性由命令PROPERTY设置各向同性渗透系数，孔隙率，Biot系数等。 对于各向同性渗流，渗透系数由perm关键词定义。对于各向异性流，三个主方向的渗透系数用关键词k1，k2，k3来定义，主方向可以用关键词fdip，fdd，和frot来定义。 如果没有设置，Biot系数默认为1,孔隙率默认为0.5。节点流体属性可以用INITIAL命令来设置。这些属性包括，流体密度，流体体积模量，Biot模量，流体抗拉极限和饱和度。流体密度可以用全局命令WATER命令设置。 有一点需要注意流体的压缩性可以通过两种方式定义（设置CONFIG fluid）1定义Biot模量和Biot系数；2定义流体体积模量和孔隙率。第一种方法考虑了土体颗粒的可压缩性，第二种方法假设为固体颗粒不可压缩。 节点上初始孔隙压力分布可以按照1.4.1中的方法定义。孔隙水压力在特定的节点上可以用FIX pp命令固定。源和汇可以用APPLY命令定义。 渗流分析由SET fluid和SOLVE命令控制。一些关键词可以用于帮助求解过程。例如，SET fluid on或者 off打开或者关闭流体渗流计算。这些命令和关键词的应用取决于流固耦合的程度。 1. 5渗流分析中涉及到的属性和单位 FLAC3D渗流计算中涉及的属性为渗透系数，流体的密度，或者Biot系数和Biot模量（颗粒可压缩土体中的渗流），或者流体体积模量和孔隙率（只适用于颗粒不可压缩的土体）。 1.5. 1渗透系数 FLAC3D中各向同性渗透系数，在文献中也称为移动系数。该系数为达西定律中压力的系数，而且与水利传导系数有关。相关公式见手册。 如果沿着节点渗透系数有所变化，计算时间步由最大的渗透系数决定。在稳定流计算问题中，最好减少渗透系数的差异以提高收敛速度。例如，对于最后的稳定状态，201的渗透系数差别与2001的渗透系数差别几乎对结果没有影响。 1.5.2土体质量 FLAC3D中涉及的土体密度有3种干密度，饱和密度和流体的密度。只有重力荷载需要考虑的时候才设置土体的重度。 如果FLAC3D设置了CONFIG fluid，那么干密度必须设置。FLAC3D会自动计算每个单元的饱和重度。 唯一一种需要设置饱和重度的情况是在计算有效应力时没有设置CONFIG fluid。 土体的密度由INITIAL density命令设置。流体密度可以功过WATER density命令设置，或者可以用INITIAL fdensity命令设置不同位置的不同流体密度。所有的密度是zone的变量。 1.5.3流体模量 1.5.3.1Biot系数和Biot模量 Biot系数定义为当孔隙压力改变时，吸收流体的体积与排出土体的流体体积的比例。该系数可以由测定排水体积模量的排水试验来决定。该系数的变化区域为，3n/2n与1之间，n是孔隙率。在土体骨架不可压缩的情况下，Biot系数为1。该值为FLAC3D的默认值。（计算公式参考手册） 可压缩骨架的计算模式可以由命令SET fluid biot on打开。Biot系数是zone的属性，用PROPERTY命令定义。Biot模量是节点变量，由INITIAL命令定义。 1.5.3.2流体体积模量 在土体骨架的压缩性可以忽略的分析中，用户可以是用薄人的Biot系数，定义Biot模量；或者直接定义流体的体积模量。 流体的“可压缩性”，是体积模量的倒数。例如，对于失稳下的纯水，其体积模量为2Gpa。在实际土体中，孔隙水含有溶解的空气气泡，这降低了水的体积模量。对于地下水问题，考虑到空气含量的不同，水的体积模量在不同的节点可能有不同的模量。 当流体的模量作为输入参数时，对于颗粒不可压缩的土体，Biot模量由程序自动计算。计算中，孔隙率，zone的属性由节点平均体积在节点上计算。Biot系数在渗流区域设置为1,忽略已经赋予的数值。 1.5.3.3流体模量和计算收敛速度 如果考虑稳定渗流，完全饱和渗流，Biot模量对于数值收敛来说不是很重要。因为，系统的响应时间和时间步的长度与Biot模量成反比。对于包含浸润线的系统，较低的体积模量可以加速收敛到稳定状态的速度。因为，饱和度的变化涉及到时间间隔，而不是Biot模量和时间间隔的乘积。在流固耦合比较强烈的系统中，比较难以估计收敛速度和计算时间。我们可以做如下评价高的模量比低的模量导致更快的收敛速度。无论如何，从数值计算的观点，没有必要将Biot模量设置为大于20倍的土体的体积模量。 1.5.3.4排水和不排水分析的流体模量 在FLAC3D中，当设置了CONFIG fluid命令，而且选择了输入流体的模量（或者Biot模量），必须为土体设置排水模量。土体的表象模量（不排水模量），由程序计算，而且随着之间更新。 不排水分析也可以不设置CONFIG fluid命令来分析。在这种情况下，直接设置土体的不排水体积模量。 1.5.4孔隙率 孔隙率是一个无量纲的数，由孔隙的体积与土体的总体积的比值定义。程序默认孔隙率为0.5.孔隙率的取值范围为01,但是必须谨慎对待小于0.2的情况。对于，孔隙率比较小的情况，表象刚度与骨架刚度相比可能相差很多，导致了计算收敛缓慢。这种情况下，可以考虑折减流体刚度。 FLAC3D中的孔隙率用于计算饱和重度和Biot模量（当流体的模量作为输入参数时）。由于更新孔隙率非常耗时，而且仅影响到瞬态流浸润线的坡度，FLAC3D不会更新孔隙率。 1.5.5饱和度 饱和度定义为流体所占的体积与所有孔隙体积的比值。在FLAC3D中，如果任何一点的饱和度不是精确的等于1，该点的孔隙水压力为0.流体中溶解的空气或者气泡，可以通过降低流体的刚度来实现，但是必须在饱和度为1的情况下实施。虽然在部分饱和区域的初始孔隙水压力为0,但是孔隙中的流体还是要考虑重量的。这些流体在自重下流动。 用户可以设置初始饱和度，但是程序为了遵守质量守恒定律，会自动更新。有一点需要指出饱和度不是一个独立的变量，其节点上的值不能固定。 1.5.7流体的抗拉极限 在细粒土中，孔隙水可以承受明显的拉力（负超孔隙水压力）。在FLAC3D中，负超孔隙水压力能够发展起来，直到发生“去饱和”现象（孔隙水从孔隙中流出，土体不再饱和）。负超孔隙水压的极限用INITIAL ftens命令设置，程序默认值为－1015。需要指出负超孔隙水压与毛细压力，电化学力是不一样的。后者可以在本构模型上增加有效应力表现出来。负超孔隙水压与土体是否由颗粒组成没有关系。负超孔隙水压只是由于包含水的土体体积膨胀所引起。 1. 6渗流边界条件，初始条件，源与汇 默认的边界条件是不透水边界，所有的节点上流体可以流动的。一般情况，如果孔隙水压力固定，流体会从外部边界进入或者流出。如下所示两种边界条件 1. 孔隙压力自由这是一个不可渗透边界，而且是默认选项。节点上的流体与外界没有交换。孔隙压力与饱和度根据式1.69和式1.70计算。该项计算依赖于当前饱和度的值和流体是否已经受拉破坏。 2. 固定了的孔隙水压力这是一个可渗透边界，沿着该边界流体可以流到模型区域外。只有当孔隙压力设置为零时，饱和度才会变化。否则，饱和度必须为1.孔隙水压力的设置值不能小于流体的抗拉极限。 如上所示，有几种边界条件的组合是不可能实现的。例如，节点的孔隙水压力设定值小于先前设定的流体的抗拉极限。FLAC3D在执行任何一个step之前会纠正这个错误。孔隙水压力可以通过FIX pp固定。或者，用APPLY pp命令设置外部或者内部的边界条件。需要指出，如果不是设定外边界面上的孔隙水压力，必须用interior关键词。用APPLY命令的一个好处就是，可以用history控制孔隙水压力。 单相渗流的边界条件可用APPLY命令应用到单个range或者几个range，节点，zone的表面或者zone。APPLY pwell命令将预定义的流入或者流出量应用到边界的节点上。如果加了interior关键词，该条件可应用到内部节点。 固定孔隙压力的节点的作用有点类似于井。在节点上的渗流量没有现成的命令可以利用去监测。然而，FISH中的gp_flow记录了节点的不平衡流量；可以写一个简单的程序监控节点的流量。 孔隙压力，孔隙率，饱和度和流体属性的初始分布可以用INITIAL命令或者PROPERTY命令定义。如果设置了重力加速度，以上变量的初始分布必须和由重力加速度和节点上的水的密度，饱和度和孔隙率引起的重力梯度一致。如果初始的孔隙压力与静水压力不匹配，以开始计算zone内就产生渗流。因此，在开始模拟的时候必须先运算几个step以检验这种可能性是否存在。 如果模型中包含有接触面，有效应力将会沿着接触面初始化。例如，当节点上的应力初始化是孔隙水压力的存在将会计入接触面的应力。有一点需要注意如果接触面的一边与另一边没有接触，从一接触面到另一接触面之间的渗流不考虑阻力。