有限元法在水文求参及矿坑涌水量预测中应用.docVIP

有限元法在水文求参及矿坑涌水量预测中应用 　　摘 要：本文简述有限元法数学原理及其推导，通过运用某矿山大口径群孔抽水试验数据，划定计算水文区域的界线及性质，确定水文工程结点（井、观测孔）及内插虚拟结点的水文参数（初始水位、单元性质及顶底板标高），运用ALGOL60编制优选参数程序及矿坑涌水量预测程序进行数值模拟，逐个优选参数，确定参数最终值，预测矿坑涌水量。 　　关键词：有限元法 边界条件 单元剖分 参数优选 矿坑涌水量预测 　　 　　一、数学模型 　　有限元法（finite elementmethod）是1956年研究飞机结构的计算时提出来的。1960年，Clough进一步处理了平面弹性问题，并第一次提出了“有限单元法”，使人们认识到它的功效。我国著名力学家、教育家徐芝纶院士（河海大学教授）首次将有限元法引入我国，对它的应用起了很大的推动作用。自从1969年以来，某些学者在流体力学中应用加权余数法中的迦辽金法（Galerkin）或最小二乘法等同样获得了有限元方程。其数学原理及其推导简述如下： 　　描述非均质各向异性二维不稳定流的基本方程是： 　　式中：H―承压水头或潜水面标高（米）； 　　TX―X方向的导水系数（米2/日）； 　　TY―Y方向的导水系数（米2/日）； 　　Q―入渗或蒸发量（米/日）； 　　μ―承压含水层的储水系数或潜水含水层的给水度（无量纲） 　　在一类边界上有； 　　在二类边界上有； 　　f，g是已知函数。 　　（x，y）―表示域内的点；（Xb，Yb）―表示边界的点； 　　Tb―表示弱透水边界的导水系数；―表示水头沿边界法线方向的导数； 　　1、根据变分原理，解上述非稳定流的微分方程等价于求下述泛函的变分问题。 　　此方程下一时刻的解： 　　对时间导数用中心差分来近似，重复这一过程可求出水头H在空间和时间上的分布。 　　为解上述变分问题把域D用三角形单元离散化，h用三角线性插值离散化为结点上的值h1，h2，…hm，于是泛函I=I（h1，h2，…hm），求极小值问题就归结为解线代数方程组（i=1，2，…m）的解。 　　2、方程组的系数矩阵用变带宽存储，方程组的解法用直接法把系数矩阵分为上三角阵和下三角阵。为系数矩阵，是对角线矩阵，的元素。 　　令{g}=[L]{}[L1]=[L]T[L0] 　　为下三角阵，（i=1,2…n） 　　先向前代入确定{g}，再向后回代确定{} 　　3、解正问题的方法 　　先给出一组参数初值然后解正，比较计算水头与实测水头是否一致，若差别较大，修改参数后重新解正，使计算水头与实测水头的差达到足够小。 　　考虑目标函数 　　其中：R1，R2……是待求参数；t-表示比较时刻的号码； 　　r-表示观测点的号码；Wt，r-权因子； 　　与―分别表示第t个比较时刻第r个比较点上的计算水位与观测水位。 　　目标函数E是计算水头与实测水头拟合程度和度量。 　　在优选参数中，用试估法和斐波那西分数法相结合。即人工和机器相结合的方法逐个优选参数。 　　二、流程 　　根据矿区水文地质条件确定各参数的初值及上下限代入方程解正，有观测孔的结点将某一时刻的计算水位与该时刻的实测水位进行比较，输出计算水位实测水位之差及误差平方和，然后对每一参数逐一进行修改。修改是根据每个观测孔的计算水位与实测水位之差并考虑参数的修改情况对水位分布的影响，再结合对水文地质条件的认识来进行分析，找出变化差异的原因，从而修改不均质分区，参数变化范围条件。这样多次反复就加深了对水文地质条件的认识，使对参数的取值、边界和“天窗”位置的选择不均质分区等更趋合理。流程如右图。参数优选程序及水量预测程序均用ALGOL60编制。 　　三、计算方案及数据整理 　　1、边界条件 　　某矿山通过群孔抽水试验与恢复水位试验，暴露了矿区供水与隔水边界条件。根据抽水试验等水位线下降漏斗可反映出矿区内部水位呈等幅下降，透水性强，南面与西面水位下降曲线梯度逐渐变陡，透水性弱，东与北面各观测孔均未受影响，为隔水边界。根据恢复水位试验分析，矿区的西与西南面水位恢复较快，显示出供水边界现象，南面与东面水位恢复迟缓，由东面隔水和南部透水性相对弱造成。矿区内经222天后，水位才恢复到抽水前的水位，恢复速度十分缓慢，反映出矿区内存在一个弱透水边界现象。 　　2、“天窗”补给 　　在大口径群孔抽水试验中，对第四系各含水层进行了观测。在天窗部位的QK003号孔及CK706号孔，在坡积层的观测中得到了系统的越流水位。计算的天窗注入量用CK706孔坡积层的观测水位与基岩动水位的差与“天窗”面积及越流系数的积来计算。 　　3、单元剖分 　　本次计算在拟合过程中分三个时段：第一时段是抽水0-3天的拟合