自由边圆板固有频率与模态形状系数的MATLAB解析计算

简介：面向结构动力学与振动分析领域研究者的自由边圆板模态计算工具，提供基于Python实现的脚本代码。该资源依据Zagrai与Donskoy提出的弹性边缘支撑均匀圆板分析方法，通过求解特征方程得到无量纲系数lambda_mn与模态形状系数C_mn，进而计算自由边缘圆板的固有频率和振型。代码已通过Itao、Amabili等文献实验数据验证，适合用于教学演示或科研辅助计算。包体共4个文件，包括Python主程序circular_plate_free_edge.py、两组不同泊松比（v=0.33和v=0.35）的lambda计算结果dat数据，以及1份README说明文档，整体仅54KB，轻量易读。目前已有1314人学习下载。使用者可直接运行脚本并替换材料参数，快速获得前若干阶自由边圆板的模态参数，同时通过dat文件对比不同泊松比下的结果，方便检验计算一致性。对于需要开展圆板振动特性研究的初学者，该资源提供了从理论公式到数值实现的完整参考路径。 圆板类的结构，恐怕是所有机械件里最容易跟振动较劲的一种。硬盘盘片转速一上去就开始抖，圆锯片切硬料的时候尖叫，音箱振膜高频段失真，背后全是同一类问题：这块板在哪几阶频率上会共振。我这次梳理的自由边圆板振动计算，就是给你一套用MATLAB把自由边界条件下圆板固有频率和模态形状系数完整算出来的代码和思路。输入材料参数、半径、厚度，程序自动扫描特征方程，输出各阶固有频率，还能还原成三维振型图，整个过程不需要几何建模，也不需要网格。

这个题目听起来是振动力学教材里的经典题，但真正动手写代码会碰到不少坑。贝塞尔函数组合出来的特征方程高度非线性，根不是肉眼看出来的；边界条件稍不注意就列错，公式里抄错一个符号，整个行列式结果就废掉。所以我这篇文章不打算只丢一个脚本，而是把模型怎么建、方程怎么来、代码怎么组织、调试时踩过哪些坑，一层层说清楚。不管你是做课程设计、研究生课题，还是在做盘片类产品结构校核，这份东西都能直接用得上。

1. 项目解读：这份代码到底在算什么

1.1 自由边圆板的实际工程场景

自由边圆板，工程上最常见的例子就是硬盘盘片和光盘。盘片中心被主轴电机夹持，外缘完全悬空，边界条件就是“自由”——没有支撑，也没有外力。圆锯片、砂轮片、扬声器纸盆、光学镜片，甚至MEMS麦克风里的薄膜，在简化模型时都可以抽象成一块外缘自由的圆板。这些结构最怕什么？怕工作转速或者激励频率落在某阶固有频率上，一共振就是噪声、磨损、疲劳断裂。

所以算固有频率是一件事，算出模态形状系数是另一件事。固有频率告诉你“什么转速会出事”，振型告诉你“出事的时候板子是怎么弯的”。比如最低阶弹性模态通常对应两节径的蝶形变形，这时候外缘上下大幅摆动，对动平衡和噪声的影响就完全不一样。标题里提到的“circular-plate-free-edge”这个搜索词，在振动分析里基本就是指向这类问题的关键词组合。

1.2 为什么选解析解而不是直接上有限元

直接拿ANSYS或者Abaqus建个圆板模型，几秒钟也能算出频率。那为什么还要写MATLAB解析解？原因有三个。

第一是快。有限元每次改厚度、改半径都要重建模型重新划分，解析解只需要改两个输入参数，一秒钟出结果，做参数扫描的时候差别巨大。第二是物理图像清晰。解析解能直接告诉你频率跟厚度、半径、弹性模量、密度之间的幂次关系，这比有限元跑一百个工况更有工程指导价值。第三是基准验证。有限元算完心里没底的时候，用解析解一对比，能快速确认边界条件有没有加错、网格够不够密。

当然，解析解只能处理理想圆板。偏心夹持、环板、变厚度这种就得靠有限元，或者后面第6部分讲的扩展方法。两者不是替代关系，是配合关系。

2. 数理模型与特征方程推导

2.1 控制方程和振型函数

薄板横向自由振动满足经典四阶偏微分方程：

D ∇⁴w + ρh ∂²w/∂t² = 0

其中 D = Eh³/[12(1-ν²)] 是弯曲刚度，E是弹性模量，ν是泊松比，h是厚度，ρ是密度。这个方程假设板厚远小于半径，变形符合直法线假设，也就是Kirchhoff薄板理论。

对方程做时间谐波分离，w = W(r,θ)e^{iωt}，在极坐标下求解空间部分。圆板是各向同性的，振型沿圆周方向自然按 cos(nθ) 或 sin(nθ) 展开，n对应节径数。代入方程后发现，径向函数R(r)的解是四类贝塞尔函数的线性组合：第一类贝塞尔函数J_n、第二类Y_n、第一类修正贝塞尔函数I_n、第二类修正K_n。

对于实心圆板，中心点处Y_n和K_n发散，物理上不可能，所以这两项必须去掉。振型函数就简化成：

R(r) = A·J_n(βr/a) + C·I_n(βr/a)

这里引入无量纲频率参数β，a是半径。最终固有圆频率和β是平方关系：

ω = β²·√[D/(ρh·a⁴)]

这个式子是整个求解的核心。算出β，频率立刻就出来了。后续所有工作，本质都是在解β。

2.2 自由边边界条件怎么列

自由边意味着边界上既没有弯矩，也没有等效横向剪力。具体到圆板外缘 r=a 处，要满足两个条件：

第一，径向弯矩Mr等于零。第二，等效横向剪力Vr等于零。

为什么用“等效横向剪力”而不是直接让剪力Qr等于零？因为薄板理论里，边界上的扭矩和剪力是耦合的，Kirchhoff把二者合并成一个等效横向剪力，才能在边界上恰好给出两个标量条件。这也是板振动和梁振动在边界处理上最大的区别，刚接触这个领域的人特别容易在这里栽跟头。

把振型函数R(r)代入这两个边界方程，自然会得到一组关于未知系数A和C的齐次线性方程。为了让A和C有非零解，这个方程组的系数行列式必须等于零，这就是频率方程（特征方程）。

2.3 频率方程的行列式形式

先看最简单的轴对称情况 n=0。此时振型与角度θ无关，边界条件的推导可以完整手推出来，结果非常漂亮。

利用贝塞尔函数的递推关系和微分方程，径向弯矩条件化为：

A·[J₁(β)(1-ν)/β - J₀(β)] + C·[I₀(β) - I₁(β)(1-ν)/β] = 0

等效剪力条件更简洁，直接得到：

A·J₁(β) + C·I₁(β) = 0

把这两个方程写成矩阵形式，行列式等于零：

2(1-ν)/β · J₁(β)·I₁(β) - [J₀(β)·I₁(β) + I₀(β)·J₁(β)] = 0

这个方程就是n=0模态的频率方程，程序里按这个式子写不会错。对于n≥1的情况，推导过程类似但表达式更长，还要考虑扭矩项，所以代码里我统一用数值组装2×2矩阵的方式来处理，每个矩阵元素由贝塞尔函数及其递推关系现场计算，避免手抄公式抄错。对照文献的话，Leissa的NASA SP-160《Vibrations of Plates》第二章有完整的系数表，可以直接校验。

3. MATLAB代码实现：求解流程与核心函数

3.1 程序框架怎么搭

整个程序按四个文件拆，主脚本、频率行列式函数、求根函数、振型生成函数。主脚本只负责输入参数和展示结果，具体算法全部封装成函数，这样换一组参数不用动代码逻辑。

主脚本输入块很简单：

% 材料参数 E = 70e9; % 弹性模量，Pa（铝合金） rho = 2700; % 密度，kg/m^3 nu = 0.33; % 泊松比 a = 0.1; % 半径，m h = 0.002; % 厚度，m n_modes = [0 1 2 3]; % 需要计算的节径数 D = E*h^3/(12*(1-nu^2)); fprintf('弯曲刚度 D = %.4f N.m\n', D);

3.2 频率行列式函数

频率行列式是计算的核心。n=0直接用推导好的方程，n≥1用通用矩阵组装。代码结构大致如下：

function detVal = freeEdgeDet(beta, n, nu) % 自由边圆板频率行列式 if n == 0 J0 = besselj(0, beta); J1 = besselj(1, beta); I0 = besseli(0, beta); I1 = besseli(1, beta); detVal = 2*(1-nu)/beta*J1*I1 - (J0*I1 + I0*J1); else % n>=1 通式，利用贝塞尔递推关系组装 2x2 矩阵后求行列式 [m11, m12, m21, m22] = freeEdgeMatrix(beta, n, nu); detVal = m11*m22 - m12*m21; end end

这里有一个容易被忽略的细节：贝塞尔函数在MATLAB里分别是besselj、bessely、besseli、besselk，函数名属于老牌数值库，精度和可靠性都经过大量验证，放心用。修正贝塞尔函数I和K随自变量增长一个暴涨一个骤减，后面会专门讲溢出问题。

3.3 求根策略：不要直接拿fzero乱试

这是这个项目最关键的工程经验。

频率行列式det(β)是极度振荡的函数，直接叫fzero并给一个初值，十有八九会飞到一个莫名其妙的根上去。原因是行列式振荡剧烈，初值稍微偏离一点，牛顿迭代就不知道该收敛到哪一边了。

我的做法是分段扫描加括号细化。先在0到β_max区间均匀取几千个点，计算行列式，找到所有变号的区间，然后以每个变号区间作为fzero的括号约束，精确求根。代码很简单：

function betas = scanRoots(n, nu, betaMax, N) if nargin < 4, N = 5000; end x = linspace(1e-4, betaMax, N); F = arrayfun(@(b) freeEdgeDet(b, n, nu), x); dF = diff(sign(F)); idx = find(dF ~= 0); betas = zeros(1, length(idx)); for k = 1:length(idx) a = x(idx(k)); b = x(idx(k)+1); betas(k) = fzero(@(t) freeEdgeDet(t, n, nu), [a b]); % 残留判断：防止扫描点太疏混入假根 if abs(freeEdgeDet(betas(k), n, nu)) > 1e-8 betas(k) = NaN; end end betas(isnan(betas)) = []; end

扫描点数量N要跟β_max匹配。经验值是保证每个振荡周期内至少有10个采样点，否则极窄的“双根”区间会被漏掉。β_max取多少可以看贝塞尔函数的性质：J₀在β=2.4048处过零点，I₀单调暴涨，自由边圆板的高阶模态间距会随着β增大逐渐趋近一个常数，实践里取β_max=80已经能覆盖工程关心的前十几阶。

3.4 频率换算与模态排序

得到β之后，按前面的平方关系换算频率：

omega = beta.^2 .* sqrt(D/(rho*h*a^4)); freq = omega / (2*pi);

这里要特别提醒，β必须是无量纲频率参数本身，不是它的平方。有些文献里把“频率参数”直接定义成λ² = ρh a⁴ω²/D，跟这里的β²是一回事，换算的时候要分清口径。全自由圆板还存在β=0的刚体模态，对应平动和刚体转动，扫描时通常会把它们列为β≈0的根，筛选时直接忽略前两个极小值即可。

4. 模态形状系数计算与振型可视化

4.1 模态形状系数的意义和求法

标题里说的“模态形状系数”，指的就是振型函数里A和C的比值。它决定了径向截面到底是J项主导还是I项主导，最终表现出什么样的变形轮廓。

求法有两种。一种是用边界条件中任何一个方程代数算比值，比如n=0时可以直接由 A·J₁ + C·I₁ = 0 得到 A/C = -I₁/J₁。但这个做法有个隐患：当J₁恰好接近零时，比值会发散。更稳妥的办法是把两个边界方程组装成系数矩阵，然后用MATLAB的null函数求零空间向量：

function [Acoef, Ccoef] = modeCoeff(beta, n, nu) M = freeEdgeMatrix(beta, n, nu); % 通用矩阵组装 v = null(M); Acoef = v(1); Ccoef = v(2); % 归一化，方便比较振型 Acoef = Acoef / sqrt(Acoef^2 + Ccoef^2); Ccoef = Ccoef / sqrt(Acoef^2 + Ccoef^2); end

null函数本质上是做奇异值分解，数值稳定性比手写代数表达式好得多。用这种方法，得到的A和C就是模态形状系数，振型的径向形状完全由这两个系数决定。

4.2 三维振型图与节线判读

得到系数之后，把径向函数和圆周方向cos(nθ)乘起来，就能画出整个面上的位移分布：

r = linspace(0, a, 100); theta = linspace(0, 2*pi, 100); [RR, TH] = meshgrid(r, theta); Rr = Acoef*besselj(n, beta*RR/a) + Ccoef*besseli(n, beta*RR/a); W = Rr .* cos(n*TH); [X, Y] = pol2cart(TH, RR); surf(X, Y, W, 'EdgeColor', 'none'); colormap(jet); colorbar;

画完图重点看两个特征：节径和节圆。节径是穿过圆心的零位移线，n=0没有节径，n=1有一条直线节径，n=2有两条互相垂直的节径形成“蝶形”变形。节圆是同心圆形状的零位移环。自由边圆板的第一阶弹性模态通常是n=2的那个蝶形模态，频率最低，工程上最危险的就是它。

4.3 结果验证：行列式残差与正交性检查

代码写完最怕的是“看起来对，实际错”。我强烈建议做两个验证。

第一个是残差检查。把求出来的每个β代回频率行列式，确认绝对值在1e-10量级。如果残差偏大，基本可以断定是漏根或者扫描精度不够，增大N重新扫。

第二个是正交性检查。理论要求不同模态之间满足质量正交关系，也就是两个不同模态位移的乘积在整个板面上积分应该为零：

∫₀^a ∫₀^{2π} ρh·Wᵢ·Wⱼ·r drdθ = 0，i≠j

用数值积分跑一圈，正交性能到1e-6以下，就说明振型函数、系数、边界条件都没有原则性错误。这一步是我个人认为所有振动计算里最值得做、也最容易被忽略的验证。

5. 常见问题与调试心得

5.1 贝塞尔函数溢出和NaN问题

修正贝塞尔函数I_n(β)在自变量大的时候会爆炸性增长，β超过150左右，besseli直接给你返回Inf，K_n直接返回0，行列式瞬间变成NaN。处理办法有两个方向。

第一个方向是用MATLAB提供的自适应缩放参数：besseli(n, β, 1)和besselk(n, β, 1)，这两个函数会额外乘上exp(-β)或exp(β)，把数值范围拉回可计算区间。第二个方向是限制搜索范围，β_max控制在100以内，工程上通常足够了。真需要算超高阶模态，建议把频率方程改写成J和I的比值形式，比如(J/I)，把暴涨因子消掉再算。

5.2 漏根、假根与刚体模态

漏根几乎都是扫描间隔太宽导致振荡区间被跳过去了。解决办法是按采样密度要求回推N，别偷懒用几百个点扫到100。

假根则出在符号变化的判断上。频率行列式在某些β处可能非常接近零但没有真正过零，浮点数误差会让sign函数误判。滤除方法就是在fzero收敛后增加一次行列式残差判断，绝对值超过1e-8一律当作假根剔除。

刚体模态是自由边圆板特有的问题。整个板自由悬浮时，它可以整体上下平移（n=0刚体模态）和刚体转动（n=1刚体模态），对应β=0。程序扫描时会在零点附近扫出一个根，这不是弹性模态，排序时可以直接丢弃。

5.3 参数单位与薄板适用范围

单位不一致是新手最常见的错误。弹性模量用GPa，密度用kg/m³，半径用mm，结果就是一通算下来频率差了好几个量级。建议全系统一使用国际单位：E用Pa，长度用m，密度用kg/m³，频率输出自然就是Hz。如果是自己封装函数，可以约定输入E的单位为GPa，但一定要在注释里写清楚，否则两个月后回来看代码绝对会懵。

薄板理论也有适用范围。经典四阶方程假设厚度远小于半径，工程经验是h/a < 0.1。超过这个比例，剪切变形和转动惯量开始显著影响高阶模态，需要换Mindlin板理论或者直接上三维有限元。遇到厚圆板，解析解只适合做定性参考。

我把这些常见问题整理成一个速查表：

6. 这个程序怎么扩展

算完实心自由边圆板之后，这套框架的价值在于能往多个方向扩展。如果遇到环形板（中心带孔），振型函数里Y_n和K_n两项要加回来，未知系数变成四个，边界条件在内外两个边界上一共四条方程，频率行列式从2×2变成4×4，求解思路完全一致。如果遇到中心夹持外缘自由的盘片，边界条件改成内孔固定、外缘自由，同样只需换一组边界方程。

加筋圆板、变厚度圆板这类情况，解析解就很难直接套了，但可以把这个程序作为基准解，用来校验有限元模型。哪怕只是算完自由边圆板后，顺手用ANSYS建个同样尺寸的模型对比前五阶频率，误差在1%以内，那你的有限元边界条件基本就是可信的。这就是解析解在实际工程里最大的作用——它不是用来替代仿真，而是用来给仿真兜底。

最后多提一句，这个程序在改参数的时候有个小技巧：把半径a、厚度h、弹性模量E这些输入参数写成一个结构体变量struct，比如plate.E = 70e9这样，频率计算和振型绘制都从同一个结构体里取数，避免多个脚本之间参数不一致。我早期做参数扫描时，就是因为主脚本和绘图脚本里各写了一遍参数，厚度改了一处忘了另一处，白白浪费了大半天排查时间。这算是这个项目里最不起眼，但最实在的一条经验了。

本文还有配套的精品资源，点击获取