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

摘要：为解决水平面二维有限元方法计算减压井时井点附近奇异区问题，在李祖贻等针对特定网格划分推导出修正井水位法的成果基础上，推导适应一般网格划分的修正公式。通过特例分析计算，验证修正公式精度和适应性，减压井计算受单元尺寸、尺寸差异影响很小，流量和井点以外节点水头值满足精度要求。针对有限元计算时涉及的井阻力情况，给出实用的修正公式处理方法。

关键词：有限元 减压井 渗流 中断面法

20

(1华南理工大学 土木工程系，广东 广州510641；2广东省水利水电科学研究院，广东 广州510610)

　　减压井是堤防防渗加固的一种常用工程措施，需要定量分析它的出水量和减压效果。多数堤防强透水层水平成层分布，满足缓变渗流条件^[1]，可以用水平面二维有限元方法。井的公式大多在缓变渗流条件下推导出的，因此在水平面有限元计算中，可将减压井设为一个节点，以出水量和井内外水位差的协调，将井的公式与其结合，解决计算问题。因井点附近是奇异区，水头分布为对数旋转面，无法用有限个平面或低阶曲面拟合，必须作特殊处理。

　　李祖贻等^[1，2]以在井周划分为4个相同的等腰直角三角形单元的特殊情况(图1) ，推出修正井水位法及修正井周单元渗透系数法。设h_w为井水位，h₀为井节点计算水头，h_a为井周节点计算水头。井水位修正的要求是：给定修正量Δh，当h₀满足h₀=h_w+Δh时，计算得井的出水量Q和井周节点水头h_a与解析解相同。由此得：

式中：q=Q/T为单位厚度流量；a为节点间距。

由于该修正公式针对特定单元划分，应用时有一定局限性。修正井周单元渗透系数法的数学推导与修正井水位法相同，只在计算中处理方式不同，两者计算出水量差不多。但当井数多，单元尺寸与井间距相比不很小时，由于单元渗透系数修正使井后区域计算得水头偏低(回升水头偏小)。因此，后续讨论仅就修正井水位法进行。

减压井还受非完整井、井壁摩阻力和动力水头等影响，井水位与滤管外砂层的平均水头不同，分别用h′_w和h_w表示，两者之差是井出水量的函数，可由井的公式^[1，3]获得。

1 按等分圆周角划分单元时的修正公式

1.1 修正公式推导 图1中的4个三角形可看成在以井点为圆心，a为半径的圆周上4等分而分割成的单元。下面进一步讨论划分任意n等分的情况。图2 所示为等分n个单元后的其中一个。为便于讨论，将坐标系平移、旋转，使井点i落在(0，0) ，j点落在x轴上，这样不影响流量计算。该三角形两相邻边长为a，夹角φ＝２π/n，井点计算水头为h₀，j、m点水头h_a，单元流量q^e按中断面法计算^[1]：

式中：Δ为单元面积；b_i，b_j，b_m，c_i，c_j，c_m为单元节点对边向y及x轴投影长度，可在一般有限元书中找到。注意到图2中流量定义与井出水相反，按井点习惯出水为正，反号后由式(2)推得：

设流向井点的水均匀分布，总流量q与q^e的关系为q=nq^e=2π=q^e/φ，整理可得

按修正的要求，同一流量下，半径为r_w的井，距井点a处水头也为h_a，解析公式为

由式(5)与式(4)可得按等分圆周角划分单元时的修正公式

显然，式(1)是式(6)在φ＝π/2或n=4时的特例。

1.2 对修正公式的讨论 修正公式含有流量，受远处单元 影响情况需进一步论证。

　　讨论最简单的情形：单元分划成放射状，如图3所示。第1圈节点距井点a，节点数n，两节点对应圆心角φ，节点水头h_a。第m圈节点距井点ma，节点数mn，两节点对应圆心角φ/m ，坐标旋转为图3所示情况下，φ角对应区域的m+1个节点坐标为(macosi/φm，masini/φm(i=0，1，…，m)，节点水头h_m。该圈与m-1圈间共(2m-1)n个单元，其中mn个单元有两个节点在m圈上，(m-1)n个单元只有一个节点在m圈上。由中断面法可求得该圈单元向井流量近似值为

式中：l_m为m圈各三角形单元中断面长度与三角形高的比值累加再乘以单元尺寸a，l₁=antan(φ/2)，

由于通过各圈单元流向井点的流量相等，递推可得

式中：M为井的影响半径R对应的节点圈数，R＝Ma。设远方水头为h_R，将式(4)、式(6)代入式(8)，消去h₀，h_a得

　　解析解流量q′=2πk(h_R-h_w)/ln(Ma/r_w)。
　　作为对比，可计算同等条件下不作修正的流量q_u，这只须在式(8)中令h₀=h_w，并将式(4)代入消去h_a，得。
　　定义流量相对误差Δq=(q-q′)/q′，可对各种不同网格划分的计算流量进行比较。
　　表1列出R=1000rw时(相当于井径02m，影响半径100m)，不同单元尺寸计算流量的 相对误差Δq。

由表可见，未作修正时，流量计算误差很大，且受单元尺寸影响很大。经修正后，算得流量几乎不受尺寸影响。当φ由大变小，流量误差随之由负变正，当φ趋向于0，误差趋向于0.53%，不大于1%。误差最小值在φ＝π/2～π/4，即4～8等分圆周时。

总流量q与q_i的关系为∑q_i=q，按修正公式要求：

整理后得

式中：A_i为三角形面积；b_i为井点的对边边长。

　　按式(12)计算，当井周单元较均匀时，误差很小，但当井周节点分布不均时 ，会有偏差。其原因在于仅用一个修正量不可能使所有节点满足式(11)。为此增加对流向 井点水量的约束条件：设其均匀分布，q_i=qφ_i/2π。
定义单元修正量：Δh_i=h₀-h_w。让每个单元独立满足流量、节点水头条件， 求出各自对应的修正量，再选用适当的权函数加权平均求平均修正量。第i单元的修正量为

　　为保证修正后尽可能满足流量条件，权函数应与单元流量相关。由于单元流量与总流量关系是夹角φ_i与圆周角之比，因此取φi/2π为权函数，加权平均得

　　该式即为适应一般网格划分的修正公式。作为特例，当所有a_i=a，φ_i=φ时，式(13)变成式(6)。

3 修正公式在有限元程序中的处理方法

　　由式(13)或式(6)可见，Δh=h₀-h_w是q=Q/T的线性函数，可表示为一般形式：h₀-h_w=CQ或Q=(h0-hw)/C，式中C是与井周节点分布相关的常数。
　　在有限元计算时可根据涉及的井阻力情况选用不同处理方法：①完整井且不计井阻力。由Q=(h₀-h_w)/C，(注意有限元计算时以流入为正)只需在总系数矩阵对应于井节点的对角元素加1/C，右端项对应于井节点的元素加h_w/C。②非完整井不计井阻力。非完整井的井水位h′w与滤管外砂层的平均水头h_w不同，两者之差也是q的线性函数，具体形式在文献^[1，3]中不完全相同，但其形式都可写为hw-h′w=DQ，因此可按完整井的方法处理，只需将C变成C+D。③一般情况。计入井阻力时，流量Q与h_w-h′_w关系为非线性函数^[1，3]，但仍有Q＝(h₀-h_w)/C。在有限元计算时须用迭代方法求解。迭代过程为：先假定流量Q，通过有限元计算h₀。用h₀代入流量Q与h_w-h′_w和h₀的非线性函数解出新的流量Q′，用新流量调整旧流量，再进入计算，直至得到满意结果。如果是两层或多层强透水层，则对每一层分别进行迭代过程，各层流量Q与各层h_w-h′_w和h₀的关系一般为 Q的多元二次方程组^[3，4]。上述处理方法已编入“水平二维有限元程序”^[4]。

4 精度验算及应用情况简介

　　用一个简单算例可验证修正公式，设减压井为完整井，透水层厚T＝5m，渗透系数k=100m/d，井半径r_w=0.1m，影响半径R＝100m，远方水头与井口高程差h_R-hw=3m，不计井阻力。由解析公式可得井出水量Q＝1364.4m³/d，到井点距离为5m、10m、20m处降深分别为1.30m、1.00m、0.70m。

　　有限元计算区域为半径R＝100m的圆形，减压井设在圆心，三角形单元，按图6所示6种井周单元划分方法。
　　计算结果列于表2。表中可见经修正后，6种划分所得流量都与解析式算得相近，相对误差在-0.59%～0.73%之间，降深误差也很小。不修正的两种情况，流量误差为77%和50%，与表1所列相符。表中②④⑤单元尺寸相同，但划法不同，所得结果完全相同。⑥井周单元差异较大，但结果仍很满意。不同单元划分Δh差异较大，这与按等分圆周角的讨论结果相同，即h₀仅是一个与出水量相关的过渡数，并无实际物理意义。井点的真实水位仍是井水位 h′_w。

　　北江大堤大量应用减压井，是验证计算方法的主要工程对象。早期的研究计算只考虑井阻力，未采用井点水头修正方法，计算结果与实际有一定差异，流量、降深偏大，应用时只能按较大安全系数折减。例如1999年石角段莲藕塘险区加固，新布置56个减压井，计算得百年一遇洪水位15.3m，井附近水头比无井时下降3m，设计时只能按水头下降1.5m考虑
。工程现已完成，2001年3次小洪水，江水位11.1m时，测压管实测水头与1997年相同江水位时的测值相比，下降了0.5～0.8m。现采用井点水头修正方法复核^[4]，江水位11.1m时，井周水头，下降0.7m，与实测相近；江水位15.3m时，井附近水头比无井时下降1.5m。

　　目前计算程序都已采用井点水头修正方法，在北江大堤强透水地基堤段防渗加固处理分析、韩江潮州枢纽两岸浸没问题研究^[5]等项目中，得到大量应用，计算结果合理，符合工程规律。

5 结 语

　　(1)用水平面二维有限元方法计算减压井时，由于井点附近是奇异区，无法用有限个平面或低阶曲面拟合，必须作特殊处理。在李祖贻等针对特定网格划分的修正井水位法成果基础上，推导适应一般网格划分的修正公式。(2)以按任意等分圆周角划分单元为例，推导出相应的修正公式。通过特例分析计算，获得对修正公式误差量的认识：①流量计算几乎不受单元尺寸影响，受等分圆周角大小影响，流量误差由负变正，相对误差不大于1%。误差最小值在4～8等分圆周时；②井周外第二圈及以外节点水头值上与理论解吻合；③井点处计算结果h₀仅是一个与出水量相关的过渡数，并无实际物理意义。(3)用加权平均方法获得适应一般网格划分的修正公式。(4)针对有限元计算时涉及的井阻力情况，给出了实用的修正公式处理方法。(5)通过一个算例，用差异很大的6种单元划分方法，验证修正公式精度和适应性。(6)此方法已在实际工程中应用，满足工程要求。

参 考 文 献：
［1］毛昶熙主编.渗流计算分析与控制［M］.北京：水利电力出版社 ，1990.
［2］李祖贻，陈平.有限元法计算井的渗流对奇异点的处理［J］.水利学报，1984，(5)：41-50.
［3］陈雨孙.单井水力学［M］.北京：中国建筑工业出版社，1977.
［4］曹洪，张挺，陈小丹，郑存灼，羽海英.堤防中多层透水地基渗流分析计算方法研究［R］.广州：广东省水利水电科学研究院，2002.5.
［5］曹洪，尹小玲.潮州供水枢纽工程库区两岸浸没分析［R］.广州：华南理工大学，2002.6.