FLAC3D流力耦合作用.docx

FLAC3D流力耦合作用1.1耦合作用简介11. 2数学模型描述21.1.1 规定和定义21.1.2 流体重量平衡方程31.1.3 流动法则41.1.4 力学结构法则41.1.5 边界及初始条件51. 3数值公式51.1.1 空间导数的有限差分近似51.1.2 质量平衡方程的节点公式61.1.3 显式有限差分公式81.1.3.1 稳定标准91.1.4 隐式有限差分公式91.1.4.1 收敛准则111.1.5 力学时间步和力学稳定性121.1.6 总应力修正121. 4流动耦合问题的属性和单位121.1.1 渗透系数131.1.2 Biot系数G和Biot模数M131.1.3 流体体积模量141.1.4 孔隙率141.1.5 密度141.1.6 流体张力限151. 5单一流动问题和耦合流动问题151.1.1 恒定孔压（用于有效应力计算）151.1.2 建立了孔压分配的单一流动计算161.1.3 非流动，力学变形产生的孔隙压力161.1.4 耦合流动和力学计算171. 6对于渗流分析的输入指导181.1.1 FLAC3D命令181.1.2 FISH变量211.7验证举例221.7.1 在限制层内的不稳定地下水流动221.7.2 单方向固结251.7.3 穿透浅含水层限制边界的井水流动291.1 耦合作用简介FLAC3D允许在饱和多孔材料中进行流体流动的瞬时模拟。

流动计算可以脱离FLAC3D中的力学计算独立进行，也可以与其他力学模型进行耦合计算，以控流流一一固耦合作用的影响，其计算具有如下特征1. 提供了在各向同性条件下的流体运动法则，也提供了在流动区域中的无渗流材料的流动零模型2. 不同的区域可以有不同的流动模型和法则3. 流体压力、流量通量、松散和不透水边界可以设定4. 流体源可以以点源和体积源插入在材料中这些流源可以是随时间变化的流入源和流出源5. 对于流体计算既提供了显式算法也提供了隐式算法模式6. 任何力学和热学模型都可以和流体模型进行耦合计算在耦合计算中，流体和固体骨架的压缩和热膨胀都是允许发生的7. 由于骨架变形产生的流一一力耦合效果通过Biot固结系数«来实现8. 流体与热之间的耦合作用可以通过热力学线性膨胀系数必和不排水热力学系数P来实现9. 流动法则是基于线性理论的，其假设恒定的材料属性并忽略对流流体与固体骨架的温度保持局部平衡非线性变形则要通过FISH函数设定孔压和材料的属性来实现1.2 数学模型描述在FLAC3D中的变形传播增量的准确描述过程在Biot半静态线性理论框架内进行下面详细介绍水-热-力耦合与热-孔压-力的方程。

这些方程可以被被应用于孔隙介质中的达西流动定律各种各样的流体-包括空气和水-都可以应用本模型包含在描述经过孔隙介质的流动的变量是孔隙压力和比流量矢量的三个组成部分这些变量是通过流体运动的重量平衡方程-达西定律和描述流动对于孔隙压力、体积应变和温度的改变的反应机制的方程联系在一起的孔隙压力和温度的影响被包含在完全水热力耦合的力学结构法则增量中假设体应变率和温度率已知，把重量平衡方程-达西定律-代入流动结构关系中，根据孔隙压力产生一个不同的方程，此方程可以根据具体的几何特性、属性、边界和初始条件解出1.2.1 规定和定义作为一个习惯符号，符号向表示笛卡儿坐标系中矢量［a］的i部分；Aj表示张量［A］的（i,j）部分而f,i用来代表与为相联系的f的偏导数部分（符号f表示变量的矢量或标量部分）爱因斯坦的求和方程只应用于有下标i,j,k的变量，下标的取值分别为1,2,3代表空间组成当用于矩阵方程时，下标可以取任何值我们用国际单位制来描述变量参数下面这个无量纲参数来描述瞬时流动特性，特性长度Lc：LcVfaf(1-1)这里Vf表示流动区域的体积,af表示流动区域的表面积流体扩散率:(1-2)K是体这里k是渗透率，M是Biot模数，口是Biot系数，a1=K+3G/4积模量，而G是排水弹性材料的剪切模量。

Biot系数考虑孔隙材料土骨架的压缩性如果口取单位1,那么土骨架被认为是无压状态的，而Biot模数M等于,这里kw是流体的体积模量，而n是孔隙压力流体的扩散率成为下面的式子：(1-3)kc=nFkw-1注意，对于单一流动不考虑耦合效应的计算，流体扩散率成为下面的式子:(1-4)c=kKI特性时间:tc(1-5)1.2.2 流体重量平衡方程对于小变形的问题，流体的重量平衡方程表述为：加-qi,i+qv=—(1-6)二t这里qi为比流量矢量单位是［m/s］,qv是体积流量源强度，单位是1sec］，而-二是每单位体积孔隙材料的流体体积或流体容量在FLAC3D中，流体体积变量的改变与孔隙压力p、力学体积应变名和温度T的改变成线性关系流体结构法则表述为：(1-7)二二」卫-h_.£T-:tM：:t：:t：:t这里M是Biot模数，单位是N/m2】，口是Biot系数，P是不排水热力系数,单位［1"C］,他们都考虑土骨架和流体的热力膨胀因素把（1-6）代入（1-7）得到：-qi,i+qV=/与（1-8）M::t这里qjv-喙+哈（1-9）1.2.3流动法则基本的法则是定义比流量矢量与孔隙压力之间关系的达西定律。

对于均质，各向同性土体和恒定流密度条件下，这个法则以下列形式给出：qi=—kb—PfXjgj］i（1-10）这里k是渗透系数，单位是［m4/Ns］,Pf是流体密度，单位［kg/m3］,gi,i=1,3是重力加速度的三个分量，单位是［m/s2］为了供将来参考，量：p-，xjgj*=-——f-^1（1-11）汗g被定义为重力头，而Pfg4定义为压力头1.2.4 力学结构法则关于孔隙压力的体应变的影响在流体力学法则中反映出来依次地，孔隙压力的改变引起力学变形的发生：孔隙固体结构方程的增量形式表述如下形式：△可十口即久=Hi；Bij,A%-△*：）（1-12）这里也是关联循环应力增量，Hj是给定的函数，上】是总应变，bT】是热应变，而与是克罗内克尔增量尤其，弹性关系可以用下式表示：%—仃0+a（p—po乐=2G佃j—鸟：计^2（%卜一4）（1-13）这里%=K-2G/3,50和p0代表初始状态1.2.5 边界及初始条件用（1-10）代替（1-8）中的qi,产生不同的流动方程，假设qV已知初始条件与给定的压力区一致边界条件通常根据孔隙压力或垂直与边界的比流量矢量来表述在FLAC3D中，定义了四个类型的条件：（1）给定孔隙压力；（2）给定垂直于边界的比流量矢量；（3）渗透边界；（4）不渗透边界。

在flac3D中，渗透边界有如下形式：qn=h（p-Pe）（1-14）这里qn是垂直于边界的比流量矢量的外垂直线方向的部分，h是渗透系数，单位[m3/N-sec],p是在边界表面的孔隙压力，而pe是在渗透层的孔压默认情况下，边界条件是不透水边界1.3数值公式在FLAC3D中，重力平衡方程式（1-8）与达西定律（1-10）,对于给定的qV两式组合在一起用基于介质离散化的有限差分方法解出，离散化区域由两层包定的比流量四面体覆盖层组成数值计算是根据重力平衡方程的节点公式进行的这些公式即为牛顿力法的节点公式它是通过用孔隙压力、比流量矢量和孔隙压力梯度分别代替速度矢量、应力张量和应变率张量来得到初等常微分方程的解系是用与时间一致的显式和隐式离散化两种模式解出其主解描述如下1.3.1 空间导数的有限差分近似根据习惯，四面体节点的指局部由1到4编号的节点面n与编号为n的节点相对上标f与面f上的独立变量的值相联系线性孔隙压力变量和恒定流体密度被假设在四面体内部；根据应用高斯离散定理得到的节点孔压值得到压力头梯度的表达方式：4（P-PfXigi,j=---s（p—PfXiginjS（）（1-15）3V11这里n《）是垂直于面l外单位矢量，S是面的表面积，V是四面体的体积。

对于数值精确性，量Xi—x1,这里x1与四面体顶点之一的坐标一致，在（1-15）的压力头公式中可以用Xi来代替，而且并不影响压力头梯度的值，根据（1-15）,有下面公式：14,,,（p-PfXigi,j=——£p*lnfSO（1-16）3V11这里节点量p"定义如下：p*=p1_Pf（x：—X；gi（1-17）1.3.2 质量平衡方程的节点公式质量平衡方程（1-8）可以被描述为：qij+b：1-18）这里b冲」型_qV（1-19）Mft它等于瞬时“体力”，①用在力学节点公式应用这种相似，节点流量，Q；m3/sjn=1,4,等于四面体的比流量和体积源强度可以表述为：(1-20)etqvVndpdtQn=Qnm一这里Q；_ qi；i（n S（n）(1-21)n4M(1-22)作用于原则上，在每个节点的重力平衡方程节点形式都是根据需要建立的,边界的等节点流量（-Q；）与节点贡献流量（QW柏勺总和为零式（1-21）的比流量矢量部分与根据流动法则（1-10）得到的压力头梯度相联系依次地，压力头梯度部分可以根据式（1-16）四面体节点孔压来表述为了节省计算时间，在用式（1-21）在四边形每个节点n遇到的网格区获得的网格区贡献流量q；时，采用网格区公式。

局部网格区的矩阵集合把节点值q；与节点压力头p*n联系起来因为这些矩阵是对称的，36个组成部分被计算；这些计算在大变形模式下每十个计算步更新一次根据网格区矩阵的定义，我们得到：Q；=MnjP*j（1-23）这里P”是对于网格区的节点孔隙压力头的局部矢量依次地，总节点值Q；通过网格区贡献的累加来获得我们可以写出每个节点的总流量公式：Q；=CnjP刷（1-24）这里C]是总矩阵，B"]节点压力头的总矢量返回到我们前面的考虑因素，我得到：-zQen+£QW=0（1-25）这里，为了符号的简单化，符号z用来代表在节点n与该节点相遇的所有网格区的贡献流量之和一个网格区的贡献由网格区内所有的四面体流量之和）应用（1-20）,（1-9）,（1-25）,我们得到:dpndtapp'、QtnhJ(1-26)这里Q；是用式（1-24）和（1-17）得到节点孔隙压力的函数，Z；pp是体积源、边界通量和点源贡献的已知流量有如下形式：n£Qnpp=—£qvV+Qw（1-27）一4而£Qthm是节点热力贡献流量定义如下：工Qthm*,V[.£[V-P（1-28）dt_14j<4J方程（1-26）是节点n的重力平衡方程的节点形式；右边Q；+Q：pp+Q，是指不平衡流量。

它由两部分组成：流体不平衡流量QTn+Z：pp;不平衡热力流量£Q，当流体保持静止时（流体计算停止关闭），不平衡流消失，孔隙压力的改变只由力学或热变形引起对于独立的流体计算，不平衡热力流量等于零在FLAC3D中，Biot模数是节点属性，利用（1-22）我们得到：这里、. n、, mVn(1-29)Vn —(1-30)(1-28)代入(1-26),我们写出:nL?Q— Q：PP1(1-31)。