孔口應(yīng)力集中問題的Cosserat連續(xù)體有限元分析

唐洪祥 管毓輝

(1大連理工大學(xué)土木工程學(xué)院,大連 116024)(2中國科學(xué)院巖土力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,武漢 430071)

當(dāng)工程結(jié)構(gòu)尺寸或材料特性產(chǎn)生突變時,會發(fā)生應(yīng)力集中現(xiàn)象．對于平面應(yīng)變狀態(tài)下圓孔和橢圓孔的應(yīng)力集中問題,已經(jīng)在經(jīng)典彈性理論下得到了解決[1]．但是在某些實(shí)驗(yàn)中,所得的應(yīng)力集中因子常常會小于經(jīng)典彈性理論解[2],這是因?yàn)閷?shí)際工程材料具有自身特有的微結(jié)構(gòu)及內(nèi)在的特征尺度,在產(chǎn)生應(yīng)力集中現(xiàn)象的部位具有較陡的應(yīng)變梯度,經(jīng)典彈性理論中所假設(shè)的均勻性和連續(xù)性條件不能得到滿足,因而也不能對其進(jìn)行準(zhǔn)確地模擬．

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

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

1 平面Cosserat連續(xù)體模型

1.1 彈性Cosserat連續(xù)體

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

u=[uxuyωz]T

(1)

式中,ux,uy為x-y平面內(nèi)的平移自由度;ωz為繞z軸旋轉(zhuǎn)的旋轉(zhuǎn)自由度．相應(yīng)地,應(yīng)變和應(yīng)力向量定義為

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

(2)

(3)

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

ε=Lu

(4)

LTσ+f=0

(5)

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

(6)

線彈性應(yīng)力-應(yīng)變關(guān)系為

σ=Dεe

(7)

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

(8)

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

1.2 基于Cosserat連續(xù)體理論的常規(guī)有限元

1.2.1 平面四邊形四節(jié)點(diǎn)單元

標(biāo)準(zhǔn)的平面四邊形四節(jié)點(diǎn)單元(u4ω4) 如圖1(a)所示,每個節(jié)點(diǎn)具有2個位移自由度和1個旋轉(zhuǎn)自由度．單元的位移與旋轉(zhuǎn)由4個節(jié)點(diǎn)的雙線性多項(xiàng)式插值逼近,即

(9)

(10)

圖1 2種Cosserat連續(xù)體單元

1.2.2 平面四邊形八節(jié)點(diǎn)單元

標(biāo)準(zhǔn)的平面四邊形八節(jié)點(diǎn)單元(u8ω8) 如圖1(b)所示,每個節(jié)點(diǎn)具有2個位移自由度和1個旋轉(zhuǎn)自由度．單元的位移與旋轉(zhuǎn)由8個節(jié)點(diǎn)的二次多項(xiàng)式插值逼近,即

(11)

(12)

2 常規(guī)情況下的孔口應(yīng)力集中問題

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

圖2 孔口應(yīng)力集中問題的計算模型

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

圖3 應(yīng)力集中因子隨Gc/G的變化規(guī)律

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

圖4 應(yīng)力集中因子隨lc/r的變化規(guī)律

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

圖5 圓孔周圍應(yīng)力和應(yīng)變分布圖

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

圖6 橢圓孔周圍應(yīng)力和應(yīng)變分布圖

3 不可壓縮狀態(tài)下的孔口應(yīng)力集中問題

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

圖7 菱形孔周圍應(yīng)力和應(yīng)變分布

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

3.1 u4ω4p單元

(13)

對平面應(yīng)變問題,有

(14)

引入獨(dú)立的膨脹場變量θ,即

(15)

式中,Ve為單元體積．

修改后的應(yīng)變向量為

(16)

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

(17)

則有

(18)

引入與獨(dú)立膨脹場對應(yīng)的壓力場p,對于Cosserat連續(xù)體,Hu-Washizu混合變分原理可表示為

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

(19)

(20)

離散的虛功表達(dá)式為

(21)

(22)

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

(23)

因此,修正應(yīng)變分量的離散形式可表示為

(24)

3.2 孔口應(yīng)力集中問題分析

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

圖8 v=0.450時孔口應(yīng)力集中因子的變化規(guī)律

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

圖9 v=0.499時孔口應(yīng)力集中因子的變化規(guī)律

4 結(jié)語

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

)

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

[2] 李又村,程赫明,周友坤. 仿真在薄板小孔應(yīng)力集中現(xiàn)象實(shí)驗(yàn)中的應(yīng)用[J]. 昆明冶金高等?？茖W(xué)校學(xué)報,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] 唐洪祥,李錫夔. 地基漸進(jìn)破壞及極限承載力的Cosserat連續(xù)體有限元分析[J]. 巖土力學(xué),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] 劉俊,黃銘,葛修潤. 考慮偶應(yīng)力影響的應(yīng)力集中問題求解[J]. 上海交通大學(xué)學(xué)報,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理論的小孔應(yīng)力集中問題的有限元分析[J]. 力學(xué)季刊,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節(jié)點(diǎn)單元與分片試驗(yàn)研究[J]. 計算力學(xué)學(xué)報,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.