孫傳智，李愛(ài)群，繆長(zhǎng)青，黎少華，喬 燕

（1.東南大學(xué) 混凝土及預(yù)應(yīng)力混凝土結(jié)構(gòu)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210096；2.宿遷學(xué)院 建筑工程系，江蘇 宿遷 223800）

粘滯阻尼器是目前高層建筑結(jié)構(gòu)振動(dòng)控制中應(yīng)用較多的結(jié)構(gòu)被動(dòng)控制裝置之一。在進(jìn)行粘滯阻尼器減震設(shè)計(jì)時(shí)，可以通過(guò)調(diào)整阻尼系數(shù)和速度指數(shù)取得不同的參數(shù)組合達(dá)到同樣的減震效果，但是不同的參數(shù)取值，阻尼器提供的阻尼力相差較大，從而對(duì)與其連接的柱內(nèi)力影響較大，同時(shí)阻尼器本身也會(huì)出現(xiàn)磨損或損壞，致使減震設(shè)計(jì)存在不能達(dá)到減震要求的隱患［1］。同時(shí)粘滯阻尼器的價(jià)格主要取決于阻尼力大小，因此在進(jìn)行粘滯阻尼器減震結(jié)構(gòu)設(shè)計(jì)時(shí)通過(guò)粘滯阻尼器參數(shù)優(yōu)化分析，既能保證結(jié)構(gòu)安全，又能盡量降低建造成本，這對(duì)于粘滯阻尼器的使用推廣尤為重要。目前已對(duì)粘彈性阻尼器的參數(shù)優(yōu)化進(jìn)行了諸多研究［2－5］，而粘滯阻尼器參數(shù)優(yōu)化研究還比較少，目前多依靠經(jīng)驗(yàn)和試算來(lái)確定粘滯阻尼器參數(shù)，工作量大。響應(yīng)面分析方法是一種基于試驗(yàn)設(shè)計(jì)理論的近似方法，最早由Box等［7］提出，其基本思想是對(duì)樣本點(diǎn)處的函數(shù)響應(yīng)計(jì)算值或試驗(yàn)值，用回歸分析法構(gòu)造一個(gè)具有明確表達(dá)形式的多項(xiàng)式來(lái)表達(dá)隱式功能函數(shù)。近年來(lái)，響應(yīng)面分析方法已被應(yīng)用到結(jié)構(gòu)的可靠度評(píng)估、模型修正和結(jié)構(gòu)優(yōu)化［8－12］等方面，其中武和全等［12］將響應(yīng)面法與試驗(yàn)設(shè)計(jì)、有限元分析計(jì)算等結(jié)合起來(lái)，對(duì)S型薄壁梁的結(jié)構(gòu)進(jìn)行抗撞性優(yōu)化設(shè)計(jì)。筆者將響應(yīng)面法與有限元分析相結(jié)合，構(gòu)造參數(shù)優(yōu)化目標(biāo)函數(shù)和約束條件，求得非線性粘滯阻尼器參數(shù)最優(yōu)解。算例應(yīng)用表明該方法在減震設(shè)計(jì)過(guò)程中具有較強(qiáng)的實(shí)用性和有效性。

1 粘滯阻尼器計(jì)算模型及減震結(jié)構(gòu)非線性分析

1.1 粘滯阻尼器計(jì)算模型

國(guó)內(nèi)外學(xué)者對(duì)粘滯阻尼器力學(xué)模型進(jìn)行了大量研究［13－15］，美國(guó)Taylor公司給出的阻尼力通用表達(dá)式為：

1.2 減震結(jié)構(gòu)非線性分析

未設(shè)置粘滯阻尼器時(shí)，原結(jié)構(gòu)的動(dòng)力平衡方程為：

式中：［M］為原結(jié)構(gòu)質(zhì)量矩陣；［C］為原結(jié)構(gòu)阻尼矩陣，采用Rayleigh阻尼；［K］為原結(jié)構(gòu)剛度矩陣；［I］為地震動(dòng)作用的位置向量；｛x｝、｝分別為節(jié)點(diǎn)的位移向量、速度向量和加速度向量；為地震動(dòng)加速度。當(dāng)結(jié)構(gòu)采用附加方式設(shè)置粘滯阻尼器時(shí)，結(jié)構(gòu)的動(dòng)力平衡方程變?yōu)椋?/p>

式中 ［Cd］為阻尼器提供的附加阻尼矩陣。

求解附加粘滯阻尼器減震結(jié)構(gòu)的控制方程時(shí)，出于研究目的和提高計(jì)算效率的考慮，可以只進(jìn)行彈性時(shí)程分析，采用快速非線性分析法進(jìn)行分析［17］。該方法是一種將模態(tài)疊加法和增量法組合使用來(lái)快速求解只含有少量非線性單元的分析方法，該方法可在保證計(jì)算精度的前提下快速提高計(jì)算速度［18］，通用有限元軟件SAP2000可以實(shí)現(xiàn)快速非線性分析［19］。

2 基于響應(yīng)面法的減震結(jié)構(gòu)粘滯阻尼器參數(shù)優(yōu)化方法

2.1 減震結(jié)構(gòu)粘滯阻尼器參數(shù)優(yōu)化數(shù)學(xué)模型

如前所述，在保證結(jié)構(gòu)安全的前提下，應(yīng)盡量控制附加阻尼器提供的阻尼力，這樣可以降低結(jié)構(gòu)的建造成本，所以在優(yōu)化設(shè)計(jì)時(shí)以所有粘滯阻尼器提供的阻尼力最小為目標(biāo)函數(shù)，以滿足結(jié)構(gòu)正常使用時(shí)的層間位移限值作為約束條件，進(jìn)行粘滯阻尼器的參數(shù)優(yōu)化，具體形式如下：

設(shè)計(jì)變量：求粘滯阻尼器阻尼系數(shù)Cdi和速度指數(shù)αi。

目標(biāo)函數(shù)：所有粘滯阻尼器提供的阻尼力之和∑Fi最小，即min（∑Fi）。

約束條件：各層層間最大位移小于限值，即max（│Δi│）＜［Δ］，［Δ］＝5mm；速度指數(shù)α滿足0.2≤αi≤0.8。

2.2 基于響應(yīng)面法減震結(jié)構(gòu)粘滯阻尼器參數(shù)優(yōu)化步驟

基于響應(yīng)面進(jìn)行減震結(jié)構(gòu)粘滯阻尼器參數(shù)優(yōu)化，首先利用響應(yīng)面法和有限元分析相結(jié)合得到各層層間位移和各阻尼器阻尼力的響應(yīng)面函數(shù)，然后進(jìn)行參數(shù)優(yōu)化。具體來(lái)說(shuō)分為6個(gè)步驟：

1）試驗(yàn)設(shè)計(jì)。為了確保響應(yīng)面函數(shù)精度，需要以概率論、數(shù)理統(tǒng)計(jì)和線性代數(shù)等為理論基礎(chǔ)進(jìn)行試驗(yàn)設(shè)計(jì)，常用的試驗(yàn)設(shè)計(jì)方法有中心復(fù)合設(shè)計(jì)方法和Box－Behnken矩陣抽樣設(shè)計(jì)方法。根據(jù)影響減震效果的主要因素，選擇粘滯阻尼器阻尼系數(shù)和速度指數(shù)作為自變量xi（i＝1，2，…，L），各阻尼器所提供的最大阻尼力和層間位移作為因變量ym（m＝1，2，…，n），然后采用 Box－Behnken矩陣抽樣法進(jìn)行試驗(yàn)設(shè)計(jì)，確定K組樣本點(diǎn)。

2）有限元計(jì)算分析。根據(jù)試驗(yàn)設(shè)計(jì)的樣本點(diǎn)參數(shù)值，利用有限元分析軟件SAP2000進(jìn)行計(jì)算，得到K組層間位移和各阻尼器最大阻尼力變化值ym。

3）參數(shù)篩選。在研究的初始階段有可能考慮了所有的參數(shù)，如果參數(shù)很多，則需要進(jìn)行參數(shù)篩選。參數(shù)篩選主要采用方差分析方法，其基本思想是將總變異平方和分解為試驗(yàn)因素效應(yīng)與隨機(jī)誤差所引起的變異平方和，由此構(gòu)造出方差分析的F檢驗(yàn)統(tǒng)計(jì)量［20］，找出顯著性參數(shù)。

4）響應(yīng)面擬合。將K組自變量及其對(duì)應(yīng)的K組因變量代入式（4），采用最小二乘法估計(jì)多項(xiàng)式系數(shù)，得到響應(yīng)面模型函數(shù)。

5）響應(yīng)面函數(shù)精度檢驗(yàn)。為了驗(yàn)證響應(yīng)面模型是否可靠，需對(duì)其進(jìn)行精度檢驗(yàn)，可采用多重?cái)M合系數(shù)R2進(jìn)行響應(yīng)面精度檢驗(yàn)，如式（5）所示，式中代表響應(yīng)面模型計(jì)算值，yj代表有限元模型計(jì)算值，代表有限元模型計(jì)算值的平均值。R2判定系數(shù)在0～1之間取值，R2值越大，則回歸模型就越接近實(shí)際情況，如果精度較低則回到第1）步重新進(jìn)行試驗(yàn)設(shè)計(jì)。

6）粘滯阻尼器參數(shù)優(yōu)化。以層間位移小于限值作為約束條件，所有粘滯阻尼器提供的阻尼力之和最小作為目標(biāo)函數(shù)，建立粘滯阻尼器參數(shù)優(yōu)化數(shù)學(xué)模型，運(yùn)用非線性規(guī)劃優(yōu)化方法進(jìn)行參數(shù)優(yōu)化。

3 算例分析

3.1 算例概況

鋼筋混凝土框架結(jié)構(gòu)，5層，層高5m，跨度6m，柱截面為400mm×400mm，梁截面為250mm×500mm。梁柱材料均為C35混凝土，每層梁柱節(jié)點(diǎn)處質(zhì)量為50kN，輸入地震波為EL－CENTRO，地面運(yùn)動(dòng)加速度峰值為140cm／s2，層間位移最大限值為5mm。有限元模型如圖1所示。

圖1 算例有限元模型

3.2 層間位移和阻尼力之和響應(yīng)面函數(shù)

采用Box－Behnken樣本抽樣法進(jìn)行試驗(yàn)設(shè)計(jì)，1～5號(hào)阻尼器的阻尼系數(shù)為xi（i＝1、3、5、7、9），速度指數(shù)為xi（i＝2、4、6、8、10）。各層層間位移為Ri（i＝1、2、3、4、5），阻尼力之和為R6，得到165組樣本點(diǎn)。

應(yīng)用數(shù)理統(tǒng)計(jì)的F檢驗(yàn)法分析所選參數(shù)對(duì)特征頻率的顯著性，計(jì)算各參數(shù)的統(tǒng)計(jì)特征量的顯著性水平P值，當(dāng)P＜0.05時(shí)，則該參數(shù)為顯著。表1為顯著性水平為顯著的參數(shù)。由表1可以看出，參數(shù)x1、x4、x6和x8對(duì)所有因變量的影響都顯著，參數(shù)x2和x3對(duì)因變量R1、R2、R3和R6的影響顯著，參數(shù)x7除了對(duì)因變量R2影響不顯著外，對(duì)其它因變量均顯著，參數(shù)x9僅對(duì)因變量R5和R6影響顯著，參數(shù)x10對(duì)因變量R3、R5和R6影響顯著。此外，部分交叉項(xiàng)和二次項(xiàng)對(duì)該模型的顯著性較參數(shù)的影響要小，但為了保證精度，在對(duì)試驗(yàn)設(shè)計(jì)計(jì)算得到的樣本數(shù)據(jù)進(jìn)行二次多項(xiàng)式響應(yīng)面回歸時(shí)，采用完全二次多項(xiàng)式函數(shù)模型。

表1 顯著性檢驗(yàn)結(jié)果

采用完全二次多項(xiàng)式對(duì)該樣本值進(jìn)行響應(yīng)面回歸即可獲得各層層間位移及各阻尼器阻尼力與各參數(shù)之間的響應(yīng)面函數(shù)模型，通過(guò)對(duì)各響應(yīng)面模型進(jìn)行多重?cái)M合系數(shù)R2判定，多重?cái)M合系數(shù)R2分別為0.9843、0.9756、0.9867、0.9621、0.9513 和0.9769，均大于0.95，能夠滿足參數(shù)優(yōu)化的精度要求。圖2所示R1殘差正態(tài)分布概率圖，各殘差點(diǎn)主體部分呈直線，說(shuō)明誤差呈正態(tài)分布。以上分析說(shuō)明模型擬合度好，利用響應(yīng)面函數(shù)模型精確能夠反映各層層間位移、粘滯阻尼器阻尼力之和與各阻尼器的阻尼系數(shù)、速度指數(shù)之間的關(guān)系，從而能夠保證后續(xù)粘滯阻尼器參數(shù)優(yōu)化的精度。

圖2 R1殘差正態(tài)分布概率圖

3.3 粘滯阻尼器參數(shù)影響分析

限于篇幅，僅以R1為例利用響應(yīng)面函數(shù)分析粘滯阻尼器參數(shù)對(duì)結(jié)構(gòu)的影響。圖3和圖4分別為1號(hào)阻尼器和3號(hào)阻尼器對(duì)R1的影響曲線圖，即參數(shù)x1和x2、x5和x6對(duì)因變量R1的影響曲線圖，對(duì)于1號(hào)阻尼器對(duì)R1的影響，從圖3可以看出，R1是隨著阻尼系數(shù)的增大而減小，隨著速度指數(shù)的減小而減小。但是，從圖4可以看出3號(hào)阻尼器對(duì)的R1影響關(guān)系復(fù)雜，在速度指數(shù)等于0.2時(shí)，R1隨著阻尼系數(shù)的增大而減小；在速度指數(shù)等于0.8時(shí)，R1隨著阻尼系數(shù)的增大而增大，反之，在阻尼系數(shù)等于30kN·s／m時(shí)，R1隨著阻尼系數(shù)的增大而增大；在阻尼系數(shù)等于150kN·s／m時(shí)，R1隨著速度指數(shù)的增大而減小。說(shuō)明如果3號(hào)阻尼器參數(shù)選擇不合適，有可能使其對(duì)R1有放大作用。

圖3 x1和x2對(duì)R1的交互影響

圖4 x5和x6對(duì)R1的交互影響

由以上分析可知，各阻尼器之間、阻尼器與結(jié)構(gòu)之間的影響關(guān)系比較復(fù)雜。阻尼器阻尼系數(shù)與速度指數(shù)的不同取值組合，對(duì)減震結(jié)構(gòu)的減震效果有個(gè)最優(yōu)解的問(wèn)題，需對(duì)阻尼器的參數(shù)進(jìn)行優(yōu)化設(shè)計(jì)，即能保證結(jié)構(gòu)安全，又使得建造成本最低。

3.4 粘滯阻尼器參數(shù)優(yōu)化

當(dāng)確定了響應(yīng)面模型函數(shù)后，以各阻尼力之和最小為目標(biāo)函數(shù)，各層間位移和各參數(shù)的取值范圍作為約束條件，建立參數(shù)優(yōu)化模型，運(yùn)用非線性規(guī)劃優(yōu)化方法進(jìn)行參數(shù)優(yōu)化，即可得到非線性粘滯阻尼器的優(yōu)化參數(shù)，優(yōu)化后的非線性粘滯阻尼器的參數(shù)如表2所示。

通過(guò)有限元分析和上文所得的響應(yīng)面函數(shù)可得1～5層層間位移和各粘滯阻尼器阻尼力之和，如表3所示。響應(yīng)面函數(shù)預(yù)測(cè)所得的層間位移和各阻尼器阻尼力之和與有限元分析所得的相差較小，能夠滿足工程精度要求。

表2 非線性粘滯阻尼器參數(shù)優(yōu)化值

表3 響應(yīng)面函數(shù)預(yù)測(cè)和有限元分析比較

圖5、6為無(wú)控結(jié)構(gòu)和減震結(jié)構(gòu)節(jié)點(diǎn)6的位移時(shí)程曲線圖和加速度時(shí)程曲線圖，從圖中可以看出，結(jié)構(gòu)由于設(shè)置了粘滯阻尼器，取得了較好的減震效果。

圖5 節(jié)點(diǎn)6位移時(shí)程比較

圖6 節(jié)點(diǎn)6加速度時(shí)程比較

對(duì)于非線性粘滯阻尼器阻尼系數(shù)和速度指數(shù)的取值，選擇如表4所示的5種不同方案進(jìn)行減震效果比較。第1種方案是基于響應(yīng)面法優(yōu)化結(jié)果；第2種方案是阻尼系數(shù)取值150kN·s／m，速度指數(shù)為0.2；第3種方案是阻尼系數(shù)取值150kN·s／m，速度指數(shù)為0.8；第4種方案是阻尼系數(shù)取值90kN·s／m，速度指數(shù)為0.5；第5種方案是阻尼系數(shù)取值30kN·s／m，速度指數(shù)為0.2。

表5為5種不同方案下的參數(shù)取值所得的層間位移和各粘滯阻尼器阻尼力之和?？梢钥闯?，方案2雖然能夠滿足層間位移限值的要求，但是各阻尼器之和太大，建造成本高；方案3、4、5雖然各阻尼器阻尼力之和較小，但是不能滿足層間位移限值的要求。

4 結(jié) 語(yǔ)

1）以一榀鋼筋混凝土框架結(jié)構(gòu)為例，進(jìn)行了基于響應(yīng)面法的非線性粘滯阻尼器參數(shù)優(yōu)化設(shè)計(jì)。結(jié)果表明采用簡(jiǎn)單的代數(shù)公式能夠精確擬合設(shè)計(jì)變量和響應(yīng)量之間的復(fù)雜關(guān)系，把復(fù)雜問(wèn)題顯式化，并且響應(yīng)面函數(shù)精度檢驗(yàn)表明，響應(yīng)面函數(shù)精度高，能夠滿足優(yōu)化精度要求。

表4 不同粘滯阻尼器參數(shù)取值方案

表5 不同粘滯阻尼器參數(shù)取值時(shí)的減震效果比較

2）選擇阻尼系數(shù)和速度指數(shù)作為自變量，結(jié)構(gòu)層間位移、阻尼力之和作為響應(yīng)值，進(jìn)行響應(yīng)面函數(shù)擬合。研究表明層間位移和粘滯阻尼力之和與非線性粘滯阻尼器參數(shù)取值呈現(xiàn)復(fù)雜的非線性關(guān)系。

3）結(jié)合響應(yīng)面函數(shù)模型，建立目標(biāo)函數(shù)和約束條件，進(jìn)行非線性粘滯阻尼器參數(shù)優(yōu)化，得到了粘滯阻尼器的參數(shù)優(yōu)化值。該方法精度高、計(jì)算量小，過(guò)程思路清晰，簡(jiǎn)單易行，避免了憑經(jīng)驗(yàn)進(jìn)行參數(shù)調(diào)整過(guò)程中的反復(fù)迭代試算。

［1］陳瑜，劉偉慶，陸偉東.線性與非線性黏滯阻尼支撐減震效果的比較［J］.南京工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，31（6）：45－50.Chen Y，Liu W Q，Lu W D.Comparison of seismicreduced effect between linear and nonlinear viscous damper bracings ［J］.Journal of Nanjing University of Technology：Natural Science Edition，2009，31（6）：45－50.

［2］李正良，李疏影，孫毅.基于遺傳算法的阻尼支撐優(yōu)化分析［J］.重慶建筑大學(xué)學(xué)報(bào)，2007，29（3）：57－60.Li Z L，Li S Y，Sun Y.Optimum analysis of damper braces based on genetic algorithm ［J］.Journal of Chongqing Jianzhu University，2007，29（3）：57－60.

［3］涂青，徐趙東，彭軍.隔減震結(jié)構(gòu)中黏彈性阻尼裝置的遺傳算法優(yōu)化分析［J］.東南大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，39（1）：73－77.Tu Q，Xu Z D，Peng J.Parametric optimization of viscoelastic device in earthquake isolation and mitigation of structures ［J］.Journal of Southeast University：Natural Science Edition，2009，39（1）：73－77.

［4］Xu Z D，Shen Y P，Zhao H T.A synthetic optimization analysis method on structures with viscoelastic dampers［J］.Soil Dynamics and Earthquake Engineering，2003，23：683－689.

［5］郭勇，孫炳楠，葉尹.多目標(biāo)優(yōu)化方法在輸電塔阻尼器布置中的應(yīng)用［J］.浙江大學(xué)學(xué)報(bào)：工學(xué)版，2006，40（10）：1755－1760.Guo Y，Sun B N，Ye Y.Multi－objective optimization approach for damper distribution of transmission towers［J］.Journal of Zhejiang University： Engineering Science，2006，40（10）：1755－1760.

［6］Tandjiria V，Teh C I，Low B K.Reliability analysis of laterally loaded piles using response surface methods［J］.Structural Safety，2000，22：335－355.

［7］Box G E P，Wilson K B.On the experimental attainment of optimum conditions［J］.Journal of the Royal Statistical Society，1951，13（1）：1－34.

［8］Lee S Y，Haldar A.Reliability of frame and shear wall structural systems II：Dynamic loading［J］.Journal of Structural Engineering，2003，129（2）：233－240.

［9］Deng L，Cai C S.Bridge model updating using response surface method and genetic algorithm ［J］.Journal of Bridge Engineering，2010，15（5）：553－564.

［10］ Bayramov F， Tasdemir C， Tasdemir M A.Optimisation of steel fiber reinforced concretes by means of statistical response surface method ［J］.Cement ＆Concrete Composites，2004，26：665－675.

［11］費(fèi)慶國(guó)，韓曉林，蘇鶴玲.響應(yīng)面有限元模型修正的實(shí)現(xiàn)與應(yīng)用［J］.振動(dòng)、測(cè)試與診斷，2010，30（2）：132－134.Fei Q G，Han X L，Su H L.Response surface based finite element model updating and its application ［J］.Journal of Vibration， Measurement and Diagnosis，2010，30（2）：132－134.

［12］武和全，楊家興，辛勇.基于響應(yīng)面法的S型薄壁梁抗撞性優(yōu)化設(shè)計(jì)［J］.機(jī)械科學(xué)與技術(shù)，2010，29（9）：1132－1136.Wu H Q，Yang J X，Xin Y.Optimal design of S－shaped rail for crashworthiness analysis based on response surface method ［J］.Mechanical Science and Technology for Aerospace Engineering，2010，29（9）：1132－1136.

［13］Makris N，Constantinou M C.Fractional－derivative Maxwell model for viscous dampers ［J］.Journal of Structural Engineering，1991，117（9）：2708－2724.

［14］葉正強(qiáng)，李愛(ài)群，程文瀼，等.采用粘滯流體阻尼器的工程結(jié)構(gòu)減振設(shè)計(jì)研究［J］.建筑結(jié)構(gòu)學(xué)報(bào)，2001，22（4）：61－66.Ye Z Q，Li A Q，Cheng W R，et al.Study on vibration energy dissipation design of structures with fluid viscous dampers ［J］.Journal of Building Structures，2001，22（4）：61－66.

［15］葉正強(qiáng)，李愛(ài)群，徐幼麟.工程結(jié)構(gòu)粘滯流體阻尼器減振新技術(shù)及其應(yīng)用［J］.東南大學(xué)學(xué)報(bào)：自然科學(xué)版，2002，32（3）：466－473.Ye Z Q，Li A Q，Xu Y L.Fluid viscous damper technology and its engineering application for structural vibration energy dissipation ［J］.Journal of Southeast University：Natural Science Edition，2002，32 （3）：466－473.

［16］Housner G W，Bergrnan L A.Structural control：past，present and future ［J］.Journal of Engineering Mechanics，1997，123（9）：73－79

［17］陳學(xué)偉，韓小雷，毛貴牛，等.粘滯阻尼器在連體高層結(jié)構(gòu)中的抗風(fēng)減振效果［J］.土木建筑與環(huán)境工程，2009，31（5）：74－80.Chen X W，Han X L，Mao G N，et al.Fluid viscous damper technology and its engineering application for structural vibration energy dissipation［J］.Journal of Civil，Architectural ＆ Environmental Engineering，2009，31（5）：74－80.

［18］Wilson E L.Three dimensional dynamic analysis for structures with emphasis on earthquake engineering［M］.Berkeley：Computers and Structures，Inc.，1997.

［19］彭俊生，羅永坤，彭地.結(jié)構(gòu)動(dòng)力學(xué)、抗震計(jì)算與SAP2000應(yīng)用［M］.成都：西南交通大學(xué)出版社，2007.

［20］呂棟雷，曹志耀，鄧寶，等.利用方差分析法進(jìn)行模型驗(yàn)證［J］.計(jì)算機(jī)仿真，2006，23（8）：46－48.Lyu D L，Cao Z Y，Deng B，et al.Model verification using variance analysis method ［J］. Computer Simulation，2006，23（8）：46－48.