我编的三维梁的线弹性VUMAT程序，求围观，求指导！！！

nuaaxhf

C
C USER SUBROUTINE VUMAT
      SUBROUTINE VUMAT (
C READ ONLY -
     *     NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL,
     *     STEPTIME, TOTALTIME, DT, CMNAME, COORDMP, CHARLENGTH,
     *     PROPS, DENSITY, STRAININC, RELSPININC,
     *     TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
     *     STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
     *     TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,
C WRITE ONLY -
     *     STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW )
C
      INCLUDE 'VABA_PARAM.INC'
C
      DIMENSION COORDMP(NBLOCK,*), CHARLENGTH(NBLOCK), PROPS(NPROPS),
     1     DENSITY(NBLOCK), STRAININC(NBLOCK,NDIR+NSHR),
     2     RELSPININC(NBLOCK,NSHR), TEMPOLD(NBLOCK),
     3     STRETCHOLD(NBLOCK,NDIR+NSHR),
     4     DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
     5     FIELDOLD(NBLOCK,NFIELDV), STRESSOLD(NBLOCK,NDIR+NSHR),
     6     STATEOLD(NBLOCK,NSTATEV), ENERINTERNOLD(NBLOCK),
     7     ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
     8     STRETCHNEW(NBLOCK,NDIR+NSHR),
     9     DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
     1     FIELDNEW(NBLOCK,NFIELDV),
     2     STRESSNEW(NBLOCK,NDIR+NSHR), STATENEW(NBLOCK,NSTATEV),
     3     ENERINTERNNEW(NBLOCK), ENERINELASNEW(NBLOCK)
C     
      CHARACTER*80 CMNAME
      DIMENSION INTV(2)
      PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0, THREE = 3.D0,
     *     THIRD = ONE / THREE, HALF = 0.5D0, TWOTHDS = TWO / THREE,
     *     OP5 = 1.5D0 )
      PARAMETER ( TEMPFINAL = 1.D2, TIMEFINAL = 1.D-2 )
C
*
*   
      E      = PROPS(1)
      XNU    = PROPS(2)
*     
      TWOMU  = E / ( ONE + XNU )
      ALAMDA = TWOMU * XNU / ( ONE - TWO * XNU )
*     
*     IF STEPTIME EQUALS TO ZERO, ASSUME THE MATERIAL PURE ELASTIC
*     AND USE INITIAL ELASTIC MODULUS
*
c        IF(STEPTIME.EQ.ZERO)THEN
c        STRESS
      IF ( STEPTIME .EQ. ZERO ) THEN     
        DO K = 1, NBLOCK
*     TRIAL STRESS
          STRAININC(K,2) = (-XNU)*STRAININC(K,1)
          STRAININC(K,3) = (-XNU)*STRAININC(K,1)
          TRACE = STRAININC(K,1)+STRAININC(K,2)+STRAININC(K,3)
          STRESSNEW(K,1) = STRESSOLD(K,1) +
     *                  TWOMU*STRAININC(K,1) + ALAMDA*TRACE
          STRESSNEW(K,2) = ZERO
          STRESSNEW(K,3) = ZERO
          STRESSNEW(K,4) = STRESSOLD(K,4) + TWOMU*STRAININC(K,4)
          STRESSNEW(K,5) = ZERO
          STRESSNEW(K,6) = ZERO
        END DO
      ELSE
*     
        DO K = 1, NBLOCK
c        write(400,*) 'COORDMP=',COORDMP(K,*)
*     TRIAL STRESS
      STRAININC(K,2) = (-XNU)*STRAININC(K,1)
      STRAININC(K,3) = (-XNU)*STRAININC(K,1)
      TRACE = STRAININC(K,1)+XNU*STRAININC(K,2)+XNU*STRAININC(K,3)
          STRESSNEW(K,1) = STRESSOLD(K,1) +
     *                  TWOMU*STRAININC(K,1) + ALAMDA*TRACE
      STRESSNEW(K,2) = ZERO
      STRESSNEW(K,3) = ZERO
      STRESSNEW(K,4) = STRESSOLD(K,4) + TWOMU*STRAININC(K,4)
      STRESSNEW(K,5) = ZERO
      STRESSNEW(K,6) = ZERO

C-----UPDATE THE SPECIFIC INTERNAL ENERGY--------
      STRESSPOWER = HALF *
     1   (  ( STRESSOLD(K,1)+STRESSNEW(K,1) )*STRAININC(K,1)
     2    + ( STRESSOLD(K,2)+STRESSNEW(K,2) )*STRAININC(K,2)
     3    + ( STRESSOLD(K,3)+STRESSNEW(K,3) )*STRAININC(K,3)
     4    + ( STRESSOLD(K,4)+STRESSNEW(K,4) )*STRAININC(K,4)*TWO  )
C
      ENERINTERNNEW(K) = ENERINTERNOLD(K)
     1                  + STRESSPOWER / DENSITY(K)
C
          END DO
        END IF
*     
*     
      RETURN
      END

nuaaxhf

首先告诉大家，我的这个程序肯定有问题，就看你能不能看出来了！

lxzzwn1004

没看出来 哪里有问题

nuaaxhf

lxzzwn1004 发表于 2013-5-25 22:21
没看出来 哪里有问题

     IF ( STEPTIME .EQ. ZERO ) THEN     
        DO K = 1, NBLOCK
*     TRIAL STRESS
          STRAININC(K,2) = (-XNU)*STRAININC(K,1)
          STRAININC(K,3) = (-XNU)*STRAININC(K,1)
          TRACE = STRAININC(K,1)+STRAININC(K,2)+STRAININC(K,3)
          STRESSNEW(K,1) = STRESSOLD(K,1) +
     *                  TWOMU*STRAININC(K,1) + ALAMDA*TRACE
          STRESSNEW(K,2) = ZERO
          STRESSNEW(K,3) = ZERO
          STRESSNEW(K,4) = STRESSOLD(K,4) + TWOMU*STRAININC(K,4)
          STRESSNEW(K,5) = ZERO
          STRESSNEW(K,6) = ZERO
        END DO
      ELSE
*   
请问这一部分需要吗？

diesure

其实我也不是很明白 不过STRAININC不应该是read only 吗。。。你怎么给他赋值了呢

nuaaxhf

diesure 发表于 2013-6-6 03:00
其实我也不是很明白 不过STRAININC不应该是read only 吗。。。你怎么给他赋值了呢 ...

For beams in space you must define the thickness strains, strainInc(*,2) and strainInc(*,3)。
这个是手册中的一句话

princewhb

nuaaxhf 发表于 2013-5-27 15:01
IF ( STEPTIME .EQ. ZERO ) THEN     
        DO K = 1, NBLOCK
*     TRIAL STRESS

需要。。

lvhongtao666

没看出来