企业官网建设流程全解析

简介：这份MATLAB代码源于Sandberg著作中的声振耦合算例，面向刚接触声振耦合有限元方法的初学者，旨在演示如何用梁单元与矩形声学单元搭建结构与声场耦合模型。代码中已加入较详细的中文注释，帮助读者理解有限元建模、声学方程离散以及结构振动与声场之间的双向作用，特别适合用于课程学习或自学入门。压缩包内共有1个文件，为m脚本（约3KB），主体是ASI_full_code.m，可直接运行并对照注释逐步拆解，无需额外复杂工程文件。目前已有223人学习下载。通过学习这份代码，初学者可以掌握在MATLAB中实现声振耦合求解的基本流程，弄清梁单元、矩形声学单元与声振耦合的关键处理思路，为后续深入结构声学分析打下基础。

1. 为什么从Sandberg的梁-声腔例子入手学声振耦合

低频NVH问题里，结构振动和封闭声腔之间是强耦合的。最直观的现象是：一块薄板背后有个空腔，板振动得越欢，腔内的声压反过来又把板推回去；这个“推回去”的力在某些频率上会明显改变结构模态频率和阻尼。如果只把结构有限元和声学有限元分开算，会漏掉这种附加质量效应和反作用力。

Sandberg那本关于声振耦合有限元的书里，恰好安排了一个用梁模拟结构、用矩形声学单元模拟二维声腔的例子。ASI_full_code.rar解压出的ASI_full_code.m就是这个例子的MATLAB实现，代码里加了中文注释。适合两类人：一类是学过有限元理论但没见过“耦合项在矩阵里怎么摆”的初学者；另一类是用商业NVH软件做分析、想验证自己耦合矩阵组装逻辑的工程师。下面按结构侧、声学侧、求解和验证四个层面拆开讲。

2. 结构侧建模：梁单元自由度约定与ASI_full_code的组装顺序

2.1 梁节点的自由度表与界面法向约定

在二维梁-声腔模型里，梁被当成一维结构放在声腔边界上，但横向弯曲会沿法线方向推动声学介质。每个节点需要三个自由度：轴向位移u、横向位移w、绕面外轴的转动θ。如果用Euler-Bernoulli梁理论，横向位移的导数给出转动，所以自由度不需要再加高阶项。

需要注意：ASI_full_code里梁位于声场上方，结构侧的法向就是全局y方向，所以只有w自由度与声压直接耦合。如果模型里梁倾斜，就必须把u和w投影到声界面法线方向；不要把θ直接和声压关联起来，这是耦合建模时最容易出错的点。

2.2 Euler-Bernoulli梁单元的刚度与质量矩阵

单元矩阵并不需要每次都从形函数重推。对均匀截面Euler-Bernoulli梁，局部坐标系下单元刚度矩阵是标准形式，质量矩阵用一致质量近似。组装时先由节点编号生成自由度向量，再按自由度数填入全局矩阵：

% 向量式梁单元自由度编号（单元连接node1-node2） dofs = [3*node1-2, 3*node1-1, 3*node1, ... 3*node2-2, 3*node2-1, 3*node2]; K(dofs, dofs) = K(dofs, dofs) + Ke; M(dofs, dofs) = M(dofs, dofs) + Me;

Ke、Me是单元刚度、质量矩阵；dofs按“轴-横-转”的顺序排列，这样和大多数教材的梁单元自由度定义一致。在稀疏矩阵上逐单元赋值效率很低，常见做法是先把每个单元的dofs和矩阵值存到三个临时数组里，最后用sparse一次性组装；对几十个单元的教学模型差异不明显，改成几百个单元后能差出一个数量级。

实际的单元矩阵如果用一致质量，长度和弯曲自由度耦合会产生一个2×2分块结构。要验证组装是否正确，可以做一个刚体平移测试：把K乘一个全部节点u=1、w=0的向量，结果应该为零向量；把M乘这个向量，得到的总质量应等于所有单元质量之和。

u_test = zeros(ns_dof,1); u_test(1:3:end) = 1; % 所有节点轴向单位位移 force = K * u_test; total_mass = u_test' * M * u_test;

如果force的模超过1e-10，说明梁单元矩阵里有坐标变换或自由度顺序的问题。这个自检在改写单元矩阵时非常有用。

2.3 边界条件与自由度压缩

固支端要约束w和θ，简支端约束w，轴向约束看模型是否需要。耦合分析中，与声学域接触的边界自由度不能压缩，否则耦合矩阵的列数对不上，组装时维度就崩了。

% 节点1固支，末端节点约束横向位移 fixed_dofs = [1 2 3]; % 节点1的u,w,theta fixed_dofs = [fixed_dofs, 3*N-1]; % 节点N的w free_dofs = setdiff(1:ns_dof, fixed_dofs);

重排后，结构矩阵被分成free和fixed两个分块，求解时只对free_dofs做分解；但写耦合矩阵时仍使用完整结构编号。setdiff返回的是排过序的索引，如果你后续还需要按原节点顺序组装声学界面，一定先把编号向量单独存好，不能在原数组上覆盖。这个细节在ASI_full_code这类教学代码里是一个很常见调试点。

另外，ASI_full_code.m从rar里解出来后，如果直接在旧版MATLAB里打开，中文注释可能显示成乱码，但这不影响运行结果。最好的做法是把.m文件另存为UTF-8编码，或者用实际使用的MATLAB版本重新读取。你也可以顺手在命令行跑一个单元级测试：构造两个节点的梁，把K和M与手算值比较，这样之后改任何参数都不会心里没底。

3. 声学侧建模：矩形声学单元、流体矩阵与耦合界面的处理

3.1 矩形声学单元的压力自由度与形函数

二维矩形声学单元每个节点只有一个压力自由度p，四个节点构成双线性单元。均匀介质中声场用Helmholtz方程描述，加权余量后得到流体刚度矩阵Kf和质量矩阵Mf。声学单元与结构单元最大的不同在于未知量是标量，所以单元矩阵是4×4，规模小很多。

形函数在局部坐标(ξ,η)下取双线性插值：N1=(1-ξ)(1-η)/4，N2=(1+ξ)(1-η)/4，N3=(1+ξ)(1+η)/4，N4=(1-ξ)(1+η)/4。做2×2高斯积分即可满足精度。ASI_full_code里采用一致质量矩阵而不是集中质量矩阵，因为集中质量近似会让高频声模态偏硬，且设置完全刚性壁边界时容易出现伪模态。初学者看到流体质量矩阵里有1/ρ0、刚度矩阵里有ρ0c0²时，往往怀疑是单位错了，其实这正是声学有限元的定义方式。

3.2 组装流体刚度与质量矩阵的MATLAB流程

下面这个片段只负责组装若干4节点矩形单元，没加边界吸收项，完整代码里会在右边界再加一层阻尼条件。为了突出声学矩阵的组装逻辑，把它独立出来：

function [Kf, Mf] = assem_fluid(nodes4, rho0, c0) for e = 1:size(nodes4,1) % 对每个四节点矩形单元 xy = squeeze(nodes4(e,:,:)); % 4x2坐标矩阵 [gp, gw] = gauss2d(2); % 2x2高斯点与权重 ke = zeros(4,4); me = zeros(4,4); for i = 1:4 [N, dN] = shape_quad(gp(i,:)); % 形函数与局部导数 J = dN * xy; % 2x2雅可比矩阵 dNxy = J \ dN; % 全局坐标导数 ke = ke + (rho0*c0^2) * (dNxy'*dNxy) * det(J) * gw(i); me = me + (1/rho0) * (N'*N) * det(J) * gw(i); end dof = [e*4-3:e*4]; % 假定的连续节点编号 Kf(dof,dof) = Kf(dof,dof) + ke; Mf(dof,dof) = Mf(dof,dof) + me; end end

ke对应声学势能，me对应声学惯性项。乘rho0c0²再乘det(J)后，尺寸正好是压力×体积，频率单位才能对得上。如果漏掉rho0c0²，求出的特征频率无量纲且数值离谱。矩形单元的det(J)通常是边长乘积的一半再乘以4，出现负值说明节点顺序反了，必须逆时针排列。ASI_full_code的网格生成函数按逆时针生成，但当你手动补网格时，这是最常踩的雷。

3.3 耦合矩阵：从结构法向速度到声学边界激励

耦合矩阵C的行数是结构自由度数，列数是声学自由度数，只在接触面上非零。物理上，结构侧法向振动速度作为声学域边界速度源；反过来，声压合力又作为结构外载荷。离散方程里会出现两项，一项是-Cp，一项是ρ0C^T u，二者不能互消，因为一个作用在结构方程，一个作用在声学方程。

% 耦合界面自由度匹配：结构节点snode与声学节点anode一一对应 Coupling(3*snode-1, anode) = -boundary_width; % 梯形分配

这段示意假设界面网格匹配，每个结构节点对应一个声学节点。如果两边网格不等，不能直接填这一个元素，先用投影矩阵做插值。常见做法是保证界面网格节点重合，这样耦合矩阵不需要额外插值，也最方便调试。一个快速自检是看Coupling的行和或列和：对封闭声腔，行和应近似等于界面长度；符号则取决于界面法向定义，如果所有频率都明显偏高或偏低，多半是这里差了一个负号。

3.4 刚性壁与阻抗边界条件的扩展

ASI_full_code原始模型里，声腔除耦合界面外都设为刚性壁。刚性壁条件的处理方式是把边界自由度和普通内部自由度一样保留，不做额外约束，因为声学有限元中刚性壁是自然边界条件。如果想模拟开口端，需要把压力释放边界自由度从总矩阵中删掉。问题在于，矩形声学单元在四角重叠处容易出现法向不唯一，删自由度时容易删错。

% 在边界单元上叠加吸声项，扩展成阻抗边界 Kd = zeros(nf, nf); for be = 1:n_boundary Kd(edof, edof) = Kd(edof, edof) + ... (1i * omega / (rho0 * c0 * Z_n)) * Mb_edge; end

这里Mb_edge是边界边的质量矩阵，Z_n是法向声阻抗率。加入该项后总刚度矩阵变成复数，特征值求解从eig(A,B)变成复广义特征值问题，计算量明显增加。教学阶段建议先不加，跑通刚性壁再接阻抗边界，否则排错时很难判断是耦合矩阵的问题还是边界矩阵的问题。

4. 把ASI_full_code.m跑通：求解流程、参数设置与频率响应

4.1 主程序的数据流与矩阵分块

ASI_full_code.m的流程很直接：先生成结构网格，再生成声学网格，分别组装K_s、M_s、K_f、M_f和耦合矩阵C，然后求解特征值或扫频响应。为了不让变量名混在一起，建议把五个矩阵单独命名，最后组成分块矩阵。标准声振耦合方程写成广义特征值问题如下：

MATLAB代码可以直接按照分块展开：

A = [Ks, -Coupling; sparse(nf, ns), Kf]; B = [Ms, sparse(ns, nf); rho0 * Coupling', Mf]; [V, lambda] = eig(A, B); freq = sqrt(diag(lambda)) / (2*pi);

eig(A,B)返回广义特征值，lambda虚部接近零时说明总矩阵组装正确。注意C的转置出现位置和符号：声学方程里耦合项是ρ0ω² C^T u，所以质量分块B(2,1)=ρ0C^T。有些教材把C定义为声压对结构的作用力而不是结构边界速度，C的符号相反。ASI_full_code原代码的C符合“结构方程-声压”这一约定，抄到别处时先确认定义。

运行前还可以加两条维度断言，快速定位矩阵写反的问题：

assert(size(Coupling,1) == ns_dof, '耦合矩阵行数应等于结构自由度数'); assert(size(Coupling,2) == nf_dof, '耦合矩阵列数应等于声学自由度数');

这两行会在矩阵尺寸不一致时直接报错，而不是等到eig阶段抛出难以理解的维度异常。对从经典教程抄矩阵组装代码的人来说，这是成本最低的防御性写法。

4.2 参数设置：材料、网格密度与频率范围

先给一组能跑出稳定模态的初始参数：梁弹性模量2.1e11 Pa、密度7800 kg/m³、矩形截面宽0.02 m、高0.01 m；声腔介质密度1.21 kg/m³、声速343 m/s。梁取8~16个单元，声腔取8×8个矩形单元。结构模态的前3阶应与手算悬臂梁理论解接近，声学模态的前几个应与矩形声腔解析解接近。

如果扫频上限取得太高，矩形声学单元会因色散误差产生明显的频率偏移，响应曲线上的尖峰位置随着网格加密不断左移或右移。判断方法是对同一模型跑32个和64个声学单元，看目标峰频率变化；变化超过2%，就继续加密。这里不要迷信“网格越多越准”，声振耦合问题里结构网格和声学网格的界面上节点必须对齐，只加密一侧反而会让耦合矩阵插值误差变大。

4.3 扫频响应求解与后处理

只想要某个频段的传递函数时，不必算全部特征值。用直接法多次求解一个大型线性方程组更省时间：

om = 2*pi*freq; for k = 1:numel(freq) Atot = [Ks - om(k)^2*Ms, -Coupling; -rho0*om(k)^2*Coupling', Kf - om(k)^2*Mf]; q(:,k) = Atot \ [Fs(:,k); zeros(nf,1)]; end % 输出最后一个声学节点声压 p_end = q(ns+end_anode, :); semilogy(freq, abs(p_end));

矩阵Atot顺序与特征值分块一致，但声学行中多了一个-ρ0ω²C^T项，等于从质量矩阵里移项。求解器直接左除，对几十阶自由度的教学模型足够快。如果Atot接近奇异，先检查结构侧是否固支：梁如果完全没有约束，A_total在低频会近似奇异，响应曲线会出随机尖峰。也可以在结构方程里加一个小阻尼项，比如C_d=0.01Ms，用实部把峰值抹平一点，便于观察趋势。

5. 用模态对比和矩阵检查验证声振耦合结果的可信度

拿到一份能跑的ASI_full_code后，第一步不是直接看声压云图，而是做一次“解耦测试”。把Coupling设为零矩阵，分别求解结构侧和声学侧的模态频率，记录前几阶；恢复耦合后重新求解，再对比频率变化。物理上声振耦合会让相近的模态频率互相靠近，低频段通常只移动几个百分点。如果某个模态频率在耦合后忽然掉了一半，或者出现复数特征值，大概率是耦合矩阵符号反了或界面法向不一致。

另一个高效验证是检查总矩阵的结构。理想无阻尼声振耦合的特征值lambda应该是正实数；如果虚部明显不为零，说明A和B组装时引入了非对称伪影。此时使用spy(A)、spy(B)看非零块分布，耦合矩阵的非零部分应只出现在界面自由度对应的行列；如果结构刚度矩阵里出现声学自由度的非零项，说明dofs索引在某个循环里被覆盖了。这个检查能定位大多数“矩阵装错位置”的bug。

最后是一个实战技巧：把单元节点顺序的检查写进程序。矩形声学单元组装前，对每个单元计算det(J)并判断正负，一旦发现负值就直接报错并输出单元编号。这个检查对网格生成和手动修改非常有用，也是教学代码最容易忽略的点。完成这些检查后，再绘制声压云图或结构振动响应，你就能确认结果可信，并开始把梁替换成板、把二维声腔替换成三维声场。

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