機(jī)載拖曳天線動力學(xué)建模與仿真

鄭小洪，韓維，賈忠湖，郭衛(wèi)剛，張原

（1.海軍航空工程學(xué)院飛行器工程系，山東煙臺 264001；2.海軍航空兵學(xué)院模擬中心，遼寧 葫蘆島 125001）

0 引言

機(jī)載甚低頻拖曳系統(tǒng)的拖曳天線在工作過程中的穩(wěn)態(tài)構(gòu)型及垂直度（指拖曳天線在鉛垂線上的投影長度與拖曳天線總長度的比值）是檢驗天線能否正常工作的主要依據(jù)，是設(shè)計人員必須著重考慮的因素。出于對拖曳天線安全使用考慮，天線最大張力也是該系統(tǒng)需要關(guān)注的因素。為保證天線的有效正常工作狀態(tài)，載機(jī)必須以一定狀態(tài)進(jìn)行穩(wěn)定盤旋飛行［1］。

早在1969年，NADC就開始編寫計算機(jī)程序，用來計算TACAMO拖曳天線的平衡狀態(tài)形狀，其所建模型的缺點是只有在良好的初始條件下才能獲得較好的收斂效果［2］。文獻(xiàn)［3］在拖索上建立局部坐標(biāo)系，并使用有限元法思想將拖纜離散成多個分段，針對任一個分段建立平衡方程，并運用中心差分法進(jìn)行求解過程。近年來，文獻(xiàn)［4-6］運用多體動力學(xué)建模方法對航空拖索進(jìn)行建模，對航空拖索構(gòu)型進(jìn)行了一定的研究。

從目前的研究成果看，大部分都集中于離散模型研究。離散模型在一定程度上能夠反映拖索的構(gòu)型，但是為了進(jìn)一步深入研究航空拖索的動力學(xué)特性，建立航空拖索的連續(xù)的動力學(xué)模型是非常必要的。本文在文獻(xiàn)［7-8］基礎(chǔ)上，建立機(jī)載甚低頻拖曳天線的動力學(xué)模型，運用Newton-Raphson迭代法對甚低頻拖曳天線的動力學(xué)模型進(jìn)行求解。

1 甚低頻拖曳天線的建模

在建立其動力學(xué)模型前，首先要選定合適的空間坐標(biāo)系，如圖1建立慣性坐標(biāo)系Ox'y'z'和旋轉(zhuǎn)直角坐標(biāo)系Oxyz，Oxy平面為海平面。旋轉(zhuǎn)坐標(biāo)系繞縱軸坐標(biāo)軸以載機(jī)旋轉(zhuǎn)角速度Ω旋轉(zhuǎn)，它們具有相同的原點。設(shè)i'，j'，k'分別為慣性坐標(biāo)系下沿Ox'軸、Oy'軸和Oz'軸方向的單位矢量；i，j，k分別為旋轉(zhuǎn)坐標(biāo)系下沿Ox軸、Oy軸和Oz軸方向的單位矢量；XD為天線末端錐袋在旋轉(zhuǎn)坐標(biāo)系下的Ox軸坐標(biāo)，HD為其Oz軸坐標(biāo)。

如圖1所示，S1和S2分別為拖曳天線自然無伸張狀態(tài)下、穩(wěn)定狀態(tài)和動態(tài)平衡狀態(tài)下甚低頻天線任意位置E的弧坐標(biāo)，故E點位移為:

式中，U（S0，T）為E點從穩(wěn)定狀態(tài)到T時刻動態(tài)平衡狀態(tài)的位移；R2為動態(tài)平衡下圖中O點到E點在T時刻的距離；R1為穩(wěn)定狀態(tài)下O點到E點的距離；U1，U2，U3為U（S0，T）旋轉(zhuǎn)坐標(biāo)系下沿Ox軸、Oy軸和Oz軸方向的分量大小。

圖1 拖曳天線回轉(zhuǎn)示意圖Fig.1 Rotary diagram of the trailing antenna

由于切向阻力很?。?］，模型中只考慮法向空氣阻力。根據(jù)牛頓第二定律，天線的動力學(xué)方程為:

式中，P為天線拉力的數(shù)值大??；ρb為天線單位長度質(zhì)量密度；ρa(bǔ)為空氣密度；A0為天線的橫截面積；CD為法向空氣阻力系數(shù)；D為天線的直徑；為法向速度。

由于機(jī)載甚低頻拖曳天線末端連著一個錐袋，因此拖曳天線系統(tǒng)的固有邊界條件為:

式中，W，M分別為錐袋重力和質(zhì)量；CDD為錐袋的空氣阻力系數(shù)；SD為錐袋的橫截面積；W'，M'分別為空氣對錐袋產(chǎn)生的浮力和由于空氣使得錐袋增加的附加質(zhì)量；VD=?UD/?T+Ωk×（UD+R1）；UD=U（0，T）；aD=Ω2k×（k×R1）+Ω2k×（k×UD）+2Ω×（?UD/?T）+ ?2UD/?T2。

2 穩(wěn)態(tài)動力學(xué)模型

在穩(wěn)態(tài)時，即位移U=0，將U=0、式（1）和式（3）代入式（2）并無因次化得:

其中:

式中，a為載機(jī)盤旋半徑。

穩(wěn)態(tài)時，天線末端錐袋的位移UD=0，將其代入式（4）并無因次化得:

令 dr/ds=τ/p，則:

由式（5）可以得到:

上述無因次化模型的邊界條件如下:

（1）在s=0處

（2）在s=L/a=l處

式中，HP為載機(jī)盤旋高度；hp為無因次盤旋高度。

由式（7）可以得到:

3 基于Newton-Raphson迭代的邊值問題求解算法

由上文可知，機(jī)載甚低頻拖曳天線穩(wěn)態(tài)動力學(xué)模型屬于兩點邊值問題。目前已有許多學(xué)者對邊值問題進(jìn)行了研究。常見的主要有離散法和初值方法。離散法最直接的就是運用差分法［9］和變分法［10-11］進(jìn)行求解；而初值方法目前最主要的就是打靶法［12］。本文利用 Newton-Raphson 迭代法［13］對打靶法試探解進(jìn)行修正，該方法具有收斂速度快等優(yōu)點?？紤]如下一類兩點邊值問題:

式中，Gβ（y（1，β））為G（y（1，β））的雅克比矩陣。

通過式（14）便可以在一定誤差下將該邊值問題轉(zhuǎn)化為初值問題（本文中取截斷誤差為10－5）。因此問題的關(guān)鍵轉(zhuǎn)化為求解Gβ（y（1，β）），而Gβ（k）（y（1，β（k）））=Gy（y（1，β（k）））yβ（k）（1，β（k）），顯然Gy（y（1，β（k）））是可以求出解析式的，對于yβ（k）（1，β（k））則需要進(jìn)一步求解。

將式（13）轉(zhuǎn)化為邊值問題得:

對式（15）關(guān)于βj求導(dǎo)，得:

由式（16）可以求出yβ（k）（1，β（k））。由此將該邊值問題轉(zhuǎn)化為初值問題，便可以運用Runge-Kutta法進(jìn)行求解。

4 數(shù)值仿真

式（10）和式（11）提供了天線末端處的6個初始條件，h0，xd作為假定初值。運用Newton-Raphson迭代法進(jìn)行求解上述常微分方程式（9）。某甚低頻天線長度L=6 184.4 m，載機(jī)飛行高度HP=5 585.5 m，載機(jī)空速為 80.253 3 m/s［14］。由于空氣密度較小，暫且不考慮空氣密度和附加質(zhì)量的影響，因此ρ=1，w=m。設(shè)φ為載機(jī)盤旋時的傾斜角，Vtrue為載機(jī)的真實空速，則載機(jī)的旋轉(zhuǎn)半徑和旋轉(zhuǎn)角速度為:

由于載機(jī)的飛行高度低于11 km，因此天線的活動區(qū)域?qū)儆趯α鲗?，其空氣密度隨高度變化為:

式中，ρ0為海平面的空氣密度。

4.1 仿真實例1

仿真條件:ρb=0.093 4 kg/m，CD=1.03，CDD=0.41，M=37.2 kg，φ=34°。圖2 和圖3 分別為天線垂直坐標(biāo)和天線張力隨天線長度分布情況的仿真結(jié)果，顯然與AD報告［14］仿真計算結(jié)果是相符的。

由圖2和圖3可知，本文建立的機(jī)載甚低頻拖曳天線動力學(xué)模型是可信的，運用Newton-Raphson迭代法對該動力學(xué)模型進(jìn)行求解也是可行的。

圖2 天線垂直坐標(biāo)分布情況Fig.2 Vertical coordinate of the antenna

圖3 天線張力分布情況Fig.3 Distribution of the antenna tension

4.2 仿真實例2

仿真條件:φ=36°，CD分別為1.03和0.69，其他條件同仿真實例1。圖4和圖5分別為天線垂直坐標(biāo)和天線張力隨天線長度分布情況的仿真結(jié)果。

圖4 天線垂直坐標(biāo)分布情況Fig.4 Vertical coordinate of the antenna

圖5 天線張力分布情況Fig.5 Distribution of the antenna tension

由圖4和圖5可知，天線的垂直度與法向空氣阻力系數(shù)相關(guān)，法向空氣阻力系數(shù)越大，天線垂直度越小。法向空氣阻力系數(shù)對天線末端的張力影響較小，但是對于天線與載機(jī)的連接點附近影響較大，法向空氣阻力系數(shù)越大，天線張力越小。

4.3 仿真實例3

仿真條件:φ=36°，M分別為40 kg和30 kg，其他條件同仿真實例1。圖6和圖7分別為天線垂直坐標(biāo)和天線張力隨天線長度分布情況的仿真結(jié)果。

圖6 天線垂直坐標(biāo)分布情況Fig.6 Vertical coordinate of the antenna

圖7 天線張力分布情況Fig.7 Distribution of the antenna tension

由圖6和圖7可知，在該飛行條件下，天線末端錐袋質(zhì)量越大，天線的垂直度越大，但是天線的張力也隨著錐袋質(zhì)量變大而明顯變大。考慮到天線的強(qiáng)度，因此不能通過單一增加錐袋質(zhì)量來增加天線的垂直度。

5 結(jié)束語

本文以機(jī)載甚低頻拖曳天線為研究對象，建立了機(jī)載甚低頻拖曳天線動力學(xué)模型，并對其穩(wěn)態(tài)動力學(xué)模型運用Newton-Raphson迭代進(jìn)行求解。通過數(shù)值計算分析了法向空氣阻力系數(shù)和天線錐袋質(zhì)量對天線張力和垂直度的影響。本文的仿真結(jié)果對機(jī)載甚低頻拖曳天線系統(tǒng)的設(shè)計有一定的指導(dǎo)作用。

［1］韓維，侯志強(qiáng).機(jī)載拖曳系統(tǒng)的建模和若干動力學(xué)問題［C］//飛行力學(xué)與飛行試驗學(xué)術(shù)交流年會論文集（2005）.西安:中國飛行試驗研究院飛行力學(xué)雜志社，2005:1-5.

［2］Huang S L.Mathematical model for long cable by orbiting aircraft［R］.NADC-AM-6849，1969.

［3］James M C，Dr Louis V S.Dynamic modeling of a trailing wire towed by an orbiting aircraft［R］.AIAA-93-3663-CP，1993.

［4］Paul W，Pavel T.Periodic solutions for flexible cable-body systems towed in circular and elliptical paths［C］//AIAA Atmospheric Flight Mechanics Conference and Exhibit.Keystone，Colorado，American Institute of Aeronautics and Astronautics，2006.

［5］Paul W，Pavel T.A study on the transitional dynamics of a towed-circular aerial cable system［C］//AIAA Atmospheric Flight Mechanics Conference and Exhibit.San Francisco，California，American Institute of Aeronautics and Astronautics，2005.

［6］Paul W，Daniel S，Pavel T.Motion planning for an aerialtowed cable system［C］//AIAA Atmospheric Flight Mechanics Conference and Exhibit.San Francisco，California，American Institute of Aeronautics and Astronautics，2005.

［7］Zhu F，Rahn C D.Stability analysis of a circularly towed cable-body system ［J］.Journal of Sound and Vibration，1998，217（3）:435-452.

［8］Zhu F，Sharma R，Rahn C D.Vibrations of ballooning elastic stings［J］.Journal of Applied Mechanics，1997，64（3）:676-683.

［9］Saadatmandi A，F(xiàn)arsangi J A.Chebyshev finite difference method for a nonlinear system of second order boundary value problems［J］.Applied Mathematics and Computation，2007，192（2）:586-591.

［10］Lu J F.Variational iteration method for solving a nonlinear system of second-order boundary value problems［J］.Computers and Mathematics with Applications，2007，54（4）:1133-1138.

［11］Geng F Z，Cui M G.Solving a nonlinear system of second order boundary value problems［J］.Mathematical Analysis and applications，2007，327（2）:1167-1181.

［12］Yuki N，Satoshi T.On the existence of multiple solutions of the boundary value problem for nonlinear second-order differential equations ［J］.Nonlinear Analysis，2004，55（6）:919-935.

［13］Jeffery J L.Number analysis and scientific computation［M］.Beijing:Tsinghua University Press，2008.

［14］James M C.Modeling and control of a trailing wire antenna towed by an orbiting aircraft［R］.ADA256450，1992.