彭從文，董龍彬，吳 群

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

0 引 言

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

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

1 Cosserat連續(xù)體模型

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

應(yīng)力、應(yīng)變分別定義為：

幾何方程為：

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分別是平面內(nèi)平移與轉(zhuǎn)動自由度；mxz,myz偶應(yīng)力；κxz,κyz為微曲率.

對于彈性材料，其本構(gòu)關(guān)系為[7]

S=DeE

(5)

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

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

(6)

(7)

式(7)中,

Cosserat連續(xù)體彈塑性本構(gòu)關(guān)系推導(dǎo)方法同經(jīng)典連續(xù)介質(zhì)力學(xué)，應(yīng)力應(yīng)變增量關(guān)系為

(8)

式(8)中，

硬化項定義為

2 UEL實現(xiàn)方法及計算流程

2.1 子程序編寫注意事項

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

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

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

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

2.2 AMATRX與RHS計算

UEL界面與ABAQUS內(nèi)核主要通過AMATRX與RHS等變量進(jìn)行數(shù)值傳遞.本文采用八節(jié)點等參單元，設(shè)單元節(jié)點位移、插值函數(shù)與位移應(yīng)變轉(zhuǎn)換矩陣分別為d、N與B，單元內(nèi)任一點位移u及應(yīng)變ε為

u=∑Nidi=Nd

(9)

ε=Bd

(10)

B=

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

(11)

對于線性問題，結(jié)合材料本構(gòu)關(guān)系，得到式(12).

(12)

Kd=f′

(13)

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

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

(14)

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

(15)

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

(16)

2.3 計算流程

計算流程如圖1所示.

圖1 UEL流程圖

3 性狀分析

3.1 有限元模型

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

3.2 計算結(jié)果

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

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

圖2 不同網(wǎng)格密度等效塑性應(yīng)變云圖

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

圖4 不同網(wǎng)格密度下應(yīng)力位移曲線

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

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

圖6 不同特征長度的影響

4 結(jié) 語

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

參考文獻(xiàn)：

[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] 楊樂, 吳德倫, 許年春. 偶應(yīng)力理論的層狀巖體洞室數(shù)值模擬[J]. 重慶建筑大學(xué)學(xué)報, 2008, 30(3):73-77.

[5] 尹雪英, 楊春和, 李銀平. 層狀鹽巖體三維Cosserat介質(zhì)擴(kuò)展本構(gòu)模型的程序?qū)崿F(xiàn)[J]. 巖土力學(xué), 2007,28(7):1415-1420,1426.

[6] 朱珍德, 秦天昊, 王士宏, 等. 基于Cosserat理論的柱狀節(jié)理巖體各向異性本構(gòu)模型研究[J]. 巖石力學(xué)與工程學(xué)報,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.