鄭小洪，韓維，賈忠湖，郭衛(wèi)剛，張?jiān)?/p>

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

0 引言

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

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

從目前的研究成果看，大部分都集中于離散模型研究。離散模型在一定程度上能夠反映拖索的構(gòu)型，但是為了進(jìn)一步深入研究航空拖索的動力學(xué)特性，建立航空拖索的連續(xù)的動力學(xué)模型是非常必要的。本文在文獻(xiàn)［7-8］基礎(chǔ)上，建立機(jī)載甚低頻拖曳天線的動力學(xué)模型，運(yùn)用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)，它們具有相同的原點(diǎ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點(diǎn)位移為:

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

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

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

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

由于機(jī)載甚低頻拖曳天線末端連著一個(gè)錐袋，因此拖曳天線系統(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)時(shí)，即位移U=0，將U=0、式（1）和式（3）代入式（2）并無因次化得:

其中:

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

穩(wěn)態(tài)時(shí)，天線末端錐袋的位移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é)模型屬于兩點(diǎn)邊值問題。目前已有許多學(xué)者對邊值問題進(jìn)行了研究。常見的主要有離散法和初值方法。離散法最直接的就是運(yùn)用差分法［9］和變分法［10-11］進(jìn)行求解；而初值方法目前最主要的就是打靶法［12］。本文利用 Newton-Raphson 迭代法［13］對打靶法試探解進(jìn)行修正，該方法具有收斂速度快等優(yōu)點(diǎn)?？紤]如下一類兩點(diǎn)邊值問題:

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

通過式（14）便可以在一定誤差下將該邊值問題轉(zhuǎn)化為初值問題（本文中取截?cái)嗾`差為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)化為初值問題，便可以運(yùn)用Runge-Kutta法進(jìn)行求解。

4 數(shù)值仿真

式（10）和式（11）提供了天線末端處的6個(gè)初始條件，h0，xd作為假定初值。運(yùn)用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ī)盤旋時(shí)的傾斜角，Vtrue為載機(jī)的真實(shí)空速，則載機(jī)的旋轉(zhuǎn)半徑和旋轉(zhuǎn)角速度為:

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

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

4.1 仿真實(shí)例1

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

由圖2和圖3可知，本文建立的機(jī)載甚低頻拖曳天線動力學(xué)模型是可信的，運(yùn)用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 仿真實(shí)例2

仿真條件:φ=36°，CD分別為1.03和0.69，其他條件同仿真實(shí)例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ī)的連接點(diǎn)附近影響較大，法向空氣阻力系數(shù)越大，天線張力越小。

4.3 仿真實(shí)例3

仿真條件:φ=36°，M分別為40 kg和30 kg，其他條件同仿真實(shí)例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é)模型運(yùn)用Newton-Raphson迭代進(jìn)行求解。通過數(shù)值計(jì)算分析了法向空氣阻力系數(shù)和天線錐袋質(zhì)量對天線張力和垂直度的影響。本文的仿真結(jié)果對機(jī)載甚低頻拖曳天線系統(tǒng)的設(shè)計(jì)有一定的指導(dǎo)作用。

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