馬 鑫，關(guān)志東，薛 斌，劉 露

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

0 引言

復(fù)合材料的性能與傳統(tǒng)金屬材料相比，具有較大的變異性。因此在飛機(jī)結(jié)構(gòu)設(shè)計中，如何經(jīng)濟(jì)安全地應(yīng)用復(fù)合材料就成為了一個非常重要的問題。這也就不可避免的觸及了復(fù)合材料設(shè)計的許用值問題。復(fù)合材料的許用值一般采用B基準(zhǔn)值，B基準(zhǔn)值的定義為：一個力學(xué)性能的限定值，在95%的置信度下，90%的性能數(shù)值群的值不低于其值。美國軍用手冊MIL-HDBK-17-1F中給出了計算B基準(zhǔn)值的方法，之后FAA根據(jù)在長期使用經(jīng)驗的基礎(chǔ)上基于以前的方法提出了最新的B基準(zhǔn)值計算方法。我國綜合上述國外研究成果，編訂了航空行業(yè)標(biāo)準(zhǔn)HB 7618-2009《聚合物基復(fù)合材料力學(xué)性能數(shù)據(jù)表達(dá)準(zhǔn)則》。

本文將B基準(zhǔn)值的計算過程通過C#語言編程實現(xiàn)，同時考慮到對兩種不同方法的需求，程序可以通過不同選項完成兩種方法的B基準(zhǔn)值計算，并且通過報告的方式輸出計算結(jié)果。

1 B基準(zhǔn)值計算方法

1.1 B基準(zhǔn)值計算流程

航空行業(yè)標(biāo)準(zhǔn)HB 7618-2009中給出的B基準(zhǔn)及計算過程流程圖如圖1所示，美國軍用手冊MILHDBK-17-1F中給出的B基準(zhǔn)值計算方法流程圖如圖2所示。

1.2 B基準(zhǔn)值計算方法的不同

1.2.1 分布擬合順序不同

MIL-HDBK-17-1F中在進(jìn)行統(tǒng)計處理時，首選威布爾分布，根據(jù)最新的研究成果，F(xiàn)AA認(rèn)為采用威布爾分布計算的結(jié)果過于保守，因此推薦優(yōu)先選擇正態(tài)分布進(jìn)行擬合優(yōu)度檢驗。HB 7618-2009中采用了FAA推薦的方法首選正態(tài)分布。

1.2.2 環(huán)境數(shù)據(jù)處理的不同

MIL-HDBK-17-1F對于多試驗環(huán)境的的B基準(zhǔn)值計算，認(rèn)為B基準(zhǔn)值與溫度為線性關(guān)系，但是這種線性關(guān)系并不總是適用的，并且無法處理不同濕度環(huán)境的數(shù)據(jù)集。HB 7618-2009對于多環(huán)境數(shù)據(jù)進(jìn)行歸一化后合并計算B基準(zhǔn)值。

1.2.3 離散系數(shù)處理的不同

MIL-HDBK-17-1F中對數(shù)據(jù)的離散系數(shù)沒有要求，HB 7618-2009采用FAA推薦的方法，F(xiàn)AA認(rèn)為復(fù)合材料性能的離散系數(shù)一般為4%～10%，低于此值則為不正常的，會導(dǎo)致較高的基準(zhǔn)值，需要對數(shù)據(jù)進(jìn)行變換將數(shù)據(jù)的離散系數(shù)擴(kuò)大為4%。

1.2.4 子體相容性檢驗的不同

在子體相容性檢驗時，采用K-AD檢驗的方法，MIL-HDBK-17-1F中認(rèn)為具有5%的錯判風(fēng)險，而HB 7618-2009中檢驗失敗時允許通過將錯判風(fēng)險改為1%來增大臨界值，增大了數(shù)據(jù)來自同一母體的可能性。

圖1 HB 7618-2009 B基準(zhǔn)值計算流程圖

表1 兩種方法計算結(jié)果對比

通過兩種方法計算B基準(zhǔn)值的結(jié)果如表1所示。其中ETW數(shù)據(jù)通過MIL-HDBK-17-1F的方法時，認(rèn)為數(shù)據(jù)來自不同的母體，而采用航標(biāo)對臨界值放大后則可以將數(shù)據(jù)合并為非結(jié)構(gòu)型數(shù)據(jù)。計算過程中，ETW數(shù)據(jù)在具有5%錯判風(fēng)險時為非結(jié)構(gòu)型數(shù)據(jù)，但在1%錯判風(fēng)險時子體相容性檢驗表明可以將其作為非機(jī)構(gòu)性數(shù)據(jù)處理。

由于HB 7618-2009方法將各個環(huán)境下數(shù)據(jù)合并后計算折減系數(shù)，所以從結(jié)果計算結(jié)果中可以看出，當(dāng)環(huán)境數(shù)據(jù)均值較高時會降低相應(yīng)基準(zhǔn)值，同時提高均值較低環(huán)境的基準(zhǔn)值，使不同環(huán)境下B基準(zhǔn)值趨于一致。由于將不同環(huán)境的數(shù)據(jù)合并后增大了樣本容量，因此采用此方法計算的B基準(zhǔn)值總體上提高了4.7%。

1.3 計算方法

程序可以按照圖1和圖2的流程圖完成不同的計算過程，計算過程使用的統(tǒng)計學(xué)方法主要包括子體相容性檢驗、異常數(shù)據(jù)檢查、擬合優(yōu)度檢驗、非參數(shù)法和方差分析（ANOVA）方法。

1.3.1 子體相容性檢驗

程序子體相容性檢驗采用k樣本的Anderson-Darling檢驗用來檢驗從中抽取的兩組或多組數(shù)據(jù)的母體是相同的假設(shè)的非參數(shù)統(tǒng)計方法，該檢驗要求各組為來自某一母體的獨立隨機(jī)樣本。

如果檢驗表明可以認(rèn)為數(shù)據(jù)來自同一母體，則可以將認(rèn)為數(shù)據(jù)屬于非結(jié)構(gòu)型數(shù)據(jù)集。一個非結(jié)構(gòu)型數(shù)據(jù)集是一個簡單隨機(jī)樣本。將結(jié)構(gòu)型數(shù)據(jù)簡化為非結(jié)構(gòu)型數(shù)據(jù)是有利的，當(dāng)不能簡化時，只能應(yīng)用結(jié)構(gòu)型的處理方法。K樣本AD檢驗統(tǒng)計量為：

錯判風(fēng)險5%時臨界值為：

圖2 MIL-HDBK-17-1F中的B基準(zhǔn)值計算流程圖

當(dāng)ADK大于ADC時認(rèn)為則可以斷定各組是從不同母體中抽取。否則，接受各組來自同一母體的假設(shè)，可以將數(shù)據(jù)作為非結(jié)構(gòu)型數(shù)據(jù)進(jìn)行處理。HB 7618-2009中提出當(dāng)子體相容性檢驗不符合同一母體假設(shè)時，可以選擇α=0.01重新進(jìn)行檢驗。

當(dāng)無法通過檢驗時，可以通過勾選1%檢驗復(fù)選框重新按照α=0.01進(jìn)行檢驗，以確實是否能將數(shù)據(jù)合并為非結(jié)構(gòu)型數(shù)據(jù)。

1.3.2 異常數(shù)據(jù)檢查

程序通過最大賦范殘差法（MNR）來對異常數(shù)據(jù)進(jìn)行篩選，這是一種識別非結(jié)構(gòu)型數(shù)據(jù)集中異常數(shù)據(jù)的篩選方法。這種方法每次只能從一個選定的數(shù)據(jù)組或者子集中檢出一個異常數(shù)據(jù)，因此，一旦檢出了一個異常數(shù)據(jù)，就必須處置這個數(shù)據(jù)，并重新進(jìn)行分析以檢測另外的異常數(shù)據(jù)。

MNR統(tǒng)計量為：

檢驗時，將上式的統(tǒng)計量與該樣本大小n的臨界值相比較，臨界值可以通過查表或者下式計算：

如果MNR統(tǒng)計量小于該臨界值，則該樣本中未檢出異常數(shù)據(jù)，如果MNR統(tǒng)計量大于該臨界值，則認(rèn)為相應(yīng)的數(shù)據(jù)值是一個異常數(shù)據(jù)。如果存在異常值，從樣本中剔除該值后再次進(jìn)行異常數(shù)據(jù)檢查，重復(fù)此過程直至檢查不出異常數(shù)據(jù)。

當(dāng)數(shù)據(jù)量較小或者分散系數(shù)較小時，可能會出現(xiàn)異常值的誤報或者漏報，所以，需要通過目視檢查和最大賦范殘差法結(jié)合來進(jìn)行異常數(shù)據(jù)檢查。

1.3.3 擬合優(yōu)度檢驗

擬合優(yōu)度檢驗采用對尾端區(qū)域差異敏感的Anderson-Darling檢驗統(tǒng)計量檢驗各種分布。Anderson-Darling檢驗將所考慮的分布累積分布函數(shù)與數(shù)據(jù)的累積分布函數(shù)進(jìn)行比較。基于Anderson-Darling檢驗統(tǒng)計量計算每個檢驗的觀測顯著性水平(OSL)：

OSL為數(shù)據(jù)來自檢驗分布的假設(shè)正確時統(tǒng)計量與計算值相等的概率，當(dāng)OSL小于或等于0.05時，拒絕數(shù)據(jù)來自該分布的假設(shè)（存在5%的誤判風(fēng)險），認(rèn)為數(shù)據(jù)不是來自被檢驗的分布。

根據(jù)擬合分布求得單側(cè)B基準(zhǔn)值容限系數(shù)kB后計算B基準(zhǔn)值：

按照MIL手冊中的要求分布擬合順序為威布爾分布、正態(tài)分布和對數(shù)正態(tài)分布，HB標(biāo)準(zhǔn)中將擬合順序調(diào)整為正態(tài)分布、對數(shù)正態(tài)分布和威布爾分布。程序為了滿足不同用戶的要求，同時給出三種分布的顯著性水平（OSL）和計算得到的B基準(zhǔn)值，用戶可以通過自己的需求安排不同的擬合順序選擇。

1.3.4 非參數(shù)法

如果檢驗表明數(shù)據(jù)的總體分布與正態(tài)分布、對數(shù)正態(tài)分布和威布爾分布均有顯著區(qū)別時，認(rèn)為總體分布是未知的，此時采用非參數(shù)方法計算基準(zhǔn)值。

對大樣本試驗數(shù)據(jù)（n〉28）時，對數(shù)據(jù)進(jìn)行排序，選取第rB個觀測值作為基準(zhǔn)值，rB可以通過查表獲得，程序采用如下公式進(jìn)行計算：

四舍五入取整。通過近似公式得出結(jié)果偏于保守。

當(dāng)樣本大小不超過28時，如果數(shù)據(jù)來自累積分布函數(shù)的對數(shù)為上凹的母體的隨機(jī)樣本，那么可以通過Hanson-Hoopmans方法計算B基準(zhǔn)值。經(jīng)驗表明復(fù)合材料的強度一般滿足該假設(shè)。基準(zhǔn)值計算公式為：

x(1)和x(r)代表最小和第r個最小數(shù)據(jù)，r和k可以通過查表得到，程序中直接讀取列表的值進(jìn)行B基準(zhǔn)值計算。

1.3.5 方差分析方法

如果子體相容性檢驗表明批次間存在明顯差異，試驗數(shù)據(jù)作為結(jié)構(gòu)型數(shù)據(jù)進(jìn)行處理，采用方差分析（ANOVA）方法來計算B基準(zhǔn)值。使用方差分析方法需要滿足一下假設(shè)：

(1)各批次的數(shù)據(jù)為正態(tài)分布；

(2)各批的批內(nèi)方差相同；

(3)批平均值為正態(tài)分布。

對第一個假設(shè)可以通過擬合優(yōu)度檢驗驗證，方差相等通過Levene檢驗來驗證，該檢驗本質(zhì)上為F檢驗。

方差分析方法在MIL手冊和HB中都有詳細(xì)介紹，當(dāng)數(shù)據(jù)批次大于3時可以通過下式計算B基準(zhǔn)值

如果只有兩個批次的數(shù)據(jù)，則需要獲得更多批次的數(shù)據(jù)，或者合并這兩個批次的數(shù)據(jù)合并采用非結(jié)構(gòu)型數(shù)據(jù)方法計算B基準(zhǔn)值。

1.3.6 樣本離散系數(shù)變換方法

FAA的研究表明，復(fù)合材料性能的離散系數(shù)一般在4%～10%，數(shù)據(jù)離散系數(shù)低于4%是不正常的，可能造成B基準(zhǔn)值偏高，所以需要通過一個簡單的方法修正樣本變異性到任意所需值，變換公式如下：

其中：

通過上式可將離散系數(shù)為CV的樣本的離散系數(shù)轉(zhuǎn)換為CV*（FAA要求為4%）。程序中，讀入文件時可以通過CV值復(fù)選框來選擇是否對數(shù)據(jù)進(jìn)行CV值檢查與變換。

1.3.7 環(huán)境數(shù)據(jù)

HB 7618-2009對環(huán)境試驗數(shù)據(jù)進(jìn)行樣本合并后歸一化計算B基準(zhǔn)值，計算過程如圖1所示。歸一化樣本后，通過下式求得各環(huán)境下的B基準(zhǔn)值容差系數(shù)：

然后可以通過式（13）得到各環(huán)境下的B基準(zhǔn)值，即折減系數(shù)：

2 程序?qū)崿F(xiàn)

2.1 編程環(huán)境

軟件通過Microsoft Visual Studio環(huán)境下的C#語言進(jìn)行編寫，C#語言作為微軟在新一代開發(fā)平臺.NET推出的、完全面向?qū)ο蟮恼Z言，憑著其簡潔、高效、模板、標(biāo)準(zhǔn)化的特性，獲得了越來越多編程者的喜愛。

2.2 程序設(shè)計

程序通過C#實現(xiàn)界面化，程序界面如圖3所示。程序界面包括5部分，分別為數(shù)據(jù)讀取顯示、異常數(shù)據(jù)檢查、子體相容性檢驗、擬合優(yōu)度檢驗以及環(huán)境數(shù)據(jù)處理。程序在分析過程中可以通過改變設(shè)置分別按照MIL-HDBK-17-1F和HB 7618-2009兩種方法進(jìn)行B基準(zhǔn)值的分析計算。

程序通過txt文件導(dǎo)入試驗數(shù)據(jù)，分析結(jié)果可以通過txt文件給出計算報告，計算結(jié)果包括樣本基本統(tǒng)計量（如均值、離散系數(shù)等）、異常數(shù)據(jù)檢驗結(jié)果、子體相容性檢驗結(jié)果、B基準(zhǔn)計算方法和結(jié)果以及環(huán)境數(shù)據(jù)處理結(jié)果。進(jìn)行擬合優(yōu)度檢驗時，給出所有分布的統(tǒng)計量及OSL，可以通過不同要求自行判斷采用的B基準(zhǔn)值。

3 程序測試及結(jié)果輸出

圖3 B基準(zhǔn)值計算程序界面

在程序投入使用之前，需要對程序的準(zhǔn)確性和可靠性進(jìn)行驗證。本文中將程序計算結(jié)果與MIL-HDBK-17-1F及FAA/DOT/AR-03/19中算例對比來驗證程序。表2中給出了在單一環(huán)境下，程序計算結(jié)果與手冊中結(jié)果的對比，其誤差均小于0.1%，誤差主要來源于計算過程中的截斷誤差以及B基準(zhǔn)值單側(cè)容限系數(shù)近似公式。如果按照手冊中的要求，對計算結(jié)果按照原始數(shù)據(jù)的有效位數(shù)進(jìn)行取整，則誤差均為0。因此，可以證明在單一環(huán)境下數(shù)據(jù)程序計算結(jié)果與算例相同。

表2 單一環(huán)境算例計算結(jié)果對比

表3分別通過excel與統(tǒng)計程序計算FAA/DOT/AR-03/19中提供的環(huán)境數(shù)據(jù)算例，計算結(jié)果對比如表3所示。計算結(jié)果表明，各個環(huán)境下所得到數(shù)據(jù)誤差均在低于萬分之一，存在誤差為計算過程中截斷誤差不同造成的。對基準(zhǔn)值按照原始數(shù)據(jù)有效位數(shù)化整后兩種方法所提供的B基準(zhǔn)值完全相同。

通過以上算例對比可以證明，B基準(zhǔn)值計算程序得到的結(jié)果符合工程精度要求，軟件工作正常，滿足使用要求。

表3 不同環(huán)境基準(zhǔn)計算對比

4 結(jié)論

通過復(fù)合材料B基準(zhǔn)值計算程序的編制，可以通過用戶選擇實現(xiàn)按照航空行業(yè)標(biāo)準(zhǔn)HB 7618-2009或美國軍方手冊MIL-HDBK-17-1F中的方法來實現(xiàn)B基準(zhǔn)值的計算，包括數(shù)據(jù)統(tǒng)計特征計算（平均值、離散系數(shù)等）、異常數(shù)據(jù)檢查、子體相容性檢驗、基準(zhǔn)值計算、方差變換、環(huán)境數(shù)據(jù)處理等功能，計算結(jié)果可以通過報告形式給出。算例計算表明，計算程序準(zhǔn)確可靠，能夠滿足使用要求。

通過復(fù)合材料B基準(zhǔn)值計算程序計算B基準(zhǔn)值方便、快捷，對用戶統(tǒng)計學(xué)知識沒有很高要求，大大提高了復(fù)合材料B基準(zhǔn)值的計算效率。

[1] MIL-HDBK-17-1F,Composite Materials Handbook[S] Volume 1:Polymer Matrix Composites Guidelines for Characterization of Structural Materials.

[2] HB 7618-200 9,聚合物基復(fù)合材料力學(xué)性能數(shù)據(jù)表達(dá)準(zhǔn)則[S].

[3] FAA/DOT/AR-03/19,Material Qualification and Equivalency for Polymer Matrix Composite Material Systems:Updated[S].

[4] H.Press.Numerial Recipes in C:The Art of Scientific Computing:Second Edition[M].Cambridge University Press.2004.

[5] 孫堅石,葉強.復(fù)合材料力學(xué)性能數(shù)據(jù)B基準(zhǔn)值計算程序[J].航空制造技術(shù),2009,增刊:P19-24.

[6] 王翔,陳新文,王海鵬,等.基于統(tǒng)計的復(fù)合材料B基準(zhǔn)值計算方法研究[J].失效分析與預(yù)防,2010,5(4):210-215.

[7] 盛驟,謝式千,潘承毅.概率論與數(shù)理統(tǒng)計[M].高等教育出版社.2001.