火箭仿真系列-稳定性与动力学分析模块介绍

稳定性与动力学分析模块详解

稳定性与动力学分析模块是其最专业的组成部分，提供了从静态稳定性到复杂动态响应的全面分析能力。以下是该模块的深入解析：

一、模块架构概览

1.1 核心分析层次

StabilityAnalysis (稳定性分析)
├── StaticStability (静态稳定性)
│ ├── NeutralPoint (中性点)
│ ├── StaticMargin (静稳定裕度)
│ └── AerodynamicCenter (气动中心)
├── DynamicStability (动态稳定性)
│ ├── Eigenanalysis (特征值分析)
│ ├── ModeAnalysis (模态分析)
│ └── DampingAnalysis (阻尼分析)
├── FlightDynamics (飞行动力学)
│ ├── 6DOFDynamics (六自由度动力学)
│ ├── CoupledDynamics (耦合动力学)
│ └── FlexibleBodyDynamics (柔性体动力学)
└── Aeroelasticity (气动弹性)
├── FlutterAnalysis (颤振分析)
└── DivergenceAnalysis (发散分析)

二、静态稳定性分析

2.1 气动中心与中性点计算

class StaticStability:
"""静态稳定性分析"""

def __init__(self, rocket):
self.rocket = rocket
self.components = rocket.aerodynamic_surfaces

def compute_aerodynamic_center(self, mach, alpha=0.0):
"""计算气动中心位置"""

# 收集所有气动面的贡献
total_lift = 0.0
total_moment = 0.0
ac_positions = []

for comp in self.components:
# 计算每个组件的升力和力矩
lift = comp.compute_lift_coefficient(mach, alpha)
moment = comp.compute_moment_coefficient(mach, alpha)
position = comp.position # 从火箭尾部算起

total_lift += lift
total_moment += moment
ac_positions.append(position)

if total_lift != 0:
# 气动中心位置 (从火箭尾部算起)
# Σ(升力 * 位置) / Σ升力
ac_position = sum([lift * pos for lift, pos in
zip([c.compute_lift_coefficient(mach, alpha)
for c in self.components], ac_positions)]) / total_lift
else:
ac_position = 0.0

return {
'position': ac_position, # 从尾部算起 (m)
'total_lift_coefficient': total_lift,
'total_moment_coefficient': total_moment,
'mach': mach,
'alpha': alpha
}

def compute_neutral_point(self, mach_range=(0.1, 5.0), n_points=50):
"""计算中性点随马赫数的变化"""

mach_numbers = np.linspace(mach_range[0], mach_range[1], n_points)
neutral_points = []

for mach in mach_numbers:
# 计算升力曲线斜率
dCL_dalpha = self.compute_lift_curve_slope(mach)

# 计算俯仰力矩曲线斜率
dCm_dalpha = self.compute_moment_curve_slope(mach)

if dCL_dalpha != 0:
# 中性点位置: x_NP = x_ref – (dCm/dα) / (dCL/dα)
# 其中x_ref是参考点位置 (通常为箭体长度25%处)
neutral_point = self.rocket.reference_position – dCm_dalpha / dCL_dalpha
neutral_points.append(neutral_point)
else:
neutral_points.append(np.nan)

return {
'mach_numbers': mach_numbers,
'neutral_points': np.array(neutral_points)
}

def compute_static_margin(self, cg_position, mach, alpha=0.0):
"""计算静稳定裕度"""

# 获取气动中心位置
ac_data = self.compute_aerodynamic_center(mach, alpha)
ac_position = ac_data['position']

# 计算距离
distance_ac_cg = ac_position – cg_position

# 静稳定裕度 (以箭体直径为参考)
static_margin = distance_ac_cg / self.rocket.radius

# 静稳定裕度 (以箭体长度为参考)
static_margin_length = distance_ac_cg / self.rocket.length

stability = {
'static_margin_calibers': static_margin, # 倍直径
'static_margin_length': static_margin_length, # 倍长度
'distance_ac_cg': distance_ac_cg, # 绝对距离 (m)
'cg_position': cg_position, # 从尾部算起 (m)
'ac_position': ac_position, # 从尾部算起 (m)
'mach': mach,
'alpha': alpha,
'is_stable': static_margin > 1.0 # 通常要求 > 1倍直径
}

return stability

2.2 压力中心计算

class PressureCenterAnalysis:
"""压力中心分析"""

def compute_pressure_center(self, mach, alpha, method='vortex_lattice'):
"""计算压力中心"""

if method == 'vortex_lattice':
return self.vortex_lattice_method(mach, alpha)
elif method == 'panel_method':
return self.panel_method(mach, alpha)
elif method == 'cfd_interpolation':
return self.cfd_interpolation(mach, alpha)
else:
return self.empirical_method(mach, alpha)

def vortex_lattice_method(self, mach, alpha):
"""涡格法计算压力中心"""

# 离散火箭表面为涡格
panels = self.discretize_rocket()

total_normal_force = 0.0
total_moment = 0.0
cp_weighted_sum = 0.0

for panel in panels:
# 计算每个面板的法向力系数
cn_panel = panel.compute_normal_force_coefficient(mach, alpha)

# 计算面板中心位置
panel_center = panel.get_center()

# 累加
total_normal_force += cn_panel
cp_weighted_sum += cn_panel * panel_center

if total_normal_force != 0:
cp_position = cp_weighted_sum / total_normal_force
else:
cp_position = 0.0

return {
'position': cp_position,
'normal_force_coefficient': total_normal_force,
'method': 'vortex_lattice',
'n_panels': len(panels)
}

def empirical_method(self, mach, alpha):
"""经验公式计算压力中心 (适用于初步设计)"""

# Barrowman 方法 (用于细长体火箭)
cp_components = []

# 鼻锥贡献
if hasattr(self.rocket, 'nose'):
cp_nose = self.rocket.nose.get_cp_position(mach)
cn_nose = self.rocket.nose.get_normal_force_coefficient(mach, alpha)
cp_components.append((cp_nose, cn_nose))

# 鳍片贡献
if hasattr(self.rocket, 'fins'):
for fin in self.rocket.fins:
cp_fin = fin.get_cp_position(mach, alpha)
cn_fin = fin.get_normal_force_coefficient(mach, alpha)
cp_components.append((cp_fin, cn_fin))

# 箭体贡献 (忽略，因为通常很小)

# 计算加权平均
if cp_components:
total_cn = sum(cn for _, cn in cp_components)
if total_cn > 0:
cp_weighted = sum(cp * cn for cp, cn in cp_components) / total_cn
else:
cp_weighted = 0.0
else:
cp_weighted = 0.0

return {
'position': cp_weighted,
'normal_force_coefficient': total_cn,
'method': 'empirical',
'n_components': len(cp_components)
}

三、动态稳定性分析

3.1 特征值分析

class DynamicStability:
"""动态稳定性分析"""

def __init__(self, rocket, flight_conditions):
self.rocket = rocket
self.conditions = flight_conditions

# 质量特性
self.mass = rocket.mass
self.inertia = rocket.inertia
self.cg_position = rocket.center_of_mass

# 气动导数
self.aero_derivatives = {}

def compute_aerodynamic_derivatives(self, mach, alpha, beta):
"""计算气动导数"""

# 纵向导数
derivatives = {
# 力系数对攻角的导数
'C_L_alpha': self.compute_CL_alpha(mach, alpha), # 升力曲线斜率
'C_D_alpha': self.compute_CD_alpha(mach, alpha), # 阻力对攻角导数
'C_m_alpha': self.compute_Cm_alpha(mach, alpha), # 俯仰力矩曲线斜率

# 阻尼导数
'C_m_q': self.compute_Cmq(mach, alpha), # 俯仰阻尼
'C_L_q': self.compute_CLq(mach, alpha), # 升力对俯仰率导数

# 侧向导数
'C_Y_beta': self.compute_CY_beta(mach, beta), # 侧力对侧滑角导数
'C_l_beta': self.compute_Cl_beta(mach, beta), # 滚转力矩对侧滑角导数
'C_n_beta': self.compute_Cn_beta(mach, beta), # 偏航力矩对侧滑角导数

# 侧向阻尼导数
'C_l_p': self.compute_Clp(mach, beta), # 滚转阻尼
'C_n_r': self.compute_Cnr(mach, beta), # 偏航阻尼
'C_Y_r': self.compute_CYr(mach, beta), # 侧力对偏航率导数
}

self.aero_derivatives = derivatives
return derivatives

def compute_longitudinal_modes(self, velocity, density, mach, alpha):
"""计算纵向模态"""

# 获取气动导数
derivatives = self.compute_aerodynamic_derivatives(mach, alpha, 0)

# 参考值
S = self.rocket.reference_area
c = self.rocket.reference_length
q = 0.5 * density * velocity**2

# 质量参数
m = self.mass
Iyy = self.inertia[1] # 俯仰惯性矩

# 构建纵向状态矩阵
A_long = np.zeros((4, 4))

# 状态变量: [u, w, q, theta]
# u: 纵向速度扰动, w: 垂向速度扰动, q: 俯仰角速率, theta: 俯仰角

# 方程系数
Xu = q * S * (derivatives['C_D_alpha'] – 2 * self.rocket.CD) / (m * velocity)
Xw = q * S * derivatives['C_L_alpha'] / (m * velocity)
Zu = q * S * (derivatives['C_L_alpha'] + 2 * self.rocket.CL) / (m * velocity)
Zw = q * S * derivatives['C_D_alpha'] / (m * velocity)
Mu = q * S * c * derivatives['C_m_alpha'] / (Iyy * velocity)
Mw = q * S * c * derivatives['C_m_alpha'] / (Iyy * velocity)
Mq = q * S * c**2 * derivatives['C_m_q'] / (Iyy * 2 * velocity)

# 填充状态矩阵
A_long[0, 0] = Xu
A_long[0, 1] = Xw
A_long[0, 2] = 0
A_long[0, 3] = -9.81 * np.cos(self.conditions.theta)

A_long[1, 0] = Zu
A_long[1, 1] = Zw
A_long[1, 2] = velocity
A_long[1, 3] = -9.81 * np.sin(self.conditions.theta)

A_long[2, 0] = Mu
A_long[2, 1] = Mw
A_long[2, 2] = Mq
A_long[2, 3] = 0

A_long[3, 2] = 1
A_long[3, 3] = 0

# 特征值分析
eigenvalues, eigenvectors = np.linalg.eig(A_long)

# 分析模态
modes = self.analyze_longitudinal_modes(eigenvalues)

return {
'state_matrix': A_long,
'eigenvalues': eigenvalues,
'eigenvectors': eigenvectors,
'modes': modes
}

def analyze_longitudinal_modes(self, eigenvalues):
"""分析纵向模态"""

modes = {}

for i, eig in enumerate(eigenvalues):
freq_natural = np.abs(eig) # 自然频率
damping_ratio = -np.real(eig) / np.abs(eig) if np.abs(eig) > 0 else 0

# 识别模态类型
if np.imag(eig) == 0:
# 实特征值 -> 非周期模式
tau = -1 / np.real(eig) if np.real(eig) < 0 else np.inf
mode_type = 'aperiodic'
description = f"衰减时间常数: {tau:.3f}s"
else:
# 复特征值 -> 振荡模式
freq_damped = np.imag(eig) / (2 * np.pi)
period = 1 / freq_damped if freq_damped > 0 else np.inf

if freq_damped < 1.0:
mode_type = 'phugoid' # 长周期模态
description = f"频率: {freq_damped:.3f}Hz, 周期: {period:.1f}s"
else:
mode_type = 'short_period' # 短周期模态
description = f"频率: {freq_damped:.3f}Hz, 周期: {period:.2f}s"

modes[f'mode_{i}'] = {
'eigenvalue': eig,
'natural_frequency': freq_natural,
'damping_ratio': damping_ratio,
'type': mode_type,
'description': description,
'is_stable': np.real(eig) < 0
}

return modes

3.2 横向-方向稳定性

def compute_lateral_directional_modes(self, velocity, density, mach, beta):
"""计算横向-方向模态"""

# 获取气动导数
derivatives = self.compute_aerodynamic_derivatives(mach, 0, beta)

# 参考值
S = self.rocket.reference_area
b = self.rocket.span
c = self.rocket.reference_length
q = 0.5 * density * velocity**2

# 质量参数
m = self.mass
Ixx = self.inertia[0] # 滚转惯性矩
Izz = self.inertia[2] # 偏航惯性矩

# 构建横向-方向状态矩阵
A_lat = np.zeros((4, 4))

# 状态变量: [v, p, r, phi]
# v: 侧向速度, p: 滚转角速率, r: 偏航角速率, phi: 滚转角

# 方程系数
Yv = q * S * derivatives['C_Y_beta'] / (m * velocity)
Yp = q * S * b * derivatives['C_Y_r'] / (2 * m * velocity)
Yr = q * S * b * derivatives['C_Y_r'] / (2 * m * velocity)

Lv = q * S * b * derivatives['C_l_beta'] / (Ixx * velocity)
Lp = q * S * b**2 * derivatives['C_l_p'] / (2 * Ixx * velocity)
Lr = q * S * b**2 * derivatives['C_l_p'] / (2 * Ixx * velocity)

Nv = q * S * b * derivatives['C_n_beta'] / (Izz * velocity)
Np = q * S * b**2 * derivatives['C_n_r'] / (2 * Izz * velocity)
Nr = q * S * b**2 * derivatives['C_n_r'] / (2 * Izz * velocity)

# 填充状态矩阵
A_lat[0, 0] = Yv
A_lat[0, 1] = Yp
A_lat[0, 2] = -(velocity – Yr)
A_lat[0, 3] = 9.81 * np.cos(self.conditions.theta)

A_lat[1, 0] = Lv
A_lat[1, 1] = Lp
A_lat[1, 2] = Lr
A_lat[1, 3] = 0

A_lat[2, 0] = Nv
A_lat[2, 1] = Np
A_lat[2, 2] = Nr
A_lat[2, 3] = 0

A_lat[3, 1] = 1
A_lat[3, 2] = np.tan(self.conditions.theta)

# 特征值分析
eigenvalues, eigenvectors = np.linalg.eig(A_lat)

# 分析模态
modes = self.analyze_lateral_modes(eigenvalues)

return {
'state_matrix': A_lat,
'eigenvalues': eigenvalues,
'eigenvectors': eigenvectors,
'modes': modes
}

def analyze_lateral_modes(self, eigenvalues):
"""分析横向模态"""

modes = {}

for i, eig in enumerate(eigenvalues):
freq_natural = np.abs(eig)
damping_ratio = -np.real(eig) / np.abs(eig) if np.abs(eig) > 0 else 0

if np.imag(eig) == 0:
# 实特征值
if np.real(eig) < 0:
mode_type = 'roll_subsidence' # 滚转收敛
tau = -1 / np.real(eig)
description = f"滚转收敛, 时间常数: {tau:.3f}s"
else:
mode_type = 'spiral_divergence' # 螺旋发散
tau = 1 / np.real(eig)
description = f"螺旋发散, 倍幅时间: {tau:.3f}s"
else:
# 复特征值
freq_damped = np.imag(eig) / (2 * np.pi)
period = 1 / freq_damped

if damping_ratio < 0.1:
mode_type = 'dutch_roll' # 荷兰滚
description = f"荷兰滚, 频率: {freq_damped:.3f}Hz, 阻尼比: {damping_ratio:.3f}"
else:
mode_type = 'oscillatory'
description = f"振荡, 频率: {freq_damped:.3f}Hz, 阻尼比: {damping_ratio:.3f}"

modes[f'mode_{i}'] = {
'eigenvalue': eig,
'natural_frequency': freq_natural,
'damping_ratio': damping_ratio,
'type': mode_type,
'description': description,
'is_stable': np.real(eig) < 0
}

return modes

四、飞行动力学仿真

4.1 六自由度动力学

class SixDOFDynamics:
"""六自由度飞行动力学"""

def __init__(self, rocket, environment):
self.rocket = rocket
self.env = environment

# 状态变量
self.state = {
'position': np.zeros(3), # 位置 (ECI系)
'velocity': np.zeros(3), # 速度 (ECI系)
'quaternion': np.array([1, 0, 0, 0]), # 姿态四元数
'angular_velocity': np.zeros(3), # 角速度 (机体系)
'mass': rocket.total_mass
}

# 质量特性
self.mass = rocket.total_mass
self.inertia = rocket.inertia
self.cg_position = rocket.center_of_mass

# 气动数据库
self.aero_db = {}

def equations_of_motion(self, t, state):
"""六自由度运动方程"""

# 解包状态
pos = state[0:3] # 位置
vel = state[3:6] # 速度
quat = state[6:10] # 四元数
omega = state[10:13] # 角速度
mass = state[13] # 质量

# 获取环境条件
altitude = pos[2]
env_cond = self.env.get_conditions(altitude, t)

# 计算空气动力
aero_forces, aero_moments = self.compute_aerodynamics(vel, quat, omega, env_cond)

# 计算推力
thrust_force, thrust_moment = self.compute_thrust(t, quat)

# 计算重力
gravity_force = self.compute_gravity(pos, mass)

# 总外力
total_force = aero_forces + thrust_force + gravity_force

# 总力矩
total_moment = aero_moments + thrust_moment

# 计算加速度
acceleration = total_force / mass

# 四元数运动学方程
q0, q1, q2, q3 = quat
omega_matrix = np.array([
[0, -omega[0], -omega[1], -omega[2]],
[omega[0], 0, omega[2], -omega[1]],
[omega[1], -omega[2], 0, omega[0]],
[omega[2], omega[1], -omega[0], 0]
])
quat_dot = 0.5 * omega_matrix @ quat

# 欧拉动力学方程
inertia_inv = np.linalg.inv(self.inertia)
omega_dot = inertia_inv @ (total_moment – np.cross(omega, self.inertia @ omega))

# 质量变化
mass_dot = self.compute_mass_rate(t)

# 组装状态导数
state_dot = np.zeros(14)
state_dot[0:3] = vel
state_dot[3:6] = acceleration
state_dot[6:10] = quat_dot
state_dot[10:13] = omega_dot
state_dot[13] = mass_dot

return state_dot

def compute_aerodynamics(self, velocity, quaternion, omega, env_cond):
"""计算空气动力和力矩"""

# 转换到机体坐标系
DCM = self.quaternion_to_dcm(quaternion)
vel_body = DCM @ velocity

# 计算相对风速
wind_body = DCM @ np.array([env_cond['wind_u'], env_cond['wind_v'], env_cond['wind_w']])
vel_rel = vel_body – wind_body

# 计算气流参数
vel_mag = np.linalg.norm(vel_rel)
alpha, beta = self.compute_flow_angles(vel_rel)
mach = vel_mag / env_cond['speed_of_sound']
q_bar = 0.5 * env_cond['density'] * vel_mag**2

# 气动系数
CD = self.compute_drag_coefficient(mach, alpha)
CL = self.compute_lift_coefficient(mach, alpha)
CY = self.compute_side_force_coefficient(mach, beta)

Cm = self.compute_pitching_moment_coefficient(mach, alpha)
Cn = self.compute_yawing_moment_coefficient(mach, beta)
Cl = self.compute_rolling_moment_coefficient(mach, beta)

# 无量纲角速率
p = omega[0] * self.rocket.radius / (2 * vel_mag) if vel_mag > 0 else 0
q = omega[1] * self.rocket.length / (2 * vel_mag) if vel_mag > 0 else 0
r = omega[2] * self.rocket.radius / (2 * vel_mag) if vel_mag > 0 else 0

# 阻尼导数
CD_q = self.compute_CD_q(mach, alpha)
CL_q = self.compute_CL_q(mach, alpha)
CY_r = self.compute_CY_r(mach, beta)

Cm_q = self.compute_Cm_q(mach, alpha)
Cn_r = self.compute_Cn_r(mach, beta)
Cl_p = self.compute_Cl_p(mach, beta)

# 总系数
CD_total = CD + CD_q * q
CL_total = CL + CL_q * q
CY_total = CY + CY_r * r

Cm_total = Cm + Cm_q * q
Cn_total = Cn + Cn_r * r
Cl_total = Cl + Cl_p * p

# 气动力 (机体坐标系)
force_body = q_bar * self.rocket.reference_area * np.array([
-CD_total,
CY_total,
-CL_total
])

# 气动力矩 (机体坐标系)
moment_body = q_bar * self.rocket.reference_area * np.array([
Cl_total * self.rocket.radius,
Cm_total * self.rocket.length,
Cn_total * self.rocket.radius
])

return force_body, moment_body

def compute_flow_angles(self, velocity_body):
"""计算攻角和侧滑角"""

u, v, w = velocity_body

if u > 0:
alpha = np.arctan2(w, u)
beta = np.arcsin(v / np.linalg.norm(velocity_body))
else:
alpha = 0.0
beta = 0.0

return alpha, beta

4.2 耦合动力学分析

class CoupledDynamics:
"""耦合动力学分析"""

def __init__(self, rocket, environment):
self.rocket = rocket
self.env = environment

# 结构模态
self.structural_modes = self.compute_structural_modes()

# 气动弹性耦合
self.aeroelastic_coupling = {}

def compute_structural_modes(self):
"""计算结构模态"""

# 简化为欧拉-伯努利梁
E = rocket.youngs_modulus
I = rocket.area_moment_of_inertia
L = rocket.length
m = rocket.mass_per_length

# 计算固有频率
frequencies = []
mode_shapes = []

for n in range(1, 6): # 前5阶模态
# 悬臂梁固有频率
beta_L = [1.875, 4.694, 7.855, 10.996, 14.137][n-1]
omega_n = (beta_L**2 / L**2) * np.sqrt(E * I / m)
freq_n = omega_n / (2 * np.pi)

# 模态振型
x = np.linspace(0, L, 100)
phi = (np.cosh(beta_L * x / L) – np.cos(beta_L * x / L) –
(np.cosh(beta_L) + np.cos(beta_L)) /
(np.sinh(beta_L) + np.sin(beta_L)) *
(np.sinh(beta_L * x / L) – np.sin(beta_L * x / L)))

frequencies.append(freq_n)
mode_shapes.append(phi)

return {
'frequencies': np.array(frequencies),
'mode_shapes': np.array(mode_shapes),
'damping_ratios': np.array([0.01, 0.01, 0.01, 0.01, 0.01]) # 典型值
}

def coupled_aeroelastic_analysis(self, mach, dynamic_pressure):
"""耦合气动弹性分析"""

# 结构质量矩阵
M = self.compute_mass_matrix()

# 结构刚度矩阵
K = self.compute_stiffness_matrix()

# 结构阻尼矩阵
C = self.compute_damping_matrix()

# 气动刚度矩阵
K_a = self.compute_aerodynamic_stiffness(mach)

# 气动阻尼矩阵
C_a = self.compute_aerodynamic_damping(mach)

# 总矩阵
M_total = M
C_total = C + C_a
K_total = K + K_a

# 特征值问题
n = len(M)
A = np.zeros((2*n, 2*n))
A[:n, n:] = np.eye(n)
A[n:, :n] = -np.linalg.inv(M_total) @ K_total
A[n:, n:] = -np.linalg.inv(M_total) @ C_total

eigenvalues, eigenvectors = np.linalg.eig(A)

# 分析稳定性
stability = self.analyze_aeroelastic_stability(eigenvalues, dynamic_pressure)

return {
'eigenvalues': eigenvalues,
'eigenvectors': eigenvectors,
'stability': stability,
'matrices': {
'M': M, 'K': K, 'C': C,
'K_a': K_a, 'C_a': C_a,
'M_total': M_total, 'C_total': C_total, 'K_total': K_total
}
}

def analyze_aeroelastic_stability(self, eigenvalues, dynamic_pressure):
"""分析气动弹性稳定性"""

stability = {
'flutter_speed': None,
'divergence_speed': None,
'is_stable': True,
'critical_modes': []
}

for i, eig in enumerate(eigenvalues):
damping = -np.real(eig) / np.abs(eig) if np.abs(eig) > 0 else 0

if damping < 0:
# 负阻尼 -> 颤振
stability['is_stable'] = False
stability['critical_modes'].append({
'mode': i,
'damping': damping,
'frequency': np.abs(eig.imag) / (2 * np.pi),
'type': 'flutter',
'critical_q': dynamic_pressure
})

if stability['flutter_speed'] is None or dynamic_pressure < stability['flutter_speed']:
stability['flutter_speed'] = dynamic_pressure

elif np.real(eig) > 0 and eig.imag == 0:
# 正实部 -> 发散
stability['is_stable'] = False
stability['critical_modes'].append({
'mode': i,
'damping': damping,
'frequency': 0,
'type': 'divergence',
'critical_q': dynamic_pressure
})

if stability['divergence_speed'] is None or dynamic_pressure < stability['divergence_speed']:
stability['divergence_speed'] = dynamic_pressure

return stability

五、稳定性指标与判据

5.1 稳定性判据

class StabilityCriteria:
"""稳定性判据"""

def __init__(self, rocket):
self.rocket = rocket

def check_static_stability(self, cg_position, mach, alpha):
"""检查静态稳定性"""

criteria = {}

# 1. 静稳定裕度
static_margin = (self.rocket.neutral_point(mach) – cg_position) / self.rocket.radius
criteria['static_margin'] = {
'value': static_margin,
'requirement': '> 1.0',
'pass': static_margin > 1.0,
'severity': 'critical' if static_margin < 0.5 else 'warning' if static_margin < 1.0 else 'ok'
}

# 2. 俯仰力矩曲线斜率
Cm_alpha = self.rocket.compute_Cm_alpha(mach, alpha)
criteria['Cm_alpha'] = {
'value': Cm_alpha,
'requirement': '< 0',
'pass': Cm_alpha < 0,
'severity': 'critical' if Cm_alpha >= 0 else 'ok'
}

# 3. 航向静稳定性
Cn_beta = self.rocket.compute_Cn_beta(mach, 0)
criteria['Cn_beta'] = {
'value': Cn_beta,
'requirement': '> 0',
'pass': Cn_beta > 0,
'severity': 'critical' if Cn_beta <= 0 else 'ok'
}

# 4. 横向静稳定性
Cl_beta = self.rocket.compute_Cl_beta(mach, 0)
criteria['Cl_beta'] = {
'value': Cl_beta,
'requirement': '< 0',
'pass': Cl_beta < 0,
'severity': 'warning' if Cl_beta >= 0 else 'ok'
}

return criteria

def check_dynamic_stability(self, eigenvalues):
"""检查动态稳定性"""

criteria = {}

for i, eig in enumerate(eigenvalues):
damping_ratio = -np.real(eig) / np.abs(eig) if np.abs(eig) > 0 else 0
frequency = np.abs(eig.imag) / (2 * np.pi)

if eig.imag == 0:
# 非周期模式
if np.real(eig) < 0:
status = 'stable'
requirement = 'real part < 0'
else:
status = 'unstable'
requirement = 'real part < 0'
else:
# 振荡模式
if damping_ratio > 0:
status = 'stable'
elif damping_ratio == 0:
status = 'neutrally stable'
else:
status = 'unstable'

requirement = f'damping ratio > 0, current: {damping_ratio:.3f}'

criteria[f'mode_{i}'] = {
'eigenvalue': eig,
'damping_ratio': damping_ratio,
'frequency_hz': frequency,
'status': status,
'requirement': requirement,
'pass': status in ['stable', 'neutrally stable']
}

return criteria

def check_aeroelastic_stability(self, flutter_speed, divergence_speed, max_q):
"""检查气动弹性稳定性"""

criteria = {}

# 颤振边界
if flutter_speed is not None:
flutter_margin = (max_q – flutter_speed) / max_q
criteria['flutter'] = {
'flutter_speed': flutter_speed,
'max_q': max_q,
'margin': flutter_margin,
'requirement': 'margin > 0.2',
'pass': flutter_margin > 0.2,
'severity': 'critical' if flutter_margin <= 0 else 'warning' if flutter_margin <= 0.2 else 'ok'
}

# 发散边界
if divergence_speed is not None:
divergence_margin = (max_q – divergence_speed) / max_q
criteria['divergence'] = {
'divergence_speed': divergence_speed,
'max_q': max_q,
'margin': divergence_margin,
'requirement': 'margin > 0.2',
'pass': divergence_margin > 0.2,
'severity': 'critical' if divergence_margin <= 0 else 'warning' if divergence_margin <= 0.2 else 'ok'
}

return criteria

六、高级分析工具

6.1 稳定性边界分析

class StabilityBoundaryAnalysis:
"""稳定性边界分析"""

def __init__(self, rocket, flight_envelope):
self.rocket = rocket
self.envelope = flight_envelope

def compute_stability_map(self, parameter_space):
"""计算稳定性地图"""

# 参数空间: 马赫数, 动压, 攻角, 质心位置等
mach_range = parameter_space.get('mach', (0.1, 5.0))
q_range = parameter_space.get('dynamic_pressure', (0, 100e3))
cg_range = parameter_space.get('cg_position', (0.3, 0.7))

n_points = 50
mach_points = np.linspace(mach_range[0], mach_range[1], n_points)
q_points = np.linspace(q_range[0], q_range[1], n_points)

stability_map = np.zeros((n_points, n_points, 4)) # 静态, 动态, 颤振, 发散

for i, mach in enumerate(mach_points):
for j, q in enumerate(q_points):
# 计算当前条件下的稳定性
stability = self.analyze_point(mach, q)

# 存储结果
stability_map[i, j, 0] = stability['static']
stability_map[i, j, 1] = stability['dynamic']
stability_map[i, j, 2] = stability['flutter']
stability_map[i, j, 3] = stability['divergence']

return {
'mach_numbers': mach_points,
'dynamic_pressures': q_points,
'stability_map': stability_map,
'parameter_space': parameter_space
}

def analyze_point(self, mach, dynamic_pressure):
"""分析单个点的稳定性"""

# 计算攻角 (假设平衡攻角)
alpha = self.compute_trim_alpha(mach, dynamic_pressure)

# 静态稳定性
cg = self.rocket.center_of_mass
static_margin = (self.rocket.neutral_point(mach) – cg) / self.rocket.radius
static_stable = static_margin > 1.0

# 动态稳定性
eig_analysis = self.rocket.compute_eigenvalues(mach, alpha, dynamic_pressure)
dynamic_stable = all(np.real(eig) < 0 for eig in eig_analysis['eigenvalues'])

# 气动弹性稳定性
aeroelastic = self.rocket.aeroelastic_analysis(mach, dynamic_pressure)
flutter_stable = aeroelastic['flutter_speed'] > dynamic_pressure * 1.2
divergence_stable = aeroelastic['divergence_speed'] > dynamic_pressure * 1.2

return {
'static': static_stable,
'dynamic': dynamic_stable,
'flutter': flutter_stable,
'divergence': divergence_stable,
'static_margin': static_margin,
'eigenvalues': eig_analysis['eigenvalues'],
'flutter_speed': aeroelastic['flutter_speed'],
'divergence_speed': aeroelastic['divergence_speed']
}

6.2 蒙特卡洛稳定性分析

class MonteCarloStability:
"""蒙特卡洛稳定性分析"""

def __init__(self, rocket, n_samples=1000):
self.rocket = rocket
self.n_samples = n_samples

# 不确定性定义
self.uncertainties = {
'cg_position': {'type': 'normal', 'mean': rocket.cg_nominal, 'std': 0.01},
'mass': {'type': 'normal', 'mean': rocket.mass_nominal, 'std': 0.02},
'inertia': {'type': 'normal', 'mean': 1.0, 'std': 0.05},
'Cp_position': {'type': 'normal', 'mean': 1.0, 'std': 0.03},
'aero_coefficients': {'type': 'uniform', 'low': 0.9, 'high': 1.1}
}

def sample_parameters(self):
"""采样不确定参数"""

samples = {}

for name, dist in self.uncertainties.items():
if dist['type'] == 'normal':
samples[name] = np.random.normal(dist['mean'], dist['std'])
elif dist['type'] == 'uniform':
samples[name] = np.random.uniform(dist['low'], dist['high'])

return samples

def run_stability_analysis(self, flight_conditions):
"""运行稳定性分析"""

results = []

for i in range(self.n_samples):
# 采样参数
params = self.sample_parameters()

# 创建带不确定性的火箭
rocket_sample = self.create_rocket_sample(params)

# 分析稳定性
stability = self.analyze_stability(rocket_sample, flight_conditions)

results.append({
'sample': i,
'parameters': params,
'stability': stability,
'is_stable': stability['overall_stable']
})

# 统计分析
stats = self.compute_statistics(results)

return {
'results': results,
'statistics': stats,
'reliability': stats['reliability']
}

def compute_statistics(self, results):
"""计算统计量"""

n_stable = sum(1 for r in results if r['is_stable'])
reliability = n_stable / len(results)

# 收集关键指标
static_margins = [r['stability']['static_margin'] for r in results]
damping_ratios = []
for r in results:
for mode in r['stability']['modes'].values():
if 'damping_ratio' in mode:
damping_ratios.append(mode['damping_ratio'])

stats = {
'reliability': reliability,
'static_margin': {
'mean': np.mean(static_margins),
'std': np.std(static_margins),
'min': np.min(static_margins),
'max': np.max(static_margins),
'percentile_5': np.percentile(static_margins, 5),
'percentile_95': np.percentile(static_margins, 95)
},
'damping_ratio': {
'mean': np.mean(damping_ratios) if damping_ratios else 0,
'min': np.min(damping_ratios) if damping_ratios else 0
}
}

return stats

七、使用示例

7.1 完整稳定性分析流程

# 1. 创建火箭和环境
from rocketpy import Rocket, Environment, Flight
from rocketpy.stability import StaticStability, DynamicStability, StabilityCriteria
import numpy as np

# 创建火箭
calisto = Rocket(
radius=0.0635,
mass=14.426,
inertia=(6.321, 6.321, 0.034),
center_of_mass=1.5,
reference_area=np.pi * 0.0635**2,
reference_length=2.0
)

# 添加气动面
calisto.add_aerodynamic_surface(nose_cone, position=1.7)
calisto.add_aerodynamic_surface(fins, position=0.2)

# 2. 静态稳定性分析
static_analysis = StaticStability(calisto)

# 分析不同马赫数下的静稳定裕度
mach_numbers = np.linspace(0.1, 3.0, 20)
static_results = []

for mach in mach_numbers:
result = static_analysis.compute_static_margin(
cg_position=calisto.center_of_mass,
mach=mach,
alpha=0.0
)
static_results.append(result)

# 3. 动态稳定性分析
flight_conditions = {
'velocity': 300, # m/s
'density': 1.225, # kg/m³
'mach': 0.88,
'alpha': 0.0,
'theta': np.radians(85) # 发射角度
}

dynamic_analysis = DynamicStability(calisto, flight_conditions)

# 纵向模态分析
longitudinal = dynamic_analysis.compute_longitudinal_modes(
velocity=flight_conditions['velocity'],
density=flight_conditions['density'],
mach=flight_conditions['mach'],
alpha=flight_conditions['alpha']
)

# 横向-方向模态分析
lateral = dynamic_analysis.compute_lateral_directional_modes(
velocity=flight_conditions['velocity'],
density=flight_conditions['density'],
mach=flight_conditions['mach'],
beta=0.0
)

# 4. 稳定性判据检查
criteria_checker = StabilityCriteria(calisto)

# 静态稳定性判据
static_criteria = criteria_checker.check_static_stability(
cg_position=calisto.center_of_mass,
mach=1.0,
alpha=0.0
)

# 动态稳定性判据
dynamic_criteria = criteria_checker.check_dynamic_stability(
longitudinal['eigenvalues']
)

# 5. 稳定性边界分析
stability_map = StabilityBoundaryAnalysis(
rocket=calisto,
flight_envelope={'mach_max': 3.0, 'q_max': 50e3}
)

parameter_space = {
'mach': (0.1, 3.0),
'dynamic_pressure': (0, 50e3),
'cg_position': (1.3, 1.7)
}

stability_boundary = stability_map.compute_stability_map(parameter_space)

# 6. 蒙特卡洛分析
mc_analysis = MonteCarloStability(calisto, n_samples=1000)
mc_results = mc_analysis.run_stability_analysis(flight_conditions)

# 7. 输出结果
print("="*60)
print("稳定性分析报告")
print("="*60)

print(f"\\n1. 静态稳定性:")
for mach, result in zip(mach_numbers, static_results):
if result['is_stable']:
status = "✓ 稳定"
else:
status = "✗ 不稳定"
print(f" 马赫数 {mach:.2f}: 静稳定裕度 {result['static_margin_calibers']:.2f}倍直径 – {status}")

print(f"\\n2. 动态稳定性:")
print(" 纵向模态:")
for mode_name, mode_info in longitudinal['modes'].items():
stability = "稳定" if mode_info['is_stable'] else "不稳定"
print(f" {mode_name}: {mode_info['type']}, 阻尼比 {mode_info['damping_ratio']:.3f}, {stability}")

print("\\n 横向-方向模态:")
for mode_name, mode_info in lateral['modes'].items():
stability = "稳定" if mode_info['is_stable'] else "不稳定"
print(f" {mode_name}: {mode_info['type']}, 阻尼比 {mode_info['damping_ratio']:.3f}, {stability}")

print(f"\\n3. 稳定性判据:")
print(" 静态稳定性:")
for name, criterion in static_criteria.items():
status = "通过" if criterion['pass'] else "失败"
print(f" {name}: {criterion['value']:.3f} ({criterion['requirement']}) – {status}")

print(f"\\n4. 蒙特卡洛分析结果:")
print(f" 可靠性: {mc_results['reliability']*100:.1f}%")
print(f" 静稳定裕度: {mc_results['statistics']['static_margin']['mean']:.2f} ± {mc_results['statistics']['static_margin']['std']:.2f}倍直径")
print(f" 5%分位静稳定裕度: {mc_results['statistics']['static_margin']['percentile_5']:.2f}倍直径")

# 8. 可视化
import matplotlib.pyplot as plt

# 静稳定裕度随马赫数变化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 子图1: 静稳定裕度
ax1 = axes[0, 0]
static_margins = [r['static_margin_calibers'] for r in static_results]
ax1.plot(mach_numbers, static_margins, 'b-', linewidth=2)
ax1.axhline(y=1.0, color='r', linestyle='–', label='最小稳定裕度 (1.0)')
ax1.set_xlabel('马赫数')
ax1.set_ylabel('静稳定裕度 (倍直径)')
ax1.set_title('静稳定裕度 vs 马赫数')
ax1.grid(True)
ax1.legend()

# 子图2: 根轨迹
ax2 = axes[0, 1]
eigenvalues = longitudinal['eigenvalues']
real_parts = [eig.real for eig in eigenvalues]
imag_parts = [eig.imag for eig in eigenvalues]
ax2.scatter(real_parts, imag_parts, c='b', s=100)
ax2.axvline(x=0, color='r', linestyle='–')
ax2.set_xlabel('实部')
ax2.set_ylabel('虚部')
ax2.set_title('特征值根轨迹')
ax2.grid(True)

# 子图3: 稳定性边界
ax3 = axes[1, 0]
X, Y = np.meshgrid(stability_boundary['mach_numbers'],
stability_boundary['dynamic_pressures']/1000)
Z = stability_boundary['stability_map'][:, :, 0] # 静态稳定性
contour = ax3.contourf(X, Y, Z.T, levels=[-0.5, 0.5, 1.5],
colors=['red', 'green'], alpha=0.5)
ax3.set_xlabel('马赫数')
ax3.set_ylabel('动压 (kPa)')
ax3.set_title('稳定性边界')
ax3.grid(True)

# 子图4: 蒙特卡洛结果
ax4 = axes[1, 1]
static_margins_mc = [r['stability']['static_margin'] for r in mc_results['results']]
ax4.hist(static_margins_mc, bins=30, edgecolor='black', alpha=0.7)
ax4.axvline(x=1.0, color='r', linestyle='–', linewidth=2, label='稳定边界')
ax4.set_xlabel('静稳定裕度 (倍直径)')
ax4.set_ylabel('频数')
ax4.set_title('蒙特卡洛分析 – 静稳定裕度分布')
ax4.legend()
ax4.grid(True)

plt.tight_layout()
plt.show()

八、验证与测试

8.1 稳定性验证

def validate_stability_models():
"""验证稳定性模型"""

validation_cases = [
{
'name': 'Standard Rocket',
'parameters': {
'mass': 10.0,
'length': 1.5,
'radius': 0.05,
'cg_position': 0.75,
'cp_position': 1.0
},
'expected': {
'static_margin': 5.0, # (1.0-0.75)/0.05
'is_stable': True
}
},
{
'name': 'Unstable Rocket',
'parameters': {
'mass': 10.0,
'length': 1.5,
'radius': 0.05,
'cg_position': 1.1,
'cp_position': 1.0
},
'expected': {
'static_margin': -2.0, # (1.0-1.1)/0.05
'is_stable': False
}
}
]

results = []

for case in validation_cases:
# 创建测试火箭
rocket = Rocket(**case['parameters'])

# 运行稳定性分析
static_analysis = StaticStability(rocket)
result = static_analysis.compute_static_margin(
cg_position=rocket.center_of_mass,
mach=0.5,
alpha=0.0
)

# 验证
static_margin_error = abs(result['static_margin_calibers'] –
case['expected']['static_margin'])
stability_match = (result['is_stable'] == case['expected']['is_stable'])

results.append({
'case': case['name'],
'expected_static_margin': case['expected']['static_margin'],
'computed_static_margin': result['static_margin_calibers'],
'static_margin_error': static_margin_error,
'stability_match': stability_match,
'pass': static_margin_error < 0.01 and stability_match
})

return results

总结

稳定性与动力学分析模块提供了：

 能够为火箭设计提供从初步稳定性评估到详细动态响应的全面分析，确保设计的火箭在各种飞行条件下都能保持稳定可控。无论您是业余火箭爱好者还是专业航天工程师，这个模块都能提供可靠的稳定性分析工具。