简介:这份资源面向航空工程领域研究人员、研究生及工程师,尤其适合关注倾转旋翼飞行器过渡段问题的专业人士,核心是XV-15倾转旋翼机的动力学建模与操纵特性分析。内容建立了考虑旋翼尾流对机翼、平尾、垂尾气动干扰的非线性动力学模型,借助MATLAB/Simulink完成配平、线性化与模态分析,并覆盖直升机模式与定翼机模式的稳定性、操纵特性以及过渡走廊确定,还延伸到航模建模验证。压缩包为1个docx文件,约56KB,内含论文复现说明与关键MATLAB代码及逐段解释,例如参数初始化、状态变量与控制输入定义、ODE动力学求解、trim配平与linmod线性化、特征根求解和响应绘图等,便于复现与二次扩展。已有144人学习。读者可据此理解过渡段气动与操纵机理,获得可运行的建模仿真脚本、图表化结果呈现以及后续高阶旋翼动力学与气动干扰模型的研究思路。
1. 倾转旋翼过渡段为什么是整条包线里最难仿真的一段
做倾转旋翼飞行器的人,几乎都会被同一段飞行卡住:短舱从 90° 转到 0° 的这一两分钟。直升机模式下旋翼是全部升力来源,固定翼模式下机翼接管升力,中间那段两边都不占优——旋翼拉力随短舱前倾,垂直分量往下掉,机翼升力要正好在同一个速度点上补上来,补早了抖振,补晚了掉高度,补错了就是俯仰发散。
XV-15 之所以几十年被反复拿出来建模,不是因为它新,而是因为它把这套耦合做成了可复现的样例:两副可倾转旋翼、大展弦比机翼、T 型尾翼,短舱角、总距、周期变距、升降舵四个通道在过渡段互相抢对方的活。动力学建模和操纵特性分析,绕不开这个平台。下面从状态量定义讲到配平、非线性仿真和过渡走廊的边界排错,代码都是能直接跑的形式。
2. XV-15 六自由度动力学建模:状态量、旋翼入流与部件气动
建模的第一步不是写方程,是决定哪些量进状态、哪些量当输入。过渡段最容易出错的地方,就是把短舱角当成了状态量——它其实是操纵输入,由飞控或飞行员按调度规律给出,模型里应该作为随时间变化的参数进入气动力计算,而不是积分出来。
2.1 机体坐标系、气流坐标系与 12 维状态向量
常用的做法是用机体轴系写刚体方程,气动力在部件本地坐标系算完再投影回机体轴。状态向量取 12 维:
| 分组 | 符号 | 含义 |
|---|---|---|
| 线速度 | u, v, w | 机体轴三轴速度,m/s |
| 角速度 | p, q, r | 机体轴滚转/俯仰/偏航角速度,rad/s |
| 姿态 | φ, θ, ψ | 欧拉角,rad |
| 位置 | x, y, z | 地轴系位置,m |
短舱角记作 β_n,0° 为固定翼模式,90° 为直升机模式。它和总距 θ0、横向周期变距 θ1s、纵向周期变距 θ1c、升降舵 δe、副翼 δa、方向舵 δr、左右油门一起构成控制向量。过渡段的关键是 β_n 的调度律,通常按指示空速或按时间线性推进,转速限制在每秒零点几度到几度之间。
提示:短舱角速率上限对配平连续性影响很大,建模时先固定速率再谈响应,不要把速率和配平混在一起调。
2.2 旋翼入流:动量理论打底,叶素理论补非均匀
旋翼是过渡段升力和推力的主要来源,入流模型选错,后面配平全歪。悬停时用动量理论就够了,前飞时诱导速度下降得很快,尤其在有迎角工况下,桨盘左右入流不对称会直接产生滚转力矩。
动量理论给的是拉力系数和诱导入流比的隐式关系:
CT = 2·λi·sqrt((λc + λi)² + μ²)
其中 μ 是前进比,λc 是爬升入流比,λi 是待求的诱导入流比。这是个一元非线性方程,用 Brent 法求根既快又稳:
import numpy as np from scipy.optimize import brentq RHO = 1.225 # 海平面空气密度 kg/m^3 R = 3.81 # XV-15 旋翼半径量级 m A = np.pi * R**2 # 单副旋翼桨盘面积 m^2 def inflow_ratio(CT, mu, lam_c, bracket=(1e-6, 0.5)): """动量理论求诱导入流比 lam_i CT = 2*lam_i*sqrt((lam_c+lam_i)^2 + mu^2) """ f = lambda lam: 2.0 * lam * np.sqrt((lam_c + lam)**2 + mu**2) - CT # 拉力为正时根一定在 (0, 0.5) 内,桨盘过载不会超过这个量级 return brentq(f, bracket[0], bracket[1], xtol=1e-10) def rotor_thrust(CT, Omega): """由拉力系数反算拉力 N,Omega 单位 rad/s""" return CT * RHO * A * (Omega * R)**2 if __name__ == "__main__": # 直升机模式悬停,桨尖马赫数量级下 CT 约 0.006~0.010 CT = 0.008 lam_i = inflow_ratio(CT, mu=0.0, lam_c=0.0) print("诱导入流比:", round(lam_i, 5)) print("旋翼拉力 N:", round(rotor_thrust(CT, Omega=61.7), 1)) # 589 rpminflow_ratio里bracket的选法很关键。λi 的物理上界由动量理论的理想极限给出,悬停时 CT = 2λi²,CT 取 0.01 时 λi 只有 0.07 左右,所以 0.5 这个上界足够宽。如果配平时报 “f(a) and f(b) must have different signs”,九成是 CT 给了负值或者 μ 写反了,先检查这两个入参。
桨盘非均匀入流用 Pitt-Peters 三状态动态入流模型补上,三个状态分别对应均匀、横向一阶、纵向一阶分量。过渡段短舱倾转过程中,桨盘迎角连续变化,动态入流能明显改善俯仰力矩的相位,稳态入流模型会把过渡段的俯仰响应算得过于超前。
2.3 机翼、机身、平尾的气动叠加与滑流修正
机翼在过渡段承担从 0 到接近全部升力的交接任务,麻烦在于旋翼滑流会扫过机翼内侧。滑流区的动压比自由来流高,局部升力增加,但同时带来低头力矩和额外的阻力。工程上通常把机翼分成滑流区和自由流区两段,分别查升力曲线,再按滑流面积比例加权。
| 部件 | 建模方式 | 过渡段主要影响 |
|---|---|---|
| 旋翼 | 叶素+动态入流 | 拉力矢量随短舱角旋转 |
| 机翼 | 分段升力线+滑流修正 | 升力交接、低头力矩 |
| 机身 | 静导数表 | 垂向阻力、俯仰阻尼 |
| 平尾 | 二维升力曲线+下洗修正 | 俯仰配平主力 |
| 垂尾 | 侧力导数 | 过渡段航向稳定性 |
滑流修正系数一般写成动压比的形式,klip = sqrt(1 + 2·T/(ρ·V²·A_slip)),其中 A_slip 是滑流扫过的机翼面积。V 很小时这个系数会爆掉,所以低速段要设下限,常见做法是把 V 钳到 5 m/s 以上再算。
2.4 六自由度刚体运动方程与状态方程组装
力与力矩汇总到机体轴后,刚体方程写成标准形式:
import numpy as np def rigid_body_deriv(state, F, M, I, m): """state = [u,v,w,p,q,r,phi,theta,psi] F, M 为机体轴合力与合力矩,I 为惯量张量对角项""" u, v, w, p, q, r, phi, theta, psi = state Ixx, Iyy, Izz = I du = F[0]/m + r*v - q*w dv = F[1]/m + p*w - r*u dw = F[2]/m + q*u - p*v dp = (M[0] - (Izz - Iyy)*q*r) / Ixx dq = (M[1] - (Ixx - Izz)*r*p) / Iyy dr = (M[2] - (Iyy - Ixx)*p*q) / Izz dphi = p + (q*np.sin(phi) + r*np.cos(phi)) * np.tan(theta) dtheta = q*np.cos(phi) - r*np.sin(phi) dpsi = (q*np.sin(phi) + r*np.cos(phi)) / np.cos(theta) return np.array([du, dv, dw, dp, dq, dr, dphi, dtheta, dpsi])注意dphi和dpsi里的 tanθ 和 cosθ。过渡段俯仰角如果超过 ±60°,这两个函数会把姿态微分放大到失真,纯欧拉角积分就不够用了,得换四元数。XV-15 正常过渡段的俯仰角一般在 ±15° 以内,欧拉角够用,但做大幅度机动仿真时一定要切四元数。
3. 过渡段配平:未知量怎么选、fsolve 怎么写、走廊边界在哪
配平是过渡段一切分析的地基。配平点找不到,后面的线性和响应都是空中楼阁。
3.1 配平残差与未知量的配对
配平的本质是解一组非线性代数方程,未知量个数必须等于残差个数。XV-15 过渡段常用的配对是六个未知量配六个残差:
| 未知量 | 作用方向 | 对应残差 |
|---|---|---|
| 总距 θ0 | 垂向力 | Fz |
| 纵向周期变距 θ1c | 轴向力/俯仰 | Fx |
| 横向周期变距 θ1s | 侧向力 | Fy |
| 俯仰角 θ | 俯仰力矩 | My |
| 滚转角 φ | 滚转力矩 | Mx |
| 升降舵 δe | 偏航/俯仰补充 | Mz |
短舱角 β_n 和飞行速度 V 是外部给定的调度参数,不参与配平求解。这个配对的好处是物理意义清楚:拉力大小由总距管,姿态由周期变距和舵面管。
3.2 用 scipy.optimize.fsolve 解过渡段配平点
import numpy as np from scipy.optimize import fsolve def aero_forces(theta0, theta1c, theta1s, theta, phi, de, V, beta_n): """简化版气动力汇总,返回机体轴合力与合力矩 真实模型里这里要调用旋翼、机翼、尾翼三套子函数""" # 旋翼拉力方向随短舱角旋转:beta_n=90° 时全垂直,0° 时全水平 T = 5.0e4 * theta0 # 总距到拉力的线性化近似 Tx = -T * np.cos(np.radians(beta_n)) # 前向分量 Tz = -T * np.sin(np.radians(beta_n)) # 垂向分量 # 机翼升力随速度平方增长,小迎角线性段 q_dyn = 0.5 * 1.225 * V**2 Lw = q_dyn * 15.7 * (0.09 * (theta + 0.02)) # 15.7 m^2 机翼面积量级 Dw = 0.02 * Lw Fx = Tx - Dw Fz = Tz - Lw + 5.0e3 My = -Tx * 1.2 + Lw * 0.4 * (1 + theta) - de * 2.0e4 Mx = -theta1s * 8.0e3 - phi * 6.0e3 Mz = de * 3.0e3 - theta1c * 4.0e3 Fy = theta1s * 6.0e3 - phi * 4.0e3 return np.array([Fx, Fy, Fz, Mx, My, Mz]) def trim(V, beta_n, x0=None): if x0 is None: x0 = [0.10, 0.0, 0.0, 0.0, 0.0, 0.0] fun = lambda x: aero_forces(x[0], x[1], x[2], x[3], x[4], x[5], V, beta_n) sol, info, ier, msg = fsolve(fun, x0, full_output=True, xtol=1e-10) return sol, ier, msg # 沿过渡段逐点配平,上一点解当前点初值(continuation) V_list = np.arange(30, 130, 10) beta_list = np.linspace(90, 0, len(V_list)) x = None for V, b in zip(V_list, beta_list): x, ier, msg = trim(V, b, x) print(f"V={V:5.1f} m/s beta={b:5.1f} deg " f"theta0={x[0]:6.3f} theta={np.degrees(x[3]):6.2f} deg ier={ier}")两个参数要说明。xtol=1e-10是为了让配平残差落到力矩量级以下,过渡段力矩对姿态角特别敏感,容差松了结果会来回跳。更关键的是x0的传递——这叫延拓法,用上一个速度点的解做初值。如果不做延拓,每个点都从零初值起解,ier会频繁返回 4(迭代不收敛),尤其是在 β_n 落在 30°~60° 这个交接区间时。
3.3 过渡走廊的三条硬约束
配平解存在不代表能飞。XV-15 的过渡走廊由三件事框住:
第一是俯仰姿态。过渡段姿态角如果超过 ±15°,机身阻力和飞行员视界都会出问题,一般把姿态作为不等式约束加进配平求解,越界就给个惩罚项把解推回来。
第二是功率。总距和转速决定了需用功率,过渡段某些速度点上需用功率会超过发动机可用功率,表现为配平旋翼转速掉转,仿真里会看到 Ω 的积分一路下滑。
第三是抖振边界。机翼在滑流和来流叠加处先失速,配平解在数学上存在,但气动导数已经失效。工程做法是配平完成后回查机翼局部迎角,滑流区局部迎角超过失速迎角就标记该配平点不可用,走廊边界由此确定。
4. 非线性仿真与操纵特性分析:从状态方程到模态
配平只给了一个瞬间,操纵特性要看的是整个过渡过程。
4.1 积分器选型:为什么不能一路用 RK4
过渡段的动力学是刚性的。旋翼入流状态的时间常数在 0.1 s 量级,机体长周期在 10 s 量级,两者差两个数量级。定步长 RK4 为了照顾快模态,步长必须压到 0.001 s,算一次 60 秒过渡要 6 万步,效率低且容易积累舍入误差。
scipy.integrate.solve_ivp里的LSODA或BDF是自适应刚性问题求解器,会自动在快慢模态之间切换步长。我的习惯是先用LSODA跑通,再用固定步长的 RK4 做一次交叉验证,如果两条曲线在力矩峰值处差异超过 5%,说明气动模型里某个导数有分段或插值不连续。
import numpy as np from scipy.integrate import solve_ivp def dynamics(t, state, m, I, beta0, beta_rate, V_cmd): """state = [u,v,w,p,q,r,phi,theta,psi,x,y,z] 短舱角按固定速率从 beta0 线性倾转""" beta_n = max(0.0, beta0 - beta_rate * t) # 度 u, v, w = state[0:3] # 简化控制器:按速度误差给总距,按姿态误差给周期变距 theta0 = np.clip(0.10 + 0.002 * (V_cmd - u), 0.0, 0.35) theta1c = np.clip(-0.05 * state[8], -0.1, 0.1) F, M = compute_forces(state, theta0, theta1c, beta_n) # 见 3.2 的力汇总 dstate = np.zeros(12) dstate[0:9] = rigid_body_deriv(state[0:9], F, M, I, m) return dstate I = (2.4e4, 1.6e5, 1.5e5) # 惯量张量对角项量级 kg·m^2 m = 5900.0 # 过渡段质量 kg sol = solve_ivp( dynamics, t_span=(0, 60), y0=y0, args=(m, I, 90.0, 1.5, 100.0), # 90° 起转,1.5°/s,目标 100 m/s method="LSODA", rtol=1e-6, atol=1e-8, dense_output=True ) print("姿态角范围 deg:", np.degrees(sol.y[8]).min(), np.degrees(sol.y[8]).max()) print("最大俯仰角速度 rad/s:", np.abs(sol.y[4]).max())beta_rate取 1.5°/s 是 XV-15 过渡段的常见量级,60 秒走完 90°。rtol和atol不能设得太松,atol放到 1e-5 以上时,低速段的速度积分会漂,过渡末端速度能差出好几米每秒。args里把短舱角初值和速率传进去,是为了不改方程就能扫不同调度律。
4.2 沿过渡轨迹逐点线性化
操纵特性分析需要的是小扰动模型。做法很简单:在配平点上对状态和控制各加一个很小的扰动,用中心差分求雅可比矩阵。
def linearize(trim_state, trim_input, V, beta_n, h=1e-5): n, m = 9, 3 # 9 个状态,3 个关键控制 A = np.zeros((n, n)) B = np.zeros((n, m)) for i in range(n): sp = trim_state.copy(); sp[i] += h sm = trim_state.copy(); sm[i] -= h A[:, i] = (f(sp, trim_input, V, beta_n) - f(sm, trim_input, V, beta_n)) / (2*h) for j in range(m): up = trim_input.copy(); up[j] += h um = trim_input.copy(); um[j] -= h B[:, j] = (f(trim_state, up, V, beta_n) - f(trim_state, um, V, beta_n)) / (2*h) return A, B A, B = linearize(x_trim, u_trim, V=80.0, beta_n=45.0) eig = np.linalg.eigvals(A) print("特征值实部排序:", np.sort(eig.real)[::-1][:4])h选 1e-5 是个平衡点。太大,非线性项污染导数;太小,浮点相减把有效位吃光。如果算出来的特征值实部大得离谱(比如几百),大概率是h取小了,改到 1e-4 再试。
过渡段最值得盯的特征值有两类:一类是俯仰短周期,随短舱角减小频率升高,从直升机模式的 0.3 Hz 附近涨到固定翼模式的 1 Hz 以上;另一类是入流的动态模态,实部很负,本身就是刚性的来源。真正的风险信号是出现实部为正的根,通常伴随滚转-偏航耦合,说明该配平点的航向稳定性裕度已经不够。
注意:特征值分析只在配平点附近有效。过渡过程是强参数时变系统,逐点特征值只能作为趋势判断,不能替代时域仿真。
5. 过渡段操纵分配与排错:几个能省下大量时间的技巧
5.1 短舱角速率和升降舵提前量的配合
过渡段最常见的失控不是发散,是姿态缓慢下沉。原因是机翼升力的建立比旋翼拉力的损失晚半个节拍。经验做法是让升降舵提前介入:在 β_n 降到 45° 之前,就让升降舵给出一个与 β_n 变化率成比例的前馈量,而不是等俯仰角误差出现再反馈。
def elevator_schedule(beta_n, beta_dot, theta_err): """前馈 + 反馈复合的升降舵指令""" k_ff = 0.030 # 前馈增益,每度短舱角 k_p = 1.2 # 姿态反馈增益 # 短舱角越小,机翼越接管升力,前馈量越要减小 ff = k_ff * (90.0 - beta_n) / 90.0 * (beta_dot / 1.5) return np.clip(ff + k_p * theta_err, -0.35, 0.35)k_ff按短舱角速率缩放,是因为倾转越快、升力交接越剧烈,需要的前馈越大。如果仿真里把beta_rate从 1.5°/s 提到 4°/s 而k_ff不动,俯仰角会在 β_n≈40° 附近出现一个明显的下沉尖峰,这正是前馈量不足的典型波形。
5.2 常见数值问题与排查方向
| 现象 | 大概率原因 | 处理方式 |
|---|---|---|
| fsolve 返回 ier=4 | 初值离解太远 | 用延拓法,从悬停配平起逐点推 |
| 姿态角在过渡末段发散 | 机翼升力曲线迎角越界 | 检查局部迎角,钳位后重算 |
| 速度积分单调漂移 | 积分器 atol 太松 | atol 收到 1e-8 以下 |
| 力矩出现周期性尖刺 | 气动表插值不连续 | 用样条替代线性插值 |
| 旋翼转速持续下滑 | 需用功率超过可用功率 | 降低倾转速率或推迟倾转起点 |
5.3 一个验证模型可信度的做法
配平解算完之后,把每个配平点的力和力矩代回残差函数,打印归一化残差(力除以重量,力矩除以参考力矩)。所有残差应该落在 1e-6 量级。如果有某个点是 1e-3 量级,说明该点的雅可比条件数很差,通常是两个配平未知量在互相打架——最常见的是 θ1c 和升降舵在低速段作用重叠。这时候改用带边界的least_squares并加上变量的物理上下限,比继续调fsolve的初值有效得多。
本文还有配套的精品资源,点击获取