基于ABAQUS的復(fù)合材料結(jié)構(gòu)材料參數(shù)識(shí)別算法及其軟件設(shè)計(jì)

1、第30卷第1期2010年2月桂林理工大學(xué)學(xué)報(bào)Journal of Guilin University of Technol ogyVol130No11Feb12010文章編號(hào):1006-544X(201001-0159-05基于ABAQU S的復(fù)合材料結(jié)構(gòu)材料參數(shù)識(shí)別算法及其軟件設(shè)計(jì)李雙蓓1,周小軍2,劉玉印1,黃立新1,郭相武1(11廣西大學(xué)土木建筑工程學(xué)院,南寧530004;21廣州市城市規(guī)劃勘測(cè)設(shè)計(jì)研究院,廣州510060摘要:結(jié)合優(yōu)化技術(shù)和有限元分析,提出了正交各向異性復(fù)合材料結(jié)構(gòu)材料參數(shù)的識(shí)別方法。通過(guò)建立以測(cè)量位移與有限元計(jì)算的相應(yīng)位移之差的平方和作為目標(biāo)函數(shù),把參數(shù)識(shí)別問(wèn)題轉(zhuǎn)化為

2、極小化目標(biāo)函數(shù)的問(wèn)題。采用Levenberg2Marquardt方法解極小化目標(biāo)函數(shù)的問(wèn)題,其中靈敏度的計(jì)算是基于復(fù)合材料的有限元離散結(jié)構(gòu)的求解方程對(duì)識(shí)別的材料參數(shù)求導(dǎo)?；谟邢拊浖嗀BAQUS,進(jìn)行了參數(shù)識(shí)別算法軟件設(shè)計(jì),編制了可以計(jì)算結(jié)點(diǎn)位移對(duì)材料參數(shù)的靈敏度的用戶單元程序。數(shù)值算例表明本文提出的方法是有效的。關(guān)鍵詞:正交各向異性復(fù)合材料;ABAQUS;參數(shù)識(shí)別;靈敏度分析;軟件設(shè)計(jì)中圖分類號(hào):T B330文獻(xiàn)標(biāo)志碼:A0引言具有獨(dú)特優(yōu)點(diǎn)的復(fù)合材料在航空航天、土木工程、機(jī)械工程、汽車工程、能源工程和醫(yī)學(xué)等眾多的領(lǐng)域得到越來(lái)越多的應(yīng)用,相應(yīng)地增加了復(fù)合材料及其結(jié)構(gòu)設(shè)計(jì)和評(píng)價(jià)的需求,而準(zhǔn)確的

3、材料參數(shù)數(shù)據(jù)是這些需求的基礎(chǔ)。然而,與各向同性材料相比,復(fù)合材料性能參數(shù)的離散度比較大,在實(shí)驗(yàn)室用標(biāo)準(zhǔn)測(cè)試方法測(cè)定的復(fù)合材料性能參數(shù)很可能與工廠制造的以及工程中使用的實(shí)際構(gòu)件的性能參數(shù)有較大的偏差1-2。在測(cè)試各向異性材料參數(shù)時(shí)還存在邊界效應(yīng)、試樣尺寸效應(yīng)及不易實(shí)現(xiàn)各向異性材料系統(tǒng)的均勻應(yīng)力狀態(tài)或均勻應(yīng)變狀態(tài)等問(wèn)題3-7。此外,對(duì)實(shí)際存在的復(fù)合材料結(jié)構(gòu),不可能從中切出一個(gè)試件,進(jìn)行參數(shù)測(cè)定8。融合測(cè)量技術(shù)、數(shù)值分析方法和優(yōu)化技術(shù)的數(shù)值-實(shí)驗(yàn)參數(shù)識(shí)別方法是獲取復(fù)合材料性能參數(shù)的有效途徑。材料參數(shù)識(shí)別屬于反問(wèn)題范疇,而正問(wèn)題是反問(wèn)題的基礎(chǔ)。正問(wèn)題必須正確解決,才能提反問(wèn)題。Kavanagh和Cl

4、ough9首次結(jié)合靜態(tài)測(cè)量值與有限元技術(shù),識(shí)別正交各向異性復(fù)合材料及其結(jié)構(gòu)的材料參數(shù)。除了測(cè)量靜態(tài)值之外,許多研究人員又研究了測(cè)量動(dòng)態(tài)值的方法進(jìn)行材料參數(shù)的識(shí)別,他們采用不同的數(shù)值分析方法作為正問(wèn)題求解方法,Ayorinde等10采用瑞利(Rayleigh方法;De2 obald等11、Moussu等12采用瑞利-里茲(Ray2 leigh2R itz方法;Cunha13和R ikards等7采用有限元法。他們的識(shí)別計(jì)算中,靈敏度計(jì)算采用差分法,計(jì)算量大,并且正問(wèn)題的數(shù)值計(jì)算沒(méi)有采用成熟的有限元商業(yè)軟件,因此材料參數(shù)識(shí)別方法在適用性方面受到一定的限制。有限元法有建立于嚴(yán)格理論基礎(chǔ)的可靠性,并且

5、對(duì)復(fù)雜幾何構(gòu)形和各種物理問(wèn)題具有很強(qiáng)的適應(yīng)性14。經(jīng)過(guò)多年的研究開發(fā),有限元法已經(jīng)發(fā)展了成熟的計(jì)算技術(shù)和高效可靠的計(jì)算商用軟件。自然地,在很多數(shù)值-實(shí)驗(yàn)參數(shù)識(shí)別方法中,有限元法被選作正問(wèn)題數(shù)值求解的方法。收稿日期:2008-10-20基金項(xiàng)目:廣西科學(xué)基金項(xiàng)目(桂科自0728005;廣西教育廳科研項(xiàng)目(桂教科研200626作者簡(jiǎn)介:李雙蓓(1963,女,博士研究生,副教授,研究方向:復(fù)合材料結(jié)構(gòu)與力學(xué)。引文格式:李雙蓓,周小軍,劉玉印,等.基于ABAQUS的復(fù)合材料結(jié)構(gòu)材料參數(shù)識(shí)別算法及其軟件設(shè)計(jì)J.桂林理工大學(xué)學(xué)報(bào),2010,30(1:159-163.本文選擇結(jié)構(gòu)的位移為目標(biāo)變量,以測(cè)量位移

6、與有限元計(jì)算的相應(yīng)位移之差的平方和作為目標(biāo)函數(shù)。采用Levenberg2M arquardt方法極小化目標(biāo)函數(shù),迭代計(jì)算得出材料參數(shù)。優(yōu)化計(jì)算中,基于復(fù)合材料的有限元離散結(jié)構(gòu)的求解方程對(duì)識(shí)別的材料參數(shù)求導(dǎo),計(jì)算得出靈敏度?；诖笮蜕逃糜邢拊浖嗀BAQUS的計(jì)算平臺(tái),編制了用戶單元CPS4-HZ單元,加載此單元使ABAQUS 可以輸出單元結(jié)點(diǎn)位移對(duì)材料參數(shù)的靈敏度矩陣。算例表明本文提出的方法是有效的。1正交各向異性復(fù)合材料的有限元求解方程在一定的載荷和邊界條件下,基于最小勢(shì)能原理可以導(dǎo)出復(fù)合材料的有限元離散結(jié)構(gòu)的求解方程K(pu(p=R,(1其中:K(p是有限元結(jié)構(gòu)剛度矩陣;u(p是有限元結(jié)構(gòu)

7、結(jié)點(diǎn)位移向量;R是有限元結(jié)構(gòu)結(jié)點(diǎn)荷載向量;p為材料參數(shù)向量。對(duì)于平面正交各向異性復(fù)合材料來(lái)說(shuō),結(jié)構(gòu)的x-y坐標(biāo)與材料主方向1-2坐標(biāo)不一致時(shí),偏軸應(yīng)力-應(yīng)變關(guān)系為15xyxy=T-1Q(p(T-1Txyxy,(2其中:T是坐標(biāo)轉(zhuǎn)換矩陣;Q(p矩陣各非零元素為Q11=E211/(E11-E22212;(3Q12=E11E2212/(E11-E22212;(4Q22=E11E22/(E11-E22212;(5Q66=G12。(6本文中材料參數(shù)向量取為p=p1p2p3p4T=E12E2212G12T。(7作為正問(wèn)題的式(1的求解采用大型商用有限元軟件ABAQUS進(jìn)行。2正交各向異性復(fù)合材料結(jié)構(gòu)材料參

8、數(shù)識(shí)別算法211目標(biāo)函數(shù)目標(biāo)函數(shù)可以定義結(jié)構(gòu)位移的計(jì)算值與測(cè)量值之誤差的平方和,即F(p=Si=1f2i(p=Si=1u i(p-u3i2,pR q,sq。(8其中:ui(p是有限元計(jì)算的位移;u3i是測(cè)量的位移;s是測(cè)量位移的數(shù)量;q是識(shí)別材料參數(shù)的數(shù)量。定義向量函數(shù)f(p=f1(pf2(pf s(pT=u1(p-u31u2(p-u32us(p-u3sT。(9 212靈敏度分析式(9中測(cè)量位移ui與材料參數(shù)p無(wú)關(guān)。f(p對(duì)材料參數(shù)的靈敏度,即f(p的雅可比矩陣,可表示為J(p=5u i(p5p m s×m,i=1,2,s;m=1,2,3,4。(10式(1的R與材料參數(shù)無(wú)關(guān),該方程兩

9、邊對(duì)材料參數(shù)pm求偏導(dǎo)數(shù),并整理得5u(p5p m=-K(p-15K(p5p mu(p,m=1,2,3,4。(11式中的u(p和K(p均可以在ABAQU S進(jìn)行正問(wèn)題分析時(shí)求得。根據(jù)有限元單元?jiǎng)偠染仃嚳傻?K(e(p5p m=+1-1+1-1BTT-15Q(p5pm(T-1TB|J|t dd,m=1,2,3,4。(12其中:B是有限元等參單元的應(yīng)變矩陣;|J|是坐標(biāo)變換的雅可比行列式;t是等參單元的厚度; 5Q(p5p m各元素表達(dá)式可以根據(jù)式(3(6求出,以5Q115E11為例5Q11E11=E211-2E11E22v212(E11-E222122。(13其他公式可以同理推出。根據(jù)有限元結(jié)構(gòu)

10、總體剛度矩陣14,求得總體剛度矩陣對(duì)材料參數(shù)的偏導(dǎo)數(shù)矩陣,即5K(p5p m=e G T5K(e(p5p m G,m=1,2,3,4。(14其中:G矩陣是與材料參數(shù)無(wú)關(guān)的單元結(jié)點(diǎn)自由度061桂林理工大學(xué)學(xué)報(bào)2010年和整體結(jié)構(gòu)結(jié)點(diǎn)自由度轉(zhuǎn)換矩陣。由式(12、(14和(11就可求出位移對(duì)材料參數(shù)的偏導(dǎo)數(shù),即靈敏度。213Le ve nbe rg2M a rqua rd t優(yōu)化方法極小化目標(biāo)函數(shù)的問(wèn)題采用Levenberg2M ar2 quardt優(yōu)化方法解決,其中第k次迭代的材料參數(shù)值可以由下列方程計(jì)算16J(p(kT J(p(k+(kI(k=-J(p(kT f(p(k;(15p(k+1=p(k

11、+(k。(16其中:I是單位矩陣;(k是迭代過(guò)程中非負(fù)的Levenberg2M arquard t調(diào)整參數(shù)。收斂準(zhǔn)則定義為J(p(kT f(p(k<1,(17和m axi p(k+1i-p(kip(ki<2。(18其中:1和2是給定的精度要求。3正交各向異性復(fù)合材料結(jié)構(gòu)材料參數(shù)識(shí)別算法的軟件設(shè)計(jì)ABAQUS是國(guó)際上最先進(jìn)的大型通用有限元軟件之一,具有強(qiáng)大的力學(xué)分析求解功能,同時(shí)為用戶提供了強(qiáng)大而又靈活的用戶子程序接口(USER S UBROUTI N E和應(yīng)用程序接口(UTI L I2 TY ROUTI N E17-19。本文材料參數(shù)識(shí)別的正問(wèn)題求解由大型商用有限元軟件ABAQUS

12、完成。針對(duì)復(fù)合材料的平面問(wèn)題,有限元計(jì)算采用平面四節(jié)點(diǎn)等參單元。ABAQUS的單元庫(kù)中本身就帶有平面四節(jié)點(diǎn)雙線性完全積分單元CPS4,但是ABAQUS自身所帶單元不能計(jì)算本文所需要的單元?jiǎng)偠染仃噷?duì)各材料參數(shù)偏導(dǎo)數(shù)矩陣,所以在正問(wèn)題計(jì)算時(shí)引入用戶自定義單元CPS4-HZ參與計(jì)算。本文以F OR2 TRAN語(yǔ)言編制了用戶自定義單元CPS4-HZ,加載此單元可使ABAQUS輸出單元結(jié)點(diǎn)位移對(duì)材料參數(shù)的靈敏度矩陣。在參數(shù)識(shí)別算法中,總剛度矩陣組裝、靈敏度矩陣J(p(k的計(jì)算以及Lev2 enberg2M arquardt優(yōu)化方法求解材料參數(shù)過(guò)程均由主程序Load CPS4-HZ完成。材料參數(shù)識(shí)別算法的

13、軟件設(shè)計(jì)流程圖參見圖1 。圖1參數(shù)識(shí)別算法流程圖Fig11Flow chart of parameter identification algorithm4算例均布?jí)毫ψ饔孟碌恼桓飨虍愋詮?fù)合材料長(zhǎng)方形板(圖2,左端固定,右端自由,板的長(zhǎng)度和高度分別為200mm和100mm,單位厚度,均布?jí)毫=20N/mm。材料彈性主方向1與x軸的夾角=0°。硼-環(huán)氧復(fù)合材料,材料參數(shù)20為E11=113000M Pa,E22=52700M Pa,v12=0145,G12=28500MPa。本文材料參數(shù)的識(shí)別采用以下數(shù)據(jù)和條件: Levenberg-M arquardt調(diào)整參數(shù)初值(0=10-14

14、;收斂準(zhǔn)則精度要求1=10-6和2=10-3;材料參數(shù)識(shí)別搜索范圍22600M PaE11565000M Pa; 10540M PaE22263500M Pa;012v1201 7;圖2均布?jí)毫ψ饔孟碌恼桓飨虍愋蚤L(zhǎng)方形板Fig12Rectangular orthotrop ic p late under unifor m l oad161第1期李雙蓓等:基于ABAQUS的復(fù)合材料結(jié)構(gòu)材料參數(shù)識(shí)別算法及其軟件設(shè)計(jì) 5700M Pa G 12142500M Pa 。測(cè)量位移點(diǎn)及其坐標(biāo):點(diǎn)A (40,0、點(diǎn)B (100,0及點(diǎn)C (50,90的水平方向及垂直方向的位移。在極小化過(guò)程中,材料參數(shù)取真

15、值,采用大型商用有限元軟件ABAQUS 進(jìn)行正問(wèn)題的有限元計(jì)算,得到測(cè)量點(diǎn)的位移計(jì)算值模擬相應(yīng)的測(cè)量位移值。材料參數(shù)識(shí)別的過(guò)程及結(jié)果見圖3。基于材料參數(shù)的先驗(yàn)信息,選擇6組離真值較遠(yuǎn)的材料參數(shù)初值,經(jīng)過(guò)若干步迭代之后,材料參數(shù)就可以識(shí)別出來(lái)。從識(shí)別的結(jié)果來(lái)看,該方法具有較高的精度和好的穩(wěn)定性。5結(jié)束語(yǔ)本文提出正交各向異性復(fù)合材料結(jié)構(gòu)材料參數(shù)的識(shí)別方法,數(shù)值算例表明該方法具有較高的精度和好的穩(wěn)定性。材料參數(shù)識(shí)別的優(yōu)化計(jì)算中,由于靈敏度的計(jì)算是根據(jù)復(fù)合材料有限元離散結(jié)構(gòu)的求解方程求導(dǎo)得到的,因此計(jì)算工作量比較小。ABAQUS 是比較成熟的大型商用有限元軟件,基于這樣的計(jì)算平臺(tái),該方法有很強(qiáng)的實(shí)用性

16、。根據(jù)材料參數(shù)的先驗(yàn)信息,確定材料參數(shù)的搜索范圍,可以提高識(shí)別的成功率。材料參數(shù)初值的選擇會(huì)影響材料參數(shù)識(shí)別結(jié)果。如果材料參數(shù)初值離真值太遠(yuǎn),會(huì)增加識(shí)別的難度甚至?xí)?dǎo)致識(shí)別失敗。因此材料參數(shù)的先驗(yàn)信息是非常重要的,在識(shí)別材料參數(shù)時(shí)應(yīng)盡可能利用材料參數(shù)的先驗(yàn)信息 。圖3不同工況條件下材料參數(shù)的識(shí)別過(guò)程Fig 13Parameter identification p r ocess under different cases參考文獻(xiàn):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王震鳴.復(fù)合材料力學(xué)和復(fù)合材料結(jié)構(gòu)力學(xué)M .北京:機(jī)械工業(yè)出版社,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桂林理工大學(xué)學(xué)報(bào)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 .北京:清華大學(xué)出版社,2003:1-4,61-72.15沈觀林,胡更開.復(fù)合材料力學(xué)M .北京:清華大學(xué)出版社,2006:47-55.16席少霖,趙鳳治.最優(yōu)化計(jì)算方法M .上海:上?？茖W(xué)技術(shù)出版社,1983:189-200.17H ibbitt,Karlson &Sorensen Inc .ABAQUS/Standard:U ser s Manual M .美國(guó):HKS 公司,2005.18H ibbit,Karlson &Sorensen I nc .ABAQUS

25、/CAE:U ser sManual M .美國(guó):HKS 公司,2005.19莊茁,張帆,岑松.ABAQUS 非線性有限元分析與實(shí)例M .北京:科學(xué)出版社,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