流体的质量守恒方程.docxVIP

流体的质量守恒方程 这是流量修正的恒定方程。若单元体积内不发生源汇,即W=0,则: divJ=0 (10) 这是整个流体体系的质量守恒方程,即连续性方程。 各向同性介质的渗透系数是一个表示多孔介质输送流体能力的标量。此时不可压缩流体的三维运动的Darcy定律为: {qx=ΚΙx=-Κ?Η?xqy=ΚΙy=-Κ?Η?yqz=ΚΙz=-Κ?Η?z(11) 对于各向异性介质,渗透系数是一个张量,三维情况下: [Κ]=[ΚxxΚxyΚxzΚyxΚyyΚyzΚzxΚzyΚzz](12) Darcy定律的表达为: q=[Κ]Ι(13) 根据连续性方程式(10)及Darey定律,可以得出:divJ=div(gradH)=0 (14) 即:ΔΗ=?2Η=?2Η?x2+?2Η?y2+?2Η?z2=0(15) 式(15)为Laplace方程,亦称不可压缩流体运动的微分方程,它描述均质各向同性的不变形介质中不可压缩流体渗流场中水头的分布H=H(x,y,z,t)。它对于稳定流动和不可压缩流体的非稳定流动都适用。满足Laplace方程的非稳定流动是由于随时间变化和边界条件引起的。 任何偏微分方程都有无穷多个解。要从众多可能解中得到所研究问题的解,除微分方程外,还要提供求解条件,包括初始条件和边界条件。 3.1.1 初始条件 是指渗流过程开始的某一特定时刻研究区域上所有点必须满足的初始条件。 例如必须给出初始时刻区域内的水头分布。设水头H=H(x,y,z,t),初始条件的表达式一般为 H(x,y,z,t)=H0(x,y,z),(x,y,z∈Ω) (16) 式(16)中,H0—已知函数;Ω—渗流区域。 3.1.2 边界条件 ① 第一类边界条件为已知水头的边界条件,在边界的所有点上水头是给定的,例如: H(x,y,z)=Hb(x,y,z), (x,y,z)∈s(17) H(x,y,z,t)=Hb(x,y,z,t),(x,y,z)∈s(18) 式中,s—三维区域的边界曲面。 ②第二类边界条件为已知通量的边界条件,即垂直于边界面的流量给定的,一般有: qn=→q?→n=-qb(x,y,z)?(x,y,z)∈s(19) 或qn=→q?→n=-qb(x,y,z,t)?(x,y,z)∈s(20) 式中,qb—边界面上沿法线方向的单位面积入量;n—边界外法线的单位矢量。 ③第三类边界条件为水头和水头的法向导数的组合在边界上为已知,即: ?Η?n+λ(x,y,z)Η=f(x,y,z)?(x,y,z)∈s(21) 式(21)中,λ、f—均为已知函数,λ称为变换系数。 3.2 二维有限元方法 首先,根据渗流运动微分方程与边界条件约束建立二维数学模型: {??x(Τ?Η?x)+??y(Τ?Η?y)=0,在Ω内Η=Ηb(x,y),在Γ1上Τ?Η?xcos(n,x)+Τ?Η?ycos(n,y)=q,在Γ2上(22) 式(22)中,Ω—计算区域;Γ1—第一类边界;Γ2—第二类边界;Hb—第一类边界上的已知水头;n—第二类边界的外法线方向;q—第二类边界上法向单宽流量,流入为正,流出为负。 采用伽辽金加权剩余法,设ˉΗ(x,y)为区域Ω上地下水头的试探解。 ˉΗ(x,y)=ne∑L=1ΝLΗL(23) 式(23)中,NL—在单元 e上分片定义的基函数;HL—单元节点上的水头值;ne—单元结点数目。令基函数为权函数,代入试函数ˉΗ。令权剩余为零,则权剩余方程为: ∫ΩΝL(L(ˉΗ-f)dΩ=Μ∑e=1∫ΩeΝL(L(ˉΗ)-f)dΩ=0(24) 式(24)把权剩余在整个区域上的积分,化为在各个单元上的积分,然后求和,这样得到一个线性方程组,求解此线性方程组,便可得到各节点的水头值。 3.2.1 基函数 三角形单元节点i,j,m按逆时针编号,则: ˉΗ(x,y)=12A[(ai+bix+ciy)Ηi+(aj+bjx+cjy)Ηj+(am+bmx+cmy)Ηm](ΝiΗi+ΝjΗj+ΝmΗm)(25) 于是,三角形单元基函数表达式为: NL=(aL+bLx+cLy)/2A, (L=i,j,m) (26) 3.2.2 渗透矩阵 由上面给出的基函数,故试函数的表达式为: ˉΗ=ne∑i=1ΝLΗL=ΝiΗi+ΝjΗj+ΝmΗm(27) 权剩余方程为: ?ΩΝL[??x(Τ?ˉΗ?x)+??y(Τ?ˉΗ?y)]dxdy=0(L=1,2,?,Ν)(28) 对式(28)进行分部积分,并应用格林公式,得: ∫Γ2ΝLqds-?DL(Τ?ˉΗ?x?ΝL?x+Τ?ˉΗ?y?ΝL?y)dxdy=0(29) 对于每一个单元,有: QLe=?Ω(Τe?Ηˉ?x?ΝL?x+Τe?Ηˉ?y?ΝL?y)dxdy-∫Γ2eΝLds(30) 根据式(25)、式(