CN109918748B - 基于随机游走的非均质含水层水流问题评估方法 - Google Patents
基于随机游走的非均质含水层水流问题评估方法 Download PDFInfo
- Publication number
- CN109918748B CN109918748B CN201910136315.1A CN201910136315A CN109918748B CN 109918748 B CN109918748 B CN 109918748B CN 201910136315 A CN201910136315 A CN 201910136315A CN 109918748 B CN109918748 B CN 109918748B
- Authority
- CN
- China
- Prior art keywords
- aquifer
- boundary
- random walk
- water flow
- grid
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于随机游走的非均质含水层水流问题评估方法，该方法首先建立随机游走密度与格林函数之间的定量关系；将高度非均质含水层离散化；利用网格随机游走方法在非均质介质中模拟大量的随机游走实现；统计随机游走在含水层中的密度；将随机游走密度转化为所求的格林函数。与传统水流问题格林函数技术相比，本发明的方法可以处理一般的、实际的含水层非均质问题，应用范围更广，适用性更强。
Description
技术领域
本发明属于水力学技术领域，具体指代一种基于随机游走的高度非均质含水层中水流问题评估方法。
背景技术
格林函数应用于地下水问题求解已超过四十年。格林函数法的特点在于，所求问题对应的格林函数已知的情况下，原问题的定解条件和源汇项组合改变时，无需再次求解，而可以直接把定解条件、源汇项代入，即可求得新解，相比于其他方法具有显著的优势。已有一些文献总结整理了目前用格林函数法可解的偏微分方程问题类型库，这些基本都是考虑介质均质或分块的情况。均质的情况比较基础和简单；对于分块的情况，需要借助伽辽金法，对每个不同参数分区引入一个基函数。而真实的含水层常常具有高度非均质性(性质处处变化)和不规则几何结构。若使用伽辽金法，则基函数过多，而且分区界面也过于复杂，无法实际操作。因此在非均质含水层中应用格林函数法非常困难。
文献表明，在均质情况下，随机游走密度与格林函数之间存在简单的正比关系(参见Deaconu和Lejay于2006年发表的“A Random Walk on Rectangles Algorithm”等文献)，由此在均质条件下，可以通过随机游走评估格林函数。然而，在非均质条件下，“是否存在类似的关系”、“如果存在具体表达式如何”这两个问题还有待回答。
发明内容
针对于上述现有技术的不足，本发明的目的在于提供一种基于随机游走的非均质含水层水流问题评估方法，以给出随机游走密度与格林函数的定量关系；本发明回避直接求解格林函数的困难，采用了随机游走模拟求得游走密度再间接计算格林函数。
为达到上述目的，本发明采用的技术方案如下：
本发明的一种基于随机游走的非均质含水层水流问题评估方法，包括步骤如下：
(1)根据所要研究的非均质含水层形状，将含水层空间用网格进行剖分离散；
(2)将含水层中第一类边界作为吸收边界，第二类边界作为反弹边界，第三类边界作为半吸收-半反弹边界；
(3)根据网格剖分，计算每个网格节点到相邻网格节点的转移概率；
(4)确定需要计算格林函数的位置，作为后续随机游走出发点；
(5)以出发点为起点，用网格随机游走法重复模拟N个随机游走，并记录每个网格节点被所有随机游走路过的次数Ni，i为网格节点编号；
(6)计算每个网格节点被随机游走路过的平均密度qi＝Ni/N；
(7)根据平均密度qi和网格剖分参数，计算出发点对应的格林函数；
(8)代入格林函数形式解中得水流方程的解。
进一步地，所述步骤(1)中非均质含水层的水流问题用方程描述如下：
其中，x为空间向量，K(x)为随位置变化的渗透系数，w(x,t)为源汇项，Ss为贮水系数，h(x,t)为待求水头，D为含水层模型空间范围，HD(x,t)、f(x,t)、b(x)、r(x,t)为已知函数，ΓD、ΓN、ΓR分别为第一、二、三类边界，n为边界外向单位法向长度。
该方程的解可以表示为h(x,t)＝HS+HR+HD+HI，其中
HI＝∫ΩG(x,t|s,0)H0(s)ds(详见文献Cole等“Heat Conduction Using Green’s Functions”一书)；其中，G(x,t|s,τ)为待求的格林函数，s、τ为积分参数。
进一步地，所述步骤(1)中的将含水层空间用网格进行剖分离散，具体包括：采用有限元法或有限差分法的离散方式，基于均匀网格或非均匀网格对含水层进行剖分离散。
进一步地，所述步骤(2)中的吸收边界会将到达该处的行走器吸收，使该次随机游走终止；反弹边界将到达该处的行走器反弹回孔隙空间；半吸收-半反弹边界是前两类边界的混合，使随机游走终止，或使其反弹。
进一步地，所述步骤(3)具体包括：当前网格节点到相邻网格节点的转移概率通过相关点之间的空间及时间距离和含水层介质性质计算得到。
具体来说，以三维问题有限差分法为例，空间网格上一个节点最多有6个相邻网格节点，加上时间上与该节点上一时刻的相邻关系，共有7个相邻网格节点，因此共需考虑向这些相邻网格节点的7个转移能力(C)和转移概率(p)；向前一时刻转移能力定义为：
其中Δt为剖分时间步长，Δx、Δy、Δz分别为当前网格节点三个方向上的剖分大小，SS为该节点的贮水系数；向空间相邻网格节点的转移能力定义为与这些点间的水力传导度(hydraulic conductance，详见MODFLOW手册)，并分别记为Cx+、Cx-、Cy+、Cy-、Cz+、Cz-；该网格节点总转移能力为Cs＝Ct+Cx++Cx-+Cy++Cy-+Cz++Cz-，而向7个相邻网格节点的转移概率分别为对应的转移能力与Cs之比；例如，x+方向上转移概率px+＝Cx+/Cs；据此计算所有网格节点到相邻网格节点的转移概率，对于稳态问题，Ct≡0。
进一步地，所述步骤(4)中的出发点即为格林函数G(x0,t0|s,τ)中的变量x0和t0。
进一步地，所述步骤(7)中的平均密度qi与格林函数G(x0,t0|si,τi)的转换关系为
(动态水流问题)及
(稳定水流问题)，si和τi为编号i的网格节点对应的空间位置和时刻。
进一步地，所述步骤(8)中，将步骤(7)中求得的数值格林函数，代入格林函数形式解的表达式即得所求水流问题的解。
本发明的有益效果：
1、本发明可以将格林函数法扩展到一般非均质含水层水流问题，使格林函数法适用于实际含水层，适用范围大大增加。
2、与现有含水层模拟工具相比，本发明的格林函数解可以快速计算实际含水层在变化环境下的响应，对真实含水层快速预警等应用具有重要意义。
3、本发明对真实的复杂含水层也适用，在求得格林函数的基础上计算效率和精度都很高，可以作为场景模拟计算中的替代模型以节省大量计算时间。
附图说明
图1为本发明执行中对二维空间和一维时间进行离散所得的三维时空网格示意图。
图2为实施例1中傍河均质潜水含水层示意图。
图3a为实施例1中t0＝10d时x0＝0.25L位置上计算的格林函数与对应解析解结果对比示意图。
图3b为实施例1中t0＝10d时x0＝0.50L位置上计算的格林函数与对应解析解结果对比示意图。
图3c为实施例1中t0＝10d时x0＝0.75L位置上计算的格林函数与对应解析解结果对比示意图。
图3d为实施例1中t0＝10d时x0＝L位置上计算的格林函数与对应解析解结果对比示意图。
图4a为实施例2中渗透系数有1～3个分区情况下x0＝0.15L位置上计算的格林函数与解析解结果对比示意图。
图4b为实施例2中渗透系数有1～3个分区情况下x0＝0.25L位置上计算的格林函数与解析解结果对比示意图。
图4c为实施例2中渗透系数有1～3个分区情况下x0＝0.35L位置上计算的格林函数与解析解结果对比示意图。
图5a为实施例3中在中心位置的格林函数解析解示意图。
图5b为实施例3中在中心位置的格林函数数值解示意图。
图6a为实施例4中的二维非均质潜水含水层示意图。
图6b为实施例4中渗透系数示意图。
图6c为实施例4中给水度空间分布示意图。
图7a为实施例4中抽水井位置上0～20d内的格林函数示意图。。
图7b为实施例4中τ＝20d时刻的格林函数空间分布示意图。
图7c为实施例4中τ＝10d时刻的格林函数空间分布示意图。
图7d为实施例4中τ＝0d时刻的格林函数空间分布示意图。
图8a为实施例4中随时间变化的抽水量示意图。
图8b为实施例4中随时间变化的降水量示意图。
图8c为实施例4中观测井#1处的水头变化示意图。
图8d为实施例4中观测井#2处的水头变化示意图。
图8e为实施例4中抽水井处的水头变化示意图。
具体实施方式
为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。
图1表示了一个三维时空(空间二维+时间一维)问题的网格离散，每一层网格对应一个离散时刻；对于稳态问题，则只有一层。x0和t0分别表示所要求的格林函数的空间位置和时刻，粗线箭头示意该三维网格上的一个随机游走实现。根据含水层介质性质(K及SS)以及网格剖分(Δx、Δy、Δt)计算所有网格节点的同层(空间)转移概率和向下层(时间)的转移概率。稳定水流问题中网格只有一层，不存在时间上多层，不需考虑时间。从(x0,t0)出发根据转移概率进行随机游走直到被边界吸收或到达最底层(t＝0)，记录每个网格节点被路过的次数；共重复随机游走N次，对每个网格节点汇总其在N次游走中被路过的总次数(例如，网格节点i被路过的总次数为Ni)。然后计算每个网格节点被路过的平均次数(平均密度)qi＝Ni/N。之后根据关系式
(动态水流问题)或
(稳定水流问题)，计算每个网格节点上的格林函数值。最后，将定解条件、源汇项代入，利用求得的格林函数计算水流问题的解。
实施例1：傍河均质潜水含水层一维非稳定流
研究对象为地下分水岭与河流之间的均质潜水含水层(左右两侧分别对应第二、三类边界)。含水层下方为水平隔水底板，上方开放有入渗补给，如图2所示。L＝100m，渗透系数为K＝1m/d，给水度μ＝0.3，含水层平均厚度
河床渗透系数Kb＝0.1m/d，河床厚度Lb＝0.1m。左侧边界为零通量边界，表示为
右侧边界为第三类边界，表示为
其中hb(t)为河水位。模拟期为(0,10d)。可以令h2＝u，对原潜水方程进行线性化后方便求解。空间离散步长Δx＝0.5m，时间离散步长Δt＝0.05d。本实施例存在传统格林函数的解析解，可以对WOG计算的格林函数数值解进行验证，这也是本例中考虑均质含水层的原因。我们对t0＝10d时x0＝0.25L、0.50L、0.75L、L四个位置的格林函数进行求解。图3a-图3d中比较了解析解(实线)与本发明所得的数值解(虚线)结果，表明本发明的计算结果与解析解相比除了极小的数值误差以外几乎完全一致，证明本发明在均质含水层中的有效性。
实施例2：具有一/二/三个K分区的承压含水层一维稳定流；
考虑一维承压含水层，两侧边界均为第一类边界。含水层渗透系数K为均质(一个分区)或非均质(二或三个分区)，三种情况如图4a-图4c中顶部子图所示。对于K的三种空间分布情况，我们计算了x0＝0.15L、x0＝0.25L、x0＝0.35L三个位置的格林函数。图4a-图4c中对比了解析解(细实线)与本发明的数值解(粗虚线)的结果，表明本发明的计算结果与解析解高度一致，证明本发明在参数分区含水层稳定流问题中的有效性。
实施例3：均质承压含水层二维非稳定流
考虑一平面二维矩形均质承压含水层中的非稳定流问题。含水层长宽分别为1440m和800m。含水层东西边界为第一类边界，南北边界为第二类边界。为了验证本发明对该问题的有效性，分别用解析解和本发明计算了t0＝10d中心位置的格林函数。图5a和图5b分别表示了这两种方法所得的格林函数结果。图中表明两种方法所得结果高度一致(注意本发明的数值格林函数有一定的数值计算误差)，证明了本发明对该问题的有效性。
为了能有精确解析解进行对比，以上三个实施例考虑的是均质或简单分区含水层问题。为了进一步说明本发明对真实含水层水流问题的普适性，我们考虑一个一般性的非均质含水层问题。
实施例4：高度非均质潜水含水层二维非稳定流
如图6a-图6c所示，一长宽分别为1000m、800m的矩形非均质含水层，东边界为已知水头边界(水头为已知函数
)，南、西、北侧为隔水边界。底部为水平隔水底板，标高为0。含水层平均厚度50m。渗透系数K的对数和给水度μ的空间分布如图6b、图6c所示。一抽水井位于(xw,yw)＝(500m,400m)，两个观测井#1和#2分别位于(x1,y1)＝(450m,400m)和(x2,y2)＝(550m,400m)。同样令h2＝u对潜水方程线性化并用格林函数法求解。图7a展示了抽水井位置上0～20d内的格林函数。图7b、图7c、图7d则分别显示了t＝20d、10d、0d三个时刻的格林函数空间分布，清楚地显示出格林函数随时间的变化以及介质非均质性对格林函数的影响。由于本例问题没有格林函数解析解可以对比，我们用算得的格林函数计算抽水井和观测井的水头，将水头与地下水领域广泛接受的MODFLOW软件算得的水头进行对比，间接验证格林函数的正确性。为了保证检验的可靠性，以下所用的源汇项具有很强的时空变化性。假设东部边界水头为50m；含水层内部初始水头也均为50m；抽水井的抽水量随时间变化如图8a所示；降雨量随时间变化如图8b所示。由此可以计算两观测井及抽水井位置上的水位随时间的变化，如图8c、图8d、图8e所示，其中格林函数法计算的水位变化用粗虚线表示，MODFLOW的计算结果用细实线表示。从图中可以看出，两种方法的结果高度一致，表明本发明计算出的格林函数是非常准确的。
本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。
Claims (6)
1.一种基于随机游走的非均质含水层水流问题评估方法，其特征在于，包括步骤如下：
(1)根据所要研究的非均质含水层形状，将含水层空间用网格进行剖分离散；
(2)将含水层中第一类边界作为吸收边界，第二类边界作为反弹边界，第三类边界作为半吸收-半反弹边界；
(3)根据网格剖分，计算每个网格节点到相邻网格节点的转移概率；
(4)确定需要计算格林函数的位置，作为后续随机游走出发点；
(5)以出发点为起点，用网格随机游走法重复模拟N个随机游走，并记录每个网格节点被所有随机游走路过的次数Ni，i为网格节点编号；
(6)计算每个网格节点被随机游走路过的平均密度qi＝Ni/N；
(7)根据平均密度qi和网格剖分参数，计算出发点对应的格林函数；
(8)代入格林函数形式解中得水流方程的解；
所述步骤(1)中非均质含水层的水流问题用方程描述如下：
其中，x为空间向量，K(x)为随位置变化的渗透系数，w(x,t)为源汇项，Ss为贮水系数，h(x,t)为待求水头，D为含水层模型空间范围，HD(x,t)、f(x,t)、b(x)、r(x,t)为已知函数，ΓD、ΓN、ΓR分别为第一、二、三类边界，n为边界外向单位法向长度。
2.根据权利要求1所述的基于随机游走的非均质含水层水流问题评估方法，其特征在于，所述步骤(1)中的将含水层空间用网格进行剖分离散，具体包括：采用有限元法或有限差分法的离散方式，基于均匀网格或非均匀网格对含水层进行剖分离散。
3.根据权利要求1所述的基于随机游走的非均质含水层水流问题评估方法，其特征在于，所述步骤(2)中的吸收边界会将到达该处的行走器吸收，使该次随机游走终止；反弹边界将到达该处的行走器反弹回孔隙空间；半吸收-半反弹边界是前两类边界的混合，使随机游走终止，或使其反弹。
4.根据权利要求1所述的基于随机游走的非均质含水层水流问题评估方法，其特征在于，所述步骤(3)具体包括：当前网格节点到相邻网格节点的转移概率通过相关点之间的空间及时间距离和含水层介质性质计算得到。
5.根据权利要求1所述的基于随机游走的非均质含水层水流问题评估方法，其特征在于，所述步骤(4)中的出发点即为格林函数G(x0,t0|s,τ)中的变量x0和t0。
6.根据权利要求1所述的基于随机游走的非均质含水层水流问题评估方法，其特征在于，所述步骤(8)中，将步骤(7)中求得的数值格林函数，代入格林函数形式解的表达式即得所求水流问题的解。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910136315.1A CN109918748B (zh) | 2019-02-25 | 2019-02-25 | 基于随机游走的非均质含水层水流问题评估方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910136315.1A CN109918748B (zh) | 2019-02-25 | 2019-02-25 | 基于随机游走的非均质含水层水流问题评估方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109918748A CN109918748A (zh) | 2019-06-21 |
| CN109918748B true CN109918748B (zh) | 2023-04-07 |
Family
ID=66962200
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910136315.1A Active CN109918748B (zh) | 2019-02-25 | 2019-02-25 | 基于随机游走的非均质含水层水流问题评估方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109918748B (zh) |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106886682B (zh) * | 2017-01-04 | 2020-01-07 | 中国环境科学研究院 | 用于单裂隙中溶质运移数值模拟的随机行走粒子追踪方法 |
| CN109284555B (zh) * | 2018-09-28 | 2023-02-14 | 南京大学 | 基于多孔介质几何形状估计渗透率的网格随机游走方法 |
-
2019
- 2019-02-25 CN CN201910136315.1A patent/CN109918748B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109918748A (zh) | 2019-06-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hunt et al. | Flow, transport, and reaction in porous media: Percolation scaling, critical‐path analysis, and effective medium approximation | |
| Hyman et al. | Predictions of first passage times in sparse discrete fracture networks using graph-based reductions | |
| Liu et al. | A fractal model for characterizing fluid flow in fractured rock masses based on randomly distributed rock fracture networks | |
| Rees et al. | A three-dimensional numerical model of borehole heat exchanger heat transfer and fluid flow | |
| Di Stefano et al. | Flow resistance equation for rills | |
| Begnudelli et al. | Conservative wetting and drying methodology for quadrilateral grid finite-volume models | |
| Hunt et al. | Unsaturated hydraulic conductivity modeling for porous media with two fractal regimes | |
| Zhu et al. | A fully coupled numerical modeling for regional unsaturated–saturated water flow | |
| Liggett et al. | Influence of the first-order exchange coefficient on simulation of coupled surface–subsurface flow | |
| Brunetti et al. | A hybrid finite volume-finite element model for the numerical analysis of furrow irrigation and fertigation | |
| Mousavi Nezhad et al. | Stochastic finite element modelling of water flow in variably saturated heterogeneous soils | |
| Hyman et al. | Characterizing the impact of fractured caprock heterogeneity on supercritical CO 2 injection | |
| WO2012021292A1 (en) | Reservoir upscaling method with preserved transmissibility | |
| Xue et al. | A fast numerical method and optimization of 3D discrete fracture network considering fracture aperture heterogeneity | |
| Bautista et al. | New results for an approximate method for calculating two-dimensional furrow infiltration | |
| CN104091065A (zh) | 一种求解浅水问题模拟间断水流数值的方法 | |
| Wang et al. | The importance of capturing topographic features for modeling groundwater flow and transport in mountainous watersheds | |
| Esmaeilpour et al. | Scale-dependent permeability and formation factor in porous media: Applications of percolation theory | |
| Fan et al. | Establishment and verification of the prediction model of soil wetting pattern size in vertical moistube irrigation | |
| Naghedifar et al. | Optimization of quadrilateral infiltration trench using numerical modeling and Taguchi approach | |
| An et al. | High-order averaging method of hydraulic conductivity for accurate soil moisture modeling | |
| CN109918748B (zh) | 基于随机游走的非均质含水层水流问题评估方法 | |
| Zhang et al. | Modeling free-surface seepage flow in complicated fractured rock mass using a coupled RPIM-FEM method | |
| He et al. | A new semi-analytical method for estimation of anisotropic hydraulic conductivity of three-dimensional fractured rock masses | |
| Xu et al. | A multipoint flux mixed finite element method for Darcy–Forchheimer incompressible miscible displacement problem |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |