孔口應力集中問題的Cosserat連續體有限元分析

唐洪祥 管毓輝

(1大連理工大學土木工程學院,大連 116024)(2中國科學院巖土力學國家重點實驗室,武漢 430071)

當工程結構尺寸或材料特性產生突變時,會發生應力集中現象．對于平面應變狀態下圓孔和橢圓孔的應力集中問題,已經在經典彈性理論下得到了解決[1]．但是在某些實驗中,所得的應力集中因子常常會小于經典彈性理論解[2],這是因為實際工程材料具有自身特有的微結構及內在的特征尺度,在產生應力集中現象的部位具有較陡的應變梯度,經典彈性理論中所假設的均勻性和連續性條件不能得到滿足,因而也不能對其進行準確地模擬．

為了正確模擬應力集中問題中較陡的應變梯度及微結構尺度效應,需要在經典彈性連續體模型中引入某種類型的特征長度參數和高階應變梯度．有效方法之一是采用引入了高階連續體結構的Cosserat微極連續體理論．對于平面二維問題,Cosserat連續體中引入了旋轉軸垂直于2D平面的旋轉自由度及與之相關的微曲率、與微曲率能量共軛的偶應力以及作為材料參數的內部長度尺寸[3-5]．目前,就筆者所知,基于Cosserat連續體理論對孔口應力集中問題的研究主要集中于圓孔應力集中問題的模擬[6-9],而對橢圓孔和菱形孔這2種在實際工程問題中不難見到的孔口應力集中問題的研究還不多見．另外,對于Cosserat連續體在不可壓縮或接近不可壓縮情況下的數值模擬研究中基本都采用了高階單元,盡管低階單元具有計算量小及不易扭曲等優點,但由于其自身局限性,很少得到應用．

本文基于Cosserat連續體理論,建立了3種類型的有限單元,即具有標準位移和旋轉自由度的平面四邊形四節點單元(u4ω4)、平面四邊形八節點單元(u8ω8)以及基于Hu-Washizu混合變分原理推導出的平面四邊形四節點單元(u4ω4p),以模擬平面應變狀態下圓孔、橢圓孔和菱形孔在一般情況和不可壓縮或幾乎不可壓縮情況下的應力集中問題,提出不同條件下合適的數值模擬方法,為實際工程問題提供借鑒和參考．

1 平面Cosserat連續體模型

1.1 彈性Cosserat連續體

在Cosserat連續體平面問題中,作為微極固體的每個材料點具有3個自由度,即

u=[uxuyωz]T

(1)

式中,ux,uy為x-y平面內的平移自由度;ωz為繞z軸旋轉的旋轉自由度．相應地,應變和應力向量定義為

ε=[εxxεyyεzzεxyεyxκzxlcκzylc]T

(2)

(3)

式中,εxx,εyy,εzz,εxy,εyx為常規意義上的應變分量;κzx,κzy分別為繞z軸沿x軸和y軸方向的微曲率;mzx,mzy分別為與κzx,κzy對應的偶應力;lc為內部長度參數．應變-位移關系和平衡方程可表示為

ε=Lu

(4)

LTσ+f=0

(5)

式中,f為外力向量；L為算子矩陣,且

(6)

線彈性應力-應變關系為

σ=Dεe

(7)

式中,εe為彈性應變向量；D為各向同性彈性模量矩陣,且

(8)

式中,λ=2Gv/(1-2v)為Lame常數,其中G,v分別為經典意義上的剪模和泊松比;Gc為Cosserat剪模．

1.2 基于Cosserat連續體理論的常規有限元

1.2.1 平面四邊形四節點單元

標準的平面四邊形四節點單元(u4ω4) 如圖1(a)所示,每個節點具有2個位移自由度和1個旋轉自由度．單元的位移與旋轉由4個節點的雙線性多項式插值逼近,即

(9)

(10)

圖1 2種Cosserat連續體單元

1.2.2 平面四邊形八節點單元

標準的平面四邊形八節點單元(u8ω8) 如圖1(b)所示,每個節點具有2個位移自由度和1個旋轉自由度．單元的位移與旋轉由8個節點的二次多項式插值逼近,即

(11)

(12)

2 常規情況下的孔口應力集中問題

考慮二維平面條件下長方形結構中的圓孔、橢圓孔和菱形孔周圍的應力集中問題(見圖2)．長方形結構中長邊邊長為320.0mm,短邊邊長為60.6mm,分別用u4ω4單元和u8ω8單元來劃分網格．圓孔直徑為19.9mm;橢圓孔的長軸長度為19.9mm,短軸長度為長軸的1/2;菱形孔的長對角線長度為19.9mm,短對角線長度為長對角線的1/2．模型的左邊界固定,右邊界施加600μPa的均布荷載,并等效為節點力施加到對應的節點上．材料參數如下:楊氏模量E=200GPa,泊松比v=0.3,則經典意義上的剪模G=7.69GPa．Cosserat剪切模量Gc和內部長度參數lc作為變量,在模擬結果中給出．

圖2 孔口應力集中問題的計算模型

將應力集中因子定義為孔頂拉應力與右邊界均布力的比值．當內部長度參數lc=r(r為圓孔半徑或橢圓孔、菱形孔的長軸半徑)時,各孔頂端應力集中因子隨Gc/G的變化情況見圖3．由圖可知,隨著Gc/G的增大,孔周圍的應力集中因子逐漸減小,應力集中現象得到弱化;當Gc/G足夠大時,應力集中因子趨于穩定．另外,橢圓孔與菱形孔的應力集中因子明顯大于圓孔,且菱形孔的應力集中因子最大,說明當孔周角度變得尖銳時應力集中因子顯著增大．u4ω4單元和u8ω8單元所得的結果較為一致．

圖3 應力集中因子隨Gc/G的變化規律

當Gc/G=0.5時,各孔頂端應力集中因子隨lc/r的變化情況見圖4．由圖可知,隨著lc/r的增大,孔周圍的應力集中因子逐漸減小,應力集中現象得到弱化;當lc/r足夠大時,應力集中因子趨于穩定．另外,對各孔的應力集中因子進行比較,也可得出與圖3類似的結論,即菱形孔的應力集中因子最大,橢圓孔次之,圓孔最小,且u4ω4單元和u8ω8單元所得到的結果基本相同．需要指出的是,當Gc/G=0或lc/r=0時,Cosserat連續體單元退化為經典連續體單元,所得結果與經典彈性理論的結果一致,即此時圓形孔和橢圓孔的應力集中因子分別為3.440和5.365．

圖4 應力集中因子隨lc/r的變化規律

當lc=r時,經典連續體(Gc/G=0)與Cosserat連續體(Gc/G=0.2)分別對應的圓孔周圍的應力應變云圖見圖5．從圖中可以看出,相對于經典彈性理論解,Cosserat連續體中的應力集中現象明顯弱化．其原因在于,采用Cosserat連續體理論來分析應力集中問題,可以反映大應變梯度和微結構對應力集中的影響．

圖5 圓孔周圍應力和應變分布圖

圖6和圖7分別給出了與圖5類似情況下橢圓孔和菱形孔周圍的應力應變云圖．由圖可知,兩者的結果均與圓孔結果類似;菱形孔周圍的應力集中現象最明顯,橢圓孔次之,圓孔最微弱．

圖6 橢圓孔周圍應力和應變分布圖

3 不可壓縮狀態下的孔口應力集中問題

對于不可壓縮材料或接近不可壓縮材料,當采用位移有限元進行分析時,某些低階單元節點無法自由變位,單元剛度會不合理地增大[10],應力反應也相應放大,導致孔口應力集中因子也會被放大．

圖7 菱形孔周圍應力和應變分布

一般而言,采用減縮積分的高階單元進行模擬是比較合適的,但高階單元的計算量較大,且較大變形所產生的扭曲會造成數值收斂上的困難．鑒于此,本文參考文獻[11]的工作,引入獨立的膨脹場與壓力場,由Hu-Washizu混合變分原理推導得到平面四邊形四節點單元(u4ω4p),以此來分析不可壓縮問題的孔口應力集中問題,并與u4ω4單元和u8ω8單元進行比較．

3.1 u4ω4p單元

(13)

對平面應變問題,有

(14)

引入獨立的膨脹場變量θ,即

(15)

式中,Ve為單元體積．

修改后的應變向量為

(16)

式中,I為單位向量,其定義為

(17)

則有

(18)

引入與獨立膨脹場對應的壓力場p,對于Cosserat連續體,Hu-Washizu混合變分原理可表示為

i=x,y,z;j=x,y,z
k=x,y,z;l=x,y,z

(19)

(20)

離散的虛功表達式為

(21)

(22)

對于u4ω4p單元,k=1,2,3,4,且

(23)

因此,修正應變分量的離散形式可表示為

(24)

3.2 孔口應力集中問題分析

當v=0.450(即接近不可壓縮)時,各種孔口應力集中因子的變化規律見圖8．由圖可見,與v=0.3的情況相比,當u4ω4單元和u8ω8單元蛻化為經典連續體 (即Gc/G=0.0或lc/r=0.0)時,利用其得到的應力集中因子均增大,且u4ω4單元的增大幅度更為明顯;而基于u4ω4p單元得到的應力集中因子基本不變．隨著Gc/G和lc/r的增加,應力集中因子明顯降低,并逐漸趨于穩定．整體上看,u8ω8單元與u4ω4p單元在模擬圓孔和橢圓孔時得到的應力集中因子基本一致,在模擬棱形孔時差別較大;而u4ω4單元在模擬3種孔時得到的應力集中因子均較大．

圖8 v=0.450時孔口應力集中因子的變化規律

當v=0.499(即幾乎不可壓縮)時,各種孔口應力集中因子的變化規律見圖9．由圖可見,與v=0.300及v=0.450的情況相比,u4ω4單元與u8ω8單元在蛻化為經典連續體 (即Gc/G=0或lc/r=0) 時,利用其得到的應力集中因子明顯增大,且u4ω4單元增加顯著,u8ω8單元略有增加;而基于u4ω4p單元得到的應力集中因子基本不變．與v=0.450的情況類似,隨著Gc/G和lc/r的增加,應力集中因子均有所降低,并各自趨于一個穩定值;整體上看,u8ω8單元與u4ω4p單元在模擬圓孔和橢圓孔時得到的應力集中因子基本一致,在模擬棱形孔時則差別較大;而u4ω4單元在模擬3種孔時得到的應力集中因子均較大．

圖9 v=0.499時孔口應力集中因子的變化規律

4 結語

本文基于Cosserat連續體理論,數值模擬了圓孔、橢圓孔以及菱形孔的應力集中問題．對于一般情況下的應力集中問題,采用2種Cosserat單元(即u4ω4單元和u8ω8單元)進行了數值模擬．針對接近不可壓縮或幾乎不可壓縮材料的孔口應力集中問題,基于Hu-Washizu混合變分原理發展了u4ω4p單元來進行數值模擬．結果表明,當孔周角度變得尖銳時,應力集中因子與應變梯度均顯著增大．與經典連續體理論可能高估應力集中因子相比,Cosserat連續體單元可以反映大應變梯度和微結構的影響,從而弱化了應力集中因子．從3種單元數值模擬效果上看,u4ω4p單元與u8ω8單元的性能較好,u4ω4單元較差．考慮到u4ω4p單元具有單元計算量小及不易受扭曲的影響等優點,其應用范圍更加廣闊．

)

[1] 西田正孝. 應力集中 [M]. 李安定，等譯.北京: 機械工業出版社,1986.

[2] 李又村,程赫明,周友坤. 仿真在薄板小孔應力集中現象實驗中的應用[J]. 昆明冶金高等專科學校學報,2008,24(1): 63-68.

Li Youcun,Cheng Heming,Zhou Youkun. Application of simulation in the experiment of the stress concentration of a plate hole [J].JournalofKunmingMetallurgyCollege,2008,24(1): 63-68. (in Chinese)

[3] de Borst R. Simulation of srtain localization:a reappraisal of the Cosserat continuum [J].EngineeringComputations,1991,8(4): 317-332.

[4] Li Xikui,Tang Hongxiang. A consistent return mapping algorithm for pressure-dependent elastoplastic Cosserat continua and modeling of strain localization [J].ComputersandStructures,2005,83(1): 1-10.

[5] 唐洪祥,李錫夔. 地基漸進破壞及極限承載力的Cosserat連續體有限元分析[J]. 巖土力學,2007,28(11): 2259-2264.

Tang Hongxiang,Li Xikui. Finite element analysis of Cosserat continuum for progressive failure and limit bearing capacity of soil foundation [J].RockandSoilMechanics,2007,28(11): 2259-2264. (in Chinese)

[6] Mindlin R D. Influence of couple-stress on stress concentrations [J].ExperimentalMechanics,1963,3(1):1-7.

[7] 劉俊,黃銘,葛修潤. 考慮偶應力影響的應力集中問題求解[J]. 上海交通大學學報,2001,35(10): 1481-1485.

Liu Jun,Huang Ming,Ge Xiurun. Solution of stress concentration problem considering influence of couple-stress [J].JournalofShanghaiJiaotongUniversity,2001,35(10): 1481-1485. (in Chinese)

[8] 趙勇,張若京. 基于Cosserat理論的小孔應力集中問題的有限元分析[J]. 力學季刊,2009,30(3): 410-414.

Zhao Yong,Zhang Ruojing. Finite element analysis of stress concentration problem based on Cosserat theory [J].ChineseQuarterlyofMechanics,2009,30(3): 410-414. (in Chinese)

[9] 張洪武,王輝. 平面Cosserat模型有限元分析的4和8節點單元與分片試驗研究[J]. 計算力學學報,2005,22(5): 512-517.

Zhang Hongwu,Wang Hui. 4- and 8-node isoparametric elements for finite element analysis of plane cosserat bodies [J].ChineseJournalofComputationalMechanics,2005,22(5): 512-517. (in Chinese)

[10] Ted Belytschko,Wing Kamliu,Brian Moran.Nonlinearfiniteelementsforcontinuaandstructures[M]. London: John Wiley and Sons,2000.

[11] Steinmann P. Theory and numerics of ductile micropolar elastoplastic damage [J].InternationalJournalForNumericalMethodsinEngineering,1995,38(4): 583-606.