變體航行器動(dòng)力學(xué)建模與仿真

顏奇民，胡俊華，陳國(guó)明，劉安，張強(qiáng)

（空軍工程大學(xué)航空工程學(xué)院，西安 710038）

變體航行器是指在航行過(guò)程中能夠改變自身結(jié)構(gòu)特征以適應(yīng)外界環(huán)境和任務(wù)需求的一類航行器，其涵蓋范圍十分廣泛，包括變體水下航行器、變體飛機(jī)、跨介質(zhì)航行器、變體航天器等，由于變體航行器可以實(shí)現(xiàn)不同構(gòu)型間的切換，其相比固定構(gòu)型航行器功能更加豐富、任務(wù)范圍更廣［1］。

變體航行器的動(dòng)力學(xué)建模和仿真往往十分復(fù)雜，航行器變體過(guò)程中，自身的各項(xiàng)參數(shù)會(huì)發(fā)生變化，如變體飛行器的氣動(dòng)參數(shù)、變體水下航行器的水動(dòng)參數(shù)，同時(shí)，航行器變體過(guò)程中自身結(jié)構(gòu)的調(diào)整也會(huì)對(duì)航行器的運(yùn)動(dòng)狀態(tài)產(chǎn)生影響；針對(duì)變體航行器的動(dòng)力學(xué)建模問(wèn)題，許來(lái)斌等［2］基于牛頓歐拉方程，推導(dǎo)了變體飛機(jī)運(yùn)動(dòng)方程并進(jìn)行了俯仰和轉(zhuǎn)向運(yùn)動(dòng)仿真；王志剛等［3］針對(duì)一種后掠變展長(zhǎng)組合變形飛行器模型，采用多剛體理論建立了描述飛行器變形過(guò)程的動(dòng)力學(xué)模型；Seigler等［4-6］建立了固定于飛機(jī)機(jī)身的機(jī)體坐標(biāo)軸系，采用笛卡兒矢量描述機(jī)體各點(diǎn)的位置，推導(dǎo)了變體飛機(jī)的一般動(dòng)態(tài)特性；陳偉和馮高鵬［7］將變體飛機(jī)視為一個(gè)質(zhì)點(diǎn)系，得到了變體飛機(jī)的力方程和力矩方程；譚駿怡等［8］結(jié)合有限元思想和分析力學(xué)方法分別分析了剛性翼和柔性翼的動(dòng)力學(xué)特性，并構(gòu)建了變體飛行耦合動(dòng)力學(xué)模型；上述研究分別建立了針對(duì)變體飛機(jī)的變體動(dòng)力學(xué)模型，但是存在計(jì)算復(fù)雜、適用范圍小等問(wèn)題。

本文在陳偉和馮高鵬［7］的研究基礎(chǔ)上，采用積分思想，將變體飛機(jī)視為由無(wú)數(shù)質(zhì)點(diǎn)組成的質(zhì)點(diǎn)系，建立固定于機(jī)體的機(jī)體坐標(biāo)系，分析每個(gè)質(zhì)點(diǎn)的運(yùn)動(dòng)，推導(dǎo)了變體航行器動(dòng)力學(xué)擴(kuò)展方程；對(duì)一種具有共形半環(huán)翼的變體航行器進(jìn)行了不同變體速率下的縱向面內(nèi)開(kāi)環(huán)仿真，針對(duì)建立的動(dòng)力學(xué)模型中部分附加力和附加力矩計(jì)算困難的問(wèn)題，根據(jù)該航行器的已知變體規(guī)律，計(jì)算了不同時(shí)間節(jié)點(diǎn)航行器的靜矩并擬合得到靜矩的一階和二階導(dǎo)數(shù)，附加力矩中非靜矩的未知項(xiàng)換算為包含靜矩的二階導(dǎo)數(shù)的已知項(xiàng)進(jìn)行計(jì)算，簡(jiǎn)化了計(jì)算過(guò)程。

1 變體航行器動(dòng)力學(xué)建模

1.1 坐標(biāo)系建立

以航行器未變體時(shí)的質(zhì)心作為原點(diǎn)建立機(jī)體坐標(biāo)系O1x1y1z1，地面坐標(biāo)系為O0x0y0z0，如圖1所示［9］。

圖1 坐標(biāo)系建立［9］Fig.1 Coordinate system setup［9］

式中：sinγ，cosγ，sin?，cos?，sinψ，cosψ分別為航行器滾轉(zhuǎn)角γ、俯仰角?和偏航角ψ的正弦和余弦值。

1.2 變體航行器運(yùn)動(dòng)方程

將航行器視為由n個(gè)質(zhì)點(diǎn)組成的整體，首先，在地面坐標(biāo)系中描述航行器的n個(gè)質(zhì)點(diǎn)的位置，表達(dá)式為

在此，引入運(yùn)算符“～”表示向量元素組成的反對(duì)稱矩陣，其規(guī)則為：若x＝xii+xjj+xkk（i，j，k表示沿標(biāo)準(zhǔn)坐標(biāo)系3個(gè)坐標(biāo)軸的單位矢量），則

式中：Ω ＝［ωxb，ωyb，ωzb］T為機(jī)體坐標(biāo)系相對(duì)于地面坐標(biāo)系的角速度在機(jī)體坐標(biāo)系中的表示；ω為角速度矢量，ωxb、ωyb、ωzb分別為該矢量沿機(jī)體坐標(biāo)系x軸、y軸、z軸的分量。

對(duì)式（2）兩端求導(dǎo)得

式（21）即普通航行器的動(dòng)力學(xué)方程，因此，本文得到的動(dòng)力學(xué)方程（20）即變體航行器的動(dòng)力學(xué)擴(kuò)展方程，可以適用于發(fā)生任意可知形變的變體航行器，如變后掠航行器、折疊翼航行器、伸縮翼航行器等。

轉(zhuǎn)動(dòng)慣量張量I的具體表達(dá)式為［13］

2 共形半環(huán)翼航行器變體描述

本文研究的共形半環(huán)翼航行器是一種具有共形半環(huán)翼布局的跨介質(zhì)航行器，其空中構(gòu)型和水下構(gòu)型分別如圖2（a）和圖2（b）所示［14-15］。

圖2 共形半環(huán)翼航行器構(gòu)型［14-15］Fig.2 Configuration of conformal semi-ring wing vehicle［14-15］

該航行器構(gòu)型主要針對(duì)水空跨越任務(wù)提出，航行器在水中航行時(shí)，共形半環(huán)翼收起，包裹于機(jī)體上，其構(gòu)型與魚(yú)雷相似，航行器沖出水面后，共形半環(huán)翼展開(kāi)為空中構(gòu)型，迅速加速獲取升力起飛，當(dāng)航行器準(zhǔn)備再次入水時(shí)，機(jī)翼在空中收起，轉(zhuǎn)換為水下構(gòu)型，并減速俯沖入水。

共形半環(huán)翼包含2層翼面，其中，下層翼為剛性翼，上層翼為可彎曲的柔性翼，上、下翼面在翼梢端點(diǎn)處連接，航行器變體時(shí)，作動(dòng)器控制下層翼轉(zhuǎn)動(dòng)，進(jìn)而帶動(dòng)上層翼發(fā)生轉(zhuǎn)動(dòng)和彈性形變，其變體過(guò)程的翼型正視圖如圖3所示。

圖3 共形半環(huán)翼航行器變體過(guò)程Fig.3 Process of variant of conformal semi-ring wing vehicle

由于航行器進(jìn)行對(duì)稱變體，為準(zhǔn)確描述航行器變體情況，取航行器正向單側(cè)剖面，如圖4所示。

圖4 航行器正向單側(cè)剖面圖Fig.4 Forward unilateral profile of vehicle

圖4中，角θ和角μ分別用以描述下層翼和上層翼的轉(zhuǎn)動(dòng)角度，且μ的值與θ的值一一對(duì)應(yīng)。

本文的仿真中，設(shè)計(jì)航行器以不同的變體速度在θ＝0和θ＝3π／8之間進(jìn)行收回或展開(kāi)，變體規(guī)律為

式中：˙ωθ為下層翼的角加速度，ωθ為下層翼的旋轉(zhuǎn)角速度。

3 航行器變體氣動(dòng)參數(shù)

采用Fluent流體力學(xué)仿真軟件對(duì)航行器的升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù)進(jìn)行了解算，得到了航行器在不同攻角、不同變體角度的氣動(dòng)系數(shù)。根據(jù)仿真結(jié)果和研究分析，共形半環(huán)翼航行器的安全攻角范圍為［-π／18，π／9］，因此，本節(jié)給出［-π／18，π／9］攻角范圍和［0，3π／4］下層翼變體角度范圍對(duì)應(yīng)的氣動(dòng)參數(shù)（圖中已將弧度折算為角度），如圖5所示。

圖5 航行器變體氣動(dòng)參數(shù)Fig.5 Aerodynamic parameters of vehicle in variant phase

4 附加力與附加力矩計(jì)算

航行器機(jī)翼變形對(duì)航行器運(yùn)動(dòng)的影響具體表現(xiàn)為由于變形產(chǎn)生的附加力F1～F4及力矩M1～M3，根據(jù)附加力和附加力矩表達(dá)式，計(jì)算附加力需要求得的物理量除基本的運(yùn)動(dòng)參數(shù)外，還包括靜矩S及其一階、二階導(dǎo)數(shù)，以及附加力矩

本文航行器僅機(jī)翼發(fā)生變體，因此航行器靜矩全部來(lái)源于機(jī)翼質(zhì)心的變化，文中采取積分方法計(jì)算機(jī)翼產(chǎn)生的靜矩，如圖4（b）所示。

分別計(jì)算上、下層翼產(chǎn)生的靜矩，可得

式中：Sy_u和Sy_d分別為上、下層翼產(chǎn)生的靜矩；yo″和yo′分別為上、下層翼圓心的縱坐標(biāo)；μo″、ro″、ρup分別為上層翼的圓心角、半徑和單位弧度的質(zhì)量；rt、ρdown分別為下層翼的半徑和單位弧度的質(zhì)量。根據(jù)航行器變體規(guī)律，這些參數(shù)均可求得。

本文在設(shè)定kθ后，可求取航行器變體總用時(shí)tw，以tw／20為步長(zhǎng)，計(jì)算不同時(shí)刻的Sy，并利用MATLAB擬合功能［16］，快速求解S·y、S··y的近似值，本文航行器機(jī)翼由完全閉合至完全展開(kāi)過(guò)程為例，其靜矩及靜矩的一階、二階導(dǎo)數(shù)變化情況如圖6所示。

圖6 靜矩及其導(dǎo)數(shù)變化曲線Fig.6 Variation curve of static moment and its derivative

5 不同變體速率開(kāi)環(huán)仿真分析

為研究開(kāi)環(huán)配平狀態(tài)下不同的航行器變體速度對(duì)航行器運(yùn)動(dòng)狀態(tài)的影響，設(shè)定航行器的飛行速度為0.8 Ma，分別對(duì)航行器下層翼角度θ＝0和θ＝3π／8兩種狀態(tài)進(jìn)行配平，得到兩組配平狀態(tài)［17］，分別記為A和B，其具體參數(shù)如表1所示。

令航行器以不同的變體速度在θ＝0和θ＝3π／8之間進(jìn)行收回或展開(kāi)，分別設(shè)置變體時(shí)間為8，4，2，1 s，仿真結(jié)果如圖7和圖8所示。

由圖7和圖8可知，在配平狀態(tài)下，航行器機(jī)翼由θ＝0收回至θ＝3π／8狀態(tài)，會(huì)導(dǎo)致航行器產(chǎn)生俯沖運(yùn)動(dòng)，速度增加，反之，航行器機(jī)翼展開(kāi)會(huì)導(dǎo)致爬升運(yùn)動(dòng)，速度降低，同時(shí)，航行器變體速度越快，運(yùn)動(dòng)狀態(tài)變化越劇烈。

圖7 機(jī)翼收回仿真結(jié)果Fig.7 Simulation results of retracting wing

圖8 機(jī)翼展開(kāi)仿真結(jié)果Fig.8 Simulation results of unfolding wing

2種變體情況都會(huì)使航行器的攻角和俯仰角發(fā)生振蕩。機(jī)翼收回時(shí)，攻角總體有所上升，俯仰角先上升，而后有下降趨勢(shì)；機(jī)翼展開(kāi)仿真中，總體趨勢(shì)與收回仿真相反。這是因?yàn)楹叫衅鳈C(jī)翼在收回和展開(kāi)過(guò)程中引起氣動(dòng)參數(shù)發(fā)生變化，在機(jī)翼收回過(guò)程中，航行器升力、阻力系數(shù)和負(fù)向俯仰力矩均減小，導(dǎo)致航行器產(chǎn)生俯沖，速度增大，俯仰角增大，從而導(dǎo)致航行器攻角迅速增大，由于航行器具有縱向靜穩(wěn)定性，攻角的增大會(huì)引起負(fù)向俯仰力矩增加，航行器俯仰角和攻角減小，且在縱向靜穩(wěn)定性的影響下，隨后航行器的攻角和俯仰角會(huì)產(chǎn)生逐步減小的振蕩。

航行器機(jī)翼收回過(guò)程和展開(kāi)過(guò)程中附加力和附加力矩變化方向相反而趨勢(shì)相似，因此只給出航行器機(jī)翼收回過(guò)程中的附加力和附加力矩的變化情況，如圖9所示。

圖9 附加力與附加力矩變化曲線Fig.9 Variation curve of additional forces and additional moments

基于本文航行器的結(jié)構(gòu)，其變體產(chǎn)生的靜矩均沿y1軸，根據(jù)附加力與附加力矩的表達(dá)式，在航行器俯仰運(yùn)動(dòng)中，F(xiàn)1與F4沿y1軸方向，F(xiàn)2與F3沿x1軸方向，M1、M2、M3均沿z1軸方向，因此直接給出各個(gè)附加力與附加力矩沿相應(yīng)軸向的值。

由圖9可得，配平狀態(tài)下航行器進(jìn)行變體時(shí)，變體速度越快，產(chǎn)生的附加力與附加力矩越大。在本次仿真中，F(xiàn)4的值遠(yuǎn)大于其他附加力，這是因?yàn)樵摲抡嬷泻叫衅鞯慕撬俣群徒羌铀俣染^小，導(dǎo)致F1～F3較小，F(xiàn)4為航行器靜矩的二階導(dǎo)，可以理解為機(jī)翼變體過(guò)程中機(jī)翼上質(zhì)點(diǎn)加速對(duì)機(jī)身產(chǎn)生的反作用力。

附加力矩中，M3的值遠(yuǎn)大于其他附加力矩，這是由于航行器的角速度、加速度較小，導(dǎo)致M1較小，同時(shí)，本文航行器變體過(guò)程中慣量矩陣的變化率很低，導(dǎo)致M2的值很小，而M3可以理解為機(jī)翼變體過(guò)程中的機(jī)翼上質(zhì)點(diǎn)加、減速產(chǎn)生反作用力偏離質(zhì)心所產(chǎn)生的力矩，變體越快，產(chǎn)生的M3力矩越大。

6 結(jié) 論

本文建立了變體航行器的動(dòng)力學(xué)拓展方程，針對(duì)變體規(guī)律已知的共形半環(huán)翼航行器，對(duì)不同時(shí)刻的靜矩及轉(zhuǎn)動(dòng)慣量進(jìn)行擬合和求導(dǎo)，求得其變體附加力和附加力矩，該方法也可以適用于其他具有已知變體規(guī)律的變體航行器的動(dòng)力學(xué)建模，包括變體飛行器、變體航天器等。

共形半環(huán)翼航行器變體飛行的開(kāi)環(huán)仿真結(jié)果表明，航行器變體產(chǎn)生的附加力和附加力矩大小與變體速率正相關(guān)，其中，由靜矩變化引起的F4和M3影響最大；由于本文航行器機(jī)翼在整機(jī)中占比小，附加力和附加力矩遠(yuǎn)小于氣動(dòng)力，引起航行器運(yùn)動(dòng)狀態(tài)變化的主要因素仍然是變體導(dǎo)致的氣動(dòng)參數(shù)的變化；對(duì)于大幅度變體的航行器，則必須考慮附加力和附加力矩對(duì)運(yùn)動(dòng)的影響。