机载拖曳天线动力学建模与仿真

郑小洪，韩维，贾忠湖，郭卫刚，张原

（1.海军航空工程学院飞行器工程系，山东烟台 264001；2.海军航空兵学院模拟中心，辽宁 葫芦岛 125001）

0 引言

机载甚低频拖曳系统的拖曳天线在工作过程中的稳态构型及垂直度（指拖曳天线在铅垂线上的投影长度与拖曳天线总长度的比值）是检验天线能否正常工作的主要依据，是设计人员必须着重考虑的因素。出于对拖曳天线安全使用考虑，天线最大张力也是该系统需要关注的因素。为保证天线的有效正常工作状态，载机必须以一定状态进行稳定盘旋飞行［1］。

早在1969年，NADC就开始编写计算机程序，用来计算TACAMO拖曳天线的平衡状态形状，其所建模型的缺点是只有在良好的初始条件下才能获得较好的收敛效果［2］。文献［3］在拖索上建立局部坐标系，并使用有限元法思想将拖缆离散成多个分段，针对任一个分段建立平衡方程，并运用中心差分法进行求解过程。近年来，文献［4-6］运用多体动力学建模方法对航空拖索进行建模，对航空拖索构型进行了一定的研究。

从目前的研究成果看，大部分都集中于离散模型研究。离散模型在一定程度上能够反映拖索的构型，但是为了进一步深入研究航空拖索的动力学特性，建立航空拖索的连续的动力学模型是非常必要的。本文在文献［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分别为拖曳天线自然无伸张状态下、稳定状态和动态平衡状态下甚低频天线任意位置E的弧坐标，故E点位移为:

式中，U（S0，T）为E点从稳定状态到T时刻动态平衡状态的位移；R2为动态平衡下图中O点到E点在T时刻的距离；R1为稳定状态下O点到E点的距离；U1，U2，U3为U（S0，T）旋转坐标系下沿Ox轴、Oy轴和Oz轴方向的分量大小。

图1 拖曳天线回转示意图Fig.1 Rotary diagram of the trailing antenna

由于切向阻力很小［7］，模型中只考虑法向空气阻力。根据牛顿第二定律，天线的动力学方程为:

式中，P为天线拉力的数值大小；ρb为天线单位长度质量密度；ρa为空气密度；A0为天线的横截面积；CD为法向空气阻力系数；D为天线的直径；为法向速度。

由于机载甚低频拖曳天线末端连着一个锥袋，因此拖曳天线系统的固有边界条件为:

式中，W，M分别为锥袋重力和质量；CDD为锥袋的空气阻力系数；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 稳态动力学模型

在稳态时，即位移U=0，将U=0、式（1）和式（3）代入式（2）并无因次化得:

其中:

式中，a为载机盘旋半径。

稳态时，天线末端锥袋的位移UD=0，将其代入式（4）并无因次化得:

令 dr/ds=τ/p，则:

由式（5）可以得到:

上述无因次化模型的边界条件如下:

（1）在s=0处

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

式中，HP为载机盘旋高度；hp为无因次盘旋高度。

由式（7）可以得到:

3 基于Newton-Raphson迭代的边值问题求解算法

由上文可知，机载甚低频拖曳天线稳态动力学模型属于两点边值问题。目前已有许多学者对边值问题进行了研究。常见的主要有离散法和初值方法。离散法最直接的就是运用差分法［9］和变分法［10-11］进行求解；而初值方法目前最主要的就是打靶法［12］。本文利用 Newton-Raphson 迭代法［13］对打靶法试探解进行修正，该方法具有收敛速度快等优点。考虑如下一类两点边值问题:

式中，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 数值仿真

式（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，因此天线的活动区域属于对流层，其空气密度随高度变化为:

式中，ρ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可知，天线的垂直度与法向空气阻力系数相关，法向空气阻力系数越大，天线垂直度越小。法向空气阻力系数对天线末端的张力影响较小，但是对于天线与载机的连接点附近影响较大，法向空气阻力系数越大，天线张力越小。

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 结束语

本文以机载甚低频拖曳天线为研究对象，建立了机载甚低频拖曳天线动力学模型，并对其稳态动力学模型运用Newton-Raphson迭代进行求解。通过数值计算分析了法向空气阻力系数和天线锥袋质量对天线张力和垂直度的影响。本文的仿真结果对机载甚低频拖曳天线系统的设计有一定的指导作用。

［1］韩维，侯志强.机载拖曳系统的建模和若干动力学问题［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，Farsangi 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.