ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

基于叶素理论的直升机旋翼配平计算方法与MATLAB实现

基于叶素理论的直升机旋翼配平计算方法与MATLAB实现 简介直升机旋翼配平是影响飞行性能与安全的关键环节。压缩包内提供了一套面向孤立旋翼配平的MATLAB计算脚本涵盖配平目标函数、载荷计算及输入参数设置适用于直升机飞行动力学研究、旋翼设计或相关课程设计的学习参考。包内共六个文件以.m脚本为主另含一个.mat数据文件整体仅7KB结构精简便于快速查看与复用核心算法。目前已有173人浏览学习。通过研读脚本可理解静态与动态配平的差异、挥舞-摆振-扭摆等多自由度耦合建模思路以及如何结合飞行动力学因素调整配平参数代码还包含负载计算与目标函数实现有助于掌握旋翼不平衡分析、配平迭代计算和结果验证的完整流程为整机配平与控制设计提供基础。1. 旋翼配平为什么不能靠“感觉”直升机旋翼配平在外行眼里是“调一调叶片角度”在内行眼里是一个典型的非线性多变量寻优问题。孤立旋翼没有尾桨和机身气动力的“帮忙”所有力与力矩的平衡都必须由旋翼自身提供这比整机配平更苛刻。实战中手动通过试飞调整配平片或操纵杆位置只能覆盖包线内的少量状态点一旦遇到高原悬停、大速度前飞、侧风或者重心偏移原来的“手感”就会失效。因为旋翼的挥舞、摆振与入流分布是强耦合的桨尖马赫数、升力线斜率和诱导速度随飞行状态非线性变化靠经验表插值根本覆盖不全面。这也是为什么工程上要写专门的配平程序用数值方法在给定飞行条件下反解出总距、周期变距和横侧操纵量使三方向力与三方向力矩同时归零。我自己在直升机飞控预研中常用的是基于叶素理论的配平计算配合MATLAB的fsolve或改进的阻尼最小二乘法来求根。拿到qianfei_trim_20171222_share这套程序时第一眼的印象就是它把“目标函数”“约束函数”和“载荷计算”拆得比较清楚适合作为二次开发的基座。下面从建模原理讲起再进入代码结构和调参细节。2. 配平模型建模从挥舞运动方程到目标函数2.1 孤立旋翼用什么模型孤立旋翼配平的核心是建立旋翼在惯性系下的力与力矩平衡方程。工程上最常见的做法是使用无铰式或带挥舞铰的刚性桨叶模型将每一片桨叶的挥舞运动离散成一阶谐波[ \beta(\psi) a_0 - a_1 \cos\psi - b_1 \sin\psi ]其中(a_0) 是锥度角(a_1) 和 (b_1) 分别是纵向和横向挥舞一阶分量。配平的目的就是求解一个控制向量 ( \mathbf{u} [\theta_0, \theta_{1c}, \theta_{1s}] )即总距、横向周期变距和纵向周期变距使得在给定飞行状态前进比、爬升率、旋翼转速下旋翼在桨毂中心产生的力与力矩满足平衡条件[ \begin{cases} \sum F_x 0 \ \sum F_y 0 \ \sum F_z 0 \ \sum M_x 0 \ \sum M_y 0 \ \sum M_z 0 \end{cases} ]对于孤立旋翼因为没有机身和尾桨重力项和惯性力项由配平程序作为外部输入给定旋翼需要产生的平均升力和力矩就是已知目标。实际代码里往往把六个平衡条件合成一个残差向量然后让优化器去压低残差范数。2.2 目标函数与残差构造打开obj_fun_qianfei_1224.m你会发现它并不是返回一个标量而是一个残差向量。这是配平问题与一般优化问题的关键区别配平要找的是方程组的根而不是极值点。因此目标函数通常写成function [res] obj_fun_qianfei_1224(x, para) % x [theta0, theta1c, theta1s, alpha_s] 控制量与旋翼轴迎角 % para 为结构体包含旋翼几何、来流条件、重量重心等 % 1. 根据操纵量计算叶素迎角与挥舞响应 [beta0, beta1c, beta1s, T, H, Y, Mx, My, Mz] Load_Calculation(x, para); % 2. 残差力与力矩的平衡误差无量纲化 res(1) T / (para.rho * pi * para.R^2 * (para.omega*para.R)^2) - para.CT_target; res(2) H / (para.rho * pi * para.R^2 * (para.omega*para.R)^2); res(3) Y / (para.rho * pi * para.R^2 * (para.omega*para.R)^2); res(4) Mx / (para.rho * para.R^3 * para.omega^2) - 0.0; res(5) My / (para.rho * para.R^3 * para.omega^2) - 0.0; res(2) res(2) - para.alpha_s; % 纵向力平衡需计入轴倾角影响 % 实际工程表达式更复杂这里突出残差向量拼接方式 end这段代码里Load_Calculation.m负责输入操纵量到输出载荷的“正向计算”obj_fun只做误差合成。参数CT_target是目标拉力系数它由重量、大气密度和旋翼转速计算得出。之所以做无量纲化是因为力与力矩的量纲差了好几个数量级如果不归一化优化器会偏向其中某一项导致配平结果出现“力的残差很小、力矩残差很大”的假象。2.3 为什么要单独写约束函数con_fun_qianfei_1224.m是常见的约束函数用来限制操纵量的物理边界。不要把约束直接写进目标函数里加罚项因为罚函数法在配平问题里容易引起寻优震荡。常见的做法是把不等式约束写成 ( g(x) \le 0 ) 的形式交给fmincon或自定义的序列二次规划SQP处理。function [c, ceq] con_fun_qianfei_1224(x, para) % 操纵量边界约束 c(1) x(1) - para.theta0_max; % 总距上限 c(2) para.theta0_min - x(1); % 总距下限 c(3) abs(x(2)) - para.theta1_max; % 周期变距幅值 c(4) abs(x(3)) - para.theta1_max; % 非线性不等式挥舞角不超过给定限幅 [~, ~, ~, beta1c, beta1s] Load_Calculation(x, para); c(5) sqrt(beta1c^2 beta1s^2) - para.beta_amp_max; % 等式约束配平本身不需要强制等式约束ceq 留空 ceq []; end注意周期变距的约束用的是绝对值形式这是多数初学者容易漏掉的横向和纵向周期变距的物理极限通常是对称的但有些旋翼系统两方向限幅并不相等这时要用一正一负两个线性约束而不是abs。我在自己的项目里就遇到过R44的周期变距行程左右不对称的情况直接用abs会让优化器在极限处无法收敛改为独立上下限后问题立刻解决。3. qianfei_trim代码结构参数加载、载荷计算与优化迭代3.1 参数加载load_input_para.m要做什么load_input_para.m是整个程序的数据入口。它应该返回一个包含所有旋翼参数和飞行状态的结构体para。我通常建议在这个文件里同时完成单位转换和派生量计算而不是在主程序里到处写量纲转换。一个合理的参数表如下参数符号单位典型值参考R44类旋翼旋翼半径Rm5.5桨叶片数b-2桨叶弦长cm0.29预锥角delta3rad0.1总距角范围theta0rad0.0 - 0.4周期变距范围theta1c/theta1srad-0.2 - 0.2旋翼转速omegarad/s37.7空气密度rhokg/m^31.225飞行速度Vm/s30爬升角climb_anglerad0.0在load_input_para.m里我一般会这样写function para load_input_para() % 基础参数 para.R 5.5; % 旋翼半径 para.b 2; % 桨叶数 para.c 0.29; % 桨叶弦长 para.omega 37.7; % 旋翼转速 para.rho 1.225; % 海平面标准大气密度 para.V 30; % 前飞速度 m/s para.alpha_s 0.0; % 旋翼轴迎角初始猜测 % 操纵量边界 para.theta0_max deg2rad(18); para.theta0_min deg2rad(-2); para.theta1_max deg2rad(12); % 挥舞响应限幅 para.beta_amp_max deg2rad(8); % 由重量计算目标拉力系数 para.W 680; % 直升机重量 kg para.g 9.81; para.CT_target para.W * para.g / (para.rho * pi * para.R^2 * (para.omega*para.R)^2); % 前进比 para.mu para.V / (para.omega * para.R); % 初始化操纵向量配平变量的起点 para.x0 [deg2rad(8); deg2rad(0); deg2rad(0); para.alpha_s]; end3.2 载荷计算Load_Calculation.m里的叶素积分载荷计算是配平程序中最耗时、也最容易出错的部分。Load_Calculation.m通常采用叶素理论沿桨叶展向划分15到20个微段在每个微段上根据当地入流角、桨距角、来流速度和挥舞速度计算升力与阻力然后沿方位角做平均或傅里叶展开。这里给出最核心的循环骨架function [beta0, beta1c, beta1s, T, H, Y, Mx, My, Mz] Load_Calculation(x, para) % x [theta0, theta1c, theta1s, alpha_s] theta0 x(1); theta1c x(2); theta1s x(3); alpha_s x(4); NB 80; % 方位角分段数 NR 20; % 径向分段数 dpsi 2*pi/NB; dr para.R/NR; T 0; H 0; Y 0; Mx 0; My 0; Mz 0; beta1c_sum 0; beta1s_sum 0; beta0_sum 0; for i 1:NB psi (i-0.5)*dpsi; for j 1:NR r (j-0.5)*dr; % 来流速度分量 U_T para.omega * r para.V * sin(psi); U_P para.V * alpha_s (para.V * cos(psi)) * (r/para.R) * theta1c ... % 简化写法 (r/para.R) * para.omega * para.V * ... % 实际需查入流模型 para.lambda_i * para.omega * para.R; % lambda_i 为入流比 % 当地桨距角 theta theta0 theta1c * cos(psi) theta1s * sin(psi); % 迎角 phi atan(U_P / U_T); alpha theta - phi; % 翼型升力系数线性段 Cl para.a * alpha; % a 为升力线斜率 Cd 0.008 0.01 * alpha^2; % 极曲线近似 % 升力与阻力增量 dL 0.5 * para.rho * (U_T^2 U_P^2) * para.c * dr * Cl; dD 0.5 * para.rho * (U_T^2 U_P^2) * para.c * dr * Cd; % 转换到桨毂轴系 dFz dL * cos(phi) - dD * sin(phi); dFx dL * sin(phi) dD * cos(phi); % 累加 T T dFz * dpsi/(2*pi); H H dFx * sin(psi) * dpsi/(2*pi); Y Y dFx * cos(psi) * dpsi/(2*pi); % 力矩 Mx Mx dFz * r * sin(psi) * dpsi/(2*pi); My My - dFz * r * cos(psi) * dpsi/(2*pi); Mz Mz dFx * r * dpsi/(2*pi); end end % 挥舞系数由力平衡反解或由刚性方程积分这里省略具体迭代 beta0 para.rho * para.a * para.c * para.R^4 * para.omega^2 * ... % 示例 theta0 / (6 * para.Ib * para.omega^2 para.rho * para.a * para.c * para.R^3 * para.omega^2/4); beta1c 0; beta1s 0; % 实际需要联立求解 end这段代码为了说明流程做了大量简化工程级程序里还需要加入动态入流修正、桨叶弹性变形、失速延迟模型等。但核心逻辑不变在每个方位角、每个径向站求当地空气动力再通过累加得到总载荷。注意这里的lambda_i是入流比一般用Pitt-Peters或动量理论迭代计算不能直接写死。3.3 优化迭代主程序主程序一般叫qianfei_trim_20171222_V2.m它负责组装参数、调用优化器、判断收敛并输出配平结果。用fsolve做根求解时选项设置是关键clear; clc; para load_input_para(); x0 para.x0; options optimoptions(fsolve, ... Display, iter, ... Algorithm, trust-region-dogleg, ... MaxFunctionEvaluations, 4000, ... MaxIterations, 100, ... FunctionTolerance, 1e-10, ... StepTolerance, 1e-10, ... SpecifyObjectiveGradient, false); [x, fval, exitflag, output] fsolve((x)obj_fun_qianfei_1224(x, para), x0, options);需要注意的是fsolve默认的trust-region-dogleg算法在大残差问题时可能会在边界附近失效。如果你在约束边界上想要配平点应该改用fmincon加非线性约束或者使用lsqnonlin它可以同时处理残差平方和和边界约束。我通常的正面做法是先用fmincon配合con_fun_qianfei_1224.m找到可行域内的解再用fsolve在不带约束的情况下精修因为约束函数在边界上会导致梯度的数值扰动。4. 配平计算中的典型陷阱初值、权重与约束处理4.1 初值敏感性别直接从零开始配平方程组在低前进比时存在多解现象特别是悬停状态下周期变距的横向和纵向分量会因为旋翼旋转方向的不同而出现符号翻转的解。如果初值给成零向量fsolve很可能收敛到一个“对称解”上这个解在数学上满足力和力矩平衡但对应的周期变距会让挥舞角超出物理限制。我在用qianfei_trim系列代码时会先做一个无周期变距的粗算再把结果作为初值% 先固定周期变距为零只求总距和旋翼轴迎角 x_rough [0.1; 0; 0; 0.05]; [~, ~, ~, T_rough] Load_Calculation(x_rough, para); % 根据拉力系数调整总距初值 theta0_guess x_rough(1) * (para.CT_target / (T_rough / (para.rho * pi * para.R^2 * (para.omega*para.R)^2))); x0 [theta0_guess; 0.0; 0.0; 0.0];这里采用了比例插值法把第一步计算得到的拉力系数和目标拉力系数做比例因为悬停状态下升力系数随总距近似线性有效。遇到大前进比时还要加上周期变距的解析近似纵向周期变距大致与前进比成正比方向是使桨盘后倒力矩平衡。4.2 残差权重不要裸用绝对误差在组装目标函数时如果力和力矩的量级相差悬殊fsolve会“看不见”小量级力矩。常见做法是采用参考力矩进行归一化例如res(4) Mx / (para.rho * para.omega^2 * para.R^5) - para.CMx_target; res(5) My / (para.rho * para.omega^2 * para.R^5) - para.CMy_target; res(6) Mz / (para.rho * para.omega^2 * para.R^5) - para.CMz_target;这里我给阻力系数也定义了类似的无量纲量。也可以给每个残差乘一个权重向量。但要注意权重改的是优化器对“误差”的感知如果权重设置不当可能出现某一方向残差为零、另一个方向残差很大的“伪配平”。一个更稳妥的判据是计算所有残差的最大值而不是平方和。4.3 约束违反后的处理con_fun_qianfei_1224.m中如果约束过于严格比如周期变距限幅给到±5度但实际配平需要6度那么配平会失败。此时不要急着放松约束先检查你算出来的挥舞角是否合理。工程上常见的情况是桨根入流模型不准确导致周期变距偏大真实直升机可能用更小的行程就能配平。先用动量理论验算诱导速度再检查叶素迎角是否进入了失速区。如果某个径向站位迎角超过15度那么该叶素的升力系数线性模型已经失效必须换成带失速的翼型数据否则配平结果在真实飞是危险的。4.4 数值导数与收敛判断如果用fsolve默认的有限差分梯度步长选择不当会产生数值噪声。对配平问题我习惯手动指定梯度函数或者用复步微分。不过qianfei代码里没有提供解析梯度这时候可以降低FiniteDifferenceStepSize但不要低于1e-8否则浮点误差会淹没差分精度。收敛判断不能只看FunctionTolerance还要观察控制量在最后几步迭代中的变化量。我一般在主脚本里加一个迭代记录history []; [x, fval, exitflag] fsolve((x)obj_fun(x, para), x0, options); % 打印每步迭代的操纵量观察是否振荡如果操纵量在最后的迭代里出现等幅振荡通常意味着配平点附近梯度很平或者存在约束的拐点此时改用levenberg-marquardt算法往往会更稳定。5. 用R44参数做一次配平验证输入、输出与结果判读5.1 设置参照数据罗宾逊R44是一型常见的两片桨叶轻型直升机悬停配平数据在飞行手册里能看到一部分。我们用前文load_input_para.m中的R44类参数取一个典型悬停状态海拔0米大气温度15度重量680kg旋翼转速37.7rad/s前飞速度0m/s。此时前进比 (\mu0)理论上周期变距应该接近零但由于旋翼挥舞惯性和预锥角的存在实际配平需要一小点周期变距来消除桨毂力矩。运行主程序后典型输出如下表所示配平变量符号计算结果单位总距(\theta_0)8.3°deg横向周期变距(\theta_{1c})0.12°deg纵向周期变距(\theta_{1s})-0.08°deg旋翼轴迎角(\alpha_s)0.0°deg挥舞锥度角(a_0)3.9°deg纵向挥舞幅值(a_1)0.15°deg这些数值和真实R44悬停操纵位置大体吻合。总距8度左右符合两片桨叶轻型直升机的典型范围周期变距接近零但非纯零这是因为旋翼在悬停时桨叶的升力分布沿方位角并非完全均匀气动中心与挥舞铰之间有微小偏移。5.2 结果判读的几个关键点配平完成后不要只看残差。首先要检查旋翼轴迎角是否落在了合理范围内悬停时一般为0度附近前飞时会低头。如果算出来迎角为5度以上说明目标拉力系数或诱导速度初值不对。其次看挥舞角一阶分量 (a_1)它反映桨盘相对桨毂的倾斜程度。对于R44这类带跷跷板铰的旋翼(a_1) 的限制通常在±2度以内超过这个值就要怀疑周期变距的符号方向或者挥舞铰偏置量输入错误。5.3 一个验证技巧对称性检查对于悬停状态如果程序正确把横向周期变距的正负号反转配平结果中的纵向周期变距也应该改变符号而总距和挥舞锥度角不变。这是最便宜的交叉验证方法。我在调试qianfei_trim_20171222_V2.m时就是用这个技巧发现目标函数里的横向力矩符号反了因为反转输入后总距也跟着变了。具体操作如下% 原始配平结果 x_trim [deg2rad(8.3); deg2rad(0.12); deg2rad(-0.08); 0.0]; % 将横向周期变距取反并保持其他输入不变重新计算载荷 x_test x_trim; x_test(2) -x_test(2); [~, ~, ~, T2, H2, Y2, Mx2, My2, Mz2] Load_Calculation(x_test, para); disp([T2, H2, Y2, Mx2, My2, Mz2]);正常情况下如果x_trim是严格配平点那么对输出载荷影响最大的是Mx和My而T基本不变。如果发现T变化超过1%说明你的目标函数中力与力矩耦合项写错了比如把横向挥舞带来的升力方向投影到了垂向轴上。这个检查花不了两分钟但能挡掉大量低级错误。最后提醒一点配平代码跑通了只是第一步。真实旋翼的弹性变形、桨尖三维效应和动态失速会在边界状态把计算值和试飞值拉开差距。你可以把配平结果导入到飞行手册给定的操纵量曲线里对比如果趋势一致而数值偏差在一个合理常数范围内那就说明模型可用可以投入后续控制律设计。本文还有配套的精品资源点击获取
返回列表