[发明专利]一种基于FLAC3D计算矿区地表最大水平变形的方法

权利要求书

1.基于FLAC3D计算矿区地表最大水平变形的方法，其特征在于，包括如下步骤：

S01、使用FLAC3D建立矿区的三维数值力学模型并计算至平衡，利用FLAC3D中的listgp pos和list gp disp输出模型节点的空间坐标与位移信息，并将其合并存储在Excel中，根据所述模型的矿区地表节点的高程坐标z=h，利用Excel的筛选功能，筛选出模型中z坐标为h的地表节点，将矿区地表节点的x坐标和y坐标与x方向位移的信息、矿区地表节点的x坐标和y坐标与y方向位移的信息分别存储为.txt的文本文件；

S02、使用Surfer软件对矿区地表节点文件进行插值：通过Surfer软件读取步骤S01得到的两个文本文件，利用Surfer软件网格选项中的数据功能对矿区地表节点文件进行插值，选取克里金插值方法，以矿区地表节点网格单元的平均尺寸作为插值距离∆x，对步骤S01得到的两个文本文件进行插值，分别得到表示x方向位移与y方向位移的矿区地表节点水平位移点阵，所述水平位移点阵中各节点离散且等间距，并且插值后的结果以数据文件.dat的形式导出；

S03、利用Matlab计算矿区地表最大水平变形值：在Matlab中使用reshape命令将插值后地表节点的x坐标、y坐标、x方向的位移、y方向的位移分别以矩阵的形式存储在Matlab的XM、YM、XDM、YDM四个工作区中，

使矿区地表边界节点的水平位移点阵中各节点在X方向、Y方向、XY方向的应变值、、赋值为0；

根据Matlab的四个工作区中的坐标与位移信息，采用中心差分法计算矿区地表非边界节点的水平位移点阵中每个节点的应变值、、，并将矿区地表边界及非边界节点的水平位移点阵中每个节点的应变值组合为全矿区地表节点应变矩阵strain_tensor，利用Matlab中的eig函数，获取全矿区地表应变矩阵的特征值，并将其存储在对角矩阵B中，并求strain_tensor的特征向量构成列向量存储在矩阵A中；

然后在Matlab中使用命令index_max = find(diag(B) = max(diag(B)))，将矿区地表的最大水平变形值以对角矩阵的形式存储起来，并使用find函数查询该对角矩阵中非零元素的线性索引值并存储在index_max中；

最后使用命令cart2pol(A(1,index_max),A(2,index_max)) * 180 / pi根据最大特征值的线形索引值提取矩阵A中对应的特征向量的笛卡尔坐标，并把特征向量的笛卡尔坐标系下表示的分量转化为弧坐标表达的极角，再变成角度，进而实现矿区地表最大水平变形方向的求解；

其中，函数diag提取对角矩阵B的对角元组成特征值对角矩阵D，即不在对角线上元素全为0的方阵，函数max用于找出特征值对角矩阵D的最大特征值，命令diag(B) = max(diag(B))用于将矿区地表最大的水平变形值替换特征值对角矩阵D中的数值，cart2pol函数把二维笛卡尔坐标转换为对应的极坐标，矩阵A包括两列构成特征向量的坐标值；

S04、利用Matlab编程语言根据矿区地表节点的水平位移点阵和各节点对应的最大水平变形值绘制成等值线图，并保存Surfer软件可识别的.grd格式；

S05、在Surfer软件中打开格式为.grd的矿区地表最大水平变形等值线图，利用Surfer软件地图选项中的等值线功能对其进行后处理，从而得出表征矿区地表最大水平变形大小的云图。

2.根据权利要求1所述的基于FLAC3D计算矿区地表最大水平变形的方法，其特征在于，采用中心差分法计算矿区地表非边界节点的水平位移点阵中第i行,j列节点的、、的步骤如下：

=[XDM(i,j+1) - XDM(i,j-1)] / [XM(i,j+1) - XM(i,j-1)];

= [YDM(i+1,j) - YDM(i-1,j)]/ [YM(i+1,j) - YM(i-1,j)];

= [(XDM(i+1,j) - XDM(i-1,j) ] / [YM(i+1,j) - YM(i-1,j)]+ [YDM(i,j+1)- YDM(i,j-1)] / [XM(i,j+1) - XM(i,j-1)];

XM(i,j)、YM(i,j) 、XDM(i,j)、YDM (i,j)分别为矿区地表节点的水平位移点阵中第i行,j列节点的横坐标、纵坐标、x方向的位移、y方向的位移。

3.根据权利要求1所述的基于FLAC3D计算矿区地表最大水平变形的方法，其特征在于，所述的将矿区地表边界及非边界节点的水平位移点阵中每个节点的应变值组合为全矿区地表节点应变矩阵strain_tensor的步骤包括：将每个节点的应变值、、分别存储在空矩阵strain_xx、strain_yy、strain_xy中，在Matlab使用strain_tensor =[strain_xx, 0.5 * strain_xy; 0.5 * strain_xy, strain_yy]命令集成所述的空矩阵得到全矿区地表节点应变矩阵strain_tensor。