写在前面：
新开个人公众号一枚，长期驻守。更新随心，忙时沉淀，闲时落笔。往后会陆续分享个人科研点滴，承蒙厚爱，欢迎大家关注指教✨
手动改参数、通宵跑模型、一张张截图？这个案例教你打通Python与FLAC3D（PFC/UDEC同样适用），构建顶级期刊级别的自动化数值试验流水线！
各向异性岩体中的隧道开挖
在这个案例中，我们需要模拟一个深埋在柱状玄武岩中的圆形隧道。由于岩体存在明显的正交节理（各向异性），不同方向的节理倾角会对隧道的最终变形产生巨大影响。
我们的目标是：探究第3组节理倾角 0° 变化到180°（每隔15°为一组，共13个模型）时，隧道拱顶下沉量的变化规律。
Python充当“总指挥”(批量控制)，我们不需要写13个FLAC3D脚本，只需要用Python的Itasca模块写一个循环大脑。它会自动遍历所有角度，把参数“塞”进 FLAC3D 里，并自动命名保存。
it.command("python-reset-state false")
这个隐藏的神仙指令，保证了 FLAC3D 在重复执行
Model new时，绝不会清空我们的Python循环和内存环境！
自动化数据校验与批量出图 (后处理解放双手)
模型算完了，怎么提取数据？难道要一个个点开看吗？我们使用了一套极其优雅的回归测试（Regression Testing）工作流：
通过自动化脚本，程序会将tunnel-000.sav到tunnel-180.sav依次唤醒。
精准提取：利用内置FISH函数[gp.disp(gp.near(0,0,1))->z]精准抓取隧道拱顶的Z向下沉量。
误差自检：配合-test check real指令，程序会自动对比理论值，误差超过 1e-6即报警。
静默出图：自动导出所有角度下的位移（disp）和塑性状态（state）高清位图，尺寸和视角完全一致，直接满足 SCI 论文的插图标准！
将Python强大的生态（控制逻辑、数据处理）与FLAC3D/PFC深厚的力学内核结合，不仅能让你彻底告别机械重复的调参截屏，更能极大提升科研数据的产出效率和准确度。
不同角度下位移（第3组节理倾角0°变化到180°，每隔15°）
不同角度下塑性区
tunnel.dat;---------------------------------------------------------------------; numerical solution for a long tunnel in pre-stressed; anisotropic material - plane strain;---------------------------------------------------------------------;model newzone deletemodel large-strain offzone create radial-cylinder size 1 1 50 50 ratio 1 1 1 1.1 point 1 20 0 0 ... point 2 0 0.2 0 point 3 0 0 20 dimension 1 1 1 1zone reflect normal -1 0 0 origin 0 0 0zone reflect normal 0 0 -1 origin 0 0 0;zone cmodel assign columnar-basaltzone property flag-matrix-plastic falsezone property young=60e3 poisson=0.30zone property space-1=0.1 space-2=0.1 space-3=0.1zone property stiffness-normal-1=1200e3 stiffness-normal-2=1200e3 stiffness-normal-3=300e3 zone property stiffness-shear-1=1200e3 stiffness-shear-2=1200e3 stiffness-shear-3=50.4e3zone property dip-1 [dip1] dip-direction-1 90zone property normal-x-2 0 normal-y-2 1 normal-z-2 0zone property dip-3 [dip3] dip-direction-3 90zone property joint-cohesion-1 5.0 joint-friction-1 5 joint-dilation-1 0 joint-tension-1 1e20zone property joint-cohesion-2 1e20 joint-friction-2 0 joint-dilation-2 0 joint-tension-2 1e20zone property joint-cohesion-3 5.0 joint-friction-3 10 joint-dilation-3 0 joint-tension-3 1e20;zone initialize stress-xx -8zone initialize stress-yy -10zone initialize stress-zz -12zone gridpoint fix velocity-z range position-z -20.1 -19.9zone gridpoint fix velocity-z range position-z 19.9 20.1zone gridpoint fix velocity-x range position-x -20.1 -19.9zone gridpoint fix velocity-x range position-x 19.9 20.1zone gridpoint fix velocity-y;zone history displacement-x position 1 0 0zone history velocity-x position 1 0 0zone history unbalanced-force-y position 1 0 0zone history unbalanced-force-z position 1 0 0history interval 20model solvemodel save [savefile]
check.datmodel restore 'tunnel-000'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0012558 tolerance 1e-6 windows-test check real [gp.disp(gp.near(0,0,1))->z] -0.0012558 tolerance 1e-6 linuxplot 'disp' export bitmap filename "disp-000.png"plot 'state' export bitmap filename "state-000.png"model restore 'tunnel-015'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0012441 tolerance 1e-6plot 'disp' export bitmap filename "disp-015.png"plot 'state' export bitmap filename "state-015.png"model restore 'tunnel-030'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0011793 tolerance 1e-6 windows-test check real [gp.disp(gp.near(0,0,1))->z] -0.0011793 tolerance 1e-6 linuxplot 'disp' export bitmap filename "disp-030.png"plot 'state' export bitmap filename "state-030.png"model restore 'tunnel-045'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0010735 tolerance 1e-6plot 'disp' export bitmap filename "disp-045.png"plot 'state' export bitmap filename "state-045.png"model restore 'tunnel-060'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0009492 tolerance 1e-6 windows-test check real [gp.disp(gp.near(0,0,1))->z] -0.0009492 tolerance 1e-6 linuxplot 'disp' export bitmap filename "disp-060.png"plot 'state' export bitmap filename "state-060.png"model restore 'tunnel-075'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0008431 tolerance 1e-6plot 'disp' export bitmap filename "disp-075.png"plot 'state' export bitmap filename "state-075.png"model restore 'tunnel-090'-test check real [gp.disp(gp.near(0,0,1))->z] -0.000797 tolerance 1e-6plot 'disp' export bitmap filename "disp-090.png"plot 'state' export bitmap filename "state-090.png"model restore 'tunnel-105'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0008446 tolerance 1e-6plot 'disp' export bitmap filename "disp-105.png"plot 'state' export bitmap filename "state-105.png"model restore 'tunnel-120'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0009498 tolerance 1e-6plot 'disp' export bitmap filename "disp-120.png"plot 'state' export bitmap filename "state-120.png"model restore 'tunnel-135'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0010734 tolerance 1e-6plot 'disp' export bitmap filename "disp-135.png"plot 'state' export bitmap filename "state-135.png"model restore 'tunnel-150'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0011787 tolerance 1e-6 windows-test check real [gp.disp(gp.near(0,0,1))->z] -0.0011787 tolerance 1e-6 linuxplot 'disp' export bitmap filename "disp-150.png"plot 'state' export bitmap filename "state-150.png"model restore 'tunnel-165'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0012433 tolerance 1e-6 windows-test check real [gp.disp(gp.near(0,0,1))->z] -0.0012433 tolerance 1e-6 linuxplot 'disp' export bitmap filename "disp-165.png"plot 'state' export bitmap filename "state-165.png"model restore 'tunnel-180'-test check real [gp.disp(gp.near(0,0,1))->z] -0.0012543 tolerance 1e-6 windows-test check real [gp.disp(gp.near(0,0,1))->z] -0.0012543 tolerance 1e-6 linuxplot 'disp' export bitmap filename "disp-180.png"plot 'state' export bitmap filename "state-180.png"
master.pyimport itasca as itit.command("python-reset-state false")dip3List = [0,15,30,45,60,75,90,105,120,135,150,165,180]outfile = 'tunnel-'for dip3 in dip3List : it.command("model new") it.fish.set('dip3', dip3) it.fish.set('dip1', dip3 + 90) it.fish.set('savefile', outfile + '{:03d}'.format(dip3) + '.sav') it.command("program call 'tunnel.dat' ")it.clear_save_variables()
master.datprogram call "master.py"program call "check.dat"
