气动模型与简单飞行仿真的Simulink（MATLAB）实现

罗天宇

一个飞机爱好者又双来搞点事情……
为便于理解，请在确保理解前文的基础上阅读本文。转载请注明出处。
感谢@优雅的刺猬 一如既往地支持~
上次我们说到，关于气动模型建立，就是要建立一种数据查找机制，让六自由度方程计算过程中实时能够按需拿到对应角度的气动力和力矩的值。这里用到的是 Simulink 中的 n-D Lookup Table。前文各矩阵录入到工作空间，然后再去Simulink中的n-D Lookup Table，就可以直接填入矩阵名称了。
以 \Delta C_{\square}(\alpha,\beta,\delta e) 为例，一个典型的数据查表输入输出流如下图所示：



数据查表流程结构

这里我为了方便确认输出数据是否有误，把气动数据的升力系数和阻力系数也放了进去，以便在出现意外情况的时候判断数据是否合理，从而找出数据问题。如果建表的时候用的是自然下标那么记得把输入的角度转换成下标对应值。
将 \Delta C_{\square}(\alpha,\beta) 、 \Delta C_{\square}(\alpha,\beta,\delta a) 、\Delta C_{\square}(\alpha,\beta,\delta e) 、 \Delta C_{\square}(\alpha,\beta,\delta r) 均按照上图做好，并进行封装后，整个基于数据的气动力模型就在simulink中建立好了，如下图所示：



基于数据的气动力模型

（这里之所以说是基于数据的气动力模型，是因为还有另外一种基于公式的气动力模型，在前文提到过，大家还记得吗？）
这里的输入是 \alpha 和 \beta 以及三组舵面偏角，输出是飞机体轴系六分量力和力矩，外加升力系数和阻力系数（这两个也可以不要）。
到这里为止，simulink中的气动力模型就建立完毕了。
<hr/>接下来进入仿真环节。这一部分目前我只搭起了一个简单仿真，比较粗糙，但是已经包含了各核心内容。先上一张全图。


是的，就是大标题上方那张。接下来我会逐一说明这里面的关键模块。
还是先从源头——气动力模型开始。这里将本文开始所述的气动力模型又封装了一层。成为下图这样：


可以看到输入输出和前文是一样的。这里左边的三个舵偏用一个名为fcs_cmd的数组导入，输出的体轴系六分量力和力矩通过Mux合成一个数组。我们就以这里为起点，沿着流程继续向下，到了这里：


这个模块封装的是从系数到力和力矩的转换计算，就是如下公式：




这里一定要注意，力矩系数和力矩相互转换时，参考长度 l 针对绕不同的轴的旋转所取的量是不一样的。对于绕 Y 轴的俯仰运动， l 为平均气动弦长；对于绕 X 轴的滚转和绕 Z 轴的偏航， l 为机翼展长。我前期处理数据的时候这地方犯错导致全部重新处理了一遍。
这个地方我没有用官方的aerodynamic forces and moments模块是因为我气动力模型的数据是以质心为参考点的，没有压心或焦点与质心的力臂，所以简化了。
现在我们有了力和力矩，接下来就要把他们送入动力学方程了。这里用的是最简单的六自由度模型：


固定质量的欧拉角描述的六自由度模型，里面同时包含了动力学和运动学模型。我认为这里要明白这个模型的原理是什么才算懂了，不能只知道拖个模型过来，所以要多讲一下里面包含的内容。参数设置官方帮助已有，不再赘述。我们CTRL+U进去看看里面是什么。


这张图从中间一分为二上下两部分看，上半部分是平动部分，就是力、质量、加速度、速度、位置之间的关系。下半部分是力矩、惯量、角加速度、角速度、欧拉角之间的关系。已经写明的数学运算不再赘述。我们看一下图片中间靠下这个封装模块。输入是pqr，惯量和惯量变化率，力矩。其内部是这样的：


由于我们这次搭的模型中不牵扯到质量和惯量变化，所以选用的这个六自由度模块惯量变化率是保持零。这一部分是矩阵运算形式，对应内容为文章开头链接第一篇文章中的P10~P12。感兴趣的同学可以推导一下看是否一致。
回到最上层，现在我们有了飞机位置、飞机姿态欧拉角、飞机体轴系下的速度矢量、地坐标系到体坐标系的转换矩阵（图中的DCMbe）。这里要用这个矩阵把重力加上，同时把推力直接加上作为输入，如下图这一部分。


同理，推力的力矩也可以加到下面力矩上一起作为输入。
飞机体轴系下的速度引出来送到这里：


这个封装的模块用来求速度标量、气流偏角。公式比较简单不再写了。用官方的Incidence, Sideslip, & Airspeed模块里面那个asin除零会报错，要换成atan2，所以自己写了。速度标量传到空速表可以实时显示。
飞机姿态欧拉角送到这里：


解算出的角度送到相应的地平仪和航向仪中。
本次模型搭建没有考虑地球曲率，用的大气模型为：


注意北东地坐标系高度为Z坐标乘以-1。
当然，如果想实时显示的话可以把三个位置量和三个姿态角传送到UE4（虚幻4游戏引擎）中进行显示：


这个展开的话又是另外一个故事了。
至此整个模型搭建完毕。先到这里，大家还想看哪个模块的讲解可以留言，我会不定时对文章进行补充。
后续可以尝试给飞机加个飞控系统，我打算搞一搞看看能不能搞通……（抱歉，我又挖坑了……）
未完待续……
出处：
https://zhuanlan.zhihu.com/p/428517085
https://zhuanlan.zhihu.com/p/410423919
作者：机智的刺猬

mfy0508

能把模型分享一下吗大神？[大笑]你太优秀了，想学习一下！

雾满岚江

不敢当。学习的话还是请按照步骤走一遍吧，debug也是个很重要的学习过程呢。

流云

非常棒的学习资料[赞]

Pegasus8

感谢~开心~

freejaja

为啥不自己搭模型呀 这个是航空航天工具箱里的六自由度模型吧[大笑]

a840177949

牛掰[赞][赞][赞]，得找你好好指导我一下，我最近在头疼simulink搭建飞机动力学模型

甜心萝卜糕

打算先用matlab自带的跑通，我前面文章提到过，后面计划在UE4里面自己搭动力学方程

a840177949

这个跟着做就能搭起来一套哈，可以看看前几篇文章，一步一步过来就行了

Pegasus8

师兄，有用JSBsim搭模型吗？

甜心萝卜糕

没有，看到过这个东西，还没搞懂怎么用

shunshui

请问下计算动压rho给的速度是不是写错了？应该是V而不是Vb？

yangwei850622

这个地方的Vb是矢量，是速度在体轴系的三分量，送进去后在里面求了矢量大小的

shunshui

这里的仿真运行结果是什么样的啊，是不是会按照设置的高度和转速运行啊？

shunshui

是根据舵偏角和推力输入可以看到仪表盘飞机飞行参数变化，也就是模拟飞行驾驶了

show712416

你好，想请问一下，搭建出飞机的动力学模型，然后通过simulink自带的linear analysis tool对飞机的工作点做了线性化，请问怎么根据这个线性化的结果设计相应的pid参数啊？

A...周东

我还没有用过这个工具[大哭]

御野紘

好吧[飙泪笑]

mfy0508

你好，最近在做飞行器仿真，想问一下您在Fxyz输入处，g乘了一个7.5的gain是怎么来的。

PLA

有可能是我一开始做的时候发现飞机受重力效果不明显增加了重力加速度测试用的。我目前装有这个工程的电脑送修了，等回头电脑回来我看一下