工程p-S-N曲線的小子樣預測方法研究

劉瀟然, 孫秦, 梁珂

(西北工業大學 航空學院, 陜西 西安 710072)

p-S-N曲線是結構工程中進行疲勞可靠性分析與耐久性設計的基礎依據[1-2]，現行估計p-S-N曲線的主要方法有基于最小二乘法的線性回歸方法和極大似然法[3]。大量疲勞試驗結果表明，疲勞壽命的概率分布隨不同應力水平作用而不同，因此最小二乘法估計的p-S-N曲線不是線性無偏估計，近年來有學者提出采用異方差加權最小二乘法線性回歸技術來估計p-S-N曲線[4]。

然而上述方法大多需要成組疲勞試驗數據，由于疲勞試驗需要耗費較多的時間和財力，各應力水平下的疲勞試驗件數常少于10件，甚至只取3～5件。有學者提出在小樣本條件下預測p-S-N曲線的方法，但大多集中在不同可靠度下應力-壽命模型中參數的確定性估計方面[5]。

近年來，基于變量不確定性譜分解的近似隨機函數方法備受關注，非嵌入多項式混沌(non-intrusive polynomial chaos,NIPC)展開方法[6]是其之一，因其可采用較少試驗數據便能有效估計響應量的概率特性，被廣泛應用于結構力學、計算流體力學和可靠性分析等學科中[7]。本文基于NIPC展開方法，提出一種p-S-N曲線的小子樣預測方法，主要技術思路是將應力-壽命模型中的參數視為隨機變量，在較少應力水平下完成小子樣疲勞試驗，并運用優化算法反演壽命模型參數數值，由此獲得應力-壽命模型參數的完整小子樣數據集。將疲勞壽命變換為標準正態隨機變量，并運用所得到的小樣本模型參數數據對各參數進行關于標準壽命變量的NIPC展開，由此計算模型各參數對應試驗應力水平的大子樣數據集， 繼而完成各參數的邊緣概率分布檢驗及參數間相關系數計算。在此基礎上，運用最小二乘擬合方法獲得模型參數統計特征值關于應力水平的定量關系，由此得到應力-壽命模型的全概率描述，即本文的p-S-N曲線。為提高前述方法的預測精度，本文額外選定其他應力水平的疲勞試驗壽命對其進行了貝葉斯修正。通過鋁合金2024-T3小子樣疲勞試驗數據，檢驗了本文方法預測中等壽命區概率疲勞壽命特性的有效性。

1 多變量概率函數的NIPC展開

基于小子樣疲勞壽命試驗數據進行應力-壽命模型的概率分布預測是一個統計學的難題。利用數值大樣本擴充方法是其基本方法途徑，Bootstrap方法是早期的擴充樣本方法之一，但對于樣本數小于10的情況，仍效果不佳，且其主要目的是對統計量進行區間估計。本文的技術途徑是利用NIPC的概率函數展開法來擴充應力-壽命模型中各參數的樣本數據，進而完成基于小樣本試驗數據的p-S-N曲線預測。

NIPC是對隨機函數進行其多項式響應面展開的方法，即采用多項式形式來近似表達一未知隨機函數，其優點在于通過較少樣本數據可獲得隨機函數的概率特性。

設ξ={ξ1,…,ξn}為一組獨立同分布隨機變量,基于NIPC展開,任意隨機函數φ均可表示為隨機變量組中各變量ξk(k=1,…,n)的正交多項式組合

(1)

式中,Bm為m階多項式,公式為:

(2)

式中，φad是單變量ad階帶權正交多項式,且

(3)

假設截斷式(1)右端項至p階多項式,即用有限項來近似φ,將(1)式寫成緊湊格式為

(4)

于是,隨機函數φ分為兩部分,具有隨機性質的多項式基函數ψk(ξ)和其確定性系數Ck。(4)式中的Npc為該展開式的項數,可由下式來確定

(5)

式中,p為多項式展開最高階數,n為隨機變量個數。Npc亦為確定多項式確定性系數所需最小樣本量。

(4)式中正交多項式基函數的選擇依賴于多項式中獨立同分布隨機變量的概率分布類型。根據Askey法則[6],對應隨機變量的不同概率分布,存在不同的最優多項式基函數,如當隨機變量服從正態分布時,對應的最優多項式基函數為Hermit多項式。當帶權正交多項式的權函數與隨機變量概率密度函數形式相同時,(4)式右端的有限項展開式隨階數p的增加依指數收斂于隨機函數φ。

1.1 NIPC展開系數的確定

NIPC展開方法最重要的環節是獲得多項式混沌展開的各項系數Ck。一般情況下,利用一定樣本量(不少于由(5)式確定的數量,最佳為2倍于Npc)的隨機變量和隨機函數值就可以計算出多項式混沌展開的系數。

對于n維輸入隨機向量ξ={ξ1,…,ξn},每取一組樣本,通過試驗或數值模擬得到一個對應的輸出φ*(ξ),取Npc組樣本,由(4)式可得Npc個方程,如(6)式所示,由此即可解出各個系數值。

(6)

1.2 輸出響應的統計特性

當求出多項式各項系數之后,即可計算隨機函數φ的統計特性。由期望值定義

(7)

式中,Dφ是φ的積分域。由(4)式及多項式的正交性,將(7)式中積分域Dφ轉換到隨機變量ξ的積分域Dξ可得

(8)

(8)式說明隨機函數φ的均值為多項式混沌展開式的0階項。同理,可以得到隨機函數φ的方差

(9)

以上分析過程的計算流程如圖1所示。

圖1 NIPC展開法分析步驟

2 p-S-N曲線預測方法

在金屬材料疲勞分析中,中等壽命區最常用的確定性模型為如下的Basquin冪函數式

Sm·N=C

(10)

式中,S為給定應力比下的應力幅值,N為疲勞壽命,m和C為模型參數。對(10)式兩邊取對數得

lgN=-mlgS+lgC

(11)

在同一應力幅值下,由于材料內在的不確定性,疲勞壽命N總是隨機的。以下分析中,認為應力作用是確定量,壽命N服從對數正態分布,故(10)式中m和C也是服從某種分布的隨機變量。

對于工程分析模型,模型變量可分為以下3類:

1) 外部作用量,如(10)式中的應力幅值S,一般可認為是精確給定的;

2) 模型狀態量,如(10)式中的壽命N,可由試驗直接測得;

3) 模型內變量,如(10)式中的m和C,通常不可測,但在前述參量已知條件下可通過模型反演的方法求解。

為使(10)式的確定性應力-壽命模型概率化,本文提出以下方法與步驟予以完整實現。

2.1 m和C的概率化方法

相應于每次試驗壽命N的m和C值可通過(12)式優化反演獲得

min(lgN+mlgS-lgC)2

(12)

在中等壽命區,m取值一般在(0,10]之間,lgC取值一般在(0,30]之間[1]。(12)式即為在約束條件下的最優解問題。

從統計抽樣的觀點,對每一個試驗壽命數據Ni,由(12)式得到的mi和lgCi可認為是模型內變量的一次間接采樣。由此得到的樣本數據可分別用于建立各應力幅下m和lgC關于對數正則隨機壽命N的NIPC展開式。

將每一試驗應力水平下的小子樣試驗壽命進行對數正則化處理,即轉換為標準正態分布,再運用(12)式計算得到的m(Si)和lgC(Si) 數據進行NIPC展開,即由小樣本數據集(需滿足(5)式的要求)按(6)式計算可得各展開式的系數。由這些NIPC展開式,即可實施各應力幅下m(Si)和lgC(Si)的大容量樣本快速計算,隨后按大子樣EDF擬合優度檢驗其給定置信度下的邊緣概率分布。取應力幅Si,記由EDF檢驗得到的兩模型內變量的邊緣概率分布為:

m(Si)～EDF(ωi,ηi)

lgC(Si)～EDF(λi,θi)

(13)

式中,ωi,ηi,λi,θi為概率分布參數。

(14)

為獲得模型內變量概率分布參數及其相關系數關于應力幅的關系,運用最小二乘擬合計算可得:

(15)

由此計算步驟,即可得到任意應力幅值下m和lgC的概率分布特性。

2.2 貝葉斯更新方法

上述大樣本數據計算過程中必然存在先驗假設及模型誤差問題。首先,NIPC展開過程中,假設疲勞壽命N服從對數正態分布,且對數正則化過程中只能使用小子樣試驗數據估計其均值和標準差;其次,(15)式的擬合過程會存在擬合誤差。為適當解決前述矛盾,本文提出一種利用貝葉斯更新方法修正上述誤差的方案:

1) 另取一組應力幅下的小樣本疲勞試驗數據,按前述方法計算m和lgC的樣本數據;

2) 將獲得的樣本數據代入貝葉斯公式,對m和lgC的分布參數進行修正。在前述各應力水平下,m和lgC的概率分布類型可通過(13)式的EDF擬合優度檢驗確定,且其在此應力水平下的概率分布參數可由(15)式計算得到,作為貝葉斯更新的先驗分布參數。

假設m和lgC通過關于正態分布的EDF檢驗,可認為其服從二維正態分布N(μ,Σ),貝葉斯修正后m和lgC的均值和協方差聯合后驗分布為(16)式所示的Normal-inverse-Wishart。

p(μ,Σ|D,μ0,k0,Σ0,v0)=NIW(μ,Σ|μn,kn,Σn,vn)

(16)

式中,μ和Σ的后驗分布參數為:

(17)

式中,n為新增樣本的容量,μ0,κ0,α0和β0分別為μ和Σ的先驗分布參數;

3) 將步驟2得到的貝葉斯修正后模型參數視為原擬合(15)式的新增數據實例,重新擬合m和lgC的概率分布參數關于應力幅的關系,由此改進這些概率分布參數的擬合精度。

2.3 模型參數的抽樣方法

為預測任意應力幅值下疲勞壽命的概率特性,需對m和lgC進行大子樣抽樣。由于m和lgC之間的相關性,需將其等效地變換成獨立正態分布才能進行抽樣,目前常用的變換方法為Nataf變換,其是基于Copula函數變換的一種特例,理論上適用于任意分布的變換,但對于正態化程度較好的變量具有較高精度,適用范圍廣,其基本原理如下。

設有n維隨機變量X=(X1,X2,…Xn)T,已知各變量邊緣概率密度函數為fi(xi),累積分布函數為Fi(xi),相關系數矩陣為ρ=(ρij)n×n。定義n維標準正態隨機變量Y=(Y1,Y2,…,Yn)T,其相關系數矩陣為ρ0=(ρ0ij)n×n,相應的聯合概率密度函數為

(18)

根據等概率變換原則可得X和Y中的變量有如下關系

yi=φ-1(Fi(xi)),i=1,2,…n

(19)

根據Nataf變換理論,利用隱函數求導法則,由(19)式可推導出變量X的聯合概率密度函數為

(20)

根據相關系數定義及(19)式和(20)式可得隨機變量X各變量間的相關系數與對應的標準正態變量Y各變量間的相關系數有以下關系

(21)

當Xi和Xj的邊緣分布及相關系數ρij已知時,通過數值求解(21)式可得Yi和Yj的相關系數ρ0ij。

ρ0是一對稱正定矩陣,對其進行Cholesky分解可得下三角矩陣Γ0,左乘Γ0的逆可將相關標準正態變量Y轉化為獨立標準正態變量U:

(22)

至此,采用拉丁超立方抽樣方法對獨立標準正態變量U抽樣,再對其進行上述變換的逆過程,如(23)式所示,即可得到隨機變量X的樣本。

Y=Γ0U

(23)

3 試驗設計及驗證

為驗證上述方法的有效性,本文在中等壽命區完成了4組不同應力水平下帶孔板的疲勞試驗。由(11)式可知,雙對數坐標的壽命與應力幅滿足線性關系,故本文選取其中2組應力水平下的試驗壽命樣本,確定m和lgC概率特性,并擬合其統計特征值與應力水平的關系,然后對第3組應力水平下的模型參數進行貝葉斯更新,并擬合修正各參數統計特征值與應力水平關系的系數,最終預測第4組應力水平下的概率疲勞壽命并與試驗壽命進行對比。

疲勞試驗件構型如圖2所示,圖中所標注尺寸單位為mm,材料為鋁合金2024-T3。

圖2 疲勞試驗件構形

疲勞試驗獲得的壽命列于表1,應力比均為0.06。

表1 疲勞試驗數據

試驗機為MTS810.13±250 kN電液伺服萬能試驗機,加載頻率為20 Hz,試驗壽命為帶孔板斷裂時壽命。通過Bootstrap方法,可獲得A、C和D組應力水平下對數疲勞壽命均值μ和標準差σ的95%置信區間,結果列于表2,為后文初步驗證本文方法預測概率疲勞壽命參數的合理性做準備。

表2 各組疲勞壽命參數的區間估計

由前述小節理論分析可知,分別構建m和lgC關于疲勞壽命N的NIPC概率多項式函數,由以下2個步驟完成:

1) 選定NIPC的基函數及階數。由于疲勞壽命N服從對數正態分布,可直接變換為服從標準正態分布的隨機變量,其對應的最優基函數為Hermit多項式。工程上一般取多項式階數為3即可有效預測隨機響應的概率特性。由(5)式可知此時僅需4個疲勞壽命N的樣本。

2) 確定NIPC概率多項式函數的各項系數。將疲勞壽命N變換到標準正態空間作為多項式函數中隨機變量的基本樣本值,即(6)式左端的輸入矩陣,由(12)式求出與各疲勞壽命N對應的m和lgC的值,即(6)式右端的輸出響應矩陣。求解(6)式即可獲得構建的NIPC展開式各項系數。

由于不同應力水平下m和lgC的概率特性不同,本文選取A、C 2組試驗為例,各自隨機選取4個疲勞試驗壽命,采用NIPC方法獲得各自的m和lgC關于lgN的Hermit多項式函數。在此基礎上對疲勞壽命N進行抽樣,即可獲得m和lgC的各自大樣本,此時僅是對多項式簡單計算,而并非求解(12)式的優化方程,故節省大量計算時間。再通過m和lgC的大樣本數據對其進行相關性檢驗和關于正態分布的EDF檢驗,結果列于表3,表中Dn為上確界型統計量,W2和A2為平方差型統計量,具體定義可參見相關文獻[8]。

表3 A組和C組試驗數據對應的C和m的統計特性

注：表中,“√”表示在95%置信度下,可接受各變量服從正態分布。

由表3可知,m和lgC的均值在2個應力水平下基本相同,標準差隨著應力幅值增加而減小,相關系數隨應力幅值變化,且隨應力幅值增加由負相關變為正相關。這2個應力水平下的m和lgC的概率密度函數如圖3所示。

圖3 A組和C組應力水平下m和lgC的概率密度函數

在獲得m和lgC概率分布的基礎上,對其進行隨機抽樣,將抽樣值反代入(10)式中獲得疲勞壽命N的樣本,繼而獲得其統計特性,結果列于表4。

表4 由NIPC展開計算所得的對數壽命N的參數估計值

結合表2和表4可知,由NIPC展開法預測的對數疲勞壽命的均值和標準差均落在由Bootstrap獲得的對應的95%區間估計中,且如圖4所示,各組的試驗數據也均落在NIPC預測的疲勞壽命的上下95%分位數區間內,初步驗證了上述算法在小樣本情況下預測概率疲勞壽命的可行性。

圖4 A組和C組試驗數據及NIPC預測的疲勞壽命分布

由于m和lgC的均值不隨應力水平變化,僅擬合m和lgC的標準差σ及其相關系數ρ與應力水平的關系,如(24)式所示[14]。

σ=BlgS+A

ρ=BlgS+A

(24)

由表3中數據求解(24)式中各系數,結果列于表5。

表5 各參數關于應力幅值關系的系數

由(24)式預測B組應力水平下m和lgC的標準差及其相關系數,得到m和lgC的標準差分別為0.033 7和0.026 0,相關系數為-0.272 3。由于此結果與實際情況可能存在偏差,因此通過B組中3個試驗數據運用貝葉斯更新對其進行修正,結果如表6所示。

表6 B組貝葉斯修正后lgC和m的統計特性

結合表6的數據,對(24)式的擬合系數進行修正,結果列于表7,由線性擬合相關系數R的計算值可以看出,本文采用線性擬合是恰當的,B組試驗數據均落在由此預測的疲勞壽命的上下95%分位數區間內,如圖5所示。

圖5 B組試驗數據及NIPC預測的疲勞壽命分布

表7 修正后各參數關于應力幅值擬合系數及線性相關系數R

表7中給出的修正擬合系數,可進一步用于預測任意應力水平下的m和lgC的標準差及其相關系數,進而用于預測該應力水平下的概率疲勞壽命。對于D組試驗,將其應力幅值70.97 MPa代入(24)式即可獲得lgC的標準差為0.014 1,m的標準差為0.018 9,相關系數為0.099 1。再對m和lgC進行隨機抽樣,將抽樣值代入(10)式獲得疲勞壽命的樣本值,進而統計可得對數疲勞壽命的均值為4.781 7,標準差為0.026 5,與通過試驗數據統計的概率特性對比如表8所示,均落在由Bootstrap獲得的對應的95%區間估計內。且D組試驗數據均落在預測的疲勞壽命的上下95%分位數區間內,如圖6所示,說明了由此預測D組的概率疲勞壽命是合理的。

綜上可知,在工程應用中,在高低載2個應力水平下,各做4個試驗件,對m和lgC的標準差及其相

表8 D組疲勞壽命概率預測與試驗結果對比

圖6 D組試驗數據及NIPC預測的疲勞壽命分布

關系數關于應力幅值擬合,在中間應力水平下做3個試驗件,對m和lgC的統計特征值進行貝葉斯修正并重新擬合與應力幅值的關系。在此基礎上即可預測任意應力幅值下疲勞壽命的分布。因此,僅需3個應力水平下11個試驗數據,即可完成中等壽命區p-S-N曲線的預測。通過對D組疲勞壽命的預測,也進一步說明了此方法的合理性。

4 結 論

1) 在非嵌入多項式混沌展開法的基礎上，將應力-壽命模型中的參數概率化，應用少量疲勞壽命樣本數據可以有效預測金屬材料中等壽命區的p-S-N曲線及疲勞壽命概率分布。

2) 利用本文方法對鋁合金帶孔板的疲勞壽命進行了預估，試驗疲勞壽命均處于預測疲勞壽命分布的95%分位數區間內，表明了該方法的有效性。