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

程度面渗流有限元计算中减压井点处理摘要：为解决程度面二维有限元方法计算减压井时井点附近奇异区问题，在李祖贻等针对特定网格划分推导出修正井水位法的成果根底上，推导适应一般网格划分的修正公式。通过特例分析计算，验证修正公式精度和适应性，减压井计算受单元尺寸、尺寸差异影响很小，流量和井点以外节点水头值满足精度要求。针对有限元计算时涉及的井阻力情况，给出实用的修正公式处理方法。关键词：有限元减压井渗流中断面法20(1华南理工大学土木工程系，广东广州510641；2广东省水利水电科学研究院，广东广州510610)减压井是堤防防渗加固的一种常用工程措施，需要定量分析它的出水量和减压效果。多数堤防强透水层程度成层分布，满足缓变渗流条件[1]，可以用程度面二维有限元方法。井的公式大多在缓变渗流条件下推导出的，因此在程度面有限元计算中，可将减压井设为一个节点，以出水量和井内外水位差的协调，将井的公式与其结合，解决计算问题。因井点附近是奇异区，水头分布为对数旋转面，无法用有限个平面或低阶曲面拟合，必须作特殊处理。李祖贻等[1,2]以在井周划分为4个一样的等腰直角三角形单元的特殊情况(图1)，推出修正井水位法及修正井周单元浸透系数法。设h为井水位，h0为井节点计算水头，ha为井周节点计算水头。井水位修正的要求是：给定修正量Δh,当h0满足h0=h+Δh时，计算得井的出水量Q和井周节点水头ha与解析解一样。由此得：(1)式中：q=Q/T为单位厚度流量；a为节点间距。由于该修正公式针对特定单元划分，应用时有一定局限性。修正井周单元浸透系数法的数学推导与修正井水位法一样，只在计算中处理方式不同，两者计算出水量差不多。但当井数多，单元尺寸与井间距相比不很小时，由于单元浸透系数修正使井后区域计算得水头偏低(上升水头偏小)。因此，后续讨论仅就修正井水位法进展。减压井还受非完好井、井壁摩阻力和动力水头等影响，井水位与滤管外砂层的平均水头不同，分别用h′和h表示，两者之差是井出水量的函数，可由井的公式[1,3]获得。1按等分圆周角划分单元时的修正公式1.1修正公式推导图1中的4个三角形可看成在以井点为圆心，a为半径的圆周上4等分而分割成的单元。下面进一步讨论划分任意n等分的情况。图2所示为等分n个单元后的其中一个。为便于讨论，将坐标系平移、旋转，使井点i落在(0,0)，j点落在x轴上，这样不影响流量计算。该三角形两相邻边长为a，夹角φ＝２π/n，井点计算水头为h0，j、点水头ha，单元流量qe按中断面法计算[1]：(2)式中：Δ为单元面积；bi,bj,b,i,j,为单元节点对边向y及x轴投影长度，可在一般有限元书中找到。注意到图2中流量定义与井出水相反，按井点习惯出水为正，反号后由式(2)推得：(3)设流向井点的水均匀分布，总流量q与qe的关系为q=nqe=2π=qe/φ，整理可得(4)按修正的要求，同一流量下，半径为r的井，距井点a处水头也为ha，解析公式为(5)由式(5)与式(4)可得按等分圆周角划分单元时的修正公式(6)显然，式(1)是式(6)在φ＝π/2或n=4时的特例。1.2对修正公式的讨论修正公式含有流量，受远处单元影响情况需进一步论证。讨论最简单的情形：单元分划成放射状，如图3所示。第1圈节点距井点a，节点数n，两节点对应圆心角φ，节点水头ha。第圈节点距井点a，节点数n，两节点对应圆心角φ/，坐标旋转为图3所示情况下，φ角对应区域的+1个节点坐标为(asi/φ,asini/φ(i=0,1,…,)，节点水头h。该圈与-1圈间共(2-1)n个单元，其中n个单元有两个节点在圈上，(-1)n个单元只有一个节点在圈上。由中断面法可求得该圈单元向井流量近似值为q=kl(h-h-1)/a(6)式中：l为圈各三角形单元中断面长度与三角形高的比值累加再乘以单元尺寸a，l1=antan(φ/2)，(7)由于通过各圈单元流向井点的流量相等，递推可得(=1,2,…,)(8)式中：为井的影响半径R对应的节点圈数，R＝a。设远方水头为hR，将式(4)、式(6)代入式(8)，消去h0，ha得(9)表1不同网格划分的计算流量相对误差比拟(%)a/rπ/2π/4π/8φ→0修正不修正修正不修正修正不修正修正不修正2050100-0.404-0.403-0.39925.350.377.00.3530.3520.35019.441.965.40.4830.4820.47818.039.962.70.5310.5300.52517.639.361.9由表可见，未作修正时，流量计算误差很大，且受单元尺寸影响很大。经修正后，算得流量几乎不受尺寸影响。当φ由大变小，流量误差随之由负变正，当φ趋向于0，误差趋向于0.53%，不大于1%。误差最小值在φ＝π/2～π/4，即4～8等分圆周时。转贴于论文联盟.ll.qi=k/2aiai+1sinφi[(hai-h0)ai+1(ai+1-aisφi)+(hai+1-h0)ai(ai-ai+1sφi)](10)总流量q与qi的关系为∑qi=q，按修正公式要求：Δh=h0-h(i=1,2，……，)(11)整理后得(12)式中：Ai为三角形面积；bi为井点的对边边长。按式(12)计算，当井周单元较均匀时，误差很小，但当井周节点分布不均时，会有偏向。其原因在于仅用一个修正量不可能使所有节点满足式(11)。为此增加对流向井点水量的约束条件：设其均匀分布，qi=qφi/2π。定义单元修正量：Δhi=h0-h。让每个单元独立满足流量、节点水头条件，求出各自对应的修正量，再选用适当的权函数加权平均求平均修正量。第i单元的修正量为为保证修正后尽可能满足流量条件，权函数应与单元流量相关。由于单元流量与总流量关系是夹角φi与圆周角之比，因此取φi/2π为权函数，加权平均得(13)表2不同单元划分时有限元计算结果与解析解比拟解析式修正公式不修正公式①②③④⑤⑥①⑤Q/(3·d-1)Δh/1364.401356.51.311356.41.011357.91.211356.41.011356.41.011374.41.172407.80.002045.90.00降深r=5r=10r=201.301.000.701.010.701.311.000.691.311.010.691.311.000.691.311.000.691.341.020.701.801.251.981.511.05北江大堤大量应用减压井，是验证计算方法的主要工程对象。早期的研究计算只考虑井阻力，未采用井点水头修正方法，计算结果与实际有一定差异，流量、降深偏大，应用时只能按较大平安系数折减。例如1999年石角段莲藕塘险区加固，新布置56个减压井，计算得百年一遇洪水位15.3，井附近水头比无井时下降3，设计时只能按水头下降1.5考虑。工程现已完成，2001年3次小洪水，江水位11.1时，测压管实测水头与1997年一样江水位时的测值相比，下降了0.5～0.8。现采用井点水头修正方法复核[4]，江水位11.1时，井周水头，下降0.7，与实测相近；江水位15.3时，井附近水头比无井时下降1.5。目前计算程序都已采用井点水头修正方法，在北江大堤强透水地基堤段防渗加固处理分析、韩江潮州枢纽两岸浸没问题研究[5]等工程中，得到大量应用，计算结果合理，符合工程规律。5结语(1)用程度面二维有限元方法计算减压井时，由于井点附近是奇异区，无法用有限个平面或低阶曲面拟合，必须作特殊处理。在李祖贻等针对特定网格划分的修正井水位法成果根底上，推导适应一般网格划分的修正公式。(2)以按任意等分圆周角划分单元为例，推导出相应的修正公式。通过特例分析计算，获得对修正公式误差量的认识：①流量计算几乎不受单元尺寸影响，受等分圆周角大小影响，流量误差由负变正，相对误差不大于1%。误差最小值在4～8等分圆周时；②井周外第二圈及以外节点水头值上与理论解吻合；③井点处计算结果h0仅是一个与出水量相关的过渡数，并无实际物理意义。(3)用加权平均方法获得适应一般网格划分的修正公式。(4)针对有限元计算时涉及的井阻力情况，给出了实用的修正公式处理方法。(5)通过一个算例，用差异很大的6种单元划分方法，验证修正公式精度和适应性。(6)此方法已在实际工程中应用，满足工程要求。参考文献：［1］毛昶熙主编.渗流计算分析与控制［］.北京：水利电力出版社，1990.［2］李祖贻，陈平.有限元法计算井