機載拖曳天線動力學建模與仿真

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

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

0 引言

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

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

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

1 甚低頻拖曳天線的建模

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

如圖1所示，S1和S2分別為拖曳天線自然無伸張狀態(tài)下、穩(wěn)定狀態(tài)和動態(tài)平衡狀態(tài)下甚低頻天線任意位置E的弧坐標，故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）旋轉坐標系下沿Ox軸、Oy軸和Oz軸方向的分量大小。

圖1 拖曳天線回轉示意圖Fig.1 Rotary diagram of the trailing antenna

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

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

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

式中，W，M分別為錐袋重力和質量；CDD為錐袋的空氣阻力系數(shù)；SD為錐袋的橫截面積；W'，M'分別為空氣對錐袋產生的浮力和由于空氣使得錐袋增加的附加質量；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)動力學模型

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

其中:

式中，a為載機盤旋半徑。

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

令 dr/ds=τ/p，則:

由式（5）可以得到:

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

（1）在s=0處

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

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

由式（7）可以得到:

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

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

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

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

將式（13）轉化為邊值問題得:

對式（15）關于βj求導，得:

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

4 數(shù)值仿真

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

由于載機的飛行高度低于11 km，因此天線的活動區(qū)域屬于對流層，其空氣密度隨高度變化為:

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

4.1 仿真實例1

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

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

圖2 天線垂直坐標分布情況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分別為天線垂直坐標和天線張力隨天線長度分布情況的仿真結果。

圖4 天線垂直坐標分布情況Fig.4 Vertical coordinate of the antenna

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

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

4.3 仿真實例3

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

圖6 天線垂直坐標分布情況Fig.6 Vertical coordinate of the antenna

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

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

5 結束語

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

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