qliao 发表于 2012-7-14 23:24:13

vumat计算主应力

本帖最后由 qliao 于 2012-7-17 21:27 编辑

理论：

代码：   
      DOUBLE PRECISION :: I1,I2,I3,A,B,C,E1,E2,Fai,SP1,SP2,SP3,SP
      PARAMETER (PI=3.141592653589)
C
CCALCULATE 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
Cmake 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 发表于 2013-6-3 18:34:59

考虑的因素很全面了 防止各种出错
不过计算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 发表于 2012-7-14 23:33:57

本帖最后由 qliao 于 2012-7-17 21:28 编辑

要求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 发表于 2013-6-1 19:53:53

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

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

huorili 发表于 2017-3-22 21:17:37

好贴留名，顶

huorili 发表于 2017-3-22 21:22:06

十分感谢，游泳

陈祥华xh 发表于 2018-4-19 11:05:11

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

taojiayu瑜 发表于 2019-12-22 10:58:18

求主应变的计算公式

慌张的小铁匠 发表于 2020-2-25 12:08:33

学习了，太有帮助了

wkd123 发表于 2020-7-4 15:29:53

好帖子啊，挺

逆时针向 发表于 2021-6-23 21:18:23

顶一个，非常有用

wanghongbo 发表于 2021-7-6 17:02:44

顶一个，

枯藤老树昏鸦 发表于 2022-7-7 14:44:01

学习一下
页: [1]
查看完整版本: vumat计算主应力