高等土力学课程CamClay.docx

1、高等土力学课程CamClay基于修正剑桥模型模拟理想三轴不排水试验 两种积分算法的对比分析（CZQ-SpringGod）1、修正剑桥模型 在塑性功中考虑体积塑性应变的影响，根据屈服面一致性原则，假定屈服函数对硬化参数的偏导为0，就获得了以理想三轴不排水试验为基础的修正剑桥模型屈服函数： （1）其中，M为临界线斜率，为前期固结压力。硬化/软化法则： （2） 式中为体积塑性应变，v为比体积，为正常固结线斜率，为回弹线斜率。由于不排水屈服面推导过程是基于硬化参数偏导为0，也就是说不排水试验中硬化参数同体积塑性应变无关，屈服面不变化，而若引入硬化法则就同屈服面推导过程中的假定矛盾，因此计算时将模型处理

2、为理想塑性模型。2、显式和隐式两种积分格式 考虑应变增量驱动下，第n增量步到第n+1增量步之间的应力积分格式。显式积分格式的推导参考文献1，其中弹塑性矩阵中的塑性硬化模量H=0。 隐式积分格式推导如下： （3） （4） （5） （6） （7）在这一组方程中没有硬化规律方程表明为理想塑性，并将式（3）-（7）合并化简得到： （8）式中求解（8）式方程组即可得到n+1增量步的各个增量。两种积分格式的matlab程序分别见邮件附件文件夹camclayexp和camclayimp，显式运行主程序为camclayexp.m，而隐式运行主程序为camclayimp.m。3、数值分析 （1）修正剑桥模型的参

3、数设定： 临界线斜率：M=1.1 正常固结线斜率： =0.17回弹线斜率： =0.034初始比体积：v0 =2.12前期固结压力： =100 KPa剪切与体积模量的比值：GK=0.46155每个增量步体积模量的计算： 剪切模量G=GKK 其中固结线方程为：。（2）计算结果： 不排水有效应力路径： （a）显示算法 (b) 隐式算法 图1 不排水有效应力路径 偏应力随轴向应变的变化： （a）显示算法 （b）隐式算法 图2 偏应变随轴向应变的变化孔隙水压力随轴向应变的变化： （a）显示算法 （b）隐式算法 图3 孔压随轴向应变的变化两种算法的每个增量步同屈服面的偏移程度： （a）显式算法 （b）隐式

4、算法 图4 每个增量屈服面的偏移程度结论：两种算法在计算理想塑性修正剑桥模型时，数值解能很好地同理论屈服面符合。显示算法的误差是递增的，而隐式是收敛的。理想塑性模型的分析结果表明，经过屈服面修正后的显示算法在精度上要高于隐式算法，可能同收敛参数的设定有关，不过两者都是精确的。参考文献：1 S.W.Sloan. A.J.Abbo. D.Sheng. Refined explicit integration of elasoplastic models with automatic erro controlJ. Engineering Computations. 2001:18,121-19程序代

5、码：显式积分算法：(Explicit Integration Algorithm)% function camclayexp% Undrained condition(perfect plasticity)% initialization of parameterek=0.034; % 回弹斜率lam=0.17; % 固结斜率M=1.1; % 临界线斜率v0=2.12; % 初始比体积GK=0.46155; % 剪切与体积模量的比值pc=100; % 初始固结压力% PreliminaryS=pc pc pc 0 0 0;Pst,deviS=deviT(S);J2,J3,sJ2,q,lode=

6、invar(deviS);E=0 0 0 0 0 0;nstep=300;de1=0.0004;q1=0;dEpvol=0;devidEp=zeros(1,6);for n=1:nstep pcre(n)=pc; Sre(n,:)=S; pre(n)=Pst;qre(n)=q;q1re(n)=q1; % strain increment dE=de1 -de1/2 -de1/2 0 0 0; v1=v0-lam*log(pc); % 固结曲线 K=v1*Pst/ek; % 体积模量 G=GK*K; % 剪切模量 m=1 1 1 0 0 0; De=K*m*m+2*G*(eye(6)-m*m/3

7、); % 弹性刚度矩阵 dre(n)=det(De);var(n)=q/Pst; meanE,devidE=deviT(dE); dEvol=meanE*3; ddS=(De*dE); % 弹性应力增量 pc=harden(pc,v1,lam,ek,0); % px(1)=Pst;py(1)=q; % increment of strain: initialization YF1=ydfun(sJ2,Pst,pc,M); S1=S+ddS; Pst1,deviS1=deviT(S1); J2p,J3p,sJ2p,qp,lodep=invar(deviS1); YF2=ydfun(sJ2p,Ps

8、t1,pc,M); if YF20 %塑性加载 if YF1=0 dfp0=2*Pst-pc;dfj0=6*sJ2/M2;dfo0=0; flow0=FlowPl(deviS,dfp0,dfj0,dfo0,J2,sJ2,lode); pW=flow0*ddS; if pW0 alph=0; else disp(弹性卸载，又加载) St=S+0.2*ddS; alph=alphfun(St,ddS,pc,M); end end loop=1; S=S+alph*ddS;alphre(n)=alph; sige=(1-alph)*ddS; % 找到屈服之后的弹性预测 end end % Error

9、 control of Corrector toler=0.001;iter=0;T=0;dT=1;dpv=0; while loop=1 & Ttoler beta=max(Tp,0.1); dT=beta*dT; else Sc=sm;lare(n)=(dlam1+dlam2)/2;dAre(n)=(dA1+dA2)/2;dpre(n)=(det(Dep1)+det(Dep2)/2; dpvc=dpv+(dk1+dk2)/2; % 塑性体积应变 % stress correction S,dpv=correctfun(Sc,pc,v1,lam,ek,M,dpvc,De,iter,0); %

10、 0 为无硬化% S=sm;% dpv=dpvc; % T=T+dT; beta=min(Tp,2); %必须输入数组做参数 dT=beta*dT; dT=min(dT,1-T); end % record of iteration ps=Sc; px(iter+2),pds=deviT(ps); aa,bb,cc,py(iter+2),dd=invar(pds); iter=iter+1; if iter10 disp(too much iteration) break end end disp(coverged iteration: ,num2str(iter)% px=;py=; % n

11、ext increment disp(increment: ,num2str(n) Pst,deviS=deviT(S); J2,J3,sJ2,q,lode=invar(deviS); dvpc(n)=dpv;% pc=harden(pc,v1,lam,ek,dpv); % 含有固结过程 fre(n)=q2/M2+Pst*(Pst-pc); %fre1(n)=q-M*sqrt(-p.*(p-pc);end隐式算法：（Implicit Integration Algorithm）% function camclayimp% Undrained condition(perfect plastici

12、ty)% initialization of parameterek=0.034; % 回弹斜率lam=0.17; % 固结斜率M=1.1; % 临界线斜率v0=2.12; % 初始比体积GK=0.46155; % 剪切与体积模量的比值pc=100; % 初始固结压力% PreliminaryS=pc pc pc 0 0 0;Pst,deviS=deviT(S);J2,J3,sJ2,q,lode=invar(deviS);E=0 0 0 0 0 0;nstep=300;de1=0.001;q1=0;for n=1:nstep pcre(n)=pc; Sre(n,:)=S; pre(n)=Pst

13、;qre(n)=q;q1re(n)=q1;% q1-q % strain increment dE=de1 -de1/2 -de1/2 0 0 0;% v1=v0-lam*log(pc); % 固结曲线 K=v0*Pst/ek; % 体积模量 G=GK*K; % 剪切模量 meanE,devidE=deviT(dE); dEvol=meanE*3; % Elastic predictor Pst1=Pst+K*dEvol; for i=1:6 deviS1(i)=deviS(i)+2*G*devidE(i); end J2t,J3t,sJ2t,qt,lodet=invar(deviS1); p

14、c1=pc; q1=qt; % Yieldf=qt2/M2+Pst1*(Pst1-pc1); % Plastic corrector iter=0; toler=1e-3; dEpvol=0; devidEp=zeros(1,6); dpl=0; % record px(2)=Pst1;px(1)=Pst; py1(2)=q1;py1(1)=q; % py2(2)=q1;py2(1)=q; % iteration of residule% Pst1=Pst; while Yieldf0 res(1)=Pst1-Pst-K*dEvol+K*dpl*(2*Pst1-pc1);% res(2)=pc

15、1-pc*exp(v1/(lam-ek)*dpl*(2*Pst1-pc1); % hardening rule res(2)=qt2/(M*(1+6*G*dpl/M2)2+Pst1*(Pst1-pc1); % disp(num2str(iter), interation ,num2str(dpl) resmax=getmax(res); disp(范数,num2str(resmax) if resmax=10 disp(too much interation) break end iter=iter+1; % the Jacobian Matrix of Residule Vector ntd

16、m(1,1)=1+2*K*dpl; ntdm(1,2)=K*(2*Pst1-pc1);% ntdm(1,3)=K*dpl;% ntdm(2,1)=-pc*exp(v1/(lam-ek)*dpl*(2*Pst1-pc1)*2*dpl*v1/(lam-ek);% ntdm(2,2)=-pc*exp(v1/(lam-ek)*dpl*(2*Pst1-pc1)*(2*Pst1-pc1);% ntdm(2,3)=1-pc*exp(v1/(lam-ek)*dpl*(2*Pst1-pc1)*(-v1/(lam-ek)*dpl); ntdm(2,1)=M2*(1+6*G*dpl/M2)2*(2*Pst1-pc1

17、); ntdm(2,2)=Pst1*(Pst1-pc1)*M2*2*(1+6*G*dpl/M2)*6*G/M2;% ntdm(3,3)=-Pst1*M2*(1+6*G*dpl/M2)2; BM=-res; AM=ntdm; dru=soluN(AM,BM); % dru=inv(AM)*BM % update the residule Pst1=Pst1+dru(1); dpl=dpl+dru(2);% pc1=pc1+dru(3); % record of iteration px(iter+2)=Pst+K*dEvol-K*dpl*(2*Pst1-pc1); py1(iter+2)=qt/

18、(1+6*G*dpl/M2); dEpvol=dpl*(2*Pst1-pc1); end disp(increment: ,num2str(n) px=;py1=; % next increment Pst=Pst1; q1=qt/(1+6*G*dpl/M2); deviS=deviS1/(1+dpl*6*G/M2); J2,J3,sJ2,q,lode=invar(deviS);% pc=pc1; % 有固结过程 pc=pc; % 无固结过程 S=backT(deviS,Pst); fre(n)=q12/M2+Pst*(Pst-pc); end子程序：function y=ydfun(stef

19、f,p,pc,M)-屈服函数% yield functiony=3*steff2/M2+p*(p-pc);function p,sd=deviT(s)-求张量的偏量p=(s(1)+s(2)+s(3)/3;for i=1:3 sd(i)=s(i)-p;endfor i=4:6 sd(i)=s(i);endfunction J2,J3,sJ2,q,lode=invar(DEVIA)-求偏量的不变量n=length(DEVIA);ROOT3=1.73205080757;J2=0.5*(DEVIA(1)*DEVIA(1)+DEVIA(2)*DEVIA(2)+DEVIA(3)*DEVIA(3). )+D

20、EVIA(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)*DEVIA(6)*DEVIA(6)-DEVIA(2)*DEVIA(5)*DEVIA(5)-DEVIA(3)*.DEVIA(4)*DEVIA(4);sJ2=sqrt(J2);q=ROOT3*sJ2;if J2=0 SINT3=0;elseSINT3=-3.0*ROOT3*J3/(2.0*J2*sJ2);endif (SINT31) SINT3=1 ;end

21、 if(SINT3-1) SINT3=-1 ;endlode=asin(SINT3)/3.0;function uh=harden(pc,v1,lam,ek,dpv)-硬化法则uh=pc*exp(v1/(lam-ek)*dpv);function flow=FlowPl(s,dfp,dfj2,dfo,varj2,steff,lode)-流动法则if steff=0 flow=1 1 1 0 0 0;return endroot3=sqrt(3);tant3=tan(3*lode);cos3=cos(3*lode);%dfp= % 对P偏导%dfj2= % 对sqrt(J2)偏导%dfo= %

22、对洛德角偏导c1=dfp;c2=dfj2-tant3/steff*dfo;c3=-root3*dfo/(2*cos3*steff*varj2);vn1=1/3 1/3 1/3 0 0 0;for i=1:6 vn2(i)=s(i)/(2*steff);endvn3(1)=s(2)*s(3)-s(6)2+varj2/3.0;vn3(2)=s(3)*s(1)-s(5)2+varj2/3.0;vn3(3)=s(1)*s(2)-s(4)2+varj2/3.0;vn3(4)=s(6)*s(5)-s(3)*s(4);vn3(5)=s(5)*s(4)-s(1)*s(6);vn3(6)=s(4)*s(6)-s(2)*s(5); for i=1:6 flow(i)=c1*vn1(i)+c2*vn2(i)+c3*vn3(i); end