复合材料vumat单元畸变

divinelove

仿照一个博士的论文做的，因为没有发表所以不能提供给大家，但是基本思想是用hashin准则，考虑纤维断裂，纤维挤压，机体开裂，机体激烈，机纤剪切，五种模式，子程序和inp一起附上，看那位高手可以解决一下，算到一定程度就算不下去了，而且我还想知道怎么样才可以知道是那种破坏模式导致了网格畸变，谢谢。

divinelove

自己回答一下。程序改用应变的起始破坏准则，可以算下去，具体的原因在考虑中。或许是因为有限元本身事先给出应变然后积分得到应力，应力的变化比应变更加剧烈吧。。。

divinelove

把子程序贴出来吧，急待高手提点意见，是不是因为我用的是应力判据啊，我看有用应变判据的，我主要是模仿论文里面的，人家用的是应力的，。。而且我用单个单元测试了，可以发生折减的，具体的测试，我晚籼时候发过来，

       subroutine vumat(

c Read only -

     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,

     2  stepTime, totalTime, dt, cmname, coordMp, charLength,

     3  props, density, strainInc, relSpinInc,

     4  tempOld, stretchOld, defgradOld, fieldOld,

     5  stressOld, stateOld, enerInternOld, enerInelasOld,

     6  tempNew, stretchNew, defgradNew, fieldNew,

c Write only -

     7  stressNew, stateNew, enerInternNew, enerInelasNew )

c

      include 'vaba_param.inc'

c

c 3D Orthotropic Elasticity with Hashin 3d Failure criterion

c

c The state variables are stored as:

c    state(*,1)   = material point status

c    state(*,2:7) = damping stresses

c

c User defined material properties are stored as

c  * First line:

c     props(1) --> Young's modulus in 1-direction, E1

c     props(2) --> Young's modulus in 2-direction, E2

c     props(3) --> Young's modulus in 3-direction, E3

c     props(4) --> Poisson's ratio, nu12

c     props(5) --> Poisson's ratio, nu13

c     props(6) --> Poisson's ratio, nu23

c     props(7) --> Shear modulus, G12

c     props(8) --> Shear modulus, G13

c

c  * Second line:

c     props(9)  --> Shear modulus, G23

c     props(10) --> beta damping parameter

c     props(11) --> "not used"

c     props(12) --> "not used"

c     props(13) --> "not used"

c     props(14) --> "not used"

c     props(15) --> "not used"

c     props(16) --> "not used"

c

c  * Third line:

c     props(17) --> Ultimate tens stress in 1-direction, sigu1t

c     props(18) --> Ultimate comp stress in 1-direction, sigu1c

c     props(19) --> Ultimate tens stress in 2-direction, sigu2t

c     props(20) --> Ultimate comp stress in 2-direction, sigu2c

c     props(21) --> Ultimate tens stress in 2-direction, sigu3t

c     props(22) --> Ultimate comp stress in 2-direction, sigu3c

c     props(23) --> "not used"

c     props(24) --> "not used"

c

c  * Fourth line:

c     props(25) --> Ultimate shear stress, sigu12

c     props(26) --> Ultimate shear stress, sigu13

c     props(27) --> Ultimate shear stress, sigu23

c     props(28) --> "not used"

c     props(29) --> "not used"

c     props(30) --> "not used"

c     props(31) --> "not used"

c     props(32) --> "not used"

c

      dimension props(nprops), density(nblock),

     1  coordMp(nblock,*),

     2  charLength(*), strainInc(nblock,ndir+nshr),

     3  relSpinInc(nblock,nshr), tempOld(nblock),

     4  stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr),

     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),

     6  stateOld(nblock,nstatev), enerInternOld(nblock),

     7  enerInelasOld(nblock), tempNew(*),

     8  stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr),

     9  fieldNew(nblock,nfieldv), stressNew(nblock,ndir+nshr),

     1  stateNew(nblock,nstatev),

     2  enerInternNew(nblock), enerInelasNew(nblock)

*

      character*80 cmname

*

      parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = 0.5d0 )

      real

     *    dft,dfc,dmt,dmc,ds,gg

*

      parameter(

     *     i_svd_DmgFiberT   = 1,

     *     i_svd_DmgFiberC   = 2,

     *     i_svd_DmgMatrixT  = 3,

     *     i_svd_DmgMatrixC  = 4,

     *     i_svd_DmgDels   = 5,

     *     i_svd_Strain   = 6,

c     *    i_svd_StrainXx = 6,

c     *    i_svd_StrainYy = 7,

c     *    i_svd_StrainZz = 8,

c     *    i_svd_StrainXy = 9,

c     *    i_svd_StrainYz = 10,

c     *    i_svd_StrainZx = 11,

     *     n_svd_required = 11 )

*

* Structure of property array

      parameter (

     *     i_pro_E1    = 1,

     *     i_pro_E2    = 2,

     *     i_pro_E3    = 3,

     *     i_pro_nu12  = 4,

     *     i_pro_nu13  = 5,

     *     i_pro_nu23  = 6,

     *     i_pro_G12   = 7,

     *     i_pro_G13   = 8,

     *     i_pro_G23   = 9,

*

*

     *     i_pro_sigu1t = 17,

     *     i_pro_sigu1c = 18,

     *     i_pro_sigu2t = 19,

     *     i_pro_sigu2c = 20,

     *     i_pro_sigu3t = 21,

     *     i_pro_sigu3c = 22,

     *     i_pro_sigu12 = 25,

     *     i_pro_sigu13 = 26,

     *     i_pro_sigu23 = 27 )

* Read material properties

*

      E1 = props(i_pro_E1)

      E2 = props(i_pro_E2)

      E3 = props(i_pro_E3)

      xnu12 = props(i_pro_nu12)

      xnu13 = props(i_pro_nu13)

      xnu23 = props(i_pro_nu23)

      G12 = props(i_pro_G12)

      G13 = props(i_pro_G13)

      G23 = props(i_pro_G23)

*

      xnu21 = xnu12 * E2 / E1

      xnu31 = xnu13 * E3 / E1

      xnu32 = xnu23 * E3 / E2

*

*

*

      f1t = props(i_pro_sigu1t)

      f1c = props(i_pro_sigu1c)

      f2t = props(i_pro_sigu2t)

      f2c = props(i_pro_sigu2c)

      f3t = props(i_pro_sigu3t)

      f3c = props(i_pro_sigu3c)

      f12 = props(i_pro_sigu12)

      f13 = props(i_pro_sigu13)

      f23 = props(i_pro_sigu23)

C      write(*,*)'goog',f1t,f1c,f2t,f2c,f3t,f3c,f12,

C     *       f23,E1,E2,E3,xnu12,xnu13,xnu23,G12,G23,G13

C--------------------------------------------------------------------

      do 100 k=1,nblock

      stateNew(k,6)=stateOld(k,6)+strainInc(k,1)

      stateNew(k,7)=stateOld(k,7)+strainInc(k,2)

      stateNew(k,8)=stateOld(k,8)+strainInc(k,3)

      stateNew(k,9)=stateOld(k,9)+strainInc(k,4)

      stateNew(k,10)=stateOld(k,10)+strainInc(k,5)

      stateNew(k,11)=stateOld(k,11)+strainInc(k,6)

C      write(*,*)'strainfirst',stateNew

C-------------------------------------------------

      stateNew(k,1)=stateOld(k,1) 纤维拉伸

      stateNew(k,2)=stateOld(k,2)纤维压缩

      stateNew(k,3)=stateOld(k,3)基体拉伸

      stateNew(k,4)=stateOld(k,4)基体压缩

      stateNew(k,5)=stateOld(k,5)基线剪切

C------------------------------------------------

      dft=zero

      dfc=zero

      dmt=zero

      dmc=zero

      ds=zero

C---------------------------------------------------

      if (stateNew(k,1) .EQ. one) then

         E1=E1*0.07

         E2=E2*0.07

         G12=G12*0.07

         G13=G13*0.07

      end if

      if (stateNew(k,2) .EQ. one) then

         E1=E1*0.14

         E2=E2*0.14

         G12=G12*0.14

         G13=G13*0.14

      end if

      if (stateNew(k,3) .EQ. one) then

         E2=E2*0.2

         G12=G12*0.2

         G23=G23*0.2

      end if

      if (stateNew(k,4) .EQ. one) then

         E2=E2*0.4

         G12=G12*0.4

         G23=G23*0.4

      end if

      if (stateNew(k,5) .EQ. one) then

         G12=G12*zero

         G13=G13*zero

      end if

C      write(*,*)'goog',f1t,f1c,f2t,f2c,f3t,f3c,f12,

C     *       f23,E1,E2,E3,xnu12,xnu13,xnu23,G12,G23,G13

*       Compute terms of stiffness matrix

      gg = one / ( one - xnu12*xnu21 - xnu23*xnu32 - xnu31*xnu13

     *     - two*xnu21*xnu32*xnu13 )

      C11  = E1 * ( one - xnu23*xnu32 ) * gg

      C22  = E2 * ( one - xnu13*xnu31 ) * gg

      C33  = E3 * ( one - xnu12*xnu21 ) * gg

      C12  = E2 * ( xnu12 + xnu32*xnu13 ) * gg

      C13  = E3 * ( xnu13 + xnu12*xnu23 ) * gg

      C23  = E3 * ( xnu23 + xnu21*xnu31 ) * gg

C       write(*,*)'gg',gg,C11,C22,C33,C12,C13,C23

C---------------------------------------------------------------------

      stressNew(k,1)=C11*stateNew(k,6)+C12*stateNew(k,7)

     *   +C13*stateNew(k,8)

      stressNew(k,2)=C12*stateNew(k,6)+C22*stateNew(k,7)

     *   +C23*stateNew(k,8)

      stressNew(k,3)=C13*stateNew(k,6)+C23*stateNew(k,7)

     *   +C33*stateNew(k,8)

      stressNew(k,4)=G12*two*stateNew(k,9)

      stressNew(k,5)=G23*two*stateNew(k,10)

      stressNew(k,6)=G13*two*stateNew(k,11)

C      write(*,*)'stress',stressNew

C---------------------------------------------------------------------

      if ((stateNew(k,1) .LT. one) .or. (stateNew(k,2) .LT. one)) then

      if (stressNew(k,1) .GE. zero) then

        dft=(stressNew(k,1)/f1t)**two

      else

        dfc=(stressNew(k,1)/f1c)**two

      end if

      end if

      if ((stateNew(k,3) .LT. one) .or. (stateNew(k,4) .LT. one)) then

      if (stressNew(k,2) .GE. zero) then

        dmt=(stressNew(k,2)/f2t)**two+(stressNew(k,4)/f12)**two

     *     +(stressNew(k,5)/f23)**two

      else

        dmc=(stressNew(k,2)/(two*f12))**two

     *     +(stressNew(k,2)/f2c)*((f2c/(two*f12))**2-one)

     *     +(stressNew(k,4)/f12)**two+(stressNew(k,5)/f23)**two

      end if

      end if

      if ((stateNew(k,5) .LT. one)) then

      if (stressNew(k,1) .GE. zero) then

        ds=(stressNew(k,4)/f12)**two+(stressNew(k,6)/f13)**two

      else

        ds=(stressNew(k,1)/f1c)**two

     *       +(stressNew(k,4)/f12)**two+(stressNew(k,6)/f13)**two

      end if

      end if

C      write(*,*)'df',dft,dfc,dmt,dmc,ds

C------------------------------------------------------------------

      if (dft .GE. one) then

        stateNew(k,1)= one

      end if

      if (dfc .GE. one) then

        stateNew(k,2)= one

      end if

      if (dmt .GE. one) then

        stateNew(k,3)= one

      end if

      if (dmc .GE. one) then

        stateNew(k,4)= one

      end if

      if (ds .GE. one) then

        stateNew(k,5)= one

      end if

C      write(*,*)'statenew',statenew

  100  continue

      return

      end

ido88

divinelove 发表于 2012-6-3 12:25
把子程序贴出来吧，急待高手提点意见，是不是因为我用的是应力判据啊，我看有用应变判据的，我主要是模仿论 ...

算不下去应该是单元发生畸变，也即变形过大；所以算的不收敛了就算不下去了。把失效单元删了就可以继续算下去了，具体的那种失效形式可以通过你定义的几个SDV量看到。比如你定义了SDV1-SDV5，在后处理中就可以到了。

divinelove

ido88 发表于 2012-6-3 20:25
算不下去应该是单元发生畸变，也即变形过大；所以算的不收敛了就算不下去了。把失效单元删了就可以继续算 ...

不可以删掉的，因为从实验看没有伤到纤维，只有机体的损伤。感觉算不下去还有其他的原因，我在找了，一句一句的找，nnd

ido88

divinelove 发表于 2012-6-3 21:16
不可以删掉的，因为从实验看没有伤到纤维，只有机体的损伤。感觉算不下去还有其他的原因，我在找了，一句 ...

你是怀疑子程序出错了？一句一句地找？如果在单胞上能通过就说明它是没有问题的..

divinelove

w又有问题了，这个是我用上面的子程序做的冲击模型，这个是sdv4显示的是基体的拉伸损伤，为什么周围部分会使间断的不是我想象的那样连续的。

上表面也是一圈一圈的，

这个是接触力曲线

不知道真么回事，是我的子程序编的有问题？还是？？？

divinelove

zhanghaidong 发表于 2012-6-24 20:18
应该还是子程序的问题，楼主解决没有

还没有，你也试试看

zhanghaidong

应该还是子程序的问题，楼主解决没有

qliao

楼主我下了你的程序，非常感谢！
你程序的下面这部分有问题：
if (dft .GE. one) then
        stateNew(k,1)= one
      end if

      if (dfc .GE. one) then
        stateNew(k,2)= one
      end if

      if (dmt .GE. one) then
        stateNew(k,3)= one
      end if

      if (dmc .GE. one) then
        stateNew(k,4)= one
      end if

      if (ds .GE. one) then
        stateNew(k,5)= one
      end if

比如：假定单元1积分点1在第10个增量步结束时满足dft=1。在第11个增量步对其刚度进行折减，其应力水平大幅降低df=0。
在第11个增量步，除stateNew(k,1)=stateOld(k,1)外，你的程序将不会对stateNew(k,1)额外赋值，abaqus不会将stateNew(k,1)的值保存到第12个增量步的stateOld(k,1)中
（除非你在当前增量步对其赋值），第12个增量步单元1积分点1上stateOld(1)将取0,将不会对刚度进行折减，刚度将恢复到损伤前的水平
（即你的刚度折减只在第11个增量步有效，并没有继承下来），这个可能就是你的损伤区域和未损伤区域相间的原因

divinelove

qliao 发表于 2012-7-13 21:07
楼主我下了你的程序，非常感谢！
你程序的下面这部分有问题：
if (dft .GE. one) then

      大概明白你说的意思了，但是我的程序开始先把Old值写入到New里面了，
     stateNew(k,1)=stateOld(k,1) 纤维拉伸
      stateNew(k,2)=stateOld(k,2)纤维压缩
      stateNew(k,3)=stateOld(k,3)基体拉伸
      stateNew(k,4)=stateOld(k,4)基体压缩
      stateNew(k,5)=stateOld(k,5)基线剪切
     后面进行判断，如果dft大于1就进行替换，小于则保持Old的值，感觉是连续的啊。

divinelove

divinelove 发表于 2012-7-13 21:51
大概明白你说的意思了，但是我的程序开始先把Old值写入到New里面了，
     stateNew(k,1)=stateOld ...

另外我看到程序的前面有问题，就是如果纤维拉伸和压缩同时破坏，会在原来的基础上折减两次，
if (stateNew(k,1) .EQ. one) then
         E1=E1*0.07
         E2=E2*0.07
         G12=G12*0.07
         G13=G13*0.07
      end if
      if (stateNew(k,2) .EQ. one) then
         E1=E1*0.14
         E2=E2*0.14
         G12=G12*0.14
         G13=G13*0.14
会让这个值越来越小，是一个原因，我现在的程序是纤维破坏值折减一次，选哪个严重的折减，这样子可以算下去，但是和实验结果差距很大，所以应该还有其他的问题，也在考虑，很是头疼。。。还是多多交流，一但调试成功，想留在网站上，希望大家不要再和我一样头疼，呵呵

divinelove

上传一下最新的程序，比原来的效果好，但是还是与实验结果会很大，希望大家多交流：
主程序依然贴在下面，
程序修改了一下几个方面：
1.刚度折减，认为纤维，集体，剪切，同时发生，选择最严重的，而不是连续折减，
2刚度矩阵，改用全部写法，不再改用对称形式，
另外不明白的是这个泊松比，是应该放到折减后的模量计算，还是折减前的模量计算，另外本人还没有见到系统介绍这种折减方法的理论，如果哪位有的话，可以放上来分享一下。谢谢
      subroutine vumat(
c Read only -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2  stepTime, totalTime, dt, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
c Write only -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
c
      include 'vaba_param.inc'
c
c 3D Orthotropic Elasticity with Hashin 3d Failure criterion
c
c The state variables are stored as:
c    state(*,1)   = material point status
c    state(*,2:7) = damping stresses
c
c User defined material properties are stored as
c  * First line:
c     props(1) --> Young's modulus in 1-direction, E1
c     props(2) --> Young's modulus in 2-direction, E2
c     props(3) --> Young's modulus in 3-direction, E3
c     props(4) --> Poisson's ratio, nu12
c     props(5) --> Poisson's ratio, nu13
c     props(6) --> Poisson's ratio, nu23
c     props(7) --> Shear modulus, G12
c     props(8) --> Shear modulus, G13
c
c  * Second line:
c     props(9)  --> Shear modulus, G23
c     props(10) --> beta damping parameter
c     props(11) --> "not used"
c     props(12) --> "not used"
c     props(13) --> "not used"
c     props(14) --> "not used"
c     props(15) --> "not used"
c     props(16) --> "not used"
c
c  * Third line:
c     props(17) --> Ultimate tens stress in 1-direction, sigu1t
c     props(18) --> Ultimate comp stress in 1-direction, sigu1c
c     props(19) --> Ultimate tens stress in 2-direction, sigu2t
c     props(20) --> Ultimate comp stress in 2-direction, sigu2c
c     props(21) --> Ultimate tens stress in 2-direction, sigu3t
c     props(22) --> Ultimate comp stress in 2-direction, sigu3c
c     props(23) --> "not used"
c     props(24) --> "not used"
c
c  * Fourth line:
c     props(25) --> Ultimate shear stress, sigu12
c     props(26) --> Ultimate shear stress, sigu13
c     props(27) --> Ultimate shear stress, sigu23
c     props(28) --> "not used"
c     props(29) --> "not used"
c     props(30) --> "not used"
c     props(31) --> "not used"
c     props(32) --> "not used"
c

      dimension props(nprops), density(nblock),
     1  coordMp(nblock,*),
     2  charLength(*), strainInc(nblock,ndir+nshr),
     3  relSpinInc(nblock,nshr), tempOld(nblock),
     4  stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr),
     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(*),
     8  stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv), stressNew(nblock,ndir+nshr),
     1  stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock)
*
      character*80 cmname
*
      parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = 0.5d0 )
      real

     *    dft,dfc,dmt,dmc,ds,gg
*
      parameter(
     *     i_svd_DmgFiberT   = 1,
     *     i_svd_DmgFiberC   = 2,
     *     i_svd_DmgMatrixT  = 3,
     *     i_svd_DmgMatrixC  = 4,
     *     i_svd_DmgDels   = 5,
     *     i_svd_Strain   = 6,
c     *    i_svd_StrainXx = 6,
c     *    i_svd_StrainYy = 7,
c     *    i_svd_StrainZz = 8,
c     *    i_svd_StrainXy = 9,
c     *    i_svd_StrainYz = 10,
c     *    i_svd_StrainZx = 11,
     *     n_svd_required = 11 )
*
* Structure of property array
      parameter (
     *     i_pro_E1    = 1,
     *     i_pro_E2    = 2,
     *     i_pro_E3    = 3,
     *     i_pro_nu12  = 4,
     *     i_pro_nu13  = 5,
     *     i_pro_nu23  = 6,
     *     i_pro_G12   = 7,
     *     i_pro_G13   = 8,
     *     i_pro_G23   = 9,
*
*
     *     i_pro_sigu1t = 17,
     *     i_pro_sigu1c = 18,
     *     i_pro_sigu2t = 19,
     *     i_pro_sigu2c = 20,
     *     i_pro_sigu3t = 21,
     *     i_pro_sigu3c = 22,
     *     i_pro_sigu12 = 25,
     *     i_pro_sigu13 = 26,
     *     i_pro_sigu23 = 27 )
* Read material properties
*
      E1 = props(i_pro_E1)
      E2 = props(i_pro_E2)
      E3 = props(i_pro_E3)
      xnu12 = props(i_pro_nu12)
      xnu13 = props(i_pro_nu13)
      xnu23 = props(i_pro_nu23)
      G12 = props(i_pro_G12)
      G13 = props(i_pro_G13)
      G23 = props(i_pro_G23)
*
      f1t = props(i_pro_sigu1t)
      f1c = props(i_pro_sigu1c)
      f2t = props(i_pro_sigu2t)
      f2c = props(i_pro_sigu2c)
      f3t = props(i_pro_sigu3t)
      f3c = props(i_pro_sigu3c)
      f12 = props(i_pro_sigu12)
      f13 = props(i_pro_sigu13)
      f23 = props(i_pro_sigu23)
      E1x = E1
      E2x = E2
      E3x = E3
      G12x = G12
      G13x = G13
      G23x = G23
C      write(*,*)'goog',f1t,f1c,f2t,f2c,f3t,f3c,f12,
C     *       f23,E1,E2,E3,xnu12,xnu13,xnu23,G12,G23,G13
C--------------------------------------------------------------------
      do 100 k=1,nblock
      stateNew(k,6)=stateOld(k,6)+strainInc(k,1)
      stateNew(k,7)=stateOld(k,7)+strainInc(k,2)
      stateNew(k,8)=stateOld(k,8)+strainInc(k,3)
      stateNew(k,9)=stateOld(k,9)+strainInc(k,4)
      stateNew(k,10)=stateOld(k,10)+strainInc(k,5)
      stateNew(k,11)=stateOld(k,11)+strainInc(k,6)
C      write(*,*)'strainfirst',stateNew
C-------------------------------------------------
      stateNew(k,1)=stateOld(k,1)
      stateNew(k,2)=stateOld(k,2)
      stateNew(k,3)=stateOld(k,3)
      stateNew(k,4)=stateOld(k,4)
      stateNew(k,5)=stateOld(k,5)
      stateNew(k,12)=stateOld(k,12)
      stateNew(k,13)=stateOld(k,13)      
C------------------------------------------------
      dft=zero
      dfc=zero
      dmt=zero
      dmc=zero
      ds=zero
C---------------------------------------------------
      if (stateOld(k,13) .EQ. one)  then
       if (stateOld(k,4) .EQ. one) then
         E2x=E2*0.4
         G12x=G12*0.4
         G23x=G23*0.4
       end if
       if (stateOld(k,3) .EQ. one) then
         E2x=E2*0.2
         G12x=G12*0.2
         G23x=G23*0.2
       end if
      end if
      if (stateOld(k,12) .EQ. one) then
        if (stateOld(k,2) .EQ. one) then
         E1x=E1*0.14
         E2x=E2*0.14
         G12x=G12*0.14
         G13x=G13*0.14
        end if
       if (stateOld(k,1) .EQ. one) then
         E1x=E1*0.07
         E2x=E2*0.07
         G12x=G12*0.07
         G13x=G13*0.07
        end if
       end if
      if (stateOld(k,5) .EQ. one) then
         G12x=G12*zero
         G13x=G13*zero
        end if
      xnu21 = xnu12 * xE2 / xE1
      xnu31 = xnu13 * xE3 / xE1
      xnu32 = xnu23 * xE3 / xE2
c      write(*,*)'goog',f1t,f1c,f2t,f2c,f3t,f3c,f12,
c     *       f23,E1,E2,E3,xnu12,xnu13,xnu23,G12,G23,G13,
c     *        xnu21,xnu31,xnu32
*       Compute terms of stiffness matrix
      gg = one / ( one - xnu12*xnu21 - xnu23*xnu32 - xnu31*xnu13
     *     - two*xnu21*xnu32*xnu13 )
      C11  = E1x * ( one - xnu23*xnu32 ) * gg
      C22  = E2x * ( one - xnu13*xnu31 ) * gg
      C33  = E3x * ( one - xnu12*xnu21 ) * gg
      C12  = E2x * ( xnu12 + xnu32*xnu13 ) * gg
      C13  = E3x * ( xnu13 + xnu12*xnu23 ) * gg
      C23  = E3x * ( xnu23 + xnu21*xnu13 ) * gg
      C21  = E1x * ( xnu21 + xnu31*xnu23 ) * gg
      C31  = E1x * ( xnu31 + xnu21*xnu32 ) * gg
      C32  = E2x * ( xnu32 + xnu12*xnu31 ) * gg
c       write(*,*)'gg',gg,C11,C22,C33,C12,C13,C23
C---------------------------------------------------------------------
      stressNew(k,1)=C11*stateNew(k,6)+C12*stateNew(k,7)
     *   +C13*stateNew(k,8)
      stressNew(k,2)=C21*stateNew(k,6)+C22*stateNew(k,7)
     *   +C23*stateNew(k,8)
      stressNew(k,3)=C31*stateNew(k,6)+C32*stateNew(k,7)
     *   +C33*stateNew(k,8)
      stressNew(k,4)=G12*two*stateNew(k,9)
      stressNew(k,5)=G23*two*stateNew(k,10)
      stressNew(k,6)=G13*two*stateNew(k,11)
C      write(*,*)'stress',stressNew(k,2),stressNew(k,4),stressNew(k,5)
      s11=stressNew(k,1)
      s22=stressNew(k,2)
      s33=stressNew(k,3)
      s12=stressNew(k,4)
      s23=stressNew(k,5)
      s13=stressNew(k,6)
C---------------------------------------------------------------------
      if (s11 .GE. zero) then
        dft=(s11/f1t)**two
      else
        dfc=(s11/f1c)**two
      end if
      if (s22 .GE. zero) then
        dmt=(s22/f2t)**two+(s12/f12)**two
     *     +(s23/f23)**two
      else
        dmc=(s22/(two*f12))**two
     *     +(s22/f2c)*((f2c/(two*f12))**2-one)
     *     +(s12/f12)**two+(s23/f23)**two
      end if
      if (s11 .GE. zero) then
        ds=(s12/f12)**two+(s13/f13)**two
      else
        ds=(s11/f1c)**two
     *       +(s12/f12)**two+(s13/f13)**two
      end if
C      write(*,*)'df',dft,dfc,dmt,dmc,ds
C------------------------------------------------------------------
      if (dft .GE. one) then
        stateNew(k,1)= one
        stateNew(k,12) = one
      end if
      if (dfc .GE. one) then
        stateNew(k,2)= one
        stateNew(k,12) = one
      end if
      if (dmt .GE. one) then
        stateNew(k,3)= one
        stateNew(k,13) = one
      end if
      if (dmc .GE. one) then
        stateNew(k,4)= one
        stateNew(k,13) = one
      end if
      if (ds .GE. one) then
        stateNew(k,5)= one
      end if
C      write(*,*)'statenew',statenew(k,3),statenew(k,4)
  100  continue
      return
      end

qliao

divinelove 发表于 2012-7-13 22:08
上传一下最新的程序，比原来的效果好，但是还是与实验结果会很大，希望大家多交流：
主程序依然贴在下面，
...

      abaqus自带的二维hashin采用基于断裂能的折减方案，你可以参考abaqus analysis user's manual——21.3.3 damge evolution and element removal for fiber-reinforced composites。
      example manual——1.4.6 failure of blunt notched fiber metal lamnates 采用linde准则。采用一种基于应变的折减方案，不需要人为给定折减系数，也不需要测定断裂能，但好像也针对二维问题。

qliao

divinelove 发表于 2012-7-13 22:08
上传一下最新的程序，比原来的效果好，但是还是与实验结果会很大，希望大家多交流：
主程序依然贴在下面，
...

不一定要通过对模量和泊松比进行折减来折减刚度
那两个例子都是直接对刚度矩阵中相应的项进行折减

divinelove

qliao 发表于 2012-7-14 14:24
不一定要通过对模量和泊松比进行折减来折减刚度
那两个例子都是直接对刚度矩阵中相应的项进行折减 ...

嗯，首先谢谢你的答复，
现在做低速冲击这块，有两个方法，一直是直接折减弹性模量，就是我上面的那个方法，还有一个就是你说的，基于损伤力学的思想，我刚刚开始的时候用的是直接折减弹性模量的方法，从现有的资料看，没有一个完整的，而且程序一直有问题，中途就放弃了，改到了用损伤力学的方法，很高兴的说，基本上已经完成了vumat的编写，只是现有的损伤刚度矩阵说法不一，结果也会有不一样，在调试中。另外想，既然很多人都完成了基于弹性模量的方法折减，我当然也想尝试着做，也是我毕业设计的一部分。。。。

greenteasong

divinelove 发表于 2012-7-14 16:04
嗯，首先谢谢你的答复，
现在做低速冲击这块，有两个方法，一直是直接折减弹性模量，就是我上面的那个方 ...

楼主，现在这个问题解决没有啊？

greenteasong

我现在用简单的最大应力准则判断破坏，做准静态分析，破坏之前是纯线弹性的，破坏之后按线性degradation,采用刚度折损的方法
现在测试一个单元没有问题，两个单元也可以，而多个单元测试时发现，一旦进入破坏，就会有一个（几个）单元应变急速增加，导致单元distortion过大，
我同样试了Abaqus自带的材料模型，用显式计算并没有发现这个问题。请问楼主，这个主要是由没有加入粘度系数导致的吗？
或者楼主知道abaqus自带材料模型在下降段是怎么处理的吗，这个问题比较着急，希望得到你尽快回答~谢谢

coobe

谢谢啦 ，感觉在这里长知识特别快

zyhnwpu

divinelove 发表于 2012-6-8 15:02
自己回答一下。程序改用应变的起始破坏准则，可以算下去，具体的原因在考虑中。或许是因为有限元本身事先给 ...

这个是因为冲击的局部区域应力集中较大，导致用应力做判据无法模拟真实情况，而应变不会出现这种情况，它的变化比较平滑。