《FLAC/FLAC 3D基础与工程实例》全部命令流1gen zone bri p0 0 0 0 p1 10 0 0 p2 0 10 0 p3 0 0 10 &p4 15 15 0 p5 0 15 15 p6 15 0 10 p7 20 20 20 &size 10 10 10 rat 1.0 0.9 1.1 group brick_1gen zone bri p0 20 0 0 p1 add 10 0 0 p2 add 0 20 0 p3 add 0 0 15 &size 10 10 10 rat 1.0 0.9 1.1 group brick_2gen zone bri p0 40 0 0 edge 10 size 10 10 10 rat 1.0 0.9 1.1 group brick_3plot sur2-1newgen zon bri size 3 3 3model elasprop bulk 3e8 shear 1e8ini dens 2000fix z ran z -.1 .1fix x ran x -.1 .1fix x ran x 2.9 3.1fix y ran y -.1 .1fix y ran y 2.9 3.1set grav 0 0 -10solveapp nstress -10e4 ran z 3 x 1 2 y 1 2hist gp vel 0 0 3hist gp vel 0 3 3plo hist 1 redplo add hist 2 bluesolve3-1;-------------------------------------工程信息;Project Record Tree export;Title:Simple test;---------------------------------计算第一步;... STATE: STATE1 ....configgrid 10,10model elastic group 'User:Soil' notnullmodel elastic notnull group 'User:Soil'prop density=1500.0 bulk=3E6 shear=1E6 notnull group 'User:Soil'fix x y j 1fix x i 1fix x i 11set gravity=9.81history 999 unbalancedsolvesave state1.sav;----------------------------------计算第二步;... STATE: STATE2 ....initial xdisp 0 ydisp 0initial xvel 0 yvel 0model null i 4 7 j 8 10group 'null' i 4 7 j 8 10group delete 'null'history 1 xdisp i=4, j=11solvesave state2.sav;--------------------------------绘图命令;*** plot commands ****;plot name: syyplot hold grid syy fill;plot name: Unbalanced forceplot hold history 999;plot name: gridplot hold grid magnify 20.0 lred grid displacement;plot name: Xdis-Aplot hold history 1 line5-1new; ===============================; 定义球体半径和半径方向上单元网格数; ===============================def parmrad=10.0rad_size=5endparm; ===============================; 建立八分之一球体外接立方体网格; ===============================gen zone pyramid p0 rad 0 0 p1 rad 0 rad p2 rad rad 0 p30 0 0 &p4 rad rad rad size rad_size rad_size rad_size group 1 gen zone pyramid p0 0 rad 0 p1 rad rad 0 p2 0 rad rad p3 0 0 0 &p4 rad rad rad size rad_size rad_size rad_size group 2 gen zone pyramid p0 0 0 rad p1 0 rad rad p2 rad 0 rad p3 0 0 0 &p4 rad rad rad size rad_size rad_size rad_size group 3; ==================================; 利用FISH语言将内部立方体节点调整到球面; ==================================def make_spherep_gp=gp_headloop while p_gp#null; 获取节点点坐标值：P=(px,py,pz)px=gp_xpos(p_gp)py=gp_ypos(p_gp)pz=gp_zpos(p_gp)dist=sqrt(px*px+py*py+pz*pz)if dist>0 then; 节点位置调整maxp=max(px,max(py,pz))k=(maxp/rad)*(rad/dist)gp_xpos(p_gp)=k*pxgp_ypos(p_gp)=k*pygp_zpos(p_gp)=k*pzend_ifp_gp=gp_next(p_gp)end_loopendmake_sphere; ===============================; 利用镜像生成完整球体网格; ===============================gen zone refgen zone ref dip 90gen zone ref dip 90 dd 90; ===============================; 显示球体网格; ===============================plot surfpl set back whpl bl gr 5-3ngen zon bri size 1 1 2group soil ran z 1 1group rock ran z 0 1expgrid 1.flac3d6-1gen zon bri size 3 3 3model mohrprop bu 3e6 sh 1e6 coh 10e3 fric 15fix z ran z -.1 .1fix x ran x -.1 .1fix x ran x 2.9 3.1fix y ran y -.1 .1fix y ran y 2.9 3.1ini dens 2000hist unbalset grav 10solve elasticsave 6-1.sav6-2rest 6-1.savini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 solvesave 6-2.sav6-3rest 6-1.savini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 hist id=2 gp zdis 1 1 3hist id=3 gp zdis 1 1 2hist id=4 gp xdis 1 1 3hist id=5 gp xdis 1 1 3hist id=6 zone szz 1 1 3hist id=7 zone szz 1.5 1.5 2.5hist id=8 zone sxz 1.5 1.5 2.5solvesave 6-3.sav6-4rest 6-3.savset log onset logfile 6-2.logprint zone stressprint gp disset log off6-5rest 6-1.savini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 plot set rot 20 0 30plot con szz ou on magf 10plot add hist 1set movie avi step 1 file 6-5.avimovie startsolvemovie finish7-1ngen zon bri size 1 1 2model elasprop bulk 3e7 shear 1e7fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000set grav 0 0 -10solveplo con sz7-2ngen zon bri size 1 1 2 model mohrprop bulk 3e7 shear 1e7 c 1e10 f 15 tension 1e10 fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000set grav 0 0 -10solveprop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0solveplo con sz7-3ngen zone brick size 1 1 2model mohrprop bulk 3e7 shear 1e7 coh 10e3 fri 15 ten 0fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000set grav 0 0 -10solve elasplo con sz7-4newgen zone brick size 1 1 2model mohrprop bulk 3e7 shear 1e7 coh 10e3 fri 15 ten 0fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000ini szz -40e3 grad 0 0 20e3 ran z 0 2ini syy -20e3 grad 0 0 10e3 ran z 0 2ini sxx -20e3 grad 0 0 10e3 ran z 0 2set grav 0 0 -10solveplo con sz7-5ngen zon bri size 1 1 2model mprop bulk 3e7 shear 1e7 c 10e10 f 15 ten 1e10 fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000 ran z 0 1ini dens 1500 ran z 1 2ini szz -35e3 grad 0 0 20e3 ran z 0 1ini syy -17.5e3 grad 0 0 10e3 ran z 0 1ini sxx -17.5e3 grad 0 0 10e3 ran z 0 1ini szz -15e3 grad 0 0 15e3 ran z 1 2ini syy -7.5e3 grad 0 0 7.5e3 ran z 1 2ini sxx -7.5e3 grad 0 0 7.5e3 ran z 1 2ini pp 10e3 grad 0 0 -10e3 ran z 0 1set grav 0 0 -10solveplo con sz7-6ngen zon bri size 1 1 2model mprop bulk 3e7 shear 1e7 c 10e10 f 15 ten 1e10 fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000 ran z 0 2ini szz -50e3 grad 0 0 20e3 ran z 0 1ini syy -30e3 grad 0 0 10e3 ran z 0 1ini sxx -30e3 grad 0 0 10e3 ran z 0 1ini pp 30e3 grad 0 0 -10e3 ran z 0 2app nstress -10e3 ran z 2 set grav 0 0 -10solveplo con sz7-7newgen zone brick p0 0 0 0 p1 60 0 0 p2 0 60 0 p3 0 0 90 & p4 60 60 0 p5 0 60 90 p6 60 0 150 p7 60 60 150 & size 6 6 10model elaspro bulk 10e10 she 10e10ini den 2500apply sxx -1e9 grad 0 0 1.1111111e7 range x -.1 .1 apply sxx -1e9 grad 0 0 6.6666666e6 range x 59.9 60.1 apply syy -1e9 grad 0 0 8.3333333e6 range y -.1 .1 apply syy -1e9 grad 0 0 8.3333333e6 range y 59.9 60.1 apply szz -1e8 grad 0 0 8.3333333e5 ran z 0 120set grav 0 0 -10step 30000ini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0plo cont szz7-8newgen zone brick p0 0 0 0 p1 60 0 0 p2 0 60 0 p3 0 0 90 & p4 60 60 0 p5 0 60 90 p6 60 0 150 p7 60 60 150 & size 6 6 10model elaspro bulk 10e10 she 10e10ini den 2500ini sxx -1e9 grad 0 0 1.1111111e7 range x -.1 .1ini sxx -1e9 grad 0 0 6.6666666e6 range x 59.9 60.1ini syy -1e9 grad 0 0 8.3333333e6 range y -.1 .1ini syy -1e9 grad 0 0 8.3333333e6 range y 59.9 60.1ini szz -1e8 ran z -.1 .1fix x y z ran z -.1 .1set grav 0 0 -10solveini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0plo cont szz7-9newgen zone brick p0 0 0 -50 p1 27.5 0 -50 p2 0 5 -50 p3 0 0 -10 size 8 1 10 group claygen zone brick p0 27.5 0 -50 p1 100 0 -50 p2 27.5 5 -50 p3 27.5 0 -10 ratio 1.1 1 1 size 12 1 10 group claygen zone brick p0 0 0 -10 p1 27.5 0 -10 p2 0 5 -10 p3 0 0 0 ratio 1 1 0.8 size 8 1 4 group soilgen zone brick p0 27.5 0 -10 p1 100 0 -10 p2 27.5 5 -10 p3 27.5 0 0 ratio 1.1 1 0.8 size 12 1 4 group soilgen zone brick p0 0 0 0 p1 27.5 0 0 p2 0 5 0 p3 0 0 5 p4 27.5 5 0 &p5 0 5 5 p6 20 0 5 p7 20 5 5 size 8 1 5 group damfix x y z ran z -49.9 -50.1fix x ran x -.1 .1fix x ran x 99.9 100.1fix ymodel mohr ran z -50 0model null ran z 0 5prop bulk 7.8e6 shear 3.0e6 coh 10e10 tension 1e10 ran group soilini dens 1500 ran group soilprop bulk 3.91e6 shear 1.5e6 coh 10e10 tension 1e10 ran group clayini dens 1800 ran group clayset grav 0 0 -9.8hist id=1 unbalsolveprop bulk 7.8e6 shear 3.0e6 coh 10e3 fric 15 ran group soilprop bulk 3.91e6 shear 1.5e6 coh 20e3 fric 20 ran group claysolvesave elastic.savini xdis 0 ydis 0 zdis 0 ;将节点位移清零ini xvel 0 yvel 0 zvel 0 ;将节点速度清零hist id=2 gp zdis 0 0 0 ;记录地基顶部中心点的沉降hist id=3 gp zdis 27.5 0 0 ;记录路基坡脚处的沉降hist id=4 gp xdis 27.5 0 0 ;记录路基坡脚处的水平位移model elastic ran z 0 1 ; ;激活0 m ~ 1 m的单元prop bulk 7.8e6 shear 3.0e6 ran z 0 1ini dens 1500 ran z 0 1solve ;按软件默认精度求解save fill-1.savmodel elastic ran z 1 2prop bulk 7.8e6 shear 3.0e6 ran z 1 2ini dens 1500 ran z 1 2solvesave fill-2.savmodel elastic ran z 2 3prop bulk 7.8e6 shear 3.0e6 ran z 2 3ini dens 1500 ran z 2 3solvesave fill-3.savmodel elastic ran z 3 4prop bulk 7.8e6 shear 3.0e6 ran z 3 4ini dens 1500 ran z 3 4solvesave fill-4.savmodel elastic ran z 4 5prop bulk 7.8e6 shear 3.0e6 ran z 4 5ini dens 1500 ran z 4 5solvesave fill-5.savpau;plo bl gr;gen zone brick p0 0 0 0 p1 100 0 0 p2 0 5 0 p3 0 0 5 sizegen zone brick p0 0 0 -50 p1 27.5 0 -50 p2 0 5 -50 p3 0 0 -10 size 8 1 10 group claygen zone brick p0 27.5 0 -50 p1 100 0 -50 p2 27.5 5 -50 p3 27.5 0 -10 ratio 1.1 1 1 size 12 1 10 group clay gen zone brick p0 0 0 -10 p1 27.5 0 -10 p2 0 5 -10 p3 0 0 0 ratio 1 1 0.8 size 8 1 4 group soilgen zone brick p0 27.5 0 -10 p1 100 0 -10 p2 27.5 5 -10 p3 27.5 0 0 ratio 1.1 1 0.8 size 12 1 4 group soil gen zone brick p0 0 0 0 p1 27.5 0 0 p2 0 5 0 p3 0 0 5 p4 27.5 5 0 &p5 0 5 5 p6 20 0 5 p7 20 5 5 size 8 1 5 group dam7-10set log on ;打开log记录set logfile 1.log ;设置记录文件名为：1.log restore fill-1.sav ;调用保存的文件print gp dis range id 517 any id 533 any ;输出两个节点的变形值restore fill-2.savprint gp dis range id 517 any id 533 anyrestore fill-3.savprint gp dis range id 517 any id 533 anyrestore fill-4.savprint gp dis range id 517 any id 533 anyrestore fill-5.savprint gp dis range id 517 any id 533 anyset log off ;关闭log记录8-1def abcabc = 1 + 2 * 3abcd = 1.0 / 2.0endabcprint fish8-2def abcif aa < 0 thenabc = 0.0elseabc = 2.0 * aaendifendabc8-3def abcloop aa (1, 2.5)commandprint aaendcommandendloopend8-4newgen zon bri size 3 3 3model elasticprop bu 3e7 sh 1e7ini dens 2000fix x y z ran z -.1 .1fix x ran x -.1 .1fix x ran x 2.9 3.1fix y ran y -.1 .1fix y ran y 2.9 3.1set grav 10solveini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0save 8-4.sav8-5rest 8-4.savdef E_modifyp_z = zone_headd_k = 704d_n = 0.38d_pa = 101325.0 ;//标准大气压loop while p_z # nullsigma_3 = -1.0 * z_sig1(p_z)E_new = d_k * d_pa * (sigma_3 / d_pa) ^ d_n z_prop(p_z,'young') = E_newp_z = z_next(p_z)endloopendE_modify8-6rest 8-5.savtable 1 name load_settlementdef add_loadp_gp = gp_near(2,1,3)loop n (1,5)app_load = n * (-1000e3)file_name = '7-6_add_step' + string(n) + '.sav'commandapp nstress app_load ran z 2.9 3.1 x 1 2 y 1 2solvesave file_nameendcommandxtable(1,n) = -1.0 * app_loadytable(1,n) = gp_zdisp(p_gp)endloopendadd_loadsave 8-6.sav8-7rest 8-6.savdef find_max_dispp_gp = gp_headmaxdisp_value = 0.0maxdisp_gpid = 0loop while p_gp # nulldisp_gp = sqrt(gp_xdisp(p_gp) ^ 2 + gp_ydisp(p_gp) ^ 2 + gp_zdisp(p_gp) ^ 2)if disp_gp > maxdisp_valuemaxdisp_value = disp_gpmaxdisp_gpid = gp_id(p_gp)endifp_gp = gp_next(p_gp)endloopendfind_max_dispprint maxdisp_value maxdisp_gpidrest 8-6.savconfig zextra 1def get_sigma_difp_z = zone_headloop while p_z # nullsigma_dif = z_sig3(p_z) - z_sig1(p_z)z_extra(p_z,1) = sigma_difp_z = z_next(p_z)endloopendget_sigma_difplot con zextra 19-1 ;---------------------------------------------------;; 移来移去法接触面的建立;---------------------------------------------------ngen zone radcyl p0 (0,0,0) p1(8,0,0) p2 (0,0,-5) p3 (0,8,0) &p4 (8,0,-5) p5 (0,8,-5) p6 (8,8,0) p7 (8,8,-5) &p8 (.3,0,0) p9 (0,.3,0) p10 (.3,0,-5) p11 (0,.3,-5) &size 3 10 6 15 ratio 1 1 1 1.15gen zone radcyl p0 (0,0,-5) p1 (8,0,-5) p2 (0,0,-8) p3 (0,8,-5) &p4(8,0,-8) p5 (0,8,-8) p6 (8,8,-5) p7 (8,8,-8) &p8 (.3,0,-5) p9 (0,.3,-5) p10 (.3,0,-8) p11 (0,.3,-8) &size 3 6 6 15 ratio 1 1 1 1.15 fillgen zone reflect dd 270 dip 90group clay;interface 1 face range cylinder end1 (0,0,0) end2 (0,0,-5.1) radius .31 &cylinder end1 (0,0,0) end2 (0,0,-5.1) radius .29 not interface 2 face range cylinder end1 (0,0,-4.9) end2 (0,0,-5.1) radius .31;pausegen zone cyl p0 (0,0,6) p1 (.3,0,6) p2 (0,0,1) p3 (0,.3,6) &p4 (.3,0,1) p5 (0,.3,1) &size 3 10 6gen zone cyl p0 (0,0,6.1) p1 (.3,0,6.1) p2 (0,0,6) p3 (0,.3,6.1) &p4 (.3,0,6) p5 (0,.3,6) &size 3 1 6gen zone reflect dd 270 dip 90 range z 1 6.1group pile range z 1 6.1pauseini z add -6.0 range group pilesave pile_geom.sav9-2;导来导去法;---------------------------------------------ngen zone brick size 3 3 3group 2 range x 1 2 y 1 2 z 1 2group 1 range gr 2 notsave 1.savdel ran group 2 notinterface 1 face range x 1 y 1 2 z 1 2interface 1 face range x 2 y 1 2 z 1 2interface 1 face range x 1 2 y 1 z 1 2interface 1 face range x 1 2 y 1 2 z 1interface 1 face range x 1 2 y 1 2 z 2rest 1.savdel ran group 2expgrid 1.fac3dimpgrid 1.flac3dmodel elafix x y z ran z 0ini den 2000set grav 0 0 -10interface 1 prop kn 20e6 ks 20e6 coh 10e3 fri 15app nstr -200e3 ran x 0 1 y 1 2 z 3solve9-3;切割模型法;------------------------------ngen zone brick size 3 3 3group 1 range x 1 2 y 1 2 z 2 3group 2 range group 1 notgen separate 1int 1 wrap 1 2int 1 maxedge 0.5plo int red9-4; Create Material Zonesgen zone brick size 5 5 5 &p0 (0,0,0) p1 (3,0,0) p2 (0,3,0) p3 (0,0,5) &p4 (3,3,0) p5 (0,5,5) p6 (5,0,5) p7 (5,5,5) gen zone brick size 5 5 5 p0 (0,0,5) edge 5.0group Material; Create Bin Zonesgen zone brick size 1 5 5 &p0 (3,0,0) p1 add (3,0,0) p2 add (0,3,0) &p3 add (2,0,5) p4 add (3,6,0) p5 add (2,5,5) &p6 add (3,0,5) p7 add (3,6,5)gen zone brick size 1 5 5 &p0 (5,0,5) p1 add (1,0,0) p2 add (0,5,0) &p3 add (0,0,5) p4 add (1,6,0) p5 add (0,5,5) &p6 add (1,0,5) p7 add (1,6,5)gen zone brick size 5 1 5 &p0 (0,3,0) p1 add (3,0,0) p2 add (0,3,0) &p3 add (0,2,5) p4 add (6,3,0) p5 add (0,3,5) &p6 add (5,2,5) p7 add (6,3,5)gen zone brick size 5 1 5 &p0 (0,5,5) p1 add (5,0,0) p2 add (0,1,0) &p3 add (0,0,5) p4 add (6,1,0) p5 add (0,1,5) &p6 add (5,0,5) p7 add (6,1,5)group Bin range group Material not; Create named range synonymsrange name=Bin group Binrange name=Material group Material; Assign models to groupsmodel mohr range Materialmodel elas range Bingen separate Materialinterface 1 wrap Material Bin range plane ori 0 0 0 normal 1 -1 0 aboveinterface 2 wrap Material Bin range plane ori 0 0 0 normal 1 -1 0 belowint 1 maxedge 0.55int 2 maxedge 0.55; Assign propertiesprop shear 1e8 bulk 2e8 fric 30 range Materialprop shear 1e8 bulk 2e8 range Binini den 2000int 1 prop ks 2e9 kn 2e9 fric 15int 2 prop ks 2e9 kn 2e9 fric 15; Assign Boundary Conditionsfix x range x -0.1 0.1 any x 5.9 6.1 anyfix y range y -0.1 0.1 any y 5.9 6.1 anyfix z range z -0.1 0.1 Bin; Monitor historieshist unbalhist gp zdisp (6,6,10)hist gp zdisp (0,0,10)hist gp zdisp (0,0,0); Settingsset largeset grav 0,0,-10; Cyclingstep 4000save bin.sav9-5;;------------------------------------------------------------------ rest pile_geom.savmodel mohr range group clayprop bulk 8.333e7 shear 3.846e7 coh 30000 fric 0 range group claymodel elas range group pile prop bulk 8.333e7 shear 3.846e7 range group pile interface 1 prop kn 1e8 ks 1e8 fric 20 coh 30000 interface 2 prop kn 1e8 ks 1e8 fric 20 coh 30000;ini dens 1230 range group clayini dens 1230 range group pilemodel null range z -0.1 0.15;fix z range z -8.1 -7.9fix x range x -8.1 -7.9fix x range x 7.9 8.1fix y range y -.1 .1fix y range y 7.9 8.1set grav 0 0 -10ini szz 0. grad 0 0 12300. range z -5.5 0.ini szz 17600 grad 0 0 15500 range z -8 -5.5ini sxx 0. grad 0 0 5271.4 range z -5.5 0.ini sxx 7542.86 grad 0 0 6642.86 range z -8 -5.5ini sxx add 31428.6 grad 0 0 5714.3 range z -8 -5.5ini syy 0. grad 0 0 5271.4 range z -5.5 0.ini syy 7542.86 grad 0 0 6642.86 range z -8 -5.5ini syy add 31428.6 grad 0 0 5714.3 range z -8 -5.5; water density 1000water table origin 0,0,-5.5 normal 0 0 -1ini dens 1550 range z -8 -5.5hist unbal;solve rat 1.e-6save pile0.sav;model elas range group pile prop bulk 13.9e9 shear 10.4e9 range group pileini dens 2500 range group pile call find_add.fissolve rat 1.e-6save pile1.savrest pile1.sav ;调用保存文件ini state 0ini xdis 0.0 ydis 0.0 zdis 0.0 ;位移清零apply szz -0.4e6 range z 0.05 0.15 group pile ;桩顶加第一级荷载solvesave app0.4.savprint gp disp range id 1 ;输出第一级荷载下的桩顶位移，假定桩顶中心的id号为1apply szz -0.6e6 range z 0.05 0.15 group pile ;桩顶加第二级荷载solvesave app0.6.savprint gp disp range id 1 ;输出第二级荷载下的桩顶位移…………………………………………;依次加载，直到桩破坏;-------------------------------------------------------------------- ;速度加载法rest pile1.savini state 0ini xdis 0 ydis 0 zdis 0def zs_top ;检测桩顶竖向荷载ad = top_headzftot = 0.0loop while ad # nullgp_pnt = mem(ad+1)zf = gp_zfunbal(gp_pnt)zftot = zftot + zfad = mem(ad)endloopzs_top = zftot / 0.1414endfix z range z 0.05 .15 group pile ;固定桩顶速度，用速度来确定位移def rampwhile_steppingif step < ncut thenudapp = float(step) * udmax / float(ncut) elseudapp = udmaxendifad = top_headloop while ad # nullgp_pnt = mem(ad+1)gp_zvel(gp_pnt) = udappad = mem(ad)endloopendhist gp zdis 0,0,0hist gp zvel 0,0,0hist zs_tophist zone szz 0,0,-.1set mech damp combset udmax = -1e-8 ncut 30000step 225000save pile2.sav;;--------------------------------------------------------------------;位移控制法def solve_stepsloop n (1,21)save_file = string(n) + '-step.sav'commandstep 40000save save_filepri zone stress ran id 2381 a id 2361 a id 2341 a ;输出桩顶网格单元的应力endcommandendloopendsolve_steps ;-----------------------------------------------------------------------------10-1sel beam beg 0 0 0 end 2 0 0 nseg 2sel beam beg 2 0 0 end 4 0 -1 nseg 3sel beam id=2 beg 4 0 -1 end 5 0 -2 nseg 2plot sel geo id on nod on scale 0.04plot ad ax10-2sel node id=1 0 0 0sel node id=2 2 0 0sel node id=3 4 0 -1sel node id=4 5 0 -2sel beamsel id=1 cid=1 node 1 2sel beamsel id=1 cid=2 node 2 3sel beamsel id=1 cid=3 node 3 4plot sel geo id on nod on scale 0.04plot ad ax10-3def set_valsptA = 25.0 * sin( 40.0*degrad )ptB = 25.0 * cos( 40.0*degrad )endset_valsgen zone cylinder p0=( 0.0, 0.0, 0.0 ) &p1=( ptA, 0.0, ptB ) &p2=( 0.0, 25.0, 0.0 ) &p3=( 0.0, 0.0, 25.0 ) &p4=( ptA, 25.0, ptB ) &p5=( 0.0, 25.0, 25.0 ) &size=(1, 2, 2)sel shell id=5 range cylinder end1=(0.0, 0.0,0.0) &end2=(0.0,25.0,0.0) radius=24.5 notplot blo groplot ad sel geom black black cid on scale=0.03plot ad axpaudelete ; delete all zonessel node init zpos add -25.010-4newgen zone brick size 6 8 8model mohrprop bulk 1e8 shear 0.3e8 fric 35prop coh 1e10 tens 1e10set grav 0 0 -9.81ini dens 1000fix x range x -0.1 0.1fix x range x 5.9 6.1fix y range y -0.1 0.1fix y range y 7.9 8.1fix z range z -0.1 0.1hist n 5hist unbalset mech force 50solvesave beam-brace0.sav;prop coh 1e3 tens 1e3model null range x 2 4 y 2 6 z 5 10set largeini xdis 0 ydis 0 zdis 0sel beam begin=( 2, 4, 8) end=( 4, 4, 8) nseg=2sel beam prop emod=2.0e11 nu=0.30sel beam prop XCArea=6e-3 XCIz=200e-6 XCIy=200e-6 XCJ=0.0hist gp zdisp 4 4 8solvesave beam-brace1.sav;plot create GravVplot set plane dip 90 dd 0 origin 3 4 0plot set rot 15 0 20plot set center 2.5 4.2 4.0plot set cap size 25plot add cont disp plane behind shade onplot add sel beam force fxplot add sel geom black black node=off shrinkfac=0.0 plot add axesplot show 10-5;非全长锚固、预紧力锚杆（锚索）模拟;方法1、通过删除－建立link链接来模拟托盘gen zone radtun p0 0,0,0 p1 25,0,0 p2 0,50, 0 p3 0,0,25 size 4 25 4 10 dim 4 4 4 4 ratio 1 1 1 1.1 fill gen zone reflect normal 1 0 0 ori 0 0 0gen zone reflect normal 0 0 1 ori 0 0 0mo mohrpro bulk 2.2e9 she 1.3e9 fric 30 coh 1.3e6 ten 1.5e5ini dens 2000fix x range x -25.1,-24.9fix x range x 24.9 25.1fix y range y 49.9 50.1fix z range z -25.1 -24.9fix z range z 24.9 25.1sel cable id=1 beg 0, 0, 0 end 0 ,29, 0 nseg 10sel cable id=1 beg 0,29,0 end 0,35,0 nseg 6sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &gr_coh 1 gr_k 1 gr_per 0.0785 range cid 1,10sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &gr_coh 10e5 gr_k 2e7 range cid 11,17sel delete link range id 1sel link id=100 1 target zonesel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=rigid yrdir=rigid zrdir=rigid range id 100sel cable id=1 pretension 60e3 range cid 1,10step 2000sav 10-5.sav10-6;非全长锚固、预紧力锚杆（锚索）模拟;方法2、通过设置极大锚固剂参数模拟托盘gen zone radtun p0 0,0,0 p1 25,0,0 p2 0,50, 0 p3 0,0,25 size 4 25 4 10 dim 4 4 4 4 ratio 1 1 1 1.1 fillgen zone reflect normal 1 0 0 ori 0 0 0gen zone reflect normal 0 0 1 ori 0 0 0mo mohrpro bulk 2.2e9 she 1.3e9 fric 30 coh 1.3e6 ten 1.5e5ini dens 2000fix x range x -25.1,-24.9fix x range x 24.9 25.1fix y range y 49.9 50.1fix z range z -25.1 -24.9fix z range z 24.9 25.1sel cable id=1 beg 0, 0, 0 end 0 ,29, 0 nseg 10sel cable id=1 beg 0,29,0 end 0,35,0 nseg 6sel cable prop emod 2e10 ytension 310e3 xcarea 0.0004906 &gr_coh 1 gr_k 1 gr_per 0.0785 range cid 2,10sel cable prop emod 2e10 ytension 310e3 xcarea 0.0004906 &gr_coh 10e5 gr_k 2e7 range cid 11,17sel cable prop emod 2e10 ytension 310e3 xcarea 0.0004906 &gr_coh 10e8 gr_k 2e10 range cid 1,1sel cable id=1 pretension 60e3 range cid 1,10step 2000sav 10-6.sav10-7;非全长锚固、预紧力锚杆（锚索）模拟;方法3：借助别的结构单元（如liner单元）来模拟托盘gen zone radtun p0 0,0,0 p1 25,0,0 p2 0,50, 0 p3 0,0,25 size 4 25 4 10 dim 4 4 4 4 ratio 1 1 1 1.1 fillgen zone reflect normal 1 0 0 ori 0 0 0gen zone reflect normal 0 0 1 ori 0 0 0mo mohrpro bulk 2.2e9 she 1.3e9 fric 30 coh 1.3e6 ten 1.5e5ini dens 2000fix x range x -25.1,-24.9fix x range x 24.9 25.1fix y range y 49.9 50.1fix z range z -25.1 -24.9fix z range z 24.9 25.1sel cable id=1 beg 0, 0, 0 end 0 ,29, 0 nseg 10sel cable id=1 beg 0,29,0 end 0,35,0 nseg 6 sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &gr_coh 1 gr_k 1 gr_per 0.0785 range cid 1,10sel cable id=1 prop emod 2e10 ytension 310e3 xcarea 0.0004906 &gr_coh 10e5 gr_k 2e7 range cid 11,17sel liner range y=-.1, .1 x=-1,1 z=-1,1sel liner PROP iso=( 25e9, 0.15) thick=0.1 ; concrete sel liner PROP cs_nk=8e8 cs_sk=8e8 &cs_ncut=0.0 cs_scoh=0.0 cs_scohres=0.0 cs_sfric=0.0 sel delete link range id 1sel link id=100 1 target node tgt_num 18sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=rigid yrdir=rigid zrdir=rigid range id 100sel cable id=1 pretension 60e3 range cid 1,10step 2000sav 10-7.sav10-8ntitle Structure_dynamic_analysis_lakewaterconfig dynsel pile id=1 beg 0 0 0 end 0 0 1sel pile prop dens 2400 &Emod 1.0e10 Nu 0.3 XCArea 0.3 &XCJ 0.16375 XCIy 0.00625 XCIz 0.01575 &Per 2.8 &CS_sK 1.3e11 CS_sCoh 0.0 CS_sFric 10.0 &CS_nK 1.3e11 CS_nCoh 0.0 CS_nFric 0.0 CS_nGap offdef f1whilesteppingf0=10000*sin(10*dytime)np = nd_headloop while np # nullif nd_pos(np,1,3)=1nd_apply(np,1)=f0endifnp = nd_next(np)endloopendsel node fix x y z xr yr zr ran id=1sel set damp combinedplo cre pileplo current pileplo set back black fore white mag 0.8plo add sel geo id on sca .04 magf 1e3plo add sel fapp lgreen magf 1e3plo add sel pile mom my lblue lred magf 1e3 axe yel set movie avi step 100 file pile.avimovie startsol age 1movie finish11-1newconf dyn ;打开动力计算功能gen zone brick size 10 5 10mod elasmod null range x=0,5 z=5,10 ;删除部分网格fix z range x=-.1 .1 z=.1 10.1 ;设置静力边界条件fix z range x=9.9,10.1 z=.1 10.1fix y range y=-.1 .1fix y range y=4.9 5.1prop bulk 2e8 shear 1e8 ;设置土体参数prop bulk 4e9 shear 2e9 range x=5,6 z=5,10 ;设置墙体参数（土体参数的20倍）ini dens 2000 ;设置密度def setup ;动荷载中的变量赋值freq = 1.0omega = 2.0 * pi * freqold_time = clockendsetup ;执行变量赋值def wave ;定义动荷载函数wave = sin(omega * dytime) ;定义动荷载变量endapply xvel = 1 hist wave range z=-.1 .1 ;施加动荷载apply zvel = 0 range z=-.1 .1hist gp xvel 5,2,0hist gp xvel 5,2,10hist gp zvel 5,2,10hist dytimedef tim ;估算程序运行的时间tim = 0.01 * (clock - old_time)end set dyn multi on ;设置动态多步solve age 1.0print tim ;输出计算时间print dyn ;输出动力计算相关信息save mult1.sav11-2newconfig dyngen zone brick size 1,1,50model elasprop shear 1e7 bulk 2e7ini dens 1000def setupomega = 2.0 * pi * freqpulse = 1.0 / freqendset freq=4.0setupdef waveif dytime > pulsewave = 0.0elsewave = 0.5 * (1.0 - cos(omega * dytime))endifendrange name bottom z=-.1 .1fix z range z=.5 55 ;将上部网格都施加数值向约束apply dquiet squiet range bottomapply sxz -2e5 hist wave syz 0.0 szz 0.0 range bottom ;-2e5的系数来源于的值apply nvel 0 plane norm 0,0,1 range bottomhist gp xvel 0,0,0hist gp xvel 0,0,25hist gp xvel 0,0,50hist dytimehist waveplot create hhhplot add hist 1 2 3 vs 4plot showsolve age 211-3new;第一步：静力计算阶段config dynset dyn offgen zone brick size 6 3 2gen zone brick size 2 3 2 p0 0 0 2gen zone brick size 2 3 2 p0 4 0 2gen zone wedge size 1 3 2 p0 2 0 2gen zone wedge size 1 3 2 p0 4 3 2 p1 3 3 2 p2 4 0 2 p3 4 3 4 &p4 3 0 2 p5 4 0 4model elasticprop bulk 66667 shear 40000ini dens 0.0025set grav 0 0 -10fix x range x -0.01 0.01fix x range x 5.99 6.01fix y range y -0.01 0.01fix y range y 2.99 3.01fix z range z -0.1 0.1hist unbalsolvesave 11-3_1.sav;第二步：动力计算阶段set dyn ondef iniwaveper = 0.01endiniwavedef wavewave = 0.5 * (1.0 - cos (2*pi*dytime/per))endfree x y z ran z -0.1 0.1 ;去掉模型底部原有的静力条件apply nquiet squiet dquiet ran z -0.1 0.1 ;静态边界条件apply dstress 1.0 hist wave ran z -0.1 0.1 ;加动力荷载apply ff ;施加自由场边界条件group ff_cornergroup ff_side ran x 0 6group ff_side ran y 0 3group main_grid ran x 0 6 y 0 3set dyn time = 0 ;设置动力计算从0s开始hist reset ;清空已有的历史信息hist unbalhist dytime ; 主体网格hist gp xvel 2 1 0hist gp xvel 2 1 5.0 ; 柱体网格hist gp xvel -1 -1 0hist gp xvel -1 -1 5 ;平行于y方向的二维自由场网格hist gp xvel -1 0 0hist gp xvel -1 0 5.0;平行于x方向的二维自由场网格hist gp xvel 2 -1 0hist gp xvel 2 -1 5.0solve age 0.015save 11-3_2.sav11-4conf dygen zone brick size 3,3,3model elasprop bulk 1e8 shear 0.3e8ini dens 1000fix z range z -.1,.1set dyn=on, grav 0 0 -10, hist_rep=1hist gp zdisp 3.0,1.5,3.0hist dytimeplot create hhplot add his 1 vs 2save 11-4_1.sav;保存文件，在后续计算中需调用该文件cyc 15011-5;（1）质量分量和刚度分量共同作用rest 11-4_1.savset dyn damp rayleigh 1 22.8solve age=0.2titlevertical displacement versus time (mass & stiffness damping)plot showpause;（2）只有质量分量rest 11-4_1.sav。