薄板的后屈曲損傷分析與疲勞壽命預(yù)估

肖浩，胡偉平*，張淼，孟慶春

(1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191;2.中國空間技術(shù)研究院，北京100094)

在航空航天、管道運(yùn)輸?shù)裙こ填I(lǐng)域，結(jié)構(gòu)的穩(wěn)定性問題一直受到重點(diǎn)關(guān)注.結(jié)構(gòu)的屈曲問題早在18世紀(jì)由Euler和Lagrange提出，在這之后，Duncan［1］對不均勻柱進(jìn)行了研究，為屈曲臨界載荷推導(dǎo)了封閉形式的表達(dá)式.此后，對于封閉解的研究較少.1999 年，Elishakoff和 Rollot［2］才發(fā)現(xiàn)了新的封閉解.但工程結(jié)構(gòu)的實(shí)際屈曲載荷明顯要低于由理論解給出的屈曲臨界載荷，當(dāng)時的理論無法解釋這一現(xiàn)象.

1941 年，von Kármán 和 Tsien［3］根據(jù)非線性大撓度理論，提出了板殼后屈曲分析的一般方法以及非線性屈曲理論.Donnell和 Wan［4］于1950年將非線性大撓度分析推廣到非完善圓柱薄殼.Koiter［5］用攝動法研究了彈性結(jié)構(gòu)的初始后屈曲性態(tài)，導(dǎo)出了臨界壓力與缺陷參數(shù)之間的關(guān)系，并提出了初始缺陷敏感度的概念.1968年，Stein［6］提出了非線性前屈曲一致理論，解釋了經(jīng)典線性理論與實(shí)驗(yàn)之間差異的原因.非線性屈曲分析考慮了結(jié)構(gòu)屈曲前的變形，計(jì)算結(jié)果更符合實(shí)際情況.

以von Kármán大撓度理論為依據(jù)，對于單一的屈曲和后屈曲問題已有大量的研究成果，但在實(shí)際工程應(yīng)用中，當(dāng)薄壁構(gòu)件受到循環(huán)載荷作用時，可能會出現(xiàn)屈曲問題與疲勞問題的耦合［7-8］，即結(jié)構(gòu)可能在反復(fù)進(jìn)入后屈曲狀態(tài)下發(fā)生疲勞破壞，或是在產(chǎn)生一定的疲勞損傷累積后發(fā)生屈曲或后屈曲破壞.

目前，國內(nèi)外很多學(xué)者主要研究構(gòu)件在循環(huán)載荷下屈曲后屈曲過程中模態(tài)的改變［9-10］，以及臨界載荷和極限載荷隨著循環(huán)次數(shù)的增加而產(chǎn)生的變化［11-12］，很少有研究將屈曲分析與疲勞損傷分析進(jìn)行耦合，并進(jìn)一步預(yù)估構(gòu)件的疲勞壽命.

本文基于連續(xù)損傷力學(xué)理論與方法，研究結(jié)構(gòu)在后屈曲情況下的損傷以及在疲勞載荷作用下屈曲與疲勞損傷的耦合特性.首先采用有限元方法，通過線性屈曲分析得到屈曲臨界載荷和屈曲模態(tài)，進(jìn)而采用大撓度理論，將線性屈曲的一階屈曲模態(tài)作為初始擾動，進(jìn)行薄板的非線性屈曲分析，得到屈曲臨界載荷.其次，根據(jù)損傷力學(xué)理論與方法建立薄板材料在單次加載過程中的損傷演化方程以及參數(shù)識別方法，并根據(jù)非線性屈曲分析結(jié)果進(jìn)行后屈曲損傷分析.最后，考慮疲勞載荷的作用，基于損傷力學(xué)理論，采用有限元數(shù)值方法進(jìn)行求解，考慮每次加載引起的疲勞損傷與后屈曲應(yīng)力應(yīng)變場的耦合作用，通過反復(fù)迭代計(jì)算，給出結(jié)構(gòu)疲勞壽命.本研究為結(jié)構(gòu)在后屈曲情況下的疲勞特性分析提供了一種新的分析方法和實(shí)現(xiàn)途徑，可以應(yīng)用于飛行器中某些特殊受力狀態(tài)下薄壁構(gòu)件的疲勞壽命分析.

1 薄板的線性屈曲分析

板的線性屈曲分析是一種最常規(guī)的分析方法，可以給出結(jié)構(gòu)的各階屈曲臨界載荷以及相應(yīng)的屈曲模態(tài).

1.1 薄板的有限元模型

考慮一四邊簡支且在某兩對邊受均勻壓力的方形薄板，板邊長a=40mm，板厚h=0.6mm，材料為LC9CgS1［13］，其幾何形狀及載荷如圖1所示.

圖1 結(jié)構(gòu)幾何模型Fig.1 Geometry model of the structure

載荷大小為Px=40 N/mm，四邊簡支板的邊界條件為

式中w為結(jié)構(gòu)變形撓度.

在ANSYS中建立薄板的有限元模型.選擇Shell181單元，該單元適用于薄到中等厚度的殼結(jié)構(gòu)，單元有4個節(jié)點(diǎn)，每個節(jié)點(diǎn)有6個自由度，適用于結(jié)構(gòu)的幾何非線性分析.

1.2 特征值屈曲分析

建立有限元模型后，在ANSYS軟件中進(jìn)行特征值屈曲分析.計(jì)算得到薄板的一階屈曲臨界載荷為34.6 N/mm.

通過特征值屈曲分析能夠給出結(jié)構(gòu)屈曲模態(tài)，能為結(jié)構(gòu)設(shè)計(jì)提供參考，并為后續(xù)的非線性屈曲分析提供初始擾動.

2 板結(jié)構(gòu)的非線性后屈曲分析

在工程實(shí)際中，由于線性屈曲分析忽略了屈曲前變形的影響，導(dǎo)致過高地估計(jì)了結(jié)構(gòu)的臨界載荷.因此，為了提高分析的精度，更多地采取非線性屈曲分析方法進(jìn)行結(jié)構(gòu)穩(wěn)定性分析.

2.1 非線性屈曲分析方法

非線性屈曲分析的目的是求解結(jié)構(gòu)從穩(wěn)定平衡過渡到不穩(wěn)定平衡的臨界載荷，結(jié)構(gòu)的載荷臨界點(diǎn)可以分為兩種類型:分叉臨界點(diǎn)和極值臨界點(diǎn)［14］.

分叉臨界點(diǎn)的特征是:結(jié)構(gòu)在基本載荷-位移平衡路徑Ⅰ的附近還存在另一分叉平衡路徑Ⅱ.當(dāng)載荷達(dá)到臨界值Pcr時，如果結(jié)構(gòu)或載荷受到一微小擾動，載荷-位移曲線將沿分叉平衡路徑發(fā)展，其載荷-位移曲線如圖2(a)所示.

圖2 屈曲載荷臨界點(diǎn)Fig.2 Critical point of buckling load

極值臨界點(diǎn)的特征是:當(dāng)載荷達(dá)到臨界值時，如果載荷位移有微小變化，將分別發(fā)生位移的跳躍和載荷的快速下降，其載荷-位移曲線如圖2(b)所示.

利用有限元軟件進(jìn)行結(jié)構(gòu)的非線性屈曲分析，需要在施加沿板面內(nèi)壓縮載荷的同時，對結(jié)構(gòu)施加一橫向小擾動，通過失穩(wěn)點(diǎn)的載荷-位移曲線計(jì)算屈曲臨界載荷.

2.2 非線性屈曲分析結(jié)果

以特征值屈曲分析得到的第一階屈曲模態(tài)作為初始幾何擾動進(jìn)行非線性屈曲分析，得到失穩(wěn)點(diǎn)的離面位移隨外載荷的變化曲線，如圖3所示.從該曲線可以判斷出非線性屈曲的臨界載荷為31.5 N/mm.當(dāng)外載荷達(dá)到非線性屈曲臨界載荷時，結(jié)構(gòu)進(jìn)入后屈曲，屈曲前和屈曲后對應(yīng)結(jié)構(gòu)的兩個不同的平衡路徑，平衡路徑的轉(zhuǎn)換也是因?yàn)榻Y(jié)構(gòu)上微小擾動的緣故.圖4為結(jié)構(gòu)失穩(wěn)部位的von Mises應(yīng)力隨外載荷的變化曲線，可以看出，微小擾動在前屈曲階段對結(jié)構(gòu)的應(yīng)力影響很小，當(dāng)外載荷達(dá)到屈曲臨界載荷后，結(jié)構(gòu)發(fā)生屈曲并進(jìn)入另一個平衡狀態(tài)，此時結(jié)構(gòu)的應(yīng)力隨外載荷增加而顯著增大.

圖3 失穩(wěn)點(diǎn)離面位移隨外載的變化曲線Fig.3 Variation curve of deflection of instability point with external load

圖4 失穩(wěn)點(diǎn)的von Mises應(yīng)力變化曲線Fig.4 Variation curve of von Mises stress of instability point

3 板結(jié)構(gòu)的后屈曲損傷分析

通過板的后屈曲分析可知，當(dāng)板所受載荷超過屈曲臨界載荷時，板進(jìn)入后屈曲狀態(tài)，此時結(jié)構(gòu)失穩(wěn)部位的應(yīng)力迅速增加.當(dāng)板承受最大載荷超過屈曲臨界載荷的循環(huán)載荷時，板在每一次載荷循環(huán)中都將經(jīng)歷一次后屈曲，因此為了預(yù)估板的疲勞壽命，首先需要進(jìn)行板在后屈曲狀態(tài)下的損傷分析.下面擬采用連續(xù)損傷力學(xué)理論與方法，建立結(jié)構(gòu)在后屈曲狀態(tài)下的損傷分析模型.

3.1 損傷演化模型

根據(jù)連續(xù)損傷力學(xué)理論，在疲勞載荷作用下，材料的損傷將引起材料有效承載面積的下降.在各向同性損傷的情況下，引入損傷度D描述這一物理過程［15］，令

在單向受力狀態(tài)下，建立材料單元的力平衡方程:

對于受損傷的材料單元，其真實(shí)應(yīng)變?yōu)?/p>

在三維受力情況下，對于各向同性材料，根據(jù)應(yīng)變等價(jià)原理，得到含損傷線彈性體本構(gòu)關(guān)系為

式中，λ和μ為拉梅常數(shù);σij和εij分別為應(yīng)力張量和應(yīng)變張量分量.

單軸受力情況下，以應(yīng)變?yōu)榛咀兞浚瑩p傷材料應(yīng)變能密度為

根據(jù)熱力學(xué)原理，引入損傷驅(qū)動力［16］:

單軸受力情況下，有

對處于復(fù)雜受力狀態(tài)下的結(jié)構(gòu)而言，引入損傷力學(xué)等效應(yīng)變:

式中εm為平均應(yīng)變.

采用如下形式的損傷演化方程［17］:

式中，Ymax為損傷力學(xué)驅(qū)動力峰值;Yth為損傷力學(xué)驅(qū)動門檻值;N為載荷循環(huán)次數(shù);mk為材質(zhì)參數(shù).式(12)可以寫成如下應(yīng)變形式的表達(dá)式:

式中，β和εth為材質(zhì)參數(shù);εmax為最大損傷力學(xué)等效應(yīng)變.式(13)需要通過材料標(biāo)準(zhǔn)件疲勞試驗(yàn)結(jié)果識別.

3.2 參數(shù)識別

3.1節(jié)已經(jīng)建立了材料在疲勞載荷作用下的損傷演化方程，方程中包含有β，mk，εth等材質(zhì)參數(shù)，這些參數(shù)需要根據(jù)材料標(biāo)準(zhǔn)件的疲勞試驗(yàn)結(jié)果來確定.

對于光滑試樣，根據(jù)單軸受力的本構(gòu)關(guān)系可得

式中，D0，1為應(yīng)力集中系數(shù)KT=1試件的初始損傷;σ0，max為交變載荷引起的最大應(yīng)力;σ1，th為KT=1時材料發(fā)生疲勞損傷的應(yīng)力門檻值.

對損傷演化方程積分，可得

式中，Nf為裂紋萌生壽命;m為材質(zhì)參數(shù).式(15)可以寫為

式(16)即為KT=1時材料的裂紋萌生壽命和最大應(yīng)力之間的關(guān)系.根據(jù)這一應(yīng)力-壽命表達(dá)式，結(jié)合材料的疲勞試驗(yàn)結(jié)果，就可以確定式(16)中的材質(zhì)參數(shù) β，mk，D0，1及 σ1，th.LC9CgS1材料的靜力力學(xué)性能如表1所示.當(dāng)KT=1，應(yīng)力比R=0.1時，LC9CgS1材料標(biāo)準(zhǔn)件的疲勞壽命如表2所示.

表1 LC9CgS1材料的靜力力學(xué)性能參數(shù)Table1 Static mechanical properties of LC9CgS1

表2 KT=1，R=0.1時 LC9CgS1 材料的疲勞試驗(yàn)結(jié)果Table2 Fatigue test results of LC9CgS1 when KT=1，R=0.1

根據(jù)以上疲勞試驗(yàn)結(jié)果并假設(shè)實(shí)驗(yàn)結(jié)果中離中值壽命曲線最遠(yuǎn)點(diǎn)的實(shí)驗(yàn)點(diǎn)對應(yīng)的材料初始損傷度為0，可以得到損傷演化方程中的材質(zhì)參數(shù)如表3所示.

表3 KT=1，R=0.1時LC9CgS1的材質(zhì)參數(shù)Table3 Material parameters of LC9CgS1 when KT=1，R=0.1

計(jì)算結(jié)果與實(shí)驗(yàn)點(diǎn)的對比如圖5所示.

圖5 KT=1，R=0.1時計(jì)算結(jié)果與實(shí)驗(yàn)點(diǎn)對比曲線Fig.5 Comparison of experimental points and calculation results when KT=1，R=0.1

從圖5可以看出KT=1，R=0.1時擬合的材質(zhì)參數(shù)是合理的.

4 板結(jié)構(gòu)疲勞壽命預(yù)估

當(dāng)板受到的最大載荷超過屈曲臨界載荷的疲勞載荷作用時，每一次載荷循環(huán)中，板都要進(jìn)入后屈曲狀態(tài)，造成局部材料的損傷，這一損傷又會降低下一次加載時板的屈曲臨界載荷，這種屈曲與疲勞損傷的耦合特性難以用閉合形式的方程描述.因此，基于損傷力學(xué)理論，采用有限元數(shù)值解法，通過反復(fù)迭代的后屈曲分析和疲勞損傷分析來預(yù)估板的疲勞壽命.

4.1 損傷力學(xué)-有限元數(shù)值解法

利用ANSYS軟件中的APDL語言控制求解過程，將ANSYS軟件計(jì)算得到的應(yīng)力場代入損傷演化方程，計(jì)算結(jié)構(gòu)的損傷場，給出各單元損傷度;根據(jù)單元損傷度計(jì)算各單元的等效彈性模量，在ANSYS軟件中給各單元重新賦值，再次計(jì)算下一個載荷作用下的應(yīng)力場.如此迭代計(jì)算，在ANSYS平臺中，基于損傷力學(xué)理論，采用有限元數(shù)值解法對構(gòu)件進(jìn)行疲勞壽命預(yù)估.

具體做法如下:

1)利用ANSYS軟件計(jì)算結(jié)構(gòu)在某一載荷循環(huán)下的應(yīng)力應(yīng)變場.

2)根據(jù)上一步計(jì)算得到的應(yīng)力應(yīng)變場，通過APDL語言對ANSYS軟件進(jìn)行二次開發(fā)，計(jì)算每個單元的損傷力學(xué)等效應(yīng)力.

3)利用APDL語言，在ANSYS軟件中計(jì)算每個單元在ΔN次載荷循環(huán)之后，損傷度的增量:

式中，βk為材質(zhì)參數(shù);εe，e，max為單元的最大損傷力學(xué)等效應(yīng)變;εk，th為損傷力學(xué)應(yīng)變門檻值.

那么，可以得到當(dāng)前各個單元的損傷度:

4)由上述計(jì)算得到的各個單元的損傷度可知每個單元的等效彈性模量的變化量，在ANSYS軟件中重新賦予各個單元彈性模量，在上次計(jì)算結(jié)果基礎(chǔ)上，進(jìn)行下一個加載循環(huán)的應(yīng)力應(yīng)變分析.

5)重復(fù)步驟2)～步驟4)的過程，直到某一單元損傷度達(dá)到或超過1，此時的ΔN的累積值即為疲勞裂紋萌生壽命.

需要說明的是，載荷循環(huán)次數(shù)步長ΔN需通過不斷減小ΔN的值來進(jìn)行收斂性驗(yàn)證.

損傷力學(xué)有限元法計(jì)算流程如圖6所示.

圖6 損傷力學(xué)有限元法計(jì)算流程圖Fig.6 Calculation flow chart of damage mechanics with finite element method

4.2 計(jì)算結(jié)果

給定板受到的疲勞載荷為:最大載荷為40 N/mm，載荷比 R=0.1.由于最大載荷大于板的屈曲臨界載荷31.5 N/mm，因此每一次加載循環(huán)中板都將進(jìn)入后屈曲階段.

薄板在后屈曲狀態(tài)下，失穩(wěn)部位局部發(fā)生損傷，這一損傷會導(dǎo)致板的總體剛度矩陣發(fā)生變化，從而影響下一次加載時的屈曲臨界載荷以及結(jié)構(gòu)位移.圖7顯示了不同損傷階段失穩(wěn)部位的載荷位移曲線的變化.從圖7可以看出，板的屈曲臨界載荷隨著損傷的增加而減小，失穩(wěn)點(diǎn)的離面位移隨著損傷的增加而增大.

圖7 不同損傷度下結(jié)構(gòu)失穩(wěn)點(diǎn)位移變化曲線對比Fig.7 Comparisons of deflection changing curves of instability points in different damage degrees

通過計(jì)算，板在上述疲勞載荷作用下的疲勞壽命為6×105次，失效的單元位于失穩(wěn)點(diǎn)，失穩(wěn)點(diǎn)周圍的區(qū)域是損傷比較大的區(qū)域，因?yàn)殡S著損傷的增大，結(jié)構(gòu)的屈曲臨界載荷會逐漸變小，失穩(wěn)區(qū)域會擴(kuò)大，這一部分區(qū)域損傷也會因此而變大.

4.3 初始幾何缺陷對結(jié)構(gòu)裂紋萌生壽命的影響

在前面考慮后屈曲的疲勞壽命分析中，為了模擬屈曲這一物理過程，在每一次屈曲分析時都對結(jié)構(gòu)施加了擾動.

采用的施加擾動方法是引入初始幾何缺陷，即將第一階線性屈曲位移模態(tài)作為初始缺陷施加在板上.然而這里存在一個問題，即作為初始缺陷的一階線性屈曲位移模態(tài)的幅值如何選取，不同幅值對板的疲勞壽命的影響程度如何.

因此計(jì)算了不同缺陷幅值比(幅值比=幾何缺陷幅值/板厚，無量綱)情況下板的疲勞壽命，如圖8所示.

圖8 裂紋萌生壽命隨幾何缺陷幅值比的變化曲線Fig.8 Changing curve of crack initiation life with different geometric imperfection amplitude ratios

從圖8中可以看出，初始幾何缺陷的大小對板的疲勞壽命影響很大，因此在預(yù)估實(shí)際結(jié)構(gòu)疲勞壽命時，應(yīng)首先測量結(jié)構(gòu)的幾何外形，根據(jù)測量結(jié)果盡量準(zhǔn)確估計(jì)結(jié)構(gòu)的初始幾何缺陷，以便得到更加準(zhǔn)確的疲勞壽命分析結(jié)果.

5 結(jié)論

基于損傷力學(xué)理論構(gòu)建了薄板后屈曲狀態(tài)下的裂紋萌生壽命預(yù)估方法，為工程結(jié)構(gòu)考慮后屈曲損傷的疲勞壽命分析提供了一種新的研究方法和手段.

1)首先對板進(jìn)行了線性屈曲分析，得到了屈曲臨界載荷以及屈曲位移模態(tài)，為后續(xù)的非線性屈曲分析提供了位移擾動形式.

2)以線性屈曲分析得到的一階屈曲位移模態(tài)作為板的初始幾何缺陷，對板進(jìn)行非線性屈曲分析，得到了板的非線性屈曲載荷，其值要低于特征值屈曲分析得到的屈曲臨界載荷.

3)建立了板后屈曲損傷分析模型，并根據(jù)板材料的標(biāo)準(zhǔn)件疲勞試驗(yàn)結(jié)果確定損傷演化方程中的各材質(zhì)參數(shù).

4)建立了考慮后屈曲過程的板的疲勞壽命分析方法，通過損傷力學(xué)-有限元數(shù)值解法，可以考慮加載循環(huán)過程中板的后屈曲與疲勞損傷的耦合作用.

5)分析初始幾何缺陷的大小對板裂紋萌生壽命的影響，結(jié)果表明，初始幾何缺陷對結(jié)構(gòu)的疲勞壽命有很明顯影響，因此，一方面在制造過程中應(yīng)盡量提高加工精度，減小幾何缺陷;另一方面，在預(yù)估結(jié)構(gòu)疲勞壽命時應(yīng)根據(jù)實(shí)際測量結(jié)果正確估計(jì)幾何缺陷的大小，以便給出更符合實(shí)際情況的壽命預(yù)估結(jié)果.

References)

［1］ Duncan W J.Galerkin’s method in mechanics and differential equations，ADA951862［R］.London:Aeronautical Research Council，1937.

［2］ Elishakoff I，Rollot O.New closed-form solutions for buckling of a variable stiffness column by Mathematica［J］.Journal of Sound and Vibration，1999，224(1):172-182.

［3］ von Kármán T，Tsien H.The buckling of thin cylindrical shells under axial compression［J］.Journal of the Aeronautics Science，1941，8(6):303-312.

［4］ Donnell L H，Wan C C.Effect of imperfections on buckling of thin cylinders and columns under axial compression［J］.Journal of Applied Mechanics-transactions of the ASME，1950，17(1):73-83.

［5］ Koiter W T.The stability of elastic equilibrium，AD0704124［R］.Standford:Department of Aeronautics and Astronautics of Standford University，1970.

［6］ Stein M.Some recent advances in the investigation of shell buckling［J］.AIAA Journal，1968，6(12):2339-2345.

［7］ Kong C，Bang J，Sugiyama Y.Structural investigation of composite wind turbine blade considering various load cases and fatigue life［J］.Energy，2005，30(11):2101-2114.

［8］ Rhead A T，Butler R，Hunt G W.Post-buckled propagation model for compressive fatigue of impact damaged laminates［J］，International Journal of Solids and Structures，2008，45(16):4349-4361.

［9］ 李世榮，劉平.彈性地基上Euler-Bernoulli梁的熱屈曲模態(tài)躍遷特性［J］.應(yīng)用力學(xué)學(xué)報(bào)，2011，28(1):90-94.Li S R，Liu P.Thermal buckling mode transition characteristics of Euler-Bernoulli beam on the elastic foundation［J］.Journal of Applied Mechanics，2011，28(1):90-94(in Chinese).

［10］ 賀爾銘，張釗，胡亞琪，等.單軸壓縮載荷下夾層梁結(jié)構(gòu)屈曲及皺曲模擬研究［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2012，30(5):668-674.He E M，Zhang Z，Hu Y Q，et al.Buckling and wrinkling study of sandwich beam under uniaxial compression load［J］.Journal of Northwestern Polytechnical University，2012，30(5):668-674(in Chinese).

［11］ 左志亮，蔡健，朱昌宏.帶約束拉桿L形鋼管混凝土短柱的偏壓試驗(yàn)研究［J］.東南大學(xué)學(xué)報(bào):自然科學(xué)版，2010，40(2):346-351.Zuo Z L，Cai J，Zhu C H.Bias experimental study on the binding L-shaped concrete bars filled with steel tubes［J］.Journal of Southeast University:Natural Sciences，2010，40(2):346-351(in Chinese).

［12］ 夏盛來，何景武，海爾瀚.大展弦比機(jī)翼屈曲及后屈曲分析［J］.機(jī)械強(qiáng)度，2012，33(6):907-912.Xia S L，He J W，Hai E H.Buckling and post-buckling analysis of the high aspect ration wing［J］.Mechanics Strength，2012，33(6):907-912(in Chinese).

［13］ 吳學(xué)仁.飛機(jī)結(jié)構(gòu)金屬材料力學(xué)性能手冊:靜強(qiáng)度疲勞/耐久性［M］.北京:航空工業(yè)出版社，1996:225-353.Wu X R.Manual of metal material mechanics performance of the plane structure:static fatigue strength/durability［M］.Beijing:Aviation Industry Press，1996:225-353(in Chinese).

［14］ 王勖成.有限單元法［M］.北京:清華大學(xué)出版社，2003:649-651.Wang X C.The finite element method［M］.Beijing:Tsinghua University Press，2003:649-651(in Chinese).

［15］ Lemaitre J，Desmorat R.Engineering damage mechanics:ductile，creep，fatigue and brittle failures［M］.Berlin:Springer，2005:3-6.

［16］ 張行，崔德渝，孟慶春.斷裂與損傷力學(xué)［M］.北京:北京航空航天大學(xué)出版社，2006:329-331.Zhang X，Cui D Y，Meng Q C.Fracture and damage mechanics［M］.Beijign:Beijing University of Aeronautics and Astronautics Press，2006:329-331(in Chinese).

［17］ 張淼.構(gòu)件疲勞壽命預(yù)估的多變量損傷力學(xué)方法［D］.北京:北京航空航天大學(xué)，2009.Zhang M.Multivariate fatigue damage mechanics methods of fatigue life prediction of the structure［D］.Beijing:Beijing University of Aeronautics and Astronautics，2009(in Chinese).