vumat计算主应力

qliao

理论：

代码：     
      DOUBLE PRECISION :: I1,I2,I3,A,B,C,E1,E2,Fai,SP1,SP2,SP3,SP
      PARAMETER (PI=3.141592653589)
C
C  CALCULATE THE STRESS INVARIANTS
      I1=stressNew(k,1)+stressNew(k, 2)+stressNew(k,3)
      I2=stressNew(k,1)*stressNew(k,2)+stressNew(k,2)*
     *    stressNew(k,3)+stressNew(k,1)*stressNew(k,3)-
     *    stressNew(k,4)**2-stressNew(k,5)**2-stressNew(k,6)**2
      I3=stressNew(k,1)*stressNew(k,2)*stressNew(k,3)
     *    -stressNew(k,1)*stressNew(k,6)**2-stressNew(k,2)
     *    *stressNew(k,5)**2-stressNew(k,3)*stressNew(k,4)**2
     *    +2*stressNew(k,4)*stressNew(k,5)*stressNew(k,6)
C
C calculate principal stress
      A=I1/3.0
      E1=I1**2-3.0*I2
      IF(E1.GT.0.0001)THEN
         B=SQRT(E1)
      ELSE
C         WRITE(*,*) 'Warning: B is less than 0.01,set to 0'
         B=0.0
      ENDIF
C      
      IF(B.GT.0)THEN
         E2=(2.0*I1**3-9.0*I1*I2+27.0*I3)/B**3/2.0
         IF(E2.GT.1.0)THEN
C          WRITE(*,*) 'Warning:cos(Fai)>1,E2= ', E2, 'adjust to 1'
            E2=1.0
         ELSEIF(E2.LT.-1.0)THEN
C           WRITE(*,*) 'Warning:cos(Fai)<-1,E2= ', E2, 'adjust to -1'
             E2=-1.0
         ENDIF
         Fai=ACOS(E2)/3.0
         SP1=A+2.0*B*COS(Fai)/3.0
         SP2=A+2.0*B*COS(Fai+2.0*PI/3.0)/3.0
         SP3=A+2.0*B*COS(Fai+4.0*PI/3.0)/3.0
      ELSE
         SP1=A
         SP2=A
         SP3=A
      ENDIF
C
C  make sure SP1>=SP2>=SP3   
      IF(SP1.LT.SP2)THEN
         SP=SP1
         SP1=SP2
         SP2=SP
      ENDIF
      IF(SP1.LT.SP3)THEN
         SP=SP1
         SP1=SP3
         SP3=SP
      ENDIF
      IF(SP2.LT.SP3)THEN
         SP=SP2
         SP2=SP3
         SP3=SP
      ENDIF

diesure

考虑的因素很全面了 防止各种出错
不过计算I3的时候有一个小错误 在explicit里s23，s13是stress里的第五个和第六个 颠倒一下就能求出和abaqus一样的解了 不然影响第二第三主应力的值
      I3(k)=stressNew(k,1)*stressNew(k,2)*stressNew(k,3)
     1    -stressNew(k,1)*stressNew(k,5)**2
     2    -stressNew(k,2)*stressNew(k,6)**2
     3    -stressNew(k,3)*stressNew(k,4)**2
     4    +2*stressNew(k,4)*stressNew(k,5)*stressNew(k,6)

qliao

要求E1.GT.0.0001，避免了对负数开根号以及分母接近0的情况
E1<0.0001时B<0.01，直接忽略2.0*B*COS(Fai)/3.0，2.0*B*COS(Fai+2.0*PI/3.0)/3.0，2.0*B*COS(Fai+4.0*PI/3.0)/3.0，误差不超过0.01MPA
通过下面的代码避免了ACOS的自变量大于1（计算误差引起的，理论上是不会大于1的）      
IF(E2.GT.1.0)THEN
     WRITE(*,*) 'Warning:cos(Fai)>1,E2= ', E2, 'adjust to 1'
     E2=1.0
ELSEIF(E2.LT.-1.0)THEN
     WRITE(*,*) 'Warning:cos(Fai)<-1,E2= ', E2, 'adjust to -1'
      E2=-1.0
ENDIF
通过这些处理避免了由计算误差引起的所有可能的奇异情况。

6.10中，vumat可以调用vsprinc计算主应力，但是本人始终没有调用成功。要是有人成功了，麻烦介绍一下经验

lxzzwn1004

qliao 发表于 2012-7-14 23:33
要求E1.GT.0.0001，避免了对负数开根号以及分母接近0的情况
E1

请问你调用成功了吗？
我需要用到最大和最小主应变，不知道你有什么看法能够实现！？

huorili

好贴留名，顶

huorili

十分感谢，游泳

陈祥华xh

楼主可以将程序源文件分享一下嘛

taojiayu瑜

求主应变的计算公式

慌张的小铁匠

学习了，太有帮助了

wkd123

好帖子啊，挺

逆时针向

顶一个，非常有用

wanghongbo

顶一个，

枯藤老树昏鸦

学习一下