6.属性赋值(Properties)
我们首先定义一组具有密度(2700 kg/m3),杨氏模量(50 GPa)和泊松比(0.2)的块体属性。 以下命令中的值以MPa为单位给出。
; --- MAT 1 : rock ---
; density = 2700 kg/m3 = 0.0027e6 kg/m3
; E=50 GPa = 50e3 MPa , Poisson's ratio=0.2
prop mat 1 dens 0.0027 ymod 50e3 pratio 0.2
然后,我们定义法向刚度(10 GPa),剪切刚度(2 GPa)和摩擦角(25°)的节理材料。
; --- JMAT=1 : rock joints ---
prop jmat 1 jkn 10000 jks 2000 jfric 25
使用 change
命令将这些属性分别分配给块体和节理:
; --- assign material numbers ---
change mat 1
change jmat 1
7.初始应力(In-situ stress)
我们假设隧道处于200 m的深度,因此在z = 0时,z向下方向的原位应力为5.4 MPa。 假设水平应力为垂直应力的一半,则x和y方向的原位应力为2.7 MPa(压缩时)。 坐标系如下,位于隧道中部:
使用 insitu stress
设置初始应力,会有一个随深度变化的应力梯度,该梯度使用zgrad
关键字指定。 梯度是在垂直方向上每米应力的变化。 应力随坐标的上升而减小,并且3DEC中的应力压为负。
xx 应力和 yy 应力将从地面 z = 200m 的 0 MPa 变化到 z = 0 m(位于地面以下200m)的 2.7 MPa,而 zz 应力将从地面z = 200m的 0 MPa变化至 z = 0 m(位于地面以下200m)的5.4 MPa。具体计算可参考下面的命令注释:
; --- insitu stress state ---
; assume tunnel at 200 m depth
; vertical stress: szz=ρg(z-200) = (0.0027*g)*(z-200)
; at z=200: szz=0 ; z=0: szz=-5.4
; z-gradient of szz: 0.027
; (positive: less compression going up)
; horizontal sxx=syy=0.5*szz
;
insitu stress -2.7 -2.7 -5.4 0 0 0 &
zgrad 0.0135 0.0135 0.027 0 0 0
使用以下命令将重力应用于模型:
; gravity
grav 0 0 -10
与原位应力相关的是边界应力条件。 为了在深度上模拟块体模型,将5.13 MPa的向下垂直应力施加到模型顶部(z = 10)。 施加应力时,通过指定速度为0将块体的侧面和底部固定在适当的位置。
; --- boundary conditions for insitu stress state ---
; top of model (z=10): szz=-0.027*190=-5.13
bound stress 0 0 -5.13 0 0 0 range z 10.0
; bottom
bound zvel 0 range z -10.0
; sides
bound xvel 0 range x -10.0
bound xvel 0 range x 10.0
bound yvel 0 range y -10.0
bound yvel 0 range y 10.0
8.记录监测(Monitoring)
使用以下命令记录检测不平衡力以及x,y和z方向的位移。 默认情况下,每计算10步记录一次。
对于此模型,使用 hist ncycle
关键字将每次记录的循环数设置为 1。 历史记录会按照定义顺序自动分配 ID值(即不平衡力= 1,x位移= 2,y位移= 3,z位移= 4)。
; --- histories to monitor convergence ---
hist ncycle 1 unbal
; top of model
hist xdis 0 0 10 ydis 0 0 10 zdis 0 0 10
9.初始平衡(Initial state)
使用 solve
命令运行模型,直到模型达到平衡状态。 平衡状态定义为模型中平均不平衡力与施加力之比小于1e-5 的状态(此比例可以使用Solve ratio
命令更改,默认比例为1e-5)。 达到稳定状态后,将其另存为 “tun_c.3dsav”。
solve
save tun_c
执行以上所有的命令。 要查看位移的历史记录,请通过 File → New Item → Plot 创建一个新绘图窗口。 将会出现提示,要求输入名称。 输入“Hist”,然后单击 “OK”。 在控制面板 List 选项卡下,找到 Charts 选项卡,然后双击 History ,将出现一个空白图表。
通过单击下图中的加号可以将记录的变量添加到此绘图中。 默认情况下,横坐标为循环计算的次数,当然这可以自己更改的。
也可以将位移记录绘制出来:
10.无支护开挖(Excavate Unsupported)
第一种情况是模拟无支护开挖隧道。 创建一个名为 “unsupported.3ddat” 的新数据文件(快捷键Ctrl + N)。 使用 restore
命令恢复 tun_c.3dsav 保存的平衡状态。
使用 reset
命令将记录的位移,时间和所有历史记录重置为零。 将历史记录重新设置为每10步记录一次(这将加快计算速度)。
;--- File to simulate usupported excavation ---
restore tun_c.3dsav
reset disp time hist
his ncycle 10
首先使用 hide
命令隐藏所有块体,然后在x,y和z方向上寻找一个较小的范围,可以识别由步骤2中指定的三个节理切割形成的楔形体。 然后,我们将定位的楔形块体分配给名为 “wedge” 的组。
hide
seek range xr -0.3 .46 zr 3.8 4.5 yr -1.0 1.0
group block "wedge"
为了记录楔形体的垂直位移,我们将z位移历史点指定为靠近楔形的顶点之一。 在这种情况下,我们使用点(0,0,3.9)。 通过指定接近所需顶点的点坐标,3DEC将自动找到最近的顶点。 由于我们重置了历史记录,因此第一个定义的历史记录现在是z位移,可以使用 hist label
命令将其标记为 “Vertical Displacement”。使用 seek
命令使隧道内部可见。
hist zdis 0 0 3.9
hist label 1 'Vertical Displacement'
seek
现在可以执行命令。 使用在前面FISH函数中定义的y坐标,可以从隧道的起点开始逐步移动开挖隧道,并在四个阶段中沿y轴逐渐移动。 回想一下,区域5是隧道。 使用命令指定两个范围将给出两个范围的交集。
下面的命令确保仅挖掘隧道的各个部分。在每个阶段之间,执行数千次循环以模拟与挖掘相关的时间,并允许模型尝试达到平衡。 在最后阶段之后,将执行10,000个循环,以显示历史记录所示的楔形的向下位移。 请注意,由于楔块会掉落,因此直到楔块撞到隧道地面后,模型才会慢慢求解平衡。 这将花费很长时间,因此我们使用 cycle
命令指定了多个循环计算次数,而不是使用 solve
。
; delete interior blocks in stages
delete range group tunnel y @stage0,@stage1
cycle 1000
delete range group tunnel y @stage1,@stage2
cycle 2000
delete range group tunnel y @stage2,@stage3
cycle 4000
delete range group tunnel y @stage3,@stage4
cycle 10000
save tun_x.3dsav
运行模型后,通过上面绘制历史记录的方法,可以绘制出垂直位移图:
注意:这里使用了命令delete
。 如果您打算用材料来填充隧道(例如回填),则应改用excavate
命令。 如果希望能够绘制出开挖的材料,但又不想将其重新引入模型中,请使用remove
命令。
可以看到在前两个开挖阶段几乎没有位移,但是在第三部分开挖之后,楔块掉了下来。 这种楔形破坏也可以在模块模型中看到。 切换到 “Blocks” 视图窗口。 为了最好地可视化块体破坏,我们可以如下图所示将块模型的一部分隐藏在一个节理上方。 在控制面板的 Plot Items 窗格中,单击 BGroup 项目旁边的范围图标。 在下面 Elements 中双击 Plane。 点击 Plane 在 Attributes 选项卡中,输入0 0 5.7作为原点,输入60作为倾角,输入130作为倾角方向。 将 Kind 设置为 Below。
设置完成后,点击OK,可以得到切片视图。通过单击 Joints 图项旁边的眼睛图标隐藏节理。 通过选择 List 选项卡并选择 Blocks → Vectors → Displacement 来添加位移向量进行显示。 单击 Arrow Head 复选框以显示矢量方向。 也可以像上面在块体上的操作一样将矢量隐藏在平面上方,但是在这种情况下,这实际上不是必需的。 现在该图应如下所示:
Comments | 14 条评论
请教个问题,有个地方没有理解透。
在边界荷载那个地方,前面已经有insitu stress,而且可以计算模型上方的垂向应力,为什么还要使用boundary stress 重复施加,怎么理解,谢谢。
insitu stress -2.7e6 -2.7e6 -5.4e6 0 0 0 &
zgrad 1.35e4 1.35e4 2.7e4 0 0 0
;
;Gravity
gravity 0 0 -10.0
;
;—-boundary conditions for insitu stress state—
;top of model (z=10): szz=-2700g190=-5.13 MPa
bound stress 0 0 -5.13e6 0 0 0 range z 10.0
@xuxueliang2021 insitu stress是生成整个场的初始应力,模型是20×20×20这么大,但是要模拟200米的深度,需要构造一个应力场,使得整个模型区域都处于这个应力场中。而求解一个问题必须要设置边界条件,所以这里在模型顶部施加一个压力来设置上覆岩体的重力,然后才能进行求解
这个地方不好理解。
从insitu stress的应力场赋值过程可以看见,它是基于坐标系的原点计算的。
建模时已经设定了模型相对原点的位置,也就是能直接能推算出模型顶部所受的垂向力。
可原模型代码又重新在顶部施加一个垂向压力,是不是重复了?
@xuxueliang2021 insitu stress可以看成初始条件,boun stress是边界条件,两个并不重复。边界条件是计算问题必需的,而初始条件是给模型一个初始状态,大概可以这样理解吧
嗯嗯,这样理解也不错。感谢。
只是有的实例代码,里面只有insitu stress,没有boun stress,苦恼过。
关于地应力有一些疑问还请博主解答,
施加地应力命令insitu sxxo syyo szzo sxyo sxzo syzo 这里sxxo、syyo两个是不是可以理解成共同形成水平应力呢,szzo理解成竖直应力,sxyo、sxzo、syzo应该怎么和前面三个应力对照理解呢
目前在研究一个400400400m的模型,向上一直建到地表,资料中提到埋深较大天然地应力场以进水平构造应力为主,并且给了研究区域的最大主应力值和方位角,是不是就可以通过sxxo和syyo来形成研究区域的应力场和方位角呢,然后其他不重点研究的区域的水平应力是不是可以和研究区域的水平应力设置成一样呢,模型的竖直应力是不是可以忽略呢,就简单的考虑重力。
模型的初始应力施加后,是否还要在模型四周和边界施加应力条件呢
问题比较多,还望博主解答
补充一下上一条

sxyo、sxzo、syzo是不是可以理解成xy、xz、yz这每两个方向合成的地应力
地表的应力一般是不是0呢
@qifa insitu stress sxxo syyo szzo sxyo sxzo syzo xgradient … 用来初始化应力,以应力梯度的方式施加。这里带o的应力(sxxo…)是你坐标原点处的应力,跟你的坐标系有关,原点在地表就是0,在地表以下的话就要算一下原点的地应力值。根据剪应力互等,一个单元体有6个应力分量,就是这里的6个值,可参考弹性力学的单元体的应力分量,如果你只考虑水平和竖向没有剪应力,sxyo等即可设为0,只加3个主要方向的。随深度变化的应力可以按下面计算,具体可看帮助文档的介绍:
sxx = sxxo + (sxxx)x + (sxxy)y + (sxxz)z
@qifa syy = syyo + (syyx)x + (syyy)y + (syyz)z
@qifa szz = szzo + (szzx)x + (szzy)y + (szzz)z
@qifa sxy = sxyo + (sxyx)x + (sxyy)y + (sxyz)z
@qifa sxz = sxzo + (sxzx)x + (sxzy)y + (sxzz)z
@qifa syz = syzo + (syzx)x + (syzy)y + (syzz)z
@qifa 如果是以主应力施加,可以用insitu principal …指定主应力大小和方位
竖向应力一般不可忽略的,加了重力加速度就相当于加了竖向应力,这个还要注意前面原点应力值的计算,确认施加的应力梯度是否和重力加速度算出的应力一致。
边界条件肯定是要施加的,不然无法计算,静力问题就加零速度边界固定就行了。
solve 命令运行时显示未设置模数,无法循环,是哪里出现问题了啊
@927585067 你是按照例子的命令运行出问题了吗?
我是自己设置了好多岩层,然后开挖,已经解决啦,感谢
模拟采场开挖,但是zz stress图是空白是怎么回事?
博主您好,请问怎样在一条垂直向下的隧道的内壁上添加初始应力呢?
我参考手册试写了这样的命令
insitu stress 3e6 range cylinder end1 0 0 -3 end2 0 0 3 radius 0.15
运行没有报错,但是我不知道怎么在plot里显示这个应力添加的结果,所以不知道这样用是否正确,也不知道手册里对range cylinder命令所解释的“圆柱体的一端(end1)位于(x1, y1, z1),另一端(end2)位于(x2, y2, z2),半径为r。如果指定了r2,则范围为从r到r2的圆柱体环。”意思是选取整个圆柱体还是圆柱体的内侧面,非常期待博主的解答
(=・ω・=)
@248238376 添加的初始应力我不知道能不能在边界里显示,我原来写过检查边界条件的方法,你试试能不能显示:https://geomatlab.com/3dec-check-condition/. 至于range cylinder是选面还是体,我还不清楚