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

用戶定義的，與摩擦模型相關的屬性值的數量