技术领域
[0001] 本
发明涉及
地下水计算技术领域，特别是指一种
地下水流量的计算方法。
背景技术
[0002] 局部均衡在很多地下水流数值计算过程中都是非常重要的流场性质。Zhang等指出，在两相流研究中，局部均衡发挥着重要作用。Sun和Liu在研究地下水流场时发现，在
流线绘制与地下水运动规律的研究中，局部均衡是保障计算结果
精度所必需的。Dogrul和Kadir分析了局部均衡性质在地下水资源管理中的重要性。Cass Miller在总结地下水资源的数值模拟时更是将局部均衡和含水介质的非均质性与不确定性、非承压水流模拟并称为地下水流模拟的三大挑战，并指出局部
质量守恒性质对于运移研究是至关重要的。
[0003] 基于三
角形的连续Galerkin有限单元法通常被认为是局部不均衡的数值方法，其主要体现为：在流量利用水
力坡度直接计算的前提下，每一个三角形单元不一定是均衡的。对于非稳定流问题，或者存在垂向上的补给或排泄时，每一个三角形单元往往是不均衡的。为了克服这一问题，国内外学者提出了控制体积有限单元法、混合有限单元法、不连续Galerkin有限单元法和局部均衡的Galerkin有限单元法。这些
算法都能够保持部分三角形单元的局部均衡性质。但是，这些方法的提出也带来了一些新的问题，从而影响了其在工程中的应用，例如：混合有限单元法和不连续Galerkin有限单元法计算效率远远低于连续Galerkin有限单元法等。
[0004] 解决连续Galerkin有限单元法局部均衡问题的另一思路是后处理算法，国内外学者提出了大量的通过修正流量的计算方法来求得满足三角形单元，甚至任意局部区域均衡的流量值。Cordes和Kinzelbach讨论了不存在单元垂向补排、不存在抽注水井、不考虑计算误差情况的稳定情况的局部均衡域及其相应的流量计算方法，但是因为其假设条件太过严格，不适应于实际工程项目。Hughes证明了连续Galerkin有限单元法任意均衡区域的均衡性，并提出了满足三角形单元均衡的流量计算方法，但是利用其算法计算得到了单元边界流量值是不唯一的。Berger和Howington讨论了一维连续Galerkin有限单元法的局部均衡域及其流量计算方法，提出了流量计算与连续Galerkin有限单元法的计算方程保持一致的重要性，但是他们并没有将其方法推广到二维和三维情况，属于理论上的探讨。划定基于三角形Galerkin有限单元法的局部均衡域并给出与之相应的流量计算方法是水文地质数值计算领域最迫切的课题之一，它严重限制着地下水流场
可视化、污染物运移、地下水管理以及
多相流等等相关研究的计算精度。
发明内容
[0005] 有鉴于此，本发明的目的在于提出一种地下水流量的计算方法，以提高地下水流量的计算精度。
[0006] 基于上述目的，本发明提供的地下水流量的计算方法包括以下步骤：
[0007] 1)利用三角形连续Galerkin有限单元法求解地下水流数值模型，形成原始三角形单元网格并获得地下水的
水头场数据；
[0008] 2)将所述原始三角形单元网格进行细化，从而构造局部均衡域；
[0009] 3)计算所述局部均衡域的相关流量值。
[0010] 在本发明的一些
实施例中，所述将原始三角形单元网格进行细化的步骤包括：
[0011] 连接原始三角形单元三条边的中点，从而将原始三角形单元分割为四个子三角形；
[0012] 在所述四个子三角形中，令完全处于原始三角形单元内部的子三角形为一类局部均衡域，令另围绕原始三角形网格某
顶点的子三角形单元所组成的区域为二类局部均衡域。
[0013] 在本发明的一些实施例中，所述相关流量值包括与原始三角形网格顶点相关的流量值、与原始三角形网格边相关的流量值、与原始三角形单元区域相关的流量值和局部均衡域的边界流量值。
[0014] 在本发明的一些实施例中，所述与原始三角形网格顶点相关的流量值包括定义在原始三角形顶点的井流量值；
[0015] 所述与原始三角形网格边相关的流量值包括定义在原始三角形边或原始三角形相邻顶点的边界流量值；
[0016] 所述与原始三角形单元区域相关的流量值包括单元储释水量及单元垂向补
排量；
[0017] 所述局部均衡域的边界流量值包括原始三角形边中点连接边的流量值。
[0018] 在本发明的一些实施例中，所述井流量的计算为：
[0019] Wp＝Wi
[0020] 其中，Wp为二类局部均衡域的井流量值，Wi为原始三角形网格顶点的井流量值。
[0021] 在本发明的一些实施例中，所述边界流量的计算包括：
[0022] 如果边界流量值直接定义在原始三角形网格边上，则：
[0023]
[0024] 其中，Bp为二类局部均衡域的边界流量值，Be为原始三角形网格边的边界流量值；
[0025] 如果边界流量值定义在原始三角形网格边的相邻顶点上并由单宽流量表示，则：
[0026]
[0027] 其中，Bp为二类局部均衡域的边界流量值，B1为二类局部均衡域所围绕原始三角形网格顶点的单宽流量值，B2为原始三角形网格边另一顶点的单宽流量值，l为该原始三角形网格边的长度。
[0028] 在本发明的一些实施例中，所述二类局部均衡域的储释水量的计算为：
[0029]
[0030] 其中，Sp为二类局部均衡域的储释水量，积分区间为二类局部均衡域，Se为原始三角形单元e的储释水系数；
[0031] 所述一类局部均衡域储释水量的计算为：
[0032]
[0033] 其中，Ss为一类局部均衡域的储释水量，积分区间为一类局部均衡域，Se为原始三角形单元e的储释水系数。
[0034] 在本发明的一些实施例中，所述二类局部均衡域的垂向补排量的计算为：
[0035]
[0036] 其中，Ep为二类局部均衡域的垂向补排量，εe为原始三角形单元的单元补给强度，Δe为原始三角形单元e的面积，e为二类局部均衡域所处的原始三角形单元；
[0037] 所述一类局部均衡域的垂向补排量计算为：
[0038]
[0039] 其中，Es为一类局部均衡域的垂向补排量，εe为原始三角形单元的单元补给强度，Δe为原始三角形单元e的面积，e为一类局部均衡域所处的原始三角形单元。
[0040] 在本发明的一些实施例中，所述原始三角形边中点连接边的流量计算为：
[0041]
[0042] 其中，a，b，c为原始三角形单元e的三个顶点并逆
时针排序，Aa表示边ab和边ac中点连接边的流量值，以流向顶点a为正； 为利用原始三角形水力坡度根据达西定律计算的原始三角形单元e内的流速向量； 表示边cb的向量；
[0043] 边ba和边bc中点连接边的流量计算为：
[0044]
[0045] 边ca和边cb中点连接边的流量计算为：
[0046]
[0047] 从上面所述可以看出，本发明提供的地下水流量的计算方法基于三角形连续Galerkin有限单元法求解的地下水流问题，通过连接原始三角形单元三条边的中点，将原始三角形单元分割成四个子三角形单元，划定了两种局部均衡域，即：完全处于原始三角形内部的子三角形单元和由围绕某原始三角形网格顶点的三角形子单元组成的区域，本发明在连续Galerkin有限单元法计算的地下水水头场的
基础上，提出了能够保证两种均衡域质量守恒的流量计算方法，从而能够提高基于三角形连续Galerkin有限单元法的地下水流速场及其相关计算的精度。
附图说明
[0048] 图1为本发明实施例的地下水流量的计算方法的流程示意图；
[0049] 图2为本发明实施例的研究区域模型示意图；
[0050] 图3为本发明实施例的平行于研究区域模型的南、北边界的断面示意图；
[0051] 图4为原始三角形单元细化示意图；
[0052] 图5为本发明实施例的一类局部均衡域构造示意图；
[0053] 图6为本发明实施例的二类局部均衡域构造示意图；
[0054] 图7为本发明实施例的一类局部均衡域和二类局部均衡域构造示意图；
[0055] 图8为本发明实施例的模型中选择的二类局部均衡域
位置及各边编号示意图；
[0056] 图9为本发明实施例的通过平行于模型南、北边界的断面流量与断面位置关系示意图。
具体实施方式
[0057] 为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。
[0058] 参照附图1，为本发明提供的地下水流量的计算方法的流程示意图。作为本发明的一个实施例，所述地下水流量的计算方法包括：
[0059] 步骤101：利用三角形连续Galerkin有限单元法求解地下水流数值模型，形成原始三角形单元网格并获得地下水的水头场数据。
[0060] 水头场是指研究区域内任意一点对应一个水头值。在本发明中是指任意三角形网格顶点对应一个水头值。研究区域内其他点可以通过插值得到，所以相当于知道了研究区域内任意点的水头值。
[0061] 研究区域的模型如图2所示，其研究区域为纵、横100m的方形，北边界为水头值为100m的定水头边界，南边界为水头值为50m的定水头边界，东、西边界均为隔水边界。研究区域中有两个弱透水区域，其位置如图2所示。研究区域的地下水流动规律为由北补给边界迳流至南排泄边界，并且大部分水流绕过弱透水区域流动。由于东、西边界为隔水边界，通过任意平行于南、北边界的断面流量值均相同。以研究区域西南角为坐标原点，图3为平行于研究区域模型的南、北边界的断面示意图，断面由其纵轴坐标标记。
[0062] 在该步骤中，利用三角形连续Galerkin有限单元法求解该研究区域的地下水流数值模型，从而形成原始三角形单元网格并获得地下水的水头场数据。
[0063] 步骤102：将所述原始三角形单元网格进行细化，从而构造局部均衡域。
[0064] 具体地，所述将原始三角形单元网格进行细化的步骤包括：
[0065] 如图4所示，连接原始三角形单元三条边的中点，从而将原始三角形单元分割为四个子三角形；在所述四个子三角形中，令完全处于原始三角形单元内部的子三角形为一类局部均衡域，参见图5，令围绕原始三角形网格某顶点的子三角形单元所组成的区域为二类局部均衡域，参见图6。
[0066] 如图7所示，其为本发明实施例的一类局部均衡域和二类局部均衡域构造示意图，图中的灰色区域为一类局部均衡域，白色区域为二类局部均衡域。在该步骤中，将原始三角形单元网格进行细化，得到一类局部均衡域和二类局部均衡域。因此，如图8所示，虚线0、1、2、3、4围成的区域即为所述研究区域模型中的二类局部均衡域位置。
[0067] 步骤103：计算所述局部均衡域的相关流量值，所述相关流量值包括与原始三角形网格顶点相关的流量值、与原始三角形网格边相关的流量值、与原始三角形单元区域相关的流量值和局部均衡域的边界流量值。
[0068] 其中，所述与原始三角形网格顶点相关的流量值包括定义在原始三角形顶点的井流量值；
[0069] 所述与原始三角形网格边相关的流量值包括定义在原始三角形边或原始三角形相邻顶点的边界流量值；
[0070] 所述与原始三角形单元区域相关的流量值包括单元储释水量及单元垂向补排量；
[0071] 所述局部均衡域的边界流量值包括原始三角形边中点连接边的流量值。
[0072] 具体地，所述相关流量值的计算包括以下步骤：
[0073] 1)将与原始三角形网格顶点相关的流量完全分配到围绕该顶点的二类局部均衡域，例如井流量的计算为：
[0074] Wp＝Wi
[0075] 其中，Wp为二类局部均衡域的井流量值，Wi为原始三角形网格顶点的井流量值。
[0076] 2)将与原始三角形网格边相关的流量平均分配到相邻的两个二类局部均衡域，例如边界流量的计算。
[0077] 如果边界流量值直接定义在原始三角形网格边上，则：
[0078]
[0079] 其中，Bp为二类局部均衡域的边界流量值，Be为原始三角形网格边的边界流量值。
[0080] 如果边界流量值定义在原始三角形网格边的相邻顶点上并由单宽流量表示，则：
[0081]
[0082] 其中，Bp为二类局部均衡域的边界流量值，B1为二类局部均衡域所围绕原始三角形网格顶点的单宽流量值，B2为原始三角形网格边另一顶点的单宽流量值，l为该原始三角形网格边的长度。
[0083] 3)将与原始三角形区域相关的流量值平均分配到整个原始三角形单元，例如单元储释水量的计算为：
[0084]
[0085] 其中，Sp为二类局部均衡域的储释水量，积分区间为二类局部均衡域，Se为原始三角形单元e的储释水系数。
[0086] 一类局部均衡域储释水量的计算与此类似，只是积分区间为一类局部均衡域，即：
[0087]
[0088] 其中，Ss为一类局部均衡域的储释水量，积分区间为一类局部均衡域，Se为原始三角形单元e的储释水系数。
[0089] 与原始三角形区域相关的另一流量值——垂向补排量(降雨入渗、
蒸发、越流补排等)的计算为：
[0090]
[0091] 其中，Ep为二类局部均衡域的垂向补排量，εe为原始三角形单元的单元补给强度，Δe为原始三角形单元e的面积，e为二类局部均衡域所处的原始三角形单元。
[0092] 一类局部均衡域的垂向补排量计算为：
[0093]
[0094] 其中，Es为一类局部均衡域的垂向补排量，εe为原始三角形单元的单元补给强度，Δe为原始三角形单元e的面积，e为一类局部均衡域所处的原始三角形单元。
[0095] (4)所述局部均衡域的边界流量，即原始三角形边中点连接边的流量计算为：
[0096]
[0097] 其中，a，b，c为原始三角形单元e的三个顶点并逆时针排序，Aa表示边ab和边ac中点连接边的流量值，以流向顶点a为正； 为利用原始三角形水力坡度根据达西定律计算的原始三角形单元e内的流速向量； 表示边cb的向量。
[0098] 其他中点连接边的流量计算与此类似，边ba和边bc中点连接边的流量计算为：
[0099]
[0100] 边ca和边cb中点连接边的流量计算为：
[0101]
[0102] 采用本发明提供的地下水流量的计算方法计算得到的通过任意平行于南、北边界断面的流量值与断面的坐标值关系如图9中方
块线所示，由地下水数值模拟
软件Feflow计算的断面流量与断面坐标值的关系如图9中三角形线所示。这一对比表明，本发明提供的计算方法大大提高了地下水流量计算的精度。
[0103] 为了验证两类局部均衡区的均衡性质，从理论上检验了一类均衡区的均衡性；任意选择围绕某一顶点的二类均衡区，通过实例检验了二类均衡区的均衡性，即：
[0104] 1.由于一类局部均衡域没有井流量和研究区域边界流量，其均衡性质可以通过以下公式检验：
[0105]
[0106] 2.二类局部均衡域的均衡性可以选择以上模型中围绕某一原始三角形网格顶点的二类局部均衡域，如图8所示。该模型中，没有考虑垂向补排量、储释水量、井流量，选择的顶点处于研究区域内部，所以也不涉及研究区域的边界补给量，只与局部均衡域的边界流量有关，其通过二类局部均衡域边界各边的流量如表1所示，其各边编号(0、1、2、3、4)如图8所示。
[0107] 表1
[0108]
[0109] 由此可见，本发明提供的地下水流量的计算方法基于三角形连续Galerkin有限单元法求解的地下水流问题，通过连接原始三角形单元三条边的中点，将原始三角形单元分割成四个子三角形单元，划定了两种局部均衡域，即：完全处于原始三角形内部的子三角形单元和由围绕某原始三角形网格顶点的三角形子单元组成的区域，本发明在连续Galerkin有限单元法计算的地下水水头场的基础上，提出了能够保证两种均衡域质量守恒的流量计算方法，从而能够提高基于三角形连续Galerkin有限单元法的地下水流速场及其相关计算的精度。
[0110] 所属领域的普通技术人员应当理解：以上任何实施例的讨论仅为示例性的，并非旨在暗示本公开的范围(包括
权利要求)被限于这些例子；在本发明的思路下，以上实施例或者不同实施例中的技术特征之间也可以进行组合，并存在如上所述的本发明的不同方面的许多其它变化，为了简明它们没有在细节中提供。因此，凡在本发明的精神和原则之内，所做的任何省略、
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
