Flac问题讲解.doc

分差法
岩土工程结构的数值解是建立在满足基本方程（平衡方程、几何方程、本构方程）和边界条件下推导的。由于基本方程和边界条件多以微分方程的形式出现，因此，将基本方程近假发改用差分方程（代数方程）表示，把求解微分方程的问题改换成求解代数方程的问题，这就是所谓的分差法。差分法由来已久，但差分法需要求解高阶代数方程组，只有在计算机的出现，才使该法得以实施和发展。
FLAC3D（Fast Lagrangian Analysis of Continua）由美国Itasca公司开发的。目前，FLAC有二维和三维计算程序两个版本，二维计算程序V3.0以前的为DOS版本，V2.5版本仅仅能够使用计算机的基本内存64K），所以，程序求解的最大结点数仅限于2000个以内。1995年，FLAC2D已升级为V3.3的版本，其程序能够使用护展内存。因此，大大发护展了计算规模。FLAC3D是一个三维有限差分程序，目前已发展到V2。1版本。
FLAC3D的输入和一般的数值分析程序不同，它可以用交互的方式，从键盘输入各种命令，也可以写成命令（集）文件，类似于批处理，由文件来驱动。因此，采用FLAC程序进行计算，必须了解各种命令关键词的功能，然后，按照计算顺序，将命令按先后，依次排列，形成可以完成一定计算任务的命令文件。 FLAC3D是二维的有限差分程序FLAC2D的护展，能够进行土 质、岩石和其它材料的三维结构受力特性模拟和塑性流动分析。调整三维网格中的多面体单元来拟合实际的结构。单元材料可采用线性或非线性本构模型，在外力作用下，当材料发生屈服流动后，网格能够相应发变形和移动（大变形模式）。FLAC3D采用的显式拉格朗日算法和混合-离散分区技术能够非常准确发模拟材料的塑性破坏和流动。由于无须形成刚度矩阵，因此，基于较小内存空间就能够求解大范围的三维问题。FLAC3D采用ANSI C++语言编写的。

FLAC3D有以下几个优点：
1 对模拟塑性破坏和塑性流动采用的是“混合离散法“。这种方法比有限元法中通常采用的“离散集成法“更为准确、合理。

2 即使模拟的系统是静态的，仍采用了动态运动方程，这使得FLAC3D在模拟物理上的不稳定过程不存在数值上的障碍。

3 采用了一个“显式解“方案。因此，显式解方案对非线性的应力-应变关系的求解所花费的时间，几互与线性本构关系相同，而隐式求解方案将会花费较长的时间求解非线性问题。面且，它没有必要存储刚度矩阵，这就意味着；

采用中等容量的内存可以求解多单元结构；
模拟大变形问题几互并不比小变形问题多消耗更多的计算时间，因为没有任何刚度矩阵要被修改。 当然，它也存在以下几个不足之处：

1 对于线性问题的求解，FLAC3D比有限元程序运行得要慢；因此，当进行大变形非线性问题或模拟实际可能出现不稳定问题时，FLAC3D是最有效的工具。

2 用FLAC3D求解时间取决于最长的自然周期和最短的自然周期之比。但某些问题对模型是无效的。

接上面的介绍：Flac3D中为岩土工程问题的求解开发了特有的本构模型，总共包含了10种材料模型：

1. 开挖模型null

2. 3个弹性模型（各向同性，横观各向同性和正交各向同性弹性模型）

3 6个塑性模型（Drucker-Prager模型、Morh-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型和修正的cam粘土模型）。

Flac3D网格中的每个区域可以给以不同的材料模型，并且还允许指定材料参数的统计分布和变化梯度。 还包含了节理单元，也称为界面单元，能够模拟两种或多种材料界面不同材料性质的间断特性。节理允许发生滑动或分离，因此可以用来模拟岩体中的断层、节理或摩擦边界。FLAC3D中的网格生成器gen，通过匹配、连接由网格生成器生成局部网格，能够方便地生成所需要的三维结构网格。还可以自动产生交叉结构网格（比如说相交的巷道），三维网格由整体坐标系x,y,z系统所确定，不同于FLAC程序是由行列方式确定。这就提供了比较灵活的产生和定义三维空间参数。

FLAC3D做计算分析的一般步骤：

与大多数程序采用数据输入方式不同，FLAC采用的是命令驱动方式。命令字控制着程序的运行。在必要时，尤其是绘图，还可以启动FLAc用户交互式图形界面。为了建立FLAC计算模型，必须进 行以下三个方面的工作：
1. 有限差分网格

2. 本构特性与材料性质

3. 边界条件与初始条件

完成上述工作后，可以获 得模型的初始平衡状态，也就是模拟开挖前的原岩应力状态。然后，进行工程开挖或改变边界条件来进行工程的响应分析，类似于FLAC的显式有限差分程序的问题求解。与传统的隐式求解程序不同，FLAC采用一种显式的时间步来求解代数方程。进行一系列计算步后达到问题的解。

在FLAC中，达到问题所需的计算步能够通过程序或用户加以控制，但是，用户必须确定计算步是否已经达到问题的最终的解。
最后，进行结果的分析与总结，这就是用户自己的事情了，呵呵，
在手册中的图中有显示shell单元
但我自己做的时候不知用何命令才能在网格图中将的shell，beam等单元显示出来，高手指点！！
plot add sel geom

另外，还可在上面命令的后面加：cid=on node=off .....
再问一下，shell是否只可以加在网格的表面？怎么在网格内部加不上？
可以加上了，但是不知道厚度t如何表现。
thickness啊

gen surface有什么用？
生成表面
显示模型中用到

set dyn=on, grav 0 0 -10, hist_rep=1

是在动力分析中的问题，请问最后一个变量是什么含义，
好像只能是整数，查了Fish说明也没有看到解释，
请赐教！
hist_rep n
Histories are sampled every n timesteps. The default is n = 10. The
synonym HISTORY nstep may also be used. 每10个步时记录一次

请教:
(1) 谁用过 GENERATE surface xarc命令生成过圆弧面，能否贴出来让小弟学一下。我曾按照manual作过一些尝试但没有成功。

(2) 根据说明我曾生成过圆形隧道模型，但如何在隧道内壁上加正应力，请各位高手指点。
加正应力应该是用nstress命令
我没有做过曲面的例子，只用过规则的圆柱面

请问各位大侠， grout stiffness 是什么意思？
1.gr_k grout stiffness per unit length, kg [F/L2]
2.Ft = kg* ut
3.The grout annulus is assumed to behave as an elastic-perfectly plastic solid. As a result of relative shear displacement,
ut , between the tendon surface and the borehole surface, the shear force, Ft , mobilized
per length of cable is related to the grout stiffness, kg — i.e.,
4.可见k_g是grout annulus与zone（soil）之间的剪切刚度。（单位剪切位移引起的应力）
5.详细参考 cable 一节，有图形表示，很清楚。

写了几个命令流文件，但只要所分析问题的单元数一多，计算就很难收敛，总是在e-3或e-4之间跳来跳去，都迭代几万步了
请教大侠，怎么加快收敛速度？
一般一万五千单元左右要算多久？
flac 中的求解控制通常有
1）set mech force 50
solve
2）solve ratio 1e-7
3）step 30000
等几种，不知这几种怎么选择？每种有什么特点？？

加快收敛速度，要让单元的尺寸相差不大，另外，ratio恒定，说明产生了塑性流动或者破坏了，看看约束条件对不对！！材料模型对不对，材料参数是否太离谱？？

flac 中的求解控制通常有
1）set mech force 50
solve
2）solve ratio 1e-7
3）step 30000
第一种是控制不平衡力，是给定具体的值，必须迭代到这个值才结束运算。这是绝对的值，个人认为只有针对具体问题才使用；
第二种要重先设定，程序默认1e－5，－7也太小了，迭代的步数更多，更精确，一般情况下，个人认为－5就够了；
第三种是给定迭代步数，不至于30000吧，一般3000吧。但这种设定不一定能达到1e－5，所以要在有把握的情况下使用。
至于在循环中，solve与step n到底有何差别，本人也不甚清楚。有时某一步可能需要7000步才能达到－5次方，而其他的循环只需要1000步就达到了，这种情况下究竟该如何处理？请高手指点！
solve是通过控制不平衡力或比例的大小决定计算步
step ，cyc等命令则是由作者控制计算迭代步，在无法用step估计计算是否收敛情况下，多用solve，step自身的优势可以应用于逐步加载等方面。
加大bulk，shear的量纲试一试，还有就是检查一下边界条件！

;产生一个体，这个体里面包含隧道
gen zone radbrick &
p0 (0,0,0) p1 (10,0,0) p2 (0,10,0) p3 (0,0,10) &
size 3,5,5,7 &
ratio 1,1,1,1.5 &
dim 1 4 2 fill
;设置模型服从摩尔库仑屈服准则
mod moh
;定义bulk模量，剪切模量，摩擦角，粘聚力，和tens
prop bulk 1e8 shear 3e8 fric 35 coh 1e3 tens 1e3
;施加约束
fix x range x -0.1 0.1
fix z range z -0.1 0.1
fix y range y 9.9 10.1
;施加重力
set grav 0 10 0
;设置密度
ini dens 1000
;设置初始应力，这个有什么用呢？？？？
ini syy 0.0 grad 0 -10000 0
ini sxx 0.0 grad 0 -5000 0
ini szz 0.0 grad 0 -5000 0
apply sxx 0.0 grad 0 -5000 0 range x 9.9 10.1
apply szz 0.0 grad 0 -5000 0 range z 9.9 10.1
;挖坑的命令，使用fish语言
def digHole
command
range name trench x 0 1 y 0 4 z 0 2
;挖了
model nul range trench
set large
hist gp xdisp 1 0 0
endcommand
end
;表面蓝色显示
plot surf blue
;停顿
pause
;开挖
digHole
;设置大变形
set large
;记录位移变化
hist gp xdisp 1,0,0
; ... and then calculate the response
;设置计算步为1000
step 1000
;我不知道上述解释是否明白，是否对头？？

cjiao2000 wrote:
;产生一个体，这个体里面包含隧道
gen zone radbrick &
p0 (0,0,0) p1 (10,0,0) p2 (0,10,0) p3 (0,0,10) &
size 3,5,5,7 &
ratio 1,1,1,1.5 &
dim 1 4 2 fill
;设置模型服从摩尔库仑屈服准则
mod moh
;定义bulk模量，剪切模量，摩擦角，粘聚力，和tens
prop bulk 1e8 shear 3e8 fric 35 coh 1e3 tens 1e3
;施加约束
fix x range x -0.1 0.1
fix z range z -0.1 0.1
fix y range y 9.9 10.1
;施加重力
set grav 0 10 0
;设置密度
ini dens 1000
;设置初始应力，这个有什么用呢？？？？
ini syy 0.0 grad 0 -10000 0
ini sxx 0.0 grad 0 -5000 0
ini szz 0.0 grad 0 -5000 0
别的软件很少可以加地应力的，怎么说没用呢？这可是Flac的专长呀，象ansys也没的办法的。
apply sxx 0.0 grad 0 -5000 0 range x 9.9 10.1
apply szz 0.0 grad 0 -5000 0 range z 9.9 10.1
;挖坑的命令，使用fish语言
def digHole
command
range name trench x 0 1 y 0 4 z 0 2
;挖了
model nul range trench
set large
hist gp xdisp 1 0 0
endcommand
end
;表面蓝色显示
plot surf blue
;休息
pause
;开挖
digHole
;设置大变形
set large
;这个好像是施加位移
我认为是记录坐标是（1，0，0，）点的位移历史记录吧，不是施加位移吧。权当抛砖引玉了，各位大侠不要见笑
hist gp xdisp 1,0,0
; ... and then calculate the response
;设置计算步为1000
step 1000
;我不知道上述 解释是否明白，是否对头？？
首先确定位移是z向，还是x向

hist gp xdisp 1,0,0
有时候会出现位移都为0的情况,怎么回事?明明那是有位移的啊?
是啊，还有向上得哦！！！
;设置初始应力，这个有什么用呢？？？？
ini syy 0.0 grad 0 -10000 0
ini sxx 0.0 grad 0 -5000 0
ini szz 0.0 grad 0 -5000 0
别的软件很少可以加地应力的，怎么说没用呢？这可是Flac的专长呀，象ansys也没的办法的。
设置初始应力做什么哦？？？
一般情况下围岩地应力包括自重应力、构造应力、采掘应力。你所描速的情况应是在挖掘深度不大的情况下，只考虑土体的自重应力对开挖的影响。在古典压力理论中，认为地层的水平压力的侧压系数为u/1-u，这里取0.5。apply描述的模型的边界条件。

hist gp xdisp 1,0,0
是用于记录点（1，0，0）在开挖过程中的位移情况，模拟的目的就是要研究目标位置的位移或者应力变化情况。
set large
设置大变形的话,应该是按大变形理论吧!

最近做了算例，发现将cable单元的pretension设置为0或者100没有什么区别，大家看一下是怎么回事啊？

;generate the modle
title
'Vertical Anchored Gravity Wall'
;
set mecha ratio 1e-5
;set large

;Zone

gen zone radtunnel p0 1.5,0,-20 p1 1.5,20,-20 p2 0,0,-20 p3 1.5,0,0 &
p4 0,20,-20 p5 0,0,0 p6 1.5,20,0 p7 0,20,0 &
p8 1.5 6 -20 p9 1.5 0 -10.5 p10 0 6 -20 p11 0 0 -10.5 p12 1.5 6 -10.5 p13 0 6 -10.5 &
size 12,3,19,20 ratio 1,1,1,1 dim 9.5 5 9.5 5 fill

;Material
;For balance
model mohr
prop density = 1.80e3 bulk = 12.5e6 shear = 6.5e6 coh 10.0e3 tens 50e10 &
friction = 10.0 dilation = 0.0 range z -20.0 -18.0
prop density = 2.00e3 bulk = 20.0e6 shear = 12.0e6 coh 15.0e3 tens 50e10 &
friction = 25.0 dilation = 5.0 range z -18.0 -10.0
prop density = 2.00e3 bulk = 44.4e6 shear = 33.3e6 coh 0.0e3 tens 50e10 &
friction = 30.0 dilation = 10.0 range z -10.0 0.0

;Boundary conditon
fix x range x -0.1 0.1
fix x range x 0.4 1.5
fix y range y -0.1 .1
fix y range y 19.9 20.1
fix z range z -0.1 0.1

;Set gravitional acc
set grav 0,0,10

;Initial stress
ini szz -4.0e5 grad 0 0 -2.0e4
ini sxx -2.0e5 grad 0 0 -1.00e4
ini syy -2.0e5 grad 0 0 -1.00e4

;Initial State,Solve for balance
solve
ini state 0
ini xdisp 0.0 ydisp 0.0 zdisp 0.0

;Reinstall the soil
model mohr
prop density = 1.80e3 bulk = 12.5e6 shear = 6.5e6 coh 10.0e3 tens 58e3 &
friction = 10.0 dilation = 0.0 range z -20.0 -18.0
prop density = 2.00e3 bulk = 20.0e6 shear = 12.0e6 coh 15.0e3 tens 32.6e3 &
friction = 25.0 dilation = 5.0 range z -18.0 -10.0
prop density = 2.00e3 bulk = 44.4e6 shear = 33.3e6 coh 0.0e3 tens 5e3 &
friction = 30.0 dilation = 10.0 range z -10.0 0.0

;Set the surcharge
apply szz 10e3 range y 6 20 z -20.1 -19.9
solve

;Begin to excavate
;Step 1: 1.5m off
mod null range y 0 6.0 z -18.5 -20.0
sel liner range y 5.9 6.1 z -18.5 -20.0
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete
sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=4e6 cs_scoh=4e4 cs_scohres=2e6 cs_sfric=20.0
solve
save step1.sav

;Step 2: set the soilnail,and 1.5m off
sel cable id=1 begin=(0.75, 6.0, -18.5) end=(0.75, 15.5, -17.5) nseg=12 pretention 50e3
sel cable prop xcarea=5.1e-3 emod=200e9 ytens=210e6 &
gr_k=6e6 gr_coh=40e3
sel cable prop gr_per=0.314 gr_fric=25.0
mod null range y 0 6.0 z -17.0 -18.5
sel liner range y 5.9 6.1 z -17.0 -18.5
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete
sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=4e6 cs_scoh=4e4 cs_scohres=2e6 cs_sfric=20.0
solve
save step2.sav

;Step 3: set the soilnail,and 1.5m off
sel cable id=2 begin=(0.75, 6.0, -17.0) end=(0.75, 15.5, -16.0) nseg=12 pretention 50e3
sel cable prop xcarea=5.1e-3 emod=200e9 ytens=210e6 &
gr_k=6e6 gr_coh=40e3
sel cable prop gr_per=0.314 gr_fric=25.0
mod null range y 0 6.0 z -15.5 -17.0
sel liner range y 5.9 6.1 z -15.5 -17.0
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete
sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=4e6 cs_scoh=4e4 cs_scohres=2e6 cs_sfric=20.0
solve
save step3.sav

;Step 4: set the soilnail,and 1.5m off
sel cable id=3 begin=(0.75, 6.0, -15.5) end=(0.75, 15.5, -14.5) nseg=12 pretention 50e3
sel cable prop xcarea=5.1e-3 emod=200e9 ytens=210e6 &
gr_k=6e6 gr_coh=40e3
sel cable prop gr_per=0.314 gr_fric=25.0
mod null range y 0 6.0 z -14.0 -15.5
sel liner range y 5.9 6.1 z -14.0 -15.5
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete
sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=4e6 cs_scoh=4e4 cs_scohres=2e6 cs_sfric=20.0
solve
save step4.sav

;Step 5: set the soilnail,and 1.5m off
sel cable id=3 begin=(0.75, 6.0, -14.0) end=(0.75, 15.5, -13.0) nseg=12 pretention 50e3
sel cable prop xcarea=5.1e-3 emod=200e9 ytens=210e6 &
gr_k=6e6 gr_coh=40e3
sel cable prop gr_per=0.314 gr_fric=25.0
mod null range y 0 6.0 z -12.5 -14.0
sel liner range y 5.9 6.1 z -12.5 -14.0
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete
sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=4e6 cs_scoh=4e4 cs_scohres=2e6 cs_sfric=20.0
solve
save step5.sav

;Step 6: set the soilnail,and 1.5m off
sel cable id=3 begin=(0.75, 6.0, -12.5) end=(0.75, 15.5, -11.5) nseg=12 pretention 50e3
sel cable prop xcarea=5.1e-3 emod=200e9 ytens=210e6 &
gr_k=6e6 gr_coh=40e3
sel cable prop gr_per=0.314 gr_fric=25.0
mod null range y 0 6.0 z -11.5 -12.5
sel liner range y 5.9 6.1 z -11.5 -12.5
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete
sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=4e6 cs_scoh=4e4 cs_scohres=2e6 cs_sfric=20.0
solve
save step6.sav

;Step 7: set the soilnail,and 1.5m off
sel cable id=3 begin=(0.75, 6.0, -11.0) end=(0.75, 15.5, -10.0) nseg=12 pretention 50e3
sel cable prop xcarea=5.1e-3 emod=200e9 ytens=210e6 &
gr_k=6e6 gr_coh=40e3
sel cable prop gr_per=0.314 gr_fric=25.0
mod null range y 0 6.0 z -10.5 -11.5
sel liner range y 5.9 6.1 z -10.5 -11.5
sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete
sel liner PROP cs_nk=8e8 cs_sk=8e8 &
cs_ncut=4e6 cs_scoh=4e4 cs_scohres=2e6 cs_sfric=20.0
solve
save step7.sav

flac怎么作一个剖面！就是我建了一个模型，怎么看一个切面的视图。

pl set plane ori 0 0 0 nor 1 0 0
pl con disp plane

在所建立的模型中有两个命令；就是print gp position 与print zone gp。
不知道这两个命令中的gp是否有联系？
按找手册中的解释，print gp position 是用来输出gp的位置坐标的；
print zone gp是用来输出与zone相关的gp的id值的。
这样我就有一个上面的疑问，两个命令中的gp都有id值，是不是不同命令中的gp的id值相同的话，他们的坐标就可以看作是一样的啊？

是
一个gp有且只有一个id

谢谢dynamax兄。
这下我就明白了，
如果要画等值线的话，先输出gp的坐标来，将其另存为一个文件，再输出zone的应力值，将与zone对应的gp的id与前面输出的gp的id值对应起来，找到gp的坐标。用excel打开上面的文件进行编辑，最后用surfer画等值线。
这样的办法太笨，不知到有没有更好的办法啊？
欢迎前来讨论。

你的思路很对！
要简化操作，
如果你对fish熟悉，可以直接按surfer的文件格式输出
如果会用tecplot，可以参考前面我那篇
【原创】用tecplot对flac3d计算结果进行后处理
surfer我不熟，原则应该是一样的
看了你的哪个原创，发现里面有这样几个语句不太理解，
IO_READ = 0
IO_WRITE = 1
IO_FISH = 0
IO_ASCII = 1
其中的0，1具体指的是什么啊？

我发现一般zone中显示应力的gp的id往往不是全部模型的gp，那么如果想要用fish语句输出这些gp的坐标 ，该如和判断呢？
不知道在fish语句中如何将zone中的gp的坐标给出。
可以将gp的坐标和应力值存到一个数据文件，然后用surfer调用数据文件即可，不须在做额外的编辑

比如这个命令： print zone s输出了所有区的应力，可是要如果把它保存到指定写字板呢

set log on
set logfile *.log
print zone s
依次执行以上命令便可以实现你的目的。

其实简单的理解就是：
dd：映射平面法向与y轴正向夹角（顺时针）
dip：映射平面法向与z轴负向夹角
dip 是倾角，就是plane和水平面的夹角
dd 是dip direction ，倾角的方向，注意不是走向
就是你用来看倾角的那根线，在水平面上的投影，与y轴正向的交角
从y轴正向开始算
那位知道flac3D与ansys中单元不等比例剖分中网格长度与比例系数之间的计算公式。

s=[a(1-q^n-1)]/1-q
a 第一个网格长度
q 比例系数
不知版主问的是不是这个公式
等比公式
在做边坡稳定性分析建模的时候,一定会遇到非自立边坡,此时生成模型的时候有几种方法.大家用什么方法?
1.用不同的zone brick 拼装起来;
2.用plane above (该命令"切"出来的斜坡呈锯齿状,很难看,除非网格很细,否则与实际偏离较大)
3.使用fish,手册上有一个例子通过改变各个grid 的坐标来实现隧道的球形开挖面,依此类推,一定可以修改grid的坐标实现斜坡.

用fish或用zone brick 拼装生成非90度的基坑边坡，能给个例子吗？

下面就是我编的，请大家帮我看看，到底哪里有问题，怎么老是有问题。（我观察到第一步开挖后，图上显示只挖掉了一排网格，我要挖的可是两排。在第二次开挖前，shellSEL竟然只布置上去2个！）
我改了好多次都是类似结果，不知到底我哪弄错了!

gen zone radbrick p1 60 0 0 p2 0 60 0 p3 0 0 50 size 20 20 10 32 &
dim 20 20 10 rat 1 1 1 1.05 fill

model mohr
prop bulk 1.25e7 shear 5.77e6 fric 20 coh 1.5e4 range z 0.0 50.0

;边界条件
fix x range x -0.01 0.01
fix y range y -0.01 0.01
fix x y z range z 49.99 50.01
fix x y range x 59.99 60.01
fix x y range y 59.99 60.01
;初始条件
ini density 1.9e3 range z 0.00 50.00
set gravity 0 0 10
ini szz 0 grad 0 0 1.9e4 range z 0.00 50.00
ini sxx 0 grad 0 0 7.6e3 range z 0.00 50.00
ini syy 0 grad 0 0 7.6e3 range z 0.00 50.00

set large
hist id 1 unbal
plot hist 1
set mech ratio=1e-4
solve
save j_ini.sav
ini xdis=0 ydis=0 zdis=0

group lod1 range x 20 30 y 0 30 z -0.01 0.01
group lod2 range x 0 30 y 20 30 z -0.01 0.01
group lod range group lod1 any group lod2 any
apply szz 10e3 range group lod
model null range x 0 20 y 0 20 z 0.0 2.0
step 2000
save j_eca1.sav

sel shell id=1 range x 19.9 20.1 y 0 20 z 0.0 2.0
sel shell id=2 range x 0 20 y 19.9 20.1 z 0.0 2.0
sel node fix x yr zr range x -0.01 0.01
sel node fix y xr zr range y -0.01 0.01
sel shell id=1 prop iso=(26.6e9,0.25) thick=0.1
sel shell id=2 prop iso=(26.6e9,0.25) thick=0.1

model null range x 0 20 y 0 20 z 2.0 4.0
step 2000
save j_eca2.sav

sel shell id=3 range x 19.9 20.1 y 0 20 z 2.0 4.0
sel shell id=4 range x 0 20 y 19.9 20.1 z 2.0 4.0
sel node fix x yr zr range x -0.01 0.01
sel node fix y xr zr range y -0.01 0.01
sel shell id=3 prop iso=(26.6e9,