水平面渗流有限元计算中减压井点处理20131030310467631.doc

水平面渗流有限元计算中减压井点处理 1??按等分圆周角划分单元时的修正公式1.1??修正公式推导??图1中的4个三角形可看成在以井点为圆心，a为半径的圆周上4等分而分割成的单元。下面进一步讨论划分任意n等分的情况。图2?所示为等分n个单元后的其中一个。为便于讨论，将坐标系平移、旋转，使井点i落在(0,0)?，j点落在x轴上，这样不影响流量计算。该三角形两相邻边长为a，夹角φ＝２π/n，井点计算水头为h0，j、m点水头ha，单元流量qe按中断面法计算[1]： (2) 式中：Δ为单元面积；bi,bj,bm,ci,cj,cm为单元节点对边向y及x轴投影长度，可在一般有限元书中找到。注意到图2中流量定义与井出水相反，按井点习惯出水为正，反号后由式(2)推得： (3) 设流向井点的水均匀分布，总流量q与qe的关系为q=nqe=2π=qe/φ，整理可得 (4) 按修正的要求，同一流量下，半径为rw的井，距井点a处水头也为ha，解析公式为 (5) 由式(5)与式(4)可得按等分圆周角划分单元时的修正公式 (6) 显然，式(1)是式(6)在φ＝π/2或n=4时的特例。1.2??对修正公式的讨论??修正公式含有流量，受远处单元 影响情况需进一步论证。????讨论最简单的情形：单元分划成放射状，如图3所示。第1圈节点距井点a，节点数n，两节点对应圆心角φ，节点水头ha。第m圈节点距井点ma，节点数mn，两节点对应圆心角φ/m，坐标旋转为图3所示情况下，φ角对应区域的m+1个节点坐标为(macosi/φm,masini/φm(i=0,1,…,m)，节点水头hm。该圈与m-1圈间共(2m-1)n个单元，其中mn个单元有两个节点在m圈上，(m-1)n个单元只有一个节点在m圈上。由中断面法可求得该圈单元向井流量近似值为 qm=klm(hm-hm-1)/a (6) 式中：lm为m圈各三角形单元中断面长度与三角形高的比值累加再乘以单元尺寸a，l1=antan(φ/2)， (7) 由于通过各圈单元流向井点的流量相等，递推可得 ???? (m=1,2,…,M) (8) 式中：M为井的影响半径R对应的节点圈数，R＝Ma。设远方水头为hR，将式(4)、式(6)代入式(8)，消去h0，ha得 (9) ????解析解流量q′=2πk(hR-hw)/ln(Ma/rw)。????作为对比，可计算同等条件下不作修正的流量qu，这只须在式(8)中令h0=hw，并将式(4)代入消去ha，得。????定义流量相对误差Δq=(q-q′)/q′，可对各种不同网格划分的计算流量进行比较。????表1列出R=1000rw时(相当于井径02m，影响半径100m)，不同单元尺寸计算流量的 相对误差Δq。 表1??不同网格划分的计算流量相对误差比较(%) a/rw π/2 π/4 π/8 φ→0 修正 不修正 修正 不修正 修正 不修正 修正 不修正 2050100 -0.404-0.403-0.399 25.350.377.0 0.3530.3520.350 19.441.965.4 0.4830.4820.478 18.039.962.7 0.5310.5300.525 17.639.361.9 由表可见，未作修正时，流量计算误差很大，且受单元尺寸影响很大。经修正后，算得流量几乎不受尺寸影响。当φ由大变小，流量误差随之由负变正，当φ趋向于0，误差趋向于0.53%，不大于1%。误差最小值在φ＝π/2～π/4，即4～8等分圆周时。 图4为井点附近节点水头与理论解(对数曲线)比较(R＝1000rw,a=50rw,φ=π/4)。图中可见修正后的井周外第二圈节点水头与理论曲线吻合很好。井点处计算结果h0与理论曲线相差较大。由公式的定义即知，h0仅是一个与出水量相关的过渡数，并无实际物理意义。未经修正的计算结果则表现为井点准确，井点外节点水头明显偏低。因此，井点修正的真正意义在于，以放弃井点水头准确而换取出水量及井周节点水头的准确。2??一般网格划分的修正公式????对更一般的情况，井点附近剖分成m个三角形单元，第i单元位于井点处夹角为φi，所有单元夹角组成圆周角∑φi=2π。