Abaqus二次开发之Fric子程序

1、主要功能

主要是建立界面剪应力与剪切位移的关系

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

当Abaqus中提供的经典库仑摩擦模型的扩展版本限制过多，并且需要更复杂的接触面之间剪切作用传递的定义时，可以使用；即可以定义库伦摩擦以外的复杂摩擦

当接触相互作用属性模型包含用户子程序定义的摩擦力，将被调用在接触对从面上的点和接触单元的积分点（仅当接触点闭合时上

必须提供接触面之间剪切相互作用的完整定义；可以使用和更新与解决方案相关的状态变量。

2、需要被定义的变量

LM ：相对运动标志。仅当确定触点闭合时，才会调用用户子程序FRIC；也就是说，如果接触压力为正（接触点在前一次迭代中闭合）或如果接触点过度闭合（接触点是在前一个迭代中打开的）。

迭代期间，LM作为上一次迭代期间定义的值传递给子程序。在增量开始时，或者如果接触点在上一次迭代期间打开，则此变量将根据上一次增量中的接触条件传递到程序中。如果接触点滑动，LM等于0；如果接触点粘连，LM等于1；如果触点打开，LM等于2。

如果允许相对运动（由于滑动或弹性粘滞），则将LM设置为0。在这种情况下，子程序必须指定摩擦应力τ1、（三维有τ2）作为相对滑动运动γ1,（三维有γ2）、界面压力p、预定义或者状态变量的函数。此外，子程序必须定义摩擦应力相对于γ1,γ2和p的派生（导数等）,例如，在各向同性弹性粘附的情况下，∂τ1/∂γ1=∂τ2/∂γ2=kelas,  ∂τ1/∂γ1=∂τ2/∂γ2=0 ,kelas是界面的弹性刚度。

如果不允许相对运动，则将LM设置为1；通过拉格朗日乘子法实现界面处的刚性粘附条件。在这种情况下，无需更新其他变量。如果LM始终设置为1，则会创建一个“完全粗糙”的界面。当使用有限滑动、表面到表面接触公式时，不建议将LM设置为1。如果忽略摩擦，则将LM设置为2（假设无摩擦滑动）。在这种情况下，无需更新其他变量。如果LM始终设置为2，则会创建一个“完全平滑”的界面。您可以根据增量滑移信息和计算的摩擦应力来决定粘滞/滑移条件。这些数量由Abaqus/Standard传递，如下所述。为避免一般摩擦接触问题的收敛问题，将LM设置为2，如果接触点在上一增量结束时打开，则退出此程序；也就是说，如果Abaqus/Standard在调用该例程时设置LM=2，则退出该程序。

如果LM的返回值为0，则有以下变量需要定义

TAU(NFDIR) ：（NFDIR是摩擦力的方向分量个数）这些值作为摩擦应力分量的值传入，τα，并且必须更新为增量结束时的值。这里，以及本说明的其余部分，希腊语下标(α,β)指摩擦剪切方向。这些方向在接触面上的方向在Abaqus/Standard的contact Formulations中定义，α,β分别取1，2。这些方向在接触面上的定向在Contact Formulations in Abaqus/Standard.中定义。

DDTDDG(NFDIR,NFDIR)：元素为 ∂Δτα/∂Δγβ 的数组，α方向的摩擦应力相对于β方向上的相对运动的偏导数，DDTDDG(1，1)，DDTDDG(1，2)，DDTDDG(2，1)，DDTDDG(2，2)

DDTDDP(NFDIR)

DDTDDP(NFDIR) ：∂Δτα/∂Δp，α方向的摩擦应力相对于接触压力的偏导数与接触压力有关。由于这些项对刚度矩阵的贡献是不对称的，因此仅当使用不对称方程解算器时才使用它们（请参见定义分析）。

3、需要被更新的变量

DSLIP(NFDIR)：Δγslα不可恢复滑移运动（滑动）的增量。如果LM在上一次迭代中为0，则此数组在上一个迭代中作为用户定义的值传入；否则，它将为零。仅当LM的返回值为0时，才应更新数组。

此数组用于检测迭代之间的滑动反转。输出选项使用它来指示此点是粘滞还是滑动。增量收敛后，DSLIP（NFDIR）中的值累加到SLIP（NFDIR）中，并存储为塑性应变。

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

SFD：该变量应定义为摩擦耗散增量。如果调用FRIC的接触单元或接触对使用应力而不是力，则单位为单位面积的能量。对于常规应力分析，该变量仅用于输出，对其他解变量没有影响。在温度位移耦合和热电结构耦合分析中，如果使用间隙生热模型，则耗散被转换为热量。如果未定义SFD，则根据作为滑移增量DSLIP和摩擦应力TAU的乘积获得的耗散计算发热。

DDTDDT(NFDIR,2)：∂Δτα / ∂Δθ1， ∂Δτα /∂Δθ2 ,α方向的摩擦应力关于两个表面的温度的偏导数，这仅适用于耦合温度位移和耦合热电结构元件，其中摩擦应力是表面温度的函数。

PNEWDT：建议的新时间增量与使用的时间增量的比率（DTIME，见下文）。此变量允许您为Abaqus/Standard中的自动时间增量算法提供输入（如果选择自动时间增量）。

每次调用FRIC之前，PNEWDT设置为一个较大的值。

如果PNEWDT被重新定义为小于1.0，Abaqus/Standard必须放弃时间增量，并以较小的时间增量再次尝试。为自动时间积分算法提供的建议新时间增量是PNEWDT×DTIME，其中所用的PNEWDTs是对用户子程序的所有调用的最小值，这些调用允许重新定义此迭代的PNEWDT。

如果对于此迭代的所有用户子程序调用，PNEWDT的值都大于1.0，并且增量在此迭代中收敛，则Abaqus/Standard可能会增加时间增量。为自动时间积分算法提供的建议新时间增量是PNEWDT×DTIME，其中使用的PNEWDTs是此迭代中所有用户子程序调用的最小值。

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

STATEV(NSTATV)：包含用户定义的解决方案相关状态变量的数组。指定可用状态变量的数量；有关详细信息，请参阅为依赖于解决方案的状态变量分配空间。这个数组将在增量开始时传入，其中包含这些变量的值。如果任何求解相关状态变量与摩擦行为一起使用，则必须在此子例程中将其更新为增量结束时的值。

4、传递给信息的变量

DGAM(NFDIR)：如果在上一次迭代中将LM设置为0，则该值是当前增量中滑动运动的增量，Δγα

。否则，它将为零。通过与DSLIP（NFDIR）进行比较，可以确定滑动是否在此时变为粘滞，和/或是否在此时发生滑动方向反转。

TAULM(NFDIR)：如果在上一次迭代中将LM设置为1，则该值是增量结束时约束应力的当前值，ταLM。否则，它将为零。通过与临界剪切应力进行比较，可以确定这个点是否由粘结变为滑移。

PRESS：增量步结束时的接触压力，

DPRESS：接触压力增量

DDPDDH：软接触情况下的当前接触刚度（见接触压力超闭合关系）

SLIP(NFDIR)：增量开始时的总不可恢复滑动运动（滑移）γslα。此值是DSLIP（NFDIR）从先前增量的累积值。

TIME(1)：当前增量结束时的当前分析步时间值。

TIME(2)：当前增量结束时的当前总时间值。

DTIME：当前增量步时间增量

NOEL：接触单元（不是定义接触面后，接触面上的单元，而是单独建立的接触单元）的元素标签。如果定义了接触面，则传递为零。

CINAME：用户指定的与摩擦定义相关的面相互作用名称，左对齐。对于接触单元，它是与摩擦定义相关的界面（接触面）定义时给出的的元素集名称；如果为接口定义指定了可选名称，则CINAME将作为该名称传入，左对齐。

SECNAME：从面名称。如果使用接触单元，则作为空白传入。

MAINNAME：主面名称。如果使用接触单元，则作为空白传入。

NPT: 接触单元的积分点编号，如果定义了接触面，则作为0传入。

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

COORDS(MCRD)：包含此点当前坐标的数组。

TEMP(2)：一个包含两个元素的数组，分别为从面节点和与该从面节点相对的主表面的当前温度，

NFDIR：摩擦方向数。

MCRD：接触点的坐标方向数

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

5.子程序接口

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,SECNAME,MAINNAME,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,SECNAME,MAINNAME

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),

3 RCOORD(MCRD),DROT(2,2),TEMP(2),PREDEF(2,*),PROPS(NPROPS)

C 用户自己定义 LM, TAU, DDTDDG, DDTDDP,and, optionally, DSLIP, SED, SFD, DDTDDT, PNEWDT, STATEV等变量（文中有变量具体介绍）。

RETURN

END

6.案例

       若不考虑两个方向上剪应力的相互影响,则有：

​

       接触面的本构关系为克拉夫和邓肯认为肯认为剪应力和相对剪切位移之间符合双曲线关系：

​

       式中共有7个参数，K1、K2 、Rf 、n为非线性指标，通过试验确定;是接触面的界面摩擦角;yw是水的容重;pa是大气压力。

（获取更多资料和学习交流欢迎大家关注公众号冬生亦东生，分享旨在本人复习和交流，没有足够精力在很多地方回复）

​