Abaqus二次开发：Vfric摩擦子程序

冬生亦东生

已于 2023-09-04 18:02:02 修改

阅读量1.3k

点赞数 3

分类专栏： Abaqus子程序开发文章标签：算法

于 2023-08-06 23:22:33 首次发布

本文链接：https://blog.csdn.net/qq_46417931/article/details/132137402

版权

Abaqus子程序开发专栏收录该内容

12 篇文章 15 订阅

订阅专栏

0 简介

（1）用于定义接触对表面之间的摩擦行为；

（2）当经典库仑摩擦模型限制性太强且需要定义更复杂的接触面之间剪切传递时，可以使用；

（3）必须提供接触面之间剪切相互作用的完整定义；

（4）可以使用和更新与求解相关的状态变量；

（5）不能与软化的切向表面行为结合使用；

（6）不能与通用接触算法一起使用。

主要内容

1.调用方法

2.用户子程序接口

3.要定义的变量

4.传递信息的变量

5.案例讲解

1.调用方法

进入相互作用模块进行相互作用属性设置即可。

2.用户子程序接口

         subroutine vfric(
C Write only - 
     1 fTangential, 
C Read/Write - 
     2 statev,
C Read  only - 
     3 kStep, kInc, nContact, nFacNod, nSecNod, nMainNod,
     4 nFricDir, nDir, nStateVar, nProps, nTemp, nPred, numDefTfv,
     5 jSecUid,jMainUid, jConSecid, jConMainid, timStep, timGlb,
     6 dTimCur, surfInt, surfSec, surfMain, lContType,
     7 dSlipFric, fStickForce, fTangPrev, fNormal, frictionWork,
     8 shape, coordSec, coordMain, dirCosSl, dircosN, props,
     9 areaSec, tempSec, preDefSec, tempMain, preDefMain) 
C
      include `vaba_param.inc'
C
      character*80 surfInt, surfSec, surfMain
C
      dimension props(nProps), statev(nStateVar,nSecNod), 
     1 dSlipFric(nDir,nContact),
     2 fTangential(nFricDir,nContact),
     3 fTangPrev(nDir,nContact),
     4 fStickForce(nContact), areaSec(nSecNod),
     5 fNormal(nContact), shape(nFacNod,nContact),
     6 coordSec(nDir,nSecNod), coordMain(nDir,nMainNod), 
     7 dirCosSl(nDir,nContact), dircosN(nDir,nContact),
     8 jSecUid(nSecNod), jMainUid(nMainNod),
     9 jConSecid(nContact), jConMainid(nFacNod,nContact)
     1 tempSec(nContact), preDefSec(nContact,nPred),
     2 tempMain(numDefTfv), preDefMain(numDefTfv,nPred)

      user coding to define fTangential
      and, optionally, statev

      return
      end

3.要定义的变量

fTangential(nFricDir, nContact) : 该数组必须更新为局部切线方向上所有接触点的摩擦力分量的当前值。重置之前此数组将为零（无摩擦力）。三维问题中，nFricDir为2，fTangential(1, nContact)为摩擦力在t1方向分量；fTangential(2, nContact)为摩擦力在t2方向分量。该变量是需要用户编码定义的。t1和n关于全局坐标系下的余弦，分别在dirCosT1和dirCosN中提供。在零增量滑移的情况下,即ds=0,(我们选择t1为与n正交的任意方向）

4.传递信息的变量

nContact ：当前增量步处于接触状态的从面节点数

nSecNod ：从面节点数

nFricDir ：二维问题中为1；三维问题接触点处的切线方向数，即t1和t2，故为2，

nDir ：接触点处的坐标方向数。即滑移量ds是三维矢量还是二维矢量。（在三维模型中，如果接触副中的表面是二维分析刚性表面或由二维元素形成，则nDir将为两个。）

jSecUid(nSecNod) :此数组列出了从面上节点的用户定义的全局节点号。

jMainUid(nMainNod)：此数组列出了主曲面上节点的用户定义的全局节点号。如果主曲面是分析刚性曲面，则该阵列将作为虚拟阵列传入。

jConSecid(nContact)：此数组列出了处于接触的从面节点的曲面节点号。

dtimCur : 当前增量步的起始时刻到当前增量步的结束时刻，即tcuren 到 tcurent+dt的时间间隔。tcurent为当前增量步的结束时刻。

dSlipFric(nDir, nContact) ：该数组包含当前局部坐标系中每个接触点在当前时间增量期间的增量摩擦滑移。这些增量滑移对应于从 tcurent-dt 到tcurent 的时间增量中的切向运动，即当前增量步产生的增量滑移。

fStickForce(nContact) : 该数组包含在每个接触点强制执行粘滞条件所需的摩擦力大小。对于运动接触，该力对应于无滑移；对于罚接触，该力取决于先前的摩擦力、罚刚度值和先前的增量滑移。惩罚刚度是自动指定的。有时，在恢复与惩罚方法相关的弹性滑动期间，粘着力将被指定为负值。这里把它理解成最大静摩擦力就行。

fNormal(nContact) ：此数组包含当前时间增量结束时应用的接触点法向力的大小（在t=tcurent时刻）

5.案例讲解

这个案例主要是根据部件信息选择性的施加摩擦力，大家主要学下代码的含义和用法即可，然后根据自己的问题自行发挥。在学习这个案例之前大家最好看看我的上一篇笔记部件、节点信息公用程序

C
C User subroutine VFRIC 
      subroutine vfric (
C Write only - 
     *     fTangential, 
C Read/Write - 
     *     statev,
C Read only - 
     *     kStep, kInc, nContact, nFacNod, nSecNod, nMainNod,
     *     nFricDir, nDir, nStateVar, nProps, nTemp, nPred, numDefTfv, 
     *     jSecUid, jMainUid, jConSecid, jConMainid, timStep, timGlb,
     *     dTimCur, surfInt, surfSec, surfMast, lContType,
     *     dSlipFric, fStickForce, fTangPrev, fNormal, frictionWork,
     *     shape, coordSec, coordMain, dircosSec, dircosN, props,
     *     areaSec, tempSec, preDefSec, tempMain, preDefMain )
C
      include 'vaba_param.inc'
C
      dimension props(nProps), statev(*), 
     1     dSlipFric(nDir,nContact),
     2     fTangential(nFricDir,nContact),
     3     fTangPrev(nDir,nContact),
     4     fStickForce(nContact), areaSec(nSecNod),
     5     fNormal(nContact), shape(nFacNod,nContact), 
     6     coordSec(nDir,nSecNod), coordMain(nDir,nMainNod), 
     7     dircosSec(nDir,nContact), dircosN(nDir,nContact),
     8     jSecUid(nSecNod), jMainUid(nMainNod),
     9     jConSecid(nContact), jConMainid(nFacNod,nContact),
     1     tempSec(nContact), tempMain(numDefTfv),
     2     preDefSec(nContact, nPred),
     3     preDefMain(numDefTfv, nPred)
C
      character*80 surfInt, surfSec, surfMast
      character*80 cpname
      parameter ( j_node = 0, zero = 0.d0 )
*
      jrcd = 0
      cpname = ' '
      
C     摩擦系数      
      xMu = props(1)
      do kcon = 1, ncontact
         locnum = 0
         
C     jConSecid(kcon) 是处于接触状态从面节点编号(这个节点编号是局部的)，jSecUid(jConSecid(kcon))由局部节点编号得到对应的全局节点编号             
         jusernode = jSecUid(jConSecid(kcon))

C     调用vgetpartinfo公用程序，根据接触点确定部件相关信息       
         
         call vgetpartinfo(jusernode, j_node, cpname, locnum, jrcd)
         if (cpname(1:5).eq.'BLOCK' .and. 
     1        (locnum.eq.101 .or. locnum.eq.102)) then
C   对于二维面摩擦力有1个分量，而接触力是2个方向     
            if ( nDir .eq. 2 ) then
C     法向力            
               fn = fNormal(kcon)
C     最大静摩擦力（实际是粘滞力）               
               fs = fStickForce(kcon)
C     求摩擦力（固定的写法）               
               ft = min ( xMu * fn, fs )               
               fTangential(1,kcon) = -ft
               
C     对于三维面摩擦力有2个分量，而接触力是3个方向                
            else if ( nDir .eq. 3 ) then
               fn = fNormal(kcon) 
               fs = fStickForce(kcon)
               ft = min ( xMu * fn, fs )
               fTangential(1,kcon) = -ft
               fTangential(2,kcon) = zero
            end if
         end if
      end do
*     
      return
      end