适于直升机飞行仿真的高原大气紊流模型

2022-09-05吉洪蕾苏俊杰陈仁良孔卫红

航空学报 2022年7期

吉洪蕾，苏俊杰，陈仁良，孔卫红

1. 重庆大学航空航天学院，重庆 400044 2. 南京航空航天大学航空学院直升机旋翼动力学国家级重点实验室，南京 210016

中国西部高原地区具有重要的国防战略地位，其海拔高、地势险峻、大风天气多，独特的地理气候环境导致大气紊流出现的概率以及紊流的强度和频率宽度都远高于一般平原地区，直升机在该地区执行近地飞行任务时受到大气紊流的干扰，飞行品质下降、甚至危及飞行安全，准确掌握该地区大气紊流对直升机运动的动态干扰特性显得十分重要且迫切。然而，目前仍然缺乏一种适于直升机飞行仿真的高原大气紊流模型。

大气紊流对航空器飞行品质的影响始终是航空界重点关注的飞行安全问题之一。经过多年发展，用于固定翼飞机飞行仿真的大气紊流建模技术已经日趋成熟。与之相比，适于直升机飞行仿真的大气紊流模拟技术则略有滞后，主要原因之一是直升机升力主要由旋翼的旋转产生，大气紊流对直升机的干扰机理迥异于固定升力体的飞机，适于固定翼飞机的大气紊流建模方法对于直升机来说不再合适。为解决这一问题，Gaonkar等提出桨叶紊流建模方法，该方法考虑各片桨叶气动中心的相关性生成大气紊流速度，并将其叠加各桨叶相对来流速度上建立用于直升机飞行仿真的大气紊流模型。该方法虽然符合大气紊流对直升机旋翼的气动干扰机理，但是其计算效率不适合实时仿真应用。Mcfarland和Duisenberg提出了旋翼叶素紊流建模方法。该方法在旋翼前方固定两个大气紊流成型滤波器，基于“冻结场”假设，随着直升机前飞则大气紊流速度向后扩展，形成覆盖整个旋翼的大气紊流速度场。吉洪蕾等拓展了旋翼叶素紊流方法，在直升机前方垂直平面内布置多个紊流成型滤波器，基于“冻结场”假设，随直升机前飞形成三维空间紊流模型。该方法由于兼顾广泛的适用性和较高的计算效率而更加适于直升机飞行仿真应用，但是其采用低精度的德莱顿紊流谱模型进行数值模拟，且没有考虑到横切面内大气紊流速度的空间相关性，因而计算精度受到限制。

除了以上基于物理机理的模型，控制等效紊流输入(Control Equivalent Turblence Input,CETI)模型提供了一种数据驱动的直升机大气紊流建模方法。其以具有特定频谱的信号叠加到操纵输入产生与大气紊流干扰等效的直升机频域响应，模型参数由飞行试验数据辨识。CETI模型受益于其基于飞行试验数据的特点，具有高计算精度且满足实时仿真要求，但是也受限于飞行试验数据缺乏，仅适合悬停、小速度飞行状态。

直升机由于其独特的飞行特点，常常在低空执行机动任务，经历的大气紊流环境常常比固定翼飞机恶劣得多，影响大气紊流统计特性的因素也多。我国西部高原地区地貌复杂多变、昼夜温差大、大风天气多的特点导致其大气紊流活动更加活跃，低空大气紊流受下垫面地形地貌、大气层结及平均风速的影响更加显著，固定翼飞机飞行品质规范给出的仅随高度变化的紊流统计参数无法反映以上因素的影响，成为限制高原大气紊流环境直升机飞行仿真及飞行品质分析的主要问题之一。

本文发展一种考虑我国西部高原地形地貌、大气层结和平均风速廓线影响的高原大气紊流模型，针对已有三维空间大气紊流模型精度低、无法考虑横切面内紊流速度空间相关性的问题，采用高精度的冯·卡门紊流谱函数推导大气紊流速度的高阶成型滤波器及其时间序列生成算法，提出空间相关的随机大气紊流速度时间序列生成方法解决横切面内紊流速度的空间相关性问题，提升大气紊流模拟的计算精度，并通过数值仿真和与试验数据对比验证模型。

1 适于飞行仿真的高原大气紊流模型

1.1 高原大气紊流模型

低空大气紊流，尤其是600 m高度以下的大气紊流受下垫面地形地貌、大气层结及平均风速的影响显著。因此，首先建立平均风速模型。平均风速是指在一定时距内风速变化的平均值。600 m高度以下大气边界层内的平均风速度受到地表摩擦力和大气层结状态的影响，风速随高度增加而增大。对于均一地貌下垫面，采用地面粗糙度表征地面摩擦力对风速的影响，指数律描述的大气边界层内风速廓线为

(1)

式中：为高度的平均风速度;为参考高度(通常为10 m)的风速;为风速廓线指数。

在高于600 m的大气边界层外，可以认为风速随高度不再变化，即

=≥600

(2)

式中：为600 m高度的平均风速度。

风速廓线指数受到下垫面粗糙度、大气层结状态等因素的影响，采用式(3)计算:

=+lg+(lg)

(3)

式中：、、为风速廓线指数系数，其取值随大气层结状态变化。国外研究人员测量了多个地区不同大气层结状态的风速廓线指数，采用测试数据拟合的方式获得如表1所示的风速廓线指数系数。

表1 风速廓线指数系数取值[18]Table 1 Coefficients for wind speed profile exponent[18]

大量实测数据表明，高原大气紊流速度的无量纲频谱在惯性子区满足-2/3次方衰减的卡莫夫条件，在低频段满足趋于零的极限条件，与冯·卡门紊流谱模型的基本理论一致，并且我国多个地区大气紊流速度实测频谱也与冯·卡门谱函数吻合良好。因此，采用冯·卡门紊流谱模型描述高原大气紊流速度随机变化的频域特征。基于大气紊流速度的一维空间频谱函数以及“冻结场”假设，由角频率和空间频率形式的紊流谱转化关系，可以得到角频率形式的大气紊流谱函数为

(4)

式中：()、()、()为纵、横和垂向紊流速度谱函数;,,为3个方向紊流强度;=(为速度;为积分尺度)为特征频率。在紊流强度保持一定的条件下，特征频率决定了紊流能量的频域分布，随着值的增大，大气紊流速度呈现出更多高频成分。

由式(4)的紊流谱函数可知，大气紊流强度和积分尺度是大气紊流速度随机变化的主要统计特性参数，通过调整两者的变化即可描述不同频率分布和幅值的大气紊流。在300 m高度以下的大气边界层内，紊流活动发展强盛，并受地面摩擦、大气层结状态等因素的影响，随高度、地面粗糙度、风速等因素变化，紊流强度和积分尺度采用国外统计的近地层变化规律来描述:

(5)

(6)

在600 m高度以上的自由层，大气紊流不再受到下垫面的影响，可以认为其强度和积分尺度不再随高度变化，紊流强度和积分尺度根据飞行品质规范给出

===01

(7)

===280

(8)

在300～600 m高度，紊流强度和积分尺度由大气边界层内取值到自由大气层取值逐渐过渡:

(9)

(10)

式中：300、300、300为300 m高度的大气紊流强度，其表达式分别为

由式(3)、式(5)和式(6)可以看出，高原大气紊流模型考虑到了下垫面地貌和大气层结状态引起的平均风速度、紊流强度和积分尺度变化及其对大气紊流统计特性的影响规律。

1.2 大气紊流速度生成算法

由于冯·卡门紊流谱函数的无理式形式，直接将其用于生成大气紊流速度序列较为困难。固定翼飞机飞行品质规范针对冯·卡门谱给出了一种三阶有理形式的拟合结果，该结果的有效无量纲频率为50 rad，对于飞行速度较高的固定翼飞机是合适的。与之相比，直升机飞行速度远低于固定翼飞机，三阶有理形式的拟合结果已无法满足直升机悬停/低速状态飞行品质评估精度的要求。采用最小二乘法对式(4)的冯·卡门谱拟合，获得四阶有理形式的大气紊流谱函数，共轭分解之后得到高阶大气紊流成型滤波器如式(11) 所示:

(11)

式中：()、()、()分别为纵、横和垂向高阶紊流成型滤波器。系数,,,分别为

针对式(11)的高阶紊流成型滤波器，采用离散变换的方式获得大气紊流速度时间序列的离散递推生成算法。由于零极点匹配法离散精度高，甚至在奈奎斯特频率附近也能较好地保持频域离散精度，采用该方法对式(11)的纵向紊流滤波器()离散后得到

Δ+3=2Δ+2+1Δ+1+0Δ+(+2+1+1+0)

(12)

式中：和Δ分别表示=Δ时刻的白噪声输入和纵向紊流速度输出。式(12)各系数分别为

且纵向紊流成型滤波器的零点和极点为

式中：=。

同理，可以得到横、垂向大气紊流速度时间序列生成算法为

Δ+4=3Δ+3+2Δ+2+1Δ+1+0Δ+

,(+3+2+2+1+1+0)

(13)

式中：Δ表示在=Δ时刻的横向或垂向紊流速度分量。式(13)各系数分别为

式中:横、垂向紊流速度成型滤波器的零点和极点为

其中：,=,。

由式(12)和式(13)的离散递推算法可知，该大气紊流速度时间序列生成算法计算简单、效率高并具备较高的紊流频谱，更加适于低空、低速飞行的直升机飞行仿真研究。

1.3 适于飞行仿真的三维空间紊流模型

基于大气紊流速度时间序列，生成覆盖整个直升机的三维空间紊流场，如图1中长方体所示。为了计算每一时刻直升机各气动面处的大气紊流速度，定义了两个坐标系：一个是地轴系，其中轴指向地球正北方向为正，轴指向地球正东方向为正，轴指向地心为正；另一个是紊流坐标系，其中轴指向空速的反方向，轴垂直于与轴指向机身左侧为正，轴垂直于平面指向机身上方为正。紊流场长度为，宽度为2，高度为，固定于机体并随直升机以速度前飞。紊流场前表面与旋翼桨尖平面相切，且垂直于直升机前飞速度。基于“冻结场”假设，将=×组大气紊流成型滤波器等间隔置于紊流场前表面，每套成型滤波器生成3个分量的大气紊流速度时间序列。当直升机以速度前飞时，由于大气紊流相对地面静止，所以置于前表面的紊流滤波器生成的大气紊流时间序列将以一定的延迟时间向直升机尾部拖出，布满整个大气紊流场，形成适于直升机飞行仿真的三维空间大气紊流场。

图1 三维空间紊流模型及坐标系定义Fig.1 Three-dimensional turbulence model and definition of coordinate axes

由于式(12)和式(13)的大气紊流速度时间序列生成算法是基于冯·卡门一维紊流谱函数获得的，因此在沿纵向的紊流速度时间序列满足纵向相关性，但无法满足各横截面内的相关性要求。随机统计理论早已证明，基于协方差矩阵Cholesky因子分解的线性变换方法对多维正态分布的相关变量随机数序列生成效果最好。因此，采用协方差矩阵Cholesky因子分解变换法将表面上相互独立的紊流成型滤波器相互关联。

假设置于表面的紊流成型滤波器分为行和列均匀排列。任意两个互相关的成型滤波器生成的纵向紊流速度序列之间的相关系数为

(14)

由式(14)可以得到所有组成型滤波器的相关系数矩阵=[(,)]，对其进行Cholesky分解可得下三角矩阵，即

(15)

Δ=Δ

(16)

式中：Δ为×组相互关联的大气紊流速度随机序列向量;Δ为×组独立紊流速度序列组成的向量。其表达式为

同理，可以生成×组互相关的横、垂向紊流速度时间序列，进而形成×组互相关的大气紊流成型滤波器。

将以上互相关的大气紊流成型滤波器布于大气紊流场前表面上，结合“冻结场”假设,建立的三维空间大气紊流场沿空间3个方向的分布均满足冯·卡门模型的空间相关性要求。

2 集成大气紊流的直升机飞行动力学模型

将上述三维空间大气紊流模型集成到已有高阶非线性直升机飞行动力学模型以开展飞行仿真研究，详细的飞行动力学模型描述及验证参见文献[12-13,25]。大气紊流速度主要通过干扰直升机各气动面的相对气流速度影响直升机气动载荷及运动。由地轴系到紊流坐标系的坐标变换矩阵为

(17)

(18)

(19)

式中：、、为空速在地轴系的纵、横和垂向分量。

以机身为例计算直升机各气动面在大气紊流坐标系的坐标。考虑到大气紊流场固定于直升机重心，机身在紊流坐标系中的坐标为

(20)

(21)

采用类似方法可计算出旋翼各叶素、平尾、垂尾和尾桨处大气紊流速度，采用线性叠加方式将其计入到直升机各气动面的相对气流速度上，形成集成高原大气紊流的直升机飞行动力学模型。以4片桨叶的黑鹰直升机为例，飞行动力学模型可以表示为

(22)

3 模型验证

3.1 高原大气紊流模型验证

采用四川省理塘县气象观测站实测的风速廓线和大气紊流度验证高原大气紊流模型，气象站所在区域海拔高度为3 932 m，地表为高原草甸，下垫面粗糙度约为0.05 m。

图2为理塘观测站实测风速廓线与模型预测结果的对比。观测数据分别在春季某天的02时、12时、18时和20时测得，预测风速廓线的大气层结状态分别为弱不稳定和中性状态，参考高度10 m 的平均风速分别为6.5 m/s、4.9 m/s、4.7 m/s 和3.1 m/s。从图2中可以看出，上午由于日照的逐渐增加，大气呈现出弱不稳定状态；下午到晚上由于日照逐渐减弱，大气层结状态趋于稳定，处于中性状态。按照给定的下垫面粗糙度和大气层结状态，风速廓线模型能够捕捉到高原大气层结状态和离地高度对平均风速的影响，预测结果可以与实测风速吻合的很好。

图2 理塘模型预测风速廓线与实测结果对比Fig.2 Comparison between predicted and measured wind speed profiles in Litang

大气紊流速度虽然呈现出很强的非定常性，但是其统计结果具有特定的分布规律。图3为理塘观测站实测大气紊流度与模型预测值的对比结果。从图3中可以看出，大气紊流度的观测结果存在多变性，尤其是在10 m以下高度，由于受到局部地貌和大气状态的干扰，大气紊流度的实测结果散布范围较大；随着高度增加，同一高度处的大气紊流度实测值趋于一致，呈现出随高度降低的规律。虽然观测结果存在多变性，但模型预测结果位于各散点中间，尤其是对直升机气动特性影响最大的垂向紊流速度分量，模型能够较好地预测大气紊流度的统计值及其变化趋势。由于大气紊流的随机性，采用具有代表性的统计特性参数开展大气紊流对直升机飞行仿真和飞行品质的影响研究具有重要参考价值。

图3 理塘模型预测大气紊流度与实测结果对比Fig.3 Comparison between predicted and measured turbulence intensities in Litang

生成时域大气紊流速度时间序列，计算其功率谱密度并与冯·卡门理论谱对比验证高原大气紊流模型的紊流速度生成算法。验证以四川理塘和西藏日喀则两地为例，其中，日喀则平均海拔超过4 000 m，地表主要是凹凸不平的丘陵地带，下垫面粗糙度约为0.3 m。

图4为仿真紊流速度频谱与冯·卡门理论频谱的对比结果。其中，参考高度10 m平均风速度为11.5 m/s，大气层结状态为中性，离地高度为12 m。在同样参考平均风速和大气层结状态的条件下，由于地形地貌的差别，日喀则地区粗糙度(=0.3 m)更高，低空大气紊流的幅值和频率宽度都远大于下垫面粗糙度更低(=0.05 m)的理塘地区大气紊流。仿真得到的大气紊流速度频谱与理论冯·卡门频谱模型吻合良好，证明大气紊流仿真算法能够准确捕捉高原地形地貌变化对大气紊流速度特性的影响。

图4 大气紊流速度仿真频谱与冯·卡门理论频谱对比Fig.4 Comparison between simulated and theoretical von Krmn atmospheric turbulence velocity spectra

3.2 三维空间大气紊流模型验证

以黑鹰直升机为例生成西藏日喀则地区上空的三维空间大气紊流验证仿真模型。黑鹰直升机旋翼半径为8.18 m，长度为19.76 m，高度为4.05 m。采用80个空间相关的紊流成型滤波器布于图1前表面生成三维空间大气紊流场，其中横向分布个数=40，垂向分布个数=2，参考高度10 m的平均风速度=11.5 m/s，日喀则地区的下垫面粗糙度=0.3 m，大气紊流场所在高度=12 m。

图5为某一时刻旋翼所在平面的垂向大气紊流速度场Δ。从图中可以看出，虽然每个空间位置的垂向大气紊流速度呈现出很强的随机性，但是整个大气紊流速度场表现出一定的空间相关性。

图5 垂向紊流速度场仿真结果Fig.5 Simulated vertical turbulence velocity field

图6为垂向大气紊流速度场不同滤波器生成的时域紊流速度功率谱密度与冯·卡门谱的对比。从图6中可以看出，虽然空间相关紊流成型滤波器生成的时域紊流速度序列并不完全相同(如图5所示)，但其频谱都与理论值一致，证明各空间相关紊流成型滤波器全部满足冯·卡门模型的纵向相关性要求。

图6 空间相关垂向紊流速度仿真频谱与理论谱对比Fig.6 Comparison between simulated and theoretical spectra of spatially collected vertical turblence velocity

图7为垂向大气紊流速度场沿横向的相关系数与理论值的对比。从图7中可以看出，由各紊流成型滤波器生成的大气紊流速度时间序列满足空间相关性的理论要求，证明了空间相关紊流成型滤波器生成方法的准确性。

图7 垂向紊流速度相关系数仿真结果验证Fig.7 Validation of simulated vertical turbulence velocity correlation coefficients

3.3 直升机对大气紊流的频域响应验证

以UH-60A黑鹰直升机为例模拟直升机对大气紊流的响应，计算其功率谱密度，并与飞行试验数据对比验证集成大气紊流的直升机飞行动力学模型。其中，飞行试验为UH-60A黑鹰直升机在机库后方大气紊流环境中保持12 m高度悬停，平均风速度为11.6 m/s，垂向大气紊流强度为1.37 m/s，机库高度为12 m。飞行仿真的参考高度10 m平均风速度=11.5 m/s，根据机库高度，下垫面粗糙度设为0.3 m，大气紊流场由80个空间相关的紊流成型滤波器生成，其中横向分布个数=40，垂向分布个数=2。飞行仿真中采用驾驶员模型实现直升机在大气紊流干扰环境的姿态稳定和位置保持，驾驶员模型的详细设计参见文献[25]，姿态稳定回路的穿越频率为2.5 rad/s。

图8为飞行仿真的直升机对大气紊流频域响应结果与飞行试验数据的对比。由于仅有3个方向角速度和垂向速度的飞行试验数据，因此与试验的对比也仅限于直升机对大气紊流响应的角速度和垂向速度。从图8中可以看出，虽然由于大气紊流的干扰直升机的响应特性随时间随机变化，但直升机对大气紊流的频域响应是固定分布的，模型预测的直升机角速度和垂向速度功率谱密度与飞行试验数据吻合良好，集成大气紊流的飞行动力学模型能够有效捕捉直升机对大气紊流的频域响应特性。