pfc三轴压缩

newSETrandom;随机生成plowallploaddball;----------------------------------------------------defmake_walls;生成墙wextend=0.1hextend=0.1w_stiff=1.6e10;墙的刚度，取球的1.1倍_width=-width_x0=_width*(1.0+wextend);width是长方体的宽度的一半_y0=_width*(1.0+wextend)_z0=0_x1=width*(1.0+wextend)_y1=_width*(1.0+wextend)_z1=0_x2=width*(1.0+wextend)_y2=width*(1.0+wextend)_z2=0_x3=_width*(1.0+wextend)_y3=width*(1.0+wextend)_z3=0commandwallid=1kn=w_stiffface(_x0,_y0,_z0)(_x1,_y1,_z1)(_x2,_y2,_z2)&(_x3,_y3,_z3);要按右手法则建立,而且pfc里面的坐标轴是一定的，要按照里面的来end_command_x0=_width*(1.0+wextend)_y0=_width*(1.0+wextend)_z0=height;height是长方体的高度_x1=_width*(1.0+wextend)_y1=width*(1.0+wextend)_z1=height_x2=width*(1.0+wextend)_y2=width*(1.0+wextend)_z2=height_x3=width*(1.0+wextend)_y3=_width*(1.0+wextend)_z3=heightcommandwallid=2kn=w_stiffface(_x0,_y0,_z0)(_x1,_y1,_z1)(_x2,_y2,_z2)&(_x3,_y3,_z3)end_command_x0=width_y0=_width*(1.0+wextend)_z0=height*(1.0+hextend)_x1=width_y1=width*(1.0+wextend)_z1=height*(1.0+hextend)_x2=width_y2=width*(1.0+wextend)_z2=-hextend*height_x3=width_y3=_width*(1.0+wextend)_z3=-hextend*heightcommandwallid=3kn=w_stiffface(_x0,_y0,_z0)(_x1,_y1,_z1)(_x2,_y2,_z2)&(_x3,_y3,_z3)end_command_x0=_width_y0=_width*(1.0+wextend)_z0=height*(1.0+hextend)_x1=_width_y1=_width*(1.0+wextend)_z1=-hextend*height_x2=_width_y2=width*(1.0+wextend)_z2=-hextend*height_x3=_width_y3=width*(1.0+wextend)_z3=height*(1.0+hextend)commandwallid=4kn=w_stiffface(_x0,_y0,_z0)(_x1,_y1,_z1)(_x2,_y2,_z2)&(_x3,_y3,_z3)end_command_x0=_width*(1.0+wextend)_y0=width_z0=-hextend*height_x1=width*(1.0+wextend)_y1=width_z1=-hextend*height_x2=width*(1.0+wextend)_y2=width_z2=height*(1.0+hextend)_x3=_width*(1.0+wextend)_y3=width_z3=height*(1.0+hextend)commandwallid=5kn=w_stiffface(_x0,_y0,_z0)(_x1,_y1,_z1)(_x2,_y2,_z2)&(_x3,_y3,_z3)end_command_x0=_width*(1.0+wextend)_y0=_width_z0=height*(1.0+hextend)_x1=width*(1.0+wextend)_y1=_width_z1=height*(1.0+hextend)_x2=width*(1.0+wextend)_y2=_width_z2=-hextend*height_x3=_width*(1.0+wextend)_y3=_width_z3=-hextend*heightcommandwallid=6kn=w_stiffface(_x0,_y0,_z0)(_x1,_y1,_z1)(_x2,_y2,_z2)&(_x3,_y3,_z3)end_commandend;----------------------------------------------------defassemble;颗粒生成s_stiff=0.0;初始切向刚度n_stiff=1.4e10;初始法向刚度tot_vol=4*height*width^2.0;长方体的体积rbar=0.5*(rlo+rhi);rlo最小球半径、rhi最大球半径num=int((1.0-poros)*tot_vol/(4.0/3.0*pi*rbar^3));poros越大，球数量越少mult=1.5;初始半径放大系数，越大，球数量越多rlo_0=rlo/multrhi_0=rhi/multcommandgenid=1,numrad=rlo_0,rhi_0x=_width,widthy=_width,widthz=0,heighttries10000000;方形的话就不用过滤器了;上面球的半径是缩小了mult倍的propdens=2700ks=s_stiffkn=n_stiff;dens即是岩石颗粒的密度end_commandii=out(string(num)+'particleswerecreated')sum=0.0;获得有效孔隙度bp=ball_head;是所有球里面的首地址loopwhilebp#null;当bp不是null的时候sum=sum+4.0/3.0*pi*b_rad(bp)^3;sum是前面所有的球的体积的总和bp=b_next(bp);下一个球的地址end_looppmeas=1.0-sum/tot_vol;1-球的总体积/圆柱体的体积mult=((1.0-poros)/(1.0-pmeas))^(1.0/3.0);将总的孔隙率调成1-poroscommandiniradmulmult;将球的半径放大mult倍cycle30000;执行1000个时步的循环计算propks=1.4e10fric0.5cycle3100;执行250个时步的循环计算end_commandend;----------------------------------------------------;----------------------------------------------------macrozero'inixvel0yvel0zvel0xspin0yspin0zspin0';宏，现在定义了后面就可以直接引用了xspin可能是x方向的角速度SETheight=0.1width=0.025rlo=0.001rhi=0.002poros=0.385make_wallsassemblezerosavett_ass.SAV;fname:triax_2.DATServo-controlandinitialstressstate-triaxsamplerestt_ass.SAV;restorecompactedassembly;----------------------------------------------------defget_ss;确定墙的平均应力和应变;之前的cycle可能产生了一定的位移，下面要除掉zdif=w_z(wadd2)-w_z(wadd1);w_z(wadd1)指的是墙1在之前产生的位移ydif=w_y(wadd5)-w_y(wadd6);所有的都是正的减去负的xdif=w_x(wadd3)-w_x(wadd4)new_width1=2*width+xdif;x方向新的宽度，这个时候就不是宽度的一半了，是整个宽度new_width2=2*width+ydif;y方向新的宽度new_height=height+zdif;除掉已经产生的位移，是新的高度wsxx=0.5*(w_xfob(wadd4)-w_xfob(wadd3))/(new_width2*new_height);x方向的应力wsyy=0.5*(w_yfob(wadd6)-w_yfob(wadd5))/(new_width1*new_height);ywszz=0.5*(w_zfob(wadd1)-w_zfob(wadd2))/(new_width1*new_width2);zwexx=2.0*xdif/(2*width+new_width1);x向应变参看手册186页3.16weyy=2.0*ydif/(2*width+new_width2);y向应变wezz=2.0*zdif/(height+new_height);z向应变wevol=wexx+weyy+wezz;体积应变end;----------------------------------------------------defget_gain;选取轴向和径向获得的伺服参数alpha=0.5;松弛参数count=0avg_stiff=0cp=contact_head;找到上、下墙上接触的平均数量loopwhilecp#nullifc_gobj2(cp)=wadd1count=count+1avg_stiff=avg_stiff+c_kn(cp)end_ififc_gobj2(cp)=wadd2count=count+1avg_stiff=avg_stiff+c_kn(cp)end_ifcp=c_next(cp)end_loopncount=count/2.0avg_stiff=avg_stiff/countgz=4*alpha*width^2.0/(avg_stiff*ncount*tdel);可以参看187页公式3.22，是增益参数count=0avg_stiff=0cp=contact_head;找到左、右墙上接触的平均数量loopwhilecp#nullifc_gobj2(cp)=wadd3count=count+1avg_stiff=avg_stiff+c_kn(cp)end_ififc_gobj2(cp)=wadd4count=count+1avg_stiff=avg_stiff+c_kn(cp)end_ifcp=c_next(cp)end_loopncount=count/2.0avg_stiff=avg_stiff/countgx=2*alpha*width*height/(avg_stiff*ncount*tdel);可以参看187页公式3.22，是增益参数count=0avg_stiff=0cp=contact_head;找到前、后墙上接触的平均数量loopwhilecp#nullifc_gobj2(cp)=wadd5count=count+1avg_stiff=avg_stiff+c_kn(cp)end_ififc_gobj2(cp)=wadd6count=count+1avg_stiff=avg_stiff+c_kn(cp)end_ifcp=c_next(cp)end_loopncount=count/2.0avg_stif