【原创】一个地震波探测溶洞的例子,请多指教!

sl3913

无量纲化处理了一下,边界给的是YEE吸收边界!
  
/BATCH,LIST
/VERIFY,EV129-1S
/CONFIG,NRES,5000  
*CFOPEN,disp,dat
*dim,vp,,3
*dim,vS,,3
! 表土层深度OH1LAM=0.4 $ 异常体埋深Odlam=2 $异常体直径 ORTlam=0.5 $ 异常体与表土层弹模比EX13=0.8 ! OXY=0.9
OH1LAM=0 $ Odlam=1 $ ORTlam=0.2 $ EX13=0.8 ! OXY=0.9
/PREP7
/TITLE,AMA,EV129-1S,FLUID129,TRANS ANALYSIS
ET,1,PLANE42                  ! structural element
ET,2,SURF153                    ! acoustic infinite line element
ET,3,PLANE42        ! structural element
! material properties
!*dim,UX1,array,1200   ! Create array parameter
!*dim,Uy1,array,1200   ! Create array parameter
dx=1
  
a=40
ENDTIME=0.2 $ DST=0.0005
  
!nt=1
  
absorblayer=2
  
lamda=10  !地震波基长
XDEPTHlam=1.5  !研究区深度与基长之比
XDEPTH=XDEPTHlam*lamda!研究区深度
YWIDTHlam=20  !研究区长度与基长之比
YWIDTH=YWIDTHlam*lamda!研究区长度
  
OYYCH=YWIDTH/2  !洞穴Y坐标
  
EX12=17926/224
DENS12=2600/1930
NUXY12=0.306/0.421
EX13=1.4/2.24
DENS13=1500/1930
NUXY13=0.421/0.421
mu13=10000/38000
DAMP12=0.8
DAMP13=1.5
  
EX1=2.24e8  !介质一为表层介质
DENS1=1930
NUXY1=0.421
!SONC1=860
mu1=38000
DAMP1=0.0001
  
EX2=EX1*EX12  !介质二为基岩
DENS2=DENS1*DENS12
NUXY2=NUXY1*NUXY12
DAMP2=DAMP12*DAMP1
!SONC2=1460
EX3=EX1*EX13  !介质三为岩溶洞穴
DENS3=DENS1*DENS13
NUXY3=NUXY1*NUXY13
mu3=MU1*MU13
DAMP3=DAMP13*DAMP1
  
EX4=ex1  !边界高吸收介质，对应于表层介质
DENS4=dens1
NUXY4=nuxy1
!SONC4=860
mu4=mu1
DAMP4=DAMP1
  
EX5=EX2  !边界高吸收介质，对应于基岩
DENS5=DENS2
NUXY5=NUXY2
DAMP5=DAMP2
  
!计算各种介质的P、S波速度
  
  lame1=Ex1*nuxy1/(1+nuxy1)/(1-2*nuxy1) $ lame2=EX1/2/(1+NUXY1) $ vp(1)=((lame1+lame2)/dens1)**.5 $ vs(1)=(lame2/dens1)**.5
  lame1=Ex2*nuxy2/(1+nuxy2)/(1-2*nuxy2) $ lame2=EX2/2/(1+NUXY2) $ vp(2)=((lame1+lame2)/dens2)**.5 $ vs(2)=(lame2/dens2)**.5
  lame1=Ex1*nuxy3/(1+nuxy3)/(1-2*nuxy3) $ lame2=EX3/2/(1+NUXY3) $ vp(3)=((lame1+lame2)/dens3)**.5 $ vs(3)=(lame2/dens3)**.5
  
*VSCFUN,vMAX,MAX,VP
  
DT=DX/VMAX/3
  
*IF,DT,GE,DST,THEN
DT=DST
*ENDIF
  
totalnt=ENDTIME/dt $ t0=1/a/2 $ f=600*a  
*IF,DST/DT,LT,NINT(DST/DT),THEN
RECORDSTEP=NINT(DST/DT)-1
*ELSE
  recordstep=NINT(DST/DT)
*ENDIF
nt=1/a/dt
  
MP,EX,1,EX1
MP,DENS,1,DENS1
MP,NUXY,1,NUXY1
MP,MU,1,mu1
mp,damp,1,damp1
MP,DENS,2,DENS2
MP,EX,2,EX2
MP,NUXY,2,NUXY2
MP,DENS,3,DENS3
!MP,EX,3,EX3
MP,NUXY,3,NUXY3
MP,MU,3,mu3  
MP,EX,4,EX4
MP,DENS,4,DENS4
MP,NUXY,4,NUXY4
MP,MU,4,mu4
mp,damp,4,damp4
MP,DENS,5,DENS5
MP,EX,5,EX5
MP,NUXY,5,NUXY5
mp,damp,5,damp5
! create square
RECTNG,0,XDEPTH,0,YWIDTH,  !生成研究区
ESIZE,dx,0,    !单元尺寸为2*2
MSHAPE,0,2D  
MSHKEY,1
amesh,all
  
!检波器位置
  
bounodemin=node(0,0,0) $ bounodemax=node(dx,YWIDTH,0)
  
!设定表土层和岩溶洞穴的位置与大小
!在此考虑洞穴埋深变化的影响
  
!OH1LAM=0.4 $ Odlam=2 $ ORTlam=0.5 $ EX13=0.8 ! OXY=0.9
  
esel,all
mat,1 $ real,1 $ type 1
emodify,all
  
OH1O=OH1lam*lamda  !表土层厚度
  
BURY=ODlam*lamda  !洞穴埋深  
  
ORTO=ORTlam*lamda  !洞穴直径
  
MP,EX,3,EX1*EX13!设定洞穴充填介质的弹模
  
NSEL,S,LOC,x,OH1O+dx,XDEPTH   !选择下伏岩层对应的节点
esln
TYPE,3   
MAT,2      !下伏岩层为介质2
emodif,all     !确认修改
  
ALLSEL,ALL  
NSEL,S,LOC,Y,YWIDTH   
NSEL,A,LOC,X,0      !选择边界节点
NSEL,A,LOC,X,XDEPTH
NSEL,A,LOC,Y,0
ESLN,S         !边界节点单元  
!TYPE,   2         !边界节点单元
!MAT,       2
!REAL,       3   
!ESURF         !边界节点单元定义为无限超声边界
ALLSEL,ALL  
!边界节点单元
NSEL,S,LOC,Y,YWIDTH-absorblayer,YWIDTH   
NSEL,A,LOC,Y,0,absorblayer   
NSEL,R,LOC,X,oh1o,XDEPTH  
nsel,a,loc,x,XDEPTH-absorblayer,XDEPTH
ESLN,S         !边界节点单元  
MAT,5
EMODIFY,ALL
*if,oh1o,gt,0,then
NSEL,S,LOC,Y,YWIDTH-absorblayer,YWIDTH  
NSEL,A,LOC,Y,0,absorblayer    
NSEL,r,LOC,X,0,oh1o-dx
esln,s
mat,4
EMODIFY,ALL
*endif
  
*if,orto,gt,0,then
LOCAL, 11, 1, BURY, OYYCH  !定义局部坐标系，柱体
CsYS,11    
NSEL,S,LOC,x,0,ORTO/2    !选择岩溶洞穴对应的节点
ESLN        !岩溶洞穴对应单元
MAT,3        !岩溶洞穴为介质三  
EMODIFY,ALL!确认修改
CSYS,0
*endif
  
ALLSEL,ALL
  
!*DO,OYX,1.5,0,-0.3
/prep7
!OYX溶洞到炮点的距离与地震波基长之比
!OY=OYYCH-OYX*LAMDA!!炮点的Y坐标
!RY1=OYych+OYX/2  !检波点1的Y坐标
!ry2=oyych-oyx/2
oxlam=0.1  !炮点的埋深与基长之比
OXO=oxlam*lamda  !炮点和检波点的埋深，X坐标
!EXPLONODE=NODE(OXO,OY,0)
!RECEINODEZ=NODE(OXO,RY1,0)
!RECEINODER=NODE(0,RY1,0)
FINISH   
  
!定义边界UX、UY数组，用于存放临近边界两层节点数据
*dim,UXL,ARRAY,2,XDEPTH/DX+1
*DIM,UXR,ARRAY,2,XDEPTH/DX+1
*DIM,UXB,ARRAY,2,YWIDTH/DX+1
*dim,UYL,ARRAY,2,XDEPTH/DX+1
*DIM,UYR,ARRAY,2,XDEPTH/DX+1
*DIM,UYB,ARRAY,2,YWIDTH/DX+1
  
/SOLU  
!定义和读入时程曲线
NSUBST,1, , ,1 !1个子步  
!OUTRES,NSOL,100 !输出每个子步的结果  
!OUTPR,BASIC,recordstep,
ANTYPE,TRANS !时程分析  
!TRNOPT,REDUC
LUMPM,0  
!*   
ALLSEL,ALL   
EPLOT
   
*do,i,1,TOTALNT
  
*if,i,gt,1,then    !获得临近边界的两层节点的UX、UY，输入到相应的数组中
*DO,LRX,1,XDEPTH/DX+1
XLR=(LRX-1)*DX
BNODE1=NODE(XLR,0,0) $ BNODE2=NODE(XLR,DX,0) $ BNODE3=NODE(XLR,DX*2,0)
TNODE1=NODE(XLR,YWIDTH,0) $ TNODE2=NODE(XLR,YWIDTH-DX,0) $ TNODE3=NODE(XLR,YWIDTH-DX*2,0)  
*GET,UXL(1,LRX),NODE,BNODE2,U,X $ *GET,UYL(1,LRX),NODE,BNODE2,U,Y
*GET,UXL(2,LRX),NODE,BNODE3,U,X $ *GET,UYL(2,LRX),NODE,BNODE3,U,Y
*GET,UXR(1,LRX),NODE,TNODE2,U,X $ *GET,UYR(1,LRX),NODE,TNODE2,U,Y
*GET,UXR(2,LRX),NODE,TNODE3,U,X $ *GET,UYR(2,LRX),NODE,TNODE3,U,Y
*ENDDO
*DO,LRX,1,YWIDTH/DX+1
XLR=(LRX-1)*DX
BNODE1=NODE(XDEPTH,XLR,0) $ BNODE2=NODE(XDEPTH-DX,XLR,0) $ BNODE3=NODE(XDEPTH-DX*2,XLR,0)
*GET,UXB(1,LRX),NODE,BNODE2,U,X $ *GET,UYB(1,LRX),NODE,BNODE2,U,Y
*GET,UXB(2,LRX),NODE,BNODE3,U,X $ *GET,UYB(2,LRX),NODE,BNODE3,U,Y
*ENDDO
*ENDIF
  
*if,i,le,nt+1,then  !在NT时间前延续脉冲子波
*DO,LRX,1,YWIDTH/DX+1   !施加边界条件
XLR=(LRX-1)*DX
explonode=NODE(0,XLR,0)
ac=-2*0.0001*(i*dt-t0)*exp(-f*(i*dt-t0)**2)!脉冲子波，马在田P22
D,EXPLONODE,Ux, ac
*if,i,eq,nt+1,then  !在NT+1时刻，删除作用在节点上的位移
DDELE,EXPLONODE,Ux
*endif
*enddo
*endif
  
TIME,i*DT
solve
  
*DO,LRX,1,XDEPTH/DX+1   !施加边界条件
XLR=(LRX-1)*DX
BNODE1=NODE(XLR,0,0) $ BNODE2=NODE(XLR,DX,0) $ BNODE3=NODE(XLR,DX*2,0)
TNODE1=NODE(XLR,YWIDTH,0) $ TNODE2=NODE(XLR,YWIDTH-DX,0) $ TNODE3=NODE(XLR,YWIDTH-DX*2,0)  
*GET,UXL1,NODE,BNODE1,U,X $ *GET,UYL1,NODE,BNODE1,U,Y
*GET,UXL2,NODE,BNODE2,U,X $ *GET,UYL2,NODE,BNODE2,U,Y
*GET,UXL3,NODE,BNODE3,U,X $ *GET,UYL3,NODE,BNODE3,U,Y
*GET,UXR1,NODE,TNODE1,U,X $ *GET,UYR1,NODE,TNODE1,U,Y
*GET,UXR2,NODE,TNODE2,U,X $ *GET,UYR2,NODE,TNODE2,U,Y
*GET,UXR3,NODE,TNODE3,U,X $ *GET,UYR3,NODE,TNODE3,U,Y
  
*IF,XLR,GT,OH1O,THEN
CUX=VP(2) $ CVY=VS(2)
*ELSE IF,XLR,EQ,OH1O,THEN
CUX=(VP(2)+VP(1))/2 $ CVY=(VS(1)+VS(2))/2
*ELSE
CUX=VP(1) $ CVY=VS(1)
*ENDIF
  
UXL1=UXL1+UXL2-UXL(1,LRX)+CUX*DT/DX*((UXL2-UXL1)-(UXL(2,LRX)-UXL(1,LRX)))
UYL1=UYL1+UYL2-UYL(1,LRX)+CVY*DT/DX*((UYL2-UYL1)-(UYL(2,LRX)-UYL(1,LRX)))
  
D,BNODE1,Ux,UXL1 $ D,BNODE1,UY,UYL1
  
UXR1=UXR1+UXR2-UXR(1,LRX)-CUX*DT/DX*((UXR1-UXR2)-(UXR(1,LRX)-UXR(2,LRX)))
UYR1=UYR1+UYR2-UYR(1,LRX)-CVY*DT/DX*((UYR1-UYR2)-(UYR(1,LRX)-UYR(2,LRX)))
  
D,TNODE1,Ux,UXR1 $ D,TNODE1,UY,UYR1
*ENDDO
  
*DO,LRX,1,YWIDTH/DX+1
XLR=(LRX-1)*DX
BNODE1=NODE(XDEPTH,XLR,0) $ BNODE2=NODE(XDEPTH-DX,XLR,0) $ BNODE3=NODE(XDEPTH-DX*2,XLR,0)
*GET,UXB1,NODE,BNODE1,U,X $ *GET,UYB1,NODE,BNODE1,U,Y
*GET,UXB2,NODE,BNODE2,U,X $ *GET,UYB2,NODE,BNODE2,U,Y
*GET,UXB3,NODE,BNODE3,U,X $ *GET,UYB3,NODE,BNODE3,U,Y
CUX=VP(2) $ CVY=VS(2)
UXB1=UXB1+UXB2-UXB(1,LRX)-CUX*DT/DX*((UXB1-UXB2)-(UXB(1,LRX)-UXB(2,LRX)))
UYB1=UYB1+UYB2-UYB(1,LRX)-CVY*DT/DX*((UYB1-UYB2)-(UYB(1,LRX)-UYB(2,LRX)))
  
D,BNODE1,Ux,UXB1 $ D,BNODE1,UY,UYB1
*ENDDO
  
*if,mod(i,recordstep),eq,0,then
*DO,LRX,1,YWIDTH/DX+1   !在每一时间步读出检波点UX、UY输入到DISP.DAT文件中
XLR=(LRX-1)*DX
Nnode=NODE(0,XLR,0)   
*GET,ury,node,nnode,u,y
*GET,urx,node,NNODE,u,x
*VWRITE,nnode,i,Urx,ury
(2f6.0,2e15.6)   
*enddo
  
*DO,LRX,1,YWIDTH/DX+1     !在每一时间步读出检波点UX、UY输入到DISP.DAT文件中
XLR=(LRX-1)*DX
Nnode=NODE(DX,XLR,0)
*GET,uzy,node,NNODE,u,y
*GET,uzx,node,NNODE,u,x
*VWRITE,nnode,i,Uzx,uzy
(2f6.0,2e15.6)  
*enddo
  
*DO,LRX,1,YWIDTH/DX+1     !在每一时间步读出检波点UX、UY输入到DISP.DAT文件中
XLR=(LRX-1)*DX
Nnode=NODE(OH1O,XLR,0)
*GET,uzy,node,NNODE,u,y
*GET,uzx,node,NNODE,u,x
*VWRITE,nnode,i,Uzx,uzy
(2f6.0,2e15.6)  
*enddo
  
*DO,LRX,1,YWIDTH/DX+1     !在每一时间步读出检波点UX、UY输入到DISP.DAT文件中
XLR=(LRX-1)*DX
Nnode=NODE(BURY,XLR,0)
*GET,uzy,node,NNODE,u,y
*GET,uzx,node,NNODE,u,x
*VWRITE,nnode,i,Uzx,uzy
(2f6.0,2e15.6)  
*enddo
*endif
  
*enddo !进入下一时间步
  
fini  
save
  
*CFCLOSE
/POST1   
*do,i,1,TOTALNT
*if,i,eq,NINT(nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt1x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt1y, jpg,  
*endif
*if,i,eq,NINT(2*nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt2x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt2y, jpg,  
*endif
*if,i,eq,NINT(3*nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt3x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt3y, jpg,  
*endif
*if,i,eq,NINT(4*nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt4x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt4y, jpg,  
*endif
*if,i,eq,NINT(5*nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt5x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt5y, jpg,  
*endif
*if,i,eq,NINT(6*nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt6x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt6y, jpg,  
*endif
*if,i,eq,NINT(7*nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt7x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt7y, jpg,  
*endif
*if,i,eq,NINT(8*nt),then
set,i $ plnsol,U,X,0,1
/image, save, nt8x, jpg,  
set,i $ plnsol,U,y,0,1
/image, save, nt8y, jpg,  
*endif
  
*enddo !进入下一时间步

sl3913

给点积分呗!

Jeffjia

你把你的结果图形化，我想看看怎么个探测的，结果如何，让大家看看，我给你老哥申请积分！

duke3486

好象不能运行

sl3913

对不起，我最近没有上网，没能及时发贴，让您久等了！
很高兴收到你的回复！
首先，计算没问题，肯定能行的，我算了一些例子，现在看命令流文件写得不太好，用APDL功能就方便多了，可根据实际情况改。
至于，结果图像，就是地表各点的位移矢量历时曲线，我用matlab生成的。有明显的绕射特征，应该是洞穴的反映！可以试着改变洞穴位置，绕射点位置也会随着变化！
祝：健康