发明内容

本发明针对在不断地下水开采过程中，地下水位下降和岩土体上覆荷载增大引发的地面沉降问题，提供一种考虑岩土体上覆荷载和地下水开采双重效应下的渗流场和 应力场全耦合模拟方法，该方法能够真实有效评价预测未来地面沉降演化过程。

为了实现上述目前，本发明采用以下技术方案：一种荷载和地下水开采引发地面沉降和地裂缝三维变参数全耦合模拟计算方法，包括以下步骤：

步骤1：根据研究区范围、水文地质条件和工程地质条件，确定研究区的范围、 高程、垂向地层厚度和岩性条件，确定各地层的分层特征，建立研究区三维水文地质 和工程地质概念模型；

步骤2：对水文地质和工程地质概念模型进行有限元网格剖分，将概念模型剖分为结构单元网，建立有限元网格模型，对有限元网格模型进行参数分区，导出各分区 有效网格的节点、单元信息，分别选取有限元网格模型中的各类边界单元和节点，分 别导出边界单元和节点信息，制作出不同的边界类型；

步骤3：将步骤2中的不同类型的边界赋予不同的属性值，将模型的边界参数写 入到导出的各个边界的单元和节点信息中，生成边界条件；

步骤4：利用步骤2中导出的所有有效节点信息，将所有节点的空间坐标位置和 初始水位写入节点信息中，生成初始水位条件；

步骤5：利用步骤2中导出的模型所有有效单元和节点信息，筛选出所有地下水 开采中的抽水井和观测井的单元和节点信息，将抽水井概化为排水子结构，并将抽水 井流量写入抽水井的节点信息中，最终将观测井水位值写入观测井节点信息中，生成 抽水井条件、观测井条件；

步骤6：利用步骤2中导出的有效节点信息，筛选出承载荷载的岩土体节点信息，将所有节点的附加荷载写入节点信息中，生成岩土体荷载条件，其中，各节点的地面 荷载由上覆荷载与附加应力系数共同确定；

步骤7：利用步骤2中导出的所有有效节点和参数分区信息，将各个参数分区中 具有非线性变化的渗透系数、孔隙度、给水度/贮水率数据输入参数分区中，生成水文 地质条件，将各个参数分区具有非线性变化的弹性模量、泊松比写入参数分区中，生 成工程地质条件；

步骤8：将步骤3至步骤7中的文件导入以Biot固结理论为基础的地下水系统(GWS)，完成数值模型的建立，利用伽辽金加权余量法对时间离散，采用预处理共 轭梯度法(PCG)对已建立的数值模型进行有限元求解；

步骤9：将步骤8中计算得到的地下水位变化与实际抽水试验的水位和水文长观孔水位进行对比：若计算值与观测值的差值大于水位变化值的10％，返回步骤8调整 参数直到结果符合，保存模型；最终，利用调整完成的模型预测渗流场的水位和地面 沉降量。

进一步的，其特征在于：步骤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}为 已知结点水头对出口未知结点水头的流量。

进一步的，其特征在于：步骤6中对岩土体上覆荷载进行了精确的描述，引入了 附加应力系数，最终的岩土体上覆荷载计算公式为：

P＝P_d·a_sr

式中：P_d为地面建筑荷载；a_sr为附加应力系数,计算方式如下：

式中，l为矩形荷载的长边，b为矩形荷载的短边，x、y、z为计算点坐标值，a_sr为 附加应力系数。

进一步的，其特征在于，步骤7中生成的孔隙度、渗透系数、弹性模量和泊松比 参数呈非线性变化特征，参数变化方程如下：

式中：n为孔隙度；n₀为初始孔隙度；ΔP为孔隙水压力变化量；K_s为多孔介质骨 架固体颗粒的体积弹性压缩模量；ε_v为体应变；K为渗透系数；K₀为初始渗透系数； E_t为切线弹性模量；P_a为大气压强，k为围压在100KPa时的初始弹性模量，m为弹性 模量与固结压力曲线的斜率，在1.0～0.2之间；c为土体的黏聚力；为土体的内摩擦 角；R_f为破坏比；G₁为围压为100KP_a时的初始泊松比；D、F为土体试验常数；ε₁为 轴向应变；σ₁为最大主应力；σ₃为最小主应力。

进一步的：步骤8中对包含Biot固结理论的GWS系统在考虑岩土体上覆荷载和 孔隙水压力减小的共同作用时，进行如下改进：

式中：G为剪切模量，E为弹性模量，ν为泊松比；w_x,w_y,w_z为x,y,z 方向上的位移；k_x,k_y,k_z分别为三个方向上的渗透系数，γ_w为水的重度；γ为岩土重度；

数值模型的地应力、位移和孔隙水压力的初始条件如下：

式中：σ_x和σ_z为土体的初始水平向和垂向应力；γ为岩土重度；z为计算点深度； K₁为静止侧压力系数；w₀(x,y,z)为研究区域内初始位移；u₀(x,y,z)为研究区域内初始孔 隙水压力；

数值模型的孔隙水压力边界条件如下：

u(x,y,z)|_Γ1＝u_s

式中：u(x,y,z)为各点孔隙水压力；u_s为水头边界Γ₁上的已知孔隙水压力；

数值模型的流量的边界条件如下：

式中：K为渗透系数；H为水头；n₁为法向向量；q_L为边界Γ₂上的已知单位面积 流量；

数值模型的自由面边界条件如下：

u＝γ_wz₁

式中：u为自由面孔隙水压力；γ_w为水的重度；z₁为自由面所在高程；q为通过自 由面边界单位面积流量；μ为土体给水度；θ为自由面外法线方向与垂线的交角；

数值模型的位移边界条件如下：

式中，为位移边界Γ₄上三个方向的已知位移。

进一步的：步骤2中还包括研究区网格剖分完成后，检查岩土体网格的质量，合 格则开启步骤3，不合格需重新剖分直至合格。

有益效果：与现有技术相比，本发明的技术方案具有以下技术效果：

1、本发明的计算方法利用Biot固结理论、地下水渗流理论，建立了地面沉降三 维流固耦合数值模型。计算方法充分考虑了岩土体上覆荷载和地下水开采耦合作用下 对岩土体骨架性质的影响，地面沉降过程中岩土体变形规律，不同于定参数中渗流场 和应力场的单一叠加，从地下水开采引发的岩土体性质改变的角度实现了应力场与渗 流场的完全耦合。本发明为研究地下水开采过程中的渗流场演化和地面沉降演化规律 以及流固耦合作用机理提供依据。

2、本发明的计算方法精细描述了岩土体性质和渗流场规律，有效保证了岩土体上覆荷载和地下水开采时的质量守恒，剖分网格质量较高，提高了模拟计算结果的精度。