abaqus中UMAT子程序

九鸢

朋友们，我在计算子程序时发现，直接运行显示NAN，但是Cnf是有计算出来的值的；如果将Cnf直接设为计算出来的值，计算没有问题，下面是我的核心程序，UMAT开头格式已省略，全文源代码网址Fcode Pastebin

      integer i, j, k, l
      real(KIND=8)  ZERO, ONE, TWO, THREE, A0, n, q1, q2, v, ABFF,ABFF1,
     1 pivot, Cnf, k0, q4, Smid, S0, DTIME1, k1, A, E1, Ac, Bc, kapa0,
     2 kapa, AA, choose1, choose2, choose3, AAA, dtime2, AAAA, Smid1,
     3 SS, SSS, SSSS, SSSSS, ssssss, SA, Smid2,sssssss,ACA,ABFF2,
     4 Cnv, Taomiu, Amiu, Cnv1,abff4,abff5,timenext,alphal,alphar,
     5 alphamid,na,Cnv2, VV, aging,elastic, GC, Cnf1
      parameter(ZERO=0.0d0, TWO=2.0d0, ONE=1.0d0, ACA=1.0d1, n=1.0d-1,
     1 A0=0.0d0, THREE=3.0d0, na=1.9d0)
      real(kind=8), dimension(ntens) :: dstrainv
      real(kind=8), dimension(ntens) :: dstrainf
      real(kind=8), dimension(ntens) :: gamamiu
      real(kind=8), dimension(ntens,ntens) :: Iddsdde
      real(kind=8), dimension(ntens,ntens) :: IIddsdde
      real(kind=8), dimension(ntens,2*ntens) :: aug
      real(kind=8), dimension(ntens) :: DSTRESS1
      real(kind=8), dimension(ntens) :: DSTRESS2
      real(kind=8), dimension(ntens) :: STRESS2
      real(kind=8), dimension(ntens) :: STRESS1
      real(kind=8), dimension(2*ntens) :: ABC
      real(kind=8), dimension(2*ntens) :: ABD
      real(kind=8), dimension(2*ntens) :: ABE
      real(kind=8), dimension(2*ntens) :: ABF
      real(kind=8), dimension(ntens,ntens) :: Dve1
      real(kind=8), dimension(ntens,ntens) :: DD
      real(kind=8), dimension(ntens,ntens) :: II
      real(kind=8), dimension(ntens,2*ntens) :: aug1
      real(kind=8), dimension(ntens,ntens) :: Dve
      real(kind=8), dimension(ntens) :: PS
      real(kind=8), dimension(ntens) :: SB
      LOGICAL flok

C material properties
      q1 = props(1)     
      v = props(2)      
      q2 = props(3)      
      Ac = props(4)  
      Bc = props(5)
      k0 = props(6)
      q4 = props(7)
      S0 = props(8) * 1.0d0  
      kapa0 = props(9)
      E1 = props(10)
      
C
      if (KINC.EQ.ONE) then
          DTIME1 = ZERO
          do i = 1, ntens
              DSTRESS1(i) = ZERO
              STRESS1(i) = ZERO
          end do
      else
          DTIME1 = STATEV(1) / (2.4d1 * 3.6d3)
          Smid = STATEV(14)
          do i = 1, ntens
              k1 = i + 1.0d0
              STRESS1(i) = STATEV(k1)
          end do
          do i = 1, ntens
              k1 = i + 7.0d0
              DSTRESS1(i) = STATEV(k1)
          end do
      end if
      
C G matrix
      do i = 1, ntens
        do j = 1, ntens
          Iddsdde(i,j) = 0.0d0
        end do
      end do
      do i = 1, ndi
        do j = 1, ndi
          Iddsdde(i,j) = -v
        end do
        Iddsdde(i,i) = 1.0d0
      end do
      do i = ndi+1, ntens
        Iddsdde(i,i) = TWO * (ONE + v)
      end do
      dtime2 = dtime / (2.4d1 * 3.6d3)
      timenext = time(1) / (2.4d1 * 3.6d3) + 5.0d-1 * dtime2
      do i = 1,ntens
          STRESS2(i) = stress(i)
          DSTRESS2(i) = DSTRESS1(i)
      end do

C viscous
      choose2 = 1
      if (choose2.EQ.1) then
          do i = 1, ntens
              dstrainf(i) = 0.0d0
          end do
          Cnf = ZERO
          if (KINC.EQ.ONE) then
              Smid = S0
          else
              Smid = Smid
          end if
          Smid1 = ((1.0d0 + 2.0d0 * k0 * dtime / 2.4d1 / 3.6d3 * Smid)
     1 ** 5d-1 - 1.0d0) / (k0 * dtime / 2.4d1 / 3.6d3)
          Smid = Smid1
          Cnf1 = 5.0d-1 * q4 * k0 * Smid * dtime2
          do i =1, ntens
              A = 0.0d0
              do j = 1, ntens
                  SB(i) = STRESS(i)
                  A = A + Iddsdde(i,j) * stress(j)
              end do
              dstrainf(i) = q4 * k0 * Smid * A * dtime2
          end do
      else
          Cnf = 0
          do i = 1,ntens
              dstrainf(i) = 0
          end do
      end if
      Cnf = Cnf1
      
C Damage
      elastic = 0
      choose3 = 0
      if (choose3.EQ.1) then
          AA = 0.0d0
          CALL SPRINC(STRAN,PS,2,NDI,NSHR)
          do i = 1, ndi
              AA = AA + 1.0d-4
          end do
          kapa = kapa0 + AA ** 5.0d-1
          d = 1.0d0 - kapa0 * (1.0d0 - Ac) / kapa - Ac /
     1 (exp(Bc * (kapa - kapa0)))
      else
          d = 0
      end if

C E
      if (choose3.EQ.1) then
          AAA = 1.0d0
      else
          AAA = 0.0d0
      end if
      if (elastic.EQ.1) then
          q1 = q1
      else
          q1 = 0.0d0
      end if
      do i = 1, ntens
        do j = 1, ntens
            GC = Iddsdde(i,j) * (Cnv + Cnf + q1 + AAA / E1)
            IIddsdde(i,j) = GC
        end do
      end do
      do i = 1, ntens
        do j = 1, ntens
            Iddsdde(i,j) = IIddsdde(i,j)
        end do
      end do
      
C Dve1
      if (choose3.EQ.1) then
          do i = 1, ntens
              do j = 1,ntens
                  II(i,j) = 0
              end do
          end do
          do i = 1,ntens
              II(i,i) = 1
          end do
          do i = 1,ntens
              do j =1,ntens
                  Dve1(i,j) = ((1 - d) * E1 * II(i,j) + d *
     1 Iddsdde(i,j)) / ((1 - d) * E1)
              end do
          end do
      else
          do i = 1,ntens
              do j = 1,ntens
                  Dve1(i,j) = 0.0d0
              end do
          end do
      end if
              
C Inversion of the matrix Iddsdde
      aug(:,1:TWO*ntens) = 0.0d0
      aug(:,1:ntens) = Iddsdde(:,1:ntens)
      aug(:,ntens+1:TWO*ntens) = 0.0d0
      do k = 1, ntens
          aug(k,k+ntens) = 1.0d0
      end do
      do i = 1, ntens
          pivot = aug(i,i)
          do l = 1, 2*ntens
              aug(i,l) = aug(i,l) / pivot
          end do
          do j = 1, ntens
              if (j /= i) then
                  pivot = aug(j,i)
                  do k = 1, 2*ntens
                      aug(j,k) = aug(j,k) - pivot * aug(i,k)
                  end do
              end if
          end do
      end do
      DD(:,1:ntens) = aug(:,ntens+1:TWO*ntens)
      
C Inversion of the matrix Dve1
      if (choose3.EQ.1) then
          aug1(:,1:TWO*ntens) = 0.0d0
          aug1(:,1:ntens) = Dve1(:,1:ntens)
          aug1(:,ntens+1:TWO*ntens) = 0.0d0
          do i = 1, ntens
              aug1(i,i+ntens) = 1.0d0
          end do
          do i = 1, ntens
              pivot = aug1(i,i)
              do j = 1, 2*ntens
                  aug1(i,j) = aug1(i,j) / pivot
              end do
              do j = 1, ntens
                  if (j /= i) then
                      pivot = aug1(j,i)
                      do k = 1, 2*ntens
                          aug1(j,k) = aug1(j,k) - pivot * aug1(i,k)
                      end do
                  end if
              end do
          end do
          Dve(:,1:ntens) = aug1(:,ntens+1:TWO*ntens)
      else
      end if

C Reall Dve
      do l = 1,ntens
          do k = 1,ntens
              ddsdde(l,k) = 0
          end do
      end do
      if (choose3.EQ.1) then
          do l = 1,ntens
              do j = 1,ntens
                  do k = 1,ntens
                      ddsdde(l,j) = ddsdde(l,j) + Dve(l,k) * DD(k,j)
                  end do
              end do
          end do
      else
          do l = 1,ntens
              do j = 1,ntens
                  ddsdde(l,j) = ddsdde(l,j) + DD(l,j)
              end do
          end do
      end if
      
C STATEV
      k1 = 1
      STATEV(k1) = dtime
      do i = 1, ntens
          k1 = i + 1
          STATEV(k1) = stress(i)
      end do
      do i = 1, ntens
          k1 = i + 7
          do j = 1, ntens
              STATEV(k1) = ddsdde(i,j) * (dstran(j) -
     1 dstrainv(j) - dstrainf(j))
          end do
      end do
      STATEV(14) = Smid

C Stress increment evaluation
      do l = 1, ntens
          do k = 1, ntens
              stress(l) = stress(l) + ddsdde(l,k) * (dstran(k) -
     1 dstrainv(k) - dstrainf(k))
          end do
      end do

九鸢

全文网址http://p.fcode.cn/_gonbP