代码 subroutine vusdfld( c Read only - * nblock, nstatev, nfieldv, nprops, ndir, ns

"1.#QNAN" 是浮点数运算中出现的特殊值，表示“不是一个数”，出现这样的输出通常意味着程序在执行过程中出现了某些数学计算上的错误或者非法操作，导致了浮点数的异常值。我也出现了同样的问题，正在找bug...

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

浏览：1268 回答：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