new;网格建立;;;;;;;;;;;;;;;;;;;;;;;;;gen zone brick p0 90 0 -30 p1 202 0 -30 p2 90 4 -30 p3 90 0 0 size 112 4 30 ratio 1 1 1gen zone brick p0 90 0 -30 p1 90 0 0 p2 90 4 -30 p3 0 0 -30 size 30 4 25 ratio 1 1 1.1gen zone brick p0 90 0 -30 p1 0 0 -30 p2 90 4 -30 p3 90 0 -75 size 25 4 18 ratio 1.1 1 1.1gen zone brick p0 90 0 -30 p1 90 0 -75 p2 90 4 -30 p3 202 0 -30 size 18 4 112 ratio 1.1 1 1 gen zone brick p0 202 0 -30 p1 292 0 -30 p2 202 4 -30 p3 202 0 0 size 25 4 30 ratio 1.1 1 1 gen zone brick p0 202 0 -30 p1 202 0 -75 p2 202 4 -30 p3 292 0 -30 size 18 4 25 ratio 1.1 1 1.1;分组;;;;;;;;;;;;;;;;;;;;;;;;;;group 1 range x 90 110 y 0 4 z -30 0group 1 range x 180 202 y 0 4 z -30 0group 2 range group 1 not;建立连续墙单元;;;;;;;;;;;;;;;;;;;;;;;;;;gen separate 1gen merge 1e-4 range x 90 110 y 0 4 z -30.1 -29.9gen merge 1e-4 range x 180 202 y 0 4 z -30.1 -29.9attach face range x 89.99 90.01 y 0.0 4.0 z -29.9 0attach face range x 109.99 110.01 y 0.0 4.0 z -29.9 0attach face range x 179.99 180.01 y 0.0 4.0 z -29.9 0attach face range x 201.99 202.01 y 0.0 4.0 z -29.9 0sel liner id 1 crossdiag group 2 range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1sel liner id 2 crossdiag group 2 range x 109.9 110.1 y -0.1 4.1 z -30.1 0.1sel liner id 3 crossdiag group 2 range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1sel liner id 4 crossdiag group 2 range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1sel liner id 1 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &cs_nk=4e9 cs_sk=4e9 &cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1sel liner id 2 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &cs_nk=4e9 cs_sk=4e9 &cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &range x 109.9 110.1 y -0.1 4.1 z -30.1 0.1sel liner id 3 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &cs_nk=4e9 cs_sk=4e9 &cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1sel liner id 4 prop isotropic (2.0e10, 0.20) thickness 0.8 density 2.5e3 &cs_nk=4e9 cs_sk=4e9 &cs_ncut=4e7 cs_scoh=4e7 cs_scohres=0 cs_sfric=20.0 &range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1;定义支撑结构;;;;;;;;;;;;;;;;;;;;;;;def struct_install1loop i(1,3)structx_zz=-1.0*5.0*(i-1)structx_xx0=90.0structx_xx1=110.0structx_yy=2.0commandsel beam id=2 begin (structx_xx0,structx_yy,structx_zz) end (structx_xx1,structx_yy,structx_zz) nseg=10sel beam id=2 prop dens=0.000 emod=1.0e-6 nu=0.2 &xcarea=0.80 xcj=10.94e-2 xciy=6.67e-2 xciz=4.27e-2 ydirection=(0 0 -1) ;1000x800endcommandendloopendstruct_install1def struct_install2loop i(1,3)structx_zz=-1.0*5.0*(i-1)structx_xx0=180.0structx_xx1=202.0structx_yy=2.0commandsel beam id=3 begin (structx_xx0,structx_yy,structx_zz) end (structx_xx1,structx_yy,structx_zz) nseg=11sel beam id=3 prop dens=0.000 emod=1.0e-6 nu=0.2 &xcarea=0.80 xcj=10.94e-2 xciy=6.67e-2 xciz=4.27e-2 ydirection=(0 0 -1) ;1000x800endcommandendloopendstruct_install2;建立结构单元分组;;;;;;;;;;;;;;;;;;;;;;;;;;;sel group linerwall range sel linersel group struct1 range sel beam x (90.0 110.0) z (-0.1 0.1)sel group struct2 range sel beam x (90.0 110.0) z (-5.1 -4.9)sel group struct3 range sel beam x (90.0 110.0) z (-10.1 -9.9)sel group struct4 range sel beam x (180.0 202.0) z (-0.1 0.1)sel group struct5 range sel beam x (180.0 202.0) z (-5.1 -4.9)sel group struct6 range sel beam x (180.0 202.0) z (-10.1 -9.9);删除beam单元的linksel dele link range sel beam z (-30 0);建立liner间的节点间的刚性linkdef merge_link0node_num=0node_pnt0 = nd_headloop while node_pnt0 # null ;寻找总节点数，注：不能自己任生成node，程序缺省的方式为连续生成无不连续node_num = node_num+1node_pnt0 = nd_next(node_pnt0)endloopnode_num_minus1 = node_num-1link_id=30000loop ii (1,node_num_minus1)node_pnt1 = nd_find(ii)xxa = nd_pos(node_pnt1,2,1)yya = nd_pos(node_pnt1,2,2)zza = nd_pos(node_pnt1,2,3)ii_plus1 = ii+1loop jj (ii_plus1,node_num)node_pnt2 = nd_find(jj)xxb = nd_pos(node_pnt2,2,1)yyb = nd_pos(node_pnt2,2,2)zzb = nd_pos(node_pnt2,2,3)node_dist = sqrt((xxa-xxb)^2+(yya-yyb)^2+(zza-zzb)^2)dist_tol = 1e-1if node_dist <= dist_tol thenlink_pnt1 = nd_link(node_pnt1)link_pnt2 = nd_link(node_pnt2);if link_pnt1 # null then; temp1 = lk_delete(link_pnt1);endifif link_pnt2 # null thentemp2 = lk_delete(link_pnt2)endiflink_id = link_id+1command ;生成新link(6自由度全固结),大的node的id作为target node,小的node的id作为source node，需注意不同情况下的灵活调整sel set link node_tol=dist_tolsel link id=link_id jj target = node tgt_num =ii ;指定link的ID;sel link ii target = node tgt_num = jj ;不指定link的id，自动生成sel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=rigid yrdir=rigid zrdir=rigid range id=link_idendcommandendifendloopendloopendmerge_link0;设置土层材料参数;;;;;;;;;;;;;;;;;;;;;;;;;;;;def b_s_modb_mod =e_mod/(3.0*(1.0-2.0*p_ratio))s_mod =e_mod/(2.0*(1.0+p_ratio))endmodel elasticset e_mod 100e6set p_ratio 0.3b_s_modprop bu=b_mod sh=s_modini dens 1800 range z -75 0def ini_szzszz0=0szzgrad=1800*10commandini szz add szz0 grad 0 0 szzgrad range z -75 0endcommandendini_szzdef ini_sxx_syypnt=zone_headloop while pnt # nullval=k0*z_szz(pnt)z_sxx(pnt)=valz_syy(pnt)=valpnt=z_next(pnt)endloopendset k0=0.50ini_sxx_syy;定义边界处的结构边界条件;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;cyc 0sel node local xdir=(0,1,0) ydir=(0,0,1) range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1 sel node local xdir=(0,1,0) ydir=(0,0,-1) range x 109.9 110.1 y -0.1 4.1 z -30.1 0.1 sel node local xdir=(0,1,0) ydir=(0,0,1) range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1 sel node local xdir=(0,1,0) ydir=(0,0,-1) range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1sel node fix lsys range x 89.9 90.1 y -0.1 0.1 z -30.1 0.1sel node fix lsys range x 89.9 90.1 y 3.9 4.1 z -30.1 0.1sel node fix lsys range x 109.9 110.1 y -0.1 0.1 z -30.1 0.1sel node fix lsys range x 109.9 110.1 y 3.9 4.1 z -30.1 0.1sel node fix lsys range x 179.9 180.1 y -0.1 0.1 z -30.1 0.1sel node fix lsys range x 179.9 180.1 y 3.9 4.1 z -30.1 0.1sel node fix lsys range x 201.9 202.1 y -0.1 0.1 z -30.1 0.1sel node fix lsys range x 201.9 202.1 y 3.9 4.1 z -30.1 0.1sel node fix x yr zr range x 89.9 90.1 y -0.1 0.1 z -30.1 0.1sel node fix x yr zr range x 89.9 90.1 y 3.9 4.1 z -30.1 0.1sel node fix x yr zr range x 109.9 110.1 y -0.1 0.1 z -30.1 0.1sel node fix x yr zr range x 109.9 110.1 y 3.9 4.1 z -30.1 0.1sel node fix x yr zr range x 179.9 180.1 y -0.1 0.1 z -30.1 0.1sel node fix x yr zr range x 179.9 180.1 y 3.9 4.1 z -30.1 0.1sel node fix x yr zr range x 201.9 202.1 y -0.1 0.1 z -30.1 0.1sel node fix x yr zr range x 201.9 202.1 y 3.9 4.1 z -30.1 0.1sel node fix y range x 89.9 90.1 y 0.0 4.0 z -0.1 0.1sel node fix y range x 109.9 110.1 y 0.0 4.0 z -0.1 0.1sel node fix y range x 179.9 180.1 y 0.0 4.0 z -0.1 0.1sel node fix y range x 201.9 202.1 y 0.0 4.0 z -0.1 0.1;set plot meta;plot set rot 20 0 30 ba wh color=on cent=(10 20 0) mag=3.81;set outp node_local_sys.wmf;plot add sel geom black red link=off node=off id=off shrink=0 scale=0.03 nodesys=on range group linerwall any group struct1 any;pl ha;固定边界条件;;;;;;;;;;;;;;;;;;;;;;;;;;fix x range x -0.1 0.1fix x range x 291.9 292.1fix y range y -0.1 0.1fix y range y 3.9 4.1fix x y z range z -75.1 -74.9set grav 0,0,-10solvesave elas.sav;删除侧面内外土体间的连接约束;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;attach delete range x 89.99 90.01 y 0.0 4.0 z -29.9 0attach delete range x 109.99 110.01 y 0.0 4.0 z -29.9 0attach delete range x 179.99 180.01 y 0.0 4.0 z -29.9 0attach delete range x 201.99 202.01 y 0.0 4.0 z -29.9 0;在墙内土体的外侧建立接触面;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;interface 1 face range group 1 x 89.99 90.01 y 0.0 4.0 z -29.9 0interface 2 face range group 1 x 109.99 110.01 y 0.0 4.0 z -29.9 0interface 3 face range group 1 x 179.99 180.01 y 0.0 4.0 z -29.9 0interface 4 face range group 1 x 201.99 202.01 y 0.0 4.0 z -29.9 0interface 1 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数interface 2 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数interface 3 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数interface 4 prop kn=4e8 ks=4e8 tens=5e3 coh=0.0 fric=20 ;接触面参数interface 1 maxedge=1interface 2 maxedge=1interface 3 maxedge=1interface 4 maxedge=1;interface 1 prop kn=4e8 ks=4e8 tens=1e10 sbratio=100;plot set ba wh;pl ske interface red blue attach cyan green;set outp interface_attachment.wmf;pl ha;重新定义连续墙参数;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;sel liner id 1 prop isotropic (2.0e10, 0.20) &cs_nk=4e9 cs_sk=4e9 &cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 89.9 90.1 y -0.1 4.1 z -30.1 0.1sel liner id 2 prop isotropic (2.0e10, 0.20) &cs_nk=4e9 cs_sk=4e9 &cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 109.9 110.1 y-0.1 4.1 z -30.1 0.1sel liner id 3 prop isotropic (2.0e10, 0.20) &cs_nk=4e9 cs_sk=4e9 &cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 179.9 180.1 y -0.1 4.1 z -30.1 0.1sel liner id 4 prop isotropic (2.0e10, 0.20) &cs_nk=4e9 cs_sk=4e9 &cs_scoh=4e7 cs_scohres=0.0 cs_sfric=0.0 range x 201.9 202.1 y -0.1 4.1 z -30.1 0.1;重新定义墙底约束条件;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;def redef_wall_end_link1node_pnt = nd_headlink_id=100000loop while node_pnt # nullnode_id = nd_id(node_pnt)xx = nd_pos(node_pnt,2,1)yy = nd_pos(node_pnt,2,2)zz = nd_pos(node_pnt,2,3)link_pnt = nd_link(node_pnt)dist_x = sqrt((xx-90.0)^2+(zz+30.0)^2)if dist_x <=dist_tol thenif link_pnt # null thentemp1 = lk_delete(link_pnt)\link_id = link_id+1commandsel set link node_tol = dist_tolsel link id=link_id node_id target zonesel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_idendcommandendifendifnode_pnt = nd_next(node_pnt)endloopendredef_wall_end_link1def redef_wall_end_link2node_pnt = nd_headlink_id=150000loop while node_pnt # nullnode_id = nd_id(node_pnt)xx = nd_pos(node_pnt,2,1)yy = nd_pos(node_pnt,2,2)zz = nd_pos(node_pnt,2,3)link_pnt = nd_link(node_pnt)dist_x = sqrt((xx-110.0)^2+(zz+30.0)^2)dist_tol = 1e-1if dist_x <=dist_tol thenif link_pnt # null thenif yy < 85.0temp1 = lk_delete(link_pnt)\link_id = link_id+1commandsel set link node_tol = dist_tolsel link id=link_id node_id target zonesel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_idendcommandendifendifendifnode_pnt = nd_next(node_pnt)endloopendredef_wall_end_link2def redef_wall_end_link3node_pnt = nd_headlink_id=200000loop while node_pnt # nullnode_id = nd_id(node_pnt)xx = nd_pos(node_pnt,2,1)yy = nd_pos(node_pnt,2,2)zz = nd_pos(node_pnt,2,3)link_pnt = nd_link(node_pnt)dist_x = sqrt((xx-180.0)^2+(zz+30.0)^2)dist_tol = 1e-1if dist_x <=dist_tol thenif link_pnt # null thentemp1 = lk_delete(link_pnt)\link_id = link_id+1commandsel set link node_tol = dist_tolsel link id=link_id node_id target zonesel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_idendcommandendifendifnode_pnt = nd_next(node_pnt)endloopendredef_wall_end_link3def redef_wall_end_link4node_pnt = nd_headlink_id=250000loop while node_pnt # nullnode_id = nd_id(node_pnt)xx = nd_pos(node_pnt,2,1)yy = nd_pos(node_pnt,2,2)zz = nd_pos(node_pnt,2,3)link_pnt = nd_link(node_pnt)dist_x = sqrt((xx-202.0)^2+(zz+30.0)^2)dist_tol = 1e-1if dist_x <=dist_tol thenif link_pnt # null thentemp1 = lk_delete(link_pnt)\link_id = link_id+1commandsel set link node_tol = dist_tolsel link id=link_id node_id target zonesel link attach xdir=rigid ydir=rigid zdir=rigid xrdir=free yrdir=free zrdir=free range id=link_idendcommandendifendifnode_pnt = nd_next(node_pnt)endloopendredef_wall_end_link4;剑桥模型;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;model cam-clay;cam-clay模型则不需定义弹性模量(E、G、K)等参数，自动计算;cam-clay模型中需确定8个模型参数(①-⑧)，手册property中的初始比体积cv(v0)和shear 无须给定def install_proppnt=zone_headloop while pnt # nullabs_sxx=abs(z_sxx(pnt)) ;|sxx|abs_syy=abs(z_syy(pnt)) ;|syy|abs_szz=abs(z_szz(pnt)) ;|szz|p0=(abs_sxx+abs_syy+abs_szz)/3.0;cam-clay模型中p、q均须为正值，p0由初应力场确定，故cam-clam定义模型参数前须先已知初应力p0_effective=p0-z_pp(pnt) ;p0';q0=sqrt(((abs_sxx-abs_syy)^2+(abs_syy-abs_szz)^2+(abs_szz-abs_sxx)^2)*0.5)q0=sqrt(((abs_sxx-abs_syy)^2+(abs_syy-abs_szz)^2+(abs_szz-abs_sxx)^2)*0.5+3.0*((z_sxy(pnt)) ^2+(z_sxz(pnt))^2+(z_syz(pnt))^2))z_prop(pnt,'mm')=6.0*sin(fai*degrad)/(3.0-sin(fai*degrad)) ;①注三角函数中需将角度转化为弧度temp1=q0/(z_prop(pnt,'mm')*p0_effective)pc0=p0_effective*(1.0+temp1^2)*OCR ;先期有效固结压力，用于确定屈服面v0=1.0+_e0z_prop(pnt,'cam_cp')=p0_effective ;★重要参数，否则不能正确计算有效应力，提示出错"Mean effective pressure is negative"z_prop(pnt,'mpc')=pc0 ;②z_prop(pnt,'poisson')=p_ratio ;③z_prop(pnt,'lambda')=_lambda ;④z_prop(pnt,'kappa')=_kappa ;⑤z_prop(pnt,'mp1')=_mp1 ;⑥z_prop(pnt,'mv_l')=v0+_lambda*ln(2.0*_cu/(z_prop(pnt,'mm')*_mp1))+(_lambda-_kappa)*l n(2.0) ;⑦z_prop(pnt,'bulk_bound')=100*40e6 ;⑧;z_prop(pnt,'bulk_bound')=100*(s_mod+4.0/3.0*s_mod) ;弹性体模上界Kmax;自动确定Kmax时会出现“property bad”错误提示;因为弹性上界对计算结果无影响，在不提示Kmax太小的性况下，取值越小计算收敛越快pnt=z_next(pnt)endloopendset p_ratio=0.25 fai=34.5 _lambda=0.14 _kappa=0.012 _mp1=1e3 _e0=1.2 _cu=10e3 OCR=1.0 ;模型所需参数install_propsolvesave model_cam.sav。