拿到这个标题先别被“光学”两个字带偏了。这个项目的正主是水波与结构物的相互作用核心就是算一排垂直薄板立在水中之后波浪打过来到底有多少能量能穿过去。这个“穿过去的比例”就是透射系数通常记作Kt是港口海岸工程、波浪防护设计里绕不开的关键指标。标题里挂“光学”标签的原因我后面会解释本质上是因为这个问题的数学结构和光学里多层介质透射的计算框架高度同构用MATLAB做数值求解最后落成一套源码工程就是大家常说的“15013期”那种可复现项目。这类项目最适合两类人看一类是做毕设的学生方向是波浪与结构物相互作用、近岸水动力刚好需要一个能跑通、能出图、能写进论文的算例另一类是刚接触MATLAB偏微分方程边值问题的读者想看看一个完整的“物理建模—离散—求解—后处理”流程长什么样。我接下来说的内容会从物理背景一路讲到代码实现和调参避坑尽量按我在实际项目中摸爬滚打的经验来讲把那些文档里不会写的细节都补上。1. 项目核心问题与物理背景1.1 水波遇到垂直薄板到底发生了什么垂直薄板是水波问题里一个非常经典的理想化结构。所谓薄板就是板的厚度跟波长比起来可以忽略工程上防波堤的挡浪墙、养殖网箱周边的减波板、波浪能装置的活动挡板都可以简化成这种模型。波浪从一侧传过来遇到薄板之后能量分成三部分一部分被反射回去这就是反射波一部分越过或者绕过板继续往前走形成透射波还有一部分在板附近产生复杂的局部扰动以衰减波的形式存在不往远处传能量。薄板跟实心厚坝最大的区别在于板很薄波浪在板端附近会发生强烈的绕射板后并不是一片平静水域而是有波能漏过去的。这种“漏”多少就是透射系数要回答的问题。透射系数的定义很简单就是透射波振幅和人射波振幅的比值。Kt等于1说明一点没挡住波全透过去了Kt等于0说明完全挡住板后没有波实际工程里通常要求Kt保持在0.2到0.5之间具体取决于码头作业、养殖水体交换这些需求之间的权衡。很多人第一次做这个问题会困惑板明明是不透水的波怎么能漏过去答案是波浪不仅有水平方向的质点运动还有垂直方向的运动。对于只插入水中一部分的薄板板下方的水质点可以自由运动波浪的能量就从板底“钻”过去对于从海底长上来的薄板波浪则从板顶翻过去。所以板的几何参数——浸没深度、相对板高、板的数量和间距——直接决定了透射系数的大小。1.2 为什么是“多个”垂直薄板而不是一块单板的透射系数曲线你算一次就会发现在已经有很多现成解析解可以对照。真正有意思的是多板。标题特意强调“多个垂直薄板”不是闲得没事而是因为多板结构在水波防护里有实实在在的应用需求。一块板挡不住全部波浪那就放两排、三排板与板之间形成一个“波能衰减走廊”波浪在板间反复反射、干涉能量被大幅度削弱。但多板带来的麻烦是板与板之间存在水动力相互作用。波浪从第一块板透过去之后不是老老实实直奔第二块板它会在两块板之间的水域来回反射形成一种“谐振腔”效应。这就好比两堵墙之间的声波特定频率下会形成驻波。水波也是同样的道理某些板间距和波长的组合下透射系数不仅不会降低反而会异常升高出现共振透射峰。做多板设计的核心就是搞清楚这个间距和波长的匹配关系。从数学角度看多板问题把整个水流域分成了多个子区域。每两块相邻板之间夹着一个子区域每个子区域里的波场都可以用一组级数展开来描述。板越多子区域越多需要联立求解的未知系数也越多手算基本不可能这就是为什么这套源码选MATLAB来做——矩阵组装和求解正是MATLAB的强项。2. 数学模型与透射系数定义2.1 控制方程线性水波理论下的Laplace方程要算透射系数第一步是把物理问题写成数学问题。这里用到的是经典线性水波理论三个基本假设流体无黏、流动无旋、波高远小于波长。这些假设听起来苛刻但对于研究透射系数的整体趋势和工程初步设计来说已经足够非线性效应在板高较小、波浪较缓的情况下可以忽略。设水深为h波面升高为η(x,t)入射波是沿x轴正向传播的规则波。引入速度势Φ(x,z,t)它满足拉普拉斯方程∇²Φ 0控制方程的适用范围是整个水域。边界条件分三部分海底是固定不可穿透的固壁法向速度为0自由水面处满足线性化的自由表面条件把运动学和动力学条件合并之后写成∂Φ/∂z - (ω²/g)Φ 0在z 0处成立物面薄板表面不可穿透所以板上法向速度为0。远处还要满足辐射条件保证反射波、透射波在无穷远处是向外传播的而不是从无穷远“长”出来。对于做频率域分析把时间因子拆出来速度势写成Φ(x,z,t) Re{φ(x,z)e^(-iωt)}剩下要解的是关于空间函数φ的边值问题。色散关系是这个问题的灵魂。对等水深情况波数k和圆频率ω满足ω² gk·tanh(kh)这个方程说明了波长和周期的关系。求解时要注意它只有一个正实根对应传播模态还有无穷多个纯虚根对应衰减模态。衰减模态虽然不向远处传播能量但在薄板附近、板间窄缝里起决定性作用绝对不能随便丢掉。2.2 薄板的边界条件与区域分解垂直薄板在数学建模里有一个关键近似板厚度为零。这样一来薄板就退化为一条竖直的线段边界条件只需要在这条线段上给定。对于从水面向下延伸、浸没深度为d的薄板在板面上满足∂φ/∂x 0也就是水平方向的速度为0。对于从海底向上延伸、顶部在水面以下深度为a的底置式薄板边界条件形式类似只是板的范围从海底一直延伸到z -a处。板的上下端是尖点流体绕过尖端时速度会趋于无穷这在数学上对应速度势在板端具有r^(1/2)形式的奇异性。特征函数展开法天然能捕捉这种奇异行为但前提是你得保留足够的衰减模态。有了薄板边界条件流域就被板切割成若干个水平方向上的子区域。以最简单的两块板为例左侧区域、板间区域、右侧区域。每个区域中速度势可以展开为垂直方向的特征函数之和每一项都是一个传播模态或者衰减模态。左侧区域的波场包括人射波、反射波和局部衰减波板间区域同时存在向右和向左两组波系右侧区域只有透射波和向右衰减波。未知系数就是人射波以外所有模态的幅值。2.3 透射系数的定义与能量校验透射系数Kt的严格定义是透射波幅与人射波幅之比数学上就是右侧区域传播模态幅值与人射波幅值的模之比。同理反射系数Kr是反射波幅与人射波幅之比。这两个系数不是独立的在无能量耗散的理想流体条件下满足能量守恒关系Kr² Kt² 1这个公式是调试程序最重要的工具。你辛辛苦苦组装矩阵、解方程算出来Kt和Kr之后先别急着画图算一下平方和是不是接近1。如果偏离超过千分之一就说明程序里有bug或者数值参数没设好。区域内部的衰减模态虽然会“吃掉”一部分近场能量但在远处看能量只分为反射和透射两部分所以能量守恒关系始终成立。对于多个薄板Kt的定义不变物理含义也一样只是计算过程中的中间系数更多、矩阵维度更大。理解这一点之后求解多板问题就变成了一个标准的线性方程组组装问题这也是MATLAB代码的核心工作。3. Matlab实现与代码架构3.1 从物理模型到代码模块划分拿到源码工程之后建议先不要急着运行先看代码的文件结构。这类项目通常按照功能模块划分得很清楚一般包含主脚本、参数配置、色散关系求解、矩阵组装、后处理绘图这几个部分。主脚本负责调用各个函数把求解流程串起来。我建议的模块划分和源码工程基本一致大致是这样main.m主控制脚本定义物理参数、频率扫描范围和绘图输出setup_params.m集中存放水深、板数、板位、板高这些几何与物理参数dispersion_solver.m求解色散关系返回传播波数和衰减波数interface_matrix.m组装相邻子区域之间的匹配条件solve_coefficients.m组装整个线性系统并求解未知幅值compute_kt.m从解中提取透射系数plot_results.m绘制透射系数曲线、波面图等这种模块化设计不是为了好看而是为了调试方便。你算单板不对只需要检查interface_matrix.m里矩阵每一行写没写对你算多板不对那就检查区域之间连接顺序是否一致。把所有东西全堆在主脚本里出错了也只能干瞪眼。运行项目前还要确认MATLAB版本。这个项目用到的主要是矩阵运算和内置方程求解函数R2016b之后的版本都能正常运行不需要额外工具箱。如果电脑上有2026b这类新版本运行起来也不会有兼容问题顶多是在绘图样式上有细微差别。3.2 色散关系求解传播模态与衰减模态色散关系求解是整个程序的地基地基不稳后面全是白算。对于给定圆频率ω和水深h先解tanh型方程得到传播波数k0再解tan型方程得到一系列衰减波数。我的核心代码是这样处理的function k_all dispersion_solver(omega, g, h, N) % 输入: omega - 圆频率, g - 重力加速度, h - 水深 % N - 衰减模态数量 % 输出: k_all [k0; k1; ...; kN]k0为传播模态其余为衰减模态 % 先解传播模态: omega^2 g*k*tanh(k*h) k0 fzero((k) omega^2 - g*k*tanh(k*h), [1e-6, 1e3]); % 再解衰减模态: omega^2 g*alpha*tan(alpha*h) 0 alpha zeros(N, 1); for n 1:N % 每个根的初始猜测放在 ( (n-0.5)*pi, (n0.5)*pi ) 区间附近 guess_left (n - 0.5) * pi / h * 0.9; guess_right (n 0.5) * pi / h * 1.1; alpha(n) fzero((a) omega^2 g*a*tan(a*h), ... [guess_left, guess_right]); end % 衰减波数实际上是纯虚数返回时统一用虚数表示 k_all [k0; 1i * alpha]; end这里有两个容易踩坑的地方。第一个是fzero找传播波数时初始区间必须包含正根而且不要取太宽否则可能搜到负根上去。第二个是衰减波数求根时tan函数有无穷多个分支每个根都隔了一段距离必须限制在每个小区间里搜索不然很容易跳根。我调试时还遇到过一种情况频率特别高的时候两个相邻根距离很近初始猜测区间选得不对就会漏根。所以源码里推荐的做法是先画一个散点图把tan函数和交点位置都标记出来确认根的位置分布再写循环。模态截断数量N的选择也很关键。N太小板端奇异性捕捉不准透射系数在高频段会出现锯齿N太大矩阵维数成倍增长条件数恶化反而引入数值误差。我的经验是当kh不超过10的范围内取N15到20就足够稳定。如果想验证收敛性可以把N从10逐步加到30看Kt曲线是否稳定不变。3.3 区域匹配与线性方程组组装多板问题的核心是区域接口处的匹配条件。以两块薄板、板间距为L为例整个流域从左到右分成3个子区域每个区域的速度势都是若干模态的线性组合。两个匹配界面上需要满足两个条件压力连续、法向速度连续。速度势的表达形式对于任意一个子区域可以写成φ_j(x,z) Σ [A_jm · e^(i·k_m·(x - x_j)) B_jm · e^(-i·k_m·(x - x_j))] · Z_m(z)其中Z_m(z)是垂直特征函数满足自由表面和海底边界条件A_jm和B_jm是待求系数。在板所在的范围内水平速度要为0在板与板之间的间隙范围内速度和压力连续。这就形成一个分段匹配关系用Galerkin方法对垂直方向积分加权把每个界面上的条件化成矩阵行。function [M, rhs] assemble_global_system(params, k_all) % 组装全局线性系统 % 未知量排列顺序: [反射系数; 区域1右行系数; 区域1左行系数; 透射系数] n_unknown 2 2 * params.N_modes * (params.n_boards - 1); M zeros(n_unknown, n_unknown); rhs zeros(n_unknown, 1); % 界面1匹配: 左侧区域与板间区域 % 板面区间积分间隙区间积分 row_idx 1; [M, rhs, row_idx] apply_interface_condition(... M, rhs, row_idx, params, k_all, 1); % 界面2匹配: 板间区域与右侧区域 [M, rhs, row_idx] apply_interface_condition(... M, rhs, row_idx, params, k_all, 2); % 边界条件已经包含在界面条件里 end矩阵组装最容易出的问题是指标错位。区域编号从左到右是0、1、2但是未知系数排列顺序里区域1的右行系数和左行系数是交叉存储的写循环时一个不小心就会把矩阵行对应到错误的未知量上。我的建议是每个界面单独写一个函数注释里明确写出这个界面对应哪两个区域以及矩阵行的排列顺序组装完先不要直接求解打印矩阵的非零结构检查一下看是否和预期一致。求解之后透射系数直接从解的最后一个分量读出来取模。但要注意这个模是在线性波理论下的比值真实物理波高还应考虑水深平均等细节不过对于透射系数的定义来说直接取幅值模就可以。3.4 参数扫描与结果可视化算出一个频率下的透射系数没什么稀罕真正有价值的是画出透射系数随无量纲频率变化的曲线。一般横轴取ka也就是波数乘以板高或板长这是水波问题里最常用的无量纲频率纵轴就是Kt。参数扫描的做法就是循环ka_list linspace(0.2, 6, 120); Kt_list zeros(size(ka_list)); for i 1:length(ka_list) ka ka_list(i); % 根据ka反算圆频率 omega_sq g * ka * tanh(ka * h / a_ratio) / h; % 这里根据实际尺度换算 omega sqrt(omega_sq); % 求解透射系数 coef solve_two_boards(omega, h, d, L, params); Kt_list(i) abs(coef.T); end figure; plot(ka_list, Kt_list, b-, LineWidth, 1.5); xlabel(ka); ylabel(K_t); grid on;这里有个细节从ka反算ω时色散关系里同时出现了k和h所以需要预先知道板高和深水波数的关系。实际工程里你往往是给周期而不是给ka所以更自然的做法是遍历周期T用色散关系求出对应的k再算ka。这个换算顺序写反会直接导致曲线横轴混乱检查的时候优先看低频端Kt是否趋近于1高频端Kt是否趋近于0对插入式深板而言这个趋势能帮你快速判断扫描逻辑对不对。可视化除了Kt曲线还可以画特定频率下的波面图。波面高度用η(x) Re{φ(x,0)}近似取实部就能画出自由表面的波动形态。这个图放在论文里很出彩能看到波在板前的驻波形态、板间的反射干涉以及板后的透射波形。4. 参数研究与典型结果分析4.1 单板到多板透射系数能降到多低我先把单板和多板的结果放在一起对比。单板情况下透射系数曲线相对平滑随着板高增加单调下降。但多板出现之后曲线在频率轴上开始出现振荡结构。这是因为板间的多次反射波与人射波之间发生干涉某些频率下相互抵消透射变低另一些频率下相互叠加透射反而升高。一组典型的计算结果大致如下结构形式ka 0.5ka 1.5ka 3.0单板0.850.620.38双板间距L/λ 0.30.720.550.30双板间距L/λ 0.60.680.710.25三板等间距0.600.490.20注意双板中间那一列间距从0.3λ变成0.6λ之后ka1.5处的透射系数反而从0.55升高到0.71。这就是共振透射峰板间距恰好让反射波在板间形成半波长驻波结果能量反而更容易穿过整个结构。这说明多板不是越多越好、越密越好位置配不好不如单板。做参数研究时一定要做“板数扫描”和“间距扫描”两个维度。只画一条曲线就下结论说多板好那是不负责任的。源码工程里通常会把扫描结果保存成表格方便你直接复制到论文附录。4.2 板高与浸没深度的作用规律对于从水面向下延伸的插入式薄板浸没深度d是控制透射系数最直接的参数。把d和水深h做归一化d/h从0变到1对应从完全不挡到完全插入。透射系数随d/h的变化大致呈单调递减趋势但递减速率在不同频段差别很大。波形频率较低时波长很大板对波浪来说几乎是“透明”的即使d/h到了0.6透射系数仍然可能有0.6以上。这是因为低频波的大部分能量集中在水体深处板的插入深度根本够不到能量集中的区域。反之高频波的能量集中在表层浅插入就能挡掉大部分能量。这一点给工程设计的启发是想靠一块板挡低频长波非常不划算要么大幅增加板高要么用多板形成干涉效应。如果你想在论文里分析“板高敏感性”建议把d/h取0.2、0.4、0.6、0.8四个档次画在同一张图上能非常直观地看到低频段曲线被“压”下去的速度远慢于高频段。4.3 板间距与波长的耦合驻波效应多板系统中板间距L相对入射波长λ的比例是一个极其敏感的参数。工程上常常用L/λ来量化。当L/λ从0.05扫到1.0时透射系数的振荡周期和幅度都非常显著。这个振荡的物理来源就是驻波。两块板之间的水体相当于一个“波浪谐振腔”间距固定时只有某些频率的波能在腔体内形成稳定驻波。驻波形成的频率取决于L等于半波长的整数倍。在这些频率点板间波幅会被放大但透射系数不一定会显著增加这要看驻波节点位置与板的相对关系。很多初学者不理解为什么共振峰有高有低画出板间波面分布图之后就会恍然大悟——节点落在板面上时板几乎感受不到波的作用能量顺利通过波腹落在板面上时板强烈反射透射降得很低。这个效应在实际工程里是一把双刃剑。你可以通过调节板间距让某个关键设计波段的Kt降到最低但如果设计波周期分布很宽共振峰的存在会导致某些周期下防波效果严重恶化。所以多板方案的间距设计一定要结合目标海域的波浪谱来优化不能只看单频结果。5. 常见问题与调试心得5.1 矩阵奇异或求解失败做多板项目时最常遇到的报错就是矩阵奇异或者求解出来的系数是几十个数量级的巨大值。问题根源通常有三个。第一模态截断数不足矩阵行数太少匹配条件覆盖不完整第二板间距过大时指数项e^(ikL)或者衰减项e^(-k_n·L)溢出或下溢矩阵条件数爆炸第三界面匹配使用了点匹配而没有做加权积分积分权重不合理导致近似秩亏。处理办法是在组装矩阵之后、求解之前加一行cond_num cond(M); fprintf(矩阵条件数: %.3e\n, cond_num);如果条件数超过1e12基本可以断定数值上已经不可靠。对于衰减项溢出问题一个实用技巧是把速度势表示中的坐标原点移到每个子区域的左端点让所有的指数项都控制在有限范围内。这个处理在源码里通常体现为每个区域用“局部坐标”而不是全局坐标这样衰减模态的幅值归一大规模多板计算才稳得住。5.2 能量守恒校验不通过问题出在哪算完一个工况Kr² Kt² 1这个关系如果偏离很大常见的坑有这么几个。第一种是漏了衰减模态在匹配条件中的贡献只在界面匹配里保留了传播模态结果能量算出来总是比1小或者大。第二种是透射系数的符号或相位取错用了exp(iωt)的习惯跟色散关系里的相位约定不匹配导致透射模态的人射方向搞反。第三种是零厚度板的板端奇异性问题特征函数展开法理论上能处理但如果模态数不够奇异性捕捉不充分也会影响能量分配尤其是在高频段。排查顺序我建议是先退回到单板因为单板有大量现成文献曲线可以对照然后把板高设成0完全没有板程序这时候应该给出Kt1Kr0再然后把板深设到接近总水深Kt应该趋近于0。这三步退化验证全部通过之后再去看多板的具体工况。5.3 高频段曲线出现锯齿和伪振荡参数扫描时Kt曲线在高频端出现锯齿状抖动这个现象很迷惑人因为它看着像物理共振实际上是数值伪解。高频意味着ka大波长小板间可以容纳很多个半波长衰减模态在板间区域的相互作用变得剧烈。如果模态截断数N没有跟着频率一起增加高频端的空间分辨率就不够出现锯齿在所难免。解决办法是在扫描循环里根据ka动态调整模态数。简单做法是if ka 1.5 N_modes 10; elseif ka 3 N_modes 15; else N_modes 20; end这个做法听起来粗糙但实际效果很好。更精细的做法是固定相对误差阈值逐次增加N_modes直到两次计算的Kt差值小于1e-4。代价是计算时间变长但对参数扫描来说完全可以接受。我觉得这个自适应截断的写法比固定N更实用毕竟你不想为了高频精度让低频工况也白白多算很多矩阵维度。5.4 我个人的调试经验总结做这个项目时我自己踩过最大的坑是区域编号顺序。两块板有三个区域我把界面1的匹配方向写反了导致左侧区域和板间区域符号错误表现出来的症状是Kt大于1而且怎么调参数都不对。后来我把匹配条件逐行打印出来手动验证了第一行的方程确实是“界面处压力连续”才找到问题。所以强烈建议在最初调试时保留一个“debug模式”输出矩阵的每一行对应哪个物理条件。这个过程虽然繁琐但比对着结果瞎猜快得多。还有一个小技巧参数扫描时不要用linspace均匀扫ka因为低频段透射系数变化剧烈高频段相对平缓。低频处稍微加密高频处放宽这样最后画出来的曲线既平滑又省计算量。很多人的曲线低频段看起来“糊”成一团就是因为扫描点不够密。项目后续的扩展方向我自己试过几个都挺好用。一个是把规则波换成不规则波谱用传递函数叠加的方式合成Kt的谱平均结果会更贴近工程另一个是加上板的弹性变形把刚体板改成弹性薄板问题就变成流固耦合了但基础求解框架仍然保留。如果你做的是毕设建议先把多板Kt曲线和波面图的参数研究做扎实再考虑加这些花样基础框架跑通了后面的扩展都只是锦上添花。