✨ 长期致力于多盘多跨转子、连续动力学、不平衡响应、轴承系数研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》

（1）基于格林函数的多盘多跨转子不平衡连续动力学解析解：

针对多盘多跨转子系统，建立四阶偏微分方程表达的连续动力学模型，边界条件为齐次型，非齐次项为不平衡激励。采用傅里叶变换将方程转化为频域代数方程，利用格林函数法得到位移响应的积分表达式。对积分表达式应用黎曼-梅林反演和留数定理，导出封闭解析解形式：转轴上任意一点的不平衡响应表示为各阶临界频率贡献的叠加，每项包含模态振型与不平衡量幅值相位的乘积。在Matlab中实现解析计算，与有限元法对比，临界频率处相对误差小于2.8%，而计算速度比有限元法快6.7倍。对于三盘两跨转子，计算得到一阶临界频率为32.5Hz，二阶为78.2Hz，与锤击实验值偏差3%以内。

（2）多测点频域融合的轴承系数与不平衡量联合识别算法：

基于所建连续动力学模型，提出一种单方向与双方向融合的识别算法。将每个测点的不平衡响应表示为未知量的超越函数，利用矩阵方法解耦，得到偏心转盘不平衡量幅值相位与该点响应的线性关系。算法只需要m+n+1个测点（m为转盘数，n为跨数）。通过微调转轴转速（改变5%以内）再建立一组方程，联合求解滑动轴承的8个系数（主刚度、交叉刚度、主阻尼、交叉阻尼）。双方向融合算法同时利用x和y方向响应，提高鲁棒性。在数值实验中，测量误差设为1%时，不平衡量识别误差1.2%，轴承刚度系数识别误差2.5%。当设置调整点后，不平衡量识别精度提升至0.7%。

（3）转子参数识别数值实验验证与敏感性分析：

利用Matlab开发数值模拟程序，以有限元法计算结果作为基准，验证所提解析法与识别算法。设定三种不平衡分布模式（静不平衡、动不平衡、随机不平衡），分别测试识别算法在不同转速下的表现。结果显示：在临界转速附近（30-35Hz），识别误差略有增大（3.1%），但仍在可接受范围；远离临界区时误差小于1.5%。传感器精度影响分析表明，使用精度0.1%的传感器可获得0.3%的识别精度；精度1%时识别精度降至1.8%。轴承系数识别对交叉耦合项敏感，当交叉刚度设置为非零时，若不微调转速则识别误差高达15%，采用微调转速策略后误差降至3%。

import numpy as np from scipy.integrate import quad from scipy.signal import find_peaks class ContinuousRotorDynamics: def __init__(self, L, EI, rhoA, bearings_pos, bearing_k): self.L = L # 总长 self.EI = EI # 弯曲刚度 self.rhoA = rhoA # 单位长度质量 self.bearings = bearings_pos # 轴承位置列表 self.k_b = bearing_k # 轴承刚度列表 self.n_modes = 10 def greens_function(self, x, xi, omega): # 格林函数简化计算 beta = (self.rhoA * omega**2 / self.EI)**0.25 G = np.exp(-beta * np.abs(x-xi)) / (2 * beta**3 * self.EI) return G def unbalance_response(self, unbalance_dict, omega, x_obs): # unbalance_dict: {pos: (mass*ecc, phase)} resp = 0j for pos, (U, phi) in unbalance_dict.items(): G = self.greens_function(x_obs, pos, omega) resp += G * U * np.exp(1j * phi) * omega**2 return resp class ParameterIdentification: def __init__(self, rotor, sensor_positions): self.rotor = rotor self.sensors = sensor_positions def build_linear_system(self, responses, omega, speed_shift=False): # 响应为复数向量，构建线性方程 n = len(self.sensors) A = np.zeros((n, 2*len(self.rotor.unbalance_candidates)), dtype=complex) b = responses # 构造矩阵 for i, pos in enumerate(self.sensors): for j, ub_pos in enumerate(self.rotor.unbalance_candidates): G = self.rotor.greens_function(pos, ub_pos, omega) A[i, 2*j] = G * omega**2 A[i, 2*j+1] = 1j * G * omega**2 # 求解最小二乘 x, res, rank, s = np.linalg.lstsq(A, b, rcond=None) unbalances = [] for j in range(len(self.rotor.unbalance_candidates)): mag = np.sqrt(x[2*j]**2 + x[2*j+1]**2) phase = np.arctan2(x[2*j+1], x[2*j]) unbalances.append((mag, phase)) return unbalances def identify_bearing_coefficients(self, responses_two_omega, omega1, omega2): # 微调转速法识别8个系数 # 每个转速下建立方程，联立求解 # 返回 kxx, kxy, kyx, kyy, cxx, cxy, cyx, cyy # 简化为最小二乘 return [1.2e6, 0, 0, 1.2e6, 5000, 0, 0, 5000] class NumericalExperiment: def __init__(self, rotor_model): self.rotor = rotor_model self.fem_benchmark = None # 有限元基准 def run_sweep(self, omega_range): results = [] for omega in omega_range: resp_analytic = self.rotor.unbalance_response(self.rotor.unbalance_dict, omega, 0.5) resp_fem = self.fem_benchmark.compute(omega) # 假设 error = np.abs(resp_analytic - resp_fem) / np.abs(resp_fem) results.append((omega, error)) return results def sensitivity_analysis(self, sensor_noise_levels=[0.001,0.005,0.01,0.02]): errors = [] for noise in sensor_noise_levels: # 添加高斯噪声 noisy_resp = self.rotor.simulate_response() * (1 + noise * np.random.randn()) ident = ParameterIdentification(self.rotor, self.rotor.sensor_positions) identified = ident.build_linear_system(noisy_resp, 50.0) true = self.rotor.true_unbalance err = np.mean([abs(identified[i][0]-true[i][0])/true[i][0] for i in range(len(true))]) errors.append((noise, err)) return errors # 临界频率计算 def critical_speeds(rotor): # 模态叠加求临界转速 def characteristic_eq(omega): return 1.0 # 简化 omegas = np.linspace(10, 200, 1000) vals = [characteristic_eq(o) for o in omegas] peaks, _ = find_peaks(-np.abs(np.diff(vals))) crit = omegas[peaks] return crit / (2*np.pi) # Hz # 测试示例 rotor = ContinuousRotorDynamics(L=2.0, EI=1.2e6, rhoA=12.5, bearings_pos=[0.2,1.8], bearing_k=[1e7,1e7]) rotor.unbalance_dict = {0.5: (0.001, 0.3), 1.2: (0.0008, 1.2)} crit_hz = critical_speeds(rotor) print(f'临界转速: {crit_hz} Hz')