！！！歡迎留言讨論、轉發與分享！
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–1FRIC
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
使用者定義的，與摩擦模型相關的屬性值的數量
END