ABAQUS用VUSDFLD子程序模拟结果输出存在1.#QNAN,是什么原因? 30
代码
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