基于神经网络模型的舰面流场仿真算法

摘要

飞行甲板表面流场是影响直升机飞行安全的重要因素，在直升机实时飞行仿真中，流场速度分布特性直接影响直升机飞行动力学仿真的精确度。文章以纳维-斯托克斯方程为基础，利用BP神经网络算法，研究了实时确定流场速度分布的方法，该方法可用于甲板流场实时仿真，提高直升机飞行仿真的精度。

正文

《1 舰面流场的特点》

1 舰面流场的特点

非航空型舰船的飞行甲板, 一般位于舰船尾部。飞行甲板前方紧靠着钝体结构的机库。在舰船航行过程中, 前方相对气流掠过机库上方进入甲板区域。由于机库的外形为钝体结构, 气流掠过机库后, 产生陡壁体效应, 使稳定的气流遭到破坏, 并在机库后方和甲板上方产生下冲, 在甲板不同位置分别形成无风区、涡流区、稳定流动区等。这种紊乱的气流给机体附加了一个随时间变化的气动载荷^[1], 影响到机体有关结构的寿命, 增加了保持直升机航迹和姿态角的难度。直升机骤然进入该区域, 驾驶员感觉直升机突然被“吸”向机库, 容易造成撞机事故。

《2 舰面流场研究方法》

2 舰面流场研究方法

研究直升机的舰上起降特性, 必须研究甲板流场的流动特性。目前一般采用风洞试验和计算的方法确定流场速度和压力分布。风洞试验费用较高, 并且测量点有限, 对于紊乱流动区域难以测得准确结果。流场计算一般采用纳维-斯托克斯 (N-S) 方程。计算求解甲板流动特性虽然费用较低, 但是在计算过程中首先要划分计算网格, 然后进行有限元迭代计算, 因而计算耗时较长, 而且即使网格划分得再密, 最终得到的也是离散的速度场分布。进行直升机舰上飞行仿真时, 直升机的飞行轨迹是一条连续曲线, 这就需要实时得到甲板上方任意点的流场特性。因此单独依靠N-S方程并不能满足实时仿真要求。

神经网络技术的发展为解决这一问题提供了有效工具。神经网络以大规模模拟并行处理为主, 具有很强的鲁棒性和容错性, 善于联想、概括、类比和推广, 任何局部的损伤都不会影响整体结果, 并且具有很强的自学习能力, 适于求解流场类的非线性问题。

为了满足在实时仿真中直升机飞行特性计算的需要, 本文尝试以N-S方程为基础, 利用神经网络的误差反向传播算法, 仿真计算甲板流场任一点的流场特性。

《3 纳维-斯托克斯方程的有限元解法》

3 纳维-斯托克斯方程的有限元解法

在直升机着舰过程, 主要是机库上方的下冲气流影响直升机旋翼的气动力分布, 因此计算甲板对称平面内的下冲流场, 使用二维N-S方程即可^[2]。二维定常粘性不可压N-S方程为

$\begin{array}{l} u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \frac{l_{x_{1}}}{F_{R}^{2}} - \frac{\partial p}{\partial x} + \frac{1}{R_{e}} (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}), (1) \\ u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = \frac{l_{x_{2}}}{F_{R}^{2}} - \frac{\partial p}{\partial y} + \frac{1}{R_{e}} (\frac{\partial^{2} v}{\partial x^{2}} + \frac{\partial^{2} v}{\partial y^{2}}), (2) \\ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 。 (3) \end{array}$ $\begin{array}{l} u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \frac{l_{x_{1}}}{F_{R}^{2}} - \frac{\partial p}{\partial x} + \frac{1}{R_{e}} (\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}}), (1) \\ u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = \frac{l_{x_{2}}}{F_{R}^{2}} - \frac{\partial p}{\partial y} + \frac{1}{R_{e}} (\frac{\partial^{2} v}{\partial x^{2}} + \frac{\partial^{2} v}{\partial y^{2}}), (2) \\ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 。 (3) \end{array}$

上述方程组中, 式 (1) 和式 (2) 为动量方程, 式 (3) 为连续方程。其中R_e为雷诺数, F_R为弗洛德数, l_x₁=G_x/g, l_x₂=G_y/g为重力场沿x, y方向的余弦。 u为沿x轴方向的气流速度, v为沿y轴方向的气流速度, p为流场内的气流压力。

式 (1) 和式 (2) 用伽辽金加权余数法进行加权处理, 并将计算区域分为n个, 在每个区域里引进型函数N描述速度及压力。使用格林公式给其中的二阶项降解后, 即得x方向离散表达式:

$\begin{array}{l} \sum_{1}^{n^{e}} {\int_{A^{e}} [Ν_{i} Ν_{k} u_{k} \frac{\partial Ν_{j}}{\partial x} u_{j} + Ν_{i} Ν_{k} v_{k} \frac{\partial Ν_{j}}{\partial y} u_{j} + Ν_{i} \frac{\partial Μ_{l}}{\partial x} p_{l} - \\ Ν_{i} \frac{l_{x_{1}}}{F_{R}^{2}} + \frac{1}{R_{e}} (\frac{\partial Ν_{i}}{\partial x} \frac{\partial Ν_{j}}{\partial x} u_{j} + \frac{\partial Ν_{i}}{\partial y} \frac{\partial Ν_{j}}{\partial y} u_{j})] d A^{e} - \\ \int_{Γ_{1}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial Ν_{j}}{\partial n} u_{j} d Γ - \int_{Γ_{2}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial u_{j}}{\partial n} d Γ} = 0 。 (4) \end{array}$ $\begin{array}{l} \sum_{1}^{n^{e}} {\int_{A^{e}} [Ν_{i} Ν_{k} u_{k} \frac{\partial Ν_{j}}{\partial x} u_{j} + Ν_{i} Ν_{k} v_{k} \frac{\partial Ν_{j}}{\partial y} u_{j} + Ν_{i} \frac{\partial Μ_{l}}{\partial x} p_{l} - \\ Ν_{i} \frac{l_{x_{1}}}{F_{R}^{2}} + \frac{1}{R_{e}} (\frac{\partial Ν_{i}}{\partial x} \frac{\partial Ν_{j}}{\partial x} u_{j} + \frac{\partial Ν_{i}}{\partial y} \frac{\partial Ν_{j}}{\partial y} u_{j})] d A^{e} - \\ \int_{Γ_{1}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial Ν_{j}}{\partial n} u_{j} d Γ - \int_{Γ_{2}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial u_{j}}{\partial n} d Γ} = 0 。 (4) \end{array}$

式中l=1, …, 4;j=1, …, 8; k=1, …, 8; i=1, …, n_POIN (节点总数) ;Γ $_{1}^{e}$ $_{1}^{e}$ 为元素的本质边界, 在该边界上, 给定边界点的变量值;Γ $_{2}^{e}$ $_{2}^{e}$ 为元素的自然边界, 在该边界上, 给定变量的一阶导数值。

同理, 可得y方向的离散表达式:

$\begin{array}{l} \sum_{1}^{n^{e}} {\int_{A^{e}} [Ν_{i} Ν_{k} u_{k} \frac{\partial Ν_{j}}{\partial x} v_{j} + Ν_{i} Ν_{k} v_{k} \frac{\partial Ν_{j}}{\partial y} v_{j} + Ν_{i} \frac{\partial Μ_{l}}{\partial y} p_{l} - \\ Ν_{i} \frac{l_{x_{2}}}{F_{R}^{2}} + \frac{1}{R_{e}} (\frac{\partial Ν_{i}}{\partial x} \frac{\partial Ν_{j}}{\partial x} v_{j} + \frac{\partial Ν_{i}}{\partial y} \frac{\partial Ν_{j}}{\partial y} v_{j})] d A^{e} - \\ \int_{Γ_{1}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial Ν_{j}}{\partial n} v_{j} d Γ - \int_{Γ_{2}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial v_{j}}{\partial n} d Γ} = 0 。 (5) \end{array}$ $\begin{array}{l} \sum_{1}^{n^{e}} {\int_{A^{e}} [Ν_{i} Ν_{k} u_{k} \frac{\partial Ν_{j}}{\partial x} v_{j} + Ν_{i} Ν_{k} v_{k} \frac{\partial Ν_{j}}{\partial y} v_{j} + Ν_{i} \frac{\partial Μ_{l}}{\partial y} p_{l} - \\ Ν_{i} \frac{l_{x_{2}}}{F_{R}^{2}} + \frac{1}{R_{e}} (\frac{\partial Ν_{i}}{\partial x} \frac{\partial Ν_{j}}{\partial x} v_{j} + \frac{\partial Ν_{i}}{\partial y} \frac{\partial Ν_{j}}{\partial y} v_{j})] d A^{e} - \\ \int_{Γ_{1}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial Ν_{j}}{\partial n} v_{j} d Γ - \int_{Γ_{2}^{e}} \frac{1}{R_{e}} Ν_{i} \frac{\partial v_{j}}{\partial n} d Γ} = 0 。 (5) \end{array}$

对连续方程使用压力型函数进行加权:

$\sum_{1}^{n^{e}} \int_{A^{e}} Μ_{i} (\frac{\partial Ν_{j}}{\partial x} u_{j} + \frac{Ν_{j}}{\partial y} v_{j}) d A^{e} = 0 。 (6)$ $\sum_{1}^{n^{e}} \int_{A^{e}} Μ_{i} (\frac{\partial Ν_{j}}{\partial x} u_{j} + \frac{Ν_{j}}{\partial y} v_{j}) d A^{e} = 0 。 (6)$

对计算区域内的每个节点, 均可离散出以上三个表达式, 将全部方程装配为矩阵形式:

$A λ = F + B, (7)$ $A λ = F + B, (7)$

其中

$(λ_{i i}) = [\begin{array}{l} (u_{i}) \\ (p_{i}) \\ (v_{i}) \end{array}]$ $(λ_{i i}) = [\begin{array}{l} (u_{i}) \\ (p_{i}) \\ (v_{i}) \end{array}]$

, A, F, B为相关系数, (u_i) , (p_i) , (v_i) 为网格节点的水平速度分量、压力值和垂直速度分量。

把流场划分为若干网格, 用有限元法迭代计算每个网格节点的速度分量和压力, 即可得到甲板表面的速度场分布框架。以该速度场为基础, 运用神经网络算法, 即可得到甲板表面任一点的速度分量。

《4 用于流场仿真的神经网络算法》

4 用于流场仿真的神经网络算法

BP算法是一种很有效的神经网络算法^[3], 如图1所示, 一般采用三层神经网络结构:输入层、隐层、输出层。BP模型把训练样本的输入输出问题转变为一个非线性优化问题, 使用优化算法中的梯度下降法, 通过迭代运算求解问题的最优解。通过加入隐层神经元使优化问题的可调参数增加, 从而得到更精确的解。

《图1》

图1神经网络模型

Fig.1 Model of neural network

BP算法的基本求解思路为:首先使用已知输入向量和输出向量组对网络进行训练, 训练过程中通过不断调整网络连接权值, 使网络的实际输出和期望输出一致, 然后固定此时的权值, 在网络输入层加上待求解的输入向量, 即可在输出层得到计算结果。

使输入层的输入模式为X_p, 权值为w, 各神经元阈值用θ表示, 则基本训练过程为:

步1 置各权重和阈值的初始值w_ij (0) , θ_ij (0) 为小的随机数值。

步2 提供训练样本:输入向量X_p (p = 1, 2, 3, …, P) 和目标输出d_p (p =1, 2, 3, …, P) , 对每个p进行步3～步5操作。

步3 计算网络的实际输出及隐层单元状态:

$o_{p j} = δ_{j} ((Ν e t)_{p j}) = f_{j} (\sum w_{j i} o_{i} + θ_{j}) 。$ $o_{p j} = δ_{j} ((Ν e t)_{p j}) = f_{j} (\sum w_{j i} o_{i} + θ_{j}) 。$

步4 计算训练误差:

$\begin{array}{l} δ_{p j} = o_{p j} (1 - o_{p j}) (t_{p j} - o_{p j}), (输出层) \\ δ_{p j} = o_{p j} (1 - o_{p j}) \sum δ_{k} w_{j k} 。 (隐含层) \end{array}$ $\begin{array}{l} δ_{p j} = o_{p j} (1 - o_{p j}) (t_{p j} - o_{p j}), (输出层) \\ δ_{p j} = o_{p j} (1 - o_{p j}) \sum δ_{k} w_{j k} 。 (隐含层) \end{array}$

步5 修正权值和阈值:

$\begin{array}{l} θ_{j} (t + 1) = θ_{j} (t) + η δ_{j} + ξ (θ_{j} (t) - θ_{j} (t - 1)) ‚ \\ w_{j i} (t + 1) = w_{j i} (t) + η δ_{j} o_{p i} + \\ α (w_{j i} (t) - w_{j i} (t - 1)) 。 \end{array}$ $\begin{array}{l} θ_{j} (t + 1) = θ_{j} (t) + η δ_{j} + ξ (θ_{j} (t) - θ_{j} (t - 1)) ‚ \\ w_{j i} (t + 1) = w_{j i} (t) + η δ_{j} o_{p i} + \\ α (w_{j i} (t) - w_{j i} (t - 1)) 。 \end{array}$

式中η为学习速率;ξ为常数, 其值的大小决定了过去权值的变化对目前权值变化的影响程度;α为动量项系数。

步6 判断权值指标是否满足精度要求, 若满足, 则转到步7;若不满足, 则转到步2, 继续循环。对每个p均进行上述操作。

步7 停止。

以流场网格节点的坐标值为输入模式, 以纳维-斯托克斯方程计算得到的网格节点的速度值作为目标输出, 对上述神经网络进行训练, 然后输入甲板任一点的坐标值, 即可得到该点的速度值。

《5 算例及分析》

5 算例及分析

《5.1N-S方程计算结果讨论》

5.1N-S方程计算结果讨论

依据上述原理, 对某舰的飞行甲板下冲流场进行了仿真计算。计算区域沿甲板纵向对称面为一长方形区域, 坐标原点在机库大门顶端, y方向向上为正, 即y > 0的区域为来流入口区域; x沿舰船纵向轴线向后为正。区域高度及长度均用无量纲的机库高度和长度表示。下边界为飞行甲板表面, 前边界为机库门所在的铅垂面, 计算高度取7个机库高度, 以此高度的气流流动为自然流动, 即不受机库存在的影响;向后长度取22个机库长度。在这样的区域内计算, 可以认为甲板处的气流流动不受计算边界变动的影响, 计算区域网格是不均匀划分的。

由于神经网络进行仿真时, 要使用N-S方程的计算结果作为样本对神经网络进行训练, 因此N-S方程计算结果的准确性对最后的仿真结果具有重要影响。为此, 将N-S方程的计算结果与试验结果进行了对比。该涡流区的长度, N-S计算结果与试验结果的比较如表1所示, 两者符合较好。由于着舰格栅中心到机库大门的距离为2.31, 因此, 当气流从舰船正前方吹来时, 直升机起降区域正处于紊流区内。

表1涡流区范围

Table 1 Range of vortex region

《表1》

比较项	涡流区长度 (无量纲值)	与计算值相对偏差/%
计算	3.01
风洞试验	2.88	5.2
烟风洞试验	2.83	7.1

在计算过程中, 着舰位置上方几个特定点下冲气流的垂直速度值与试验值的比较如表2所示, 其中5.4 m高度的值差别比较大, 是由于考虑到网格划分的因素而用5.2 m处的计算点值近似代替5.4 m处的计算值造成的。

表2着舰点上方下冲气流垂直速度值比较 (来流速度为18 m/s)

Table 2 Compare of vertical speed above the landing point (free-flow velocity =18 m/s)

《表2》

距甲板高度 / m	试验结果 /m·s^-1	计算结果 /m·s^-1	相对误差 / %
2.4	2.43	2.25	7.4
5.4	2.98	3.35	12.4
7.9	3.17	2.95	6.9

从上述比较可以看出, 用N-S方程计算得到的结果具有较高的精度, 因而用其作为神经网络模型的训练样本是可行的。

《5.2使用神经网络进行流场仿真》

5.2使用神经网络进行流场仿真

首先用N-S方程对上述流场进行有限元迭代, 用该迭代结果作为神经网络计算的输入样本, 对网络进行训练。输入层为两神经元, 输入向量为待计算点的x, y方向坐标;输出层为两神经元, 输出向量为计算点的水平速度分量和垂直速度分量。对于甲板表面的任意一点, 将其坐标值输入神经网络, 即可计算得到该点的气流速度分量, 从而可以确定处于该点的旋翼气动力特性, 并可以得到该点的气流变化情况对旋翼气动特性的影响。

图2所示为机库后方、甲板上方的气流流动情况。图中箭头方向代表该点的气流速度方向。从图2中可以看出, 气流从机库上方 (y > 0位置) 掠过后, 在机库后方形成下冲, 甲板区域气流出现倒流, 并呈旋涡状。而直升机在甲板区域起降时, 正处于该涡状流动内, 随着直升机在甲板表面的移动, 受甲板涡流的影响, 直升机旋翼气动特性不断变化, 其操纵特性也随之变化^[4], 增加了飞行员的操纵难度。

《图2》

图2飞行甲板流场

Fig.2 Airflow above the flightdeck

使用神经网络进行仿真时, 计算了甲板区域着舰点上方部分点的速度分布, 计算结果与实验对比如表3所示。

表3着舰点上方下冲气流垂直速度神经网络部分计算结果与实验结果比较 (来流速度为18m/s)

Table 3 Compare of vertical speed above the landing point of test with the output of neural network (free-flow velocity =18 m/s)

《表3》

距甲板高度 / m	试验结果 /m·s^-1	计算结果 /m·s^-1	相对误差 / %
2.4	2.43	2.09	14
5.4	2.98	3.34	12.1
7.9	3.17	2.76	13

由于实验数据有限, 为了充分验证神经网络计算结果的准确性, 在前面已经验证N-S方程计算结果准确性的基础上, 又将神经网络的计算结果与N-S方程的计算结果做了比较, 见表4和表5。

《6 结论》

6 结论

从上述比较可以看出, 神经网络对流场的模拟结果与实验结果及N-S方程的模拟结果具有较好的一致性, 而且计算时间大为缩短。因此, 为满足仿真需要, 可以首先将流场划分为比较大的网格, 使用N-S方程对流场进行初步模拟, 其结果作为神经网络的训练样本对神经网络进行训练, 发挥神经网络可以模拟任意点流场的优势, 将训练后的神经网络模型结合到直升机实时仿真程序中。

表4N-S方程与神经网络水平速度计算比较 (来流速度为18 m/s)

Table 4 Compare of horizontal speed calculation of N-S equations with the output of neural network (free-flow velocity=18 m/s)

《表4》

坐标x, y (无量纲)	N-S方程 /m·s^-1	BP模型 /m·s^-1	相对误差 /%
-1.38, 0.75	10.940	11.180	2.2
-1.05, 1.0	13.520	13.650	1.0
-0.72, 1.25	15.440	15.500	0.4
0.829, -0.096	0.745	0.747	0.3
1.373, -0.35	0.218	0.247	13.3
2.459, -0.452	0.518	0.509	1.7
3.0, -0.096	1.325	1.393	5.1

表5N-S方程与神经网络垂直速度计算比较 (来流速度为18 m/s)

Table 5 Compare of vertical speed calculation of N-S equations with the output of neural network ( free-flow velocity=18 m/s)

《表5》

坐标x, y (无量纲)	N-S方程 /m·s^-1	BP模型 /m·s^-1	相对误差 /%
-1.27, 1.5	0.55 E-2	0.540 E-2	1.8
-1.05, 1.0	0.189 E-1	0.198 E-1	0.5
-0.36, 0.25	0.265 E-1	0.234 E-1	11.7
0.829, -0.401	0.179 E-0	0.207 E-0	15.6
1.918, -0.305	0.571 E-1	0.558 E-1	2.3
1.918, -0.226	0.493 E-1	0.558 E-1	13.2
0.829, -0.305	0.185 E-0	0.202 E-0	9.2

展示更多