孔口应力集中问题的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.