【后处理探讨】FLAC实现计算结果对称复制

FLAC3D小技巧

6月前浏览7870

摘要：本文对FLAC3D实现计算结果对称复制进行了探讨。

1.问题来源

在③群中，有小伙伴提出了如下问题：

图1-1 问题来源

基于以上背景，本文对FLAC3D计算结果的对称复制进行了讨论，其中计算案例来自《厚壁圆筒问题径向、切向应力求解及Extra云图显示》。

案例背景
厚壁圆筒问题径向、切向应力求解及Extra云图显示

2.正文及代码

2.1 初始计算

在《厚壁圆筒问题径向、切向应力求解及Extra云图显示》一文中，我采取建立完整圆筒模型进行计算。在本文中，我仅建立1/4圆环进行计算，并将结算结果对称复制为整个圆环。

图2-1 计算模型

计算模型位移、径向应力、切向应力代码如下：

import itasca as itit.command("python-reset-state false")# 模型参数r1 = 50  # 内半径/mmr2 = 100  # 外半径/mmth = 5  # 轴向厚度/mms_r = 30  # 径向网格数量s_th = 1  # 轴向网格数量s_t = 30  # 环向网格数量E = 2e5  # 弹性模量/MPav = 0.3  # 泊松比p = -10  # 内孔压力/MPa
# 建模计算it.command("""          model new          model largestrain off          zone create cylindrical-shell point 0(0,0,0) point 1({0},0,0) point 2(0,{1},0) point 3(0,0,{0}) ...                                       point 4({0},{1},0) point 5(0,{1},{0}) point 8({2},0,0) ...                                       point 9(0,0,{2}) point 10({2},{1},0) point 11(0,{1},{2}) size {3} {4} {5}          zone group 'zone1'          zone cmodel assign elastic          zone property young {6} poisson {7}          zone face skin         zone face apply velocity-normal 0 range group 'west2' or 'bottom'          zone face apply stress-normal {8} range group 'west1'          model solve          """.format(r2, th, r1, s_r, s_th, s_t, E, v, p))

# 圆心坐标xc = 0.0zc = 0.0
# 遍历求解各单元径向、切向应力for z in it.zone.list():    xz = z.pos_x()    zz = z.pos_z()    r = math.sqrt((xz-xc)**2+(zz-zc)**2)    sin = zz/r    cos = xz/r    sxx = z.stress()[0][0]    szz = z.stress()[2][2]    sxz = z.stress()[2][0]    sig_r = sxx*cos**2+szz*sin**2+2*sxz*cos*sin    sig_t = sxx*sin**2+szz*cos**2-2*sxz*cos*sin    z.set_extra(1, sig_r)    z.set_extra(2, sig_t)    z.set_extra(3, z.stress()[0][0])

模型位移、径向应力、切向应力计算结果如下所示。与上文计算结果对比可知：采用1/4模型进行计算，其结果与采用完整模型计算几乎无差别。

图2-2 位移云图（1/4模型）

图2-3 径向应力云图（1/4模型）

图2-4 切向应力云图（1/4模型）

2.2 单元（zone）计算结果的对称复制

对于单元（zone），在计算完成后对模型进行对称复制，其部分结算结果会随着单元的复制而自动复制，如应变增量（strain-increment）。但结果与建立全模型时的结果不一致，见下图

1/4模型-对称复制单元

全模型

图2-5 XX应变云图结果对比

但有一部分计算结果压根就不会被复制，如下图的stress-xx云图：

1/4模型-对称复制单元

全模型

图2-6 XX应力云图

在我查看Extra云图时，发现定义的Extra云图均随着模型的对称复制而自动复制，数值仅有细微偏差（个人猜测是由于对称前后的单元插值环境不同以及边界条件不同导致的）。

图2-7 径向应力云图（1/4模型-对称复制单元）

图2-8 切向应力云图（1/4模型-对称复制单元）

于是，我有了如下猜测：只需要把想要复制的单元计算结果转移到Extra云图，即可实现单元计算结果的对称复制。以图2-6~图2-7为例，在计算完成后，将XX应变与XX应力写入到Extra云图后再对模型进行对称复制，代码如下：

# 遍历求解各单元xx单元应变、应力for z in it.zone.list():    z.set_extra(3,z.strain()[0][0])    z.set_extra(4, z.stress()[0][0])

运行上述代码并对称复制模型后，得到的xx单元应变、应力云图如下所示：

图2-9 XX应变云图（1/4模型-对称复制单元-Extra云图）

图2-9 XX应变云图（1/4模型-对称复制单元-Extra云图）

由图2-8~图2-9可以看出，采用将单元计算结果导入Extra云图的方式可以实现单元计算结果的对称复制。

说明1：为了方便，本文代码以Python编写，但也可以采用fish进行实现，相关函数框架可查询往期发布，包括但不限于：

遍历及其框架教学
FLAC3D 6.0新手向——遍历详解

指针获取教学
FLAC3D6.0各类指针获取详解

说明2：云图标题可通过Legend中的Title进行修改，如下所示：

图2-10 修改Extra云图标题

2.3 网格点（gridpoint）计算结果的对称复制

对于网格点（gridpoint），在计算完成后对模型进行对称复制，其结算结果和Extra均不会随着单元的复制而自动复制，以位移为例，结果见下图。

1/4模型-Extra云图

1/4模型-合位移云图

图2-11 合位移云图

因此，为了实现网格点计算结果的对称复制，就需要人为赋值。要实现这个功能，就需要我们回忆一下一点点高中知识：如何求解一个点关于一个面的对称点坐标。

对于绝大多数情况而言，对称平面都是YOZ/XOY/XOZ面或与其平行的平面。如左右对称复制模型时，对称平面为YOZ平面，此时假设原模型中的一点P(x1,y1,z1)，显然其关于YOZ平面（X=x0）的对称点P'为（x1-2*x0,y1,z1）。

图2-12 对称点坐标关系示意图

因此，我们只需要计算出对称点坐标，并将原模型各个网格点的计算结果赋给对称点即可。效果如下所示。两种计算方式下位移值有所不同，初步考虑是全模型、1/4模型的边界条件不同导致的。

1/4模型-对称复制单元-Extra云图

全模型

图2-13 合位移云图

为了使上述代码更具普适性，考虑模型关于非YOZ/XOY/XOZ平面对称的情况（比如本例建立1/8模型进行计算）。此时，同样仅需要采用Python求取原网格点关于对称平面的对称点坐标然后赋值即可。

来源：FLAC3D小技巧

位移差？应力差？利用fish语言实现结果文件间运算。

在使用FLAC3D计算完成后，有时需要对两个sav文件进行操作，以获得两个结果文件间的孔隙压力差云图、位移差云图以及应力差云图等。由于FLAC在restore结果文件时会自动清空当前内存里的所有fish变量，因此无法直接实现采用fish函数提取变量值再相减的思路。为解决这一问题，需要将第一次提取的变量值写出到本地文件，restore第二个结果文件后，再读入本地文件并执行相关操作，操作时可选用array，table以及map存储并写出文件。本文基于以上思路，使用array完成两个结果文件间的孔隙压力差云图。======案例演示======案例模型为20*20*20的立方体，分组情况见图1。首先对整个模型施加初始孔压1e3，保存结果为a.f3sav，孔压云图见图2；然后将a组的孔压归零，b组孔压固定为1e3，保存结果为b.f3sav，孔压云图见图3。案例演示的目的是得到结果a减去结果b的孔隙压力差云图，结果见图4。图1 模型分组图2 结果a的孔压云图图3 结果b的孔压云图图4 孔隙压力差云图命令流如下： ;建立模型并储存结果文件model newzone create brick size 20 20 20 group 'a'zone group 'b' range position-x 5 15 position-z 5 15zone gridpoint initialize pore-pressure 1e3 model save 'a'zone gridpoint initialize pore-pressure 0 range group 'a'model save 'b';设置文件读写参数&获取数组所需sizemodel restore 'a.f3sav'fish def setup io_read = 0 io_write = 1 io_ascii = 1 arr_size = gp.numend@setup;提取结果a中的孔隙压力值并写出到文件fish def pp_out array id1(arr_size) array pp1(arr_size) n = 0 loop foreach gp1 gp.list n = n+1 id1(n) = gp.id(gp1) pp1(n) = gp.pp(gp1) endloop status = file.open('id_list',io_write,io_acsii) status = file.write(id1,string(arr_size)) status = file.close status = file.open('pp_list',io_write,io_acsii) status = file.write(pp1,string(arr_size)) status = file.closeend@pp_out;设置文件读写参数&获取数组所需sizemodel restore 'b.f3sav'fish def setup io_read = 0 io_write = 1 io_ascii = 1 arr_size = gp.numend@setup;读取结果a的孔隙压力值并与结果b的孔隙压力值相减fish def pp_minus array id(arr_size) array pp(arr_size) status = file.open('id_list',io_read,io_acsii) status = file.read(id,arr_size) status = file.close status = file.open('pp_list',io_read,io_acsii) status = file.read(pp,arr_size) status = file.close loop i(1,arr_size) gp = gp.find(id(i)) gp.pp(gp) = pp(i) - gp.pp(gp) endloopend@pp_minusmodel save 'c'注：命令流中文件读写参数的含义请自行查阅help文件。 ======总结======进行结果文件间的操作时，需要将前一个结果文件的数据写出到本地文件，以避免清空。进行文件的写出操作时，可采用array/table/map，写出文件的格式可以为txt/csv等，也可以像本文一样不指定后缀。应注意的是，采用csv格式存储时，会有存储行数的上限值，据反馈为104万行，本文命令流已通过800W+行数据检验，且10秒内完成孔压相减。读者可在此框架内进行修改，完成不同结果文件间不同变量的加减乘除操作。来源：FLAC3D小技巧

【后处理探讨】FLAC实现计算结果对称复制

1.问题来源

2.正文及代码

附件

位移差？应力差？利用fish语言实现结果文件间运算。