基于蒙特卡洛的發(fā)動(dòng)機(jī)推力線橫移數(shù)值仿真

房 雷，張萬里

(1中國(guó)空空導(dǎo)彈研究院，河南洛陽(yáng) 471009;2航空制導(dǎo)武器航空科技重點(diǎn)實(shí)驗(yàn)室，河南洛陽(yáng) 471009)

0 引言

固體火箭發(fā)動(dòng)機(jī)的零部件在加工、裝配過程中，會(huì)產(chǎn)生尺寸和形位公差;在發(fā)動(dòng)機(jī)工作過程中，存在燃燒、燃?xì)饬鲃?dòng)以及噴管燒蝕的不均勻性，這些因素都不可避免的造成發(fā)動(dòng)機(jī)推力與彈軸形成偏移和夾角，影響導(dǎo)彈控制精度或火箭彈射程、橫偏[1]。

由于燃燒、燃?xì)饬鲃?dòng)以及噴管燒蝕的不均勻性難以控制和測(cè)量，目前在工程中一般都對(duì)上述因素進(jìn)行簡(jiǎn)化，采用噴管擴(kuò)張段軸線與發(fā)動(dòng)機(jī)軸線的平移和夾角來表征推力線橫移與偏斜，通過測(cè)量相關(guān)幾何參數(shù)，計(jì)算得到橫移與偏斜的具體結(jié)果。張春富[2]提出了一種激光跟蹤儀測(cè)量發(fā)動(dòng)機(jī)推力線的方法;來平安[3]在實(shí)測(cè)數(shù)據(jù)的基礎(chǔ)上，建立了一種推力線橫移和偏斜的模擬計(jì)算方法。

如果在發(fā)動(dòng)機(jī)設(shè)計(jì)階段即得到橫移與偏斜的大小，對(duì)于優(yōu)化結(jié)構(gòu)設(shè)計(jì)、合理分配零部件之間的公差、提高零件的加工工藝性有重要的意義。

發(fā)動(dòng)機(jī)推力線橫移一般由零部件的同軸度誤差和裝配間隙造成，對(duì)推力線橫移的研究，可轉(zhuǎn)化為對(duì)發(fā)動(dòng)機(jī)燃燒室和噴管等零部件公差設(shè)計(jì)的研究。傳統(tǒng)的公差設(shè)計(jì)有兩種方法，極值法和概率法[4]，由于零件尺寸公差和裝配方向在實(shí)際生產(chǎn)中都是隨機(jī)分布的，如果采用極值法進(jìn)行推力線橫移分析，不能反映實(shí)際橫移大小，很可能造成設(shè)計(jì)不滿足指標(biāo)要求的結(jié)果，從而需要提高設(shè)計(jì)精度，造成加工成本增加。

文中采用極坐標(biāo)法建立推力線橫移的數(shù)學(xué)模型，并用蒙特卡洛法進(jìn)行統(tǒng)計(jì)計(jì)算，在完成設(shè)計(jì)后即可得到推力線橫移的統(tǒng)計(jì)分析結(jié)果，為評(píng)估、優(yōu)化發(fā)動(dòng)機(jī)結(jié)構(gòu)設(shè)計(jì)、降低零件加工成本提供了依據(jù)。

1 數(shù)學(xué)模型

1.1 影響推力線橫移的因素分析

典型的發(fā)動(dòng)機(jī)噴管結(jié)構(gòu)如圖1所示，由圖可見影響推力線橫移的因素包括兩部分:在發(fā)動(dòng)機(jī)軸線和噴管擴(kuò)張段軸線之間，各配合面間安裝間隙的尺寸誤差以及各連接環(huán)節(jié)的同軸度誤差。本例中，影響推力線橫移的因素包括:殼體尾端配合面軸線與殼體軸線的同軸度誤差(因素1);殼體與噴管配合面之間的安裝間隙(因素2);擴(kuò)散段軸線與噴管配合面軸線同軸度誤差(因素3)。噴管結(jié)構(gòu)越復(fù)雜，在擴(kuò)散段和發(fā)動(dòng)機(jī)軸線之間的連接、裝配環(huán)節(jié)越多，影響推力線橫移的因素也就越多。

圖1 發(fā)動(dòng)機(jī)噴管結(jié)構(gòu)示意圖

一般在發(fā)動(dòng)機(jī)設(shè)計(jì)中，尺寸誤差只控制配合面實(shí)際尺寸的變化量在設(shè)計(jì)給定的極限尺寸之內(nèi)，不控制形位公差;而同軸度只控制被測(cè)要素的形位公差，與實(shí)際尺寸無關(guān)。因此，對(duì)于這兩類誤差，屬于按獨(dú)立原則設(shè)計(jì)的零件要素，在建立推力線橫移的尺寸鏈過程中，均可作為組成環(huán)計(jì)入。

1.2 推力線橫移的數(shù)學(xué)模型

以圖1所示的典型發(fā)動(dòng)機(jī)噴管結(jié)構(gòu)為例，在不考慮誤差引起推力線偏斜的情況下，影響橫移的尺寸鏈關(guān)系可用圖2表示。

O點(diǎn)表示發(fā)動(dòng)機(jī)軸線的投影，A點(diǎn)表示燃燒室后端配合面軸線的投影，B點(diǎn)表示噴管配合面軸線的投影，C點(diǎn)表示擴(kuò)散段軸線的投影。用線段 OA、AB、BC分別代表因素1、因素2、因素3引起的橫移大小，發(fā)動(dòng)機(jī)推力線橫移大小用線段OC表征?？梢杂孟蛄康姆椒ū硎綩C的大小以及相對(duì)于發(fā)動(dòng)機(jī)軸線的方位。

由圖2，將OA、AB、BC用3個(gè)向量表示，即:

向量 OA=(ρ1cosθ1，ρ1sinθ1)

向量 AB=(ρ2cosθ2，ρ2sinθ2)

向量 BC=(ρ3cosθ3，ρ3sinθ3)

圖2 影響推力線橫移的因素與發(fā)動(dòng)機(jī)軸線間的關(guān)系

式中:ρ1、ρ3表示同軸度公差產(chǎn)生的橫移的大小，ρ2為尺寸誤差產(chǎn)生橫移的大小，θ1、θ2、θ3分別表示不同誤差產(chǎn)生的橫移相對(duì)于基準(zhǔn)的方位角。

向量OC=OA+AB+BC=

設(shè)OC(ρ，θ)，C點(diǎn)在以 O點(diǎn)為極點(diǎn)的極坐標(biāo)下，則有:

根據(jù)上式可確定ρ，取ρ≥0。

若 ρ1cosθ1+ ρ2cosθ2+ ρ3cosθ3= 0，ρ1sinθ1+ρ2sinθ2+ ρ3sinθ3＞ 0，則 θ = π/2，反之，ρ1sinθ1+ρ2sinθ2+ ρ3sinθ3＜ 0，則 θ =3π/2。

若 ρ1cosθ1+ ρ2cosθ2+ ρ3cosθ3≠ 0，則:

由式(2)、式(3)可計(jì)算向量 OC的 ρ、θ，得到推力線橫移的大小、方位。同理，如果有更多因素影響推力線橫移，可參照上述方法進(jìn)行計(jì)算。

1.3 參數(shù)的取值及分布形式

得到了推力線橫移表達(dá)式，只要知道 ρ1、ρ2、ρ3、θ1、θ2、θ3等參數(shù)值，即可計(jì)算出橫移大小和方向角。

在本例中，ρ1、ρ3為同軸度誤差，對(duì)于同軸度這種上、下偏差為對(duì)稱分布的形位公差，無論是增環(huán)還是減環(huán)，對(duì)封閉環(huán)的影響都是相同的，在計(jì)算過程中，它的大小表示為0±T/2(其中T為同軸度公差值)，概率分布為瑞利分布。ρ2為尺寸公差，由設(shè)計(jì)公差確定，取值范圍為裝配間隙的下偏差 /2～上偏差 /2，概率分布為正態(tài)分布。θ1、θ2、θ3為裝配方位誤差，概率分布為均勻分布，取值范圍0～2π。

2 計(jì)算方法

蒙特卡洛(Monte Carlo)法[5]亦稱為隨機(jī)模擬(randomation)方法，有時(shí)也稱作隨機(jī)抽樣(random sampling)技巧或統(tǒng)計(jì)試驗(yàn)(statistical testing)方法，是以概率統(tǒng)計(jì)理論為基礎(chǔ)的一種方法。

它的基本思想是為了求解數(shù)學(xué)、物理、工程技術(shù)以及生產(chǎn)管理等方面的問題，首先建立概率模型或隨機(jī)過程，使它的參數(shù)等于問題的解;然后通過對(duì)模型或過程的觀察或抽樣試驗(yàn)來計(jì)算所求參數(shù)的特征，最后給出所求解的近似值，而解的精確度可用估計(jì)標(biāo)準(zhǔn)誤差來表示。

在應(yīng)用蒙特卡洛方法進(jìn)行公差分析的過程中，一般有如下的步驟[4]:

1)確定封閉環(huán)與組成環(huán)的尺寸、公差及其分布特征和傳遞系數(shù)，并確定分布函數(shù);

2)根據(jù)模擬精度要求確定模擬次數(shù)N;

3)根據(jù)各組成環(huán)的分布規(guī)律進(jìn)行隨機(jī)抽樣，從而得到各組成環(huán)尺寸的隨機(jī)數(shù)列;

4)在得到各組成環(huán)的模擬值后，代入到尺寸鏈方程中，計(jì)算得到封閉環(huán)的數(shù)值;

5)重復(fù)3)、4)步驟N次，即可得到封閉環(huán)尺寸的N個(gè)子樣。

文中采用MATLAB軟件進(jìn)行編程，計(jì)算推力線橫移大小、方位角，通過程序?qū)崿F(xiàn)蒙特卡洛模擬計(jì)算，隨機(jī)數(shù)由MATLAB軟件中相應(yīng)的函數(shù)產(chǎn)生，分別由normrnd、raylrnd、unifrnd函數(shù)表征正態(tài)分布、瑞利分布、均勻分布。

3 算例

某發(fā)動(dòng)機(jī)設(shè)計(jì)為:殼體尾部配合面與殼體軸線的同軸度為0.05mm;殼體尾端與噴管的安裝間隙為0.043～0.293mm;擴(kuò)散段軸線與噴管配合面軸線的同軸度為0.2mm，采用MATLAB編程計(jì)算，相關(guān)計(jì)算參數(shù)見表1。計(jì)算次數(shù)為200000次，統(tǒng)計(jì)結(jié)果見圖3、圖 4。

表1 參數(shù)設(shè)置

圖3 推力線橫移大小仿真統(tǒng)計(jì)結(jié)果

圖4 推力線橫移方位角仿真統(tǒng)計(jì)結(jié)果

由圖可見，橫移大小呈瑞利分布，方位角呈均勻分布。在置信水平0.9973時(shí)，該發(fā)動(dòng)機(jī)推力線橫移不大于0.189mm，小于極值法計(jì)算的0.2715mm。與極值法相比，蒙特卡洛法的計(jì)算結(jié)果更符合實(shí)際加工、裝配情況;在設(shè)計(jì)要求確定時(shí)，通過計(jì)算可以方便的評(píng)估設(shè)計(jì)結(jié)果，也可以根據(jù)仿真結(jié)果進(jìn)行公差分配，適當(dāng)增大組成環(huán)公差，降低加工成本。

4 結(jié)論

通過對(duì)影響固體火箭發(fā)動(dòng)機(jī)推力線橫移的因素進(jìn)行分析，建立了推力線橫移的數(shù)學(xué)模型，給出了基于極坐標(biāo)法的數(shù)學(xué)表達(dá)式，采用蒙特卡洛法對(duì)推力線橫移進(jìn)行了仿真計(jì)算。結(jié)果表明，采用該方法能夠方便、迅速完成推力線橫移的計(jì)算，可以在滿足設(shè)計(jì)要求的前提下使零件的加工誤差得到合理分配，為在發(fā)動(dòng)機(jī)設(shè)計(jì)階段優(yōu)化結(jié)構(gòu)設(shè)計(jì)參數(shù)、評(píng)估設(shè)計(jì)結(jié)果提供依據(jù)。

[1]劉文，張衛(wèi)華，周珞晶，等.固體火箭發(fā)動(dòng)機(jī)推力偏斜和橫移對(duì)自旋火箭彈射程和橫偏影響[J].彈箭與制導(dǎo)學(xué)報(bào)，1999，19(2):35 －38.

[2]張春富，唐文彥，李慧鵬，等.激光跟蹤儀在固體火箭發(fā)動(dòng)機(jī)推力線測(cè)量中的應(yīng)用[J].固體火箭技術(shù)，2007，30(6):548－551.

[3]來平安，田維平，余貞勇，等.固體火箭發(fā)動(dòng)機(jī)推力線橫移和偏斜模擬計(jì)算[J].固體火箭技術(shù)，1999，22(1):35－38.

[4]王晶，石宏，黃笑飛，等.基于蒙特卡羅模擬法的航空發(fā)動(dòng)機(jī)裝配公差分析[J].沈陽(yáng)航空工業(yè)學(xué)院學(xué)報(bào)，2010，27(4):8－11.

[5]徐鐘濟(jì).蒙特卡羅方法[M].上海:上?？茖W(xué)技術(shù)出版社，1985.