MATLAB模型开发：从特征值提取到状态空间系统构建

# MATLAB模型开发：从特征值提取到状态空间系统构建 ## 1. 数据准备与模型开发启动 首先，从 `actrl.eig` 文件中提取特征值以及 UX 和 UY 特征向量条目，并将其存储在 MATLAB 的 `.mat` 文件 `actrl_eig.mat` 中。完成这一步后，就可以将 ANSYS 结果读入 MATLAB，开始开发简化模型。 ## 2. MATLAB 代码 `act8.m` 详细解析 ### 2.1 代码概述 代码的执行流程如下： 1. 读取 `actrl_eig.mat` 文件中的 ANSYS 模型特征值和特征向量结果，涵盖所有 50 个模式。 2. 依据图 17.12 中的角度，定义径向和周向的 VCM 力分量。 3. 提示用户选择是否对所有模式使用相同的阻尼值（均匀阻尼），或者每个模式使用不同的值（非均匀阻尼）。 4. 根据用户选择计算相应的增益（直流增益或峰值增益），并对增益进行排序和绘图，以确定需要保留的重要模式。 5. 构建两个状态空间系统，一个包含所有 50 个模式，另一个包含排序后简化的模式数量。 6. 绘制 50 模式响应图，可选择头部 0 或头部 1，并叠加各个模式的贡献。 7. 提示用户输入采样频率（默认 20 kHz），使用 MATLAB 的 “c2d” 命令创建原始连续系统的离散模型，并绘制离散频率响应图，与原始连续频率响应图进行叠加。 8. 使用简化、排序后的模式计算频率响应，截断不太重要的模式，并使用 “modred” “mdc” 选项。 ### 2.2 输入与自由度定义 代码的第一部分从 `actrl_eig.mat` 文件中读取特征值/特征向量数据，并明确所使用的自由度。原始 ANSYS 模型约有 21000 个自由度，通过仅定义所需频率响应的自由度，可将 MATLAB 模型所需的自由度数量减少到 12 个，具体如下表所示： | 自由度编号 | 节点编号 | 方向 | 位置 | | --- | --- | --- | --- | | 1 | 22 | ux - 径向 | 顶部磁头间隙 | | 2 | 10022 | ux - 径向 | 底部磁头间隙 | | 3 | 24061 | ux - 径向 | 线圈 | | 4 | 24066 | ux - 径向 | 线圈 | | 5 | 24082 | ux - 径向 | 线圈 | | 6 | 24087 | ux - 径向 | 线圈 | | 7 | 22 | uy - 周向 | 顶部磁头间隙 | | 8 | 10022 | uy - 周向 | 底部磁头间隙 | | 9 | 24061 | uy - 周向 | 线圈 | | 10 | 24066 | uy - 周向 | 线圈 | | 11 | 24082 | uy - 周向 | 线圈 | | 12 | 24087 | uy - 周向 | 线圈 | 以下是相关代码： ```matlab clear all; hold off; clf; % load the Block Lanczos .mat file actrl_eig.mat, containing evr - the modal matrix, % freqvec - the frequency vector and node_numbers - the vector of node numbers % for the modal matrix load actrl_eig; [numdof,num_modes_total] = size(evr); freqvec(1) = 0; % set frequency of rigid body mode to zero xn = evr; ``` ### 2.3 强迫函数定义与直流增益计算 生成用于增益计算的特征值平方向量（单位为 rad/sec）。如同之前章节讨论的刚体模式直流增益计算，这里再次使用频率响应计算中定义的最低频率来计算刚体模式的低频增益。 强迫函数的定义步骤如下： 1. 以图 17.12 为参考，定义四个线圈节点的强迫函数分量。在每个线圈施加单位力，并均匀分配到四个节点上。 2. 将每个线圈节点的力分解为径向和周向（x 和 y）分量。 3. 定义每个线圈节点在物理坐标下的线圈力，同时将磁头节点（自由度 1、2、7 和 8）的 ux 和 uy 力项设为零。 该模型属于 “SI” 或单输入模型，因为相同的力施加到所有四个线圈节点，输入矩阵 “b” 仅需一个单列向量。 直流增益和峰值增益的计算公式如下： \[ \frac{F}{z_{ji}z_{nki}\omega_{ki}^2} \] \[ \frac{F}{2\zeta_i z_{ji}z_{ki}} \] 其中，\(z_{nji}z_{nki}\) 为留数，是第 i 个特征向量的第 j 行（输出）和第 k 行（施加力）项的乘积除以第 i 个模式特征值的平方，\(\zeta_i\) 为第 i 个模式的阻尼。 以下是相关代码： ```matlab omega2 = (2*pi*freqvec)'.^2; % convert to radians and square freqlo = 501; freqhi = 25000; flo=log10(freqlo) ; fhi=log10(freqhi) ; f=logspace(flo,fhi,300) ; frad=f*2*pi ; n24061fx = 0.25*sin(9.1148*pi/180); n24061fy = 0.25*cos(9.1148*pi/180); n24066fx = 0.25*sin(15.1657*pi/180); n24066fy = 0.25*cos(15.1657*pi/180); n24082fx = -0.25*sin(15.1657*pi/180); n24082fy = 0.25*cos(15.1657*pi/180); n24087fx = -0.25*sin(9.1148*pi/180); n24087fy = 0.25*cos(9.1148*pi/180); f_physical = [ 0 0 n24061fx n24066fx n24082fx n24087fx 0 0 n24061fy n24066fy n24082fy n24087fy ]; force = f_physical'*xn; head = input('enter "0" default for head 0 or "1" for head 1 ... '); if isempty(head) head = 0; end zeta_type = input('enter "1" to read in damping vector (zetain.m) … or "enter" for uniform damping ... '); if (isempty(zeta_type)) zeta_type = 0; zeta_uniform = input('enter value for uniform damping, … .005 is 0.5% of critical (default) ... '); if (isempty(zeta_uniform)) zeta_uniform = 0.005; end zeta_unsort = zeta_uniform*ones(num_modes_total,1); gainstr = 'dc gain'; else zetain; % read in zeta_unsort damping vector from zetain.m file gainstr = 'peak gain'; end if length(zeta_unsort) ~= num_modes_total error(['error - zetain vector has ',num2str(length(zeta_unsort)), … ' entries instead of ',num2str(num_modes_total)]); end if zeta_type == 0 % dc gain gain_h0 = abs([force(1)*xn(8,1)/frad(1) ... force(2:num_modes_total).*xn(8,2:num_modes_total) ... ./omega2(2:num_modes_total)]); gain_h1 = abs([force(1)*xn(7,1)/frad(1) ... force(2:num_modes_total).*xn(7,2:num_modes_total) ... ./omega2(2:num_modes_total)]); elseif zeta_type == 1 % peak gain gain_h0 = abs([force(1)*xn(8,1)/frad(1) ... force(2:num_modes_total).*xn(8,2:num_modes_total) ... ./((2 ```