弹塑性求位移pde.doc

文档介绍：――-―――――――――――弹塑性求位移pde――-―――――――――――**已知弹性矩阵、上步总应力、上步塑性内变量,上步增量应变,求解本步应力增量、塑性矩阵。$c6callgetstr(p,str,de,dv,d,dp,1,prag,prager)……**=+$c6doi=1,6$c6dv(i)=dv(i)+str(i)$c6enddo……**方程右端项load=+[u_i]*f_i-[evp_i]*dv_i――-―――――――――――+――-―――――――――――**已知、、,,求解、。含义:X(塑性材料参数),U(,),DE(),DV(),D(),DP(),KKK(塑性内变量选取),F1(屈服函数值),FSUB(屈服函数)SUBROUTINEGETSTR(X,U,DE,DV,D,DP,KKK,F1,FSUB)implicitreal*8(a-h,o-z)DIMENSIONH(4),U(4),DE(3),DV(3),D(3,3),DP(3,3),V(4),DU(3),X(*),P(12),DF(4)EXTERNALF,fsub**以下计算针对二维N=4NDF=3E=x(2)***试探应力增量DO40I=1,NDFC==1,NDFC=C+D(I,J)*DE(J)20CONTINUEDU(I)=C40CONTINUE**DO100I=1,NV(I)=U(I)IF()V(I)=V(I)+DU(I)100CONTINUE**F1=F1=FSUB(X,V)**F1<=0,为弹性状态加载或塑性卸载,,=0,结束计算IF(.-E)THENDO150I=1,NDFDV(I)=DU(I)150CONTINUERETURNENDIF**F0=F0=FSUB(X,U)**F1>0,F0<0,为弹性状态进入塑性状态,>0R=(.-E*10)R=F0/(F0-F1)**DO200I=1,NDFV(I)=U(I)+R*DU(I)CONTINUE**CALLGETDP(X,V,D,DP,KKK,FSUB)**DO400I=1,NDFV(I)=U(I)+DU(I)C==1,NDFC=C+DP(I,J)*DE(J)300CONTINUEDV(I)=CCONTINUE**F1>0,F0<0,为弹性状态进入塑性状态,**F1>0,F0=0,为从一种塑性状态到另一种塑性状态,,DD==1,NDFIF(.-E*10)THENV(I)=V(I)-DV(I)ELSEV(I)=V(I)-(1.-R)*DV(I)ENDIFCONTINUE**求CALLGETDF(X,V,DF,FSUB)DO450I=1,NDFDV(I)=DF(I)450CONTINUER=**由、、作为初值调用root反复调整比例因子,将应力拉回屈服面上。拉回方法采用法向拉回。Root采用割线法求解,将应力点B拉回点C位置。循环调整比例因子~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~**含义:R(比例因子),DD(割线法第二个点位置),E(迭代精度),p(、、、塑性材料参数),F(函数),fsub(屈服函数),k(控制迭代结束标识)4