前端时间在做一个近岸消浪设施的可行性评估，甲方给了一个很具体的问题：一排垂直插在水里的薄板，浪打过来之后还能剩多少？这个"剩多少"在工程上就是透射系数。当时我第一反应是直接套单板公式，但对方强调"一排板、多块板"，这就没法偷懒了——多块垂直薄板之间的多重反射会把问题变得相当有意思。后来我用Matlab把基于特征函数展开法的求解流程完整跑通了，也就是题目里说到的这套针对"水波在多个垂直薄板下的透射系数"的数值计算工程。
这篇文章就把这个项目的核心内容拆开讲：从透射系数的物理定义，到多板问题的数学建模，再到Matlab里怎么把连续方程离散成能解的线性方程组。文末会结合我调试过程中踩过的坑，给出源码的使用建议和扩展方向。无论你是做港口海岸工程、波浪能装置设计，还是纯粹想学特征函数展开法在Matlab里的落地方式，这期内容都能省下你不少摸索时间。
1. 透射系数这概念，为什么工程上绕不开它
1.1 一块竖在水里的板，挡不住全部波浪
先说个反直觉的事实：你把一块薄板竖直插在波浪传播路径上，哪怕板从水面一直延伸到很深，波浪照样能绕过去——因为波浪的能量不仅在水面，它沿着水深分布。对线性波而言，波能主要集中在自由表面附近，但衰减模态在近水面区域也有不可忽略的贡献。单纯一块部分浸没的板，一般只能消掉一部分能量，剩下的就从板底绕过去了。
这就是为什么实际问题里会出现"多个垂直薄板"这种布置。一块不够，就摆一排，让波浪在板与板之间反复反射、干涉，把能量尽量"困住"。透射系数Kt定义很简单：
Kt = Ht / Hi
其中Hi是入射波高，Ht是透射波高。Kt=1代表波浪完全无阻碍通过，Kt=0代表完全挡住。实际工程中，Kt能压到0.3以下就已经算很有效的消浪措施了。
1.2 反射系数Kr：另一个必须同时盯住的量
波浪碰到薄板，除了透射还有反射。反射系数Kr = Hr / Hi，Hr是反射波高。在无能量损耗的理想流体模型里，入射能量最终只有两条出路：透过去或者弹回来。所以有：
Kr² + Kt² ≈ 1
这个关系式看着简单，却是后面验证数值程序是否写对的"照妖镜"。如果你的程序算出来的Kr² + Kt²明显偏离1，不用怀疑，肯定是代码哪里出了问题——要么边界条件配错了，要么特征函数截断项数不够，要么矩阵求解精度崩了。
我想提醒一点：别只盯着Kt。反射波往回传，同样会对上游结构物产生作用力。在港口布置里，如果你前面还有其它设施，反射波带来的附加波高往往是设计者容易漏掉的地方。所以算透射系数时，顺手把反射系数一起算出来，是基本素养。
1.3 为什么关心"多个"板而不是单板
单板问题有解析解，很多教材里都能查到结果。但多板的困难在于：板与板之间的距离不是随便给的，波在两块板之间来回反射，形成一种类似法布里-珀罗干涉的效应。某些波长下，板间反射波同相叠加，透射极小；某些波长下反相抵消，透射反而变大。这就是后面要细说的"布拉格共振"现象。
多板系统的透射系数与板数N、板间距S、板浸没深度d、水深h、入射波频率ω全都相关。工程上设计消浪结构，本质上是调参优化Kt。你没法靠直觉猜，必须得有一套能快速扫描参数的数值工具——这正是这套Matlab源码存在的意义。
2. 理论建模：多板透射问题的数学语言
2.1 线性势流理论的框架
这个问题走的是经典线性波理论路线。假设流体无黏、不可压缩，流动无旋，那么存在速度势Φ，满足拉普拉斯方程：
∂²Φ/∂x² + ∂²Φ/∂z² = 0
对简谐波，把时间因子exp(-iωt)分离掉，得到空间速度势φ(x,z)。边界条件包括：
- 自由水面条件（z=0）：∂φ/∂z = ω²/g · φ
- 海底条件（z=-h）：∂φ/∂z = 0
- 板面条件（x=x_j, -d<z<0）：∂φ/∂x = 0
- 板下开口区（-h<z<-d）：压力和法向速度连续
- 无穷远辐射条件：远离结构时只有向外传播的波
这里板从水面z=0向下延伸到深度d，水深h大于d，板底与海底之间留出通道。这对应的是最常见的"悬挂式部分浸没板"工程布置。如果你想模拟从海底向上生长的板，只需要把板面条件的z区间改成-h到-h+d，完全可以用同一套代码框架处理。
2.2 特征函数展开法的核心思想
在均匀水深条件下，沿水深方向的模态是已知的。把速度势在某一水平区间内展开成级数：
φ(x,z) = [A₀·exp(ik₀x) + B₀·exp(-ik₀x)]·Z₀(z) + Σ[Aₙ·exp(kₙx) + Bₙ·exp(-kₙx)]·Zₙ(z)
其中k₀是传播模态的波数（实数），kₙ(n=1,2,...)是衰减模态的衰减系数（正实数，对应特征值纯虚数）。Z₀(z)和Zₙ(z)分别是垂向特征函数：
Z₀(z) = cosh(k₀(z+h)) / cosh(k₀h) Zₙ(z) = cos(kₙ(z+h)) / cos(kₙh)
传播模态的波数k₀由色散关系决定：
ω² = g·k₀·tanh(k₀h)
衰减模态的kₙ满足：
ω² = -g·kₙ·tan(kₙh)
这两条方程是所有实现里最先要写对的地方。色散关系求根出错，后面全盘皆输。
2.3 多板问题的求解策略：分区间匹配
为什么要强调"多个"板？因为每块板都是一个散射体。如果只有一块板，你只需要把流场分成板前、板下、板后三个区间，写匹配条件就行。但多块板沿x方向排布时，相邻板之间的区域既是前一块板的"板后区"，又是后一块板的"板前区"。波浪在中间区域反复反射，你必须把所有区间的待定系数一起解出来。
具体做法是把计算域沿x方向切成N+1个水平区间（N块板）：
- 区间0：x<x₁，入射区
- 区间j：x_j < x < x_{j+1}，各板间区域
- 区间N：x>x_N，透射区
在每个区间内，速度势都写成2.2节形式的特征函数展开。只不过入射区里要求A₀=1（归一化入射波幅），透射区里不能有内向反射波（B系数全部为零）。这样每个区间的未知系数就是各自的传播模态幅值和衰减模态幅值。然后在每个板的位置x=x_j处，把上段板面条件和下段开口区的连续条件逐一写上，最终组装成一个复系数线性方程组：
M·X = b
解出X之后，透射系数Kt直接就是最右区间透射波的传播模态幅值系数，反射系数Kr就是最左区间反射波的幅值系数。这是整套Matlab程序的主线逻辑。
3. Matlab实现:核心流程与关键代码逻辑
3.1 程序主流程架构
我拿到或写一套这类源码时，第一件事是看它的主函数结构。正常来说，完整的Matlab工程应该包含以下功能模块：
- 参数定义模块：水深h、板数N、板间距S、板浸没深度d、入射周期T或频率ω
- 色散方程求解函数：给定ω和h，返回k₀以及前M个衰减模态的kₙ
- 矩阵组装函数：按板位置逐块写匹配条件，填充系数矩阵和右端项
- 线性方程组求解：直接调用Matlab的矩阵左除（x = M \ b即可，内部会选好求解器）
- 后处理函数：计算Kt和Kr，绘制透射系数随波长或频率变化的曲线
主程序典型运行流程可以概括成下面这段伪代码结构：
% 参数设置 h = 10; % 水深(m) N = 4; % 薄板数量 S = 4; % 板间距(m) d = 5; % 板浸没深度(m) T = 6; % 波动周期(s) omega = 2*pi/T; M = 30; % 截断模态数 % 求解色散方程获得波数和特征函数 [k0, kn] = solveDispersion(omega, h, M); % 组装线性方程组 [Amat, bvec] = assembleSystem(N, S, d, h, k0, kn, omega); % 求解 X = Amat \ bvec; % 提取透射系数和反射系数 Kt = abs(X(end)); Kr = abs(X(1));实际源码里，assembleSystem是最大也是最容易写错的部分，后面我会专门讲它的内部逻辑。
3.2 色散方程求根的数值技巧
先停留一下，因为色散关系求根是第一个拦路虎。传播模态的k₀用fzero就能稳定找到，关键在于衰减模态。
衰减模态的kₙ满足ω² = -g·kₙ·tan(kₙh)。这个方程有无穷多个正实根，分别落在区间((n-1/2)π/h, (n+1/2)π/h)附近。求根的可靠做法是：在(0, (M-1/2)π/h)范围内以很小的步长（比如h/1000）采样，检测函数符号变化，然后用fzero精确定位每个根。千万别直接用fzero乱猜初值，很容易漏根或者跳到同一个根上。
我见过很多人在这一步偷懒，结果特征函数展开不完整，匹配条件处处对不上，得到的Kt曲线出现莫名其妙的尖刺。这里值得多写两行代码老老实实扫描。
function kn = solveEvanescentModes(omega, h, M) g = 9.81; % omega^2 = -g*k*tan(k*h) 的正实根 kn = zeros(M, 1); % 在整个有理区间内扫描 kmax = (M - 0.5) * pi / h; xs = linspace(1e-6, kmax, 20000); rhs = @(k) -g * k .* tan(k * h) - omega^2; vals = rhs(xs); cnt = 0; for i = 1:length(xs)-1 if vals(i) * vals(i+1) < 0 cnt = cnt + 1; if cnt > M, break; end kn(cnt) = fzero(@(k) -g*k*tan(k*h) - omega^2, ... [xs(i), xs(i+1)]); end end if cnt < M error('截断模态数M超过可解析根范围'); end end这里有个细节：衰减模态的特征函数是cos(kₙ(z+h))，它在z=-h处天然满足海底不可穿透条件，在自由水面处通过色散关系自动满足自由水面条件。特征函数写对了，边界条件就省掉一半麻烦。
3.3 匹配条件的矩阵化：把物理条件翻译成线性代数
每块板处有两类条件要写：
第一类是板上段的刚性板面条件，0>z>-d，要求∂φ/∂x=0。但注意，板两侧的速度势并不连续，因为薄板本身可以承受压差。所以对左区间和右区间，这个条件要分别施加。在配置法实现中，通常的做法是在板面的高度范围内取若干配点，要求每个配点上两个相邻区域的水平速度分别等于0。这会带来两组方程，分别对应板左侧和板右侧的流场。
第二类是板下开口区的连续条件，-h<z<-d，要求速度势连续、水平速度连续。同样取配点，写连续性方程。
把板面条件和开口区条件各自离散后，每一块板贡献的方程个数大于未知系数个数是不行的，所以通常配合"加权余量"或者"最小二乘匹配"来保证矩阵方正。最简单稳妥的办法是伽辽金法：对每个条件乘以对应特征函数并在区间内积分，把微分/连续条件投影到特征函数空间里。虽然写起来略繁琐，但数值稳定性极好。我在后处理里通常把积分匹配后的残差打印出来检查，残差小于1e-6就说明匹配质量很高。
如果不想自己手推积分，也可以用纯配点法：在板面和开口区均匀撒点，每个点写一条方程，方程数量超过未知数时用最小二乘解。但要注意，配点法对配点位置敏感，板边缘附近（z=-d）的速度奇异性会导致局部振荡。我自己的经验是伽辽金加权积分法稳健得多，宁可多写几行积分代码。
3.4 大矩阵求解的注意事项
所有区域的待定系数装进一个向量后，方程组的规模大约是：未知数数量 = (N+1) × (2M+2)，M取30、N取4时，大概是340个左右。这对Matlab来说是小菜一碟，直接用左除就行。但有几个隐藏风险：
- 系数矩阵行与行之间如果量纲差太大（比如特征函数值从1e-3到1e3量级），求解精度会下降。解决办法是给矩阵做行均衡，或者把特征函数归一化到最大值为1。
- 板间距S非常小的时候，相邻板之间的衰减模态幅值会很大，矩阵条件数剧增。这时候要么增加截断模态数M，要么用更高精度的quad积分处理重叠积分。数值上出现警告提示矩阵接近奇异时，先检查是不是S太小或M不足。
- 高频情况下（波长短、kh大），需要更多的衰减模态才能收敛。建议在扫描频率时动态调整M，比如按kh自动给一个M = round(5 + 10 * kh / pi)的经验公式。
4. 参数影响规律：板数、间距、浸没深度如何改变Kt
4.1 板数N：挡浪效率不是线性叠加的
跑参数扫描时有一个很有意思的结论：透射系数随板数增加而下降，但下降幅度是递减的。单板Kt可能在0.6左右，加到3块板能降到0.35，再加到5块可能只降0.05。原因是多板系统里，能量主要被前几块板反射掉了，后面的板面对的波高已经很小，贡献自然有限。
所以工程上不是板越多越好。板数再多，成本上涨明显，Kt收益却越来越小。一般情况下，3到5块板是性价比较高的区间，具体数值要看目标频段。
4.2 板间距S：布拉格共振带来的透射低谷
板间距的影响比板数更微妙。当板间距S与入射波半波长λ/2满足特定关系时，各板反射波同相叠加，反射率出现峰值，相应地透射系数出现低谷。这就是周期性结构中的布拉格共振。
用这套程序扫描S时，你会在Kt随S变化的曲线上看到明显的周期性低谷。一个实用结论是：如果你希望某个设计波龄的透射系数最小，就应该把板间距取在那个波龄对应半波长附近。反过来，如果你希望透射系数比较稳定、对频率不敏感，就避开布拉格共振区间。
我调试时遇到过一个很有意思的现象：当板间距取到接近一个波长时，Kt反而可能比单板还大。这是因为板间反射在某些相位下互相抵消，相当于把"墙"拆成了"透明"结构。所以无脑加密板距是不可取的，必须靠数值扫描设计。
4.3 浸没深度比d/h：越深越好，但要尊重边际递减
板浸没深度d决定波浪可以从板底绕过的通道大小。d越接近水深h，板底与海底的缝隙越窄，透射越弱。但d/h接近1时，板底缝隙里的流速急剧增大，局部能量损失和结构受力也会显著上升。
从Kt的变化曲线看，d/h从0.2增加到0.5时，Kt下降很明显；但从0.7增加到0.9，Kt的降低幅度就小很多了。这个规律与边界元法里"缝隙流"的经典结论一致。做初步设计时，我习惯把d/h定在0.5到0.7之间，再配合板间距调优，基本上能得到不错的消浪效果。
下面是三种参数影响趋势的对比总结，方便快速查阅：
| 参数 | 增大时Kt的变化 | 效应性质 | 工程建议 |
|---|---|---|---|
| 板数N | 下降，边际递减 | 反射增强 | 常用3~5块 |
| 板间距S | 振荡，存在布拉格低谷 | 干涉效应 | 按目标波长调优S |
| 浸没深度d/h | 下降，边际递减 | 绕流通道变窄 | 建议0.5~0.7 |
5. 数值验证与调试：能量守恒是诚实裁判
5.1 能量守恒检查怎么做
程序写完后别急着画Kt曲线，先算能量守恒。前文提到理想流体中Kr² + Kt²应等于1。在Matlab里，你提取出Kt和Kr后，直接算：
residual = abs(Kt^2 + Kr^2 - 1);这个残差如果大于1e-4，就必须排查原因。我自己的调试经历中，能量残差过大的元凶通常有三个：色散方程漏根导致特征函数展开不完整；截断模态数M不够导致匹配条件不精确；板面条件里的水平速度方向写反。其中漏根最隐蔽，因为Kt曲线看起来形状正常，但能量守恒就是差那么一点。
5.2 截断模态数的收敛性测试
衰减模态的截断数M怎么选？一个稳妥做法是对同一工况连续算三组：M、2M、4M。如果Kt的变化小于0.1%，就认为收敛了。我在标准水深h=10m、周期T=5s的工况下测试，M=20时Kt已经比较稳，M=40时结果几乎不再变化。但如果板间距特别小（S<0.2h），近场衰减模态的影响显著增强，这时候M需要跟S联动，否则会看到Kt随M震荡。
另外要提醒的是，特征函数展开法的收敛性在高频端会变差。kh>3之后，衰减模态的衰减长度变短，板边缘奇异性对结果的影响范围反而变大。这时候建议把M提升到50甚至60，并且把板边缘附近的配点加密，或者在伽辽金积分里做端点加权处理。
5.3 极限情形验证：N=1与已知结果对照
收敛性测试通过后，还有一个简单有效的验证方式：把板数设成1，与文献里的单板结果对比。单板的部分浸没问题在不少教材和论文里有数据表可查，也有简化公式可用来做粗对照。
再就是极限板长验证：如果把d取得非常小（比如d/h=0.01），板几乎消失，Kt应该趋近于1、Kr趋近于0。如果程序在这个极限条件下算出来的Kt不是1，那说明板面边界条件施加方式有问题。反过来，把d取得接近h（d/h=0.99），Kt应该很小，而且接近底部固定式障碍物的已知结果。这两个极限一夹，中间区域的数值可信度就很高了。
我在验证时还会顺手画一张不同M值下Kt随周期变化的叠合图，检查是否有"毛刺"出现。如果曲线在某段频率出现抖动，多半是矩阵条件数恶化，先把M加大试试，不行就检查配点是否恰好落在z=-d的奇异点附近。
6. 源码使用建议与二次开发路线
6.1 拿到源码后第一步做什么
这套源码工程拿到手解压之后，按我的习惯会先做三件事：第一，跑一遍默认参数，确认能正常出结果；第二，把能量守恒残差打印出来，确认在1e-6量级；第三，把N改成1跑一遍单板工况，和理论值对照。这三步走完，才能确认当前机器上的Matlab版本和源码完全兼容。
参数修改集中在主程序的头部变量区。需要特别留意的是板间距S的输入格式：如果你的板不是均匀布置，而是每块板位置不同，源码里可能提供的是一个位置数组而不是单一S。我曾经遇到过用户把栅栏均匀间距S当成任意位置数组传进去，导致板位置错乱的情况。改代码前先看清楚变量定义，这是项目里最划算的注意力投入。
6.2 我遇到的三个高频报错及处理方法
第一个报错是色散方程求根函数报"没有找到足够多的特征根"。这通常不是算法问题，而是截断模态数M设置超过了当前水深和频率组合下可解析出的根数量范围。最简单的处理是把扫描区间上限加大，或者减少M。第二个报错是矩阵求解时出现NaN或Inf，多半是特征函数计算时cosh(kh)溢出，kh超过85左右双精度就危险了，需要改成cosh(k(z+h))/cosh(kh)的等价指数形式计算。第三个报错是绘图阶段Kt曲线出现负值或复数幅值异常，这种情况十有八九是匹配条件里某个指数项的符号写反，建议逐项检查exp(kₙ(x-x_j))和exp(-kₙ(x-x_j))的系数配对。
这里写一段我常用的稳定计算特征函数的技巧：
% 稳定计算cosh(k0*(z+h))/cosh(k0*h) % 避免大kh时溢出 Z0 = exp(k0*z) .* (1 + exp(-2*k0*(z+h))) ./ (1 + exp(-2*k0*h));这个写法在kh较大时数值表现远比直接cosh稳定。
6.3 扩展方向：别停在理想模型上
线性势流模型是无黏、无能量耗散的理想近似，实际工程里还存在波浪破碎、涡旋脱落的能量损失，这些都会让真实Kt比你算出来的更低。如果你后续想把模型做得更贴近实际，可以考虑几个扩展方向。
第一个方向是给薄板加"能量耗散"项，通过在匹配条件中加入一个与速度成比例的阻尼系数来近似模拟板面的粗糙度和波浪破碎损失，这样Kr² + Kt²会小于1，更符合实测结果。第二个方向是分析倾斜板，把板的几何位置换成斜面，匹配条件投影到斜面局部坐标系下就可以，特征函数展开框架完全能复用。第三个方向是接入柔性板或弹性板，此时板面条件变成线弹性振动方程，需要把结构方程和流体方程联立求解，复杂度上了一个台阶，但Matlab同样能处理。
我个人认为，对绝大多数工程预研场景来说，线性多板模型已经足够回答"挡浪效率大概是多少"这个核心问题。精细化的粘性修正，往往留给后续水槽实验或者CFD精确评估阶段再上。
回到开头那个甲方的问题，我用这套程序扫完参数后给出的建议是：四块板、板间距取周期对应的半波长附近、浸没深度比控制在0.6。这样在目标波谱主频附近Kt能压到0.3以下。对方按这个方案做初步设计，后续物模实验的实测结果也基本印证了数值趋势。这个项目让我印象最深的一点是：多板透射问题看起来不过是单板问题的重复叠加，实际算起来才知道干涉效应有多复杂——这也正是用数值方法预先探索参数空间的价值所在。
