Cosserat彈塑性模型在ABAQUS中的數值實施

彭從文，董龍彬，吳 群

(1.長江大學 城市建設學院，湖北 荊州 434023; 2 荊州市城市規劃設計研究院，湖北 荊州 434000；3.深圳中廣核工程設計有限公司，廣東 深圳 518031)

0 引 言

偶應力理論是微極理論的一個特例，Cosserat兄弟最先提出了完整的偶應力理論[1]，Toupin[2]，Mindlin等[3]對該理論作了進一步的發展和完善.Cosserat理論引入了旋轉自由度和相應的微曲率，引入了與微曲率能量共軛的偶應力、以及具有“特征長度”意義的尺度參數.該理論可以較好地處理網格敏感性和控制方程失去橢圓性的問題，近年來，由于細觀力學、非均質力學的發展，Cosserat理論重新受到關注，逐漸成為研究熱點之一.

數值方法是重要的研究手段之一，為了提高計算精度與效率，基于大型通用數值計算平臺的二次開發方法得到了廣泛的應用.目前,許多數值計算平臺沒有內嵌Cosserat計算模型，關于Cosserat模型的二次開發還不多見[4-6]. ABAQUS是目前最流行、功能最強的商用有限元軟件之一，該軟件可以進行結構靜、動力分析，具有強大的非線性計算能力、豐富的材料庫及良好的擴充功能，自1997年進入我國以來，越來越多的國內企業和研究機構采用ABAQUS作為產品研發和科學研究的工具.本文采用ABAQUS的用戶接口程序，研究壓力相關彈塑性Cosserat連續 體模型的用戶子程序UEL實施方法.

1 Cosserat連續體模型

考慮Cosserat連續體平面問題，每個材料點有三個自由度.

應力、應變分別定義為：

幾何方程為：

ui,j=uj,i-eijkωk

(1)

κij=ωj,i

(2)

靜力平衡方程為：

σij,j+fi=0

(3)

mij,i+eijkσik+qj=0

(4)

式(1)～(4)中，fi、qj分別為體積力與體積力偶；eijk為排列算子；ux,uy,ωz分別是平面內平移與轉動自由度；mxz,myz偶應力；κxz,κyz為微曲率.

對于彈性材料，其本構關系為[7]

S=DeE

(5)

式(5)中，G,v,a,l分別是材料的剪切模量、泊松比、COSSERAT材料參數及特征長度.

對于彈塑性Cosserat材料，采用基于Drucker-Prager屈服準則的彈塑性Cosserate連續體模型，其屈服函數與流動勢函數分別為[8]

(6)

(7)

式(7)中,

Cosserat連續體彈塑性本構關系推導方法同經典連續介質力學，應力應變增量關系為

(8)

式(8)中，

硬化項定義為

2 UEL實現方法及計算流程

2.1 子程序編寫注意事項

(1)ABAQUS提拱了兩類單元自定義方法.一類是線性單元，可以通過結果文件或INP文件直接給出，不需要編寫UEL；另一類就是通用單元，通過UEL子程序定義；

(2)UEL子程序中更新變量與分析問題類別有關.同一個模型中可能遇到不同的分析步，如地應力平衡、靜力分析、攝動步分析等，因此，編寫UEL時要區別處理；

(3)UEL允許自定義荷載.包括集中荷載、均布荷載及彎矩等.其中，對于均布荷載，須定義荷載標志號；

(4)自定義單元在ABAQUS/CAE中不可見.若想在ABAQUS/CAE中顯示自定義單元變形圖，可以將ABAQUS標準單元與自定義單元綁定，同時將標準單元材料參數設為小值.UEL子程序中所有輸出變量均通過SDV寫入結果文件(.fil、.dat)，其分量在ABAQUS/CAE中不可見.

2.2 AMATRX與RHS計算

UEL界面與ABAQUS內核主要通過AMATRX與RHS等變量進行數值傳遞.本文采用八節點等參單元，設單元節點位移、插值函數與位移應變轉換矩陣分別為d、N與B，單元內任一點位移u及應變ε為

u=∑Nidi=Nd

(9)

ε=Bd

(10)

B=

將式(9)、(10)代入Cosserat介質虛功方程(11)進行方程離散.

(11)

對于線性問題，結合材料本構關系，得到式(12).

(12)

Kd=f′

(13)

對于非線性問題，采用Newton-Raphson方法，將式(13)改寫為

ψ(d)=K(d)d-f′

(14)

設ψ(d)為具有一階導數的連續導數，初始近似值為d(0)，第i次迭代的近似值為d(i).將函數ψ(d)在d(i)處展開，保留線性項，忽略高階項得:

(15)

d(i+1)=d(i)+Δd(i)

(16)

2.3 計算流程

計算流程如圖1所示.

圖1 UEL流程圖

3 性狀分析

3.1 有限元模型

模型幾何尺寸10×5 m2，采用三種不同網格密度，單元數分別是15×30、20×40、25×50.邊界條件為底端豎向固定，左側水平向固定.材料彈性模量25 GPa，泊松比0.3，內摩擦角35°，粘聚力1.5 MPa，軟化模量15 MPa.采用相關聯流動法則.頂部采用位移加載方式，加載量為20 mm，加載方向向下.為了觸發局剪切帶，對左下角單元弱化處理.采用高斯完全積分，四階龍格-庫塔顯式應力積分方法.

3.2 計算結果

計算模型分析了不同網格密度及特征長度的影響，計算結果如圖2～6所示.

(1)局部化帶的客觀性.圖2為經典連續介質理論得到的等效塑性應變云圖，圖3與圖4分別是采用Cosserat理論計算得到的等效塑性應變云圖與應力應變曲線.由圖可知，采用Cosserat理論計算時，隨著網格密度增加，剪切帶厚度與等效塑性應變峰值基本不變.當采用經典連續介質理論計算時，計算結果有明顯的網格依賴性，隨著網格密度增加，軟化帶逐漸變窄，等效塑性應變峰值也不斷增大,計算收斂趨于弱化.

圖2 不同網格密度等效塑性應變云圖

圖3 不同網格密度等效塑性應變云圖(特征長度0.15)

圖4 不同網格密度下應力位移曲線

(2)特征長度的影響.Cosserat理論引入特征長度作為正則化機制，特征長度決定Cosserat連續體模型模擬應變局部化問題的能力并影響局部化剪切帶寬度大小.圖5與圖6分別為不同特征長度下采用Cosserat理論計算得到的等效塑性應變云圖和應力位移曲線.由圖2～6可知，隨著特征長度增大，剪切帶厚度增大，等效塑性應變峰值減小，材料軟化模量降低.

圖5 不同特征長度下等效塑性應變云圖(網格20×40)

圖6 不同特征長度的影響

4 結 語

基于ABAQUS接口程序UEL，開發了壓力相關彈塑性Cosserat連續體材料的用戶單元，并采用該單元分析了有限元網格密度及材料特征長度對材料局部化的影響.結果表明，采用Cosserat理論計算時網格密度對材料剪切帶厚度、等效塑性應變影響很小，這也在一定程度上說明本文方法的正確性.要特別說明的是，基于ABAQUS平臺進行二次開發能有效地利用現有程序代碼，減小開發工作量，縮短有限元程序開發周期，極大地提高科研工作效率.

參考文獻：

[1] Cosserat E, Cosserat F. Theorie des Corps Deformables [M].Paris: Herman et Files, 1909.

[2] Toupin R A. Elastic materials with couple stresses [J].Archive Rational Mechanics and analysis, 1962(11): 385-414.

[3] Mindlin R D, Tiersten H F. Effects of couple stresses in linear elasticity [J]. Archive Rational Mechanics and analysis, 1962(11): 415-448.

[4] 楊樂, 吳德倫, 許年春. 偶應力理論的層狀巖體洞室數值模擬[J]. 重慶建筑大學學報, 2008, 30(3):73-77.

[5] 尹雪英, 楊春和, 李銀平. 層狀鹽巖體三維Cosserat介質擴展本構模型的程序實現[J]. 巖土力學, 2007,28(7):1415-1420,1426.

[6] 朱珍德, 秦天昊, 王士宏, 等. 基于Cosserat理論的柱狀節理巖體各向異性本構模型研究[J]. 巖石力學與工程學報,2010,29(增2): 4068-4076.

[7] Sharbati E, Naghdabadi R. computational aspects of the cosserat finite element analysis of localization phenomenon[J]. Computational materials science,2006(38):303-315.

[8] HAYDAR ARSLAN. Localization analysis of granular materials in cosserat elastoplasticity -formulation and finite element Implementation [D]. Colorado: University of Colorado, 2006, 85-88.