《高等土力学课程》由会员分享，可在线阅读，更多相关《高等土力学课程（11页珍藏版）》请在金锄头文库上搜索。
1、基于修正剑桥模型模拟理想三轴不排水试验两种积分算法的对比分析（CZQ-SpringGod）1、修正剑桥模型 在塑性功中考虑体积塑性应变的影响，根据屈服面一致性原则，假定屈服函数对硬化参数的偏导为 0，就获得了以理想三轴不排水试验为基础的修正剑桥模型屈服函数：1)f (p, q)=M+p (p -匕)=0s = 一 p5 ,ij ij ijJ = 1 ss , q 七3ss22 ij ij2 ij ijM 为临界线斜率， p 为 c前期固结压力。 硬化/软化法则：2)dpvc = d 8 p -pv 九Kc式中8 p为体积塑性应变，V为比体积，九为正常固结线斜率，K为回弹线斜率。v由于不排水屈服
2、面推导过程是基于硬化参数 p 偏导为 0，也就是说不排水试验中硬化 c参数同体积塑性应变无关，屈服面不变化，而若引入硬化法则就同屈服面推导过程中的假定矛盾，因此计算时将模型处理为理想塑性模型。隐式积分格式推导如下：n+1 p = n p + K (Ae n+1A8 p )vv( 3 )n+1 Ae p = n+1A (2 n+1 p p )vc( 4 )n+1s = ns + 2G(Ae Aep)ijijijij( 5 )3n+1sAep = n+1Aj( 6 )ijM 2f ( n+1q, n+1 p) =q + n+1 p( n+1 p p ) = 0( 7)2、显式和隐式两种积分格式考虑
3、应变增量Ae驱动下，第n增量步到第n+1增量步之间的应力积分格式。显式积分 格式的推导参考文献，其中弹塑性矩阵中的塑性硬化模量H=0。M2c在这一组方程中没有硬化规律方程表明为理想塑性，并将式（3） -（7）合并化简得到：Vq 2+ n+1trialn+1 p np KAs + K - n+1 A(2 - n+1 p p ) = 0 vc6gp( n+1 p p ) M 2(1+ n+1A )2 = 0式中qtrialcM 23(ns + 2GAe )(ns + 2GAe )2ijijijij求解(8)式方程组即可得到n+1增量步的各个增量。两种积分格式的matlab程序分别 见邮件附件文件夹
4、camclayexp和camclayimp，显式运行主程序为camclayexp.m,而隐式运 彳亍主程序为 camclayimp.m。3、数值分析(1) 修正剑桥模型的参数设定：临界线斜率：M=1.1正常固结线斜率：九=0.17回弹线斜率：k =0.034初始比体积：v0 =2.12前期固结压力：p =100 KPac剪切与体积模量的比值：GK=0.46155剪切模量G=GKXKv每个增量步体积模量的计算：K = npk其中固结线方程为：v二v 九ln(np)。 0(2) 计算结果：不排水有效应力路径：liS5 何Ihvwlh 幣旦炽1 irlEir-iFlKinfhfrlrmwn sin-
5、js p KP&|、胡 d Ejjrfag cnlKallriG-啊100创也 切 如宙盘二富雷営ETnp睿血$ pam twinID ZD 3a 4D 90 6a ?Dylikl卫吁述窃crilitd Gna*5:srinfiEsm paih(b)隐式算法(a)显示算法图1不排水有效应力路径偏应力随轴向应变的变化:(a)显示算法(b)隐式算法图2偏应变随轴向应变的变化孔隙水压力随轴向应变的变化:图3孔压随轴向应变的变化两种算法的每个增量步同屈服面的偏移程度:(a)显式算法(b)隐式算法图 4 每个增量屈服面的偏移程度结论：两种算法在计算理想塑性修正剑桥模型时，数值解能很好地同理论屈服面符合。
6、 显示算法的误差是递增的，而隐式是收敛的。理想塑性模型的分析结果表明，经过屈服面修 正后的显示算法在精度上要高于隐式算法，可能同收敛参数的设定有关，不过两者都是精确 的。参考文献：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程序代码：显式积分算法：(Explicit Integration Algorithm)% function c
7、amclayexp% 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;dEpv
8、ol=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);dEv
9、ol=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 YF20%塑性加载if YF1=0 dfp0=2*Pst-pc;dfj0=6*sJ2/MA2;dfo0=0; flow0=Flow
10、Pl(deviS,dfp0,dfj0,dfo0,J2,sJ2,lode); pW=flow0*ddS;if pW0alph=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 & Ttoler beta=max(Tp,0.l); 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 为无硬化
