高等土力学课程CamClay.docx
- 文档编号:11362602
- 上传时间:2023-02-28
- 格式:DOCX
- 页数:19
- 大小:98.95KB
高等土力学课程CamClay.docx
《高等土力学课程CamClay.docx》由会员分享,可在线阅读,更多相关《高等土力学课程CamClay.docx(19页珍藏版)》请在冰豆网上搜索。
高等土力学课程CamClay
基于修正剑桥模型模拟理想三轴不排水试验
——两种积分算法的对比分析(CZQ-SpringGod)
1、修正剑桥模型
在塑性功中考虑体积塑性应变的影响,根据屈服面一致性原则,假定屈服函数对硬化参数的偏导为0,就获得了以理想三轴不排水试验为基础的修正剑桥模型屈服函数:
(1)
其中
,
,
,
,M为临界线斜率,
为前期固结压力。
硬化/软化法则:
(2)
式中
为体积塑性应变,v为比体积,
为正常固结线斜率,
为回弹线斜率。
由于不排水屈服面推导过程是基于硬化参数
偏导为0,也就是说不排水试验中硬化参数同体积塑性应变无关,屈服面不变化,而若引入硬化法则就同屈服面推导过程中的假定矛盾,因此计算时将模型处理为理想塑性模型。
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)修正剑桥模型的参数设定:
临界线斜率:
M=1.1
正常固结线斜率:
=0.17
回弹线斜率:
=0.034
初始比体积:
v0=2.12
前期固结压力:
=100KPa
剪切与体积模量的比值:
GK=0.46155
每个增量步体积模量的计算:
剪切模量G=GK×K
其中固结线方程为:
。
(2)计算结果:
不排水有效应力路径:
(a)显示算法(b)隐式算法
图1不排水有效应力路径
偏应力随轴向应变的变化:
(a)显示算法(b)隐式算法
图2偏应变随轴向应变的变化
孔隙水压力随轴向应变的变化:
(a)显示算法(b)隐式算法
图3孔压随轴向应变的变化
两种算法的每个增量步同屈服面的偏移程度:
(a)显式算法(b)隐式算法
图4每个增量屈服面的偏移程度
结论:
两种算法在计算理想塑性修正剑桥模型时,数值解能很好地同理论屈服面符合。
显示算法的误差是递增的,而隐式是收敛的。
理想塑性模型的分析结果表明,经过屈服面修正后的显示算法在精度上要高于隐式算法,可能同收敛参数的设定有关,不过两者都是精确的。
参考文献:
[1]S.W.Sloan.A.J.Abbo.D.Sheng.Refinedexplicitintegrationofelasoplasticmodelswithautomaticerrocontrol[J].EngineeringComputations.2001:
18,121-19
程序代码:
显式积分算法:
(ExplicitIntegrationAlgorithm)
%functioncamclayexp
%%Undrainedcondition(perfectplasticity)
%%initializationofparameter
ek=0.034;%回弹斜率
lam=0.17;%固结斜率
M=1.1;%临界线斜率
v0=2.12;%初始比体积
GK=0.46155;%剪切与体积模量的比值
pc=100;%初始固结压力
%%Preliminary
S=[pcpcpc000];
[Pst,deviS]=deviT(S);
[J2,J3,sJ2,q,lode]=invar(deviS);
E=[000000];
nstep=300;
de1=0.0004;
q1=0;
dEpvol=0;
devidEp=zeros(1,6);
forn=1:
nstep
pcre(n)=pc;
Sre(n,:
)=S;
pre(n)=Pst;qre(n)=q;q1re(n)=q1;
%%strainincrement
dE=[de1-de1/2-de1/2000];
v1=v0-lam*log(pc);%固结曲线
K=v1*Pst/ek;%体积模量
G=GK*K;%剪切模量
m=[111000];
De=K*m'*m+2*G*(eye(6)-m'*m/3);%弹性刚度矩阵
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;
%%incrementofstrain:
initialization
YF1=ydfun(sJ2,Pst,pc,M);
S1=S+ddS;
[Pst1,deviS1]=deviT(S1);
[J2p,J3p,sJ2p,qp,lodep]=invar(deviS1);
YF2=ydfun(sJ2p,Pst1,pc,M);
ifYF2<=0
loop=-1;
S=S1;%弹性加载,或卸载
else
ifYF2>0%塑性加载
ifYF1<0
alph=alphfun(S,ddS,pc,M);
alphc=YF2/(YF2-YF1);
end
ifYF1>=0
dfp0=2*Pst-pc;dfj0=6*sJ2/M^2;dfo0=0;
flow0=FlowPl(deviS,dfp0,dfj0,dfo0,J2,sJ2,lode);
pW=flow0*ddS';
ifpW>0
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
%%ErrorcontrolofCorrector
toler=0.001;iter=0;T=0;dT=1;dpv=0;
whileloop==1&T<1
%%firststepformodification(flowisafunction)
sig1=S;
k1=dpv;
[Pst11,deviS11]=deviT(sig1);
[J21,J31,sJ21,q1,lode1]=invar(deviS11);
%pc1=harden(pc,v1,lam,ek,k1);
pc1=pc;
dfp1=2*Pst11-pc1;dfj1=6*sJ21/M^2;dfo1=0;%重要的流动参数
flow1=FlowPl(deviS11,dfp1,dfj1,dfo1,J21,sJ21,lode1);
%dh1=dhard(pc1,v1,lam,ek);
dh1=0;%perfectplasticity(nohardening)
dA1=-Pst11*dh1*(2*Pst11-pc1);Dep1=De-De*flow1'*flow1*De/(dA1+flow1*De*flow1');
dlam1=max(flow1*sige'*dT/(dA1+flow1*De*flow1'),0);
dsig1=sige*dT-dlam1*(De*flow1')';
dk1=dlam1*(2*Pst11-pc1);
%塑性体积应变硬化
%%secondstepformodification
sig2=sig1+dsig1;
k2=k1+dk1;
[Pst12,deviS12]=deviT(sig2);
[J22,J32,sJ22,q2,lode2]=invar(deviS12);
%pc2=harden(pc,v1,lam,ek,k2);
pc2=pc;
dfp2=2*Pst12-pc2;dfj2=6*sJ22/M^2;dfo2=0;
flow2=FlowPl(deviS12,dfp2,dfj2,dfo2,J22,sJ22,lode2);
%dh2=dhard(pc2,v1,lam,ek);
dh2=0;
dA2=-Pst12*dh2*(2*Pst12-pc2);Dep2=De-De*flow2'*flow2*De/(dA2+flow2*De*flow2');
dlam2=max(flow2*sige'*dT/(dA2+flow2*De*flow2'),0);
dsig2=sige*dT-dlam2*(De*flow2')';
dk2=dlam2*(2*Pst12-pc2);
%%errorcontrol
Err=(dsig2-dsig1)/2;
sm=S+(dsig1+dsig2)/2;%modifiedstress
rer=getmax(Err);
rsm=getmax(sm);
R=max(10^(-10),rer/rsm);
Tp=0.8*(toler/R)^0.5;
ifR>toler
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;%塑性体积应变
%%stresscorrection
[S,dpv]=correctfun(Sc,pc,v1,lam,ek,M,dpvc,De,iter,0);%0为无硬化
%S=sm;
%dpv=dpvc;
%%
T=T+dT;
beta=min([Tp,2]);%必须输入数组做参数
dT=beta*dT;
dT=min([dT,1-T]);
end
%%recordofiteration
ps=Sc;
[px(iter+2),pds]=deviT(ps);
[aa,bb,cc,py(iter+2),dd]=invar(pds);
iter=iter+1;
ifiter>10
disp('toomuchiteration')
break
end
end
disp(['covergediteration:
',num2str(iter)])
%px=[];py=[];
%%nextincrement
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)=q^2/M^2+Pst*(Pst-pc);
%fre1(n)=q-M*sqrt(-p.*(p-pc));
end
隐式算法:
(ImplicitIntegrationAlgorithm)
%functioncamclayimp
%%Undrainedcondition(perfectplasticity)
%%initializationofparameter
ek=0.034;%回弹斜率
lam=0.17;%固结斜率
M=1.1;%临界线斜率
v0=2.12;%初始比体积
GK=0.46155;%剪切与体积模量的比值
pc=100;%初始固结压力
%%Preliminary
S=[pcpcpc000];
[Pst,deviS]=deviT(S);
[J2,J3,sJ2,q,lode]=invar(deviS);
E=[000000];
nstep=300;
de1=0.001;
q1=0;
forn=1:
nstep
pcre(n)=pc;
Sre(n,:
)=S;
pre(n)=Pst;qre(n)=q;q1re(n)=q1;
%q1-q
%%strainincrement
dE=[de1-de1/2-de1/2000];
%v1=v0-lam*log(pc);%固结曲线
K=v0*Pst/ek;%体积模量
G=GK*K;%剪切模量
[meanE,devidE]=deviT(dE);
dEvol=meanE*3;
%%Elasticpredictor
Pst1=Pst+K*dEvol;
fori=1:
6
deviS1(i)=deviS(i)+2*G*devidE(i);
end
[J2t,J3t,sJ2t,qt,lodet]=invar(deviS1);
pc1=pc;
q1=qt;
%%
Yieldf=qt^2/M^2+Pst1*(Pst1-pc1);
%%Plasticcorrector
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;
%%iterationofresidule
%Pst1=Pst;
whileYieldf>0
res
(1)=Pst1-Pst-K*dEvol+K*dpl*(2*Pst1-pc1);
%res
(2)=pc1-pc*exp(v1/(lam-ek)*dpl*(2*Pst1-pc1));%hardeningrule
res
(2)=qt^2/(M*(1+6*G*dpl/M^2))^2+Pst1*(Pst1-pc1);
%disp([num2str(iter),'interation',num2str(dpl)])
resmax=getmax(res);
disp(['范数',num2str(resmax)])
ifresmax disp(['Convergence: ',num2str(iter)]) break end ifiter>=10 disp('toomuchinteration') break end iter=iter+1; %%theJacobianMatrixofResiduleVector ntdm(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)=M^2*(1+6*G*dpl/M^2)^2*(2*Pst1-pc1); ntdm(2,2)=Pst1*(Pst1-pc1)*M^2*2*(1+6*G*dpl/M^2)*6*G/M^2; %ntdm(3,3)=-Pst1*M^2*(1+6*G*dpl/M^2)^2; BM=-res; AM=ntdm; dru=soluN(AM,BM); %dru=inv(AM)*BM' %%updatetheresidule Pst1=Pst1+dru (1); dpl=dpl+dru (2); %pc1=pc1+dru(3); %%recordofiteration px(iter+2)=Pst+K*dEvol-K*dpl*(2*Pst1-pc1); py1(iter+2)=qt/(1+6*G*dpl/M^2); dEpvol=dpl*(2*Pst1-pc1); end disp(['increment: ',num2str(n)]) px=[];py1=[]; %%nextincrement Pst=Pst1; q1=qt/(1+6*G*dpl/M^2); deviS=deviS1/(1+dpl*6*G/M^2); [J2,J3,sJ2,q,lode]=invar(deviS); %pc=pc1;%有固结过程 pc=pc;%无固结过程 S=backT(deviS,Pst); fre(n)=q1^2/M^2+Pst*(Pst-pc); end 子程序: functiony=ydfun(steff,p,pc,M)----屈服函数 %%yieldfunction y=3*steff^2/M^2+p*(p-pc); function[p,sd]=deviT(s)----求张量的偏量 p=(s (1)+s (2)+s(3))/3; fori=1: 3 sd(i)=s(i)-p; end fori=4: 6 sd(i)=s(i); end function[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)... )+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)*DEVIA(6)*DEVIA(6)-DEVIA (2)*DEVIA(5)*DEVIA(5)-DEVIA(3)*... DEVIA(4)*DEVIA(4); sJ2=sqrt(J2); q=ROOT3*sJ2; ifJ2==0 SINT3=0; else SINT3=-3.0*ROOT3*J3/(2.0*J2*sJ2); end if(SINT3>1) SINT3=1; end if(SINT3<-1) SINT3=-1; end lode=asin(SINT3)/3.0; functionuh=harden(pc,v1,lam,ek,dpv)-----硬化法则 uh=pc*exp(v1/(lam-ek)*dpv); functionflow=FlowPl(s,dfp,dfj2,dfo,varj2,steff,lode)-----流动法则 ifsteff==0 flow=[111000]; return end root3=sqrt(3); tant3=tan(3*lode); cos3=cos(3*lode); %dfp=%对P偏导 %dfj2=%对sqrt(J2)偏导 %dfo=%对洛德角偏导 c1=dfp; c2=dfj2-tant3/steff*dfo; c3=-root3*dfo/(2*cos3*steff*varj2); vn1=[1/31/31/3000]; fori=1: 6 vn2(i)=s(i)/(2*steff); end vn3 (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); fori=1: 6 flow(i)=c1*vn1(i)+c2*vn2(i)+c3*vn3(i); end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 高等 土力学 课程 CamClay