一种荷载和地下水开采引发地面沉降三维变参数全耦合模拟计算方法技术

本发明专利技术提出一种荷载和地下水开采引发地面沉降和地裂缝三维变参数全耦合模拟计算方法，所解决的是在土体上覆荷载和地下水开采双重作用下的地面沉降量和地裂缝的计算问题。本发明专利技术根据地下水开采的实际背景，建立地质概念模型。基于Biot固结理论，实现了渗流场与应力场全耦合，构建了用于预测地面沉降量以及地裂缝发育的三维全耦合数值模型。本计算方法考虑岩土体压缩过程中的机理，并耦合弹性、黏弹性和黏塑性的本构模型，结合岩土体参数的变化，将渗流场与应力场紧密联系，可保证能同时模拟地面荷载和地下水动态下降情况下的地面沉降量和地裂缝延伸量，为研究地下水开采过程中的渗流场演化和地面沉降演化规律以及流固耦合作用机理提供了依据。作用机理提供了依据。作用机理提供了依据。

【技术实现步骤摘要】
一种荷载和地下水开采引发地面沉降三维变参数全耦合模拟计算方法


[0001]本专利技术涉及地质环境保护及修复领域，特别是涉及一种荷载和地下水开采引发地 面沉降和地裂缝三维变参数全耦合模拟计算方法技术。

技术介绍

[0002]地面沉降是在自然条件或者人类活动影响下，地下岩土体骨架压缩而在一定范围 内呈现地面高程缓慢降低的一种现象，其严重影响着高层房屋、公路、高铁等建筑物 的安全性。地裂缝也是地面沉降的衍生现象。在众多影响因素中，地下水过度开采和 开展工程的荷载增大是引发地面沉降的最主要因素。目前针对地面沉降的预测方法在 大类上主要是基于地下水开采量和地面沉降量统计的回归预测模型和基于地下水和岩 土体骨架相耦合的部分流固耦合模型。但是，以上两种模型在研究地下水开采对地面 沉降的影响有着一定的优势，当存在岩土体上覆荷载时，其精度难以保证。传统的流 固耦合模型主要考虑的是渗流场和应力场的叠加效应，并没有充分考虑荷载和地下水 运移过程中岩土体参数变化、本构模型变化的三维全耦合模型。因此，为了更准确的 预测地下水开采过程中渗流场的演化和地面沉降量的变化，应利用考虑岩土参数变化 和本构模型变化的全耦合模型进行地面沉降量预测。

技术实现思路

[0003]本专利技术针对在不断地下水开采过程中，地下水位下降和岩土体上覆荷载增大引发 的地面沉降问题，提供一种考虑岩土体上覆荷载和地下水开采双重效应下的渗流场和 应力场全耦合模拟方法，该方法能够真实有效评价预测未来地面沉降演化过程。
[0004]为了实现上述目前，本专利技术采用以下技术方案：一种荷载和地下水开采引发地面 沉降和地裂缝三维变参数全耦合模拟计算方法，包括以下步骤：
[0005]步骤1：根据研究区范围、水文地质条件和工程地质条件，确定研究区的范围、 高程、垂向地层厚度和岩性条件，确定各地层的分层特征，建立研究区三维水文地质 和工程地质概念模型；
[0006]步骤2：对水文地质和工程地质概念模型进行有限元网格剖分，将概念模型剖分 为结构单元网，建立有限元网格模型，对有限元网格模型进行参数分区，导出各分区 有效网格的节点、单元信息，分别选取有限元网格模型中的各类边界单元和节点，分 别导出边界单元和节点信息，制作出不同的边界类型；
[0007]步骤3：将步骤2中的不同类型的边界赋予不同的属性值，将模型的边界参数写 入到导出的各个边界的单元和节点信息中，生成边界条件；
[0008]步骤4：利用步骤2中导出的所有有效节点信息，将所有节点的空间坐标位置和 初始水位写入节点信息中，生成初始水位条件；
[0009]步骤5：利用步骤2中导出的模型所有有效单元和节点信息，筛选出所有地下水 开
采中的抽水井和观测井的单元和节点信息，将抽水井概化为排水子结构，并将抽水 井流量写入抽水井的节点信息中，最终将观测井水位值写入观测井节点信息中，生成 抽水井条件、观测井条件；
[0010]步骤6：利用步骤2中导出的有效节点信息，筛选出承载荷载的岩土体节点信息， 将所有节点的附加荷载写入节点信息中，生成岩土体荷载条件，其中，各节点的地面 荷载由上覆荷载与附加应力系数共同确定；
[0011]步骤7：利用步骤2中导出的所有有效节点和参数分区信息，将各个参数分区中 具有非线性变化的渗透系数、孔隙度、给水度/贮水率数据输入参数分区中，生成水文 地质条件，将各个参数分区具有非线性变化的弹性模量、泊松比写入参数分区中，生 成工程地质条件；
[0012]步骤8：将步骤3至步骤7中的文件导入以Biot固结理论为基础的地下水系统 (GWS)，完成数值模型的建立，利用伽辽金加权余量法对时间离散，采用预处理共 轭梯度法(PCG)对已建立的数值模型进行有限元求解；
[0013]步骤9：将步骤8中计算得到的地下水位变化与实际抽水试验的水位和水文长观 孔水位进行对比：若计算值与观测值的差值大于水位变化值的10％，返回步骤8调整 参数直到结果符合，保存模型；最终，利用调整完成的模型预测渗流场的水位和地面 沉降量。
[0014]进一步的，其特征在于：步骤4中的井的处理方法为排水子结构法，排水子结构 的节点分为4层并重新编号，内层节点编号为9～12，次内层节点编号为5～8，外层节 点编号为1～4，最外层节点编号为13～16，9～12为已知水头的排水孔结点，1～8结点 为次内层与外层的未知水头结点，这些结点基于井流量形成内部传导矩阵；最外层的 13～16结点则是建立了连接内部与外部结构的外部传导矩阵，内部传导矩阵与外部传 导矩阵组成的总体传导矩阵[K
s
]及相应流量列阵q
s
，最终，排水子结构总传导矩阵 及相应流量列阵为：
[0015][0016]式中：[K
ii
]为排水子结构内部，未知水头结点间的传导矩阵；[K
ib
]、[K
bi
]分别 为内部未知水头的结点与出口结点之间的传导矩阵；[K
bb
]为排水子结构出口处，未知 水头结点间的传导矩阵；{h
i
}为排水子结构的内部已知结点水头列阵；{h
b
}为排水子 结构出口的结点水头列阵；{q
i
}为已知结点水头对内部未知结点水头的流量；{q
b
}为 已知结点水头对出口未知结点水头的流量。
[0017]进一步的，其特征在于：步骤6中对岩土体上覆荷载进行了精确的描述，引入了 附加应力系数，最终的岩土体上覆荷载计算公式为：
[0018]P＝P
d
·
a
sr
[0019]式中：P
d
为地面建筑荷载；a
sr
为附加应力系数,计算方式如下：
[0020][0021]式中，l为矩形荷载的长边，b为矩形荷载的短边，x、y、z为计算点坐标值，a
sr
为 附加应力系数。
[0022]进一步的，其特征在于，步骤7中生成的孔隙度、渗透系数、弹性模量和泊松比 参数呈非线性变化特征，参数变化方程如下：
[0023][0024][0025][0026][0027]式中：n为孔隙度；n0为初始孔隙度；ΔP为孔隙水压力变化量；K
s
为多孔介质骨 架固体颗粒的体积弹性压缩模量；ε
v
为体应变；K为渗透系数；K0为初始渗透系数； E
t
为切线弹性模量；P
a
为大气压强，k为围压在100KPa时的初始弹性模量，m为弹性 模量与固结压力曲线的斜率，在1.0～0.2之间；c为土体的黏聚力；为土体的内摩擦 角；R
f
为破坏比；G1为围压为100KP
a
时的初始泊松比；D、F为土体试验常数；ε1为 轴向应变；σ1为最大主应力；σ3为最小主应力。
[0028]进一步的：步骤8中对包含Biot固结理论的GWS系统在考虑岩土体上覆荷载和 孔隙水压力减小的共同作用本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种荷载和地下水开采引发地面沉降和地裂缝三维变参数全耦合模拟计算方法，其特征在于：包括以下步骤：步骤1：根据研究区范围、水文地质条件和工程地质条件，确定研究区的范围、高程、垂向地层厚度和岩性条件，确定各地层的分层特征，建立研究区三维水文地质和工程地质概念模型；步骤2：对水文地质和工程地质概念模型进行有限元网格剖分，将概念模型剖分为结构单元网，建立有限元网格模型，对有限元网格模型进行参数分区，导出各分区有效网格的节点、单元信息，分别选取有限元网格模型中的各类边界单元和节点，分别导出边界单元和节点信息，制作出不同的边界类型；步骤3：将步骤2中的不同类型的边界赋予不同的属性值，将模型的边界参数写入到导出的各个边界的单元和节点信息中，生成边界条件；步骤4：利用步骤2中导出的所有有效节点信息，将所有节点的空间坐标位置和初始水位写入节点信息中，生成初始水位条件；步骤5：利用步骤2中导出的模型所有有效单元和节点信息，筛选出所有地下水开采中的抽水井和观测井的单元和节点信息，将抽水井概化为排水子结构，并将抽水井流量写入抽水井的节点信息中，最终将观测井水位值写入观测井节点信息中，生成抽水井条件、观测井条件；步骤6：利用步骤2中导出的有效节点信息，筛选出承载荷载的岩土体节点信息，将所有节点的附加荷载写入节点信息中，生成岩土体荷载条件，其中，各节点的地面荷载由上覆荷载与附加应力系数共同确定；步骤7：利用步骤2中导出的所有有效节点和参数分区信息，将各个参数分区中具有非线性变化的渗透系数、孔隙度、给水度/贮水率数据输入参数分区中，生成水文地质条件，将各个参数分区具有非线性变化的弹性模量、泊松比写入参数分区中，生成工程地质条件；步骤8：将步骤3至步骤7中的文件导入以Biot固结理论为基础的地下水系统(GWS)，完成数值模型的建立，利用伽辽金加权余量法对时间离散，采用预处理共轭梯度法(PCG)对已建立的数值模型进行有限元求解；步骤9：将步骤8中计算得到的地下水位变化与实际抽水试验的水位和水文长观孔水位进行对比：若计算值与观测值的差值大于水位变化值的10％，返回步骤8调整参数直到结果符合，保存模型；最终，利用调整完成的模型预测渗流场的水位和地面沉降量。2.根据权利要求1所述的一种荷载和地下水开采引发地面沉降和地裂缝三维变参数全耦合模拟计算方法，其特征在于：步骤4中的井的处理方法为排水子结构法，排水子结构的节点分为4层并重新编号，内层节点编号为9～12，次内层节点编号为5～8，外层节点编号为1～4，最外层节点编号为13～16，9～12为已知水头的排水孔结点，1～8结点为次内层与外层的未知水头结点，这些结点基于井流量形成内部传导矩阵；最外层的13～16结点则是建立了连接内部与外部结构的外部传导矩阵，内部传导矩阵与外部传导矩阵组成的总体传导矩阵[K
s
]及相应流量列阵q
s
，最终，排水子结构总传导矩阵及相应流量列阵为：
式中：[K
ii
]为排水子结构内部，未知水头结点间的传导矩阵；[K
ib
]、[K
bi
]分别为内部未知水头的结点与出口结点之间的传导矩阵；[K
bb
]为排水子结构出口处，未知水头结点间的传导矩阵；{h
i
}为排水子结构的内部已知结点水头列阵；{h
b
}为排水子结构出口的结点水头列阵；{q
i
}为已知结点水头对内部未知结点水头的流量；{q
b
}为已知结点水头对出口未知结点水头的流量。3.根据权利要求1...

【专利技术属性】
技术研发人员：骆祖江，付延玲，王鑫，李兆，张也，成磊，李子义，
申请(专利权)人：河海大学，
类型：发明
国别省市：