vol=dpl*(2*Pst1-pc1);enddispdmcrement:?m2str(n)Dpx=[];pyl=[];%%nextincrementPsl=Pst1;ql=qt/(l+6*G*dpl/MA2);deviS=deviS1/(1+dpl*6*G/MA2);[J2,J3,sJ2,q,lode]=invar(deviS);%pc=pcl;%有固结过程pc=pc;%无固结过程S=backT(deviS,Pst);fre(n)=q1A2/MA2+Pst*(Pst-pc);end子程序:functiony=ydfun(steff,p,pc,M)-…屈服函数%%yieldfunctiony=3*steffA2/MA2+p*(p-pc);function[p,sd]=deviT(s)-…求张量的偏量p=(s(l)+s⑵+s⑶)/3;fori=l:3sd(i)=s(i)-p;endfori=4:6sd(i)=s(i);endfunction[J2,J3,sJ2,q,lode]二invar(DEVIA)■—求偏量的不变量n=length(DEVIA);R00T3=1.73205080757;J2=0.5*(DEVIA(1)*DEVIA(1)+DEVIA(2)*DEVIA(2)+DEVIA(3)*DEVIA(3)...)+DEVIA(4)*DEVIA(4)+DEVIA(5)*DEVIA(5)+DEVIA(6)*DEVIA(6);J3=DEVIA(1)*DEVIA(2)*DEVIA(3)+2*DEVIA(4)*DEVIA(5)*DEVIA(6)-...DEVIA(1)*DEV1A(6)*DEVIA(6)-DEVIA(2)*DEVIA(5)*DEVIA(5)-DEVIA(3)*...DEV1A(4)*DEVIA(4);sJ2=sqrt(J2);q=ROOT3*sJ2;ifJ2==0