5_3练习题-三角形单元塑性程序.doc

1、习题：程序及解答1.主程序2.子程序2.1 输入数据function InputDataglobal gMat gNdt gElt gEInfo gBct %材料属性% E Poisson ration gravity ,yield stress, *,gMat=85570000 0.3 0 17.32e3 0;%结点坐 标% gNdt=0 040 01.111 02.222 03.333 04.444 05.556 06.667 07.778 08.889 010 0（数量较多，不一一罗列）370 3 490 3 37090 1 3;% 单 元材料号ielt=size(gElt,1); %单

2、元总数for ie=1:ieltgEInfo(ie)=1;end%边界条件% type=1-fixed disp;(1-结点固定；0-结点自由）% type=2-point load% node x y *, *, *, *, typegBct= 83 1 1 0 0 0 0 1;84 1 1 0 0 0 0 1;85 1 1 0 0 0 0 1;86 1 1 0 0 0 0 1;87 1 1 0 0 0 0 1;88 1 1 0 0 0 0 1;89 1 1 0 0 0 0 1;90 1 1 0 0 0 0 1;47 1 1 0 0 0 0 1;42 0 -11000 0 0 0 0 2;2

3、.2 初始化变量值function AssignInitValuesglobal gMat gNdt gElt gBct gEInfoglobal gKA gFAglobal gEData gNDataiNdt=size(gNdt,1);%总结点数iElt=size(gElt,1);%总单元数gKA=zeros(2*iNdt,2*iNdt);gFA=zeros(2*iNdt,1);for ie=1:iEltiMatNum=gEInfo(ie);vdMat=gMat(iMatNum,:);gEData(ie).strs_tot=zeros(1,3);gEData(ie).strs_inc=zer

4、os(1,3);gEData(ie).strs_yds=vdMat(4);%初始屈服应力gEData(ie).strs_yield=vdMat(4);%屈服应力gEData(ie).strs_equ=0.0;%等效应力gEData(ie).strn_equ=0.0;%等效塑性应变endfor in=1:iNdtgNData(in).disp_tot=zeros(1,2);%结点in的总位移gNData(in).disp_inc=zeros(1,2);%结点in的增量位移end2.3 增量法求解塑性方程function IncreMethodglobal gNdtglobal gKA gFAgl

5、obal gEData gNData%处理荷 载边界条件（ 计算总荷载）LoopBoundaryLoad;%idivid=50;%荷载分段数Fdiv=gFA/idivid;%初始化向量nfre_max=size(gFA,1);Fcur=zeros(nfre_max,1);%开始增量计算for idiv=1:idivid%在屏幕上显示载荷分段数sDisplay=sprintf(Current load divide number:%d,idiv);disp(sDisplay);%当前增量步总荷载Fcur=Fcur+Fdiv;iNdt=size(gNdt,1);gKA=zeros(2*iNdt,2

6、*iNdt);%计 算并组 集当前步的整体刚度矩阵LoopCalcGlobalStif;%荷 载列阵 初始化为0gFA=zeros(nfre_max,1);%计 算前一步 应力的等效荷载（初始应力部分）LoopElePrevForce;%处 理位移 边界条件LoopBoundaryDisp;%用于本步计算的未平衡荷载Ferr=Fcur+gFA;%solve the equationDisp=gKAFerr;%nlen=size(gNdt,1);for in=1:nlen %叠加增量位移到各结点上inc_disp(1:2)=Disp(2*in-1:2*in);gNData(in).disp_in

7、c(1:2)=inc_disp(1:2);gNData(in).disp_tot(1:2)=gNData(in).disp_tot(1:2)+inc_disp(1:2);end %for in%计 算单元 应力（考虑塑性返回）LoopEleStress;end% for idiv2.4 求三角形单元的应变矩阵function Bm,A=Tri3BMtr(xe)%获取三角形 单 元的三个结点坐标值xm=xe;xm(4:5,1:2)=xe(1:2,1:2);%计算形函数系数（ai,bi,ci)for in=1:3ai(in)=xm(in+1,1)*xm(in+2,2)-xm(in+2,1)*xm(

8、in+1,2); %ai/aj/ambi(in)=xm(in+1,2)-xm(in+2,2); %bi/bj/bmci(in)=xm(in+2,1)-xm(in+1,1); %ci/cj/cmenddelta=sum(ai); %单元面积*2A=delta/2;if(delta0.999*dYieldStrs)%当前积分点进入塑性Dp=DMtrPlas(ie,1);Dep=De-Dp;end%计算应变 矩阵ENodes=gElt(ie,:);%获取单元结点xe=gNdt(ENodes(:),:);%获取结点坐标Bm,A=Tri3BMtr(xe);%计算刚 度矩阵Ke=Bm*Dep*Bm*A;2

9、.7 组集单元刚度矩阵function Ke=Tri3EleStif(ie)%global gMat gNdt gEltglobal gEData%计算弹 性矩阵De=DMtrElas(ie);Dep=De; %赋值弹塑性矩阵%dEquStrs=gEData(ie).strs_equ(1,1); %等效应力dYieldStrs=gEData(ie).strs_yield(1,1);%当前屈服应力if(dEquStrs0.999*dYieldStrs)%当前积分点进入塑性Dp=DMtrPlas(ie,1);Dep=De-Dp;end%计算应变 矩阵ENodes=gElt(ie,:);%获取单元结

10、点xe=gNdt(ENodes(:),:);%获取结点坐标Bm,A=Tri3BMtr(xe);%计算刚 度矩阵Ke=Bm*Dep*Bm*A;2.8 处理外力边界条件function LoopBoundaryLoadglobal gBct gKA gFAiBct=size(gBct,1);for ib=1:iBctiNode=gBct(ib,1);iType=gBct(ib,8);if iType=2 %点荷载边界条件for ind=1:2irow=(iNode-1)*2+ind;gFA(irow)=gFA(irow)+gBct(ib,ind+1);endendend2.9 处理位移边界条件fu

11、nction LoopBoundaryDispglobal gElt gBct gKA gFAiBct=size(gBct,1);for ib=1:iBctiNode=gBct(ib,1);iType=gBct(ib,8);if iType=1 %固定点边界条件for ind=1:2irow=(iNode-1)*2+ind;if(gBct(ib,ind+1)=1)gKA(irow,:)=0;gKA(:,irow)=0;gKA(irow,irow)=1.0;gFA(irow)=0.0;endendendend2.10 在屏幕上绘制图形function PlotResultsglobal gNdt

12、 gElt global gEData gNDatahold ondFau=50; %变形图放大系数%绘制变 形前网格trisurf(gElt,gNdt(:,1),gNdt(:,2),zeros(size(gNdt,1),1)view(2);axis equal;axis off; axis tight;ele=size(gElt,1);%获取变 形后数据iNdt=size(gNdt,1);gNdisp=zeros(iNdt,2);for in=1:iNdtgNdisp(in,1:2)=gNData(in).disp_tot(1:2);endDDisp=gNdt+gNdisp*dFau;pau

13、se(3.0);trisurf(gElt,DDisp(:,1),DDisp(:,2),zeros(size(DDisp,1),1);view(2);axis equal; axis off; axis tight;for ie=1:eleif gEData(ie).strn_equ(1,1)0xe=gElt(ie,:);x=mean(gNdt(xe),1);y=mean(gNdt(xe),2);plot(x,y,bo);endend2.11 计算应力偏量function vdSij=CalPartialStress(vdStr)%输入数据% vdStr-参与计算的 总应力%输出数据% vdSi

14、j-计算得到的偏 应力值%计算给 定点的 应力偏量Sijistr=size(vdStr,1);if istr0.999*dPrevYieldStrs)%前一步进入塑性if(dEquStrsdPrevEquStrs)%增量步后也进入塑性dMep=1.0;endelse%前一步 为弹 性if(dEquStrsdPrevYieldStrs)%增量步后进入塑性dMep=(dEquStrs-dPrevYieldStrs)/(dEquStrs-dPrevEquStrs);endend%当前步到达弹 性极限时的总应力vdStrB=strs_prev+(1-dMep)*inc_strs;vdStrInc=ze

15、ros(3,1);%塑性应力增量dPlasStrn=0.0;if dMep0.0%部分进入塑性gEData(ie).strs_tot(1,1:3)=vdStrB(1:3,1);%重新赋值单元总应力为弹性极限应力Dp=DMtrPlas(ie,1);%当前步B 点的塑性矩阵vdStrIncE=De*inc_stra*dMep;%应力增量的弹 性部分vdStrIncP=Dp*inc_stra*dMep;%应力增量的塑性部分vdStrInc=vdStrIncE-vdStrIncP;dPlasStrn=0.0;%计算等效塑性应变增量dPlasStrn=CalPlasStrain(ie,vdStrB,vd

16、StrIncE);%计算等效塑性应变增量end%塑性调 整后的 总应力vdStrA=vdStrB+vdStrInc;vdSij=CalPartialStress(vdStrA);dJ2=0.5*(vdSij(1)*vdSij(1)+vdSij(2)*vdSij(2)+vdSij(4)*vdSij(4)+vdSij(3)*vdSij(3);dEquStrs=sqrt(3)*sqrt(dJ2);%求等效应力dRback=1.0;dCurYieldStrs=dPrevYieldStrs;if(dEquStrsdPrevYieldStrs)dHp=0.0;%对理想弹塑性dCurYieldStrs=dC

17、urYieldStrs+dHp*dPlasStrn;dRback=dCurYieldStrs/dEquStrs;endvdStrC=vdStrA*dRback;%-gEData(ie).strs_tot(1,1:3)=vdStrC(1:3,1);%重新赋值gEData(ie).strs_yield(1,1)=dCurYieldStrs;%当前步的屈服应力gEData(ie).strs_equ(1,1)=dEquStrs*dRback;%当前步的等效应力gEData(ie).strn_equ(1,1)=gEData(ie).strn_equ(1,1)+dPlasStrn;%总的等效塑性应变2.1

18、3 计算并组集单元初应力的等效荷载列阵function LoopElePrevForce%global gElt gNdtglobal gFAglobal gEData%组集单 元初应 力的等效荷载列阵elen=size(gElt,1);for ie=1:elen%计 算应变 矩阵ENodes(:,1)=gElt(ie,:);%获取单元结点xe=gNdt(ENodes(:),:);%获取结点坐标Bm,A=Tri3BMtr(xe);%获 取给定 积分点的单元总应力vdStr(:,1)=gEData(ie).strs_tot(1,:);%计 算等效荷 载force=Bm*vdStr*A;%组 集f

19、or il=1:3inod=ENodes(il);gFA(2*inod-1:2*inod,1)=gFA(2*inod-1:2*inod,1)-force(2*il-1:2*il);end %for ilend2.14 计算单元等效塑性应变function dPlasStrn=CalPlasStrain(ie,vdStrB,vdStrIncE)global gMat gEInfo%求B 点 时等效 应力vdSij=CalPartialStress(vdStrB);dJ2=0.5*(vdSij(1)*vdSij(1)+vdSij(2)*vdSij(2)+vdSij(4)*vdSij(4)+vdSi

20、j(3)*vdSij(3);dEquStrsB=sqrt(3)*sqrt(dJ2);%-%求B 点 时塑性 应变增量和应 力增量偏量的比例系数dr%第一步：求B 点时对应 的流动矢量ad%按表1 计算矩 阵C（Mises准则）vdC=0.0;sqrt(3);0.0;%按公式（1.27）计算矩阵A1 、A2、A3vdA1=1.0;1.0;0.0;1.0;vdA2=0.5*vdSij/sqrt(dJ2);vdA2(3)=vdSij(3)/sqrt(dJ2);vdA3=zeros(4,1);vdA3(1)=vdSij(2)*vdSij(4)+dJ2/3.0;vdA3(2)=vdSij(1)*vdSi

21、j(4)+dJ2/3.0;vdA3(3)=-2.0*vdSij(4)*vdSij(3);vdA3(4)=vdSij(1)*vdSij(2)-vdSij(3)*vdSij(3)+dJ2/3.0;%按公式（1.28）计算流动矢量advdAd1=vdA1*vdC(1)+vdA2*vdC(2)+vdA3*vdC(3);istr=3;vdAd=zeros(3,1);for i=1:istrvdAd(i)=vdAd1(i); %只用前3项 ;且须定义vdAd 为列向量，否则如果没有vdAd=zeros(3,1)这个语句，单凭vdAd(i)=vdAd1(i)得到的vdAd默认为行向量，虽然vdAd1为列向量

22、end%第二步：求B 点时对应 的计算参数AdHp=0.0;%对理想塑性材料%第三步：求弹性矩阵DeDe=DMtrElas(ie);%第四步：求材料处在塑性阶段时的应变增量inc_str，即从B 点到D点的应变增量（等于C点到A点的弹性应变 增量*应力调整系数）inc_str=inv(De)*vdStrIncE;%第五步求出dr（根据公式1.20）dr=vdAd*De*inc_str/(dHp+vdAd*De*vdAd);%-%根据公式（1.32）等效塑性应变增量dPlasStrn=dr*vdAd*vdStrB/dEquStrsB;计算单元等效塑性应变程序如上所述，变形图见右图所示 ： gEData(1).strn_equans =0故单元1的位移为 0。