陳爾康, 廖欣, 高長生, 荊武興
(1.哈爾濱工業(yè)大學(xué) 航天學(xué)院, 黑龍江 哈爾濱 150001; 2.上海機(jī)電工程研究所, 上海 201109)
旋轉(zhuǎn)導(dǎo)彈在飛行過程中繞自身縱軸低速自旋,由此可實(shí)現(xiàn)單通道控制,以簡化控制系統(tǒng),并減小發(fā)動機(jī)推力偏心、質(zhì)量偏心和外形工藝誤差等干擾的影響[1]。但彈體自旋也使得旋轉(zhuǎn)導(dǎo)彈的縱軸繞速度矢量作錐形運(yùn)動。由于不收斂的錐形運(yùn)動會降低導(dǎo)彈射程和精度[2],其穩(wěn)定性得到了眾多專家和學(xué)者的關(guān)注,目前已建立起基于剛體假設(shè)的較完備理論[3]。但為了提高速度、機(jī)動性和射程[4],旋轉(zhuǎn)導(dǎo)彈大多數(shù)具備較大的長徑比和推重比,導(dǎo)致彈體具有較大的彈性,從而引起了彈性變形自由度和剛體運(yùn)動自由度之間的耦合(簡稱剛彈耦合)[5]。文獻(xiàn)[5-6]的研究表明剛彈耦合效應(yīng)對旋轉(zhuǎn)導(dǎo)彈的軌跡和錐形運(yùn)動有較大影響。因此,剛彈耦合效應(yīng)是旋轉(zhuǎn)導(dǎo)彈研究中必須考慮的因素,而傳統(tǒng)的氣動彈性[7]和飛行力學(xué)研究并不涉及這類問題。
飛行動力學(xué)建模是彈性旋轉(zhuǎn)導(dǎo)彈動力學(xué)和制導(dǎo)控制問題研究的基礎(chǔ)[8]。一些學(xué)者在這方面做了初步研究,建立了考慮彈性的姿態(tài)動力學(xué)方程,以考察彈性對姿態(tài)運(yùn)動特性的影響[9-10],但未建立完整的飛行動力學(xué)模型。事實(shí)上,彈性旋轉(zhuǎn)導(dǎo)彈的運(yùn)動包括變化幅度較小的彈性變形自由度和變化幅度較大的剛體運(yùn)動自由度,因此應(yīng)選擇能準(zhǔn)確描述這兩種運(yùn)動自由度的坐標(biāo)系來完成建模。為研究方便,通常選擇隨導(dǎo)彈運(yùn)動的體坐標(biāo)系為參考系[11],常用的主要有平均軸系和準(zhǔn)坐標(biāo)系。平均軸系又稱為蒂塞朗(Tisserand)坐標(biāo)系,導(dǎo)彈相對于該坐標(biāo)系的線動量和角動量始終為0,實(shí)現(xiàn)了剛體運(yùn)動自由度和彈性變形自由度間的解耦[12-13]。文獻(xiàn)[14]在平均軸系下利用拉格朗日方程推導(dǎo)了彈性飛行器飛行動力學(xué)模型,該模型假設(shè)慣量張量不受彈性變形的影響,因此形式較為簡單,廣泛應(yīng)用于各類彈性飛行器的研究[15-18]。雖然在平均軸系下建模能夠簡化動力學(xué)方程,但其定義中的限制條件在真實(shí)情況下很難滿足,而且氣動載荷也很難在平均軸系下給出[19]。針對這些問題,文獻(xiàn)[20]提出在準(zhǔn)坐標(biāo)系下建立動力學(xué)模型。該模型中準(zhǔn)坐標(biāo)系的原點(diǎn)始終位于未變形飛行器的質(zhì)心,各坐標(biāo)軸的指向也始終保持不變,因此能較好地描述剛彈耦合。但質(zhì)心不與坐標(biāo)原點(diǎn)重合會產(chǎn)生額外的重力矩,并使得該坐標(biāo)系下的平動速度和轉(zhuǎn)動角速度并非真實(shí)值(“準(zhǔn)速度”問題[8])。
針對上述兩種坐標(biāo)系的不足,文獻(xiàn)[8]提出了一種兼具二者特點(diǎn)的新坐標(biāo)系——瞬態(tài)坐標(biāo)系。該坐標(biāo)系既能夠正確描述彈性變形對質(zhì)心位置的影響,又能夠克服準(zhǔn)坐標(biāo)系的“準(zhǔn)速度”問題,因此能更好地體現(xiàn)空氣動力學(xué)、結(jié)構(gòu)動力學(xué)與飛行動力學(xué)之間的相互影響。但文獻(xiàn)[8]以非自旋飛行器為對象,并未考慮旋轉(zhuǎn)導(dǎo)彈的特點(diǎn),也未建立導(dǎo)彈所受力和力矩的計(jì)算模型,難以直接用于彈性旋轉(zhuǎn)導(dǎo)彈的研究。因此,本文在瞬態(tài)坐標(biāo)系下利用拉格朗日方程建立了彈性旋轉(zhuǎn)導(dǎo)彈的動力學(xué)模型。該模型的形式與剛體導(dǎo)彈類似,便于應(yīng)用。最后利用數(shù)值仿真將該模型與其他簡化模型進(jìn)行了對比,發(fā)現(xiàn)該模型能夠更好地描述彈性旋轉(zhuǎn)導(dǎo)彈的動力學(xué)特性。
考慮到實(shí)際情況,為建模方便,做出如下假設(shè):
1) 不考慮地球的曲率和自轉(zhuǎn);
2)不考慮彈體的軸向變形,且橫向變形是小量,可用一系列正交模態(tài)描述;
3)旋轉(zhuǎn)導(dǎo)彈為軸對稱體;
4)旋轉(zhuǎn)導(dǎo)彈在飛行過程中保持小攻角,并采用準(zhǔn)定常氣動力假設(shè)。
采用如圖1所示的3個坐標(biāo)系描述彈性旋轉(zhuǎn)導(dǎo)彈的運(yùn)動,分別定義如下:
1)地面坐標(biāo)系A(chǔ)xyz,記為I. 坐標(biāo)系原點(diǎn)A位于發(fā)射點(diǎn)。Ax軸在水平面內(nèi),指向目標(biāo)為正;Ay軸位于鉛垂平面內(nèi),豎直向上為正;Az軸與Ax軸、Ay軸構(gòu)成右手直角坐標(biāo)系。
2)準(zhǔn)坐標(biāo)系O0x0y0z0,記為B′. 準(zhǔn)坐標(biāo)系的坐標(biāo)原點(diǎn)O0始終位于未變形導(dǎo)彈的質(zhì)心。O0x0軸與未變形彈體縱軸重合,指向頭部為正;O0y0軸位于未變形導(dǎo)彈縱向?qū)ΨQ平面內(nèi)并垂直于O0x0軸,指向上方為正;O0z0軸與O0x0軸、O0y0軸構(gòu)成右手直角坐標(biāo)系。準(zhǔn)坐標(biāo)系用于計(jì)算彈體的彈性變形。
3)瞬態(tài)坐標(biāo)系Obxbybzb,記為B. 坐標(biāo)原點(diǎn)Ob始終位于彈性導(dǎo)彈的質(zhì)心,各坐標(biāo)軸指向始終固定。Obxb軸與未變形彈體的縱軸平行,指向頭部為正;Obyb軸位于未變形導(dǎo)彈縱向?qū)ΨQ平面內(nèi)并垂直于Obxb軸,指向上方為正;Obzb軸與Obxb軸、Obyb軸構(gòu)成右手直角坐標(biāo)系。瞬態(tài)坐標(biāo)系的坐標(biāo)軸指向與準(zhǔn)坐標(biāo)系相同,但原點(diǎn)始終位于飛行器的質(zhì)心,從而使得該坐標(biāo)系可以準(zhǔn)確地描述彈性變形對質(zhì)心位置的影響,并克服“準(zhǔn)速度”問題。
將地面坐標(biāo)系先繞Ay軸轉(zhuǎn)動ψ角,再繞Obzb軸或O0z0軸轉(zhuǎn)動?角,最后繞Obxb軸轉(zhuǎn)動γ角,則兩坐標(biāo)系會重合。由于瞬態(tài)坐標(biāo)系和準(zhǔn)坐標(biāo)系的各坐標(biāo)軸指向是相同的,地面坐標(biāo)系到瞬態(tài)坐標(biāo)系的轉(zhuǎn)換矩陣與地面坐標(biāo)系到準(zhǔn)坐標(biāo)系的轉(zhuǎn)換矩陣相同。地面坐標(biāo)系到瞬態(tài)坐標(biāo)系的轉(zhuǎn)換矩陣L(ψ,?,γ)如(1)式所示:
(1)
彈體的彈性變形可在準(zhǔn)坐標(biāo)系中利用固有振型和模態(tài)坐標(biāo)描述[13]:
(2)
式中:j0和k0分別為準(zhǔn)坐標(biāo)系O0y0軸和O0z0軸的單位向量;Φi(x0)為第i階振型,x0為軸段微元的軸向坐標(biāo),下文為表示方便將Φi(x0)記為Φi;ηy0,i(t)和ηz0,i(t)分別為O0y0軸和O0z0軸方向上的第i階廣義坐標(biāo),與瞬態(tài)坐標(biāo)系Obyb軸和Obzb軸方向上的廣義坐標(biāo)相同,下文為表示方便只記為ηyb,i和ηzb,i;n為振動模態(tài)的總階數(shù)。
因此軸段微元在準(zhǔn)坐標(biāo)系下的位置矢量為
(3)
式中:i0為準(zhǔn)坐標(biāo)系O0x0軸的單位向量。
由(3)式可計(jì)算準(zhǔn)坐標(biāo)系下導(dǎo)彈質(zhì)心的位置矢量ΔC為
(4)
式中:m為導(dǎo)彈質(zhì)量。
(5)
則(4)式可改寫為
(6)
由此可見,λi可以表征彈性變形引起質(zhì)心改變的大小。
綜上所述,由于準(zhǔn)坐標(biāo)系與瞬態(tài)坐標(biāo)系坐標(biāo)軸指向相同,軸段微元在二者下的位置矢量只相差矢量ΔC:
(7)
2.2.1 推力與推力矩
由(2)式可知彈體尾部的變形角為
(8)
FT=Tib-TΔτy0jb-TΔτz0kb,
(9)
式中:T為推力大小。而推力作用點(diǎn)的位置矢量RT(t)為
(10)
式中:xcg為導(dǎo)彈質(zhì)心距頭部的距離。
由(9)式和(10)式可計(jì)算推力矩:
MT=RT×FT=MTxbib+MTybjb+MTzbkb=(TΔτy0zT0-TΔτz0yT0)ib+(TzT0+TΔτz0xT0)jb+(-TyT0-TΔτy0xT0)kb,
(11)
式中:MTxb、MTyb和MTzb分別為推力矩在瞬態(tài)坐標(biāo)系3個坐標(biāo)軸上的分量;xT0、yT0和zT0分別為發(fā)動機(jī)噴口在準(zhǔn)坐標(biāo)系下的三軸坐標(biāo)。
2.2.2 氣動力與氣動力矩
在彈性變形影響下,t時(shí)刻軸向位置xb處軸段微元攻角α(xb,t)和側(cè)滑角β(xb,t)的計(jì)算公式如(12)式所示:
(12)
(13)
u、v和w分別為導(dǎo)彈質(zhì)心速度在瞬態(tài)坐標(biāo)系下的3個分量。
導(dǎo)彈所受氣動力和力矩可由軸段微元所受氣動載荷積分得到,氣動力系數(shù)計(jì)算公式如(14)式所示:
(14)
式中:Cxb為軸向氣動力系數(shù);CNyb和CNzb為側(cè)向氣動力系數(shù);d0(xb)、dA(xb)分別為單位長度的零升阻力系數(shù)和誘導(dǎo)阻力系數(shù);n(xb)為單位長度的側(cè)向氣動力導(dǎo)數(shù)。
將(12)式代入(14)式并忽略2階小量,可得氣動力系數(shù)的具體表達(dá)式為
(15)
(16)
同理可得氣動力矩系數(shù)CMyb和CMzb的表達(dá)式為
(17)
(18)
2.2.3 馬格努斯力矩
與氣動力系數(shù)計(jì)算類似,馬格努斯力矩系數(shù)的計(jì)算公式如(19)式所示:
(19)
(20)
lm為單位長度的馬格努斯力系數(shù)。
2.2.4 廣義力
廣義力可分為沿Obyb軸和Obzb軸兩個方向:
(21)
式中:Qyb,i和Qzb,i為第i階廣義力的分量;q為動壓;S為導(dǎo)彈的特征面積;Fyb,j(t)和Fzb,j(t)為第j個集中載荷的分量;xj為第j個集中載荷作用點(diǎn)距導(dǎo)彈頭部的距離;r為集中載荷的數(shù)量;CQyb,i和CQzb,i為廣義力系數(shù),
(22)
(23)
為簡化推導(dǎo)過程,采用廣義坐標(biāo)系下的拉格朗日方程[19-20]建立動力學(xué)方程如下:
(24)
(25)
(26)
D為阻尼耗散能。
依次計(jì)算動能、勢能和阻尼耗散能并代入(24)式,即可得到彈性旋轉(zhuǎn)導(dǎo)彈的動力學(xué)方程。
3.2.1 導(dǎo)彈動能
如圖2所示,軸段微元在瞬態(tài)坐標(biāo)系下的位置矢量為
R=R0+p.
(27)
對(27)式等號兩側(cè)進(jìn)行微分,可得
(28)
由(28)式即可計(jì)算導(dǎo)彈的動能K:
(29)
式中:J0為旋轉(zhuǎn)導(dǎo)彈無彈性變形時(shí)的慣量張量,
(30)
Ixx0、Itt0分別為未變形導(dǎo)彈縱向和橫向的轉(zhuǎn)動慣量;d為軸段微元在瞬態(tài)坐標(biāo)系下位置矢量的側(cè)向分量;ρ為軸段微元在瞬態(tài)坐標(biāo)系下位置矢量的軸向分量。
將各分量代入(29)式,可得導(dǎo)彈動能表達(dá)式:
(31)
式中:ωxb為導(dǎo)彈自旋角速度;σi為與振型有關(guān)的常數(shù),
(32)
Mi為第i階模態(tài)質(zhì)量,
(33)
(31)式等號右側(cè)的第1項(xiàng)是旋轉(zhuǎn)導(dǎo)彈的平動動能;第2項(xiàng)~第8項(xiàng)是考慮了質(zhì)量特性變化的導(dǎo)彈轉(zhuǎn)動動能;第9項(xiàng)~第12項(xiàng)是彈性變形和導(dǎo)彈轉(zhuǎn)動耦合產(chǎn)生的動能;第13項(xiàng)~第16項(xiàng)是彈性變形引起的動能。
3.2.2 導(dǎo)彈勢能與阻尼耗散能
導(dǎo)彈勢能U包含重力勢能Ug和彈性勢能Ue,
U=Ug+Ue,
(34)
式中:
(35)
Rx、Ry和Rz為質(zhì)心位置矢量在慣性坐標(biāo)系下的3個分量,g和g分別為引力加速度矢量和引力加速度大?。?/p>
(36)
ωi為第i階模態(tài)的固有角頻率。
導(dǎo)彈阻尼耗散能的表達(dá)式為
(37)
式中:μi為第i階模態(tài)的臨界阻尼系數(shù)。
3.2.3 動力學(xué)方程
將(31)式、(34)式和(37)式代入(24)式,即可得到動力學(xué)方程。根據(jù)(24)式,動力學(xué)方程可以分為3組:質(zhì)心平動動力學(xué)方程、繞質(zhì)心轉(zhuǎn)動動力學(xué)方程和彈性變形動力學(xué)方程。其中質(zhì)心平動動力學(xué)方程如(38)式所示:
(38)
繞質(zhì)心轉(zhuǎn)動動力學(xué)方程如(39)式~(41)式所示:
(39)
(40)
(41)
彈性變形動力學(xué)方程如(42)式和(43)式所示:
(42)
(43)
由(38)式可知,由于瞬態(tài)坐標(biāo)系的原點(diǎn)始終位于彈性旋轉(zhuǎn)導(dǎo)彈瞬時(shí)質(zhì)心,建立在該坐標(biāo)系下的質(zhì)心平動動力學(xué)方程與剛體導(dǎo)彈的質(zhì)心平動動力學(xué)方程形式一致,從而給后續(xù)研究帶來了方便。
由(39)式~(43)式可知,在繞質(zhì)心轉(zhuǎn)動動力學(xué)方程與彈性變形動力學(xué)方程中,剛體運(yùn)動自由度與彈性變形自由度相互耦合在一起。而文獻(xiàn)[14]建立的動力學(xué)方程忽略了所有耦合項(xiàng),實(shí)現(xiàn)了繞質(zhì)心轉(zhuǎn)動動力學(xué)方程與彈性變形動力學(xué)方程間的解耦:繞質(zhì)心轉(zhuǎn)動動力學(xué)方程只包含剛體運(yùn)動自由度,彈性變形動力學(xué)方程中只包含彈性變形自由度。這表明建立在瞬態(tài)坐標(biāo)系下的動力學(xué)方程考慮了彈性變形對質(zhì)量特性的影響,能夠完整地描述剛彈耦合特性。這是建立在平均軸系和準(zhǔn)坐標(biāo)系下的動力學(xué)模型無法做到的,因此稱本文建立的模型為慣性耦合完整模型。
事實(shí)上,基于不同假設(shè)對該模型進(jìn)行簡化,即可得到其他文獻(xiàn)使用的模型。
由(5)式可知,λi表征彈性變形對質(zhì)心位置的影響,σi表征彈性變形與剛體旋轉(zhuǎn)耦合引起的動能增量。平均軸系中的建模忽略了這兩種耦合,因此有
(44)
將(44)式代入動力學(xué)方程(40)式~(43)式,則繞質(zhì)心運(yùn)動動力學(xué)方程簡化為
(45)
(46)
而結(jié)構(gòu)變形動力學(xué)方程簡化為
(47)
(48)
(45)式~(48)式組成的簡化動力學(xué)模型與文獻(xiàn)[6-7]中的動力學(xué)模型一致,只保留了一部分耦合項(xiàng),本文稱為慣性耦合簡化模型。
進(jìn)一步地,忽略彈性變形對質(zhì)量特性的影響,可得到與剛體導(dǎo)彈動力學(xué)模型形式一致的彈性旋轉(zhuǎn)導(dǎo)彈動力學(xué)模型,如(49)式~(52)式所示:
(49)
(50)
(51)
(52)
該模型忽略了所有耦合項(xiàng),因此彈性變形自由度與剛體運(yùn)動自由度在形式上是解耦的,二者之間的耦合僅體現(xiàn)在力和力矩的計(jì)算上。本文稱該模型為非慣性耦合模型。
綜上所述,建立在瞬態(tài)坐標(biāo)系下的慣性耦合完整模型考慮了彈性變形對質(zhì)量特性的影響,能夠全面描述剛彈耦合效應(yīng)。而且慣性耦合完整模型利用微分方程描述彈性旋轉(zhuǎn)導(dǎo)彈的運(yùn)動,只增加了有限數(shù)量的彈性變形自由度,便于工程應(yīng)用。
某旋轉(zhuǎn)導(dǎo)彈的前2階固有頻率分別為43 Hz和106 Hz,相應(yīng)的振型如圖3所示,其他參數(shù)如表1所示。
表1 旋轉(zhuǎn)導(dǎo)彈參數(shù)
下面以該旋轉(zhuǎn)導(dǎo)彈為對象進(jìn)行仿真,分析3.3節(jié)中3種模型的差別。
分別使用慣性耦合完整模型和慣性耦合簡化模型進(jìn)行仿真。自旋角速度為263.5 rad/s的仿真結(jié)果如圖4~圖8所示。由圖4~圖8可以看出,此時(shí)完整模型的仿真結(jié)果已經(jīng)開始發(fā)散,而簡化模型的結(jié)果尚未發(fā)散。此外,簡化模型和完整模型關(guān)于剛體攻角α0、剛體側(cè)滑角β0和角速度ωyb、ωzb的仿真結(jié)果存在著較明顯的差別,其中角速度的相位一致,差別主要體現(xiàn)在幅值上,而模態(tài)坐標(biāo)ηyb,1、ηzb,1則相差不大。
自旋角速度為160 rad/s的仿真結(jié)果如圖9~圖13所示。由圖9~圖13可以看出,此時(shí)完整模型和簡化模型仿真結(jié)果之間的差別明顯減小,只有剛體攻角α0、剛體側(cè)滑角β0和角速度ωyb、ωzb間存在一定的差別,模態(tài)坐標(biāo)ηyb,1、ηzb,1的仿真結(jié)果基本一致。
自旋角速度為60 rad/s的仿真結(jié)果如圖14~圖18所示。此時(shí)兩種模型的仿真結(jié)果基本一致。
以上3組仿真結(jié)果表明,在自旋轉(zhuǎn)速較大時(shí),簡化模型和完整模型關(guān)于剛體攻角、剛體側(cè)滑角和角速度的仿真結(jié)果存在較明顯的差別,其中角速度的相位一致,差別主要體現(xiàn)在幅值上,而兩種模型關(guān)于模態(tài)坐標(biāo)的結(jié)果則基本一致。表明剛彈耦合效應(yīng)對姿態(tài)運(yùn)動有一定影響,且這種影響主要體現(xiàn)在角速度的幅值上;在導(dǎo)彈運(yùn)動不發(fā)散情況下,兩種模型在彈性變形方面的差別很小,可以忽略。此外,隨著自旋角速度的增加,彈性變形自由度的頻率減小,與(42)式和(43)式的結(jié)果一致。
綜上所述可知,慣性耦合完整模型和慣性耦合簡化模型之間的差別相對較小,在彈性變形頻率較大時(shí)甚至可以忽略,而在彈性變形頻率較小時(shí)二者之間差別較大,需要加以考慮。因此在旋轉(zhuǎn)導(dǎo)彈自旋轉(zhuǎn)速較小(遠(yuǎn)小于彈性模態(tài)固有頻率)時(shí),可使用較簡單的慣性耦合簡化模型。
下面分別使用慣性耦合完整模型和無慣性耦合模型進(jìn)行仿真,結(jié)果如圖19~圖23所示。
由圖19、圖20和圖21可知,兩種模型關(guān)于剛體攻角α0、剛體側(cè)滑角β0和角速度ωyb、ωzb的仿真結(jié)果中幅值與相位均存在明顯差別。由圖22和圖23可知關(guān)于彈性變形的仿真結(jié)果基本一致。表明慣性耦合對旋轉(zhuǎn)導(dǎo)彈剛體姿態(tài)運(yùn)動的影響較大,對彈性變形的影響則很小。因此對于彈性較大的旋轉(zhuǎn)導(dǎo)彈,目前使用較多的平均軸系下的動力學(xué)方程不夠準(zhǔn)確。
為準(zhǔn)確描述彈性變形對質(zhì)心位置的影響并克服準(zhǔn)坐標(biāo)系下的“準(zhǔn)速度”問題,本文引入了瞬態(tài)坐標(biāo)系,并在該坐標(biāo)系下建立了能夠全面描述剛彈耦合效應(yīng)的彈性旋轉(zhuǎn)導(dǎo)彈動力學(xué)模型。所得主要結(jié)論如下:
1) 慣性耦合完整模型在不同假設(shè)條件下可簡化為目前使用較多的慣性耦合簡化模型和無慣性耦合模型:慣性耦合完整模型忽略彈性變形導(dǎo)致的質(zhì)心改變和剛彈耦合產(chǎn)生的動能增量可得到慣性耦合簡化模型,進(jìn)一步忽略彈性變形對質(zhì)量特性的影響可得到無慣性耦合模型。
2)慣性耦合簡化模型在導(dǎo)彈自旋轉(zhuǎn)速較低時(shí)與慣性耦合完整模型差別較小,在自旋轉(zhuǎn)速較高時(shí)的差別則不可忽略;無慣性耦合模型所忽略掉的慣性耦合項(xiàng)對角速度幅值和相位均影響較大,無法準(zhǔn)確描述彈性旋轉(zhuǎn)導(dǎo)彈的姿態(tài)運(yùn)動。
3) 慣性耦合完整模型利用微分方程描述彈性旋轉(zhuǎn)導(dǎo)彈的運(yùn)動,在只增加有限數(shù)量自由度的情況下能夠完整地描述剛彈耦合效應(yīng),便于工程應(yīng)用。
4)在旋轉(zhuǎn)導(dǎo)彈自旋轉(zhuǎn)速較小(遠(yuǎn)小于彈性模態(tài)固有頻率)時(shí),可使用較簡單的慣性耦合簡化模型。但該模型的應(yīng)用仍需經(jīng)過模型的校驗(yàn)、驗(yàn)證和確認(rèn)。