彈箭角運(yùn)動(dòng)的非線性吸引域估計(jì)

李東陽, 常思江, 王中原

(南京理工大學(xué) 能源與動(dòng)力工程學(xué)院, 江蘇 南京 210094)

0 引言

對于飛行中的彈箭,在攻角δ較小的情況下,可將作用在彈箭上的氣動(dòng)力和力矩視為攻角的線性函數(shù),并可采用sinδ≈δ、cosδ≈1的線性化假設(shè),據(jù)此形成線性化外彈道理論,該理論曾在很大范圍內(nèi)成功地預(yù)示了彈箭運(yùn)動(dòng),促進(jìn)了彈箭飛行性能設(shè)計(jì)的發(fā)展。然而,彈箭在各種條件下的實(shí)際運(yùn)動(dòng),總是具有不同程度的非線性[1]。隨著彈箭外形及其飛行條件的日益復(fù)雜化、多樣化,彈箭運(yùn)動(dòng)的非線性特征也愈發(fā)凸顯。如在大攻角情況下,作用于彈箭的氣動(dòng)力和力矩一般為攻角δ的非線性函數(shù),有的還是彈箭滾轉(zhuǎn)方位角的非線性函數(shù)(如誘導(dǎo)側(cè)向力矩),因而數(shù)學(xué)上無法利用線性理論開展研究;對于一些客觀存在的現(xiàn)象,如“行進(jìn)軍艦上發(fā)射某尾翼式低速旋轉(zhuǎn)火箭,左舷發(fā)射穩(wěn)定、右舷發(fā)射不穩(wěn)定”、“某些火箭彈飛行中出現(xiàn)錐形運(yùn)動(dòng)而致失穩(wěn)”等,線性理論也無法給出合理的解釋。因此,彈箭非線性運(yùn)動(dòng)理論應(yīng)運(yùn)而生并迅速發(fā)展。

20世紀(jì)60年代至本世紀(jì)初,一些學(xué)者就強(qiáng)非線性靜力矩[2]、非對稱再入飛行器的非線性運(yùn)動(dòng)[3-4]、轉(zhuǎn)速-攻角閉鎖與災(zāi)變性偏航[5-6]、彈箭非線性動(dòng)力學(xué)[7-8]、有控彈非線性穩(wěn)定性[9]等問題,利用非線性數(shù)學(xué)工具開展研究,建立了彈箭非線性運(yùn)動(dòng)研究的理論方法體系,取得了大量研究成果;Mc Coy等[10]通過開展靶道飛行試驗(yàn)及氣動(dòng)系數(shù)辨識(shí),獲得了大量的彈丸非線性氣動(dòng)力數(shù)據(jù)。近20年來,隨著彈箭技術(shù)的進(jìn)一步發(fā)展,彈箭非線性運(yùn)動(dòng)研究繼續(xù)受到廣泛關(guān)注。Morote 等就無控火箭共振閉鎖[11]、尾翼彈災(zāi)變性偏航控制[12-13]、尾翼彈非線性滾轉(zhuǎn)運(yùn)動(dòng)[14]、高階非線性氣動(dòng)力[15]等開展了一系列深入研究;韓子鵬等[16]對彈箭非線性運(yùn)動(dòng)的國內(nèi)外研究成果進(jìn)行了系統(tǒng)梳理,主要涉及非旋轉(zhuǎn)尾翼彈平面非線性運(yùn)動(dòng)、旋轉(zhuǎn)彈箭非線性運(yùn)動(dòng)、非線性氣動(dòng)系數(shù)獲取等方面;李臣明等[17]對遠(yuǎn)程彈箭的轉(zhuǎn)速閉鎖現(xiàn)象進(jìn)行了理論;Li等[18]、Chang等[19]和李東陽等[20]重點(diǎn)研究了彈箭非線性運(yùn)動(dòng)方程的解析解,其中文獻(xiàn)[18]在Murphy 等[21]關(guān)于轉(zhuǎn)速-攻角閉鎖解析模型的基礎(chǔ)上,進(jìn)一步完善了攻角-轉(zhuǎn)速閉鎖的解析解,文獻(xiàn)[19-20]則分別探討了利用同倫分析法、正規(guī)形法求取高精度攻角解析解的可行性,并據(jù)此開展了彈箭穩(wěn)定性分析、參數(shù)設(shè)計(jì)等相關(guān)應(yīng)用;此外,Cross等[22]利用數(shù)值仿真對制導(dǎo)炮彈的非線性穩(wěn)定性開展了分析;鐘揚(yáng)威等[23]、楊杰等[24]利用分岔理論對彈箭非線性角運(yùn)動(dòng)等開展了研究。可以預(yù)見,彈箭非線性運(yùn)動(dòng)理論將是從根本上提高彈箭性能設(shè)計(jì)水平的重要基礎(chǔ)和有力工具。

彈箭非線性運(yùn)動(dòng)的穩(wěn)定性不僅與系統(tǒng)參數(shù)(如氣動(dòng)參數(shù)、結(jié)構(gòu)參數(shù)等)有關(guān),而且與初始條件密切相關(guān),這是非線性系統(tǒng)與線性系統(tǒng)的不同之處。故要使彈箭達(dá)到期望的穩(wěn)定狀態(tài),必須對其初始條件加以限定。定義一個(gè)包含初始條件的緊集,若從該集合內(nèi)任一點(diǎn)出發(fā)的軌跡不超出該集合且最終都將收斂到系統(tǒng)的某一個(gè)平衡點(diǎn),則稱該緊集為對應(yīng)穩(wěn)定平衡點(diǎn)的吸引域[25-26]。因此,只有當(dāng)初始運(yùn)動(dòng)狀態(tài)在吸引域內(nèi)時(shí),彈箭運(yùn)動(dòng)才是穩(wěn)定的,這就使得吸引域的確定變得十分重要。在開展彈箭穩(wěn)定性設(shè)計(jì)時(shí),必須通過設(shè)計(jì)相關(guān)參數(shù),確保吸引域包含彈箭的可能運(yùn)動(dòng)狀態(tài)范圍。然而,一般來說,除一些特殊的系統(tǒng)外,非線性系統(tǒng)準(zhǔn)確的吸引域確定相當(dāng)困難[27],一般通過解析或數(shù)值計(jì)算方法進(jìn)行估計(jì)?，F(xiàn)有的吸引域估計(jì)方法,主要采用李雅普諾夫理論進(jìn)行問題的描述[28]。對于多項(xiàng)式系統(tǒng),在吸引域問題的求解上普遍采用多項(xiàng)式的平方和方法[25]。該方法的本質(zhì)是將多項(xiàng)式的正定要求松弛為平方和問題,進(jìn)而利用平方和規(guī)劃方法求解[29]。平方和規(guī)劃包含了多項(xiàng)式半正定問題的求解、平方和優(yōu)化問題與平方和可行性問題的求解,以上問題在具體的吸引域估計(jì)問題中都會(huì)遇到[26]。盡管該方法具有一定的保守性,且其計(jì)算復(fù)雜度隨著系統(tǒng)的維數(shù)和所尋求的李雅普諾夫函數(shù)階次的增大而迅速增大,以及僅適用于多項(xiàng)式系統(tǒng)的局限性,但它使得吸引域的估計(jì)變得可行,與基于蒙特卡洛軌跡積分的方法相比,平方和規(guī)劃方法給出的結(jié)果具有可靠性和絕對性[30],并給出了吸引域的解析表達(dá)式,這是該方法的優(yōu)勢所在。

因此,本文擬探討平方和方法在彈箭角運(yùn)動(dòng)非線性吸引域估計(jì)中的應(yīng)用。對此,將針對兩種情形開展吸引域估計(jì):一是無控尾翼彈平面非線性運(yùn)動(dòng),考慮高階非線性靜力矩以及二次方非線性阻尼;二是脈沖末段修正迫彈[31]的空間非線性運(yùn)動(dòng),考慮非線性法向力系數(shù)、馬格努斯力矩系數(shù)、俯仰力矩系數(shù)及赤道阻尼力矩系數(shù)等。由于平方和規(guī)劃僅適用于多項(xiàng)式系統(tǒng),故在角運(yùn)動(dòng)方程的建模中需要對幾何非線性等非多項(xiàng)式因素進(jìn)行多項(xiàng)式近似。本文研究旨在為彈箭非線性運(yùn)動(dòng)研究提供新的有效理論工具。

1 基于平方和規(guī)劃的吸引域估計(jì)方法

本節(jié)將簡要介紹基于平方和規(guī)劃的非線性吸引域估計(jì)方法。對于一個(gè)非線性自治多項(xiàng)式:

(1)

式中:x為狀態(tài)變量,x∈n;t為自變量;f(x)是n×1維的多項(xiàng)式向量場;x0為t=0時(shí)的狀態(tài)值,即初始條件。不失一般性,設(shè)原點(diǎn)是系統(tǒng)的平衡點(diǎn),于是f(0)=0,且原點(diǎn)漸近穩(wěn)定,則根據(jù)直接李雅普諾夫定理(定理1),原點(diǎn)的吸引域可定義為

{x0∈n: 若x(0)=x0則

(2)

定理1[25-26]若存在連續(xù)可微標(biāo)量函數(shù)V(x):n→和標(biāo)量γ∈+,使得

V(x)>0, ?x≠0∧V(0)=0

(3)

Ω:={x:V(x)≤γ}

(4)

Ω?{x:(?V(x)/?x)f(x)<0}∪{0}

(5)

則原點(diǎn)漸近穩(wěn)定,Ω是原點(diǎn)的確切吸引域式(2)的子集,可作為吸引域的一個(gè)估計(jì),γ為其水平集。

可見,所尋找吸引域的估計(jì),需滿足式(3)～式(5)對V(x)和γ的約束。而式(3)～式(5)中,既含有不等式約束,又含有集合包含約束。利用代數(shù)幾何學(xué),集合包含關(guān)系可以通過著名的推廣S-procedure(定理 2)[32]轉(zhuǎn)化為不等式約束。

定理2(推廣S-procedure)給出多項(xiàng)式g0(x),…,gm(x)∈R[x]和多項(xiàng)式s1(x),…,sm(x)∈Σn,R[x]代表系數(shù)為實(shí)數(shù)的多項(xiàng)式集合,Σn代表n個(gè)變量組成的平方和多項(xiàng)式集合,若

(6)

則

{x|g1(x),…,gm(x)≥0}?{x|g0(x)≥0}

(7)

因此,問題轉(zhuǎn)化為處理一系列多項(xiàng)式不等式。通常一個(gè)多元多項(xiàng)式的非負(fù)性是很難確定的,屬于非確定性多項(xiàng)式(NP)困難問題,解決的思路是利用平方和松弛方法[25],將多項(xiàng)式的非負(fù)性問題轉(zhuǎn)化為平方和問題,進(jìn)而利用平方和規(guī)劃方法[25-26]對問題進(jìn)行求解。例如,驗(yàn)證一個(gè)多項(xiàng)式s(x)∈R(x)是否為正,就等價(jià)于驗(yàn)證是否存在正定矩陣Q,使得

s(x)=ZT(x)QZ(x)∈Σn

(8)

式中:Z(x)為單項(xiàng)式向量。關(guān)于平方和方法的具體內(nèi)容,可參見文獻(xiàn)[26,32]。

綜上所述,吸引域估計(jì)問題可以轉(zhuǎn)化為以下平方和問題:

(9)

式中:l1(x)∈R(x)且l1(x)>0,一般取εxTx,ε為大于0的常數(shù)(如取為10-6);s1(x)輔助算子,為一定階數(shù)的平方和多項(xiàng)式。

為了進(jìn)一步優(yōu)化吸引域的估計(jì),可通過一個(gè)稱為形狀函數(shù)的多項(xiàng)式h(x)>0[26, 33],將原問題式(9)改寫為一個(gè)雙線性平方和規(guī)劃問題:

(10)

式中:φ為形狀函數(shù)h(x)的水平集,φ∈;l2與l1類似,同為多項(xiàng)式小量;s2與s1類似,亦為平方和輔助算子。

將上述雙線性問題解耦為兩個(gè)單線性優(yōu)化問題[26],采用迭代算法(常稱為V-s迭代算法[26-27, 33])進(jìn)行求解,給定一個(gè)已知且可行的初始李雅普諾夫函數(shù)V0和形狀函數(shù)h,并令V(x)=V0,具體步驟簡述如下:

1)固定V(x),利用二分法求得γ的最大值γ*,并記下此時(shí)的s2,即

2)固定V(x)、變換γ*,利用二分法求得φ的最大值φ*,并記下此時(shí)的s1,即

3)固定s1、s2、φ*、γ*,求解新的滿足以下約束的李雅普諾夫函數(shù)V(x),這里可根據(jù)需要通過待定系數(shù)設(shè)定V(x)的階次(如2次、4次或6次李雅普諾夫函數(shù)):

-[(?V(x)/?x)f+l2]-(γ*-V(x))s2∈Σn
(γ*-V(x))-(φ*-h)s1∈Σn
V(x)-l1∈Σn
V(0)=0

4)V(x)=V(x)/γ*;

5)重復(fù)以上步驟,得到最終的李雅普諾夫函數(shù)V(x)=V(x)/γ*和相應(yīng)的吸引域估計(jì)Ω:={x∈n:V(x)<1}。停止條件可為以下任意一種:

①遇到數(shù)值計(jì)算不可行問題;

②相鄰兩次φ的最優(yōu)值小于既定容許度,如0.001;

③預(yù)定的迭代次數(shù)i。

值得注意的是,在尋找新的李雅普諾夫函數(shù)這一步(步驟3)中,s1、s2為已知且固定,與李雅普諾夫函數(shù)的階次無直接關(guān)系。

2 考慮高階氣動(dòng)力的無控彈箭角運(yùn)動(dòng)吸引域

引起彈箭非線性運(yùn)動(dòng)的因素包括氣動(dòng)力非線性、結(jié)構(gòu)非線性及幾何非線性等,其中又以氣動(dòng)力非線性占主導(dǎo)地位[1]。Morote等[12]的研究結(jié)果表明,對于某彈箭,當(dāng)靜力矩系數(shù)關(guān)于攻角的階數(shù)達(dá)到7次時(shí)才能準(zhǔn)確地?cái)M合試驗(yàn)數(shù)據(jù);Liano等[15]也發(fā)現(xiàn),為準(zhǔn)確預(yù)測轉(zhuǎn)速閉鎖,當(dāng)攻角為12°和20°時(shí),與滾轉(zhuǎn)角相關(guān)的非線性力矩關(guān)于攻角的階數(shù)分別不能低于 5次和7次。因此,對于彈箭非線性運(yùn)動(dòng)研究,有必要引入高階氣動(dòng)力系數(shù)。將利用第1節(jié)中的方法,研究無控尾翼彈平面非線性角運(yùn)動(dòng)的吸引域估計(jì)。

2.1 無控尾翼彈非線性平面角運(yùn)動(dòng)方程

引入高階氣動(dòng)力,考慮5次方和7次方非線性靜力矩以及2次非線性阻尼對角運(yùn)動(dòng)的影響。無控尾翼彈的平面非線性角運(yùn)動(dòng)方程可描述為

δ″+(H0+H2δ2)δ′-(M0+M2δ2+M4δ4+M6δ6)δ=0

(11)

2.2 平衡點(diǎn)的穩(wěn)定性

若該角運(yùn)動(dòng)方程的平衡點(diǎn)記為δ*,其滿足:

(12)

因此,在高階非線性靜力矩系數(shù)的影響下,角運(yùn)動(dòng)可能存在除原點(diǎn)以外的其他平衡點(diǎn)。當(dāng)不考慮阻尼項(xiàng)時(shí),平衡點(diǎn)對應(yīng)的特征根λ滿足

λ2=M0+3M2δ*2+5M4δ*4+7M6δ*6

(13)

故原點(diǎn)對應(yīng)的特征根滿足λ2=M0<0,則平衡點(diǎn)為中心,穩(wěn)定性無法由線性理論確定。以3次非線性靜力矩為例,當(dāng)M0、M2同號時(shí),不存在其他平衡點(diǎn),理論上不論初始條件如何,攻角都將做幅值為δ0的擺動(dòng);當(dāng)M0、M2異號時(shí),存在一對實(shí)平衡點(diǎn)δ*2=-M0/M2,且當(dāng)M0<0時(shí),δ*為不穩(wěn)定鞍點(diǎn)。

當(dāng)考慮阻尼項(xiàng)時(shí),阻尼項(xiàng)的存在雖然不改變角運(yùn)動(dòng)方程式(11)的平衡點(diǎn)位置,但影響其穩(wěn)定性,此時(shí)特征根滿足

λ2+(H0+H2δ*2)λ-(M0+3M2δ*2+5M4δ*4+7M6δ*6)=0

(14)

為了確定攻角運(yùn)動(dòng)的穩(wěn)定性,需綜合考慮系統(tǒng)平衡點(diǎn)式(12)和阻尼項(xiàng)H0、H2的影響。原點(diǎn)的穩(wěn)定性將由線性靜力矩系數(shù)M0和線性阻尼系數(shù)H0共同確定,即當(dāng)M0<0時(shí),若H0<0,則原點(diǎn)為不穩(wěn)定平衡點(diǎn),若H0>0,則原點(diǎn)為穩(wěn)定平衡點(diǎn)。此外,文獻(xiàn)[20]表明,考慮阻尼時(shí),符號相反的H0、H2使得非線性的角運(yùn)動(dòng)系統(tǒng)出現(xiàn)了極限環(huán)。極限環(huán)的穩(wěn)定性與原點(diǎn)穩(wěn)定性的共同作用決定了攻角幅值的變化規(guī)律,如表1所示。詳細(xì)討論參見文獻(xiàn)[20]。

表1 系統(tǒng)穩(wěn)定性分析

由表1可知,對于不同符號H0、H2的組合,可能存在穩(wěn)定或不穩(wěn)定的極限環(huán)。對于原點(diǎn)為穩(wěn)定平衡點(diǎn)的系統(tǒng),設(shè)計(jì)者感興趣的是彈箭在什么初始條件下角運(yùn)動(dòng)可以收斂到原點(diǎn)。

2.3 吸引域估計(jì)結(jié)果

采用文獻(xiàn)[19]中的平面角運(yùn)動(dòng)參數(shù)開展具體的吸引域估計(jì),其中馬赫數(shù)Ma=0.6,靜力矩系數(shù)M0=-5.0×10-5,M2=-4.5×10-4,M4=-8.0×10-3,M6=0(由于缺乏數(shù)據(jù),這里取為0值),阻尼力矩系數(shù)H0=2.0×10-3,H2=-0.4。如表1所示,對于該H0、H2的具體值,原點(diǎn)為穩(wěn)定平衡點(diǎn),且存在一個(gè)不穩(wěn)定極限環(huán),相軌線圖如圖1所示。

圖1 不同初始條件下攻角運(yùn)動(dòng)相平面圖

當(dāng)角運(yùn)動(dòng)初始條件位于此極限環(huán)內(nèi)時(shí),攻角運(yùn)動(dòng)收斂到終點(diǎn);角運(yùn)動(dòng)初始條件在此極限環(huán)外時(shí),攻角運(yùn)動(dòng)發(fā)散。因此,穩(wěn)定原點(diǎn)的確切吸引域由該極限環(huán)表示(見圖1中的陰影區(qū)域),該極限環(huán)可通過數(shù)值計(jì)算并逆向積分得到。需要說明的是,并非任何系統(tǒng)都可以采用逆向數(shù)值積分方法獲得其確切吸引域。如Van de Pol系統(tǒng),由于其形成了不穩(wěn)定極限環(huán),從環(huán)內(nèi)所有點(diǎn)出發(fā)的軌跡均收斂于環(huán)內(nèi)平衡點(diǎn),故該平衡點(diǎn)的吸引域就是該極限環(huán)包圍的區(qū)域,進(jìn)而可以通過逆向積分,使軌跡最終收斂于極限環(huán),得到確切吸引域,而角運(yùn)動(dòng)方程式(11)所示系統(tǒng),本質(zhì)上恰為一Van de Pol系統(tǒng)。

取x1=δ、x2=δ′,則角運(yùn)動(dòng)方程式(11)可表示為

(15)

利用V-s迭代算法對原點(diǎn)的吸引域進(jìn)行估計(jì),選取初始李雅普諾夫函數(shù)

V0=xTPx

(16)

式(16)中的P由李雅普諾夫方程ATP+PA=-I計(jì)算得到,A=(?f/?x)|x=0,I為單位矩陣;形狀函數(shù)取為h(x)=xTx;利用V-s迭代算法,分別尋找 2次、4次和6次李雅普諾夫函數(shù)對吸引域進(jìn)行估計(jì),最終得到的李雅普諾夫函數(shù)分別記為V2、V4、V6:

V2=0.031 708δ2+134.146 533δδ′+585 313 523.361 843δ′2

(17)

V4=1.341 471e-5δ4-0.090 763δ3δ′+453 519.717 597δ2δ′2+228 473 313 060.431 2δδ′3+2.277 762e16δ′4-3.917 141e-7δ3+0.085 829δ2δ′+776.309 701δδ′2+2 696 018 997.682 03δ′3+0.027 334δ2+481.738 113δδ′+481 989 806.810 080δ′2

(18)

V6=6.535 639e9δ6-0.002 437δ5δ′-20.835 500δ4δ′2+5.912 217e8δ3δ′3+4.210 891e14δ2δ′4+0.000 19δ4δ′+1.829 569e19δδ′5+2.396 202e24δ′6-3.296 192e-10δ5-32.114 797δ3δ′2+9.010 24e6δ2δ′3-4.976 216e11δδ′4-1.172 415 2e16δ′5+1.818 161e-5δ4-6.580 880δ3δ′-7.076 783e4δ2δ′2+1.654 866e11δδ′3-8.684 899e14δ′4+2.375 432e-8δ3-0.010 868δ2δ′+1 825.386 293δδ′2+37 720 324.123 641 36δ′3+0.014 637 006 762 577 65δ2+429.817 180 706 390 3δδ′+285 374 579.309 986 1δ′2

(19)

吸引域估計(jì)結(jié)果如圖2所示。

圖2 不同階次李雅普諾夫函數(shù)估計(jì)下的吸引域

由圖2可見,隨著李雅普諾夫函數(shù)階次的增加,吸引域的近似效果逐漸提升。6階李雅普諾夫函數(shù)的結(jié)果已無限接近確切吸引域。這一結(jié)果正如文獻(xiàn)[34]所得結(jié)論,李雅普諾夫函數(shù)的階次越高,可以越準(zhǔn)確地接近確切吸引域。主要原因在于,當(dāng)李雅普諾夫函數(shù)的階次越高,可表達(dá)成的形狀也越為多樣,故高階李雅普諾夫函數(shù)有能力表示形狀復(fù)雜的吸引域。因此,基于平方和規(guī)劃的非線性吸引域估計(jì)方法可用于分析彈箭角運(yùn)動(dòng)的穩(wěn)定性,能夠定量地確定角運(yùn)動(dòng)的初始穩(wěn)定范圍。

3 脈沖末修迫彈的非線性角運(yùn)動(dòng)吸引域

隨著彈箭技術(shù)的發(fā)展,彈箭非線性運(yùn)動(dòng)研究已不局限于無控彈箭,對于有控彈箭的非線性分析也已開展過相關(guān)研究[9, 22],但限于研究工具的缺乏,這些工作尚不深入。將以脈沖控制的末修迫彈[31]為對象,開展非線性吸引域估計(jì),根據(jù)所得吸引域反求脈沖參數(shù)設(shè)計(jì)應(yīng)滿足的約束,據(jù)此提出脈沖控制參數(shù)的設(shè)計(jì)方法。

3.1 非線性角運(yùn)動(dòng)方程及其多項(xiàng)式近似

為便于描述,定義彈軸坐標(biāo)系(簡稱彈軸系),其原點(diǎn)O位于彈質(zhì)心,x軸沿彈軸指向彈體頭部為正,y軸沿水平方向,z軸在鉛錘面內(nèi)垂直于y軸指向下為正。在彈軸系內(nèi),彈體的速度和角速度分別表示為U=[uvw]T和ω=[pqr]T,u、v、w分別為彈體速度在彈軸坐標(biāo)系三軸上的分量,p、q、r分別為彈體角速度在彈軸坐標(biāo)系三軸上的分量;彈軸系相對于慣性坐標(biāo)系的角速度在彈軸系內(nèi)的投影可表示為Ωx=[pxqr]T,px表示該角速度的軸向分量。作用在彈體的橫向(即沿y軸和z軸)氣動(dòng)力和氣動(dòng)力矩可以表示為如下復(fù)數(shù)形式:

(20)

彈體的橫向運(yùn)動(dòng)在彈軸系下可表示為

(21)

(22)

不考慮馬格努斯力,氣動(dòng)力和氣動(dòng)力矩系數(shù)可表示為

Cy+iCz=-CNαξ

(23)

(24)

其中:CNα為法向力系數(shù);CMpα為馬格努斯力矩系數(shù);CMα為俯仰力矩系數(shù)。

考慮氣動(dòng)力非線性,則有

(25)

式中:CNα0、CNα2分別為線性和三次方法向力系數(shù);CMpα0、CMpα2分別為線性和三次方馬格努斯力矩系數(shù)。

將式(23)和式(24)代入橫向運(yùn)動(dòng)方程組式(21),且彈軸系下角速度px≈0,并將自變量從時(shí)間t換為無量綱弧長s。為方便起見,選擇變量

(26)

則角運(yùn)動(dòng)方程可寫為

(27)

(28)

式中:η為幾何非線性項(xiàng),η=u/U。

由幾何關(guān)系可知

(29)

則有

(30)

由氣動(dòng)力和重力沿速度方向的分量可得速度方程

(31)

其中:CD為阻力系數(shù)。

由于平方和計(jì)算工具(如SOSTOOLs、SOSOPTs、SeDuMi等[26-27, 33])只能應(yīng)用于多項(xiàng)式,故式(27)需表示為多項(xiàng)式的形式,即對幾何非線性η進(jìn)行多項(xiàng)式近似。

為此,將式(29)和式(30)泰勒展開并取1次近似,可得

(32)

式中:O(·)為高階小量表示符號。

為驗(yàn)證上述多項(xiàng)式近似系統(tǒng)代替原系統(tǒng)進(jìn)行角運(yùn)動(dòng)穩(wěn)定性分析的準(zhǔn)確性,首先要分析二者的平衡點(diǎn)和平衡點(diǎn)穩(wěn)定性是否一致,然后通過數(shù)值計(jì)算對比二者所得軌跡的一致性。

3.2 脈沖未作用時(shí)的吸引域估計(jì)結(jié)果

以文獻(xiàn)[10]中提供的某120 mm迫彈作為無控彈平臺(tái),結(jié)構(gòu)參數(shù)如表2所示。

表2 某120 mm迫彈結(jié)構(gòu)參數(shù)[10]

由于末修迫彈是在彈道末段(如距離目標(biāo)800～1 000 m)進(jìn)行修正控制,其大部分彈道為無控彈道,末段啟控點(diǎn)參數(shù)可作為末段運(yùn)動(dòng)的初始條件。假設(shè)彈體不滾轉(zhuǎn)(P=0),以初速318 m/s、射角45°對上述120 mm迫彈進(jìn)行無控彈道計(jì)算,其至目標(biāo)斜距800 m處的彈道參數(shù)如表3所示,對應(yīng)的角運(yùn)動(dòng)方程參數(shù)如表4所示。

表3 彈道末段起始位置的彈道參數(shù)

表4 彈道末段起始位置的角運(yùn)動(dòng)方程參數(shù)

此時(shí),角運(yùn)動(dòng)方程式(27)具有唯一平衡點(diǎn)x*=[0.000 2 rad 0 0 0]T。在該平衡點(diǎn)進(jìn)行線性化,線性化系統(tǒng)特征方程的兩對特征根均具有負(fù)實(shí)部,則根據(jù)線性化理論,角運(yùn)動(dòng)在平衡點(diǎn)x*的鄰域內(nèi)可穩(wěn)定收斂到x*。狀態(tài)變量的數(shù)值積分結(jié)果如圖3所示。

圖3 原模型和多項(xiàng)式近似模型的相軌跡比較

如圖3所示,多項(xiàng)式近似系統(tǒng)和原系統(tǒng)軌跡是一致的,驗(yàn)證了上述多項(xiàng)式近似系統(tǒng)代替原系統(tǒng)進(jìn)行角運(yùn)動(dòng)穩(wěn)定性分析的準(zhǔn)確性。

對于彈箭的實(shí)際角運(yùn)動(dòng)而言,所允許的攻角和角速度是有范圍的,故需要找到具有實(shí)際意義的攻角和角速度穩(wěn)定范圍,即開展吸引域估計(jì)。

假設(shè)在選定點(diǎn)附近一段時(shí)間內(nèi),彈箭速度、飛行高度、彈道傾角的變化忽略不計(jì),則氣動(dòng)系數(shù)保持不變。將非線性氣動(dòng)力式(25)和幾何非線性近似式(32)代入角運(yùn)動(dòng)方程式(27),并通過坐標(biāo)轉(zhuǎn)換將平衡點(diǎn)移至原點(diǎn),之后可利用平方和規(guī)劃估計(jì)其吸引域。

利用V-s迭代算法,分別尋找二次和四次李雅普諾夫函數(shù),最終所得函數(shù)記為V2和V4:

V2=24.063 268α2+344.850 517αα′+24.063 267β2+344.850 514ββ′+279 665.08α′2+279 665.07β′2-0.009 189α-0.065 849α′

(33)

V4=0.291 549α4+8.517 460α3α′+0.583 178α2β2+8.467 055α2ββ′+1 204.838 384α2α′2+2 960.584 841α2β′2+8.466 879αβ2α′-3 496.704 7αβα′β′+4 404.953 287αα′3-0.000 430αα′2β′+4 305.759 8αα′β′2-0.000 508αβ′3+0.291 549β4+8.517 493β3β′+2 960.576 6β2α′2-1.018 918e-6β2α′β′+1 204.859 8β2β′2+0.000 402βα′3+4 309.150 272βα′2β′+0.000 544βα′β′2+4 404.883 9ββ′3+2 692 175.590 3α′4-8.036 145e-5α′3β′+5 401 623.866 9α′2β′2+0.000 749α′β′3+2 692 218.734 4β′4-0.000 242α3+0.018 099α2α′-0.000 243αβ2+0.028 630 1αββ′-0.074 852αα′2-6.852 724e-6αα′β′-2.955 652αβ′2-0.010 157β2α′

(34)

上述李雅普諾夫函數(shù)V2、V4等,必然滿足V(x)>0。這是由于V(x)>0為多項(xiàng)式正定的條件,在平方和規(guī)劃中已將其轉(zhuǎn)化成平方和約束,即V(x)-l1(x)∈Σn(體現(xiàn)在V-s算法的第3步),故V(x)>0本身即為問題的約束條件,求得的V2、V4等自然滿足該約束。計(jì)算所得吸引域分別為

(35)

其空間表達(dá)如圖4所示,其在各相平面的軌跡(剖面圖)如圖5所示。

圖4 吸引域估計(jì)結(jié)果的空間表示

圖5 吸引域估計(jì)結(jié)果的剖面圖

如圖4和圖5所示,基于平方和規(guī)劃的吸引域估計(jì)給出了很好的結(jié)果。需要說明的是,上述計(jì)算是在假設(shè)彈箭飛行環(huán)境不變的條件下得到,因而只在相對較短的一段彈道上成立。上述計(jì)算結(jié)果表明,平方和規(guī)劃方法用于確定角運(yùn)動(dòng)穩(wěn)定范圍是有效的。值得注意的是,由于式(26)表示的氣動(dòng)力僅考慮到立方次,如果能夠建立更為準(zhǔn)確的高階氣動(dòng)力模型,將得到更加準(zhǔn)確的結(jié)果。

3.3 脈沖作用下彈體穩(wěn)定的初始條件范圍

對于脈沖末修迫彈,發(fā)射后先做無控飛行,當(dāng)飛入預(yù)定區(qū)域后(彈道末端,如距離目標(biāo)800 m斜距處)啟控,進(jìn)行彈道修正以提高對目標(biāo)的射擊精度。由于采用的控制執(zhí)行機(jī)構(gòu)為脈沖發(fā)動(dòng)機(jī),在設(shè)計(jì)脈沖控制參數(shù)時(shí),需要確定脈沖大小J、脈沖作用的軸向位置Lx(距質(zhì)心,作用在質(zhì)心前為正)以及脈沖作用的周向位置φJ(rèn)(與彈軸系y軸的夾角,即脈沖方位角)。由于脈沖作用的時(shí)間極短(通常約幾毫秒至幾十毫秒),可認(rèn)為其作用前后僅引起攻角和攻角速度的變化。

設(shè)脈沖作用前(記為t=0-)攻角狀態(tài)為x0-=[α0,α′0,β0,β′0]T,脈沖作用后(記為t=0+)x0+=[α,α′,β,β′]T,則有

x0+=x0-+xJ

(36)

式中:xJ為脈沖作用引起的角運(yùn)動(dòng)狀態(tài)增量,

(37)

若要保證脈沖作用前后的角運(yùn)動(dòng)穩(wěn)定,需同時(shí)滿足:

(38)

則穩(wěn)定的角運(yùn)動(dòng)范圍為上述兩個(gè)集合的交集Ωx0-∩Ωx0+,設(shè)其仍具有吸引域所對應(yīng)的李雅普諾夫函數(shù)之形式,則該交集可近似為

Ω∩={x0-|V((x0-+x0+)/2)<γ∩}

(39)

式中:γ∩為函數(shù)V((x0-+x0+)/2)的水平集。因此,可以方便地利用定理2處理集合包含關(guān)系,并利用平方和規(guī)劃得到盡可能大的水平集γ∩。

考慮算例:脈沖沖量大小為J=60 N·s,作用軸向位置為Lx=0.083d(即質(zhì)心前0.083d)。將上述脈沖參數(shù)值代入式(38),并利用平方和規(guī)劃通過式(39)計(jì)算不同脈沖方位角φJ(rèn)所對應(yīng)的角運(yùn)動(dòng)穩(wěn)定范圍,即交集Ω∩的水平集γ∩。計(jì)算結(jié)果表明,對于不同的φJ(rèn),水平集γ∩的大小相同,但穩(wěn)定區(qū)域Ω∩的位置隨φJ(rèn)的變化而變化,如圖6所示。

圖6 交集估計(jì)Ω∩在α′-β′平面上的位置變化

在圖6所示α′-β′平面上,當(dāng)φJ(rèn)以45°為間隔從0～360°取值時(shí),紅色曲線為不同φJ(rèn)所對應(yīng)的Ωx0+,綠色曲線為穩(wěn)定角運(yùn)動(dòng)初始條件的估計(jì)Ω∩,顯然,Ωx0+的位置隨著φJ(rèn)的變化而變化,進(jìn)而Ω∩也隨著交集的變化而變化。

從圖6還可看出,存在一定的區(qū)域,在任意脈沖方位角下,彈箭角運(yùn)動(dòng)均穩(wěn)定。不妨將這樣的區(qū)域記為Ω?,且設(shè)其仍具有所對應(yīng)的李雅普諾夫函數(shù)之形式,即

Ω?={x0-|V(x0-)<γ?}

(40)

則仍然可以通過集合包含問題優(yōu)化計(jì)算得到水平集γ?的盡可能最大值。

由于上述問題沒有精確解,為驗(yàn)證平方和規(guī)劃所得結(jié)果的正確性,下面引入另一種計(jì)算方法。由于二次李雅普諾夫函數(shù)給出的吸引域恰為橢球體,故γ∩和γ?也可通過一定的幾何關(guān)系得到。設(shè)二次李雅普諾夫函數(shù)給出的吸引域Ω2D的半徑(即半長軸長)為rΩ,并定義

(41)

為脈沖增量xJ對應(yīng)的半徑,則易知脈沖作用下穩(wěn)定角運(yùn)動(dòng)范圍Ω∩的存在條件為

max (rJk/rΩk)<1,k=1,2,3,4

(42)

方便起見,下面將max (rJk/rΩk)簡記為max (rJ/rΩ)。

對于二次李雅普諾夫函數(shù),設(shè)橢球體Ω∩和Ω?對應(yīng)的半徑分別為r∩和r?,則根據(jù)幾何關(guān)系,可得到半徑與水平集之間的關(guān)系為

γ∩=max2(r∩/rΩ),γ?=max2(r?/rΩ)

(43)

各半徑之間的關(guān)系為

(44)

r?=2r∩-rΩ

(45)

同時(shí)可得到

(46)

且

(47)

上述關(guān)系意味著,可以通過優(yōu)化計(jì)算得到γ∩后,直接由式(46)計(jì)算得到水平集γ?;另一種方法是通過已知的rJ和rΩ直接計(jì)算γ∩和γ?。

根據(jù)式(44)第1式可得,若r∩存在,也即脈沖在某個(gè)方位角φJ(rèn)下作用時(shí),存在穩(wěn)定角運(yùn)動(dòng)范圍的條件為

max (rJ/rΩ)<2

(48)

根據(jù)式(44)第2式可得,若r?存在,即存在任意方位角下均能確保彈箭角運(yùn)動(dòng)穩(wěn)定的初始條件,那么

max (rJ/rΩ)<1

(49)

對于本節(jié)算例中的二次李雅普諾夫函數(shù)V2,其半徑為

(50)

當(dāng)脈沖沖量為J=60 N·s、作用位置Lx=0.083d時(shí),脈沖增量對應(yīng)的半徑為

(51)

利用平方和規(guī)劃方法處理集合包含關(guān)系得到的γ∩,進(jìn)而根據(jù)式(46)計(jì)算得到的γ?如表5所示,同時(shí)給出了根據(jù)式(47)直接計(jì)算的γ∩、γ?值。

表5 吸引域半徑計(jì)算結(jié)果

根據(jù)表5所示結(jié)果可見,優(yōu)化計(jì)算所得估計(jì)結(jié)果與幾何關(guān)系式(44)給出的結(jié)果非常接近,由此認(rèn)為二者互相驗(yàn)證了各自的準(zhǔn)確性。

在此基礎(chǔ)上,對于已經(jīng)過驗(yàn)證的平衡點(diǎn)吸引域Ω及其半徑rΩ,若給定脈沖大小J和作用位置Lx,根據(jù)式(48)和式(49),即可判斷此時(shí)是否存在穩(wěn)定的角運(yùn)動(dòng)范圍以及對任意脈沖方位角均穩(wěn)定的角運(yùn)動(dòng)范圍,并可根據(jù)式(47)計(jì)算出穩(wěn)定角運(yùn)動(dòng)范圍的水平集,進(jìn)而由式(39)和式(40)確定出角運(yùn)動(dòng)穩(wěn)定范圍的具體位置。

3.4 應(yīng)用:脈沖參數(shù)設(shè)計(jì)

上述結(jié)論可用于脈沖參數(shù)的設(shè)計(jì),即根據(jù)吸引域估計(jì)結(jié)果反向確定出脈沖參數(shù)的設(shè)計(jì)范圍,該范圍內(nèi)的任意脈沖作用均可使彈箭保持穩(wěn)定的角運(yùn)動(dòng)。

將脈沖增量對應(yīng)的吸引域半徑表達(dá)式(41)代入式(47),可得

(52)

即

(53)

根據(jù)表5,有γ?=0.068 6,代入式(53)即可得出滿足水平集γ?=0.068 6的脈沖參數(shù)J,Lx的所有可能組合(稱為有效脈沖參數(shù)組合),如圖7曲線與縱軸包圍的區(qū)域所示,縱軸為脈沖作用在彈軸上的位置。

圖7 滿足水平集γ?=0.068 6的脈沖參數(shù)組合

圖7中曲線與縱軸包圍的區(qū)域即為脈沖參數(shù)設(shè)計(jì)的有效范圍。由此可見,對于一定的脈沖軸向作用位置,均對應(yīng)一定的脈沖大小。這些脈沖參數(shù)的組合可保證,在任意脈沖方位角φJ(rèn)作用下,只要初始角運(yùn)動(dòng)狀態(tài)在Ω?范圍內(nèi),受控后的彈箭角運(yùn)動(dòng)仍可保持穩(wěn)定。顯然,上述有效脈沖參數(shù)組合,可直接用于非線性條件下的脈沖參數(shù)設(shè)計(jì)。

此外,當(dāng)γ?∈(0,1]時(shí),有效脈沖參數(shù)組合如圖8所示。豎線左邊給出了滿足式(53)第1式的脈沖范圍,對應(yīng)顏色的兩條曲線所包圍的區(qū)域給出了滿足式(53)第2式的脈沖參數(shù)組合范圍,二者的重合區(qū)域即為γ?所對應(yīng)的有效脈沖參數(shù)組合范圍。由圖8可見,當(dāng)γ?越小,有效脈沖參數(shù)組合所對應(yīng)的范圍越大。

圖8 不同γ?下的有效脈沖參數(shù)組合

將已知的rΩ(即式(50))代入式(43),計(jì)算可得不同γ?對應(yīng)的Ω?的半徑r?,進(jìn)而確定了α、β、α′、β′各自可取到的最大值,如圖9所示。由于角運(yùn)動(dòng)方程式(27)的對稱性,α和β的最大取值相等,α′和β′的最大取值相等。對于本節(jié)算例,當(dāng)γ?=0.068 6,α、β的最大值為3.06°,α′、β′的最大值為0.028°/cal。

圖9 γ?與Ω?中攻角和攻角速度最大值的對應(yīng)關(guān)系

4 結(jié)論

本文探討了平方和規(guī)劃方法在彈箭角運(yùn)動(dòng)非線性吸引域估計(jì)中的應(yīng)用。得到以下主要結(jié)論:

1)對于考慮高階氣動(dòng)力系數(shù)的無控尾翼彈非線性平面角運(yùn)動(dòng)吸引域估計(jì),通過構(gòu)造李雅普諾夫函數(shù),得到了高精度的吸引域估計(jì)結(jié)果,且李雅普諾夫函數(shù)的階數(shù)越高,吸引域估計(jì)值與確切值越接近,驗(yàn)證了方法的可行性和有效性。

2)針對考慮了諸非線性氣動(dòng)力系數(shù)的脈沖末修迫彈空間角運(yùn)動(dòng)吸引域估計(jì),通過構(gòu)建高精度多項(xiàng)式近似模型開展平方和規(guī)劃,對于采用二次李雅普諾夫函數(shù)的情形,吸引域估計(jì)結(jié)果與直接利用橢球幾何關(guān)系所得結(jié)果一致,驗(yàn)證了方法的準(zhǔn)確性。

3)由于脈沖作用主要影響彈箭的初始運(yùn)動(dòng)狀態(tài),故根據(jù)脈沖末修迫彈角運(yùn)動(dòng)的吸引域估計(jì),可反算并確定出脈沖沖量及其軸向作用位置的有效組合范圍(滿足非線性穩(wěn)定要求),由此得到非線性條件下脈沖參數(shù)設(shè)計(jì)的有效方法。

本文研究結(jié)果表明,平方和規(guī)劃方法為彈箭角運(yùn)動(dòng)吸引域估計(jì)、非線性條件下的彈箭參數(shù)設(shè)計(jì)等提供了一個(gè)有力的工具,為后續(xù)更為深入的應(yīng)用研究奠定了基礎(chǔ)。