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

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


背景技术：

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


技术实现要素：

3.本发明针对在不断地下水开采过程中，地下水位下降和岩土体上覆荷载增大引发 的地面沉降问题，提供一种考虑岩土体上覆荷载和地下水开采双重效应下的渗流场和 应力场全耦合模拟方法，该方法能够真实有效评价预测未来地面沉降演化过程。
4.为了实现上述目前，本发明采用以下技术方案：一种荷载和地下水开采引发地面 沉降和地裂缝三维变参数全耦合模拟计算方法，包括以下步骤：
5.步骤1：根据研究区范围、水文地质条件和工程地质条件，确定研究区的范围、 高程、垂向地层厚度和岩性条件，确定各地层的分层特征，建立研究区三维水文地质 和工程地质概念模型；
6.步骤2：对水文地质和工程地质概念模型进行有限元网格剖分，将概念模型剖分 为结构单元网，建立有限元网格模型，对有限元网格模型进行参数分区，导出各分区 有效网格的节点、单元信息，分别选取有限元网格模型中的各类边界单元和节点，分 别导出边界单元和节点信息，制作出不同的边界类型；
7.步骤3：将步骤2中的不同类型的边界赋予不同的属性值，将模型的边界参数写 入到导出的各个边界的单元和节点信息中，生成边界条件；
8.步骤4：利用步骤2中导出的所有有效节点信息，将所有节点的空间坐标位置和 初始水位写入节点信息中，生成初始水位条件；
9.步骤5：利用步骤2中导出的模型所有有效单元和节点信息，筛选出所有地下水 开
采中的抽水井和观测井的单元和节点信息，将抽水井概化为排水子结构，并将抽水 井流量写入抽水井的节点信息中，最终将观测井水位值写入观测井节点信息中，生成 抽水井条件、观测井条件；
10.步骤6：利用步骤2中导出的有效节点信息，筛选出承载荷载的岩土体节点信息， 将所有节点的附加荷载写入节点信息中，生成岩土体荷载条件，其中，各节点的地面 荷载由上覆荷载与附加应力系数共同确定；
11.步骤7：利用步骤2中导出的所有有效节点和参数分区信息，将各个参数分区中 具有非线性变化的渗透系数、孔隙度、给水度/贮水率数据输入参数分区中，生成水文 地质条件，将各个参数分区具有非线性变化的弹性模量、泊松比写入参数分区中，生 成工程地质条件；
12.步骤8：将步骤3至步骤7中的文件导入以biot固结理论为基础的地下水系统 (gws)，完成数值模型的建立，利用伽辽金加权余量法对时间离散，采用预处理共 轭梯度法(pcg)对已建立的数值模型进行有限元求解；
13.步骤9：将步骤8中计算得到的地下水位变化与实际抽水试验的水位和水文长观 孔水位进行对比：若计算值与观测值的差值大于水位变化值的10％，返回步骤8调整 参数直到结果符合，保存模型；最终，利用调整完成的模型预测渗流场的水位和地面 沉降量。
14.进一步的，其特征在于：步骤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系统在考虑岩土体上覆荷载和 孔隙水压力减小的共同作用时，进行如下改进：
[0029][0030]
式中：g为剪切模量，e为弹性模量，ν为泊松比；w
x
,w
y
,w
z
为x,y,z 方向上的位移；k
x
,k
y
,k
z
分别为三个方向上的渗透系数，γ
w
为水的重度；γ为岩土重度；
[0031]
数值模型的地应力、位移和孔隙水压力的初始条件如下：
[0032][0033]
式中：σ
x
和σ
z
为土体的初始水平向和垂向应力；γ为岩土重度；z为计算点深度； k1为静止侧压力系数；w0(x,y,z)为研究区域内初始位移；u0(x,y,z)为研究区域内初始孔 隙水压力；
[0034]
数值模型的孔隙水压力边界条件如下：
[0035]
u(x,y,z)|
γ1
＝u
s
[0036]
式中：u(x,y,z)为各点孔隙水压力；u
s
为水头边界γ1上的已知孔隙水压力；
[0037]
数值模型的流量的边界条件如下：
[0038][0039]
式中：k为渗透系数；h为水头；n1为法向向量；q
l
为边界γ2上的已知单位面积 流量；
[0040]
数值模型的自由面边界条件如下：
[0041]
u＝γ
w
z1[0042][0043]
式中：u为自由面孔隙水压力；γ
w
为水的重度；z1为自由面所在高程；q为通过自 由面边界单位面积流量；μ为土体给水度；θ为自由面外法线方向与垂线的交角；
[0044]
数值模型的位移边界条件如下：
[0045][0046]
式中，为位移边界γ4上三个方向的已知位移。
[0047]
进一步的：步骤2中还包括研究区网格剖分完成后，检查岩土体网格的质量，合 格则开启步骤3，不合格需重新剖分直至合格。
[0048]
有益效果：与现有技术相比，本发明的技术方案具有以下技术效果：
[0049]
1、本发明的计算方法利用biot固结理论、地下水渗流理论，建立了地面沉降三 维流固耦合数值模型。计算方法充分考虑了岩土体上覆荷载和地下水开采耦合作用下 对岩土体骨架性质的影响，地面沉降过程中岩土体变形规律，不同于定参数中渗流场 和应力场的单一叠加，从地下水开采引发的岩土体性质改变的角度实现了应力场与渗 流场的完全耦合。本发明为研究地下水开采过程中的渗流场演化和地面沉降演化规律 以及流固耦合作用机理提供依据。
[0050]
2、本发明的计算方法精细描述了岩土体性质和渗流场规律，有效保证了岩土体上 覆荷载和地下水开采时的质量守恒，剖分网格质量较高，提高了模拟计算结果的精度。
孔水位进行对比：若计算值与观测值差距较大，返回步骤2调整参数，直到结果在误 差范围内，保存模型；若计算值与观测值在误差范围内，保存此模型。
[0061]
步骤5：用步骤4中的模型预测岩土体上覆荷载和地下水开采共同作用下引发的 渗流场的演化情况和地面沉降的演化情况。
[0062]
本发明的计算方法通过研究地下水开采和上覆荷载情况下的水文地质和工程地质 条件的变化，建立了渗流场
‑
应力场三维全耦合数值计算的模型，预测了地下水开采过 程中的地面沉降量，为有效控制地面沉降提供可参考的依据。
[0063]
以上所述仅为本发明的较佳实施例，并不限制本发明，凡在本发明的精神和原则 范围内所做的任何修改、等同替换和改进，均应包含在本发明的保护范围之内。