高等土力学课程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）]）
ifresmaxdisp（['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
升级会员 