自己开发各项同性弹塑性模型出现的问题，请高手帮忙！

liyunliang_zi

初始状态如下图所示：

计算后如下：

同样的接触定义，用dyna自带的弹塑性模型算起来没问题，而用自己开发的就不行。
把网上的给出的例子替换自己编写的用户自程序也是这个问题。

为了验证是否是接触的问题，我尝试把：
       scle=0
       if (ak2.gt.0) then
           scle=1.0
       endif
改为：
       scle=0
c       if (ak2.gt.0) then
c           scle=1.0
c      endif
即材料永远都不会屈服，算起来就没有错误了！

实在找不到错在哪了，我想监视出错的那个单元，但不知道怎么选定那个单元！

我试着调试，
找不到好办法，我用以下方式调试，加入：
write(*,*) scle, aj1, ak fac1, depi, hsv(1)
程序运行太慢了，后来加入下面代码：
hsv(3)=hsv(3)+1
if (hsv.gt.100) then
   write(*,*) scle, aj1, ak fac1, depi, hsv(1)
   hsv(3)=0
endif

但是，不知道输出的是哪个单元的信息，而且找不出错误的原因。

请各位大侠指点！
代码如下：


请各位大侠指点！

liyunliang_zi

另外，单单元调试没有问题，与dyna自带弹塑性模型计算结果一致。
关键是，要实时输出畸变单元内的应力、应变等参数，以便于调试！
但，就是找不到如何指定畸变单元，把全部单元结果输出太慢了，而且不知道，到底是哪个单元产生畸变。

liyunliang_zi

把所有想要监视的变量赋给hsv（i），然后通过后处理软件，找到指定单元的监视量。
问题已解决。
不过，还是不知道怎么去获取单元的内部ID!
lpfinvf()函数不会用。