基于ABAQUS的复合材料结构材料参数识别算法及其软件设计

1、第30卷第1期2010年2月桂林理工大学学报Journal of Guilin University of Technol ogyVol130No11Feb12010文章编号:1006-544X(201001-0159-05基于ABAQU S的复合材料结构材料参数识别算法及其软件设计李双蓓1,周小军2,刘玉印1,黄立新1,郭相武1(11广西大学土木建筑工程学院,南宁530004;21广州市城市规划勘测设计研究院,广州510060摘要:结合优化技术和有限元分析,提出了正交各向异性复合材料结构材料参数的识别方法。通过建立以测量位移与有限元计算的相应位移之差的平方和作为目标函数,把参数识别问题转化为

2、极小化目标函数的问题。采用Levenberg2Marquardt方法解极小化目标函数的问题,其中灵敏度的计算是基于复合材料的有限元离散结构的求解方程对识别的材料参数求导。基于有限元软件ABAQUS,进行了参数识别算法软件设计,编制了可以计算结点位移对材料参数的灵敏度的用户单元程序。数值算例表明本文提出的方法是有效的。关键词:正交各向异性复合材料;ABAQUS;参数识别;灵敏度分析;软件设计中图分类号:T B330文献标志码:A0引言具有独特优点的复合材料在航空航天、土木工程、机械工程、汽车工程、能源工程和医学等众多的领域得到越来越多的应用,相应地增加了复合材料及其结构设计和评价的需求,而准确的

3、材料参数数据是这些需求的基础。然而,与各向同性材料相比,复合材料性能参数的离散度比较大,在实验室用标准测试方法测定的复合材料性能参数很可能与工厂制造的以及工程中使用的实际构件的性能参数有较大的偏差1-2。在测试各向异性材料参数时还存在边界效应、试样尺寸效应及不易实现各向异性材料系统的均匀应力状态或均匀应变状态等问题3-7。此外,对实际存在的复合材料结构,不可能从中切出一个试件,进行参数测定8。融合测量技术、数值分析方法和优化技术的数值-实验参数识别方法是获取复合材料性能参数的有效途径。材料参数识别属于反问题范畴,而正问题是反问题的基础。正问题必须正确解决,才能提反问题。Kavanagh和Cl

4、ough9首次结合静态测量值与有限元技术,识别正交各向异性复合材料及其结构的材料参数。除了测量静态值之外,许多研究人员又研究了测量动态值的方法进行材料参数的识别,他们采用不同的数值分析方法作为正问题求解方法,Ayorinde等10采用瑞利(Rayleigh方法;De2 obald等11、Moussu等12采用瑞利-里兹(Ray2 leigh2R itz方法;Cunha13和R ikards等7采用有限元法。他们的识别计算中,灵敏度计算采用差分法,计算量大,并且正问题的数值计算没有采用成熟的有限元商业软件,因此材料参数识别方法在适用性方面受到一定的限制。有限元法有建立于严格理论基础的可靠性,并且

5、对复杂几何构形和各种物理问题具有很强的适应性14。经过多年的研究开发,有限元法已经发展了成熟的计算技术和高效可靠的计算商用软件。自然地,在很多数值-实验参数识别方法中,有限元法被选作正问题数值求解的方法。收稿日期:2008-10-20基金项目:广西科学基金项目(桂科自0728005;广西教育厅科研项目(桂教科研200626作者简介:李双蓓(1963,女,博士研究生,副教授,研究方向:复合材料结构与力学。引文格式:李双蓓,周小军,刘玉印,等.基于ABAQUS的复合材料结构材料参数识别算法及其软件设计J.桂林理工大学学报,2010,30(1:159-163.本文选择结构的位移为目标变量,以测量位移

6、与有限元计算的相应位移之差的平方和作为目标函数。采用Levenberg2M arquardt方法极小化目标函数,迭代计算得出材料参数。优化计算中,基于复合材料的有限元离散结构的求解方程对识别的材料参数求导,计算得出灵敏度。基于大型商用有限元软件ABAQUS的计算平台,编制了用户单元CPS4-HZ单元,加载此单元使ABAQUS 可以输出单元结点位移对材料参数的灵敏度矩阵。算例表明本文提出的方法是有效的。1正交各向异性复合材料的有限元求解方程在一定的载荷和边界条件下,基于最小势能原理可以导出复合材料的有限元离散结构的求解方程K(pu(p=R,(1其中:K(p是有限元结构刚度矩阵;u(p是有限元结构

7、结点位移向量;R是有限元结构结点荷载向量;p为材料参数向量。对于平面正交各向异性复合材料来说,结构的x-y坐标与材料主方向1-2坐标不一致时,偏轴应力-应变关系为15xyxy=T-1Q(p(T-1Txyxy,(2其中:T是坐标转换矩阵;Q(p矩阵各非零元素为Q11=E211/(E11-E22212;(3Q12=E11E2212/(E11-E22212;(4Q22=E11E22/(E11-E22212;(5Q66=G12。(6本文中材料参数向量取为p=p1p2p3p4T=E12E2212G12T。(7作为正问题的式(1的求解采用大型商用有限元软件ABAQUS进行。2正交各向异性复合材料结构材料参

8、数识别算法211目标函数目标函数可以定义结构位移的计算值与测量值之误差的平方和,即F(p=Si=1f2i(p=Si=1u i(p-u3i2,pR q,sq。(8其中:ui(p是有限元计算的位移;u3i是测量的位移;s是测量位移的数量;q是识别材料参数的数量。定义向量函数f(p=f1(pf2(pf s(pT=u1(p-u31u2(p-u32us(p-u3sT。(9 212灵敏度分析式(9中测量位移ui与材料参数p无关。f(p对材料参数的灵敏度,即f(p的雅可比矩阵,可表示为J(p=5u i(p5p m s×m,i=1,2,s;m=1,2,3,4。(10式(1的R与材料参数无关,该方程两

9、边对材料参数pm求偏导数,并整理得5u(p5p m=-K(p-15K(p5p mu(p,m=1,2,3,4。(11式中的u(p和K(p均可以在ABAQU S进行正问题分析时求得。根据有限元单元刚度矩阵可得5K(e(p5p m=+1-1+1-1BTT-15Q(p5pm(T-1TB|J|t dd,m=1,2,3,4。(12其中:B是有限元等参单元的应变矩阵;|J|是坐标变换的雅可比行列式;t是等参单元的厚度; 5Q(p5p m各元素表达式可以根据式(3(6求出,以5Q115E11为例5Q11E11=E211-2E11E22v212(E11-E222122。(13其他公式可以同理推出。根据有限元结构

10、总体刚度矩阵14,求得总体刚度矩阵对材料参数的偏导数矩阵,即5K(p5p m=e G T5K(e(p5p m G,m=1,2,3,4。(14其中:G矩阵是与材料参数无关的单元结点自由度061桂林理工大学学报2010年和整体结构结点自由度转换矩阵。由式(12、(14和(11就可求出位移对材料参数的偏导数,即灵敏度。213Le ve nbe rg2M a rqua rd t优化方法极小化目标函数的问题采用Levenberg2M ar2 quardt优化方法解决,其中第k次迭代的材料参数值可以由下列方程计算16J(p(kT J(p(k+(kI(k=-J(p(kT f(p(k;(15p(k+1=p(k

11、+(k。(16其中:I是单位矩阵;(k是迭代过程中非负的Levenberg2M arquard t调整参数。收敛准则定义为J(p(kT f(p(k<1,(17和m axi p(k+1i-p(kip(ki<2。(18其中:1和2是给定的精度要求。3正交各向异性复合材料结构材料参数识别算法的软件设计ABAQUS是国际上最先进的大型通用有限元软件之一,具有强大的力学分析求解功能,同时为用户提供了强大而又灵活的用户子程序接口(USER S UBROUTI N E和应用程序接口(UTI L I2 TY ROUTI N E17-19。本文材料参数识别的正问题求解由大型商用有限元软件ABAQUS

12、完成。针对复合材料的平面问题,有限元计算采用平面四节点等参单元。ABAQUS的单元库中本身就带有平面四节点双线性完全积分单元CPS4,但是ABAQUS自身所带单元不能计算本文所需要的单元刚度矩阵对各材料参数偏导数矩阵,所以在正问题计算时引入用户自定义单元CPS4-HZ参与计算。本文以F OR2 TRAN语言编制了用户自定义单元CPS4-HZ,加载此单元可使ABAQUS输出单元结点位移对材料参数的灵敏度矩阵。在参数识别算法中,总刚度矩阵组装、灵敏度矩阵J(p(k的计算以及Lev2 enberg2M arquardt优化方法求解材料参数过程均由主程序Load CPS4-HZ完成。材料参数识别算法的

13、软件设计流程图参见图1 。图1参数识别算法流程图Fig11Flow chart of parameter identification algorithm4算例均布压力作用下的正交各向异性复合材料长方形板(图2,左端固定,右端自由,板的长度和高度分别为200mm和100mm,单位厚度,均布压力P=20N/mm。材料弹性主方向1与x轴的夹角=0°。硼-环氧复合材料,材料参数20为E11=113000M Pa,E22=52700M Pa,v12=0145,G12=28500MPa。本文材料参数的识别采用以下数据和条件: Levenberg-M arquardt调整参数初值(0=10-14

14、;收敛准则精度要求1=10-6和2=10-3;材料参数识别搜索范围22600M PaE11565000M Pa; 10540M PaE22263500M Pa;012v1201 7;图2均布压力作用下的正交各向异性长方形板Fig12Rectangular orthotrop ic p late under unifor m l oad161第1期李双蓓等:基于ABAQUS的复合材料结构材料参数识别算法及其软件设计 5700M Pa G 12142500M Pa 。测量位移点及其坐标:点A (40,0、点B (100,0及点C (50,90的水平方向及垂直方向的位移。在极小化过程中,材料参数取真

15、值,采用大型商用有限元软件ABAQUS 进行正问题的有限元计算,得到测量点的位移计算值模拟相应的测量位移值。材料参数识别的过程及结果见图3。基于材料参数的先验信息,选择6组离真值较远的材料参数初值,经过若干步迭代之后,材料参数就可以识别出来。从识别的结果来看,该方法具有较高的精度和好的稳定性。5结束语本文提出正交各向异性复合材料结构材料参数的识别方法,数值算例表明该方法具有较高的精度和好的稳定性。材料参数识别的优化计算中,由于灵敏度的计算是根据复合材料有限元离散结构的求解方程求导得到的,因此计算工作量比较小。ABAQUS 是比较成熟的大型商用有限元软件,基于这样的计算平台,该方法有很强的实用性

16、。根据材料参数的先验信息,确定材料参数的搜索范围,可以提高识别的成功率。材料参数初值的选择会影响材料参数识别结果。如果材料参数初值离真值太远,会增加识别的难度甚至会导致识别失败。因此材料参数的先验信息是非常重要的,在识别材料参数时应尽可能利用材料参数的先验信息 。图3不同工况条件下材料参数的识别过程Fig 13Parameter identification p r ocess under different cases参考文献:1W ang W T,Kam T Y .Material characterization of lam inatedcomposite p lates via st

17、atic testing J .Composite Structures,2000,50(4:347-352.2Lecomp te D,S m its A,Sol H,et al .M ixed numerical exper 2i m ental technique f or orthotr op ic parameter identificati on using biaxial tensile tests on crucif or m speci m ens J .I nternationalJournal of Solids and Structures,2007,44(5:1643-

18、1656.3王震鸣.复合材料力学和复合材料结构力学M .北京:机械工业出版社,1991:1-150.4R ikards R,Chate A.I dentificati on of elastic p roperties ofcomposites by method of p lanning of experi ments J .Com 2posite Structures,1998,42(3:257-263.5Pedersen P,Frederiksen P S .Identification of orthotrop ic ma 2terial moduli by a combined ex

19、peri m ental/numerical method261桂林理工大学学报2010年J .Measurement,1992,10(3:113-118.6Frederiksen P S .Experi mental p r ocedure and results f or the i 2dentification of elastic constants of thick orthotr op ic p lates J .Journal of Composite Materials,1997,31(4:360-382.7R ikards R,Chate A,Gailis G .I dent

20、ificati on of elastic p r op 2erties of lam inates based on experi mental design J .I nter 2nati onal Journal of Solids and Structures,2001,38(30-31:5097-5115.8Meuwissen M H H,Oomens C W J,Baaijens F P T,et al .De 2ter m inati on of the elasto 2p lastic p r operties of alum inium using a m ixed nume

21、rical 2experi m ental methodJ .Journal of Materials Processing Technology,1998,75(1-3:204-211.9Kavanagh K T,Clough R W.Finite element app licati ons inthe characterization of elastic s olids J .I nternational Journalof Solids and Structures,1971,7(1:11-23.10Ayorinde E O.Elastic constants of thick or

22、thotr op ic compos 2ite p lates J .Journal of Composite Materials,1995,29(8:1025-1039.11Deobald L R,Gibs on R F .Deter m ination of elastic con 2stants of orthotrop ic p lates by a modal analysis/Rayleigh 2R itz technique J .Journal of Sound and V ibration,1988,124(2:269-283.12Moussu F,N ivoit M.Det

23、er m inati on of elastic constants of or 2thotr op ic p lates by a modal analysis/method of super position J .Journal of Sound and V ibration,1993,165(1:149-163.13Cunha J,Piranda J.App lication of model updating techniquesin dynam ics for the identificati on of elastic constants of com 2posite mater

24、ialsJ .Composites:Part B ,1999,30(1:79-85.14王勖成.有限单元法M .北京:清华大学出版社,2003:1-4,61-72.15沈观林,胡更开.复合材料力学M .北京:清华大学出版社,2006:47-55.16席少霖,赵凤治.最优化计算方法M .上海:上海科学技术出版社,1983:189-200.17H ibbitt,Karlson &Sorensen Inc .ABAQUS/Standard:U ser s Manual M .美国:HKS 公司,2005.18H ibbit,Karlson &Sorensen I nc .ABAQUS

25、/CAE:U ser sManual M .美国:HKS 公司,2005.19庄茁,张帆,岑松.ABAQUS 非线性有限元分析与实例M .北京:科学出版社,2005.20Chern S M ,Tuttle M E .On disp lacement fields in orthotro 2p ic lam inates containing an ellip tical hole J .Journal ofApp lied Mechanics,2000,67(3:527-539.A lgor ithm and Solftware D esi gn of M a ter i a l Param

26、eter I den ti f i ca ti onfor O rthotrop i c Co m posite Structure Ba sed on ABAQUSL I Shuang 2bei 1,ZHOU Xiao 2jun 2,L IU Yu 2yin 1,HUANG L i 2xin 1,Guo Xiang 2wu1(11School of C ivil Engineering and A rchitecture,Guangxi U n iversity,N anning 530004,China;21Guangzhou U rban P lanning &D esign S urvey R esearch Institute,Guangzhou 510060,Ch ina Abstract:A material para meter identificati on method for orthotr op ic composite structure by the app r oach of in 2tegrating op ti m izati on and finite ele ment analysis is p r oposed .The para