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

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

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

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

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

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

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

7、結(jié)點位移向量;R是有限元結(jié)構(gòu)結(jié)點荷載向量;p為材料參數(shù)向量。對于平面正交各向異性復(fù)合材料來說,結(jié)構(gòu)的x-y坐標(biāo)與材料主方向1-2坐標(biāo)不一致時,偏軸應(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作為正問題的式(1的求解采用大型商用有限元軟件ABAQUS進行。2正交各向異性復(fù)合材料結(jié)構(gòu)材料參

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

9、邊對材料參數(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ù)有限元單元剛度矩陣可得5K(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各元素表達式可以根據(jù)式(3(6求出,以5Q115E11為例5Q11E11=E211-2E11E22v212(E11-E222122。(13其他公式可以同理推出。根據(jù)有限元結(jié)構(gòu)

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

11、+(k。(16其中:I是單位矩陣;(k是迭代過程中非負(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è)計ABAQUS是國際上最先進的大型通用有限元軟件之一,具有強大的力學(xué)分析求解功能,同時為用戶提供了強大而又靈活的用戶子程序接口(USER S UBROUTI N E和應(yīng)用程序接口(UTI L I2 TY ROUTI N E17-19。本文材料參數(shù)識別的正問題求解由大型商用有限元軟件ABAQUS

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

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

14、;收斂準(zhǔn)則精度要求1=10-6和2=10-3;材料參數(shù)識別搜索范圍22600M PaE11565000M Pa; 10540M PaE22263500M Pa;012v1201 7;圖2均布壓力作用下的正交各向異性長方形板Fig12Rectangular orthotrop ic p late under unifor m l oad161第1期李雙蓓等:基于ABAQUS的復(fù)合材料結(jié)構(gòu)材料參數(shù)識別算法及其軟件設(shè)計 5700M Pa G 12142500M Pa 。測量位移點及其坐標(biāo):點A (40,0、點B (100,0及點C (50,90的水平方向及垂直方向的位移。在極小化過程中,材料參數(shù)取真

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

16、。根據(jù)材料參數(shù)的先驗信息,確定材料參數(shù)的搜索范圍,可以提高識別的成功率。材料參數(shù)初值的選擇會影響材料參數(shù)識別結(jié)果。如果材料參數(shù)初值離真值太遠(yuǎn),會增加識別的難度甚至?xí)?dǎo)致識別失敗。因此材料參數(shù)的先驗信息是非常重要的,在識別材料參數(shù)時應(yīng)盡可能利用材料參數(shù)的先驗信息 。圖3不同工況條件下材料參數(shù)的識別過程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王震鳴.復(fù)合材料力學(xué)和復(fù)合材料結(jié)構(gòu)力學(xué)M .北京:機械工業(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é)報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)化計算方法M .上海:上海科學(xué)技術(shù)出版社,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 .北京:科學(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