FEFLOW插件PINK用于地下水开采引起的地面沉降的完全水力-力学耦合分析

SoftwareX 9（2019）15原始软件出版物FEFLOW中的完全水力-力学耦合插件（FLOW+）用于地下水开采引起的地面沉降Hung Tien PhamHung，Wolfram Rühaak，Valerian Schuster，Ingo Sass德国达姆施塔特工业大学应用地球科学研究所地球科学与ar t i cl e i nf o文章历史记录：2018年7月9日收到收到修订版2018年11月7日接受2018年11月8日保留字：地下水模拟地面沉降FEFLOW Biot全耦合流体力学a b st ra ct大量开采地下水导致地下水位下降，这会增加土壤的应力并导致地面沉降。地面沉降分析一般需要基于孔隙弹性理论的数值模型FEFLOW是一个著名的商业有限元软件，是最广泛用于地下水，传热和运输建模的软件之一然而，它不能模拟地面沉降。在这里，我们介绍了一个插件，名为PINK+和预处理软件PINK，可用于完全耦合的地面沉降分析，使用孔隙弹性理论。插件和软件是使用C++和Qt-Framework开发的。Eigen和英特尔MKL库用于矩阵运算和稀疏矩阵求解器。虽然FEFLOW是封闭源代码的，但PIFLOW+和PIFLOW是开源的，可以在GitHub上使用。版权所有©2018作者.由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本V1此代码版本所用代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2018_104法律代码许可证MIT，LGPL使用的代码版本控制系统使用C++、OpenGL的软件代码语言、工具和服务编译要求、操作环境依赖性Qt-5静态，特征库，英特尔MKL如果可用，链接到开发人员文档/手册https://github.com/pham-hung/Subsidence_FEFLOW/blob/master/Documents/User_manual.pdf问题支持电子邮件pham@geo.tu-darmstadt.de软件元数据当前软件版本V1此版本可执行文件的永久链接https://github.com/pham-hung/Subsidence_FEFLOW法律软件许可证MIT，LGPL计算平台/操作系统Windows x64安装要求依赖项如果可用，请链接到用户手册-如果正式出版，请在参考列表https://github.com/pham-hung/Subsidence_FEFLOW/blob/master/Documents/User_manual.pdf问题支持电子邮件pham@geo.tu-darmstadt.de1. 动机和意义地下水开采引起的地面沉降问题在世界许多地区都是一个重要问题，*通讯作者。电子邮件地址： pham@geo.tu-darmstadt.deH.T. Pham）。https://doi.org/10.1016/j.softx.2018.11.004例如，中国上海[1]、印度尼西亚雅加达[2]、泰国曼谷[3]、墨西哥城[4]、越南湄公河三角洲[5]、美国加利福尼亚州[6]和德国东北部萨克森州[7]。通常，地面沉降模型是由水文地质学家在不存在水平变形的假设下开发的[6]。因此，地面沉降问题可以用两个非耦合模型或两步模型来解决。首先，地下水模型用于2352-7110/©2018作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx16H.T. Pham，W.吕哈克河谷Schuster等人/SoftwareX 9（2019）15==++=阿阿阿克阿克埃埃拉拉材料. 1xx3阿克3埃3拉YY3埃3阿克3拉σxz =σzx= −GBZZ+BZXYZzy埃拉获得地下水位的变化，导致有效应力的差异。然后，使用解析方法或数值模型计算沉降量。这种方法 通常用于水文地质学领域，因为易于实施，特别是使用有限差分法 -例如，具有层间存储包的MODFLOW [8]。然而，地面沉降是一个三维（3D）问题。忽略水平变形可能导致高估沉降的幅度，特别是在抽水井附近[6]。在对地下水开采引起的区域地面沉降的回顾中，Galloway和Burbey [6]指出，基于孔隙弹性理论的模型对于分析含水层设置的实际行为是必要的。孔隙弹性理论最初是由Biot [9]基于太沙基的一维固结理论发展起来的。后来，Biot和Verruijt进一步扩展了它[10，11]。FEFLOW是一个常用的地下水、热传导和输运模拟的有限元软件，然而，FEFLOW不能模拟地面沉降。 FEFLOW提供了一个API（应用程序编程接口）来开发插件。根据孔隙弹性理论，研制了一种封堵器-在，命名为“+”，不仅可以用于地面沉降，Fig. 1. 笛卡尔坐标系中具有总应力分量的代表性单元。孔隙弹性理论的系统方程包括储水方程（或地下水方程）和应力平衡方程。存储方程定义为：问题，但也一般巩固分析。而且有开发了一个名为PADER的软件包，<$ε<$H<$α +S−（kxH）−（kyH）−（kzH）=0（2）从+使用FLOW+，FEFLOW充当预处理工具。的几何、网格、输入材料参数、时间步长设置和边界条件在FEFLOW中处理。随后，XNUMX+使用这些数据来组装全局矩阵并求解方程组。最后，将来自P2P+的结果保存为ASCII文件，并使用P2P进行可视化。Python+和Python是使用C++和跨平台框架Qt开发的。源代码可以在Windows和Linux上编译。矩阵运算与稀疏矩阵其中，ε εxxεyyεzz是总体积应变;α是比奥存储率定义为：S=nCf+（α−n）Cs（3）式中n为孔隙度;Cf为水的压缩性;Cs为固体颗粒（或土壤骨架）的压缩性应力平衡方程为：求解器基于Eigen和英特尔MKL库[12，13]，保证了XNUMX+的快速和可靠。使用分析物验证了E10 +<$σx′x+<$σyx+<$σxz+布吕普α−fx=0太沙基和德勒夫问题的物理解。对于Windows最终用户，Windows+和PXP作为单个.dll和.exe提供xσxy +α−fy=0（四）文件，分别。所有第三方库都是静态链接的因此，Windows用户可以立即使用Windows+和PSDK，而无需任何配置。本文概述了孔隙弹性理论和P2O3+和P2O3+的主要有关安装、理论、验证和用户手册的详细信息将单独提供，心理材料。12. 理论我们介绍了孔隙弹性理论（或Biot详细的信息可以在Biot [9，10，14，15]，Verruijt [11]的作品中找到，并在补充说明中阿姆斯特鲁吉σxz<$σyz<$σz′z阿克斯 公司简介 公司简介+αz−fz=0其中σij和σ′ij分别为总应力和有效应力分量，fi为体积力。总应力是有效应力和孔隙压力之和σij=σi′j+αδp=σi′j+αδ（H-Ele）γw（5）其中δ是克罗内克δ。有效应力和应变之间的关系遵循胡克σ′= −（K+4G）<$u −（K−2G）<$v −（K−2G）<$wσ′= −（K+4G）<$v −（K−2G）<$u −（K−2G）<$w我们考虑一个三维的代表性元素σ′= −（K+4G）<$w −（K−2G）<$u −（K−2G）<$v（图1），其中x、 y、 z方向上的位移表示为zz3个月300万3埃什基（六）分别表示为u、v和w孔隙水压力表示为σ′=σ′= −G（n+n）;p，总压头为H。总水头与孔隙压力的关系为：xy y x′ ′埃什基（u阿克斯（见附件）其中γw是水的单位重量，Ele是高程。σ′=σ′= −G（w+u）1补充材料可在github链接www.example.com pham-hung/Subsidence_FEFLOW/blob/master/Documents/User_manual.pdf中找到https://github.com/。其中K1/Cm是体积模量;Cm是多孔介质的压缩性;G是剪切模量。2.1. 完全耦合孔隙弹性理论H= pγw+ Ele（1）;H.T. Pham，W.吕哈克河谷Schuster等人/SoftwareX 9（2019）1517∇ ∇--H阿ZZvZZvZZ−LM +PZH−tP×H0+tQ阿场;L是耦合矩阵;拉乌波夫·T 是阿2图二、 FEFLOW程序流程及VBA + Plug-In的主要功能。等式（2）和（4）是包含四个变量的全耦合方程：x、y、z方向上的三个位移和孔隙压力p（或总水头H）。当量（2）意味着总体积的变化，它有两个分量α（ε ε/ε t）和S（εH/ε t），等于离开元素的水量（k H）。前一个分量是由有效应力的变化引起的，后一个分量是由孔隙压力2.2. 非耦合方法--孔隙弹性理论的一种简单形式在传统的地下水理论中，假定没有水平变形。如果重力方向是负z方向，我们有：εxx=εyy=0（7）如果多孔介质表现为弹性材料，则：ε=ε= −mσ′= −m（σ −αp）当量（11）是FEFLOW中新的参数Ss包括水、固体颗粒和多孔介质的压缩性当量（11）只以总水头为变量，或者是Biot理论的耦合方程组的解耦形式现在解决地面沉降问题分两步走：– 步骤1：求解方程（11）获得每个时间步的总头场的结果。– 第二步：根据水头的变化计算每个时间步的竖向变形。当不存在水平变形和总垂直应力恒定时，非耦合模型是适用的。在某些情况下，当水平移动较大时（例如，靠近抽水井），这种方法会高估一维沉降[6]。2.3. 有限元逼近应用Galerkin近似方法，方程的ward积分方案（2）和（4），耦合方程= −mv[σzz−αγw（H−Ele）]（8）式中mv=1/（K+4G/ 3）为多孔介质的有限压缩系数以矩阵形式表示：[KTL]{d}={F（十二）∂εα= −αmv联系我们阿ZH+αmv（九）式中，Δt为计算时间步长;K为位移场矩阵，M和P为瞬态矩阵替换Eq。（9）在Eq.（2）给出：增量位移矢量;ΔH是增量水（α2m+S）=（kH）+m联系我们阿（十）头向量为上一计算步骤的水头如果假定总应力随时间恒定，则右侧的第二项为零，并且等式（10）变成：SsH=（kH）（11）3. +插件和P图2显示了FEFLOW中有限元模型的程序流程。FEFLOW有三个主要模块：媒体从等式（8）我们有：阿vv是增量负载向量;Q是提取速率向量。18H.T. Pham，W.吕哈克河谷Schuster等人/SoftwareX 9（2019）15= ×−=2个=– 问题编辑器：编辑模型的几何、输入参数和边界条件。– 模拟器：求解方程组– 后处理：可视化结果。在问题编辑器阶段和模拟器阶段调用了SQL+。而PSTOM是后处理阶段的替代方案。+的主要功能是：– 在问题编辑器阶段：+创建一个保存的文件夹，一个设置文件夹，用户节点和元素数据，以分配边界条件和输入参数的机械领域。在FEFLOW中以常规方式处理其他输入数据– 在模拟阶段：+收集与节点、元素、材料属性、计算时间和边界条件相关的所有信息。然后，+使用这些数据来组装全局耦合矩阵并求解方程组。最后，计算单元应力并将结果保存为ASCII文件到保存的文件夹中。Pinteger导入由Pinteger+生成的输出数据并将其可视化。Pendix可以缩放坐标，缩放变形，显示特定计算步骤的结果或制作动画。User_manual. pdf文件中提供了有关使用Windows+和PXP的分步说明3.1. 使用多线程组装全局矩阵有两种特定的技术用于提高Linux+的速度第一个是组装方程中的全局矩阵（12）使用多线程。第二种方法是在单元刚度矩阵组装成整体矩阵之前直接将Dirichlet边界图图3显示了如何并行汇编全局矩阵。可用线程数为m，元素数为noe。每个线程将在n个元素上循环以计算n个元素刚度矩阵。单元刚度矩阵的每个系数被存储为向量的元素。使用m个可用线程，创建m个向量。向量的元素是元组类型。向量元素包含行索引i、列索引j和来自元素矩阵的当所有线程完成时，所有向量被连接为单个向量，并使用此向量生成全局矩阵3.2. 对每个单元刚度矩阵应用Dirichlet边界条件全局矩阵是稀疏的。只有非零值与行索引和列索引一起存储。如果Dirichlet边界条件图三. 多线程汇编的全局矩阵算法.在组装全局矩阵之后应用，非零系数的数量这是耗时的，而不是内存效率。为了克服这个困难，Dirichlet边界条件被施加到每个单元刚度矩阵之前，它被组装成整体矩阵。由于元素矩阵是小的稠密矩阵，所以这个运算是平凡的。4. 例如在本节中，我们考虑一个有14层的含水层系统（类似于越南湄公河三角洲的情况[5]）。模型数据在表1中给出，其中kxx、kyy、kzz是水力传导率（kxx 3kzz），并且K是体积模量。所有层被认为具有相同的孔隙度n 0。6和固体颗粒的压缩性Cs1 10−10 kN/ m。系统的初始总水头为0.0 m。模拟时间为10年，时间步数为365。边界条件、模型尺寸和抽油井位置如图4所示。每口井的开采量为200m3/d节点数为23636，总自由度为94544，元件数为43484。在这个特定的例子中，对于每个计算步骤，在作者的计算机（Xeon1230v5，4核，3.4 Ghz）上，Xeon+大约需要两秒钟来组装全局矩阵，另外四秒钟用于稀疏矩阵求解器图5显示了在使用PNDG的最终分析步骤中的总水头和地面沉降幅度。地下水水头最大下降2.68m，地表最大沉降0.216m。尽管所有抽油井的抽水率相同，但远离恒定总水头边界（总水头为零）的区域沉降最大表1含水层系统的输入数据。层 类型kxx=kyykzz泊松厚度顶标高（月/日）（月/日）比（kN/m）（米）（米）阿奎塔尔星3.3e−4 1.1e−4 0.4 2991 18−92含水层8.2 2.73 0.3 56277 24− 27阿奎塔尔星1.1e−3 3.6e−4 0.4 20468 22−514含水层46.8 15.6 0.3 51587 14− 73阿奎塔尔星1.8e−2 6e-3 0.4 48611 39−876含水层59.5 19.8 0.3 106732 13− 129阿奎塔尔星5.7e−7 1.9e−4 0.4 77778 39−1398含水层11.3 3.77 0.3 116801 21− 178九、阿奎塔尔3.7e−3 1.23e−3 0.4 109546 44−19910含水层67.5 22.5 0.3 114638 18− 243阿奎塔尔星2.6e−3 8.67e−4 0.4 131827 38−26112含水层39.5 13.17 0.3 108605 21− 299阿奎塔尔星1.1e−3 3.67e−4 0.4 146751 41−32014含水层6.8 2.27 0.3 71982 59− 361H.T. Pham，W.吕哈克河谷Schuster等人/SoftwareX 9（2019）1519图四、 14层含水层系统和抽采井的几何形状和边界条件图五. 十年后的模拟结果（z坐标放大五倍）;（a）显示总水头（m），（b）显示地面沉降量（m）。5. 影响和结论FEFLOW中的插件程序PSEF+和PSEF程序为分析基于Biot理论的地面沉降问题提供了一个简单、强大、有效的工具。FEFLOW中的许多可用地下水模型现在都可以从PSEF+和PSEF中受益。当需要用这些模型进行地面沉降分析时，只需增加力学场的参数和边界条件。其他的信息都没动过。此外，应用程序的+和PPENS不仅限于地面沉降分析，他们可以用来解决其他问题，共享相同的理论。例如，软土地基的固结分析，或者它们可以扩展到支持热-水力-力学（THM）建模。未来的发展将包括扩展的ESTA+考虑非线性材料，并从其他知名的地 下 水 代 码 ， 如 ， 例 如 ， MODFLOW ， OpenGeoSys ， 和Hydrus。一个独立的版本正在计划中。确认洪范是由德国DAAD的一个项目NaWaM资助的。引用[1]徐燕生，沈世良，任东君，吴宏南.上海市地面沉降影响因素分析：基于战略环境评价的观点。Sustainability2016;8（6）：573.[2]Abidin HZ，Andreas H，Gumilar I，Fukuda Y，Pohan YE，Deguchi T.印度尼西亚雅加达地面沉降及其与城市发展的关系。Nat Hazards 2011;59（3）：1753.[3]作者：Phien-wejN，Giao PH，Nutalaya P.泰国曼谷的地面沉降 工程地质2006;82（4）：187-201。[4]达莉亚O-Z阿德里安O-G墨西哥城附近长期地面沉降的演变：回顾，实地调查和预测模拟。水资源研究2010;46（1）。[5]MinderhoudPSJ，Erkens G，Pham VH，Bui VT，Erban L，Kooi H，等.越南湄公河三角洲地下水开采25年对沉降的影响。环境研究快报2017;12（6）：064006。[6]加洛韦DL，伯比TJ。综述：地下水开采引起的区域地面沉降。HydrogeolJ2011;19（8）：1459-86.[7]沃尔克斯多弗，蒂姆.德国萨克森州东北部的地下水抽取和地面沉降。矿井水环境1999;18（1）：81[8]Leake S，D.加洛韦MODFLOW地下水模型：地下水位含水层的沉降和含水层系统压实包（SUB-WT）用户指南。2007年[9]Biot MA.三维固结理论。J Appl Phys1941;12（2）：155-64.[10]毕奥特多孔材料弹性和固结方程的通解。J Appl Mech1956;23（1）：91-6.[11] 维鲁伊特河PoroElasticity：http://geo.verruijt.net/. 2016年。[12]Intel（R）.英特尔并行工作室xe。2011年。[13]放大图片作者：Jacob B. Eigen V3。2010年。[14]BiotMA.多孔各向异性固体的弹性和固结理论J Appl Phys1955;26（2）：182-5.[15]Biot MA.多孔介质中的变形与声传播力学。J Appl Phys1962;33（4）：1482-98.