Abaqus二次开发Vfric摩擦子程序，接触行为自定义

Abaqus/Explicit analysis exited with an error，msg里刷出一行“User subroutine VFRIC is currently not supported for use with the general contact algorithm”，VFRIC白写了，模型用的是通用接触，这个子程序根本不认。VFRIC只能挂在接触对（Contact Pair）上，通用接触用不了，这一点在inp里就得配好，Interaction模块里的General Contact勾选去掉，老老实实建接触对。

VFRIC的接口参数不算少，写起来容易漏。subroutine vfric(fTangential, statev, kStep, kInc, nContact, nFacNod, nSlvNod, nMstNod, nFricDir, nDir, nStateVar, nProps, nTemp, nPred, numDefTrv, jSlvUid, jMstUid, jConslvid, jConMstid, timStep, timGlb, dTimCur, surfInt, surfSlv, surfMst, lContType, dSlipFric, fStickForce, ffangPrev, fNormal, frictionWork, shape, coordSlv, coordMst, dirCosSl, dirCosN, props, areaSlv, tempSlv, preDefSlv, tempMst, preDefMst)，include 'vaba_param.inc'这行必须有，跟Standard侧的ABA_PARAM.INC是两套东西，别搞混。dimension声明里fTangential的维度是(nFricDir, nContact)，statev的维度是(nStateVar, nSlvNod)，dSlipFric是(nDir, nContact)，这几个维度在循环里用错了就是越界，编译不报错但运行直接崩。

要用户自己填的变量不多，核心就两个。fTangential(nFricDir, nContact)是摩擦力分量，三维问题里nFricDir=2，fTangential(1,:)是t1方向的分量，fTangential(2,:)是t2方向的分量，t1和t2都在切平面内，t1沿滑移方向定义，t2=n×t1。另一个是statev，用来存跨增量步的状态量，摩擦系数、累积滑移、磨损量这些都可以塞进去。fStickForce是ABAQUS算好的“粘滞力”，调用VFRIC之前就已经在接触算法里算出来了，你给定的摩擦力大小不应该超过这个值，超了求解器会震荡。fNormal是接触点的法向压力，算库仑摩擦直接拿它乘系数就行，但注意它是在调用VFRIC之前由接触算法提供的，不是你自己算的。

有一次做一个土与管线相互作用的模型，想用VFRIC实现一个跟累积滑移量相关的动态摩擦系数，滑移越大摩擦系数越小。模型是三维的，接触对建好了，inp里也写了*FRICTION, USER，子程序编过了，一跑就崩。当时是夏天，办公室空调排水管堵了，天花板上滴答滴答往下滴水，保洁阿姨拿了个桶放在走道里接着，声音很规律。第一反应是VFRIC没被调用，在子程序开头加了write(*,*)打不动，因为Explicit不往屏幕输出，改成写文件，跑完一看调用日志是有的。又怀疑是维度搞错了，检查了一遍fTangential用的是fTangential(1,nCnt)和fTangential(2,nCnt)，循环变量是nCnt从1到nContact，没毛病。翻了半天发现statev的索引写成了statev(nCnt, i)，但声明里是statev(nStateVar, nSlvNod)，第一个维度才是状态变量编号，第二个是从面节点编号，反了。nSlvNod是几百个，nStateVar只设了5个，statev(300, 2)直接踩到内存外面去了。改过来之后跑通了，摩擦系数随滑移的衰减规律跟试验数据基本吻合，峰值摩擦力差了不到8%。

编译配置上，VFRIC是Explicit侧的，生成的是explicitU.dll，跟UMAT生成的standardU.dll不是一回事。命令行abaqus job=Job-1 user=vfric_test.for cpus=4 interactive，如果模型里同时还要用VUMAT或者DLOAD，所有子程序必须写在同一个.for文件里，分开写两个文件只认一个，另一个直接不编译。子程序之间要共享数据的话用Fortran Module或者COMMON块，Module干净一些。VFRIC里读VUMAT算出来的损伤变量，可以用COMMON块传，但注意Explicit多线程下COMMON块会互相覆盖，cpus=4的时候跑出来的结果可能跟cpus=1不一样，调通了之后再开多核。

fTangential的符号约定跟直觉有点反。t1沿滑移方向定义，摩擦力阻挠滑移，所以fTangential(1,nCnt)通常是负值，在0到负的fStickForce之间取值。fStickForce一般是正的，代表需要多大的力才能完全阻止滑移。如果算出来的fTangential(1)是正的，说明你在给系统加能量，摩擦在帮忙推物体，这个方向就错了。调试的时候把每个接触点的fStickForce、dSlipFric和算出来的fTangential一起写文件，看看符号和量级对不对。

验证方法用一个最简单的模型，刚性平面上的弹性方块，加初始速度让它滑，VFRIC里实现固定摩擦系数，跟ABAQUS内置的Penalty摩擦对一下，摩擦力和滑动距离应该基本一致。内置的对上了再换变摩擦系数的模型，逐步加复杂度。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删