[发明专利]一种采水型地裂缝数值模拟方法

权利要求书

1.一种采水型地裂缝数值模拟方法，其特征在于，所述方法采用钝化断裂带模型表征地裂缝，所述钝化断裂带模型将土体破坏分为完整、微裂隙产生和宏观裂缝产生三个阶段，通过两个临界应变和划分所述三个阶段；其中f_t为土体的抗拉强度，E为土体拉伸模量，G_f为裂缝能量释放率临界值，w_c为与土体单元尺寸有关的断裂带宽度；所述模拟方法具体步骤包括：

(1)建立开采地下水引起的地层移动模型，并根据模拟的具体情况设定初始条件和边界条件；

(2)对模拟区域进行有限元网格剖分，得到地层移动模型的有限元的离散形式：

其中K为整体劲度矩阵，K′为整体结点孔隙水压力和结点位移的耦合矩阵，为整体渗透矩阵，Δw和Δp分别为整体结点位移增量和整体结点孔隙水压力增量，F为整体结点外荷载，w_t-1和p_t-1分别为上一时间步结束时的整体结点位移和整体结点孔隙水压力，Q为源汇项，Δt为时间步长增量；

(3)考虑初始条件和边界条件，求解步骤(2)中的方程组，得到当前时间步的整体结点位移增量和整体孔隙水压力增量，从而得到各单元的应变增量和应力增量，进而得到各单元的累计应变和应力；

(4)判断是否有单元小主应变达到临界应变ε₀，若达到则该单元进入软化阶段，在软化阶段，应力随着应变的增加而减小；判断是否有单元小主应变达到临界应变ε_f，若达到则该单元产生破坏，出现垂直于小主应变方向的宏观裂缝，在宏观裂缝出现后，对该单元在裂缝法线方向上的刚度进行折减，从而修正该破坏单元的应力应变矩阵；在单元出现宏观裂缝后，将单元的应力应变矩阵修正为：

[D]^cs＝[D]^s-[D]^s[N]([D]^c+[N]^T[D]^s[N])^-1[N]^T[D]^s

其中，

[D]^cs为含裂缝土体的应力应变矩阵，[D]^s为未开裂土体的应力应变矩阵，[D]^c为裂缝的应力应变矩阵，[N]为以垂直和平行于裂缝方向建立的局部坐标系与整体坐标系之间的坐标转换矩阵，n、s、t为裂缝处局部坐标系的三个坐标轴方向，n垂直于裂缝方向，s和t平行于裂缝方向，Dⁿⁿ为裂缝的法向刚度，D^ns、D^nt为裂缝的切向刚度；法向刚度和切向刚度的计算式为：

其中，E、G分别为土的拉伸模量和剪切模量，E_s为应力应变软化段直线的斜率，为负值，β为剪力传递系数；

(5)返回步骤(2)，进行下一时间步的计算，直至所有时间步模拟结束。

2.根据权利要求1所述的一种采水型地裂缝数值模拟方法，其特征在于，所述步骤(3)中，各单元的应变增量{Δε}＝[B]{Δw}^e，应力增量{Δσ}＝[D]{Δε}，累计应变{ε}_t＝{ε}_t-1+{Δε}，累计应力{σ}_t＝{σ}_t-1+{Δσ}；其中[D]为单元的应力应变矩阵，[B]为单元应变矩阵，{Δw}^e为单元结点位移增量列向量，{ε}_t-1、{ε}_t分别为t-1时刻和t时刻的单元应变，{σ}_t-1、{σ}_t分别为t-1时刻和t时刻的单元应力。

3.根据权利要求1所述的一种采水型地裂缝数值模拟方法，其特征在于，在单元进入软化阶段后，应力随着应变的增加而减小，具体方法是：在全过程曲线上找到单元当前应变对应的应力，并与单元当前应力求差得到该单元的不平衡应力；如果该不平衡应力大于收敛的阈值，则将该单元应力调整到全过程曲线上对应的应力，并且将不平衡应力对应的荷载反向加载到该单元的结点上，计算相应的应变和应力；重复这一过程，直至单元的不平衡应力不超过收敛的阈值。