弹塑性求位移pde

1、 弹塑性求位移pde *已知弹性矩阵、上步总应力、上步塑性内变量，上步增量应变，求解本步应力增量、塑性矩阵。$c6 call getstr(p,str,de,dv,d,dp,1,prag,prager) *$c6 do i=1,6$c6 dv(i) = dv(i)+str(i)$c6 enddo *方程右端项load = +u_i*f_i-evp_i*dv_i stres.for plas.for *已知、，求解、。含义：X(塑性材料参数)，U(,)，DE()，DV()，D()，DP()，KKK(塑性内变量选取)，F1(屈服函数值)，FSUB(屈服函数)SUBROUTINE GETSTR(X,

2、U,DE,DV,D,DP,KKK,F1,FSUB) implicit real*8 (a-h,o-z) DIMENSION H(4),U(4),DE(3),DV(3),D(3,3),DP(3,3),V(4),DU(3),X(*),P(12),DF(4) EXTERNAL F,fsub*以下计算针对二维 N = 4 NDF = 3 E = x(2)*0.0001*试探应力增量 DO 40 I=1,NDF C = 0.0 DO 20 J=1,NDF C = C+D(I,J)*DE(J)20 CONTINUE DU(I) = C40 CONTINUE*DO 100 I=1,N V(I) = U(I)

3、 IF (I.LE.NDF) V(I) = V(I)+DU(I)100 CONTINUE*F1F1 = FSUB(X,V)*F1<=0，为弹性状态加载或塑性卸载，0，结束计算IF (F1.LT.-E) THEN DO 150 I=1,NDF DV(I) = DU(I)150 CONTINUE RETURN ENDIF*F0 F0 = FSUB(X,U)*F1>0，F0<0，为弹性状态进入塑性状态， >0 R = 0.0 IF (F0.LT.-E*10) R = F0/(F0-F1)* DO 200 I=1,NDF V(I) = U(I)+R*DU(I)200 CONTI

4、NUE* CALL GETDP(X,V,D,DP,KKK,FSUB)* DO 400 I=1,NDF V(I) = U(I)+DU(I) C = 0.0 DO 300 J=1,NDF C = C+DP(I,J)*DE(J)300 CONTINUE DV(I) = C400 CONTINUE*F1>0，F0<0，为弹性状态进入塑性状态，*F1>0，F00，为从一种塑性状态到另一种塑性状态， DD=0.3 DO 430 I=1,NDF IF (F0.GT.-E*10) THEN V(I) = V(I)-DV(I) ELSE V(I) = V(I) - (1.-R)*DV(I) E

5、NDIF430 CONTINUE*求 CALL GETDF(X,V,DF,FSUB) DO 450 I=1,NDF DV(I) =DF(I)450 CONTINUE R = 1.0*由、作为初值调用root反复调整比例因子，将应力拉回屈服面上。拉回方法采用法向拉回。Root采用割线法求解，将应力点B拉回点C位置。循环调整比例因子*含义：R(比例因子)，DD(割线法第二个点位置)，E(迭代精度)，p(、塑性材料参数)，F(函数)，fsub(屈服函数)，k(控制迭代结束标识)470 DO 500 I=1,N P(I) = V(I) IF (I.LT.N) P(I+N) = DV(I) P(I+N*

6、2) = X(I)500 CONTINUE CALL ROOT(R,DD,E,p,F,fsub,k) *修正 if(k.eq.1) then do 550 i=1,ndf v(i)=v(i)-(1.-r)*dv(i)550 continue *修正 call getdf(x,v,df,fsub) do 570 i=1,ndf dv(i)=df(i)570 continue r = 1.0 dd=-0.1 goto 470循环调整比例因子 end if *设迭代n次，*应力增量 DO 600 I=1,NDF V(I) = V(I) - (1.-R)*DV(I) DV(I) = V(I)-U(I)

7、600 CONTINUE*由计算最终塑性矩阵 CALL GETDP(X,V,D,DP,KKK,FSUB) F1=FSUB(X,V) RETURN END*已知、，求解。含义：X(塑性材料参数)，U(,)，D()，DP()，KKK(塑性内变量选取)，yield(屈服函数) SUBROUTINE GETDP(X,U,D,DP,KKK,yield) implicit real*8 (a-h,o-z) DIMENSION U(4),D(3,3),DP(3,3),X(*),DF(4),DDF(3) EXTERNAL YIELD N = 4 NDF = 3*求 CALL GETDF(X,U,DF,yiel

8、d)*求 DO 200 I=1,NDF C = 0.0 DO 100 J=1,NDF C = C+D(I,J)*DF(J)100 CONTINUE DDF(I) = C200 CONTINUE* A = 0.0 AM = 0.0 DO 400 I=1,NDF A = A+DF(I)*DDF(I) IF (KKK.EQ.1) AM = AM+U(I)*DF(I) IF (KKK.EQ.2 .AND. I.LE.2) AM = AM+DF(I) IF (KKK.EQ.3) AM = AM+DF(I)*DF(I)400 CONTINUE IF (KKK.EQ.3) AM = SQRT(AM)* 内变

9、量只用了塑性功 A = A - AM*DF(N)* DO 600 I=1,NDF DO 500 J=1,NDF DP(I,J) = DDF(I)*DDF(J)/A500 CONTINUE600 CONTINUE RETURN END *已知、，求解和含义：X(塑性材料参数)，U(,)，DF(,)，yield(屈服函数) SUBROUTINE GETDF(X,U,DF,yield) implicit real*8 (a-h,o-z) DIMENSION H(4),U(4),X(*),DF(4),V(4) N = 4*采用中心差商的方法求偏导，选取步长h。 DO 10 K = 1,N H(K) =

10、 x(2)*0.01 V(K) = U(K)10 CONTINUE* DO 50 K = 1,N V(K) = U(K) + H(K) F2 = yield(X,V) V(K) = U(K) - H(K) F1 = yield(X,V) DF(K) = (F2-F1)/H(K)*0.5 V(K) = U(K)50 CONTINUE RETURN END *已知、，求解值含义：X(塑性材料参数)，U(,) real*8 FUNCTION PRAGER(X,U) implicit real*8 (a-h,o-z) DIMENSION U(4),S(3),x(4) ，分别为材料的初始内聚力、初始摩擦

11、系数。为屈服面强化系数。 tgo=x(1) q0=x(4) c=x(2)+q0*U(4) PV= X(3) D = (U(1)+U(2)*(1+pv)/3. S(1) = U(1) - D S(2) = U(2) - D S(3) = U(3) UK = (u(1)+u(2)*pv-D SJ = (S(1)*2+S(2)*2+UK*2)/2.+S(3)*2 alpha=tgo/dsqrt(9+12*tgo*2) ck=3*c/dsqrt(9+12*tgo*2) PRAGER = DSQRT(SJ) + ALPHA*D*3 - CK RETURN END*已知、，求解比例因子应力在屈服面上满足，

12、求解，等价于求中的（通常称为沿法线方向拉回）。采用割线法求解，梁老师编写含义：x1 (比例因子)，d (割线法第二个点位置)，e (迭代精度)，p (、塑性材料参数)，f (函数)，fsub(屈服函数)，k(控制迭代结束标识) subroutine root(x1,d,e,p,f,fsub,k) implicit real*8 (a-h,o-z) dimension p(*) nliq=01099 k=0 nliq=nliq+1 f1 = f(x1,p,fsub) if (abs(f1).lt.e.or.nliq.ge.20) goto 999 if (f1.gt.0.0) then x2 =

13、 x1 - d else x2 = x1 + d endif f2 = f(x2,p,fsub) n12 = 010 if (f1*f2 .le. 0.0) go to 30 if (f1*f1 .lt. f2*f2) go to 20 x1 = x2 + 1.1*(x2-x1) f1 = f(x1,p,fsub) if (n12.ge.1) then x=x2 go to 85 end if n12 = 1 goto 1020 x2 = x1 + 1.1*(x1-x2) f2 = f(x2,p,fsub) if (n12.le.-1) then x=x1 go to 85 end if n1

14、2 = -1 go to 1030 if (x1 .ge. x2) go to 40 xm1 = x1 xm2 = x2 go to 5040 xm1 = x2 xm2 = x150 x3 = (x1*f2-x2*f1)/(f2-f1) if (x3 .ge. xm1) go to 60 x3 = xm160 if (x3 .le. xm2) go to 70 x3 = xm270 f3 = f(x3,p,fsub) if (abs(f3) .lt. e) go to 90 if (abs(x3).gt.1e-6.and.abs(x3-x1)/x3).lt.1.0e-6) go to 90 if (abs(x3).gt.1e-6.and.abs(x3-x2)/x3).lt.1.0e-6) go to 90 if (f2*f3 .ge. 0.0) then f2 = f3 x2 = x3 else f1 = f3 x1 = x3 end if go to 5085 continue k=1 x1=x if(d.gt.1.e17)return d=d*2. goto 1099990 continue return90 x1