罚函数法分析一类不可压缩橡胶的大变形史守峡杨嘉陵( 北京航空航天大学固体力学研究所, 北京 100083)摘要基于大变形描述, 用罚函数的概念, 引入 P en n 不变量分解的 Y eo h 型本构, 建立了适用于分析非线性不可压橡胶材料有限元列式, 结合选择降阶积分技术来处理罚因子项, 克服了罚因子通常不易选取、计算稳定性差等缺点, 改 善了计算精度. 算例说明: 位移与应力能够很好地与理论解吻合.关键词平面应变 罚函数 不可压缩性橡胶材料具有许多良好的特性 (密封性能、耐磨擦、柔韧性等) , 在航天、航空、土建、交通等领域有着广泛的应用. 研究橡胶材料的力学特性已被人们所重视, 通常从描述其变形特性的应变能密度函数出发, 一般采用M oo n ey 型橡胶本构1, 但该本构在描述橡胶材料的大变形时与实验有很大的偏差2 . 近来, Y eo h 3本构与橡胶材料的大变形实验数据吻合.橡胶材料在变形过程中, 伴随着几何非线性、物理非线性, 理论分析具有局限性, 有限 元成为常用的求解方法, 实验4 已验证, 在承受几千个大气压作用下, 其体积变化很小, 认为不可压缩的, 这又给有限元计算带来了很多困难. 自 1965 年, H e rm an n 5 关于不可压 缩弹性体的变分原理的开创性工作, 有限元分析不可压缩材料的领域取得很大进展, 在分析橡胶材料的大变形时, 罚有限元是一种重要的方法, 有广泛的应用, 显然罚因子的选取 在罚有限元中占有很重要的地位, 一直存在以下缺点: 首先找出合理的罚因子比拟困难,其次罚有限元的计算稳定性较差, 不易得到收敛的解. 本文引入 P en n 6 不变量的分解, 采 用 Y eo h 型新本构, 重点讨论在描述橡胶材料的超大变形时, 罚因子的选取及位移、应力解的影响, 并与理论解比拟, 说明计算结果是合理的.根据大量实验建议了一类新的本构, 该1 根本方程1. 1变形描述取一固定的笛卡尔坐标, 在外力作用下连续地改变其位形, 用 X i 表示变形前一物质 点的坐标 (L ag ran ge 描述) , x i , u i 分别表示变形后物质点的坐标、位移, 即x i = X i + u i变形梯度 F、变形张量 G 分别为( i = 1, 2, 3)(1) 5u i (2)F ij =G ij = F k iF k j5X j 1998206215 收稿, 1998205213 收修改稿186应用根底与工程科学学报V o l. 7取初始构形 R x 作为参考构形, 应力、应变分别取第二 P io la2K irch ho ff 应力和柯西2格林应变, 分别用 S i j 和 E ij 表示, 那么内力做虚功 U in t 表示为U in t = S i j E ij dR X(3)R X其中E ij =1 (G ij - ij )(4)25W (5)S ij =5Eij符号说明: W应变能密度函数, ijK ro n eck e 符号.1. 2平衡方程的线性化根据虚功原理, 外力在虚位移上所作的虚功 U ex t 等于内力在虚应变上所作的虚功U in t , 即物体的平衡方程为U ex t = U in t其中 表示对位移的变分, ex t 表示外力, in t 表示内力, 由 (3) 式取变分得用增量形式(6)5u i5 u k= T i j k ldR x +U in t5X j5X lR X5u p5u p5u q5u q 1 12+dR X(7)F p iF p jD ij k lF qkF q l5X j5X i5X l5X l2R X其中T ij k l = ikS j l(8)52W(9)D ij k l =5E ij 5E k lT ij k l、D ij k l 分别为初应力张量和材料张量.又外力虚功 U ex t 的变分U ex t = R ex t u dR X(10)R X最终得增量位移形式的有限元方程 K T u n = R ex t -f in t()11nn其中 n 为载荷的分级数n = B S dR Xf in tT(12)R X K T =K M+K G(13)(14) K M = B 0 + B L TD B 0 + B L d R XR X K G = B TD B G dR X(15)GR X其中 T 为转置, 矩阵 B 0 、B L 、B G 及 的具体形式都可以由 (7) 式导出. 1994-2021 China Academic Journal Electronic Publishing House. All rights reserved.N o. 2史守峡等: 罚函数法分析一类不可压缩橡胶的大变形1872 本构关系2. 1应变能密度函数Y eo h 2根据实验提出了适合于橡胶材料大应变的应变能函数, 即W = A 10 ( I 1 -3) + A 20 ( I 1 - 3) 2 + A 30 ( I 1 - 3) 3(16)(17)e- ( I - 3)+ A 10 ( I 1 - 3) + A 20 ( I 1 - 3) 2 +- 3) 3A 30 ( I 1W =1 -1其中, 、A 10、A 20、A 30 为材料系数, 除了 A 20 小于零外, 其余都大于零, 根据实验来确定,I 1、I 2、I 3 分别为变形张量的第一、第二、第三不变量, 如果材料是不可压缩的, 那么(18)I 3 1用有限元处理不可压缩橡胶材料, 通常在应变能函数中引入罚函数6H ( I 3 ) =1 ( I 3 - 1) 2H ( I 3 ) =1 ln I 3(19)或22为罚因子, ln 为自然对数, P en n 6 建议把描述橡胶材料变形特征的 I 1、I 2、I 3 ,变形和形状变形进行分解, 引入另一组不变量:按照体积I12-1 =I 2 =I 3 =(20)I 1 I 3 3I 2 I 3 3I 3进一步, 应变能密度函数记为W (I , I 3 ) =W(I +1 )( )()HI 32112. 2应力与应变的关系变形张量 G 的三个不变量I 2 =1 (Gm m G nn -1 lm n p q rG lp Gm qG n rI 1 = G k kGm nGm n )I 3 =(22)26lm n 为置换张量, 不变量的导数 5I 1 5I 2 5I 3 2 I 3G -j 15E ij =2=2 ( I 1 i j - G i j )=(23)i ji5E j5E i ji其中 G -j 1 为 G ij 的逆张量. 把式 (22)、(23) 代入 (5) 式第二 P io lo a2K irch ho ff 应力i1 I G - 1S i j = 2 A 10 + 2A 20 (I 1 - 3) + 3A 30 (I 1 - 3) 2 I - 1+ 2( I 3 - 1) I G - 1ij -3 31 i j3 ij3(24)柯西应力 ij 与 S ij 的关系1(25)i j =F k iS k l F l jI 32. 3材料张量 D ij k l据 (9) 式, 求得增量应力 S ij 与增量应变 E k l 的关系:1 I G - 11 I G - 126A 30 (I 34 I -D i j k l =3 3 2 A 20 +-3) (-) (-1 k l )ij1 ijk l334 I -1 I G - 1 )3 A 10 + 2A 20 (I - 3) + 3A(I - 3) 2 G - 1 (11 30 3-k li j -31 i j33+ (k lG -j 1 +iI 1G - 1G - l 1 ) +4 I 2G -j 1G - l 1 +3 ikI 3 ( I 3 - 1) (G -j 1G - l 1 -G - 1G - l 1 ) ikjikikj(26) 1994-2021 China Academic Journal Electronic Publishing House. All rights reserved. 188应用根底与工程科学学报V o l. 7式 (26) 中 D ij k l 为四阶张量, 共 81 个分量, 对于平面问题, 应转换成三阶矩阵, 转换十分复杂, 积分计算, 结合选择降阶积分技术1 来处理罚因子项, 防止“自锁现象7 .3数值计算以无限长厚壁圆筒受均匀内压为例1 , 这是一典型的平面应变问题, 取 14 结构 ( 如图 1 示) , 几何尺寸: 内半径 R 1 =7.0cm , 外半径 R 2= 18. 625cm , 橡胶材料本构采用 ( 16) 式, 有限元数值解与理论解8 进行比拟.Y eo h 3 本构可以反映橡胶材料的大应变, 根据 实 验 确 定 本 构 参 数 为 A 10 = 0. 373M P a、A 20 =-0. 037M P a、A 30 = 0. 005M P a, 在橡胶本构中引入P en n 不变量后, 为研究罚函数在分析橡胶材料极大变形中的应用, 计算模型如图 1 示, 网格采用两 种划分: 105、1020, 采用四结点等参元, 内压 P 等分成假设干均等载荷步, 并给定相对位移收敛精度图 1 计算模型F ig. 1 T h e m o de