[工程实例]岩土力学求解围岩特征曲线_隧道开挖

nevercry_zju

申请加精
这个例子很经典
其中用到宏rf enforce 等  见附件


! 直接法
! 岩土力学模式
! 求解围岩特征曲线
finish
/clear
toolbar
/title,Excavation of tunnel
r=5
r1=2*r
a=10*r
b=6*r
c=3*r
/PREP7
/PNUM,KP,1  
! 生成关键点
k,1,0,0
k,2,r,0
k,3,0,r
k,4,r1,0
k,5,0,r1
k,6,c,0
k,7,0,c
k,8,c,c
k,9,a,0
k,10,0,b
k,11,c,b
k,12,a,b
k,13,a,c
! 生成线
/PNUM,KP,0
/PNUM,LINE,1
l,1,2
larc,2,3,1,r
l,3,1
l,2,4
larc,4,5,1,r1
l,5,3
l,4,6
l,6,8
l,7,8
l,5,7
l,7,10
l,10,11
l,8,11
l,11,12
l,12,13
l,8,13
l,9,13
l,6,9
LPLOT  
! 生成面
/PNUM,LINE,0
/PNUM,AREA,1
al,1,2,3
al,2,4,5,6
al,5,7,8,9,10
al,9,11,12,13
al,13,14,15,16
al,8,16,17,18
! 镜像
asel,all
arsym,x,all
asel,all
arsym,y,all
APLOT
! 改变第二、四象限衬砌圆弧段的法线方向
csys,1
lsel,s,loc,y,90,180
lsel,a,loc,y,270,360
lsel,r,loc,x,r
LREVERSE,all,0  
csys,0
! 合并重合的关键点和线并压缩编号
NUMMRG,ALL
NUMCMP,ALL

! 定义衬砌单元
ET,1,BEAM3  
mp,ex,1,3e10             ! 设置弹性模量ex
mp,prxy,1,0.2            ! 设置泊松比
mp,dens,1,2500           ! 设置密度
d=0.3
r,1,d,d*d*d/12,d         ! 设置实常数
! 定义围岩单元
ET,2,PLANE42
KEYOPT,2,3,2             ! 平面应变
mp,ex,2,6e8              ! 设置弹性模量ex
mp,prxy,2,0.38           ! 设置泊松比
mp,dens,2,1900           ! 设置密度
TBDE,DP,2   
TB,DP,2,,,0
TBMODIF,1,1,80000   
TBMODIF,1,2,27
! 划分网格
type,1
real,1
mat,1
MSHAPE,0,2D
MSHKEY,1
! 划分衬砌单元
csys,1
lsel,s,loc,x,r
lesize,all,,,10
lmesh,all
csys,0
allsel,all
type,2
real,2
mat,2
! 划分内圆网格
csys,1
asel,s,loc,x,0,r
amesh,all
allsel,all
! 划分圆环网格
asel,s,loc,x,r,r1
amesh,all
csys,0
allsel,all
/PNUM,KP,1
/PNUM,LINE,1
/PNUM,AREA,1
! 划分五边形网格
AMAP,3,4,6,7,5
AMAP,9,5,7,17,15
AMAP,21,15,17,24,23
AMAP,15,23,24,6,4
! 划分左右两边的中间部分网格
lsel,s,loc,x,-a,a
lsel,u,loc,x,-c,c
lsel,u,loc,x,-a
lsel,u,loc,x,a
lsel,r,loc,y,-c,c
n1=nint(5*distkp(6,9)/distkp(6,8))
lesize,all,,,n1
asel,s,loc,x,-a,a
asel,u,loc,x,-c,c
asel,r,loc,y,-c,c
amesh,all
n1=
! 划分上下两边的中间部分网格
lsel,s,loc,y,-b,b
lsel,u,loc,y,-c,c
lsel,u,loc,y,-b
lsel,u,loc,y,b
lsel,r,loc,x,-c,c
n2=nint(5*distkp(7,10)/distkp(7,8))
lesize,all,,,n2
asel,s,loc,y,-b,b
asel,u,loc,y,-c,c
asel,r,loc,x,-c,c
amesh,all
n2=
asel,all
! 划分剩余部分网格
asel,u,loc,x,-c,c
asel,u,loc,y,-c,c
aplot
amesh,all
allsel,all
! 合并重合的关键点和线并压缩编号
NUMMRG,ALL
NUMCMP,ALL
finish

/solu
ALLSEL,ALL
ANTYPE,0
NLGEOM,1                      ! 打开大变形
TIME,1                        ! 为载荷步设置时间
NSUBST,10,0.2,0.01            ! 指定载荷步中所需要的子步数
NROPT,FULL, ,ON               ! 在静态或完全瞬态分析中，指定Newton-Raphson选项
!LUMPM,1                      ! 指定一个集中质量矩阵公式
! 施加约束
nsel,all
nsel,s,loc,x,-a
nsel,a,loc,x,a
d,all,ux,0
nsel,all
nsel,s,loc,y,-b
d,all,uy,0
! 加载
! 施加重力
ALLSEL,ALL
ACEL,0,10,0,
! 杀死梁单元
ESEL,S,TYPE,,1  
ekill,all               
! 求解
allsel,all
solve
finish

/POST1
allsel,all
! 读出数据
SET,LAST
! 画变形云图
PLDISP,1
! 读取等效节点力
csys,1
APLOT   
csys,1  
ASEL,R,LOC,X,r,r1   
APLOT  
ESLA,S  
EPLOT
LSEL,R,LOC,X,r  
LPLOT  
NSLL,R,1
NPLOT   
rf
csys,0
finish

/SOLU   
ALLSEL,ALL
ANTYPE,,REST,1,6,0  
NLGEOM,1                      ! 打开大变形
TIME,2                        ! 为载荷步设置时间
NSUBST,10,0.2,0.01            ! 指定载荷步中所需要的子步数
NROPT,FULL, ,ON               ! 在静态或完全瞬态分析中，指定Newton-Raphson选项
!LUMPM,1                      ! 指定一个集中质量矩阵公式
! 杀死梁单元
ALLSEL,ALL
ESEL,S,TYPE,,1  
ekill,all     
! 杀死衬砌内的围岩单元
csys,1
ALLSEL,ALL
ASEL,S,loc,x,0,r  
ESLA,R     
ekill,all  
csys,0  
allsel,all
! 施加等效节点力
*ask,Factor,'Please input factor : factor=',0    ! 提示用户输入释放力系数,缺省值为0
enforce,Factor                  
! 求解
allsel,all
solve
finish

/POST1
allsel,all
! 读出数据
SET,LAST
! 画变形结果图
PLDISP,1
finish

/SOLU   
ALLSEL,ALL
ANTYPE,,REST,2,6,0  
NLGEOM,1                      ! 打开大变形
TIME,3                        ! 为载荷步设置时间
NSUBST,10,0.2,0.01            ! 指定载荷步中所需要的子步数
NROPT,FULL, ,ON               ! 在静态或完全瞬态分析中，指定Newton-Raphson选项
!LUMPM,1                      ! 指定一个集中质量矩阵公式
! 激活梁单元
ALLSEL,ALL
ESEL,S,TYPE,,1  
ealive,all
! 杀死衬砌内的围岩单元
csys,1
ALLSEL,ALL
ASEL,S,loc,x,0,r  
ESLA,R     
ekill,all  
csys,0  
allsel,all  
! 除去等效节点力
enforce,1                  
! 求解
allsel,all
solve
finish

/POST1
allsel,all
! 读出数据
SET,LAST
! 画变形结果图
PLDISP,1
! 建立单元表
ETABLE,NI,SMISC,1     !单元I点轴力
ETABLE,NJ,SMISC,7     !单元J点轴力
ETABLE,MI,SMISC,6     !单元I点弯矩
ETABLE,MJ,SMISC,12    !单元J点弯矩
! 更新单元表
ETABLE,REFL
! 画轴力分布图
/TITLE,Axial  force  diagram
PLLS,NI,NJ,0.1,0
! 画弯矩分布图
/TITLE,Bending  moment diagram
PLLS,MI,MJ,-0.1,0
finish

nevercry_zju

！这是enforce宏的命令流。


!  施加等效节点力
PARRES,NEW,'RF','pav'                            ! 恢复参数
factor=ARG1
factor=1-factor
*cfopen,data,mac                                 ! 打开新文件mydata.mac(如存在则覆盖原文件)
*dim,Q,array,Nnum,3,1                            ! 定义一个nnum行3列的数组
*do,i,1,Nnum
     Q(i,1)=P(i,1)
     Q(i,2)=-P(i,2)*factor                          
     Q(i,3)=-P(i,3)*factor                           
*enddo
*vwrite,Q(1,1),Q(1,2)                            ! 施加X方向等效节点力
('F,',f8.0,',FX,',f18.4)                         ! 写入格式
*vwrite,Q(1,1),Q(1,3)                            ! 施加Y方向等效节点力
('F,',f8.0,',FY,',f18.4)                         ! 写入格式                                       
*cfclos                                          ! 关闭文件
!施加等效节点力
nsel,all
data
*msg,ui                                          ! 输出成功消息
You have enforced nodal force successfully !                                 
i=                                               ! 删除参数                                       
Nnum=
factor=
Q(1,1)=
P(1,1)=

nevercry_zju

！这是rf宏的命令流。



! 读取等效节点力
*get,Nnum,node,0,count                     ! 读取所选取节点总数
*get,Nmin,node,0,num,min                   ! 读取所选取节点中最小节点号
*dim,P,array,Nnum,3,1                      ! 定义一个nnum行3列的数组
k=nmin                                     ! 令K等于所选取节点中最小节点号
*do,i,1,nnum                                                            
       P(k,1)=k                            ! 保存当前节点号   
       k=ndnext(k)                         ! 令K等于下一个所选取节点号
*enddo

*do,i,1,Nnum
   k=P(i,1)
   Nsel,s,node,,k
   Fsum,,,
   *GET,ForX,FSUM,0,ITEM,FX                ! 读取X方向节点力
   P(k,2)=forx                             ! 保存当前X方向节点力
   *GET,ForY,FSUM,0,ITEM,FY                ! 读取Y方向节点力
   P(k,3)=fory                             ! 保存当前Y方向节点力
*enddo
i=
k=
Nmin=
forx=
fory=
PARSAV,ALL,'RF','pav'                      ! 保存参数
*msg,ui                                    ! 输出成功消息
!You have read nodal force successfully !

[ 本帖最后由 nevercry_zju 于 2006-6-5 10:30 编辑 ]

llzzpp1980

我对你得例子很感兴趣，可惜rf附件要求权限太高了。

youngzz

哥们儿，能不能联系一下？我正在搞一个模拟隧道开挖的课题，有些地方不明白，想向你请教一下。youngzz@21cn.com

laomao_0217

wo yun

deeryao

非常感谢，先看看的说

deeryao

昨天做了一下，还不错，有几个问题想请教一下
1.  我看了地表沉降曲线，楼主的帖子中的沉降是绝对的，我认为是需要减去前一个载荷步才能得到每一步开挖得到的位移，这样的话，最大位移大约是两个毫米，但是边上的位移是0.8mm左右，要是我把材料模量降低之后，边上的位移大约大约是九个毫米左右，有点大。按说周边的位移应该是比较小才对，请教楼主，这个是由于什么模型尺寸的原因还是别的什么原因呢？
2.另外在运行的时候，rf宏命令好像直接输入不可以，要用图形操作才可以？请教楼主是这样子的么？
3.楼主在第二步中杀死了梁单元，并且杀死了衬砌内部的围岩单元，并且考虑应力释放，然后施加了节点力。但是在第三步中为什么还要杀死衬砌内部的围岩单元？而且最后不知道楼主是不是笔误？应该是施加等效节点力而不是除去等效节点力吧？
4.关于楼主取节点力的宏，我原来都是用图形操作取的，先选择所需节点的外部单元，然后选择单元的节点，最后在general postproc－－nodal calcs－－然后可以计算节点力，不知道用楼主的方法和我刚才说的是不是得到一样的结果？
希望楼主帮忙解释一下！！！谢谢！！！

xyghzzj

未经作者允许就把内部程序泄漏出去，小心师兄削你！


罪过！罪过！

[ 本帖最后由 xyghzzj 于 2006-10-2 21:57 编辑 ]

xyghzzj

好好看看下面的帖子（引用）：

三维锚喷隧道模拟开挖

看到很多朋友在问关于锚索的建模问题，闲着没事，做了一个，和大家讨论
1.单元类型：锚索link8，土体SOLID45，衬砌shell181
2.为简单起见，采用均质土，模型的尺寸随便定的
3.建模的基本步骤（其实很简单，不要笑话）:
   a.先建立二维结构；
   b.用vext拖拉成三维模型，拖拉过程应分步实施，因为是分步开挖，所以先对以后需开挖的部分拖拉，再EXTRUDE围岩；
   c.找出衬砌所在的面，划分衬砌结构
   d.在源面的锚索处划分锚索单元，此后用lgen复制整个支护结构，此过程也应定义成分组
4.计算

其实计算有很多种方法，简要的说一下，献丑了
第一种方法：不考虑应力释放系数
        在这种方法很简单，杀死土体与支护同时进行，一个循环下来就能搞定
第二种方法：考虑应力释放的过程
1.第一步开挖：提取开挖部分节点的等效节点力，杀死开挖部分土体，将等效节点力同时施加在对应的节点上
2.计算，因为要考虑应力释放的过程，比如释放系数为30%，那么需要将前面施加的等效节点力改为原始值的70%,此步不激活支护单元，计算后，激活支护单元，将施加的等效节点力改为0，计算。这样就完成了第一步的开挖。
3.重复1，2步骤即可完成整个开挖模拟

当然，考虑到避免支护体系的初始变形，对应的处理方法是使用初应力文件方法。但是实质的方法还是上面提到的。
      需要指出的是：对于围岩级别比较差的情况下，锚索的作用很弱，此时可用注浆层代替。
此处仅列出两种方法，其实还有其他的方法，只是这些比较常用。希望能和大家讨论。

xyghzzj

原帖由 deeryao 于 2006-8-24 11:41 发表
昨天做了一下，还不错，有几个问题想请教一下
1.  我看了地表沉降曲线，楼主的帖子中的沉降是绝对的，我认为是需要减去前一个载荷步才能得到每一步开挖得到的位移，这样的话，最大位移大约是两个毫米，但是边上的 ...

2.可以直接输入。
3.应力全部释放。
4.一样。

xyghzzj

①叶子的离开,是风的追求,还是树的不挽留? ②Rejoicing in hope , patient in tribulation ! nevercry_zju个人资料

主页:http://nevercry.simwe.com  
Email:zengex@gmail.com
QQ:10259728   
MSN:zengex@hotmail.com  
ICQ:  
Yahoo:  
淘宝旺旺:  
支付宝:  UID:164152  
注册日期:2006-04-29 21:11:43  
上次访问:2006-10-02 10:55:32  
最后发表:2006-09-30 10:44:26  
页面访问量:960  
在线时间:总计在线67.33小时, 本月在线0.67小时,升级剩余时间 13 小时
昵称:  
用户组:初级会员  
发帖数级别:小试牛刀  
阅读权限:20  
积分:8  
技术积分:8  
帖子:107(占全部帖子的 0.01%)
平均每日发帖:0.69 帖子
精华帖:0  
性别:保密
来自:  
生日:1984-05-21  
自我介绍:读研浙大岩土. 热心学习ANSYS,Plaxis,GeoSlope,FLAC,3d-Sigma等专业软件……  
访问数:3611  
文章总数:41  
日志数:12  
图片数:7  
文件数:4  
商品数:0  
书签数:18  
省份:北京  
城市:海淀  
访问域名:http://nevercry.simwe.com

deeryao

谢谢

inkey

该附件被下载次数 2, 阅读权限 20
天,太高了!

zhizhuozhe

师兄我现在也作隧洞，我把命令流传上来给我看看，到底稳地在哪里，还有算出来的结果我不知道对不对
ini
/cle
*set,b,10
*set,h1,8
*set,th,0.5
/prep7

k,,0,0 $k,,0,h1
k,,b/2,b/2+h1
k,,b,h1 $k,,b,0


l,1,2 $larc,2,3,5,b/2
larc,3,4,1,b/2
l,4,5
l,5,1


a,1,2,3,4,5
blc4,-3*b,-2*(h1+b/2),7*b,5*(h1+b/2)
aovl,1,2

numcmp,all


et,1,beam3
et,2,plane42
keyopt,2,3,2
r,1,th,th*th*th/12,th, ,

mp,ex,1,2.6e7
mp,prxy,1,0.2
mp,dens,1,25

mp,ex,2,1.3e5
mp,prxy,2,0.32
tb,dp,2
tbdata,1,200,20,
mp,dens,2,22

mp,ex,3,1.3e6
mp,prxy,3,0.32
tb,dp,3
tbdata,1,200,20,
mp,dens,3,22

lsel,s,line,,1
lesize,all,,,20,
latt,1,1,1
lmesh,all
alls

lsel,s,line,,2
lesize,all,,,20,
latt,1,1,1
lmesh,all
alls

lsel,s,line,,3
lesize,all,,,20,
latt,1,1,1
lmesh,all
alls

lsel,s,line,,4
lesize,all,,,20,
latt,1,1,1
lmesh,all
alls
lsel,s,line,,5
lesize,all,,,20,
latt,1,1,1
lmesh,all
alls
lsel,s,line,,6,9
lesize,all,,,20
alls

mopt,split,on
mopt,split,err
mopt,split,warn
mopt,qmesh,main

asel,s,area,,3
aatt,2,,2
amesh,all
alls

asel,s,area,,1
aatt,3,,2
amesh,all
alls

lsel,s,loc,x,-3*b
lsel,a,loc,x,-3*b+7*b
dl,all,,ux,0
alls
lsel,s,loc,y,-2*(h1+b/2)
dl,all,,uy,0
alls

acel,0,9.8,0
                              
fini


/solu
antype,static
deltim,0.1,0.05,0.2,off
autots,on !使用自动时间步
pred,on !打开时间步长预测器
lnsrch,on !打开线性搜索
nlgeom,on !打开大位移效果
nropt,full !设定牛顿-拉普森选项
cnvtol,m,,0.02,0

esel,s,type,,1
ekill,all
esel,all
esel,s,live
nsle,s
nsel,inve
d,all,all,0
nsel,all
esel,all
solve

esel,s,mat,,3
ekill,all

esel,s,type,,1
ealive,all
nsle,s
ddele,all,all
esel,all
esel,s,live
nsle,s
nsel,inve
d,all,all,0
nsel,all
esel,all

solve
fini
/post1
etable,if,smisc,1
etable,jf,smisc,7
etable,im,smisc,6
etable,jm,smisc,12