直升机以其独特优越的低空、低速飞行能力,在军用及民用方面均得到了广泛应用,例如海上救援、森林救火、低空巡航监测等。随着直升机可应用领域愈加广阔,已成为航空领域中不可替代的一部分。然而,也正因其低空低速的飞行包线,以及低空任务的广泛多样性,直升机更易遭受低空大气环境影响,从而影响飞行性能甚至威胁飞行安全。通过深入研究大气扰动对飞行器的危害,主要有两种风类型影响飞行器性能:风切变与大气湍流[1]。
风切变是指平均风在一段时间或空间上的变化,湍流是指叠加在平均风上的连续随机脉动。在众风切变形式中,最危险的是水平漩涡一般在600 m高空下的微下击暴流。其具有风速变化多样性、风切变强度变化剧烈的特点,并考虑到直升机的低空低速低能量的飞行特点,微下击暴流对直升机的威胁不容忽视。大气湍流会降低直升机的飞行性能及飞行品质、增加结构载荷振荡,甚至引起驾驶员诱发振荡,增大驾驶难度,影响飞行安全。
微下击暴流作为一种最危险、较简单的空气流动,已有较多成熟的模型建立方案。主要有3种建模方式,其一是建立风场数据库,采用内插法取值,所需数据量大,且难以直观展示风切变特性,不利于定性分析;其二是利用简化的数学模型描述以适用于工程研究的工程化模型[2-3],主要作为水平距离的函数,风场三维特性不足;其三是根据流体力学和热力学规律建立并求解大气动力学方程[4-5]。结合本文研究目的,选取第3种方法完成风场的建模,生成的风场可较好地体现风场的三维空间变化特性,有利于分析风切变对直升机的本质影响。此外,也可调整建模参数达到不同风切变强度,适应性较强。
为增强风场的不均匀性,进一步分析湍流与风切变的耦合影响,在基础风场上叠加湍流流场。直升机湍流流场的建模也有诸多实现方案。McFarland[6]发展了一种可用于飞行仿真的旋翼叶素紊流仿真(Simulation of Rotor Blade Element Turbulence, SORBET)模型,该模型仅考虑了包含旋翼平面内的二维湍流场,以此为基础,吉洪蕾等[7]进一步发展了新的三维空间湍流流场,本文对该模型所使用的高斯插值算法进行了改进,以研究直升机对大气湍流的响应特性。
目前国内外对固定翼风切变威胁的研究较完善成熟[3, 8-9],且系统地提出威胁因子[10]等定量结论,广泛应用于风切变预警中,此外,Dogan和Kabamba[11]进行了飞机以不同操纵策略飞出微下击暴流的研究,并提出建议的风场逃离策略,以降低坠机事故率。然而风切变对直升机威胁分析的相关文献较少,起步较晚。高华[2]建立了三维组合风切变,并研究了不同方向风切变对直升机的影响,但未考虑风切变量级变化的情况。Liu等[12]分析了在风场不同位置及高度侧向风与垂向风对机体性能的影响,但仅考虑了机体以高速穿越风场的影响,并未考虑涵盖整个飞行包线的飞行速度在风场中的潜在威胁。在文献[4, 12-13]中,风切变速度项一般叠加于机体质心处,而未考虑由于风切变相对于机体的转动,或是通过对纵向运动方程求导加入风切变项[14]的方式,不适合于仿真分析,且涉及求导,不利于计算。
本文利用涡环法建立了三维微下击暴流的流场模型,并叠加了可用于实时仿真的三维湍流场。为提高计算精度,捕捉风切变的切变特性,选取特征点建立了含大气扰动的直升机飞行动力学模型,并配备相应的姿态保持增稳控制系统。通过对比分析直升机以不同飞行速度、从不同位置穿越风场的仿真结果,并结合理论推导,得到了直升机状态量变化与风场的关系,总结了含湍流的风切变场对直升机的潜在威胁因素。
1 微下击暴流流场模型涡环的诱导速度场与微下击暴流风场形式相似,故利用涡环原理可构造微下击暴流风场模型,图 1为建模原理示意图,地面上下对称布置强度为Γ的主涡环及镜像涡环,特别的,由于微下击暴流风场的下沉气流并不一定垂直于地面,其外流流谱具有明显的非对称性[1],涡环面与地面存在倾角。图中, Oxyz为地面坐标系,并分别建立主涡环坐标系OPxyz,其原点OP在O(xP, yP, zP)处;镜像涡环坐标系OLxyz,其原点OL在O(xP, yP, -zP)处;M为参考质点,其坐标为OP(xM, yM, zM)。风场坐标系与地面坐标系相同。
1.1 涡环诱导速度场计算为便于分析与计算,在涡环坐标系下计算诱导速度场。首先考虑一般情况,即M点距离涡环较远。设涡环半径为Rv,r1、r2分别为点M距涡环最近及最远的距离值,引入k=(r2-r1)/(r2+r1),已知当0≤k2≤1时,主涡环的流线方程可用式(1)逼近[2]:
$ {\phi _{\rm{v}}} = - \frac{\varGamma }{{2\pi }}({r_1} + {r_2})\frac{{0.788{k^2}}}{{0.25 + 0.75\sqrt {1 - {k^2}} }} $ | (1) |
由流线方程得出涡环在M点处的诱导速度vM的分量分别为
$ {{v_x} = \frac{1}{{{r_M}}} \cdot \frac{{\partial {\phi _{\rm{v}}}}}{{\partial z}} \cdot \frac{{{x_M}}}{{{r_M}}}} $ | (2) |
$ {{v_y} = \frac{1}{{{r_M}}} \cdot \frac{{\partial {\phi _{\rm{v}}}}}{{\partial z}} \cdot \frac{{{y_M}}}{{{r_M}}}} $ | (3) |
$ {{v_z} = - \frac{1}{{{r_M}}} \cdot \frac{{\partial {\phi _{\rm{v}}}}}{{\partial {r_M}}}} $ | (4) |
式中:rM为M点距涡环中心轴的距离。
接着,计算涡环中心轴线处的诱导速度。因为此时rM为零,无法通过式(2)~式(4)求解,利用涡环的位函数可推导得,在涡环的中轴线处,水平方向速率为0,垂直方向的速率为
$ {v_z} = \frac{\varGamma }{{2{R_{\rm{v}}}}} \cdot \frac{1}{{{{(1 + {{({z_M}/{R_{\rm{v}}})}^2})}^{1.5}}}} $ | (5) |
最后利用Rankine涡原理计算涡丝附近的诱导速度,将涡核看作半径为r的圆环,涡核内部的流速沿半径呈线性分布,而涡核外部仍服从流线方程。
如图 2所示,若点M在涡核内部,记点Or为过点OP、M的垂直平面与涡丝的交点,点N为平面与涡环的交点。点N的位置由点M通过定比分点公式求得,即
$ \overrightarrow {{O_r}N} = \lambda \overrightarrow {{O_r}M} \;\;\;{\kern 1pt} \lambda > 0 $ | (6) |
因为点N位于涡核的边界处,亦满足流线方程,故根据流线方程可求得点N的诱导速度vN,继而点M的诱导速度
$ {v_M} = \frac{{{v_N}}}{\lambda } $ | (7) |
地面风场分布由上下涡环综合得到。考虑到涡环绕自身中心轴(z轴)旋转不影响流场分布,可设主涡环相对地面坐标系旋转倾角为[ϕv θv 0]T,相应的镜像涡环的旋转倾角为[-ϕv -θv 0]T, ϕv与θv的定义见图 1。引入旋转矩阵L(ϕv, θv)、L(-ϕv, -θv)分别实现从主涡环坐标系到地面坐标系的转换、从镜像涡环坐标系到地面坐标系的转换。最终生成风场轴系(地面轴系)风速矢量即为两者速度场的叠加,即
$ \left[ {\begin{array}{*{20}{l}} {{W_{Mx}}}\\ {{W_{My}}}\\ {{W_{Mz}}} \end{array}} \right] = \mathit{\boldsymbol{L}}({\phi _{\rm{v}}},{\theta _{\rm{v}}})\left[ {\begin{array}{*{20}{l}} {{v_{Px}}}\\ {{v_{Py}}}\\ {{v_{Pz}}} \end{array}} \right] + \mathit{\boldsymbol{L}}( - {\phi _{\rm{v}}}, - {\theta _{\rm{v}}})\left[ {\begin{array}{*{20}{l}} {{v_{Lx}}}\\ {{v_{Ly}}}\\ {{v_{Lz}}} \end{array}} \right] $ | (8) |
式中:WMx、WMy、WMz分别为水平风、垂向风、侧向风风速。
1.3 流场建模参数设置及其三维分布根据联合机场天气研究(Joint Airport Whether Studies, JAWS)计划收集的实际微下击暴流风场强度与空间尺度的统计分析,设置模型的基本参数为:主涡环高度为610 m,涡环半径为915 m,涡核半径r=400 m,中心轴处垂直速率为12 m/s,涡环无倾角。表 1为风切变模型参数与高频特征参数范围对比,可看出,所建立的模型参数与一般的风切变特征参数相符。
风切变特征参数 | 高频参数范围 | 模型参数 |
风切变尺度/m | 1 830~3 660 | 1 800 |
最大水平风速变化 (150 m处)/(km·h-1) |
37~83 | 83.61 |
风速变化强度/(m·s-1·m-1) | 0.016 7~0.026 7 | 0.01 7 |
涡环中心截面处水平风及垂向风相比风场侧面,风速最大、风切变强度最强,图 3为随高度h变化的水平风与垂向风的分布图。由图可得,从上至下,垂向风强度逐渐较小,水平风强度逐渐增大,符合风场变化规律。
接着分析侧向风的变化趋势。图 4为300 m高度处,侧向位置y取0~1 800 m时的侧向风风速剖面分布。由图可得,中心截面处的侧向风速度为0 m/s,在y=900 m处侧向风速度最大,速度变化范围为0~10 m/s。
1.4 风切变风场接口一般计算旋翼气动力有两种方式:①仅考虑旋翼桨毂上的大气风速,即假定风速在桨盘上均匀分布;②利用叶素理论,分别计算旋翼各叶素处含大气扰动的相对来流速度。前者较为简单,但未能充分体现出风切变的切变特性;后者的计算量较大,但计算结果较为精确。
风速在一定范围内变化不大,为充分考虑其切变特性,且在精确建模且不增加计算量的前提下,参考飞机的4点模型[15],由于直升机旋翼基本覆盖了全机尺度,在旋翼上选取4个特征点计算风切变强度表征全机强度,其示意图见图 5。
旋翼桨盘上的风切变强度可等效为直升机相对空气扰动的旋转角速度[pgs qgs rgs]T,结合图 5可推导得
$ {{p_{{\rm{gs}}}} = \frac{{\partial {W_{{\rm{s}}z}}}}{{\partial y}} = \frac{{{W_{{\rm{s}}bz}} - {W_{{\rm{s}}dz}}}}{{2R}}} $ | (9) |
$ {{q_{{\rm{gs}}}} = - \frac{{\partial {W_{{\rm{s}}z}}}}{{\partial x}} = - \frac{{{W_{{\rm{s}}cz}} - {W_{{\rm{s}}az}}}}{{2R}}} $ | (10) |
$ \begin{array}{*{20}{c}} {{r_{{\rm{gs}}}} = {r_{{\rm{gs1}}}} + {r_{{\rm{gs2}}}} = - \frac{{\partial {W_{{\rm{s}}x}}}}{{\partial y}} + \frac{{\partial {W_{{\rm{s}}y}}}}{{\partial x}} = }\\ {\frac{{{W_{{\rm{s}}dx}} - {W_{{\rm{s}}bx}}}}{{2R}} + \frac{{{W_{{\rm{s}}cy}} - {W_{{\rm{s}}ay}}}}{{2R}}} \end{array} $ | (11) |
式中:各点的气流速度定义于旋翼轴系;Wsaz表示点a处气流速度的z向分量,其他相似符号以此类推;R为旋翼半径。
将旋翼轴系旋转角速度矢量转换至机体轴系可得
$ \left[ {\begin{array}{*{20}{l}} {{p_{\rm{g}}}}\\ {{q_{\rm{g}}}}\\ {{r_{\rm{g}}}} \end{array}} \right] = {\mathit{\boldsymbol{L}}_{{\rm{SB}}}}\left[ {\begin{array}{*{20}{l}} {{p_{{\rm{gs}}}}}\\ {{q_{{\rm{gs}}}}}\\ {{r_{{\rm{gs}}}}} \end{array}} \right] $ | (12) |
式中:LSB表示旋翼轴系到机体轴系的转换矩阵。
桨毂处的风速WsH由4点的风速度平均求得
$ {\mathit{\boldsymbol{W}}_{{\rm{sH}}}} = \frac{1}{4}({\mathit{\boldsymbol{W}}_{{\rm{s}}a}} + {\mathit{\boldsymbol{W}}_{{\rm{s}}b}} + {\mathit{\boldsymbol{W}}_{{\rm{s}}c}} + {\mathit{\boldsymbol{W}}_{{\rm{s}}d}}) $ | (13) |
所以,机体相对大气扰动的旋转角速度可等效表示为
$ \left[ {\begin{array}{*{20}{c}} {{p_{{\rm{rel}}}}}\\ {{q_{{\rm{rel}}}}}\\ {{r_{{\rm{rel}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {p - {p_{\rm{g}}}}\\ {q - {q_{\rm{g}}}}\\ {r - {r_{\rm{g}}}} \end{array}} \right] $ | (14) |
即考虑风切变影响的机体状态量可表示为
$ \mathit{\boldsymbol{x}} = {\left[ {\begin{array}{*{20}{l}} u&v&w&{{p_{{\rm{ rel }}}}}&{{q_{{\rm{ rel }}}}}&{{r_{{\rm{ rel }}}}}&\phi &\theta &\psi \end{array}} \right]^{\rm{T}}} $ |
机体平尾、垂尾、机身、尾桨等各部件的风速均由其在地轴系中的位置代入风场模型求得。
2 湍流模型及其三维扩展根据军用品质规范MIL-F-8785C建立符合平稳随机过程的Dryden紊流模型。利用时间序列数组生成包围机体的三维空间边界,空间内部采用插值算法计算各点湍流速度。
2.1 二维平面大气湍流流场的生成首先以旋翼所在平面为例说明二维大气湍流流场的生成。如图 6所示,构造一个覆盖整个直升机旋翼平面的长方形,且固定于机体上。长方形短边AB与旋翼桨尖平面相切,且垂直于直升机平飞速度V。长方体的宽度为旋翼直径,长度L大于直升机总长度,并平均分成N段。A、B两点处分别放置一套大气湍流滤波器,每套湍流滤波器根据Dryden模型生成离散大气紊流速度。随着直升机以水平速度V前飞,因为大气紊流相对地面静止,所以生成的大气紊流时间序列UA、UB将按照一定的采样时间分布于两侧边AD、BC上。
矩形长度L可表示为
$ L = NV\Delta t $ | (15) |
利用一长度为N的数组记录生成的湍流时间序列,在t=nΔt时刻,其与湍流分段位置的映射关系如图 6所示。湍流分段位置k与数组位置m的对应关系为
$ m = (n + 1 - k)\% N $ | (16) |
式中:“%”为取余符号。
桨叶模型采用叶素分段法,则桨叶i上叶素j在面ABCD中的位置为
$ \begin{array}{*{20}{l}} {{x_{i,j}} = R + {r_j}{\kern 1pt} {\rm{cos}}{\kern 1pt} {\kern 1pt} {\psi _i}}\\ {{y_{i,j}} = {r_j}{\kern 1pt} {\rm{sin}}{\kern 1pt} {\kern 1pt} {\psi _i}} \end{array} $ | (17) |
式中:rj为桨叶i上叶素j到桨毂轴的径向长度。
由xi, j可以计算得到桨叶气动中心所在横截线与两边交点E、F处的紊流速度分别为UA(mi, j)、UB(mi, j),且
$ {m_{i,j}} = \left( {n - \left\lfloor {\frac{{{x_{i,j}}}}{{L/N}}} \right\rfloor + 1} \right)\% N + 1 $ | (18) |
式中:“
应用反距离加权插值(Inverse Distance Weighted)算法计算各叶素气动中心处的大气紊流速率,则
$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{W}}_{{\rm{T}}i,j}} = }\\ {\quad \;\;\;\frac{{{\mathit{\boldsymbol{U}}_A}({m_{i,j}})/{{(R + {y_{i,j}})}^2} + {\mathit{\boldsymbol{U}}_B}({m_{i,j}})/{{(R - {y_{i,j}})}^2}}}{{1/{{(R + {y_{i,j}})}^2} + 1/{{(R - {y_{i,j}})}^2}}}} \end{array} $ | (19) |
图 7展示了直升机浸入湍流流场过程中各叶素及近桨尖处叶素的湍流速度变化历程。仿真中各叶素按等圆环面积法分为5段,桨叶片数取为4片。前小段浸入风场阶段,仅桨叶前缘处有湍流速度,随着进一步前飞,各叶素湍流速度分布逐渐扩宽,桨尖处湍流速度较前缘处存在一定的时间延迟,左右两端湍流速度存在明显差异,各桨叶的湍流速度由于使用插值计算,基本在由左右两端曲线包围的范围内变化。
此外,根据图 7可看出,与桨毂处的湍流速度相比,各叶素的紊流速度高频部分明显,实际上,由于旋翼的旋转,使得整个桨盘平面内的紊流速度分布是随机的,各个叶素受到的紊流扰动的高频分量会相互消减。
2.2 三维扩展为得到其他部件的湍流流场,依据二维湍流流场将其扩展至三维。在直升机底部构造与旋翼平面ABCD平行的二维湍流生成面MNOP,与旋翼面距离为zH,则两个面共同组成了包围直升机的三维长方体,且类似的,在前端M、N两点放置湍流滤波器。各点气动力计算亦采用反距离加权插值算法,以机身为例说明。
假设机身在湍流坐标系下的位置坐标为[xf yf zf]T,则该部件的湍流速度WTf可由4边的大气紊流时间序列UA、UB、UM、UN表示为
$ {\mathit{\boldsymbol{W}}_{{\rm{Tf}}}} = \frac{{\frac{{{\mathit{\boldsymbol{U}}_A}({k_{\rm{f}}})}}{{l_A^2}} + \frac{{{\mathit{\boldsymbol{U}}_B}({k_{\rm{f}}})}}{{l_B^2}} + \frac{{{\mathit{\boldsymbol{U}}_M}({k_{\rm{f}}})}}{{l_M^2}} + \frac{{{\mathit{\boldsymbol{U}}_N}({k_{\rm{f}}})}}{{l_N^2}}}}{{1/l_A^2 + 1/l_B^2 + 1/l_M^2 + 1/l_N^2}} $ | (20) |
式中:
$ {{k_{\rm{f}}} = \left( {n - \left\lfloor {\frac{{{x_{\rm{f}}}}}{{L/N}}} \right\rfloor + 1} \right)\% N + 1} $ | (21) |
$ {\left\{ {\begin{array}{*{20}{l}} {l_A^2 = {{(R + {y_{\rm{f}}})}^2} + z_{\rm{f}}^2}\\ {l_B^2 = {{(R - {y_{\rm{f}}})}^2} + z_{\rm{f}}^2}\\ {l_M^2 = {{(R + {y_{\rm{f}}})}^2} + {{({z_{\rm{H}}} - {z_{\rm{f}}})}^2}}\\ {l_N^2 = {{(R - {y_{\rm{f}}})}^2} + {{({z_{\rm{H}}} - {z_{\rm{f}}})}^2}} \end{array}} \right.} $ | (22) |
在含湍流的风切变风场模型中,假设湍流模型与微下击暴流模型相互独立。然而,实验表明,在风切变中,湍流特征长度随风切变的大小与强度以未知的方式增加,所以,本文的湍流强度设为σT=2.1 m/s,表征严重的湍流等级。
3 飞行动力学模型详细的建模过程可参考文献[16-18],本文建立了一种单旋翼带尾桨直升机通用的、精度较高的高阶非线性飞行动力学数学模型。旋翼模型采用Pitt-Peters的一阶谐波动态入流模型计算旋翼的诱导速度[19-20],通过求解挥舞运动学方程[21]计算桨叶挥舞角及挥舞角速率。由各片叶素处的来流速度计算得各叶素的迎角、侧滑角及来流马赫数,对风洞试验数据插值得到翼型的升力系数CL及阻力系数CD,继而可得各片叶素上作用的气动力,最终求得整个旋翼的气动力和气动力矩。机身、平尾、垂尾的相对来流速度均考虑了旋翼下洗、侧洗的影响,并通过风洞试验数据插值得各部件上作用的气动力及气动力矩。采用Bailey模型计算尾桨的拉力和扭矩。
模型运动学方程可表示为
$ \mathit{\boldsymbol{\dot x}} = f(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{\delta }},t) $ | (23) |
式中:x=[u v w p q r ϕ θ ψ],表示体轴系下机体各状态量;操纵杆量δ=[δcol δlon δlat δped]分别表示总距、纵横向周期变距及尾桨操纵。
3.1 模型验证将飞行力学模型计算得到直升机配平结果与稳态飞行试验数据[22]对比, 如图 8~图 10所示。用于配平的直升机总重为7 257 kg,飞行高度为1 600 m。从图中可以看出,本文计算结果与飞行试验结果吻合良好,误差基本保持在10%之内。小速度下配平计算结果与飞行试验结果相差较大,总距操纵量与需用功率均小于飞行试验数据,主要有两方面因素,一方面为在小速度时保持稳定飞行非常困难,所以飞行试验的误差不可避免;另一方面为本文所采用的动态流入模型不能有效捕捉小速度飞行时的尾迹收缩效应。尾桨操纵量与飞行试验数据保持约5%的误差。
综上,可认为本文建立的飞行力学模型满足仿真计算需要。
3.2 增稳控制系统直升机自身的不稳定性,尤其是在飞越风场时,配备飞行控制增稳系统以改善飞行品质、提高飞行安全显得尤为必要[23]。因此,仅考虑直升机飞过风场时来自风场变化导致的响应威胁,排除机体自身的不稳定性,将具有姿态保持功能的增稳控制系统集成在动力学模型中。
本文采用的飞控系统主要根据文献[16-17]构建。其为样机提供全包线全权限实时控制,包括了4个部分:内环提供角速率阻尼的增稳系统;提供纵向稳定性的俯仰偏置舵机;具有姿态保持及空速保持的飞行航迹稳定系统; 能进行初级操纵解耦的混合器。
4 飞行动力学仿真 4.1 无湍流的风切变场为研究直升机飞越微下击暴流风场的威胁因素,选取不同速度、不同高度、不同侧向位置出发作为控制变量进行仿真,以40 m/s飞行速度、300 m飞行高度、飞越风场中心的飞行状态为基准,分析各历程中姿态角(以俯仰角θ为主)、地速[Ue Ve We]T、高度的变化量,其中姿态角与地速以初始状态为基准原点考虑其差值,高度以地面为参照物,由地轴系下的z轴方向的坐标值表示,并进行了对比飞行试验。
首先进行不同飞行速度的比较,选取3个较为典型的前飞速度V,即20 m/s、40 m/s、60 m/s,分别代表低速飞行、巡航飞行以及高速飞行。以300 m飞行高度从距离风场中心1 800 m处出发从中心剖面飞越风场。飞越风场的过程中,姿态角变化不大,其中俯仰角变化幅度相对较为明显,如图 11所示。其水平、垂直地速以及高度随飞行距离变化的历程曲线如图 12所示。由图中可以看出,不同速度的地速变化曲线接近,且与图 3中水平风速的变化曲线基本重合,这说明在空速保持控制器的增稳控制下,机体基本可以较好地跟踪风速变化。换言之,直升机优异独特的水平机动性能,使得水平风切变对直升机的影响较小。
接着观察图 12中垂向速度与高度变化曲线。可以看出,在初期飞行高度变化不大时,不同前飞速度下的垂向速度变化相似,且与下降风风速相似。随着进一步的前飞,高度降低,随之遭遇的垂向风也减弱,因此机体下降速度减缓。此外,小速度飞行时对垂向下降风更敏感,也更易遭遇威胁,主要是因为相同的下降速度下,小速度飞行所需时间更长,导致下降高度更大。因此,从高度变化上看,随前飞速度增大,高度变化减小。
为进行不同飞行高度的比较,选取了3个较为典型的高度h,即150 m、300 m、450 m,分别代表水平风较弱而垂向风较强、水平风及垂向风均衡以及垂向风较强而水平风较弱的飞行条件。直升机以40 m/s的平飞速度飞越风场中心,各主要状态量变化如图 13、图 14所示。从图 13中可看出,水平风主要影响俯仰角,但俯仰角变化幅度正常,其余姿态角变化幅度更小,不予赘述。地速变化幅度在初期由于高度相差较大,因此区别较明显,随着进一步飞行,高度差异缩小,地速变化趋于一致。类似的,在高度更高的位置,垂向下降速度越大,但综合不同高度飞行的高度变化历程可以看出,高度越高,虽然遭遇更大速度下降气流,但是最终飞出风场时仍处于较高位置,换言之,在遭遇风切变时,较高的高度更安全。
最后是不同飞行侧方位的比较,结合图 4选取了3个较为典型的方位,即y=±900 m、y=0 m。y=±900 m处为左右侧向风最大侧向位置处,主要验证机体对不同方向侧向来流的对称性,以及侧向风的威胁程度。y=0 m处无侧向风作为对照基准组。机体飞行速度为40 m/s,飞行高度为300 m。
偏航角及地速变化历程如图 15和图 16所示。在侧向风较大的区域,由于距离风场中心较远,水平风及垂向风影响显著减弱,侧向风占主导地位,因此,姿态角中偏航角变化最明显,其变化历程如图 15所示。由图可看出,机体的航向稳定性较好,且随侧方位变化呈现出显著的对称性。图 16为地速变化曲线,与中心截面变化历程相比,水平及垂向地速变化显著减小,威胁降低;侧向运动速度变化范围小,且由于直升机的航向稳定性作用,侧向风影响基本可控。
所以,若遭遇微下击暴流,应向远离风场中心,垂向风更小的区域进行规避,可有效降低风场威胁。
4.2 飞行器运动与风场关系推导本节结根据直升机的动力学方程,推导含有风速及风切变项的运动学方程,从宏观理论角度,解释验证由仿真推测出的直升机状态量变化与风场的关系。
在直升机的气流轴系下,定义机体运动速度为Va,加速度为
则由动力学方程可得
$ {\mathit{\boldsymbol{\dot V}}_{\rm{a}}} + {\mathit{\boldsymbol{\omega }}_{\rm{a}}} \times {\mathit{\boldsymbol{V}}_{\rm{a}}} = \mathit{\boldsymbol{A}}/m + \mathit{\boldsymbol{g}} - ({\mathit{\boldsymbol{\dot W}}_{\rm{a}}} + {\mathit{\boldsymbol{\omega }}_{\rm{a}}} \times {\mathit{\boldsymbol{W}}_{\rm{a}}}) $ | (24) |
式中:A表示气动力项;g表示重力项。
在上文的仿真中,由于姿态保持增稳控制器的作用,水平及横向的气动力项变化与操纵有关,且纵横向耦合强,解耦困难; 另一方面,由上述仿真结果可看出,机体可跟踪横纵向风速变化,不是引起威胁的主要因素。
相反的,垂向风与其他通道耦合较小,未配备垂向通道保持的增稳控制器,机体响应体现了在裸机状态下直升机受垂向风的影响。此外,由上文仿真可推出,垂向风对机体的威胁占主导作用,因此有必要着重分析垂向通道的动力学方程。
将式(24)展开后可得
$ \begin{array}{l} {V_{{\rm{a}}z}} + {p_{\rm{a}}}{V_{{\rm{a}}x}} - {q_{\rm{a}}}{V_{{\rm{a}}y}} = g{\rm{cos}}{\theta _{\rm{w}}}{\rm{cos}}{\phi _{\rm{w}}} - {L_{\rm{a}}}/m + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {W_{{\rm{a}}z}} + {p_{\rm{a}}}{W_{{\rm{a}}x}} - {q_{\rm{a}}}{W_{{\rm{a}}y}} \end{array} $ | (25) |
式中:升力项可表示为
$ {L_{\rm{a}}} \approx T = \frac{1}{2}\rho \pi {R^2}{(\varOmega R)^2}{C_T} $ | (26) |
忽略横纵向操纵量对拉力系数CT的影响,则CT可表示为总距杆量与垂向空速的函数[24],即
$ {C_T} = F({\theta _0}) + \frac{1}{2}\kappa {a_\infty }\sigma \left( {\frac{{{V_{{\rm{a}}z}} - {W_{{\rm{a}}z}}}}{{\varOmega R}}} \right) $ | (27) |
记:
$ \begin{array}{*{20}{l}} {{P_{{\rm{ const }}}} = \frac{1}{2}\rho \pi {R^2}(\varOmega R)2F({\theta _0}) - g{\rm{cos}}{\theta _{\rm{w}}}{\rm{cos}}{\phi _{\rm{w}}}}\\ {{P_{{\rm{ varia }}}} = \frac{1}{{4m}}\kappa {a_\infty }\rho \pi {R^2}(\varOmega R)} \end{array} $ | (28) |
考虑到姿态保持器可保证机体姿态角变化量在正常平衡范围内小幅波动,故ωa×Va、ωa×Wa项可看作小量,暂时忽略其影响。
式(25)化简并移项后可得
$ {\dot V_{{\rm{a}}z}} - {\dot W_{{\rm{a}}z}} = - {P_{{\rm{ const }}}} - {P_{{\rm{ varia }}}}({V_{{\rm{a}}z}} - {W_{{\rm{a}}z}}) $ | (29) |
式(29)可看作一元非齐次线性常微分方程。求解得(Vaz-Waz)关于时间的变化函数:
$ ({V_{{\rm{a}}z}} - {W_{{\rm{a}}z}})(t) = - \frac{{{P_{{\rm{ const }}}}}}{{{P_{{\rm{ varia }}}}}} + {P_{{\rm{ coef }}}}{{\rm{e}}^{ - {P_{{\rm{ varia }}}}^t}} $ | (30) |
式中:Pcoef与初始状态量有关。
在上文的仿真中,直升机以水平飞行的配平状态飞入风场,则此时的配平状态可表示为
$ \begin{array}{*{20}{l}} {{V_{{\rm{a}}z}} = {W_{{\rm{a}}z}} = 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{m/s}}}\\ {T \approx mg{\kern 1pt} {\kern 1pt} {\rm{cos}}{\theta _{\rm{w}}}{\rm{cos}}{\phi _{\rm{w}}}} \end{array} $ | (31) |
若总距θ0保持定值,代入式(28)可得
$ {P_{{\rm{ const }}}} = 0 $ | (32) |
设初始时刻遭遇的垂向风大小为Waz0,则
$ {V_{{\rm{a}}z}}(t) = - {W_{{\rm{a}}z0}}{{\rm{e}}^{ - {P_{{\rm{ varia }}}}^t}} + {W_{{\rm{a}}z}}(t) $ | (33) |
由式(33)可知:机体垂向运动速度在经历短暂的过渡后,与当前位置处的垂向风速保持一致。
从合作用力平衡角度分析,当机体垂向运动速度与垂向风速相等时,垂向空速为0,与初始配平状态相同,此时,机体在气流轴系下作平飞运动,与4.1节仿真得到的结论相符。
同理,该分析方法也可用于纵横向通道的计算及验证。
4.3 叠加湍流的风切变场上文分析在无湍流的情况下,风切变对直升机的影响特性,由于实际的风场中必有湍流场存在,且湍流可降低飞行器飞行性能,甚至影响飞行安全,因此有必要验证湍流与风切变场共同作用的影响,分析两者是否会由于耦合作用产生更严重的威胁。
通过仿真得到的含湍流的风切变的主要状态量变化如图 17所示。飞行条件为:直升机以40 m/s的平飞速度从300 m高度、距离中心点1 800 m处从中心剖面飞越含湍流的微下击暴流风场。图 17中展示了3条曲线的变化历程, 观察仅由湍流作用的曲线变化,发现其主要激励姿态角小幅高频震荡,但基本在平衡点附近。类似的,垂向速度变化的震荡亦体现出高频的效果,而水平风速的变化则由于积分作用表现更为平缓。图中的插值拟合曲线为将同样飞行条件下,通过仿真得到的微下击暴流与湍流分别作用的时间历程曲线按相同水平位移做插值拟合并合并得到的曲线。比较该条拟合曲线与仿真得到的变化历程,可发现两者较为重合,这说明湍流与风切变的影响相互独立,无明显耦合,换言之,湍流主要影响直升机小幅短期高频的姿态角等变化,而宏观的风切变则主要诱导机体速度、高度等状态量的大幅低频变化。
综上所述,由于湍流与微下击暴流风场的作用相互独立,且风切变在机体状态量变化中占主导作用,以上对无湍流的风切变场的分析在含湍流的风场中同样适用。
5 结 论1) 为捕捉风切变的切变项,在不增加计算量的前提下,发展了可适用于直升机飞行动力学的三维风切变风场模型,并在风场中加入了三维湍流模型,提高了直升机在风切变气流场中的动态响应计算精度。
2) 分析了不同风场位置、飞行速度等直升机的响应。在增稳系统的辅助作用下,水平风及侧向风对飞行安全威胁较小。垂直气流是直升机在微下击暴流中的主要威胁来源,可导致同等幅度的机体下降速度,且与机体飞行速度无关,因此,水平飞行速度越慢,下降高度越多,坠地威胁越强。
3) 遭遇风切变时,提高飞行高度以及向侧向远离风场中心可有效降低风场对直升机威胁。提高飞行高度可增加高度裕度,降低坠地可能性;向风场侧向规避可有效减弱水平风尤其是垂向风的影响,侧向风的威胁较弱。
4) 湍流与微下击暴流风场对直升机的影响相互独立,且湍流主要引起高频小幅姿态角的震荡,总体而言对直升机的威胁次于微下击暴流对速度等状态量的作用。
[1] |
肖业伦, 金长江. 大气扰动中的飞行原理[M]. 北京: 国防工业出版社, 1993: 166-173. XIAO Y L, JIN C J. Flight theory in atmosphere turbulence[M]. Beijing: National Defense Industry Press, 1993: 166-173. (in Chinese) |
[2] |
高华.风切变对直升机飞行特性的影响[D].南京: 南京航空航天大学, 2009. GAO H. Research on helicopter flight characteristics in wind shear[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2009(in Chinese). |
Cited By in Cnki (6) | Click to display the text | |
[3] | FAA. Windshear training aid. Volume 1. overview, pilot guide, training program[R]. Washington, D.C.: Federal Aviation Administration, 1987. |
[4] | IVAN M. A ring-vortex downburst model for flight simulations[J]. Journal of Aircraft, 1986, 23(3): 232-236. |
Click to display the text | |
[5] |
高振兴.复杂大气扰动下大型飞机飞行实时仿真建模研究[D].南京: 南京航空航天大学, 2009. GAO Z X. Research on real-time flight simulation of large aircraft in complex atmospheric disturbance[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2009(in Chinese). |
Cited By in Cnki (42) | Click to display the text | |
[6] | MCFARLAND R E. Element simulation of helicopter turbulence: USA, 5860807[P]. 1999-01-09. |
Click to display the text | |
[7] |
吉洪蕾, 陈仁良, 李攀. 适用于直升机飞行力学分析的三维空间大气紊流模型[J]. 航空学报, 2014, 35(7): 1825-1835. JI H L, CHEN R L, LI P. A model of three dimension-al field atmospheric turbulence for helicopter flight dynamics analysis[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(7): 1825-1835. (in Chinese) |
Cited By in Cnki (9) | Click to display the text | |
[8] | FROST W, BOWLES R L. Wind shear terms in the equations of aircraft motion[J]. Journal of Aircraft, 1984, 21(11): 866-872. |
Click to display the text | |
[9] | BOWLES R L. Windshear detection and avoidance: Airborne systems survey[C]//29th IEEE Conference on Decision and Control. Piscataway: IEEE, 1990: 708-736. |
Click to display the text | |
[10] | PROCTOR F H, HINTON D A, BOWLES R L. A windshear hazard index[C]//9th Conference on Aviation, Range and Aerospace Meteorology. Orlando: American Meteorology Society, 2000: 482-487. |
Click to display the text | |
[11] | DOGAN A, KABAMBA P T. Escaping microburst with turbulence:Altitude, dive, and pitch guidance strategies[J]. Journal of Aircraft, 2000, 37(3): 417-426. |
Click to display the text | |
[12] | LIU T, DAI Y, HONG G. Flight dynamic simulation of helicopter forward flight through microburst wind field[J]. Advances in Mechanical Engineering, 2017, 9(2): 168781401769121. |
Click to display the text | |
[13] |
洪冠新, 庞健. 风切变场中直升机前飞状态动态响应[J]. 北京航空航天大学学报, 2005, 31(5): 524-528. HONG G X, PANG J. Helicopter dynamic response to wind shear in forward flight[J]. Journal of Beijing University of Aeronautics and Astronautics, 2005, 31(5): 524-528. (in Chinese) |
Cited By in Cnki (11) | Click to display the text | |
[14] |
赵维义, 傅百先. 低空风切变中直升机纵向运动特性分析[J]. 飞行力学, 2002, 20(1): 25-28. ZHAO W Y, FU B X. Analysis of the helicopter longitudinal motion characteristics in the low-level windshear[J]. Flight Dynamics, 2002, 20(1): 25-28. (in Chinese) |
Cited By in Cnki (9) | Click to display the text | |
[15] | ETKIN B. Turbulent wind and its effect on flight[J]. Journal of Aircraft, 1981, 18(5): 327-345. |
Click to display the text | |
[16] | HOWLETT J J. UH-60A Black Hawk Engineering Simulation Program: Volume I-mathematical model: NASA CR 166309[R]. Washington, D.C.: NASA, 1984. |
Click to display the text | |
[17] | HOWLETT J J. UH-60A Black Hawk Engineering Simulation Program: Volume II-background report: NASA CR 166310[R]. Washington, D.C.: NASA, 1988. |
[18] | HILBERT K B. A mathematical model of the UH-60 Helicopter: NASA-TM-85890[R]. Washington, D.C.: NASA, 1984. |
Click to display the text | |
[19] | GAONKAR G, PETERS D. Review of dynamic inflow modeling for rotorcraft flight dynamics[C]//27th Structures, Structural Dynamics and Materials Conference, 1986. |
Click to display the text | |
[20] | PETERS D A, HAQUANG N. Technical note:Dynamic inflow for practical applications[J]. Journal of the American Helicopter Society, 1988, 33(4): 64-68. |
Click to display the text | |
[21] | TALBOT P D, TINLING B E, DECKER W A, et al. A mathematical model of a single main rotor helicopter for piloted simulation: NASA TM-84281[R]. Washington, D.C.: NASA, 1982. |
Click to display the text | |
[22] | BALLIN M G. Validation of a real-time engineering simulation of the UH-60A helicopter: NASA TM 88360[R]. Washington, D.C.: NASA, 1987. |
Click to display the text | |
[23] |
黄子安, 韩智修, 朱家祯. 直升机三轴增稳改善风切变作用下的飞行品质[J]. 南京航空航天大学学报, 1996(3): 98-103. HUANG Z A, HAN Z X, ZHU J Z. A three-axis stability augmentation system for improving helicopter flying qualities under effect of wind shear[J]. Journal of Nanjing University of Aeronautics and Astronautics, 1996(3): 98-103. (in Chinese) |
Cited By in Cnki (3) | Click to display the text | |
[24] |
高正, 陈仁良. 直升机飞行动力学[M]. 北京: 科学出版社, 2003. GAO Z, CHEN R L. Helicopter flight dynamics[M]. Beijing: Science Press, 2003. (in Chinese) |