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

习题：程序及解答1.主程序2.子程序2.1输入数据functionInputDataglobalgMatgNdtgEltgEInfogBct%材料属性%EPoissonrationgravity,yieldstress,*,gMat=[855700000.3017.32e30];%结点坐标%gNdt=[004001.11102.22203.33304.44405.55606.66707.77808.8890100……（数量较多，不一一罗列）370349033709013];%单元材料号ielt=size(gElt,1);%单元总数forie=1:ieltgEInfo(ie)=1;end%边界条件%type=1--fixeddisp;(1--结点固定；0--结点自由）%type=2--pointload%nodexy*,*,*,*,typegBct=[831100001;841100001;851100001;861100001;871100001;881100001;891100001;901100001;471100001;420-1100000002];2.2初始化变量值functionAssignInitValuesglobalgMatgNdtgEltgBctgEInfoglobalgKAgFAglobalgEDatagNDataiNdt=size(gNdt,1);%总结点数iElt=size(gElt,1);%总单元数gKA=zeros(2*iNdt,2*iNdt);gFA=zeros(2*iNdt,1);forie=1:iEltiMatNum=gEInfo(ie);vdMat=gMat(iMatNum,:);gEData(ie).strs_tot=zeros(1,3);gEData(ie).strs_inc=zeros(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;%等效塑性应变endforin=1:iNdtgNData(in).disp_tot=zeros(1,2);%结点in的总位移gNData(in).disp_inc=zeros(1,2);%结点in的增量位移end2.3增量法求解塑性方程functionIncreMethodglobalgNdtglobalgKAgFAglobalgEDatagNData%处理荷载边界条件（计算总荷载）LoopBoundaryLoad;%idivid=50;%荷载分段数Fdiv=gFA/idivid;%初始化向量nfre_max=size(gFA,1);Fcur=zeros(nfre_max,1);%开始增量计算foridiv=1:idivid%在屏幕上显示载荷分段数sDisplay=sprintf('Currentloaddividenumber:%d',idiv);disp(sDisplay);%当前增量步总荷载Fcur=Fcur+Fdiv;iNdt=size(gNdt,1);gKA=zeros(2*iNdt,2*iNdt);%计算并组集当前步的整体刚度矩阵LoopCalcGlobalStif;%荷载列阵初始化为0gFA=zeros(nfre_max,1);%计算前一步应力的等效荷载（初始应力部分）LoopElePrevForce;%处理位移边界条件LoopBoundaryDisp;%用于本步计算的未平衡荷载Ferr=Fcur+gFA;%solvetheequationDisp=gKA\Ferr;%nlen=size(gNdt,1);forin=1:nlen%叠加增量位移到各结点上inc_disp(1:2)=Disp(2*in-1:2*in);gNData(in).disp_inc(1:2)=inc_disp(1:2);gNData(in).disp_tot(1:2)=gNData(in).disp_tot(1:2)+inc_disp(1:2);end%forin%计算单元应力（考虑塑性返回）LoopEleStress;end%foridiv2.4求三角形单元的应变矩阵function[Bm,A]=Tri3BMtr(xe)%获取三角形单元的三个结点坐标值xm=xe;xm(4:5,1:2)=xe(1:2,1:2);%计算形函数系数（ai,bi,ci)forin=1:3ai(in)=xm(in+1,1)*xm(in+2,2)-xm(in+2,1)*xm(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(delta=1e-10)error('单元面积小于0！');end%计算单元的应变矩阵Bm=[bi(1)0bi(2)0bi(3)0;0ci(1)0ci(2)0ci(3);ci(1)bi(1)ci(2)bi(2)ci(3)bi(3)];Bm=Bm/delta;end2.5计算单元的弹性矩阵functionDe=DMtrElas(iele)%globalgMatgEInfo%MatNum=gEInfo(iele);Mat=gMat(MatNum,:);%获取单元材料向量E=Mat(1);%弹性模量u=Mat(2);%泊松比%弹性矩阵（平面应力问题）A1=E/(1-u*u);De=[A1A1*u0;A1*uA10;00A1*(1-u)/2];2.6计算单元的刚度矩阵functionKe=Tri3EleStif(ie)%globalgMatgNdtgEltglobalgEData%计算弹性矩阵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,:);%获取单元结点xe=gNdt(ENodes(:),:);%获取结点坐标[Bm,A]=Tri3BMtr(xe);%计算刚度矩阵Ke=Bm'*Dep*Bm*A;2.7组集单元刚度矩阵functionKe=Tri3EleStif(ie)%globalgMatgNdtgEltglobalgEData%计算弹性矩阵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,:);%获取单元结点xe=gNdt(ENodes(:),:);%获取结点坐标[Bm,A]=Tri3BMtr(xe);%计算刚度矩阵Ke=Bm'*Dep*Bm*A;2.8处理外力边界条件functionLoopBoundaryLoadglobalgBctgKAgFAiBct=size(gBct,1);forib=1:iBctiNode=gBct(ib,1);iType=gBct(ib,8);ifiType==2%点荷载边界条件forind=1:2irow=(iNode-1)*2+ind;gFA(irow)=gFA(irow)+gBct(ib,ind+1);endendend2.9处理位移边界条件functionLoopBoundaryDispglobalgEltgBctgKAgFAiBct=size(gBct,1);forib=1:iBctiNode=gBct(ib,1);iType=gBct(ib,8);ifiType==1%固定点边界条件forind=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在屏幕上绘制图形functionPlotResultsglobalgNdtgEltglobalgEDatagNDataholdondFau=50;%变形图放大系数%绘制变形前网格trisurf(gElt,gNdt(:,1),gNdt(:,2),zeros(size(gNdt,1),1))view(2);axisequal;axisoff;axistight;ele=size(gElt,1);%获取变形后数据iNdt=size(gNdt,1);gNdisp=zeros(iNdt,2);forin=1:iNdtgNdisp(in,1:2)=gNData(in).disp_tot(1:2);endDDisp=gNdt+gNdisp*dFau;pause(3.0);trisurf(gElt,DDisp(:,1),DDisp(:,2),zeros(size(DDisp,1),1));view(2);axisequal;axisoff;axistight;forie=1:eleifgEData(ie).strn_equ(1,1)0xe=gElt(ie,:);x=mean(gNdt((xe),1));y=mean(gNdt((xe),2));plot(x,y,'bo');endend2.11计算应力偏量functionvdSij=CalPartialStress(vdStr)%输入数据%vdStr----参与计算的总应力%输出数据%vdSij----计算得到的偏应力值%%计算给定点的应力偏量Sijistr=size(vdStr,1);ifistr4vdStr(4)=0;%对平面应力问题，增加1虚拟项endvdSij=zeros(4,1);dStc=(vdStr(1)+vdStr(2)+vdStr(4))/3.0;vdSij(1)=vdStr(1)-dStc;%x方向应力偏量vdSij(2)=vdStr(2)-dStc;%y方向应力偏量vdSij(3)=vdStr(3);vdSij(4)=vdStr(4)-dStc;%z方向应力偏量（即使平面应力，也有此偏量）2.11计算塑性矩阵functionDp=DMtrPlas(iele,ipoi)%计算单元的塑性矩阵%输入数据：%iele----单元编号%ipoi----积分点编号%输出数据%Dp-------计算得到的塑性矩阵%globalgMatgEltgEInfogEData%MatNum=gEInfo(iele);Mat=gMat(MatNum,:);%获取单元材料向量%---------------------------------%计算给定点的应力偏量SijvdStress(:,1)=gEData(iele).strs_tot(ipoi,:);%获取给定积分点的单元总应力vdSij=CalPartialStress(vdStress);%计算应力偏量的不变量dJ2=0.5*(vdSij(1)*vdSij(1)+vdSij(2)*vdSij(2)+vdSij(4)*vdSij(4))+vdSij(3)*vdSij(3);dSJ2=sqrt(dJ2);%计算矩阵C（Mises准则）vdC=[0.0;sqrt(3);0.0];%--------------------%计算矩阵A1、A2、A3vdA1=[1.0;1