技术领域
[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为原始三棱柱单元的高度，Kzz为原始三棱柱单元的垂向渗透系数，ha、hb、hc、ha’、hb’和hc’分别为原始三棱柱六个顶点的水头值；
[0021] 1号子三棱柱单元与5号子三棱柱单元间相邻面的流量计算为：
[0022]
[0023] 2号子三棱柱单元与6号子三棱柱单元间相邻面的流量计算为：
[0024]
[0025] 3号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0026]
[0027] 0号子三棱柱单元与3号子三棱柱单元间相邻面的流量计算为：
[0028]
[0029] 其中，Aha为0号子三棱柱单元与3号子三棱柱单元间相邻面的流量值，以3号子三棱+柱单元流向0号子三棱柱单元为正。 J 为原始三棱柱单元上表面的水
力坡度，
J-为原始三棱柱单元下表面的水力坡度， 表示边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] 将由
单层子三棱柱单元围成的三维区域对应为由三角形围成的二维区域，通过子三棱柱围成三维区域的三角形顶面流量值qvti和三角形底面流量值qvbi计算由三角形围成二维区域的三角形的垂向流量值εi，其中εi＝qvti-qvbi,i＝1,2,3,4,5…；
[0046] 令三角形围成的二维区域内某一三角形边界的流量值为Q，以此参数化所述二维区域内所有边界的流量值，利用参数化的二维区域内部边界流量值，并根据Raviart-Thomas空间插值得到参数化的子区内部流速场，由此得到参数化的所述二维区域内的所有边界的切向流速差；
[0047] 最优化所述二维区域内部边界的切向流速差，计算出最优解或最优区间；如果计算得到的是最优区间，则选取其中使得所述二维区域内部流速场旋度最小的值为最优解；将所述最优解代入所述二维区域内部边界流量的参数化公式，计算得到所述二维区域内的所有边界流量值，即为子三棱柱单元之间的垂向面的流量值。
[0048] 在本发明的一些实施例中，所述计算子三棱柱单元内部流速场的步骤包括：
[0049] 利用子三棱柱顶、底面的流量值计算得到两个流速值，再利用顶、底面的流速值，线性插值得到子三棱柱内任意一点的垂向流速分量；
[0050] 子三棱柱的横向截面为一三角形，由子三棱柱三个垂向面的流量值作为三角形三条边的流量值，利用Raviart-Thomas空间插值得到三角形内任意一点的流速值，即为子三棱柱内任意一点的横向流速分量。
[0051] 从上面所述可以看出，本发明提供的地下水流场的三维流线可视化方法是基于高精度的地下水流速场和质点追踪算法得到的，该方法对于地下水流速场的计算，不仅考虑到了单元边界垂向流速的连续性、局部均衡性，还考虑到了介质边界切向流速的变化规律，从而保障了三维流线和合理性与准确性，能够真实地刻画地下水的补、径、排等运动特征。
附图说明
[0052] 图1为本发明实施例的地下水流场的三维流线可视化方法的流程示意图；
[0053] 图2为本发明实施例的研究区域模型的横向截面示意图；
[0054] 图3为本发明实施例的研究区域模型的断面示意图；
[0055] 图4为原始三棱柱单元细化为子三棱柱单元的示意图；
[0056] 图5为本发明实施例的一类局部均衡域构造示意图；
[0057] 图6为本发明实施例的二类局部均衡域构造示意图；
[0058] 图7为本发明实施例的一类局部均衡域和二类局部均衡域的俯视图；
[0059] 图8为本发明实施例的一类局部均衡域及三角形层间面示意图；
[0060] 图9为本发明实施例的二类局部均衡域及多边形层间面示意图；
[0061] 图10为本发明实施例中由子三棱柱围成的三维区域和由三角形围成的二维区域相应边界流量的对应关系图；
[0062] 图11为本发明实施例子三棱柱单元的内部流速场构造示意图
[0063] 图12为本发明实施例的三维流线平面图；
[0064] 图13为本发明实施例的三维流线剖面图；
[0065] 图14为本发明实施例的三维流线3D视图。
具体实施方式
[0066] 为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。
[0067] 参照附图1，为本发明提供的地下水流场的三维流线可视化方法的流程示意图。作为本发明的一个实施例，所述地下水流场的三维流线可视化方法包括以下步骤：
[0068] 步骤101：利用三棱柱连续Galerkin有限单元法求解地下水流数值模型，形成原始三棱柱单元网格并获得地下水的水头场数据。
[0069] 地下水数值模型是利用
电子计算机技术，将地下水概念模型通过空间离散和时间离散而建立的近似数学模型。它是反映自然界客观过程(模拟、再现自然界发生、发展的一些水文地质、工程地质过程)和定量评价的手段。
[0070] 水头场是指研究区域内任意一点对应一个水头值。在本发明中是指任意三棱柱单元顶点对应一个水头值。研究区域内其他点可以通过插值得到，所以相当于知道了研究区域内任意点的水头值。
[0071] 研究区域的模型如图2和3所示，图2和图3分别为本发明实施例的研究区域模型的横向截面示意图和断面示意图，该研究区域的模型由两层含水层和一层弱透水层组成，上层含水层为自由水面的潜水含水层，中间为弱透水层，在该弱透水层中存在一个圆形
天窗区域，下层含水层为承压含水层。模型的南北边界为隔水边界，西边界为水头值为50米的定水头边界，东边界为水头值为46米的定水头边界。下层的承压含水层中存在一口完整井，抽水强度为2000m3/day，抽水井和天窗的
位置如图2和3所示。弱透水层的渗透系数为含水层的1％，各含水层的垂向渗透系数为水平向渗透系数的十分之一。研究区域的网格剖分为：在水平向上，研究区域剖分为3750个三角形单元；在垂向上，研究区域剖分为3层，即潜水含水层、弱透水层和承压含水层。
[0072] 在该步骤中，利用三棱柱连续Galerkin有限单元法求解该研究区域的地下水流数值模型，从而形成原始三棱柱单元网格并获得地下水的水头场数据。
[0073] 步骤102：将所述原始三棱柱单元网格进行细化，从而构造局部均衡域。
[0074] 具体地，所述将原始三棱柱单元网格进行细化的步骤包括：
[0075] 如图4所示，原始三棱柱单元的横向剖面为一个三角形(剖面三角形)，连接剖面三角形三条边的中点，从而将该剖面三角形分割为面积相等的四个子三角形；再将该剖面三角形三条边的中点的连线沿着垂向连续地无限延展，从而将该原始三棱柱单元分割为四个中间三棱柱单元。然后，连接原始三棱柱单元的三条垂向边的中点，使其构成一个三角形(中线三角形)，该中线三角形所在的平面将四个中间三棱柱单元分别分割为两个子三棱柱单元。即：所述细化步骤将原始三棱柱单元通过水平向和垂向的细化，被分割为八个子三棱柱单元。
[0076] 为便于分清这八个子三棱柱单元，本发明将其进行编号。具体地，如图4所示，按照先上后下、先外后内、上下层顺序相同的原则，对八个子棱柱单元进行编号。记a、b、c分别为原始三棱柱单元上层的三个顶点，a’、b’、c’分别为原始三棱柱单元下层的三个顶点，围绕顶点a的子三棱柱单元记为0号，围绕顶点b的子三棱柱单元记为1号，围绕顶点c的子三棱柱单元记为2号，围绕顶点a’的子三棱柱单元记为4号，围绕顶点b’的子三棱柱单元记为5号，围绕顶点c’的子三棱柱单元记为6号，不围绕原始三棱柱顶点的上层子三棱柱单元记为3号，不围绕原始三棱柱顶点的下层子三棱柱单元记为7号。
[0077] 如图5所示，图中示出了两个相邻的原始三棱柱单元，按照图4所示的方法，分别将其细化为八个子三棱柱单元，令3号子三棱柱单元与相邻的原始三棱柱单元的7号子三棱柱单元所组成的区域为一类局部均衡域。或者，也可以令7号子三棱柱单元与相邻的原始三棱柱单元的3号子三棱柱单元所组成的区域为一类局部均衡域。
[0078] 如图6所示，令围绕原始三棱柱单元顶点的子三棱柱单元所组成的区域为二类局部均衡域。图7为本发明实施例的一类局部均衡域和二类局部均衡域的俯视图，其中灰色区域为一类局部均衡域的俯视图效果，白色区域为二类局部均衡域的俯视图效果。
[0079] 步骤103：计算所述局部均衡域的边界流量值。
[0080] 具体地，所述局部均衡域的边界流量值的计算方法如下：
[0081] 1)0号子三棱柱单元与4号子三棱柱单元间相邻面的流量计算为：
[0082]
[0083] 其中，Ava为0号子三棱柱单元与4号子三棱柱单元间相邻面的流量值，以0号子三棱柱单元流向4号子三棱柱单元为正。Δe为原始三棱柱单元剖面三角形的面积，Δz为原始三棱柱单元的高度，Kzz为原始三棱柱单元的垂向渗透系数，ha、hb、hc、ha’、hb’和hc’分别为原始三棱柱六个顶点的水头值。
[0084] 2)1号子三棱柱单元与5号子三棱柱单元间相邻面的流量计算为：
[0085]
[0086] 3)2号子三棱柱单元与6号子三棱柱单元间相邻面的流量计算为：
[0087]
[0088] 4)3号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0089]
[0090] 5)0号子三棱柱单元与3号子三棱柱单元间相邻面的流量计算为：
[0091]
[0092] 其中，Aha为0号子三棱柱单元与3号子三棱柱单元间相邻面的流量值，以3号子三棱柱单元流向0号子三棱柱单元为正。 J+为原始三棱柱单元上表面的水力坡度，J-为原始三棱柱单元下表面的水力坡度， 表示边cb的向量。
[0093] 6)1号子三棱柱单元与3号子三棱柱单元间相邻面的流量计算为：
[0094]
[0095] 7)2号子三棱柱单元与3号子三棱柱单元间相邻面的流量计算为：
[0096]
[0097] 8)4号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0098]
[0099] 9)5号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0100]
[0101] 10)6号子三棱柱单元与7号子三棱柱单元间相邻面的流量计算为：
[0102]
[0103] 通过上述10个公式，可以计算得到一类、二类局部均衡域的边界流量值。
[0104] 步骤104：计算子三棱柱单元的边界流量值。
[0105] 具体地，所述计算子三棱柱单元的边界流量值的步骤包括：
[0106] 一)计算一类局部均衡域的子三棱柱单元的边界流量值：如图8所示，经过步骤103的计算，位于不同原始三棱柱单元中的一类均衡域的子三棱柱单元的水平相邻面(如图8中白色填充的三角形所示)的流量值未知。需要说明的是，图8中的7号子三棱柱单元、3号子三棱柱单元分别位于不同且相邻的原始三棱柱单元中。根据任意子三棱柱单元都均衡的原则，即流入、流出任意子三棱柱单元的流量之和为零，计算7号子三棱柱单元与3号子三棱柱单元相邻面的流量值。例如，根据公式4)、公式8)、公式9)以及公式10)即得到位于不同原始三棱柱单元中的一类均衡域的两个子三棱柱单元的相邻面的流量值。
[0107] 二)计算二类局部均衡域的子三棱柱单元之间水平面的边界流量值：如图9所示，经过步骤103的计算，位于不同原始三棱柱单元中的二类均衡域的所有子三棱柱单元的水平相邻面(如图9中白色填充的五边形所示)的流量值未知。需要说明的是，所述二类均衡域的水平相邻面不限于五边形，结合图7所示，可能是七边形、六边形、四边形等其他多边形。同样地，根据任意子三棱柱单元都均衡的原则，以及子三棱柱单元构成的二类均衡域的上层和下层区域也都是均衡的，同时结合步骤103中的公式，既可以计算得到所述位于不同原始三棱柱单元中的二类均衡域的子三棱柱单元的水平相邻面的流量值。然后，再根据所述相邻子三棱柱单元的相邻面面积和等效渗透系数对所述二类均衡域水平相邻面的流量值进行摊派，计算得到所述水平相邻面上各个子三棱柱的流量值。此时，二类局部均衡域的上层、下层分别相当于由三角形围绕某一顶点组成的一个子区(即图9中的五边形)。
[0108] 三)计算相邻的二类均衡域的子三棱柱单元之间的垂向面的流量值，具体步骤如下：
[0109] 首先，将由子三棱柱围成三维区域的各边界面流量值对应到二维情况，即由三角形围成二维区域各边界线的流量和三角形面上的垂向流量值，如图10所示：由单层子三棱柱围成三维区域的五个矩形侧面的流量值qh1…qh5对应为由三角形围成二维区域的五条边界线的流量值q1…q5；通过子三棱柱围成三维区域的五个三角形顶面流量值qvt1…qvt5和五个三角形底面流量值qvb1…qvb5计算由三角形围成二维区域的五个三角形的垂向流量值，即ε1…ε5，其中εi＝qvti-qvbi,i＝1,2,3,4,5。因此，计算子三棱柱单元之间的垂向面的流量值即计算三角形围成的二维区域内的边界流量值。需要说明的是，εi是二维情况下三角形的垂向流量值，可以理解为三维情况下子三棱柱垂向方面的净(入/出)流量值。
[0110] 其次，假设三角形围成的二维区域内某一三角形边界的流量值为Q，以此参数化所述二维区域内所有边界的流量值，利用参数化的二维区域内部边界流量值，并根据Raviart-Thomas空间插值得到参数化的子区内部流速场。由此可得，参数化的所述二维区域内的所有边界的切向流速差。
[0111] 然后，最优化所述二维区域内部边界的切向流速差，计算出最优解或最优区间。如果计算得到的是最优区间，则选取其中使得所述二维区域内部流速场旋度最小的值为最优解。将最优解代入所述二维区域内部边界流量的参数化公式，计算得到所述二维区域内的所有边界流量值，即子三棱柱单元之间的垂向面的流量值。
[0112] 步骤105：计算子三棱柱单元内部的流速场。
[0113] 由于步骤103和步骤104计算得到了子三棱柱所有面的流量值，在此基础上，计算子棱柱单元内部的流速场。
[0114] 具体地，所述计算子三棱柱单元内部流速场的步骤包括：
[0115] 计算子三棱柱单元内任意一点的垂向流速分量：如图11所示，利用子三棱柱顶、底面的流量值(qv1和qv2)可以计算得到两个流速值，再利用顶、底面的流速值，线性插值可以得到子三棱柱内任意一点的垂向流速分量。
[0116] 计算子三棱柱单元内任意一点的横向流速分量：如图11所示，子三棱柱的横向截面为一三角形，由子三棱柱三个垂向面的流量值作为三角形三条边的流量值，利用Raviart-Thomas空间插值可以得到三角形内任意一点的流速值，即子三棱柱内任意一点的横向流速分量。
[0117] 步骤106：利用质点追踪算法形成水流场的三维流线。
[0118] 质点追踪算法是根据人为给定的或按照某种原则自动布置的流线初始点的位置及该位置的流速适量的方向与大小，计算质点的运动规律，按照向前质点追踪至汇或按照向后质点追踪至源为止，所得到的质点的轨迹线即为地下水三维流线。
[0119] 在本实施例中，布置流线初始点位置上层潜水含水层的西部边界处，垂向上位于潜水含水层的中心位置，共计24个流线初始点。按照向前质点追踪至汇为止，生成24条流线，描述了地下水的潜水补给到排泄的地下水运动过程。生成的三维流线如图12-14所示，图12-14分别为地下水三维流线平面图、地下水三维流线剖面图和地下水三维流线3D视图。三个不同角度的视图模式清晰描绘了地下水的运动规律特征，即：潜水补给主要由天窗区域向下层承压含水层补给，也有一部分通过中间弱透水层向下层承压含水层补给，并最终由抽水井排泄。抽水井严重改变了倾斜地下含水层的运动规律，形成了扫把状的降落漏斗。
虽然，抽水井为下层承压水含水层完整井，但是由于天窗的导水作用，改变了原始的地下水以
层流为主的运动规律，形成了地下水的三维运动规律。
[0120] 由此可见，本发明提供的地下水流场的三维流线可视化方法是基于高精度的地下水流速场和质点追踪算法得到的，该方法对于地下水流速场的计算，不仅考虑到了单元边界垂向流速的连续性、局部均衡性，还考虑到了介质边界切向流速的变化规律，从而保障了三维流线和合理性与准确性，能够真实地刻画地下水的补、径、排等运动特征。
[0121] 所属领域的普通技术人员应当理解：以上任何实施例的讨论仅为示例性的，并非旨在暗示本公开的范围(包括
权利要求)被限于这些例子；在本发明的思路下，以上实施例或者不同实施例中的技术特征之间也可以进行组合，并存在如上所述的本发明的不同方面的许多其它变化，为了简明它们没有在细节中提供。因此，凡在本发明的精神和原则之内，所做的任何省略、
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
