罚函数法分析一类不可压缩橡胶的大变形

罚函数法分析一类不可压缩橡胶的大变形史守峡杨嘉陵北京航空航天大学固体力学研究所,北京100083摘要基于大变形描述,用罚函数的概念,引入PENN不变量分解的YEOH型本构,建立了适用于分析非线性不可压橡胶材料有限元列式,结合选择降阶积分技术来处理罚因子项,克服了罚因子通常不易选取、计算稳定性差等缺点,改善了计算精度算例表明位移与应力能够很好地与理论解吻合关键词平面应变罚函数不可压缩性橡胶材料具有许多良好的特性密封性能、耐磨擦、柔韧性等,在航天、航空、土建、交通等领域有着广泛的应用研究橡胶材料的力学特性已被人们所重视,通常从描述其变形特性的应变能密度函数出发,一般采用MOONEY型橡胶本构1,但该本构在描述橡胶材料的大变形时与实验有很大的偏差2近来,YEOH3本构与橡胶材料的大变形实验数据吻合橡胶材料在变形过程中,伴随着几何非线性、物理非线性,理论分析具有局限性,有限元成为常用的求解方法,实验4已验证,在承受几千个大气压作用下,其体积变化很小,认为不可压缩的,这又给有限元计算带来了很多困难自1965年,HERMANN5关于不可压缩弹性体的变分原理的开创性工作,有限元分析不可压缩材料的领域取得很大进展,在分析橡胶材料的大变形时,罚有限元是一种重要的方法,有广泛的应用,显然罚因子的选取在罚有限元中占有很重要的地位,一直存在以下缺点首先找出合理的罚因子比较困难,其次罚有限元的计算稳定性较差,不易得到收敛的解本文引入PENN6不变量的分解,采用YEOH型新本构,重点讨论在描述橡胶材料的超大变形时,罚因子的选取及位移、应力解的影响,并与理论解比较,表明计算结果是合理的根据大量实验建议了一类新的本构,该1基本方程11变形描述取一固定的笛卡尔坐标,在外力作用下连续地改变其位形,用XI表示变形前一物质点的坐标LAGRANGE描述,XI,UI分别表示变形后物质点的坐标、位移,即XIXIUI变形梯度F、变形张量G分别为I1,2,315UI2FIJGIJFKIFKJ5XJ1998206215收稿,1998205213收修改稿186应用基础与工程科学学报VOL7取初始构形RX作为参考构形,应力、应变分别取第二PIOLA2KIRCHHOFF应力和柯西2格林应变,分别用SIJ和EIJ表示,则内力做虚功UINT表示为UINTSIJEIJDRX3RX其中EIJ1GIJIJ425W5SIJ5EIJ符号说明W应变能密度函数,IJKRONECKE符号12平衡方程的线性化根据虚功原理,外力在虚位移上所作的虚功UEXT等于内力在虚应变上所作的虚功UINT,即物体的平衡方程为UEXTUINT其中表示对位移的变分,EXT表示外力,INT表示内力,由3式取变分得用增量形式65UI5UKTIJKLDRXUINT5XJ5XLRX5UP5UP5UQ5UQ112DRX7FPIFPJDIJKLFQKFQL5XJ5XI5XL5XL2RX其中TIJKLIKSJL852W9DIJKL5EIJ5EKLTIJKL、DIJKL分别为初应力张量和材料张量又外力虚功UEXT的变分UEXTREXTUDRX10RX最终得增量位移形式的有限元方程KTUNREXTFINT11NN其中N为载荷的分级数NBSDRXFINTT12RXKTKMKG1314KMB0BLTDB0BLDRXRXKGBTDBGDRX15GRX其中T为转置,矩阵B0、BL、BG及的具体形式都可以由7式导出19942014CHINAACADEMICJOURNALELECTRONICPUBLISHINGHOUSEALLRIGHTSRESERVEDHTTP/WWWCNKINETNO2史守峡等罚函数法分析一类不可压缩橡胶的大变形1872本构关系21应变能密度函数YEOH2根据实验提出了适合于橡胶材料大应变的应变能函数,即WA10I13A20I132A30I1331617EI3A10I13A20I13233A30I1W11其中,、A10、A20、A30为材料系数,除了A20小于零外,其余都大于零,根据实验来确定,I1、I2、I3分别为变形张量的第一、第二、第三不变量,如果材料是不可压缩的,则18I31用有限元处理不可压缩橡胶材料,通常在应变能函数中引入罚函数6HI31I312HI31LNI319或22为罚因子,LN为自然对数,PENN6建议把描述橡胶材料变形特征的I1、I2、I3,变形和形状变形进行分解,引入另一组不变量按照体积I121I2I320I1I33I2I33I3进一步,应变能密度函数记为WI,I3WI1HI321122应力与应变的关系变形张量G的三个不变量I21GMMGNN1LMNPQRGLPGMQGNRI1GKKGMNGMNI32226LMN为置换张量,不变量的导数5I15I25I32I3GJ15EIJ22I1IJGIJ23IJI5EJ5EIJI其中GJ1为GIJ的逆张量把式22、23代入5式第二PIOLOA2KIRCHHOFF应力I1IG1SIJ2A102A20I133A30I132I12I31IG1IJ331IJ3IJ324柯西应力IJ与SIJ的关系125IJFKISKLFLJI323材料张量DIJKL据9式,求得增量应力SIJ与增量应变EKL的关系1IG11IG126A30I34IDIJKL332A2031KLIJ1IJKL334I1IG13A102A20I33AI32G111303KLIJ31IJ33KLGJ1II1G1GL14I2GJ1GL13IKI3I31GJ1GL1G1GL1IKJIKIKJ2619942014CHINAACADEMICJOURNALELECTRONICPUBLISHINGHOUSEALLRIGHTSRESERVEDHTTP/WWWCNKINET188应用基础与工程科学学报VOL7式26中DIJKL为四阶张量,共81个分量,对于平面问题,应转换成三阶矩阵,转换十分复杂,积分计算,结合选择降阶积分技术1来处理罚因子项,避免“自锁”现象73数值计算以无限长厚壁圆筒受均匀内压为例1,这是一典型的平面应变问题,取14结构如图1示,几何尺寸内半径R170CM,外半径R218625CM,橡胶材料本构采用16式,有限元数值解与理论解8进行比较YEOH3本构可以反映橡胶材料的大应变,根据实验确定本构参数为A100373MPA、A200037MPA、A300005MPA,在橡胶本构中引入PENN不变量后,为研究罚函数在分析橡胶材料极大变形中的应用,计算模型如图1示,网格采用两种划分105、1020,采用四结点等参元,内压P等分成若干均等载荷步,并给定相对位移收敛精度图1计算模型FIG1THEMODELOFCACULATION为10E3,当罚因子75E4时,我们给出两种不同载荷下的变形形貌如图2示当受内压P为0128MPA时,内壁位移随内压P变化曲线如图3所示,变形后的位移、应力沿径向的分布如图4图6所示,数值结果分别与理论解比较,其中R形参考柱体初始构图2不同内压P作用的大变形形貌FIG2THEMAPOFLARGEDEFORMATIONFORDIFFERENTINNERPRESSURES图3表明,随内压载荷的增加、橡胶圆筒的内壁位移也增加,并与理论分析十分接近,19942014CHINAACADEMICJOURNALELECTRONICPUBLISHINGHOUSEALLRIGHTSRESERVEDHTTP/WWWCNKINETNO2史守峡等罚函数法分析一类不可压缩橡胶的大变形189并与YEOH3实验一致,反映了橡胶材料的极大位移,径向的最大应变达到412,图4至图6反映了沿圆柱径向方向的位移、应力的变化,从图中看出,都与理论解吻合图3内压P2位移U的曲线图4位移U2R的曲线DISPLACEMENTRADIUSCURVEFIG3PRESSUREDISPLACEMENTCURVEFIG4应力RR2R的曲线应力2R的曲线图5图6FIG5RADIUSSTRESS2RADIUSCURVEFIG6CIRCLESTRESS2RADIUSCURVE31罚因子的选取在橡胶本构中引入PENN不变量后,用罚函数在分析橡胶材料极大变形时,提高了有限元解的精度,罚因子在较宽的范围内与理论解接近,因此便于工程上应用,而未采用PENN不变量分解的本构计算表明,选取合适的罚因子困难,计算稳定性较差,位移解易发散当然,罚因子取值也不能太大,不然方程造成“闭锁”现象,计算严重失真,罚因子在19942014CHINAACADEMICJOURNALELECTRONICPUBLISHINGHOUSEALLRIGHTSRESERVEDHTTP/WWWCNKINET190应用基础与工程科学学报VOL725E3至25E6之间,位移与应力与理论解相近,罚因子对内壁位移及应力的影响见表1,应力以单元1为例,未采用PENN不变量分解的本构计算见表2罚因子对位移、应力的影响引入PENN不变量分解表1TABLE1THEINFLUENCESOFPENALTYFACTORTODISPLACEMENTANDSTRESSPENNSINVARIANTDECOMPSIONISUSED罚因子内壁位移CM外壁位移CM单元1形心点的CAUCHY应力MPA网格I3误差误差RR误差误差UR1UR225E425E525E620528820411620401203509209313877013768713759204512212900933009310091018121043142313542856442806343631433610341012100410525E425E525E620684220560620552404102002414010013900113888005002804300961009730098310525032144198643760743823001010909510281003100021020理论解20601013940100951442441表2罚因子对位移、应力的影响未引入PENN不变量分解TABLE2THEINFLUENCESOFPENALTYFACTORTODISPLACEMENTANDSTRESSPENNSINVARIANTDECOMPSIONISNTUSED罚因子内壁位移CM外壁位移CM单元1形心点的CAUCHY应力MPA网格I3UR1误差UR2误差RR误差误差25E32022971851365172070064410525E42029271531366311980064825E5计算发散3233184329124276132153350981098625E320449007613781811400866102025E42044530771379531040086625E5计算发散89389342217542665545835609270996理论解20601013940100951442441从表二得出,在橡胶的YEOH本构中,变形张量的不变量如果未采用PENN不变量的分解,罚因子太小,位移和应力与理论解误差比较大,材料的不可压缩形不能得到较好的满足罚因子较大,计算不稳定,得不到收敛的位移解反之,橡胶本构中引入PENN不变量得分解,罚因子在一较宽的范围内,所得到的位移、应力与理论解十分接近,同时材料的不可压缩形得到较好的满足,更重要的是,我们可以在25E3至25E6之间选取任意一罚因子,都能得到稳定收敛的位移解、应力解结论1罚函数在分析橡胶材料的大变形时是有效的,特别是在橡胶的本构中引入PANN不变量的分解,改善了位移解的稳定性,收敛性,并且与理论解的误差较小2在橡胶的本构中引入PENN不变量的分解,罚因子在25E3至25E6范围内任意选取,都能得到相近的解,合理的罚因子容易找到3罚有限元是分析橡胶平面应变问题或轴对称问题的重要方法之一,本文讨论19942014CHINAACADEMICJOURNALELECTRONICPUBLISHINGHOUSEALLRIGHTSRESERVEDHTTP/WWWCNKINET4NO2史守峡等罚函数法分析一类不可压缩橡胶的大变形191的平面应变问题与相应的平面应力相比要复杂,后者的罚因子可以通过应力消去7,本文的方法可以应用到轴对称、空间问题,因此本文分析具有代表性,建议在橡胶的本构中应引入PENN不变量的分解参考文献12ODENJTFINITEELEMENTOFNONLINEARCONTINUAMCGRAW2HILLNEWYORK,1972SCHARNHORSTT,PIANTHHFINITEELEMENTANALYSISRUBBER2LIKEMATERIALSBYMIXEDMODELINTJNUMERMETHODSENGNG,1978,12665676YEOHOHCHARACTERIZATIONOFELASTICPROPERTIESOFCARBONBLACKFILLEDRUBBERVULCANIZATESRUBBERCHEMICALANDTECHNOLOGY,1990,63792805CHENJSAPRESSUREPROJECTIONMETHODFORNEARLYINCOMPRESSIBLERUBBERHYPERELASTICITYASMEJOURNALOFAP2PLIEDMECHANICS,1996,63862868HERRMANNLRELASTICITYEQUATIONFORINCOMPRESSIBLEANDNEARLYINCOMPRESSIBLEMATERIALSBYAVARIATIONALTHEO2REMAIAAJOURNAL,1965,318961900PENNRWVOLUMECHANGESACCOMPANYINGTHEEXTENSIONOFRUBBERTRANSOFTHESOCOFRHEA,1970,144509517CESCOTTOS,FONDERGAFINITEELEMENTAPPROACHFORLARGESTRAINSOFNEARLYINCOMPRESSIBLERUBBER2LIKEMATERIALSINTJSOLIDSSTRUCT,1988,15589605RIVLINRSLARGEELASTICDEFORMATIONSOFISOTROPICMATERIALSIFUNDAMENTALCONCEPTSPHILTRANSSOC,1948,240259490345678PLANESTRAININCOMPRESSIBLERUBBERFORLARGEDEFORMATIONBYPENALTYFEMSHISHOUXIAYANGJIALINGSOLIDMECHANICSRESEA