FLAC3D模型某平面应力的suffer后处理

pdf版本在百度文库就有：

http://wenku.baidu.com/link?url=y714SR3EK9CO0WpYqK-iz1MCsXFrgCekl1PBqrKH31ylc_gWpYsjOQ17hyT5hwWnNqxuvNGVyL4OeLCgAsjtpOO-0IoA7f3bJruXSBiy2ka

FLAC^3D模型应力的Surfer可视化后处理

王志强¹，朱丙龙¹，马奇飞¹，司永达¹，王云云¹，赵位位²

（1.中国矿业大学（北京）资源与安全工程学院，北京100083；2. 山东省济宁市兖州区新兖镇牛楼社区杨庄煤矿，272000）

摘要：FLAC^3D是用于工程力学计算的显式有限差分程序，主要用来解决有限单元法不能企及的大变形问题、复杂非线性行为和不稳定系统等问题，被广泛应用于岩土工程、土木工程和采矿工程领域，但其应力的形象表达等后处理环节还有待提高，而Surfer软件能够将三维数据简便快速的生成二维、三维图形。针对FLAC^3D模型分析结果可视化不丰富的问题，采取用FLAC^3D内置fish语言提取模型相应变量数据并以Surfer可以读取的格式组织起来，用Surfer进行模型应力的可视化后处理。通过具体实例表明，给出的fish语言代码能够方便正确的提取数据，Surfer软件生成的底板应力可视化图形能够将应力大小与高低对应起来，较形象的表达FLAC^3D分析数据。

关键词： FLAC^3D；Surfer；可视化；fish语言；底板应力

中图分类号：TD315       文献标识码：A

Visible post-process of stresses in FLAC^3D model with Surfer

WANG Zhi-qiang, ZHU Bing-long, MA Qi-fei, SI Yong-da, WANG Yun-yun, ZHAO Wei-wei

(1.School of Resource and Safety Engineering, China University of Mining and Technology, Beijing 100083, China; 2. Yangzhuang coal mine in Niulou community Xinyan town Yanzhou county Jining city in Shandong province, 272000)

Abstract：FLAC^3D is a explicit finite difference program，which were used to calculate engineering mechanics. It also be used in geotechnical engineering, civil engineering and mining engineering field to solve the large deformation problem, the complicated nonlinear behavior and unstable system, etc. which finite element method cannot solve, but the post-processing such as image expression of stress need to be improved. Surfer software could generate 2d and 3d graphics from 3d data simply and rapidly. For the problem of little abundant visualization of FLAC^3D model analysis results, using fish language of FLAC^3D extracted the variable data that Surfer needed and organized that data which Surfer can read, Surfer could process that stresses of FLAC^3D model visible. The given fish code could extract stresses easily and correctly through a specific example, and the size of stress could correspond to height in the stress graphics generated by Surfer, which expressed the FLAC^3D analysis data more readily relate.

Kewwords：FLAC^3D； Surfer； visualization； fish language; floor stresses

 

 

科学计算可视化是发达国家20世纪80年代后期提出并发展起来的一个新的研究领域。有限差分计算结果数据的图形化显示，即FLAC^3D模型应力的Surfer后处理属于科学计算的可视化，^[1]就是即所谓的FLAC^3D分析的后处理。目前FLAC^3D后处理主要通过FLAC^3D软件自带后处理功能、TECPLOT软件处理和MATLAB处理出图。FLAC^3D自带后处理功能处理最方便，但所出云图不形象也不易具体量化，切片云图不能同时出现两个，如图1所示，功能较简单。

 

图1 FLAC3D所出应力切片云图示例

Fig.1 the FLAC3D plane of szz stress nephogram

TECPLOT软件处理后有各种应力和位移变量可以显示，可以同时做多个切片，如图2所示，功能丰富，但一般仅能以颜色和数值区分数据大小。

 

图2 tecplot所出应力多切片示例

Fig.2 the contour line of tecplot plane of szz stress

MATLAB软件虽然能够通过图形高低表示数据大小，如图3所示，但软件本身较大，对计算机硬件要求较高，另外将层面非矩形FLAC^3D模型数据放入矩阵过程复杂，使用范围受限。^[2]

 

图3 MATLAB光滑处理后工作面底板szz应力图

Fig.3 the floor szz mess graph of coal face in the MATLAB after smoothing

找到一种像MATLAB一样能够以网格高低表示应力大小，且不要求FLAC^3D模型为长方体、所出图形同样直观形象、操作较简便的方法将有利于科研人员分析数据。

1  Surfer可视化处理FLAC^3D分析数据

Surfer可视化处理FLAC^3D分析数据时，首先需要知道Surfer能处理怎样的数据格式，然后据此提取FLAC^3D模型中相应数据，最后Surfer软件读入数据并处理。具体流程图如图4所示。

 

图4 Surfer可视化FLAC3D分析数据流程图

Fig.4 the flow diagram of Surfer visualize the FLAC3D analysis data

1.1  Surfer读入的数据格式

Surfer软件是一套二维和三维图形绘制软件，相比MATLAB而言具有占用内存小，处理计算量大的优点。^[3]。Surfer强大的插值功能和绘制图件能力，已经使它成为处理xyz三维数据的首选软件，它能迅速地将离散点的测量数据通过插值转换为连续的数据曲面，进一步简单、正确地把xyz数据转换成丰富多彩的等值线图、粘贴图、影像图、地貌晕渲图、矢量地图、线框图和表面图等。^[4-5]

Surfer能够读入文本文档中存放的成列的数据。比如第一列为x坐标，第二列为y坐标，第三列为相应x、y坐标处的属性值，在此为应力值，第一行的三列分别为数据的题头。^[6]例如：

Xposition(m) Yposition(m) SZZ(Mpa）

2.50E+00    5.00E+00   7.13E+00

7.50E+00    5.00E+00   7.21E+00

1.25E+01    5.00E+00   7.37E+00

……

当然，文本文档中数据也可超过三列，在Surfer中可以自由选取不同列显示不同项目。

1.2  FLAC^3D模拟结果的数据提取

fish语言作为FLAC^3D的内嵌语言，是一种嵌入式编程语言，实现了使用者对软件的完美控制。^[7]当FLAC^3D计算结果保存为（*.sav）后，可以通过FLAC^3D内置的fish语言提取想要Surfer表达的变量，并写入文本文档中。

在FLAC^3D中，各种应力保存在单元体中，位移保存在节点中。通过fish语言保留字z_xcen(pnt)、z_ycen(pnt)、z_szz(pnt)可以提取pnt指向的单元体的x、y坐标和szz应力，通过pnt = zone_head和pnt = z_next(pnt)遍历不同单元体。^[8]提取FLAC^3D模型单元体各种应力的代码如下，FLAC^3D软件3.0版本读入（*.sav）之后在读入以下列代码为内容的文本文档即可在同目录下生成文本文档1xyand9stresses2sufer8.dat供Surfer读入处理出图。

; Flac3d Mesh to suffer

;fist rest *.sav use flac3d，then call this file

def ini_mesh2suf

;;;;Edit the tec_range to set plot range

    command

;;;;define the range of plot in use of two plane

range name Lay2 plane nor 0 0 1 ori 0 0 20.5 below

range name Lay1 plane nor 0 0 1 ori 0 0 19.5 above

        ;group haitang range group 6 any group 7 any

        ;model null range group haitang not

;ran name tec_range group haitang  Lay2 Lay1; the group and the range between two plane codetermine that whether the zone plot or not

;ran name tec_range group haitang;the group control whether the zone plot or not

;ran name tec_range Lay2 Lay1;the range between two plane codetermine that whether the zone plot or not

        ran name tec_range Lay2 Lay1 ;

    endcommand

array buf1(1)

        file1='1xyand9stresses2sufer8.dat'

end

;;;; If plotit = 1, plot the zone

def plot_test

    plotit = 0;;;;don't plot by default

    if z_model(pnt) = 'NULL' then;

        plotit = 0

    endif   

    if inrange('tec_range',pnt) = 1 then;;;;if the zone in the group tec_range,the inrange will equivalent to 1，and the zone will plot

        plotit = 1

    endif

end

def GET_STRESS_

  status=open(file1,1,1)

    buf1(1)=' X(m) Y(m) SIG1(MPa) SIG2(MPa) SIG3(MPa) SXX(MPa) SYY(MPa) SZZ(MPa) SXY(MPa) SYZ(MPa) SZX(MPa)'

status=write(buf1,1)

    pnt = zone_head

    loop while pnt # null

buf1(1) = ' '

          plot_test

    if plotit = 1 then

buf1(1) = buf1(1) + string(z_xcen(pnt)) +' '

buf1(1) = buf1(1) + string(z_ycen(pnt)) +' '

       buf1(1) = buf1(1) + string(z_sig1(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_sig2(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_sig3(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_sxx(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_syy(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_szz(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_sxy(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_syz(pnt)*(-0.000001)) + ' '

       buf1(1) = buf1(1) + string(z_sxz(pnt)*(-0.000001)) + ' '

status=write(buf1,1)

endif

pnt = z_next(pnt)

    end_loop

end

def Wait

ini_mesh2suf

GET_STRESS_

status=close

ii = out('Successfully Write Zone Stress Data Into Suffer Format File: '+ file1 )

end

Wait

上述平面夹逼法获取某层面全部单元体各种应力代码中，主要包含了ini_mesh2suf、plot_test、GET_STRESS_和Wait四个自定义函数，其中Wait为程序代码入口即主函数，主函数被首先执行，之后依据顺序依次调取执行。注意上述代码仅适用于FLAC^3D3.0版本，FLAC^3D5.0版本对fish语言语法做了一定修改，需要修改才能正确适用于5.0版本，其余代码同样如此。^[9]

在主函数中首先执行初始化函数ini_mesh2suf，初始化函数的作用是限定FLAC^3D模型数据提取范围，定义缓存数组buf1以及命名数据存放文本文档的名称。提取范围一般为某一层面单元体，层面或平或曲。对于采矿工程中水平煤层开挖模拟时，底板应力可以直接使用上述代码中方法，即用两平面夹逼限定提取单元体范围，此范围要将单元体重心包括在两夹逼平面内。若煤层存在一定角度，只需要修改由点法式定义的两平面即可，在此建议建立FLAC^3D倾斜煤层模型时先在AutoCAD中画出平面图，模型各关键点及各层岩层法向量即可轻松获取。由于实际建模中模型层面可能会有起伏，两平面未必能精确限定所提取单元体范围，若某层面单元体可以由组（group）精确限定，则可以尝试上述代码中被注释掉的组的方法。需要注意的是status=write(buf1,1)需要放在“loop-end loop”循环体内，否则单元体数目较多时就会出现错误，文献[2]中代码需要如此改进。

接下来是plot_test函数，其作用是判定单元体是否在要提取数据的单元体范围之内，若在此范围之内，则提取数据共Surfer显示。GET_STRESS_是获取单元体应力函数，由于FLAC^3D中压应力为负且单位为Pa，为了在Surfer中显示方便，将九种应力都乘以-0.000001变为正值且单位为MPa。^[10]

若初始化函数ini_mesh2suf中两种方法仍不能够精确限定所要提取单元体范围，由于FLAC^3D建模时单元体编号按层累加，所以可以根据单元体编号来获取相应单元应力。此方法首先需要用鼠标双击所需提取层面单元体最靠近坐标原点的表面以获取单元体编号id，根据编号用pnt=find_zone(id)获得单元体指针pnt以便遍历层面内所有单元。根据z_id()函数，由指针pnt可获得单元体编号z_id(pnt)，同理，可双击要获取层面末单元体即距离原点最远的单元体表面获得末单元体编号，然后根据首末单元体编号范围遍历单元体应力。改动后代码如下：

def ini_mesh2suf

       array buf1(1)

        file1='1idxyand9stresses2sufer8.dat'

end

def GET_STRESS_

  status=open(file1,1,1)

    buf1(1)=' Zoneid X(m) Y(m) Z(m) SIG1(MPa) '

       status=write(buf1,1)

    pnt=find_zone(id);id（43513） is the first zone's id which in the range

   loop while z_id(pnt)<=ID;ID（49728） is the end zone's id which in the range

       buf1(1) = ' '

       buf1(1) = buf1(1) + string(z_id(pnt)) +' '

       buf1(1) = buf1(1) + string(z_xcen(pnt)) +' '

       buf1(1) = buf1(1) + string(z_ycen(pnt)) +' '

       buf1(1) = buf1(1) + string(z_zcen(pnt)) +' '

           buf1(1) = buf1(1) + string(z_sig1(pnt)*(-0.000001)) + ' '

       status=write(buf1,1)

       pnt = z_next(pnt)

    end_loop

end

def Wait

ini_mesh2suf

GET_STRESS_

status=close

ii = out('Successfully Write Zone Stress Data Into Suffer Format File: '+ file1 )

end

Wait

此处仅输出了单元体编号、对应编号单元体重心的x、y坐标和第一主应力四列，若想输出其他应力变量，仅需参照平面夹逼法获取单元体各种应力代码添加所需应力即可。若需要某层面中某列测线上单元体数据以供数据化出图，仅需在限定范围处加上相应的if条件语句进行限定即可。

1.3  Surfer软件读入数据并处理

首先打开Surfer8.0软件，菜单栏网格-->数据…-->在打开对话框中选择要网格化的单元体坐标及应力数据文本文档（1xyand9stresses2sufer8.dat），点击打开会弹出网格化数据对话框，如图5所示。数据列X选择列A即x坐标列，数据列Y选择列B即y坐标列，数据列Z选择列C即要表达坐标xy处的属性，比如需要最大主应力SIG1(MPa)分布图，那么就按图5所示选择，若需要竖直压应力SZZ分布图，那么选择相应的SZZ(MPa)即可。

 

图5 Surfer8.0网格化数据对话框

Fig.5 the grid data dialog of Surfer8.0

此时查看数据可以得到如图6所示的对话框，有图6可以看出A、B两列为单元体坐标，C~K为单元体九种应力。

 

图6 单元体数据文本文档读入sufer

Fig.6 zone's data that stored in text documents be readed in sufer

另一处需要设置的选择为网格化方法。网格化就是把以xyz数据表示的文件格式表示的、通常是不规则分布的原始数据点，经过数学处理，构筑一个规则的空间矩形网格的过程。原始数据的不规则分布会造成缺失数据的“空洞”，网格化则用外推或内插的算法填充了这些“空洞”。^[11]Surfer共有精确插值和平滑插值两大类共12种插值方法^[12]。根据FLAC^3D模型提取信息数据量大、数据点密集、要求精确的特点，初步选定加权反距离法、带线性三角剖分法、最小曲率法、径向基本函数法中的多重二次曲面法和最近邻点插值法。具体插值方法的最终选取根据有效性评价进行。

选择z列数据和网格化方法后，点击确定，即可在同目录下生成网格文件（1idxyand9stresses2sufer8.grd）。然后点击右侧工具栏中线框图按钮会弹出对话框，选择网格文件就会生成三维网格图。

2  工程实例

主要通过屯兰矿18203工作面瓦斯聚集区和煤柱合理性留设的模拟首采面推进120m实例演示上述提取FLAC^3D模型层面单元体应力供Surfer可视化处理的fish代码接口程序的实用性，并说明Surfer可视化后处理的优点。

2.1  FLAC^3D模型的建立

取屯兰矿18203工作面处围岩建模，首采面工作面长度设为100m，沿走向（y轴）推进，接续面沿倾向推进，工作面长度为150m。因为屯兰矿18203工作面开采煤层为近水平煤层，为了简化模拟建立水平岩层模型进行模拟，x为工作面倾向方向，y为工作面推进的走向方向，z为重力方向。8号煤埋深284m，模型最上部埋深254m，上覆岩层平均破断角按80度应留18.51m边界，此处留30m边界，可以基本消除模型边界效应，模型尺寸为296m(x)×210m(y)×131.5m(z)，整个模型内单元体边长在x和y方向分别为4m和2.5m。计算模型共划分有304584个单元，318750个节点，各层名称、层序见图7所示，由于层数较多，左面窗口未能全部列出。

 

图7 计算机数值模拟模型

Fig.7 Computer numerical simulation model

采用摩尔库伦准则并施加边界条件达到平衡。图8为工作面岩y轴推进120m时底板应力切面。

 

图8 z=20m水平切面szz应力图

Fig.8 The SZZ stress of horizontal plane when z equivalent 20 with sketch

可见FLAC3D自带的后处理功能切出的平面云图直观形象性不强。

2.1  Surfer可视化处理有效性评价

FLAC^3D运行接口程序，会在同目录下生成1xyand9stresses2sufer8.dat，Surfer读入之，分别按下面五种插值方法将数据网格化，为了计算方便需要在交叉网格中将使用的样点值和估计值输出，在网格化时都取1000个数据点，分别计算平均估计误差百分比PAEE、相对均方差RMSE、均方根预测误差RMSPE，误差越小的插值方法就越好。^[13]另外并不是网格点越多越好，网格过密,则网格分析值连续性好，但有时会导致失真。^[14]

表1不同插值方法的有效性评价指标计算

Table5.1 The evaluation index result of effectiveness of different interpolation methods

模型	PAEE(%)	RMSE	PMSPE
加权反距离法	12.35	0.0505	1.0631
带线性三角剖分法	1.71	0.007765	0.4024
最小曲率法	3.14	0.01252	0.5312
径向基本函数法	0.74	0.003129	0.2575
最近邻点插值法	10.04	0.04251	0.9530

由表1可以看出在用Surfer可视化处理应力图时径向基本函数法最好，其次依次是带线性三角剖分法、最小曲率法、最近邻点插值法，加权反距离法最次。

 

图9加权反距离法显示工作面底板SZZ

Fig.9 The SZZ stress of floor using Inverse Distance to a Pow-er Function

 

图10 带线性插值的三角剖分方法显示工作面底板SZZ

Fig.10 The SZZ stress of floor using Triangulation with Linear Interpolation Function

 

 

图11 最小曲率法显示工作面底板SZZ

Fig.11 The SZZ stress of floor using Minimum Curvature Function

 

图12 径向基本函数法中的多重二次曲面法显示工作面底板SZZ

Fig.12 The SZZ stress of floor using Multiquadratic of Radial Basis Function

 

图13 最近邻点插值法显示工作面底板SZZ

Fig.13 The SZZ stress of floor using Nearest Neighbor Function

以上都是对工作面底板竖直压应力SZZ的可视化处理，将应力大小与位置高低对应了起来，应力分布直观形象。图13中网格最低处为采空区位置，周围最高处为采空区四周压应力升高区域，而升高区外侧为原始应力区，应力升高区的范围及应力集中系数也可有直观印象。

 

 

 

 

图14 x=102测线底板单元各种应力图

Fig.14 all stress of zone which center at x direction equivalent to 102

用if条件语句和单元体重心x坐标提取函数z_xcen(pnt)限定x坐标范围提取x=102这条测线上单元体应力，并用origin处理，如图14所示，由坐标数值可知30m边界是足够的，工作面前方40m左右支承压力为零，应力增高系数约为0.73。由此可见，结合Surfer和origin出图进行分析即形象又精确。当然，origin可以在一张图中处理不同推进度下某一种应力图，然后对比不同推进度下的多个Surfer三维网格图，这样更有实际意义。有效性评价主要是评价样点估计值对样点值的逼真程度，选取哪种插值方法与样点值的密度、多少、样点值的属性值（例如本论文中的各种应力）水平和连续性有关，对于fish语言提取的其他8个应力，这几个性质几乎相同，仅应力大小水平和连续性有些差别，有效性评价同样适用。

3  结论与展望

1、论文中给出了FLAC^3D和Surfer软件数据传输的程序接口即一段fish代码。通过工程实例证明此代码能够将FLAC^3D单元体信息提取并供Surfer读入出图处理。并且提取单元体信息的fish代码提供了两平面夹逼、组和单元体编号三种限定范围的提取方法。

2、通过Surfer可视化处理有效性评价，最终选定径向基本函数法中的多重二次曲面法最适合工作面底板应力SZZ的可视化处理。直观经验来看此方法确实最贴近实际且美观，且同样适用于单元其他应力的显示。

3、需要提取某条测线上单元体应力以供origin或Excel量化处理时，只需在根据编号提取单元体应力的fish语言中限定某条测线范围。

4、研究地表沉陷需要提取某条测线上节点位移以供origin或Excel处理，只需在根据编号提取单元体应力的fish语言中限定某条测线范围及改变相应的提取函数，变提取单元应力为提取节点位移。可以预见，本论文思路同样适用于地表沉陷研究时用fish语言提取模型节点应力并用Surfer可视化后处理，并进行有效性评价选出最适合的插值方法。^[15]

参考文献：

[1]      唐泽圣.三维数据场可视化.北京:清华大学出版社,1999；1-8

[2]      王志强. FLAC3D模型interface法向应力的MATLAB处理[J] 中州煤炭,2015,03;37-44.

[3]      杨柏宁. 水下地形三维建模Surfer与MATLAB比较分析[J]. 上海国土资源,2012,01:75-78.

[4]      宋明艺,张春灌. 借助Surfer软件实现快速绘制平面等值线图[J]. 工程地球物理学报,2009,02:244-246.

[5]      王旭春,黄福昌,张怀新,张连贵. 开采沉陷可视化工程分析设计系统[J]. 矿山测量,2002,03:22-25+4.

[6]      赵方,周荣福,耿咏梅,汪丽媛. Surfer与AutoCAD结合在绘制等值线中的应用[J]. 地理空间信息,2011,02:46-48+8.

[7]      孙书伟，林杭，任连伟. FLAC3D在岩土工程中的应用. 中国水利水电出版社,2011,6;14.

[8]      彭文斌. FLAC3D实用教程.机械工业出版社,2013,1;69-70.

[9]      美国Itasca公司. FLAC^3D3.10 UserManual2004

[10]  陈育民,俆鼎平. FLAC/FLAC3D基础与工程实例 中国水利水电出版社;95-103.

[11]  张弛. 基于Surfer软件的地形图快速绘制方法[J]. 水利技术监督,2009,03:41-43.

[12]  CHEN Huan-huan, LI Xing, DING Wen-xiu. Twelve Kinds of Gridding Methods of Surfer 8. 0 in Isoline Draw-ing[J]. Chinese Journal of Engineering Geophysics, 2007,4(1): 52-57. (陈欢欢,李星,丁文秀. Surfer 8.0等值线绘制中的十二种插值方法[J].工程地球物理学报, 2007,4(1): 52-57.)

[13]   武俊红,汪云甲. 基于Surfer的煤矿等值线空间插值方法有效性评价[J]. 中国矿业,2007,01:108-110.

[14]  连志鸾. Surfer二次开发实现加密雨量图自动显示与输出[J]. 气象科技,2006,02:220-224.

[15]  马天勤. Surfer软件在开采沉陷可视化中的应用[J]. 矿业工程,2008,03:48-50.

 

 

 

基金项目：中央高校基本科研业务费专项资金项目(2011QZ06);国家自然科学基金青年科学基金项目(51404270);国家级大学生创新训练项目(201211413014)

作者简介：王志强（1980－），男，内蒙古呼伦贝尔人，副教授，博士，2009年毕业于中国矿业大学（北京），现任资源与安全工程学院采矿工程系教师。