ABAQUS用VUSDFLD子程序模拟结果输出存在1.#QNAN,是什么原因? 30

浏览:1261 回答:1

代码

subroutine vusdfld(
 c Read only -
     *   nblock, nstatev, nfieldv, nprops, ndir, nshr,
     *   jElem, kIntPt, kLayer, kSecPt,
     *   stepTime, totalTime, dt, cmname,
     *   coordMp, direct, T, charLength, props,
     *   stateOld,
c Write only -
     *   stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension props(nprops),
     *          jElem(nblock), coordMp(nblock, *),
     *          direct(nblock, 3, 3), T(nblock,3,3),
     *          charLength(nblock),
     *          stateOld(nblock, nstatev),
     *          stateNew(nblock, nstatev),
     *          field(nblock, nfieldv)
      character*80 cmname
C
C Local arrays from vgetvrm are dimensioned to
C maximum block size (maxblk)
C
      parameter( nrData=6 )
      character*3 cData(maxblk*nrData)
      dimension rData(maxblk*nrData), jData(maxblk*nrData)
      dimension StatusOld(NBLOCK,1),StatusNew(NBLOCK,1)
     
      PARAMETER(ZERO=0.0,ONE=1.0,TWO=2.0,THREE=3.0)
      DIMENSION STRESSNEW(NDIR+NSHR)    
C
C     INPUT CONSTANT PARAMETERS
C
        EMOD=80000
        ENU=0.33
        DC1=5.0992
        DC2=0.7247
        DC3=0.3611
        EA=301.0
        EQ=205.0
        EB=18.0
       
        EBULK3=EMOD/(ONE-TWO*ENU)
        EG2=EMOD/(ONE+ENU)
        EG=EG2/TWO
        EG3=THREE*EG
        ELAM=(EBULK3-EG2)/THREE
      
      DO K0=1,NBLOCK
       
C
C  INITIAL FAILURE FLAGS FROM STATE VARIABLES
C
       DAMAGE=stateOld(K0,2)
       CALL VGETVRM('S',RDATA,JDATA,CDATA,JSTATUS)
         STRESSNEW(1)=RDATA(K0*1)
         STRESSNEW(2)=RDATA(K0*2)
         STRESSNEW(3)=RDATA(K0*3)
         STRESSNEW(4)=RDATA(K0*4)
         IF(NSHR.GT.1) THEN
         STRESSNEW(5)=RDATA(K0*5) 
         STRESSNEW(6)=RDATA(K0*6)
         END IF
      
       CALL VGETVRM('PEEQ',RDATA,JDATA,CDATA,JSTATUS)
       EQPLAS=RDATA(K0)
      
       IF(NSHR.EQ.1) THEN
         POTENIAL=0.50*((STRESSNEW(1)-STRESSNEW(2))**2+
     1   (STRESSNEW(2)-STRESSNEW(3))**2+(STRESSNEW(3)-STRESSNEW(1))**2
     2   +6*STRESSNEW(4)**2)
       ELSE
         POTENIAL=0.50*((STRESSNEW(1)-STRESSNEW(2))**2+
     1   (STRESSNEW(2)-STRESSNEW(3))**2+(STRESSNEW(3)-STRESSNEW(1))**2
     2   +6*(STRESSNEW(4)**2+STRESSNEW(5)**2+STRESSNEW(6)**2))
       END IF
         EFS=SQRT(POTENIAL)
      
         HARD=0.0
       IF(EQPLAS.EQ.0.0) THEN
         SYIELD=EA
       ELSE
         HARD=EQ*EB*EXP(-1*EB*EQPLAS)
         SYIELD=EA+EQ*(1-EXP(-1*EB*EQPLAS))
       END IF
      
         RHS=EFS-SYIELD
C
C  JUDGE IF YIELD
C
       IF(RHS.GT.0.0) THEN
         DEQPLAS=RHS/(EG3+HARD)
        
         STRESSM=(STRESSNEW(1)+STRESSNEW(2)+STRESSNEW(3))/3
       
        HJ2=((STRESSNEW(1)-STRESSNEW(2))**2+(STRESSNEW(2)-STRESSNEW(3))
     1      **2+(STRESSNEW(3)-STRESSNEW(1))**2+6.0*STRESSNEW(4)**2)/6.0
        HJ3=(STRESSNEW(1)-STRESSM)*(STRESSNEW(2)-STRESSM)*
     1    (STRESSNEW(3)-STRESSM)-(STRESSNEW(3)-STRESSM)*STRESSNEW(4)**2
        THETA=ACOS(3.0**1.5*HJ3/(2.0*HJ2**1.50))/3
       
        CR11=COS(THETA)
        CR12=SIN(THETA)
        CR21=-1.0*SIN(THETA)
        CR22=COS(THETA)
       
        SP1=CR11**2*STRESSNEW(1)+CR21**2*STRESSNEW(2)+
     1      2*CR11*CR21*STRESSNEW(4)
        SP2=CR12**2*STRESSNEW(1)+CR22**2*STRESSNEW(2)+
     1      2*CR12*CR22*STRESSNEW(4)
        SP3=STRESSNEW(3)
       
        HYSTRESS=(SP1+SP2+SP3)/3.0
     
        TNH=HYSTRESS/EFS
       
        SP1=MAX(SP1,SP2,SP3)
        SP3=MIN(SP1,SP2,SP3)
       
        DDAMAGE=(((SP1-SP3)/EFS)**DC1+(TNH-1.0/3.0))**DC2*DEQPLAS
     1        /DC3
       
      
       END IF
      
       DAMAGE=DAMAGE+DDAMAGE
       stateNew(K0,2)=DAMAGE
       
       IF(stateNew(K0,2).GE.1) THEN
          stateNew(K0,1)=0.0
       ELSE
          stateNew(K0,1)=1.0
       END IF
        
     
      END DO  
      RETURN
      END

邀请回答 我来回答

全部回答

(1)
默认 最新
沉默的果冻
"1.#QNAN" 是浮点数运算中出现的特殊值,表示“不是一个数”,出现这样的输出通常意味着程序在执行过程中出现了某些数学计算上的错误或者非法操作,导致了浮点数的异常值。我也出现了同样的问题,正在找bug...
2月18日
评论 点赞

没解决?试试专家一对一服务

换一批