Abaqus子程序系列：FRIC（定义接触表面的摩擦行为）

abaqus用户子程序fric，在接触分析中，定义复杂的摩擦模型，或者在热力耦合分析中，定义摩擦生热时，潜力巨大。这里先将子程序相关的基础知识，进行了整理。后续会更新基于子程序的相关应用案例。

1.概述：

用户子程序FRIC对应于关键字*FRICTION（定义一个摩擦模型。用于将摩擦特性引入表面接触模型中，来控制接触表面、接触对或连接器单元的切向接触行为。），以及交互界面里的接触属性中切向行为的所有内容（除了用户自定义外，abaqus中可以定义5种类型的摩擦行为（摩擦公式），每个公式中，主要是定义三方面的内容：摩擦因子，剪切应力，弹性滑动（可以恢复的滑动位移））。

用户子程序FRIC：

• 可用于定义接触面间的摩擦行为；

• 当Abaqus中提供的经典库仑摩擦模型的扩展版本限制太严格，或者需要在接触面间定义更复杂的切向应力时，可以使用；

• 当接触属性模型包含用户子程序定义的摩擦时，当接触点闭合时，接触对的从属表面上的节点或者接触单元的积分点会调用子程序；

  每个增量步里的每次迭代，接触对中，从表面上，处于接触闭合状态的节点，会调用子程序。

• 必须提供接触切向行为的完整定义；

• 可以用来更新解相关的状态变量（例如常见的：弹性能量密度，摩擦生热等）；

2.子程序基础介绍

子程序接口

SUBROUTINE FRIC(LM,TAU,DDTDDG,DDTDDP,DSLIP,SED,SFD,
1 DDTDDT,PNEWDT,STATEV,DGAM,TAULM,PRESS,DPRESS,DDPDDH,SLIP,
2 KSTEP,KINC,TIME,DTIME,NOEL,CINAME,SLNAME,MSNAME,NPT,NODE,
3 NPATCH,COORDS,RCOORD,DROT,TEMP,PREDEF,NFDIR,MCRD,NPRED,
4 NSTATV,CHRLNGTH,PROPS,NPROPS)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CINAME,SLNAME,MSNAME
C
DIMENSION TAU(NFDIR),DDTDDG(NFDIR,NFDIR),DDTDDP(NFDIR),
1 DSLIP(NFDIR),DDTDDT(NFDIR,2),STATEV(*),DGAM(NFDIR),
2 TAULM(NFDIR),SLIP(NFDIR),TIME(2),COORDS(MCRD),
1.1.8–1
FRIC
3 RCOORD(MCRD),DROT(2,2),TEMP(2),PREDEF(2,*),PROPS(NPROPS)
user coding to define LM, TAU, DDTDDG, DDTDDP, and, optionally, DSLIP, SED, SFD, DDTDDT, PNEWDT, STATEV
RETURN
END

必须要定义变量：

1. LM

相对滑动标记。仅在确定接触点的接触状态处于closed时，用户子程序FRIC才会被调用；即两种情形：接触压力是正的，接触点在前一次迭代中，接触状态是closed；接触点处于过度闭合状态（干涉），接触点在前一次迭代中处于open状态。

在迭代过程中，LM的值会首先传递到子程序，此时的值为前一个迭代中定义的值。在一个增量开始时，或者如果接触点在前一个迭代中是open状态，这个变量将是否传递进子程序中，取决于前一个增量中的接触条件。如果接触点滑动，LM = 0；如果接触点是粘着的，LM = 1；如果接触点是open状态，LM等于2。

如果这次迭代中，相对滑动是允许的（由于滑移或弹性粘着），设置LM等于0。在这种情况下，子程序中必须将切向摩擦应力定义为相对滑动、界面接触压力p和其他预定义或用户定义的状态变量的函数。此外，子程序必须定义切向摩擦应力相对于接触压力p的倒数(即切向应力的变化与接触压力的变化之间的关系)。

如果这次迭代中，不允许相对滑动，设置LM等于1。界面处的刚性粘着条件，采用拉格朗日乘子方法施加。在这种情况下，不需要更新其他变量。如果LM总是设置为1，会生成一个“非常粗糙”的界面。当使用有限滑动，surface to surface的接触公式时，不建议将LM设为1。

如果摩擦被忽略（假设无摩擦的滑动），在这种情况下，不需要更新其他变量。如果LM总是设置为2，则生成一个“完美光滑”的界面。

也可以根据增量滑动量的信息，以及计算得到的摩擦应力来决定粘着/滑动状态（摩擦力计算得到摩擦应力，然后求得等效摩擦应力，去与临界剪切应力比较）。这些变量也会从abaqus/standard中传递到子程序中。

为了避免通常的摩擦接触问题中出现收敛问题，如果在前一个增量结束时，该接触点已打开，将LM设置为2，并且退出子程序，也就是说，如果Abaqus/Standard 设置 LM=2，当它调用这个子程序时，只需退出这个子程序。

2. 如果LM返回的值是0

数组TAU（NFDIR）

在增量开始时，这个值传入子程序，其值是摩擦切向应力分量，必须在增量结束时更新这个值。

将等效剪切应力与临界剪切应力相比较，可以判断时粘着，还是滑动。

二维数组DDTDDG（NFDIR,NFDIR）

方向b上的，摩擦切向应力相对于方向a上的相对运动的偏导数。

二维数组DDTDDP（NFDIR）

方向a上的摩擦应力对接触压力的偏导数。由于这些项对刚度矩阵产生非对称影响，只有在使用非对称方程求解器时才使用它们。

可以被更新的变量：

DSLIP（NFDIR）

不可恢复滑动位移(滑移)的增量。如果在先前的增量步中，LM是0，这个数组传递进子程序，其值作为前一次迭代中用户定义的值。否则，其值将是0。仅仅当LM的返回值是0时，这个数组才会被更新。

这个数组用于检测迭代之间的滑动反转。输出选项使用它来指示该点是粘着状态还是滑动状态。当一个增量步收敛时， DSLIP(NFDIR)中的值会在累积在SLIP(NFDIR)中，作为塑性应变存储。

如图，一开始不发生滑动（粘着），随着摩擦剪切力增大，到达临界剪切应力（u*p），开始滑动，滑动后，随着滑动速度降为0，摩擦剪切应力会降低。从而由滑动又变成粘着状态。随着变形的进行，粘着后，还会发生反向滑动。

其中t是滑动方向，这是由接触公式决定的（N-S,S-S,小滑动，有限滑动）。U是位移的增量。例如一个圆盘的旋转。U=速度（ω*r）*Dtime（增量步的时间）.

SED

该变量在增量开始时作为弹性能量密度传入，并应在增量结束时更新为弹性能量密度。此变量仅用于输出，对其他解变量没有影响。

SFD

该变量应定义为摩擦耗散增量。如果调用FRIC的接触单元或接触对称使用应力而不是力，单位是每单位面积的能量。对于常规应力分析，该变量仅用于输出，对其他解变量没有影响。在温度-位移耦合和热-电-结构耦合分析中，如果采用间隙生热模型，可将耗散转化为热量（也是fric子程序，比较重要的应用方向）。如果SFD没有定义，可以通过滑移增量DSLP、摩擦切向应力TAU的乘积得到的耗散来计算生热（生热（摩擦耗散、摩擦做功）=摩擦切向力（=摩擦因子*法向压力）*滑动位移）。

DDTDDT（NFDIR，2）

a方向上的摩擦切向应力,对两个表面温度的偏导数。这只对于温度-位移耦合和热电结构耦合单元是需要的，这些单元中，摩擦剪切应力是表面温度的函数。

PNEWDT

建议的新的时间增量步尺寸与现在使用的时间增量步尺寸的比值。如果使用自动时间增量步算法，会允许这个变量。

每次调用FRIC子程序之前，这个值会设置成一个非常大的值。

如果PNEWDT被重新定义，小于1，abaqus/standard将放弃此时的时间增量步尺寸，尝试用更小的时间增量尺寸。为自动时间积分算法提供的建议的新时间增量等于PNEWDT乘以DTIME（此时的时间增量步尺寸），其中使用的PNEWDT��所���调用的用户子程序（该用户子程序允许在此迭代中重新定义PNEWDT值）里定义的最小值。

如果对于这个迭代中，所有调用的子程序中，PNEWDT是一个给定的远大于1的值，并且这个增量步在这个迭代收敛，abaqus/standard会增加时间增量步尺寸。为自动时间积分算法提供的新的时间增量为PNEWDT *DTIME，其中使用的PNEWDT是所有调用的用户子程序（该用户子程序允许在此迭代中重新定义PNEWDT值）里定义的最小值。

如果在分析过程中没有选择自动时间增量算法，PNEWDT大于1.0的值将被忽略，小于1.0的PNEWDT值将导致作业终止。

STATEV（NSTATV）

包含用户定义的解相关状态变量的数组。可以指定可用状态变量的数量。这个数组将在增量步开始时传入，其中包含这些变量的值。如果将任何与解相关的状态变量与摩擦行为一起使用，则必须在增量结束时，在此子程序中更新它们的值。

仅仅为了传递信息进子程序的变量

DGAM（NFDIR）

如果在先前的迭代中，LM设置为0，它的值是这个增量步中的滑动增量。否则，为0。通过与DSLIP(NFDIR)进行比较，可以确定是否在此时，这个点由粘着变成滑动状态，或者在此时是否发生了滑动方向反转。

TAULM（NFDIR）

如果先前的迭代中，LM设置为1，该值是增量步结束时约束应力(等效剪切应力)的当前值。否则，将是0。通过与临界剪应力的比较，可以确定此时这个点是否由粘着变成了滑动。

有些量是对于增量步更新的，有些量是对于增量步中的每次迭代更新的。如上面的DGAM和TAULM是每次迭代中更新的，而DSLIP是增量步中更新的。

PRESS

增量步结束时的接触压力

DPRESS

接触压力的增量

DDPDDH

在软接触的情况下，此时的接触刚度

SLIP（NFDIR）

在增量步开始时，总的不可恢复滑动。这个值是从前面的增量步中，DSLIP(NFDIR)的累计值。

KSTEP

分析步号

KINC

增量步号

TIME(1)

增量步结束时，分析步的时间值

TIME（2）

增量步结束时，总时间的值

DTIME

这个增量步的时间尺寸

NOEL

接触单元的单元编号。如果定义的是接触表面，传递进来的是0

CINAME

与摩擦定义关联的用户定义的表面的名称，左对齐。对于接触单元，它是与摩擦定义相关联的界面定义所给定的单元集名称。如果给界面定义分配了一个可选名称，那么CINAME将作为这个名称传入，左对齐。

SLNAME

从表面名称。如果使用接触单元，传入空格。

NSNAME

主表面名称。如果使用接触单元，传入空格。

NPT

接触单元的积分点数量。如果接触表面被定义，传入0。

NODE

与这个接触点相关的用户定义的全局从节点号（或根据部件实例的组装定义的模型的内部节点编号)。如果使用面对面接触公式，则对应于约束的主要从节点。如果从接触单元调用，则作为零传递。

COORDS(MCRD)

包含这个点，此时坐标的数组

RCOORD(MCRD)

如果主表面被定义为刚性表面，则传递此数组，其中包含刚性表面上对应点的当前位置和方向的坐标。

DROT(2,2)

旋转增量矩阵。对于与三维刚性表面的接触，这个矩阵表示相对于刚性表面的表面方向的增量旋转。这样做是为了在这个子程序中，将向量值或张力值状态变量进行正确的旋转。在调用FRIC之前，应力和滑动分量已经旋转了这个量。该矩阵作为单位矩阵传入二维和轴对称接触问题。

TEMP（2）

从节点和对应的主表面，此时的温度

PREDEF（2，NPRED）

一个数组，包含当前增量步结束时的所有用户指定场变量的一对值(分析开始时的初始值和分析期间的当前值)。从接触对调用FRIC，一对值中的第一个值对应于从节点，第二个值对应于主表面上最近的点。如果从一个大滑动接触单元调用FRIC，则PREDEF(1,NPRED)对应于该单元积分点处的值，PFREDEF(2,NPRED)对应于对应表面上的最近点。如果从一个小滑动接触单元调用FRIC，则PREDEF(1,NPRED)对应于第一侧积分点处的值，PFREDEF(2,NPRED)对应于该单元对应面上积分点处的值。

NFDIR

摩擦方向的数量

MCRD

接触点处坐标方向的数量

NPRED

预定义场变量的数量

NSTATV

用户定义状态变量的数量

CHRLNGTH

接触表面的面特征尺寸，可以用于定义最大允许的弹性滑动

PROPS（NPROPS）

用户定义的属性值数组，用来定义接触表面间的摩擦行为

NPROPS

用户定义的，与摩擦模型相关的属性值的数量