最近博主因为科研的需要，接触到abaqus子程序的二次开发，其中的vuinter是博主主要使用的子程序，它可以定义两表面之间的相互作用，比如力和热，功能比较强大。但可惜的是vuinter才出现不久，网上的资料甚少，所以博主在用的时候遇到了很多坑。因此开个专题来讲这个子程序。由于博主也是第一次使用，难免会有错误之处，请谅解，也欢迎大家一起讨论。官方手册放着了[2017abaqus用户手册](SIMULIA Online User Assistance (mit.edu))

VUINTER

(一)vuinter子程序介绍

VUINTER是定义表面相互作用的一个子程序，可自定义力和热的相互作用，比如物体表面摩擦，芯片之间的接触等等,能实现像VFRIC、VFRICTION等其他子程序的功能，但比它们功能更强大，值得一提的是根据官方手册上介绍，具有类似作用，比如定义表面作用的各种子程序之间不能同时调用，但如果一个是改变材料属性的vumat，另一个是vuinter的话就可以。

abaqus子程序是用Fortran语言进行编译，难度并不大，但要非常注意语言的格式，标点符号，中英文等等都会导致数据出错，得到错误的仿真结果。首先介绍vuinter使用的界面：

这个界面并不是不能动的，很多人在刚开始学的时候往往不敢动，但其实只要按照格式，会发现改变界面会方便许多，比如定义一个用户自定义变量或者矩阵。

(二)vuinter子程序的stess变量

由于博主并没有使用全部的变量，因此主要介绍关于力的stress变量，而且这个变量很坑

stress(nDir, nSlvNod)

这个表示上一个时间增量上表面的节点力（注意这里stress并不是应力的意思，这里是一个坑），它是必须要被定义的变量。我们可以在上面的conventions for stress中更深入了解这一变量。正的并且nDir为1的stress代表进入表面的压强（与局部法向方向相反），模拟压缩等效果，负的代表与局部法向方向相同的压强，可模拟粘附力等效果。nDir为2或3（二维只有2）代表给表面节点施加剪力（！！！注意是剪力，单位是N）。如果想亲自验证一下，可更改官方的例子。选择三维施加法向位移的例子

vuinter3d_n.inp

在Fortran程序中删除多余的编程，替换stress(i_str_S11,kSlv) = 50.d0

    subroutine vuinter( sfd, scd, spd, svd, 
*     stress, fluxSlv, fluxMst, sed, statev,
*     kStep, kInc, nFacNod, nSlvNod, nMstNod, nSurfDir,
*     nDir, nUSdv, nProps, NumTemp, NumExfv, numDefTfv,
*     jSlvUid, jMstUid, jConMstid, timStep, timGlb,
*     dTimCur, surfInt, surfSlv, surfMst,
*     rdisp, drdisp, drot, stiffDflt, condDflt,
*     shape, coordSlv, coordMst, alocaldir, props,
*     areaSlv, tempSlv, dtempSlv, exfvSlv, dexfvSlv, 
*     tempMst, dtempMst, exfvMst, dexfvMst )

  include 'vaba_param.inc'

 character*80 surfInt, surfSlv, surfMst

  dimension props(nProps), statev(nUSdv,nSlvNod), 
*     drot(2,2,nSlvNod), sed(nSlvNod), sfd(nSlvNod),
*     scd(nSlvNod), spd(nSlvNod), svd(nSlvNod),
*     rdisp(nDir,nSlvNod), drdisp(nDir,nSlvNod), 
*     stress(nDir,nSlvNod), fluxSlv(nSlvNod),
*     fluxMst(nSlvNod), areaSlv(nSlvNod),
*     stiffDflt(nSlvNod), condDflt(nSlvNod),
*     alocaldir(nDir,nDir,nSlvNod), shape(nFacNod,nSlvNod),
*     coordSlv(nDir,nSlvNod), coordMst(nDir,nMstNod), 
*     jSlvUid(nSlvNod), jMstUid(nMstNod),
*     jConMstid(nFacNod,nSlvNod), tempSlv(nSlvNod),
*     dtempSlv(nSlvNod), exfvSlv(NumExfv,nSlvNod),
*     dexfvSlv(NumExfv,nSlvNod), tempMst(nSlvNod),
*     dtempMst(nSlvNod), exfvMst(NumExfv,nSlvNod),
*     dexfvMst(NumExfv,nSlvNod)
 
 parameter ( i_prp_GapCutOff = 1,
*            i_prp_YoungsModulus = 2,
*            i_prp_PoissonsRatio = 3,
*            i_prp_InitYield = 4,
*            i_prp_HardenMod = 5,
*            i_prp_ThickInter = 6,
*            i_prp_IfcCond = 7 )

 parameter ( i_usv_CompletedInit = 1,
*            i_usv_BondStatus = 2,
*            i_usv_CurYield = 3 )

 parameter ( i_str_S11 = 1, i_str_S12 = 2, i_str_S13 = 3 )

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

 if ( statev(i_usv_CompletedInit,1) .eq. zero ) then

    do kSlv = 1, nSlvNod
       statev(i_usv_CompletedInit,kSlv) = one
       gapInit = -rDisp(1,kSlv)
       if ( gapInit .le. props(i_prp_GapCutOff) ) then
          statev(i_usv_BondStatus,kSlv) = one
          statev(i_usv_CurYield,kSlv) = props(i_prp_InitYield)
       end if
     end do
 end if
         xMu = half * props(i_prp_YoungsModulus) /
*                (one + props(i_prp_PoissonsRatio))

  do kSlv = 1, nSlvNod
     if ( statev(i_usv_BondStatus,kSlv) .eq. one ) then
        stress(i_str_S11,kSlv) = 50.d0

     else

       stress(i_str_S11,kSlv) = zero

     end if
  end do

 return
 end