h=180				!大壩高度
fini_days=30			!完成大壩澆筑的總的天數
ht=h/fini_days			!每天完成的大壩高度
h1=100				!上游水位高度
h2=80				!下游水位高度
h3=100				!壩體不同混凝土材料的分界線，為100m
gm=1000				!gm為水的容重
fia1=90-atan(0.75)*180/3.1415926	!計算下游斜面夾角
fia2=90-fia1			!計算坐標旋轉角度
/prep7
et,1,plane42			!輔助劃分網格之用
et,2,solid65			!混凝土單元
et,3,solid45			!巖石單元	
keyopt,2,5,2			!輸出solid65單元的節點應力

!1號材料
mp,ex,1,2.85e10			!100m以下混凝土的材料特性 
MP,PRXY,1,0.167 
TB,CONC,1,1,9, 
TBDATA,,0.3,1,1.96e6,22e6 	!張開剪切傳遞系數為0.3，抗拉強度fc=1.96e6pa,抗壓強度fs=22e6pa,
mp,dens,1,2400

!2號材料
mp,ex,2,2.6e10			!100m以上混凝土的材料特性 
MP,PRXY,2,0.167 
TB,CONC,2,1,9, 
TBDATA,,0.3,1,1.2e6,17.5e6 	!張開剪切傳遞系數為0.3，抗拉強度fc=1.2e6pa,抗壓強度fs=17.5e6pa,
mp,dens,2,2400

!3號材料
mp,ex,3,2.9e10
mp,PRXY,3,0.3
mp,dens,3,2600

!創建幾何模型
!創建關鍵點
k,1
k,2,0.9*0.75*h+0.1*h
k,3,0.1*h,0.9*h
k,4,0,0.9*h
!創建線
l,1,2
l,2,3
l,3,4
l,4,1
!創建面
al,1,2,3,4
RECTNG,0,0.1*h,0.9*h,h				!畫壩頂矩形
RECTNG,-1.5*h,0,-2*h,0				!壩基左面矩形
RECTNG,0,0.9*0.75*h+0.1*h,-2*h,0		!壩基正下方矩形
RECTNG,0.9*0.75*h+0.1*h,(2+0.9*0.75)*h,-2*h,0	!壩基右面矩形
RECTNG,-1.5*h,(2+0.9*0.75)*h,-2*h,h		!總矩形
AOVLAP,all					!將面單元進行粘貼布爾操作
/pnum,line,1					!打開面，線的號碼開關
/pnum,area,1
nummrg,all					!合并重復元素
numcmp,all
aplot

!劃分3號面的網格,壩體的頂部
lsel,s,,,3,5,2
lesize,all,,,5
lsel,s,,,12,13,1
lesize,all,,,2
amesh,3

!劃分1號面的網格，壩體的下部
lsel,s,,,2,4,2
lesize,all,,,18
lsel,s,,,1
lesize,all,,,5
amesh,1
eplot

!劃分2號面的網格，壩體的正下方基巖
lsel,s,,,11
lesize,all,,,5
lsel,s,,,9,10,1
lesize,all,,,8,4
amesh,2

!劃分4號面的網格
lsel,s,,,14
lesize,all,,,5,4
lsel,s,,,7
lesize,all,,,8,4
lsel,s,,,6
lesize,all,,,5,0.25
amesh,4
eplot

!劃分5號面的網格
lsel,s,,,15,16,1
lesize,all,,,8,4
lsel,s,,,8
lesize,all,,,8,0.25
amesh,5
eplot

!連接4號線和13號線，并劃分7號面的網格
lsel,s,,,4,13,9
LCCAT,all
lsel,s,,,19
lesize,all,,,5,4
lsel,s,,,20,21,1
lesize,all,,,20
amesh,7
eplot

!連接2號線和12號線，并劃分6號面的網格
lsel,s,,,2,12,10
lccat,all
lsel,s,,,17
lesize,all,,,8,4
lsel,s,,,18,22,4
lesize,all,,,20
amesh,6
eplot

!刪除前面連接的線元素
allsel
lsel,r,lcca
ldele,all

!拉伸成壩體單元，采用solid65單元和1號材料，隨后改變h3以上壩體的材料屬性
EXTOPT,ESIZE,8,0,   			!拉伸的份數     
type,2
mat,1
VEXT,1,3,2,,,-0.75*h			!將1號和3號面拉伸成一半的模型
/view,1,1,1,1         
vplot
eplot

!壩體底部單元，拉伸成巖土單元，采用solid45號單元和3號材料
EXTOPT,ESIZE,8,0,        
type,3
mat,3
VEXT,2,,,,,-0.75*h
VEXT,4,5,1,,,-0.75*h
eplot

!拉伸生成大壩兩側的巖石單元
allsel
EXTOPT,ESIZE,5,4,        
EXTOPT,ACLEAR,1				!刪除面單元
type,3
mat,3
VEXT,1,7,1,,,h		
/pnum,mat,1
eplot

!利用對稱方法完成整個有限元模型
local,11,0,,,-0.75*h			!定義局部坐標系
csys,11					!激活局部坐標系
/PSYMB,CS,1				!顯示不同的局部坐標系符號
DSYS,11					!顯示局部坐標系
NSYM,Z,20000,all
ENSYM,30000, ,20000,all
nummrg,all			!合并重復節點和單元
numcmp,all
eplot

!改變材料屬性
esel,s,mat,,1			!選擇材料編號為1的單元
eplot
nsle,s				!選擇單元上的節點
nplot
nsel,r,loc,y,h3+1,h+1		!選擇h3以上的單元，以便轉換材料屬性
nsel,r,loc,z,-0.75*h+0.1,0.75*h-0.1
nplot
esln,s
eplot
MPCHG,2,all,			!改變材料屬性為2號材料
/REPLOT 
allsel
eplot
save,dam_model,db		!保存模型文件
finish

/solu
!施加邊界條件
csys,0
dsys,0
nsel,s,loc,x,(2+0.9*0.75)*h
nsel,a,loc,x,-1.5*h
nplot
d,all,ux
allsel

nsel,s,loc,z,h
nsel,a,loc,z,-(1+1.5)*h
nplot
d,all,uz
allsel

nsel,s,loc,y,-2*h
nplot
d,all,uy
allsel
gplot

!施加重力
ACEL,0,9.8,0,

!選擇上游壩面單元，施加水平靜水壓力
ESEL,S,TYPE,,2
eplot
nsel,s,loc,x,0
nsel,r,loc,z,-1.5*h+0.1,-0.1
nsel,r,loc,y,0.1,h1-5		!若采用(h1-0.1)將導致選擇過多的單元,
				!使得施加的面分布荷載出現負數的情況
esln,s
nplot
eplot
/PSF,PRES,NORM,2,0,1		!設定將顯示壓力的方向
sfgrad,pres,0,y,0,-gm		!給定載荷梯度為-gm
sfe,all,2,pres,,gm*h1
/replot

!下游壩面單元，施加水平靜水壓力
allsel
ESEL,S,TYPE,,2
nsle,s
eplot
local,12,0,0.9*0.75*h+0.1*h,,,fia2
csys,12
/PSYMB,CS,1
DSYS,12

nsel,s,loc,y,0.1,h2/sin(0.75)-25	!若采用(h2/sin(0.75)-1)將導致選擇過多的單元,
					!使得施加的面分布荷載出現負數的情況
nsel,r,loc,z,-1.5*h+0.1,-0.1
nsel,u,loc,x,-1000,-2
esln,s
nplot
eplot
/PSF,PRES,NORM,2,0,1		!設定將顯示壓力的方向
sfgrad,pres,0,y,0,-gm		!給定載荷梯度為-gm
sfe,all,4,pres,,gm*h2
/replot

!施加上游水位h1產生的滲透壓力
!同時施加下游水位h2產生的浮托力
!選擇第一部分壩底單元
csys,0
dsys,0
/psymb,cs,1
allsel
nsel,s,loc,y,0
nsel,r,loc,z,-135+0.1,0.1
esln,s
esel,r,type,,2
eplot
SFCUM,PRES,ADD			!設置載荷是疊加的
sfe,all,5,pres,,gm*h2		!下游水位產生的浮托力
p0=gm*h1/(0.9*0.75*h+0.1*h)
sfgrad,pres,0,x,0,-p0
sfe,all,5,pres,,gm*h1
/replot
!選擇第二部分壩底單元
allsel
nsel,s,loc,y,0
nsel,r,loc,z,-270+0.1,-135-0.1
esln,s
esel,r,type,,2
eplot
p0=0
sfgrad,pres,0,x,0,-p0		!令斜率為0
sfe,all,3,pres,,gm*h2		!下游水位產生的浮托力
p0=gm*h1/(0.9*0.75*h+0.1*h)
sfgrad,pres,0,x,0,-p0
sfe,all,3,pres,,gm*h1
SFCUM,PRES,repl			!設置載荷是替代的

!檢查壩底的載荷是否正確
nsel,s,loc,y,0
esln,s
esel,r,type,,2
/replot
save,dam_static_solu,db
allsel
outres,all,all
autots,on
nsubst,20
solve
save,dam_static_rst,db

/post1
esel,s,type,,2

set,1,last			!獲取第一載荷步，最后子步的結果
PLNSOL,U,X
PLNSOL,U,Y
PLNSOL,EPTO,1,0,1		!畫第一主應變
PLNSOL,S,1,0,1			!畫第一主應力
PLNSOL,S,3,0,1			!畫第三主應力
/DEVICE,VECTOR,1
PLCRACK,0,0			!顯示開裂位置
!察看在各個積分點上的開裂情況
ETABLE,11,NMISC,53		!顯示開裂的情況，包括裂縫的張開，閉合，壓碎
ETABLE,22,NMISC,60
ETABLE,33,NMISC,67
ETABLE,44,NMISC,74
ETABLE,55,NMISC,81
ETABLE,66,NMISC,88
ETABLE,77,NMISC,95
ETABLE,88,NMISC,102
PLETAB,33,AVG   
PLETAB,44,AVG   
finish