下载文档原格式
第一章绪论一、土力学的研究对象土土体土：天然的地质材料。
岩石：经过风化、搬运/迁移、沉积变成了土。
土是第四纪沉积物，由岩石碎块、矿物颗粒、粘土矿物组成的松散集合体。
土的基本性质：非均质，不连续，各相异性，抗拉强度低，（tension weak）松散性，孔隙性，多相性，在渗流压力下的破碎性，力学压缩性，渗透性。
土力学的研究内容：1、土的工程特性。
2、土工建筑物的变形固结和稳定性。
学科特点：综合性强、经验性强、地区性强（区域土、特殊土）。
土质学是从地质学的角度出发研究土的组成成分、成因、变形机理、强度及其相互关系，并以求能进一步改善土质。
土力学是从工程力学的角度，通过实验来建立物理方程和分析工程特性，即，由控制方程得到土体的应力分布、变形及稳定性。
土力学发展简史沈珠江先生指出现代土力学应该由一个模型、三个理论和四个分支组成，一个模型是指土的本构模型；三个理论是指非饱和土固结理论、液化破坏理论和逐渐破坏理论；四个分支是指理论土力学、计算土力学、试验土力学和应用土力学。
液化破坏理论：动态液化、静态液化、稳定状态稳态强度。
二、土的变形与强度特性1、一般连续介质材料的变形特征（1）、弹性线性弹性、非线性弹性，所谓弹性就是说卸载后没有残余变形，加卸载都是同一路径即沿原曲线回到原点。
弹性的特点：①、加卸载同径，无残余变形 ②、应力应变一一对应③、线弹性时叠加原理成立 ④、与应力路径及应力历史无关σ=E ε；τ=G τ；γ=E/2(1+μ)。
σij p (平面应力) εV (体积应变) εijq (广义剪应力)γ(剪切应变)由上图知：对于弹性材料，剪应力与体积应变无关，而正应力与剪切应变也无关；即平面应力p 于广义剪应变γ无关，广义剪应力q 与体积应变εV 无关。
三向应力状态下的广义胡克定律为：εX = ［σX — γ ( σY +σZ )］/E γxy = τXY /G 体积变形模量（Bulk Modulus ）：m v vpK σεε==, 3m v m K K σεε==。
05年一、一粘土试样在三轴仪中，施加有效应力300kPa 下等向固结，固结完成后等向卸载至有效应力50kPa （状态B ）。
在此基础上进行常规排水压缩试验（侧向应力保持不变），使q ′＝100kPa （状态C ）。
然后试样在不排水条件下 加载至破坏（状态D ）。
利用修正剑桥模型分析上述试验。
土的材料参数为：λ＝0.16，κ＝0.04，Γ＝3.0 和M ＝1.0。
求 （1） 在e~lnp 坐标下绘出正常固结和临界状态线和A →B →C →D 路径； （2） 估计在A →B →C →D 试验过程中q 的峰值； （3） 估计土样破坏时q 的极限值和孔隙比e （状态D ）。
（1） 已知：λ＝0.16，κ＝0.04，Γ＝3.0 和M ＝1.0。
N 与Γ之间的关系为：08.3ln2)04.016.0(0.32ln )(=-+=-+Γ=k N λ 正常固结线（NCL ）为： p p N e '-='--=ln 16.008.2ln )1(λ； 临界状态线（CSL ）为： p p e '-='--Γ=ln 16.00.2ln )1(λ ； 卸载回弹线（κ线）为：p p p k p N e '-='--='-'---=ln 04.040.1ln 04.068.008.2 ln ln )()1(0κλ修正剑桥模型为：0/222='+''-'M q p p p c，则 030022='+'-'q p p据此可以绘出试样的应力路径（见图1）。
（2）B 点，650/300/=='cc p p 为重超固结土； BC 线的斜率为3.0, 方程为q=3.0(p-50) 所以C 点的3.83='cp kPa ； 代入030022='+'-'q p p ，可解得A →B →C →D 试验过程中q 的峰值为4.134'='cq kPa（3）土体破坏时，见图1中D 点D 点的e 值与C 点的相同，C 点的P=83.3kPa ，代入p e ln 04.040.1-= 求得e=1.22，同理，将e=1.22代入p e ln 16.0000.2-= 求得q=130.3 kPa解法2解：（1）已知：λ=0.16，κ=0.04，Γ＝3.0和M＝1.0，p’0=300kPa由N＝Γ＋(λ-κ) ln2得N＝Γ＋(λ-κ) ln2=3.0+（0.16-0.04）ln2≈3.083∴正常固结线(NCL)方程：v=1+e=N-λln p’，即e= N-1-λln p’= 3.083-1-0.16ln p’ =2.083-0.16 ln p’；临界状态线（CSL）方程：v=1+e=Γ-λln p’，即e=Γ-1-λln p’=3.0-1-0.16ln p’=2.0-0.16ln p’；卸载回弹线（SL ）方程：v =1+e = v κ-κln p’= N -(λ-κ)ln p’0-κln p’，即e = N -1-(λ-κ)ln p’0-κln p’=3.083-1-(0.16-0.04)ln300-0.04ln p’ =1.399-0.04ln p’；修正剑桥模型屈服函数为：2222f=M 0p M p p q ''''-+=，代入已知参数得屈服面方程为：223000p p q '''-+=状态A ：在有效应力300kPa 下等向固结，所以300A p kPa '=，0A q '=，根据正常固结线(NCL)方程可得：e =2.083-0.16 ln p A ’=1.17。
基于修正剑桥模型模拟理想三轴不排水试验——两种积分算法的对比分析（CZQ-SpringGod ）1、修正剑桥模型在塑性功中考虑体积塑性应变的影响，根据屈服面一致性原则，假定屈服函数对硬化参数的偏导为0，就获得了以理想三轴不排水试验为基础的修正剑桥模型屈服函数：22(,)()0c q f p q p p p M =+-= （1） 其中3kkp σ=，ij ij ij s p σδ=-，212ij ij J s s =，q =M 为临界线斜率，c p 为前期固结压力。
硬化/软化法则：p c v c dp v d p ελκ=- （2） 式中p v ε为体积塑性应变，v 为比体积，λ为正常固结线斜率，κ为回弹线斜率。
由于不排水屈服面推导过程是基于硬化参数c p 偏导为0，也就是说不排水试验中硬化参数同体积塑性应变无关，屈服面不变化，而若引入硬化法则就同屈服面推导过程中的假定矛盾，因此计算时将模型处理为理想塑性模型。
2、显式和隐式两种积分格式考虑应变增量ε∆驱动下，第n 增量步到第n+1增量步之间的应力积分格式。
显式积分格式的推导参考文献[1]，其中弹塑性矩阵中的塑性硬化模量H=0。
隐式积分格式推导如下：11()n n n p v v p p K εε++=+∆-∆ （3）111(2)n p n n v c p p ε+++∆=Λ⋅- （4） 12()n n p ij ij ij ij s s G e e +=+∆-∆ （5）1123n ijp n ij s e M ++∆=Λ （6）111112(,)()0n n n n n c qf q p p p p M +++++=+-= （7）在这一组方程中没有硬化规律方程表明为理想塑性，并将式（3）-（7）合并化简得到：1112112122(2)06()(1)0n n n n v c n n n trial c p p K K p p G q p p p M Mε++++++⎧--∆+⋅Λ⋅-=⎪⎨+-+Λ=⎪⎩ （8） 式中3(2)(2)2n n trial ij ij ij ij q s G e s G e =+∆+∆ 求解（8）式方程组即可得到n+1增量步的各个增量。
两种积分格式的matlab 程序分别见邮件附件文件夹camclayexp 和camclayimp ，显式运行主程序为camclayexp.m ，而隐式运行主程序为camclayimp.m 。
3、数值分析（1）修正剑桥模型的参数设定：临界线斜率：M=1.1正常固结线斜率：λ=0.17回弹线斜率：κ=0.034初始比体积：v 0 =2.12前期固结压力：c p =100 KPa剪切与体积模量的比值：GK=0.46155每个增量步体积模量的计算：nv K p κ= 剪切模量G=GK ×K其中固结线方程为：0ln()n v v p λ=-。
（2）计算结果：不排水有效应力路径：（a ）显示算法 (b) 隐式算法图1 不排水有效应力路径偏应力随轴向应变的变化：（a）显示算法（b）隐式算法图2 偏应变随轴向应变的变化孔隙水压力随轴向应变的变化：（a）显示算法（b）隐式算法图3 孔压随轴向应变的变化两种算法的每个增量步同屈服面的偏移程度：（a）显式算法（b）隐式算法图4 每个增量屈服面的偏移程度结论：两种算法在计算理想塑性修正剑桥模型时，数值解能很好地同理论屈服面符合。
显示算法的误差是递增的，而隐式是收敛的。
理想塑性模型的分析结果表明，经过屈服面修正后的显示算法在精度上要高于隐式算法，可能同收敛参数的设定有关，不过两者都是精确的。
参考文献：[1] S.W.Sloan. A.J.Abbo. D.Sheng. Refined explicit integration of elasoplastic models with automatic erro control[J]. Engineering Computations. 2001:18,121-19程序代码：显式积分算法：(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]=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:nsteppcre(n)=pc;Sre(n,:)=S;pre(n)=Pst;qre(n)=q;q1re(n)=q1;%% strain incrementdE=[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); % 弹性刚度矩阵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: initializationYF1=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);if YF2<=0loop=-1;S=S1; % 弹性加载,或卸载elseif YF2>0 %塑性加载if YF1<0alph=alphfun(S,ddS,pc,M);alphc=YF2/(YF2-YF1);endif YF1>=0dfp0=2*Pst-pc;dfj0=6*sJ2/M^2;dfo0=0;flow0=FlowPl(deviS,dfp0,dfj0,dfo0,J2,sJ2,lode);pW=flow0*ddS';if pW>0alph=0;elsedisp('弹性卸载，又加载')St=S+0.2*ddS;alph=alphfun(St,ddS,pc,M);endendloop=1;S=S+alph*ddS;alphre(n)=alph;sige=(1-alph)*ddS; % 找到屈服之后的弹性预测endend%% Error control of Correctortoler=0.001;iter=0;T=0;dT=1;dpv=0;while loop==1 & T<1%% first step for modification (flow is a function)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; % perfect plasticity (no hardening)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);% 塑性体积应变硬化%% second step for modificationsig2=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);%% error controlErr=(dsig2-dsig1)/2;sm=S+(dsig1+dsig2)/2; % modified stressrer=getmax(Err);rsm=getmax(sm);R=max(10^(-10),rer/rsm);Tp=0.8*(toler/R)^0.5;if R>tolerbeta=max([Tp,0.1]);dT=beta*dT;elseSc=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); % 0 为无硬化% S=sm;% dpv=dpvc;%%T=T+dT;beta=min([Tp,2]); %必须输入数组做参数dT=beta*dT;dT=min([dT,1-T]);end%% record of iterationps=Sc;[px(iter+2),pds]=deviT(ps);[aa,bb,cc,py(iter+2),dd]=invar(pds);iter=iter+1;if iter>10disp('too much iteration')breakendenddisp(['coverged iteration: ',num2str(iter)])% px=[];py=[];%% next incrementdisp(['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隐式算法：（Implicit Integration Algorithm）% function camclayimp%% 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]=invar(deviS);E=[0 0 0 0 0 0];nstep=300;de1=0.001;q1=0;for n=1:nsteppcre(n)=pc;Sre(n,:)=S;pre(n)=Pst;qre(n)=q;q1re(n)=q1;% q1-q%% strain incrementdE=[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 predictorPst1=Pst+K*dEvol;for i=1:6deviS1(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);%% Plastic correctoriter=0;toler=1e-3;dEpvol=0;devidEp=zeros(1,6);dpl=0;%% recordpx(2)=Pst1;px(1)=Pst;py1(2)=q1;py1(1)=q;% py2(2)=q1;py2(1)=q;%% iteration of residule% Pst1=Pst;while Yieldf>0res(1)=Pst1-Pst-K*dEvol+K*dpl*(2*Pst1-pc1);% res(2)=pc1-pc*exp(v1/(lam-ek)*dpl*(2*Pst1-pc1)); % hardening rule 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)])if resmax<tolerdisp(['Convergence: ',num2str(iter)])breakendif iter>=10disp('too much interation')breakenditer=iter+1;%% the Jacobian Matrix of Residule Vectorntdm(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'%% update the residulePst1=Pst1+dru(1);dpl=dpl+dru(2);% pc1=pc1+dru(3);%% record of iterationpx(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);enddisp(['increment: ',num2str(n)])px=[];py1=[];%% next incrementPst=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子程序：function y=ydfun(steff,p,pc,M)----屈服函数%% yield functiony=3*steff^2/M^2+p*(p-pc);function [p,sd]=deviT(s)----求张量的偏量p=(s(1)+s(2)+s(3))/3;for i=1:3sd(i)=s(i)-p;endfor i=4:6sd(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)...)+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;if J2==0SINT3=0;elseSINT3=-3.0*ROOT3*J3/(2.0*J2*sJ2);endif (SINT3>1)SINT3=1 ;endif(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==0flow=[1 1 1 0 0 0];returnendroot3=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/3 1/3 1/3 0 0 0];for i=1:6vn2(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:6flow(i)=c1*vn1(i)+c2*vn2(i)+c3*vn3(i); end。
