飞机典型模态特性仿真实验报告.doc

飞行力学实验报告飞机典型模态特性仿真 航空科学与工程学院飞行力学实验班课程实验飞机典型模态特性仿真实验报告ii目录一、实验的目的、内容11、实验目的12、实验内容13、实验步骤1二、理论模态及参数31、纵向模态的计算31.1计算方法31.2计算结果31.3计算结果分析42、纵向模态参数计算42.1计算方法42.2计算结果52.3计算结果分析53、横、航向模态的计算53.1计算方法53.2计算结果63.3计算结果分析64、横、航向模态参数计算64.1计算方法64.2计算结果64.3结果分析7三、非实时仿真实验81、纵向模态仿真81.1仿真方法81.2仿真结果图81.3作图法计算纵向模态参数101.4仿真纵向模态参数列表141.5仿真纵向模态参数与理论纵向模态参数的对比与误差分析142、横向模态仿真152.1仿真方法152.2仿真结果图152.3作图法计算横、航向模态参数172.4仿真横、航向模态参数列表192.5仿真横、航向模态参数与理论横、航向模态参数的对比与误差分析19四、修改静稳定性导数211、纵向模态211.1长周期模态211.2短周期模态222、横、航向模态252.1滚转模态252.2螺旋模态262.3荷兰滚模态28五、实时仿真实验301、升降舵受阶跃扰动302、副翼受阶跃扰动31六、放宽静稳定性32七、思考题33附录一 理论纵向模态计算程序35附录二 理论纵向模态参数计算程序36附录三 理论横、航向模态计算程序37附录四 理论横、航向模态计算程序38附录五 仿真纵向模态参数计算程序39附录六 仿真横、航向模态参数计算程序40图目录图 1 默认气动参数下的w变化图8图 2 受到扰动前50s内w的变化图（主要表现为短周期模态）9图 3 默认气动参数下的q变化图9图 4 默认气动参数下的变化图10图 5 默认气动参数下的u变化图10图 6 做出稳态值11图 7 取点计算短周期模态参数13图 8 默认参数下p的变化图15图 9 扰动初期p值的波动16图 10 默认参数下v的变化图16图 11 默认参数下的变化图17图 12 默认参数下r的变化图17图 13 升降舵受阶跃扰动后飞行状态图30图 14 升降舵受阶跃扰动后纵向运动参数变化图31图 15 副翼受到阶跃扰动后飞行状态图31表目录表 1 纵向模态特征值3表 2 纵向各阶模态3表 3 纵向模态参数5表 4 横向模态特征值6表 5 横向各阶模态6表 6 横向模态参数6表 7 长周期模态参数计算取点过程表12表 8 仿真与理论模态对比14表 9 荷兰滚模态参数计算取点过程表18表 10 Xu=-0.003时纵向运动参数变化图21表 11 Xu=0时纵向运动参数变化图22表 12 Zw=0.4时纵向运动参数变化图23表 13 Zw=0.43时纵向运动参数变化图24表 14 Lp=0时横、航向运动参数变化图25表 15 Lp=0.2时横、航向运动参数变化图26表 16 Nv=0.0014时横、航向运动参数变化图27表 17 Nv=0.002186时横、航向运动参数变化图27表 18 Yv=0时横、航向运动参数变化图28表 19 Yv=0.05时横、航向运动参数变化图295一、实验的目的、内容1、实验目的飞机运动模态是比较抽象的概念，是学习飞行力学过程中的重难点。本实验针对飞机运动模态的问题，通过计算机动态仿真和人-机飞行仿真实验平台上的驾驶员在环仿真实验，体会飞机响应与模态特性的关系，加深对飞机运动模态特性的理解。2、实验内容（1）纵向摸态特性实验计算某机在某状态下的短周期运动、长周期运动的模态参数；进行时域的非实时或实时仿真实验，操纵升降舵激发长、短周期运动模态，并由结果曲线分析比较模态参数；放宽飞机静稳定性，观察典型操纵响应曲线，并通过驾驶员在环实时仿真体验飞机的模态特性变化。（2）横航向模态特性实验计算某机在某状态下的滚转、荷兰滚、螺旋模态参数；进行时域仿真计算，操纵副翼或方向舵，激发滚转、荷兰滚等运动模态，并由结果曲线分析比较模态参数。3、实验步骤（1）模态特性分析。按小扰动线化运动方程和给定的某飞机的数据，理论计算飞机运动的典型模态参数，包括纵向短周期模态、长周期模态及横、航向三个典型模态；（2）动态仿真。对上述模态分析结果，分别进行纵向和横、航向的实时或非实时动态仿真，观察飞机的动态响应。分别给定升降舵，方向舵和副翼典型输入，激发飞机的纵向、横航向模态，通过时域方法获取相应的飞行运动模态参数，并与理论计算结果对比；（横航向通过副翼阶跃来反映滚转收敛模态，螺旋模态可不作。）尝试改变某些气动参数，观察飞机的模态特性变化。（3）人机闭环实验。通过修改参数放宽飞机静稳定性，进行人机闭环实时操纵，通过跟踪纵向俯仰指令体验飞机模态特性与飞行品质的关系。（4）整理撰写实验报告。二、理论模态及参数1、纵向模态的计算1.1计算方法由公式X=AX计算模态及特征值。其中A矩阵为 -0.006868 0.01395 -0.09055-0.3151 0-32.2000 773.98000 0.0001187 -0.00102600 -0.4285 0.0000010X为uwq用Matlab编程计算。在Matlab中输入A矩阵。用eig函数计算，得到A矩阵的特征值及特征向量。特征向量的列向量即为各个模态。程序见附录一。1.2计算结果表 1 纵向模态特征值1234-0.3719 + 0.8875i-0.3719 - 0.8875i-0.0033 + 0.0672i-0.0033 - 0.0672i表 2 纵向各阶模态A1A2A3A40.0211 + 0.0166i0.0211 - 0.0166i-0.9983-0.99830.99960.9996-0.0573 + 0.0097i-0.0573 - 0.0097i-0.0001 + 0.0011i-0.0001 - 0.0011i-0.0001 - 0.0000i-0.0001 + 0.0000i0.0011 - 0.0004i0.0011 + 0.0004i0.0001 + 0.0021i0.0001 - 0.0021i1.3计算结果分析从表1中看出，特征值为两对共轭复数。1,2实部、虚部远大于3,4。故1,2对应短周期模态，3,4对应长周期模态。相应的，由表看出A1和A2各项共轭，A3和A4各项共轭。其中A1,2的w较大，因此，短周期模态对应迎角的迅速变化；A2，3的u变化较大，因此，长周期模态表现为速度的变化。2、纵向模态参数计算2.1计算方法模态参数包括角频率，周期T，无阻尼自然频率n，阻尼比，半衰时间t1/2（或倍增时间t2），半衰时振荡次数N1/2（或倍增时振荡次数N2）。已知特征值=abi，w，T求解公式如下：=bT=2解方程组n1-2=n=|a|可得wn及。再根据公式t2=ln2wnN2=t2T得到t2及N2。用Matlab编程计算，程序见附录二。2.2计算结果表 3 纵向模态参数wTwnt1/2N1/2短周期模态 0.88757.07930.96230.38651.86360.2632长周期模态0.067293.45650.067280.04905210.71602.25472.3计算结果分析从表3中可以看出，长周期模态的周期远大于短周期模态周期。同时，长周期模态自然频率及阻尼比均较小。3、横、航向模态的计算3.1计算方法横、航向模态的计算方法和纵向模态是一样的，也用公式X=AX。此时，A矩阵为-0.0558 0.00000-774.000032.2 -0.0039 -0.43420.41360000.00 0.00110-0.0061 -0.1458000.000.0000001.000000.0000000.00X矩阵为vpr用Matlab编程计算，程序见附录三。3.2计算结果表 4 横向模态特征值1234-0.0330 + 0.9465i-0.0330 - 0.9465i-0.5625-0.0073表 5 横向各阶模态A1A2A3A4-1.0000-1.0000-0.99720.98210.0019 - 0.0032i0.0019 + 0.0032i-0.0367-0.0014-0.0001 + 0.0011i-0.0001 - 0.0011i0.00210.0078-0.0035 - 0.0019i-0.0035 + 0.0019i0.06520.18803.3计算结果分析从表4可以看出，横向模态有一对共轭复根1,2、一个大实根3和一个小实根4。共轭复根对应荷兰滚模态，大实根对应滚转模态，小实根对应螺旋模态。4、横、航向模态参数计算4.1计算方法横、航向模态参数计算方法同纵向模态参数。由于滚转模态及螺旋模态为一阶系统，因此，只计算了荷兰滚模态参数。Matlab程序见附录四。4.2计算结果表 6 横向模态参数wTwnt2N2荷兰滚模态0.94656.63800.94710.0348420.99723.16324.3结果分析从表6中可以看出，荷兰滚模态振荡较快，周期为6.6s。经过约21s后衰减为原来的一半。三、非实时仿真实验1、纵向模态仿真1.1仿真方法积分步长为0.01s，积分方法为欧拉法，记录参数为， u，w，q及时间t，各气动参数均为默认值。令升降舵受到幅值为1的阶跃扰动，进行非实时仿真。1.2仿真结果图仿真结果如下：图 1 默认气动参数下的w变化图纵向模态分为短周期模态和长周期模态。受到扰动后，短周期模态迅速收敛，长周期模态收敛缓慢。从图1中可以明显地看出短周期模态和长周期模态的收敛情况。受到扰动后，w迅速减小到-1240左右，然后又迅速增大，出现几次跳跃式的振荡。十几秒后，w开始表现为缓慢收敛的振荡，到1000s时，振荡幅值已经很小。刚受到扰动时，w主要表现为短周期模态。大约20s后，短周期模态收敛，w表现为长周期模态。图2是受扰动前50s内w变化图。图 2 受到扰动前50s内w的变化图（主要表现为短周期模态）图 3 默认气动参数下的q变化图图3是受阶跃扰动时，q的变化图。刚受到扰动时，q表现为短时间内的大幅振荡。几十秒后，短周期模态收敛，q开始表现为长周期模态。随着时间的增加，q收敛于0。图 4 默认气动参数下的变化图图4是的变化图。变化图中，短周期模态表现不明显。长周期模态收敛于1左右。图 5 默认气动参数下的u变化图图5中可以看出，受到扰动后，u变化非常大，从0振荡至5000左右。短周期模态表现不明显，长周期模态收敛于约2700。1.3作图法计算纵向模态参数由于u、变化图中，短周期模态表现不明显，只能选择w或q变化图计算模态参数，这里选用w。1.3.1长周期模态参数的计算由于w的振荡总体上表现为长周期模态，因此长周期模态有较多可选择的数据。首先根据图像做出w的稳态值。在图中，读出w的8个峰值及与峰值对应的时间。根据各峰值对应的时间可求出周期T，可由公式w=2T求出角频率。再将10个峰值与稳态值做差，得到10个幅值。根据8个幅值求出半衰时间t1/2。解方程组wn1-2=wt1/2=ln2wn得到wn及。最后，由公式N1/2=t1/2T求出半衰时振荡次数。作图过程如下：图 6 做出w变化图的稳态值表 7 长周期模态参数计算取点过程表1.3.2短周期模态的计算由于短周期模态收敛很快，可选择的数据很少。因此用超调量Mp和峰值时间tp求短周期模态参数。在图中做出峰值wmax，稳态值w及峰值时间tp。如下：图 7 取点计算短周期模态参数由公式Mp=wmax-ww求出超调量。由公式Mp=e-1-2求出。由tp=wn1-2求出wn。由w=wn1-2T=2wt12=ln2wnN12=t1/2T求出角频率，周期，半衰时间及半衰时间内的振荡次数。用Matlab编程求出纵向模态参数，程序可见附录五。1.4仿真纵向模态参数列表wTwnt1/2N1/2短周期模态 0.91866.84000.93950.20983.51750.5142长周期模态0.067393.32860.06740.0473217.50002.33051.5仿真纵向模态参数与理论纵向模态参数的对比与误差分析表 8 仿真与理论模态对比wTwnt1/2N1/2短周期模态 仿真0.91866.84000.93950.20983.51750.5142理论0.88757.07930.96230.38651.86360.2632误差3.5%3.4%2.3%45.7%88.7%95.3%长周期模态仿真0.067393.32860.06740.0473217.50002.3305理论0.067293.45650.067280.04905210.71602.2547误差0.1%0.1%0.1%4%3.2%3.5%从图中可以看出，仿真结果总体上是比较准确的。特别是长周期模态，误差均在5%以内。对于短周期模态，w，T，wn比较准确，但是， t1/2， N1/2误差较大，经分析，原因如下：（1）在取点过程中，所取的峰值是短周期坐标与长周期的坐标的叠加。由于扰动初期，主要表现为短周期模态，因此求解过程中，忽略了长周期模态的振动，把所取的峰值认为是短周期模态的峰值。这可能是主要的误差来源。（2）由于短周期模态迅速收敛，可取的实验点很少。计算过程中，只用了一个点，具有很大的随机性，容易产生误差。2、横向模态仿真2.1仿真方法积分步长为0.01s，积分方法为欧拉法，记录参数为侧滑角，滚转角速度p，侧滑角速度r，滚转角及时间t，各气动参数均为默认值。令副翼受到幅值为1的阶跃扰动，进行非实时仿真。2.2仿真结果图图 8 默认参数下p的变化图横向模态分为滚转模态、螺旋模态和荷兰滚模态。对于稳定飞机，滚转模态迅速收敛，荷兰滚模态接着收敛，螺旋模态收敛很慢。从图8可以看出，副翼受到一阶跃相应后，p马上发生变化。在扰动初期，由于存在滚转模态和荷兰滚模态，p值一直发生小幅波动。扰动后期，小幅波动消失，表现为一阶收敛的螺旋模态。700s左右，p收敛到0。图9是0-50s内，p变化的图像，主要表现为滚转模态和荷兰滚模态。图 9 扰动初期p值的波动图 10 默认参数下v的变化图从图10中看出，受到扰动后，v发生变化，总体表现为一阶收敛的螺旋模态。700s是，v收敛至-170左右。在扰动初期，由于滚转模态和荷兰滚模态，v有小幅波动。图 11 默认参数下的变化图从图11知，主要表现为螺旋模态，滚转模态和荷兰滚模态不明显。700s，收敛到30左右。图 12 默认参数下r的变化图从图12知，r主要表现为螺旋模态，滚转模态和荷兰滚模态不明显。700s，r收敛到30左右。2.3作图法计算横、航向模态参数由于螺旋模态和滚转模态为一阶系统。因此，只计算荷兰滚模态参数。由于、r变化图中，荷兰滚模态表现不明显，只能选择p或v变化图计算，这里选用p变化图。由于荷兰滚模态和螺旋模态相互叠加，且螺旋模态对荷兰滚模态有较大的影响，因此，和纵向模态参数存在一稳态值不同，这里还需要选择相应的基准值。这里相邻的两个极值pmax和pmin求平均作为基准值p。再用（pmax-p）和（pmin-p）求出幅值。在图中做出10个峰值及其对应的时间，可以求出10个幅值。通过各个峰值对应时间求得周期T，通过幅值求出半衰时间t1/2。其他参数的求解同长周期模态。表 9 荷兰滚模态参数计算取点过程表取点后，用Matlab计算出荷兰滚模态参数，程序见附录六。2.4仿真横、航向模态参数列表wTwnt1/2N1/2荷兰滚模态 0.93936.68890.93970.02780.93933.96182.5仿真横、航向模态参数与理论横、航向模态参数的对比与误差分析wTwnt1/2N1/2荷兰滚模态 仿真0.93936.68890.93970.027826.50003.9618理论0.94656.63800.94710.0348420.99723.1632误差0.1%0.1%0.1%20.2%26.2%25.3%从图中可以看出，仿真结果总体上是比较准确的。特别是长周期模态，误差均在5%以内。对于短周期模态，w，T，wn比较准确，但是， t1/2， N1/2误差较大，经分析，原因如下：（1）在计算过程中，把两个极值的平均值作为基准，这是主要的误差来源。由于螺旋模态的影响，p变化图呈上升趋势，因此，相邻两个极值的基准的大小并相同，t坐标较大的极值，基准值也较大。取平均值作为基准的方法使幅值偏大，导致了阻尼比偏小，也对t1/2和 N1/2造成了影响。（2）在计算过程中，T是根据极值点所对应的时间得到了，因此，计算结果较为准确。w由T决定，而wn主要由w决定，因此，w和wn均较为准确。（3）取点不准确也会导致误差。四、修改静稳定性导数1、纵向模态1.1长周期模态长周期模态受w和q影响较小，因此忽略w和q。长周期模态方程变为u=Xu-g00u可见，长周期模态的稳定性主要受Xu影响。由特征方程2-Xu=0知，当Xu0时，长周期模态发散。Xu默认值为-0.006868，增大Xu可令长周期模态发散。表 10 Xu=-0.003时纵向运动参数变化图从表10的各个图中可以看出，当Xu=-0.003时，各运动参数的长周期模态振荡幅值大大增加，收敛速度减缓。当t=1000s时，各运动参数仍存在较大的振荡。但短周期模态没有受到太大影响。表 11 Xu=0时纵向运动参数变化图从表11中可以看出，当Xu=0时，各运动参数已经发散。理论上，Xu=0各参数应处于临界稳定状态，各运动参数呈等幅振荡。但由于理论分析时，忽略了w和q，因此精度受限。1.2短周期模态由于短周期模态受u和的影响较小，忽略u和。短周期模态方程变为wq=ZwZqMwMqwq特征方程为2-Zw+Mq+ZwMq-ZqMw=0因此，短周期模态稳定的条件为-Zw+Mq0ZwMq-ZqMw0可以调整Zw，Mq，Zq，Mw的值使其不满足短周期模态稳定的条件。这里改变Zw的值。Zw的默认值是-0.3151，Mq的值是-0.4285，当Zw0.4285时，短周期模态发散。表 12 Zw=0.4时纵向运动参数变化图（备注：由于q随时间变化图振荡太多，这里将t改为0-300显示）（备注：由于q随时间变化图振荡太多，这里将t改为0-300显示）从表12可以看出，当Zw=0.4时，短周期模态振荡幅度明显增大，收敛速度十分缓慢。但是，长周期模态基本不受影响，收敛速度不变。另外，u和的变化和Zw=-0.3151时几乎一样，只是受扰动初期有小幅波动，可见，短周期模态受u和的影响较小。表 13 Zw=0.43时纵向运动参数变化图从表13中，可以看出，Zw=0.43时短周期模态已经发散。但由u变化图中可以看出，此时，长周期模态仍然是收敛的。2、横、航向模态2.1滚转模态滚转模态可以近似为p=Lpp当Lp0只要改变Lv，Nr，Lr，Nv的值即可使螺旋模态发散。Lv的默认值为-0.003865，Nr的默认值为-0.1458，Lr的默认值为0.4136，Nv的默认值为0.001086。这里逐渐增大Nv使其不满足稳定性条件。表 16 Nv=0.0014时横、航向运动参数变化图从表16的各图中可以看出，当Nv为0.0014时，螺旋模态已经发散。当滚转模态和荷兰滚模态仍是收敛的。表 17 Nv=0.002186时横、航向运动参数变化图从表17的各图中可以看出，当Nv为0.002186时，螺旋模态发散剧烈。2.3荷兰滚模态荷兰滚模态方程可近似简化为vr=YvYrNvNrvr因此，特征方程为2-Yv+Nr+YvNr-YrNv=0受扰动后稳定的条件为-Yv+Nr0YvNr-YrNv0对气动参数做相应调整使其不满足上面不等式可以使荷兰滚发散。这里，增大Yv。表 18 Yv=0时横、航向运动参数变化图从表18各图中可见，当Yv=0时，荷兰滚模态收敛速度大大减小，振幅增大。对于p和v，当t=1000s时，荷兰滚模态仍有较大的波动，但其收敛趋势已经明显。从图中可知，增大Yv，螺旋模态没有受到影响。表 19 Yv=0.05时横、航向运动参数变化图从表19各图中可见，当Yv=0.05时，荷兰滚已经剧烈发散。五、实时仿真实验1、升降舵受阶跃扰动图 13 升降舵受阶跃扰动后飞行状态图图 14 升降舵受阶跃扰动后纵向运动参数变化图从图13中可以看到，当飞机升降舵受到一幅值为1的阶跃扰动后，飞机低头，迎角不断减小，直至坠毁。从图14中看出，受扰动后，u，w，q均减小。2s后，q有上升趋势，但飞机马上坠毁。2、副翼受阶跃扰动图 15 副翼受到阶跃扰动后飞行状态图从图15中可以看出，当飞机副翼受到一幅值为1的阶跃扰动后，飞机发生滚转，主要表现为滚转模态。一段时间后滚转减缓，滚转模态收敛，飞机主要表现为螺旋模态。173s后，飞机坠毁。六、放宽静稳定性放宽静稳定性，可以改变重心位置。理论上，放宽静稳定性，飞机的不稳定性增大，运动参数曲线振荡增大。实验中，将重心h的位置多次进行改动，但画出曲线与h=0时完全一样，与理论不符。实验时，还发现，改变Xu的值，实验曲线没有任何改变。但理论上，Xu=Xum，改变Xu，会导致Xu改变，必然会对稳定性产生影响。可以推测，该实验软件并未将Xu和Xu相关联，因此，稳定性没有受到影响。根据上述两种现象，可以推测，该软件在编程时，把A矩阵的各个参数设置为定值，没有与其他参数关联。只要没有改变A矩阵中各个参数的具体数值，实验结果是不会发生变化的。七、思考题1、哪些气动参数对纵向运动模态特性影响较大，试简单分析其原因；答：（1）长周期模态的稳定性主要受Xu影响。因为长周期模态受w和q影响较小，因此忽略w和q。长周期模态方程变为u=Xu-g00u可见，长周期模态的稳定性主要受Xu影响。由特征方程2-Xu=0知，当Xu0时，长周期模态发散。（2）短周期模态的稳定性主要受Zw，Mq，Zq，Mw影响。由于短周期模态受u和的影响较小，忽略u和。短周期模态方程变为wq=ZwZqMwMqwq可见，短周期模态的稳定性主要受Zw，Mq，Zq，Mw影响。由特征方程知2-Zw+Mq+ZwMq-ZqMw=0因此，短周期模态稳定的条件为-Zw+Mq0ZwMq-ZqMw0改变Zw，Mq，Zq，Mw可以使短周期模态发散。2、纵向气动舵面阶跃输入会引起飞机哪些状态参数的变化，典型的变化趋势怎样？答：纵向气动舵面的输入会引起， u，w，q等的变化。典型的变化趋势是：当受到一阶跃扰动时，飞机的状态参数， u，w，q等首先发生一个较大的改变。接着产生类似于简谐波的振荡。如果飞机使稳定的，振荡是衰减的，随着时间推移，振荡将会收敛；如果飞机是不稳定的，振荡将会发散。3、简要说明飞机线化方程的使用局限性。答：利用“小扰动”假设将飞机运动方程线化，使得解析法求解微分方程成为可能。但是，这样处理有一定的局限性。首先，进行线化处理后，求解的精度会受到限制；其次，线化方程基于“小扰动”假设，当扰动较大时，假设不成立，不能对运动方程线化。附录一 理论纵向模态计算程序function a lamda x0 = theoreticalmode()%计算a矩阵，通过a矩阵计算理论模态特征值a=zeros(4,4);x0u=-0.006868;x0w=0.01395;z0u=-0.09055;z0q=773.98;z0w=-0.3151;m0u=0.0001187;m0q=-0.4285;m0w=-0.001026;g=32.2;a(1,1)=x0u;a(1,2)=x0w;a(1,3)=0;a(1,4)=-g;a(2,1)=z0u;a(2,2)=z0w;a(2,3)=z0q;a(2,4)=0;a(3,1)=m0u;a(3,2)=m0w;a(3,3)=m0q;a(3,4)=0;a(4,1)=0;a(4,2)=0;a(4,3)=1;a(4,4)=0;x0,lamda=eig(a);End附录二 理论纵向模态参数计算程序function ws Ts ksis wns t2s n2s wl Tl ksil wnl t2l n2l = modeparameter(lamda)%短周期模态参数（用s（short首字母）后缀表示）ws=abs(imag(lamda(1,1);Ts=2*pi/ws;syms ksis wns;ksis,wns=solve(ksis*wns=0.3719,wns*sqrt(1-ksis2)=0.8875);ksis=vpa(ksis,4);wns=vpa(wns,4);t2s=log(2)/abs(real(lamda(1,1);n2s=t2s/Ts;%长周期模态参数（用l（long首字母）表示）wl=abs(imag(lamda(3,3);Tl=2*pi/wl;syms ksil wnsl;ksil,wnl=solve(ksil*wnl=0.0033,wnl*sqrt(1-ksil2)= 0.0672);ksil=vpa(ksil,4);wnl=vpa(wnl,4);t2l=log(2)/abs(real(lamda(3,3);n2l=t2l/Tl;end附录三 理论横、航向模态计算程序function a_l lamda_l x0_l =lateraltheoreticalmode()a_l=zeros(4,4);g=32.2;y0p=0;y0r=-774;y0v=-0.0558;l0p=-0.4342;l0r=0.4136;l0v=-0.003865;n0p=-0.006112;n0r=-0.1458;n0v=0.001086;a_l(1,1)=y0v;a_l(1,2)=y0p;a_l(1,3)=y0r;a_l(1,4)=g;a_l(2,1)=l0v;a_l(2,2)=l0p;a_l(2,3)=l0r;a_l(2,4)=0;a_l(3,1)=n0v;a_l(3,2)=n0p;a_l(3,3)=n0r;a_l(3,4)=0;a_l(4,1)=0;a_l(4,2)=1;a_l(4,3)=0;a_l(4,4)=0;x0_l,lamda_l=eig(a_l);end附录四 理论横、航向模态计算程序function w_n T_n ksi_n wn_n t2_n n2_n = lateralmodeparameter(lamda_l)%荷兰滚模态参数（用n(荷兰Netherlands首字母)表示）w_n=abs(imag(lamda_l(1,1);T_n=2*pi/w_n;syms ksi_n wn_n;ksi_n,wn_n=solve(ksi_n*wn_n=0.0330,wn_n*sqrt(1-ksi_n2)=0.9465);ksi_n=vpa(ksi_n,4);wn_n=vpa(wn_n,4);t2_n=log(2)/abs(real(lamda_l(1,1);n2_n=t2_n/T_n;end附录五 仿真纵向模态参数计算程序function ksi_l wn_l w_l T_l t2_l n2_l ksi_s wn_s w_s T_s t2_s n2_s = simulatemodeparameter( )%通过w图像的周期T、半衰时间t2等计算长周期模态参数y_infinite=-821;y_max=-688.5 -724.3 -748.0 -767.2 -7