地下水流场的三维流线可视化方法
【技术领域】
[0001] 本发明涉及地下水计算技术领域，特别是指一种地下水流场的三维流线可视化方 法。
【背景技术】
[0002] 由于地下水赋存于地面以下地层空隙中，人类不能够直接观察，只能通过水文地 质勘探与地下水动态监测才能揭示其赋存条件与运动规律。受勘探经费的影响，水文地质 勘探孔的布置密度有限，对区域水文地质条件的认识存在较大的误差。然而，不论是研究地 下水的赋存规律、可持续开采地下水，还是研究污染物质在地下水系统中的运移，或是地下 工程开挖作业时的突水安全研究，都需要对地下水运动所形成的流场有充分、清晰的认识。
[0003] 如何将地下水赋存运动规律和动态特征实现可视化、将地下水运动所形成的流场 完整连续且分布均匀地展现在人们眼前，使研究区水文地质条件得到较全面准确的揭示、 为水文地质研究提供科学依据已成为紧迫的研究课题之一，也是水文地质学的一国际性难 题。
[0004] 对于流场的描述有很多种方法，其中以流线为代表的线表示方法能够很好地揭示 流场的走向，是实现流带、流管、流面表示方法以及粒子动画的基础，对这种方法的研究和 使用是一个热点。基于不同的地下水流数值模拟方法，地下水流场的可视化技术也不相同。 在地下水流数值模拟方法上，有限单元法因为其灵活性与易操作性而被广泛应用。但是，基 于有限单元法计算的地下水流速场精度存在着单元边界垂向流速不连续、介质边界切向流 速加权不连续以及局部区域不均衡等三个问题，严重限制了基于有限单元法的流场可视化 技术的发展。
【发明内容】
[0005] 有鉴于此，本发明的目的在于提出一种地下水流场的三维流线可视化方法，以提 高对地下水流场认识的准确度和清晰度。
[0006] 基于上述目的，本发明提供的地下水流场的三维流线可视化方法包括以下步骤：
[0007] 包括以下步骤：
[0008] 1)利用三棱柱连续Galerkin有限单元法求解地下水流数值模型，形成原始三棱 柱单元网格并获得地下水的水头场数据；
[0009] 2)将所述原始三棱柱单元网格进行细化，从而构造局部均衡域；
[0010] 3)计算所述局部均衡域的边界流量值；
[0011] 4)计算子三棱柱单元的边界流量值；
[0012] 5)计算子三棱柱单元内部的流速场；
[0013] 6)利用质点追踪算法获得地下水流场的三维流线。
[0014] 在本发明的一些实施例中，所述将原始三棱柱单元网格进行细化的步骤包括：
[0015] 连接原始三棱柱单元的剖面三角形三条边的中点，从而将该剖面三角形分割为面 积相等的四个子三角形，将该剖面三角形三条边的中点的连线沿着垂向连续地无限延展， 从而将该原始三棱柱单元分割为四个中间三棱柱单元；然后，连接该原始三棱柱单元的三 条垂向边的中点，使其构成一个中线三角形，该中线三角形所在的平面将四个中间三棱柱 单元分别分割为两个子三棱柱单元；
[0016] 在所述八个子三角形中，记a、b、c分别为原始三棱柱单元上层的三个顶点，a'、 b'、c'分别为原始三棱柱单元下层的三个顶点，围绕顶点a的子三棱柱单元记为0号，围绕 顶点b的子三棱柱单元记为1号，围绕顶点c的子三棱柱单元记为2号，围绕顶点a'的子 三棱柱单元记为4号，围绕顶点b'的子三棱柱单元记为5号，围绕顶点c'的子三棱柱单元 记为6号，不围绕原始三棱柱顶点的上层子三棱柱单元记为3号，不围绕原始三棱柱顶点的 下层子三棱柱单元记为7号，令3号子三棱柱单元与相邻的原始三棱柱单元的7号子三棱 柱单元所组成的区域为一类局部均衡域，令围绕原始三棱柱单元顶点的子三棱柱单元所组 成的区域为二类局部均衡域。
[0017] 在本发明的一些实施例中，所述局部均衡域的边界流量值的计算方法包括：
[0018] 0号子三棱柱单元与4号子三棱柱单元间相邻面的流量计算为：
[0019]
[0020] 其中，Ava为0号子三棱柱单元与4号子三棱柱单元间相邻面的流量值，以0号子 三棱柱单元流向4号子三棱柱单元为正。△e为原始三棱柱单元剖面三角形的面积，△z为 原始三棱柱单元的高度，1^为原始三棱柱单元的垂向渗透系数，113、11 1)、11(：、113'、111)'和11(：'分 别为原始三棱柱六个顶点的水头值；
[0021] 1号子三棱柱单元与5号子三棱柱单元间相邻面的流量计算为：
[0022]
[0023] 2号子三棱柱单元与6号子三棱柱单元间相邻面的流量计算为：
[0024]
[0025] 3号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0026]
[0027] 0号子三棱柱单元与3号子三棱柱单元间相邻面的流量计算为：
[0028]
[0029] 其中，Ava为0号子三棱柱单元与3号子三棱柱单元间相邻面的流量值，以3号子 三棱柱单元流向〇号子三棱柱单元为正。v丨_./+.，J+为原始三棱柱单元上表面的水 力坡度，\ -尤、.,._?/ -，J为原始三棱柱单元下表面的水力坡度，P表示边cb的向量；
[0030] 1号子三棱柱单元与3号子三棱柱单元间相邻面的流量计算为：
[0031]
[0032] 2号子三棱柱单元与3号子三棱柱单元间相邻面的流量计算为：
[0033]
[0034] 4号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0035]
[0036] 5号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0037]
[0038] 6号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0039]
[0040] 在本发明的一些实施例中，所述子三棱柱单元的边界流量值的计算方法包括：
[0041] 根据任意子三棱柱单元都均衡的原则，计算位于不同原始三棱柱单元中的一类均 衡域的子三棱柱单元的水平相邻面的流量值；
[0042] 根据任意子三棱柱单元都均衡的原则，计算位于不同原始三棱柱单元中的二类均 衡域的所有子三棱柱单元的水平相邻面的流量值，再根据所述相邻子三棱柱单元的水平相 邻面面积和等效渗透系数对所述二类均衡域水平相邻面的流量值进行摊派，计算得到所述 水平相邻面上各个子三棱柱的流量值；
[0043] 计算相邻的二类均衡域的子三棱柱单元的垂向面的流量值。
[0044] 在本发明的一些实施例中，所述相邻的二类均衡域的子三棱柱单元之间的垂向面 的流量值的计算方法包括：
[0045] 将由单层子三棱柱单元围成的三维区域对应为由三角形围成的二维区域，通过子 三棱柱围成三维区域的三角形顶面流量值qvtl和三角形底面流量值qvbl计算由三角形围成 二维区域的三角形的垂向流量值ε;，其中εqvti-qvbi,i= 1，2, 3, 4,5···;
[0046] 令三角形围成的二维区域内某一三角形边界的流量值为Q，以此参数化所 述二维区域内所有边界的流量值，利用参数化的二维区域内部边界流量值，并根据 Raviart-Thomas空间插值得到参数化的子区内部流速场，由此得到参数化的所述二维区域 内的所有边界的切向流速差；
[0047] 最