请教各位大侠ＶＵＭＡＴ材料子程序的问题，小女子这厢有礼啦！

dwq

一个描述变化弹性模量的材料子程序在调试时没有错误，用在单个单元模型上也没有错误，可是用于我的矩形管弯曲模型上就会出现如下错误提示：
No transverse shear stiffnesses are prescribed for element set ASSEMBLY_TUBE-1__PICKEDSET8. Use *TRANSVERSE SHEAR STIFFNESS
Abaqus/Explicit Packager exited with an error - Please see the  status file for possible error messages if the file exists.
版本：ｖｉｓｕａｌ　ｓｔｕｄｉｏ２００５，ｉｎｔｅｌ　ｆｏｒｔｒａｎ　９．１，ａｂａｑｕｓ６．８．１
路过的各位大侠，还请拔刀相助啊！

dwq

怎么没有高手路过呀？

shawn2008

单元的某些刚度是由材料的弹性性能计算得到的。当采用umat或vumat时，由于abaqus不能识别材料的弹性性能参数，所以，必须在inp文件中定义这些刚度。 TRANSVERSE SHEAR STIFFNESS就属于这种情况。

定义TRANSVERSE SHEAR STIFFNESS的具体格式和数值可以查看abaqus的document，里面有详细的解释。

dwq

真是好人啊！好人一生平安呦！再次请教版主,上个问题是不是由于壳单元造成的？我把我的子程序用到实体模型上就会出现如下错误：
The executable C:\SIMULIA\6.8-1\exec\explicit.exe aborted with system error code 144. Please check the .dat, .msg, and .sta files for error messages if the files exist.  If there are no error messages and you cannot resolve the problem, please run the command "abaqus job=support information=support" to report and save your system information.  Use the same command to run Abaqus that you used when the problem occurred.  Please contact your local Abaqus support office and send them the input file, the file support.log which you just created, the executable name, and the error code.
这又是为什么呢？
还请再赐教啊··········

shawn2008

壳单元和实体单元的ndi和nshr不一样的，除非你在子程序中考虑了这些区别，否则不能通用。

比如对壳单元，张量分量一般为11, 22, 12；但对实体单元一般张量分量为11, 22, 33, 12, 23, 13。 虽然对一般金属，弹性性能都假设满足虎克定律，但在这两种分量表达下，弹性矩阵的形式是不一样的。

因此，在子程序的开发时要考虑的细节问题很多的。

如果是做冲压成形分析，像S4R单元就必须要定义TRANSVERSE SHEAR STIFFNESS
，在shell section关键词后面定义。
*Shell Section, elset=_PickedSet8, material=al
1.3, 15
*Transverse Shear
28508.8, 28508.8, 0.

aerojeff

真是好人啊！好人一生平安呦！再次请教版主,上个问题是不是由于壳单元造成的？我把我的子程序用到实体模型上就会出现如下错误：
The executable C:\SIMULIA\6.8-1\exec\explicit.exe aborted with system error cod ...
dwq 发表于 2009-11-6 20:37

这个错误是程序中数据溢出所导致的，LZ好好地调调程序吧，特别注意一下分母，根号下的那些变量

dwq

谢谢Shawn2008大侠啦！我再改一下试试。不过我不太明白的是在ABAQUS的自带例子中，Isotropic Kinematic Hardening,用于plane strain中，也没有定义transverse shear stiffness,但是却能用在我的模型计算中，这是为什么呢？下面就是我的程序，好多天还是没有调好，还请大侠一针见血，指点迷津呀......
subroutine vumat(
C Read only (unmodifiable)variables -
     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 (modifiable) variables -
     7  stressNew,stateNew,enerInternNew,enerInelasNew)
      include 'vaba_param.inc'
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(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),
     8  defgradNew(nblock,ndir+nshr+nshr),fieldNew(nblock,nfieldv),
     9  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     1  enerInternNew(nblock), enerInelasNew(nblock)
C
      character*80 cmname
C
         parameter(zero=0.,one=1.,two=2.,three=3.,EI=70.,A=15.,ES=55.,
     1  third=one/three,half=.5,twoThirds=two/three,
     2  threeHalfs=1.5)
C
      xnu=props(1)
      hn =props(2)
      hk =props(3)
      e  =props(4)

C
      if(stateold(i,1).le.0.006)then
      e=EI-A*stateold(i,1)/0.006
      HLMT=e*xnu/(one+xnu)/(one-two*xnu)
      HMU=e/two/(one+xnu)
      thremu=three*HMU
      stt=sqrt(threeHalfs)
      stt1=sqrt(twoThirds)
      stwo=sqrt(two)
      else if(stateold(i,1).gt.0.006)then
      e=ES
      HLMT=e*xnu/(one+xnu)/(one-two*xnu)
      HMU=e/two/(one+xnu)
      thremu=three*HMU
      stt=sqrt(threeHalfs)
      stt1=sqrt(twoThirds)
      stwo=sqrt(two)
      endif
C
      do 100 i = 1,nblock


C   compute  the  tiji  strain
        dstrain=strainInc(i,1)+strainInc(i,2)+strainInc(i,3)
C   compute  the  trial stress
      Tstress1=stressOld(i,1)+HLMT*dstrain+two*HMU*strainInc(i,1)
      Tstress2=stressOld(i,2)+HLMT*dstrain+two*HMU*strainInc(i,2)
      Tstress3=stressOld(i,3)+HLMT*dstrain+two*HMU*strainInc(i,3)
      Tstress4=stressOld(i,4)+two*HMU*strainInc(i,4)
     
C   compute the hydraumatic stress
        Hstress=third*(Tstress1+Tstress2+Tstress3)
C   compute the bias stress
        PTstress1=Tstress1-Hstress
        PTstress2=Tstress2-Hstress
        PTstress3=Tstress3-Hstress
        PTstress4=Tstress4
       
C   computer the longth of the bias stress
        EPTstress=sqrt(PTstress1**2+PTstress2**2+PTstress3**2+
        1  two*PTstress4**2)
C   computer the equal stress
        EQstress=stt*EPTstress
C   computer the sequent yeild stress
        HY=hk*stateold(i,1)**hn
C   computer the square longth of the sequent yeild stress
        HKF=HY**2
C Check for yield by determining the factor for plasticity,
C zero for elastic, one for yield
      
        facyld = zero
        if(  EQstress**2-HKF.ge.zero ) facyld = one
C
C Add a protective addition factor to prevent a divide by zero
C when dsmag is zero. If dsmag is zero, we will not have exceeded
C the yield stress and facyld will be zero.
       EPTstress=EPTstress+(facyld-zero)
C   computer the hard
         H=hn*hk*stateold(i,1)**(hn-1)       
C   computer the unit normal direction of the yeild surface
        Q1=PTstress1/EPTstress
        Q2=PTstress2/EPTstress
        Q3=PTstress3/EPTstress
        Q4=stwo*PTstress4/EPTstress
       
C compute the dagama
      
      term=one/(two*HMU*(one+H/thremu))
      dagama=facyld*term*(EPTstress-stt1*HY)
C   update the stress
        stressNew(i,1)=Tstress1-two*HMU*dagama*Q1
      stressNew(i,2)=Tstress2-two*HMU*dagama*Q2
      stressNew(i,3)=Tstress3-two*HMU*dagama*Q3
      stressNew(i,4)=Tstress4-two*HMU*dagama*Q4
      
C   update the state
      stateNew(i,1)=stateold(i,1)+stt1*dagama
c==============================================
C Update the specific internal energy -
      stressPower=half*(
     1  (stressOld(i,1)+stressNew(i,1))*strainInc(i,1)
     1  +(stressOld(i,2)+stressNew(i,2))*strainInc(i,2)
     1  +(stressOld(i,3)+stressNew(i,3))*strainInc(i,3)
     1  +two*(stressOld(i,4)+stressNew(i,4))*strainInc(i,4))
C
      enerInternNew(i)=enerInternOld(i)
     1  +stressPower/density(i)
C
C Update the dissipated inelastic specific energy -
      smean=third*(stressNew(i,1)+stressNew(i,2)
     1  +stressNew(i,3))
      equivStress=sqrt(threeHalfs*
     1  ((stressNew(i,1)-smean)**2
     1  +(stressNew(i,2)-smean)**2
     1  +(stressNew(i,3)-smean)**2
     1  +two*stressNew(i,4)**2))
C  computer the operator
        deqps=stt1*dagama
      plasticWorkInc=equivStress*deqps
      enerInelasNew(i)=enerInelasOld(i)
     1  +plasticWorkInc/density(i)
       
C  
  100 continue

      return
      end

shawn2008

7# dwq

横向剪切刚度定义的是横向剪切应变与剪应力之间的关系。

在plane strain状态，平面内的应力分量是\sigma_{11},  \sigma_{22} 和 \sigma_{12}，这3个分量是通过本构模型进行更新的。 而且，在plane strain条件下， 沿平面外方向上没有剪切变形。 所以，没有定义横向剪切刚度也可以运算。

但在plane stress状态， 平面内的应力分量是\sigma_{11},  \sigma_{22} 和 \sigma_{12}，这3个也是通过本构模型进行更新的，而在平面外方向上，如果有\epsilon_{23}和\epsilon_{31}存在时，这时候abaqus是不能计算\sigma_{23}和\sigma_{31}的。在这种情况下，abaqus需要用transverse shear stiffness来更新\sigma_{23}和\sigma_{31}。

3D状态下的分析，也和这个类似。

关于代码，慢慢调试，不要着急。

dwq

不愧是版主啊！太感谢您啦，我自己再试试···

dwq

还想请教版主，我的子程序可以用于弯曲计算了，可是当我把.odb文件导入回弹模型中时，却出现
User subroutine umat is missing

Abaqus/Standard Analysis exited with an error - Please see the  message file for possible error messages if the file exists.这个错误。弯曲计算用得是显示算法，回弹用隐式，难道我还要编一个用于计算回弹的UMAT吗？为什么需要UMAT呢？还请赐教啊......

shawn2008

当然需要编写用于隐式计算的umat子程序了。  因为显式计算时你用的是vumat，保存在odb中的应力和内变量都是与这个本构模型相对应的。 为了计算的一致到隐式计算时肯定要相应的umat才能保证这些量是连续的。  当然，你也可以重新指定材料的性能，在回弹时直接用纯弹性卸载。 但这个又有什么意义呢？ 精度是肯定不能保证了，但还好有限元模拟有很多参数可以调整，随便调整其他参数，总会有参数使得模拟结果与实验结果接近。  有了成果，然后论文就发表了，学位就拿到了，项目就交差了。

呵呵。

dwq

回弹过程不就是纯弹性变形嘛？变形中的弹性部分回复造成了回弹呀，我不是很明白回弹中的UMAT的作用，是用它来做弹性应力更新嘛？“你也可以重新指定材料的性能，在回弹时直接用纯弹性卸载。”是不是说我可以在Property中重新定义Elasticity和Plasticity，而不用UMAT？请版主指点！

dwq

还有就是我用VUMAT算出来的结果和实际相差很大，实在不明白为什么，哎.....
是不是子程序能算题不代表它是对的？

shawn2008

1. 子程序能运行，只能代表它能运行了，是否正确，需要进行很多种情况下的one-element实验，如果模拟结果与理论结果一致，才能说这个子程序基本正确，不过在应用于实际模拟时，有可能还会发现一些问题。 所以，一个正确的健壮的子程序，是经过很多次的验证才行的。 写出代码当然不意味着正确。

2. 如果说一个单元的卸载，我们可以确认它是纯弹性的，但在实际应用中，由于很多个单元间的几何上的约束，在回弹过程中，有可能一些单元发生了塑性变形，所以，单纯用纯弹性计算回弹是不太正确的。你可以在很多论文中发现这样的句子。

For efficiency, forming was simulated using the explicit code, ABAQUS/
explicit, with the user-defined subroutine VUMAT, while spring-back was simulated
using the implicit code, ABAQUS/standard, with the user-defined subroutine
UMAT, after forming simulation results were imported from the explicit code to the
implicit code.

dwq

版主啊，如果我只想用VUMAT，不想用UMAT了，能怎么算回弹呢？我实在是不想编程了....呵呵，这东东太累了，？“你也可以重新指定材料的性能，在回弹时直接用纯弹性卸载。”这句话有什么玄机呢？明示啊....

koukou123

达人真是好，版主很热情！赞一个～

rp_no

赞版主，很好很强大

masiwenaij@126

版主都好牛，学习了。

wuxiaotang

资料收集中，学习了。

hellaven

赞 lz 和bz