简介：面向海岸工程与波浪数值模拟研究者的SWASH模型孤立波爬升算例包，以锥形岛上的波浪爬升（Runup of solitary waves on a conical island）为对象，提供一套自制、可直接运行的输入文件与辅助脚本，适合需要复现该经典算例或入门SWASH建模的科研人员和学生。压缩包共包含5个文件，其中3个.sws格式配置文件用于设定不同工况下的模拟参数，另2个Matlab脚本分别负责生成地形（.bot）与边界（.bnd）输入文件以及模拟结果的可视化绘图。资源包大小仅5KB，轻量且结构清晰，能够省去手工编写繁琐输入文件的步骤，从地形边界生成、模型配置到后处理形成完整流程。通过该算例可掌握SWASH中自定义地形与边界文件的生成方法，理解孤立波在复杂地形上的爬升机制。目前已有301人学习下载，对从事相关模拟研究或课程学习具有直接参考价值。
1. 孤立波撞锥形岛：这份SWASH算例输入文件解决什么问题
SWASH模型在近岸水动力圈子里早就不算陌生，但真正能把孤立波打在锥形岛上、完整跑出爬高过程的算例文件却很难找。这份输入文件包就是把锥形岛地形、孤立波边界条件和测点配置全部写成可直接提交的文本，下载后配上SWASH源码就能算，省掉你自己从头搭地形、调边界条件的那几天。它适合两类人：刚接触SWASH、想找一个能跑通的数值实验当学习骨架的入门者；做海啸或风暴潮爬高研究、需要一套能改波高、改水深、改岛形的基础脚本的从业者。第一次跑通它，你得到的不是一张示意图，而是一组可量化的爬高过程数据。
2. 从选型到设定：为什么SWASH适合算孤立波爬高
2.1 可选的模式不少，为什么偏选SWASH
我最早算爬高时，第一个想到的是非线性浅水方程。浅水方程跑孤立波很顺，波峰一上岛就“啪”地铺开，干湿边界处理也快。但问题在于浅水方程没有频散项，孤立波在深水区走着走着波形就胖了，甚至生出讨厌的数值尾波。后来试过XBeach，它对近岸破碎和海岸线移动模拟很成熟，不过XBeach的默认模式是NLSW加破碎耗散，用在锥形岛这种“波从深水斜向登陆”的场景上，对水深地形和入射波长的适配反而没有Boussinesq型模式直观。
SWASH最核心的价值在于它支持扩展Boussinesq模式和非静压格式，能在保证计算效率的前提下把孤立波的频散关系描述得更准。孤立波本质上是一个“波形稳定、速度依赖波高”的浅水非线性波，它能不能在传播过程中保持形状，很大程度上取决于模式有没有频散项参与平衡。纯浅水方程没有平衡机制，波越传越失真；而用完全非静压模式跑，精度高但计算量偏大。SWASH允许你在同一个框架里切换模式，我先用Boussinesq跑工程级算例，需要校核局部细节时再换非静压，这在做爬高研究时非常顺手。
| 模式类型 | 频散能力 | 干湿边界 | 适合场景 |
|---|---|---|---|
| 非线性浅水方程 | 无频散，波形易失真 | 处理成熟 | 长波淹没、漫滩宏观过程 |
| Boussinesq型 | 弱频散，适合孤立波 | 较成熟 | 孤立波传播、爬高、岛礁绕射 |
| 非静压模式 | 频散覆盖广 | 较成熟 | 近岛复杂绕射与破碎细节 |
这个表不是要证明SWASH全面胜出，而是说明：做锥形岛孤立波爬高，你真正需要的是“保留频散、兼顾干湿边界、还能控制计算成本”的模式，SWASH正好落在区间里。
2.2 算例的物理场景与几何参数
这个算例的物理场景很清晰：一个孤立波从深水区生成，向右传播，正面撞上位于计算域中的锥形岛，然后在岛屿四周发生绕射和爬高。锥形岛的形状和尺寸完全按那组经典物理实验的比例关系重建，我拿到文件后仔细对过一遍，几何参数大概是这样的：
- 锥形岛底半径接近110m，顶部是一个直径约22.5m的平顶；
- 岛高接近62.5m，而四周静水深约32m，也就是说岛屿顶部明显出露水面；
- 计算域水平范围约640m×448m，左侧边界作为波浪入射边界；
- 入射孤立波波高约1.2m，周期特征由孤立波解析解确定。
这里要强调一点，“波高1.2m、水深32m”并不是拍脑袋定的。岛屿周围的平均水深决定孤立波的特征波长，而特征波长又决定锥形岛会不会被完全包围、背浪侧会不会出现明显的波聚焦区。波高与水深之比大概在0.04左右，属于弱非线性到中等非线性的范围，如果用线性波理论去近似，几乎复现不出背浪侧的爬高差异。所以这份算例选这个波高比例是有讲究的，实际SPH模拟时这个比例也是控制精度的关键。
文件中这些参数都不是写死的。我把2.2节列出的数值称为“基准工况”，你要复现论文量级或做敏感性分析，改起来非常直接：波高在主控文件的边界条件里改，水深在地形文件里改，网格间距在网格段落里改。唯一不建议乱动的是锥形岛的坡度，因为坡度的突变会直接影响爬高过程是否发生破碎，那是另一个物理问题。
2.3 输入文件的结构与三大模块
下载解压后你会看到一组文本文件，用任何编辑器都能打开。整个算例的逻辑可以分成三块：主控制文件负责描述计算域、网格、物理模式和输出请求；地形文件负责给每个网格点赋水深或高程；边界文件负责把孤立波“喂”进计算域。三者缺一不可，下面是我建议的阅读顺序：
| 文件 | 作用 | 需要关注的内容 |
|---|---|---|
| 主控制文件（.sws） | 整体调度，模式、网格、时间步、输出全在这里 | MODE、GRID、BOTTOM、WAVEBR、OUTPUT |
| 地形文件（.dep/.txt） | 每个网格点的水深值 | 坐标顺序、水深负值含义、岛顶高程处理 |
| 边界波文件（.bnd） | 入射孤立波的时间序列 | 时间步间距、波高、静水深 |
这三个模块之间的关系是：SWASH先读主控制文件，按网格段落生成计算网格；随后读地形文件，给每个网格点赋值水深；最后在时间推进过程中，从边界波文件里逐时刻读取入射条件，驱动孤立波向岛传播。主控制文件是大脑，地形文件是骨架，边界波文件是心脏，任何一个出错整个算例都会“心肌梗死”。
3. 拆读输入文件：网格、地形与边界条件的落地写法
3.1 计算域网格与水深：先定网格，再定时间步长
打开主控制文件，你最先看到的是计算域与网格定义。这个算例的计算域是640m×448m，网格分辨率为4m左右，这意味着横纵向网格数大约在160×112的量级，对单机并行来说压力不大。之所以选择4m而不是更粗的8m，是因为锥形岛底半径110m，一个坡面上至少要保证20到30个网格点才能把坡度描述清楚；如果用到8m，岛坡会变成“锯齿坡”，爬高过程容易在个别网格点上出现虚假堆积。
COORDINATES CARTESIAN GRID 0.0 0.0 640.0 448.0 4.0 4.0 SET LEVEL 32.0这里第一行限定坐标系为直角坐标，第二行的前四个数字是计算域左下角和右上角的坐标，后两个数字是x和y方向的网格间距。SET LEVEL定义初始静水面高程，算例里设为32m。要注意这个32m是相对基准面的水深，后续地形文件里所有数值都要以同一基准面为参考，否则地形和静水面会错位，初始时刻就出现“凭空的水位差”。
网格定完之后，时间步长不要自己瞎拍。SWASH一般用默认CFL限制自动调整时间步，但算孤立波时我会额外关注一个参数：时间步长的最大上限。算例文件里已经设置了细化的时间步控制，你能看到类似下面的配置：
TIME STEP 0.005 3600.0 CFL 0.5前一个数字是最小时间步，中间是最大时间步，最后一个数字是CFL目标值。0.5表示每个时间步内波面信息最多走过半个网格，这个值对孤立波来说是安全的。如果你想加速，先改CFL到0.7试跑，不要一上来就硬压最大时间步，否则遇到干湿边界振荡时容易直接发散。
3.2 锥形岛地形文件：从平面坐标到水深负值
地形文件是整个算例里最容易看错的部分。SWASH的地形文件不是“高程文件”，而是“水深文件”：计算域内每个网格点存一个数值，数值代表该点相对静水面的水深。静水面以下为正，静水面以上为负，所以锥形岛顶部那些出露的网格点，存的是负值，比如-30.5m表示这个点高出静水面30.5m。如果你把岛顶高程写成正的30.5，SWASH会认为那里是一个30.5m深的深坑，孤立波直接掉进去，结果完全失去物理意义。
算例自带的地形文件已经处理好了这一点，但如果你想自己重建一份锥形岛地形，可以用下面这段Python脚本作为参考。它不是算例自带的脚本，只是我平时重建地形的惯用做法：
import numpy as np dx = 4.0 x = np.arange(0.0, 640.0, dx) y = np.arange(0.0, 448.0, dx) X, Y = np.meshgrid(x, y) xc, yc = 320.0, 224.0 R_bottom = 110.0 R_top = 11.25 H_island = 30.5 h_water = 32.0 # 坡面半径随高度线性收缩 r = np.sqrt((X - xc)**2 + (Y - yc)**2) depth = np.full_like(r, h_water) # 底半径范围内的点才参与岛体构建 mask_bottom = r <= R_bottom slope = (R_bottom - R_top) / H_island # 单位高差对应的半径变化 h_above_water = (R_bottom - r[mask_bottom]) / slope depth[mask_bottom] = h_water - h_above_water # 超出静水面的点存负水深 depth[depth < 0] = depth[depth < 0] * -1.0 # 保存为SWASH可读的文本矩阵 np.savetxt("island_bathymetry.txt", depth, fmt="%.2f")这段代码的逻辑是：先铺平面网格，然后按圆锥体几何关系计算每个点到岛中心的距离，落在底半径内的点根据距离反推该点的高程；静水面以上的点强行转成负值。执行后会生成一个与网格一一对应的水深矩阵，用文本格式保存即可。要注意的是，H_island我这里写的是出露高度，不是岛的总高；算例原文件里总高约62.5m，出露高度约30.5m，这两者不要混淆。如果你改了网格间距，脚本里的dx要同步改，否则网格点数对不上，SWASH读地形时会报错。
3.3 边界条件与孤立波生成：时间序列文件不是玄学
边界条件是这份算例的精华。孤立波不能像规则波那样用周期公式简单描述，它本质上是“只有一次波峰”的瞬态波，所以最稳妥的生成方式是在入射边界上给定一个随时间变化的波面序列。算例文件里用的是左侧边界读时间序列的方式，主控制文件里的边界段落大概长这样：
BOU SIDE LEFT CCW WAVEBR 0.0 120.0 FILE 'solitary_series.bnd'含义是：左侧边界从第0秒开始，持续120秒，从solitary_series.bnd文件里读取波面数据。后面还有一行是启用边界吸波，用来吸收反射波，防止从右侧边界再弹回来干扰岛前波场。这个吸收配置不建议删除，孤立波打到岛上之后必然产生反射波，如果没有吸收，反射波会在计算域里来回振荡，最后的爬高数值就不可信了。
solitary_series.bnd文件的格式是一列时间、一列波面位移。它怎么生成？算例作者已经按孤立波解析解算好了。如果你要换波高，用下面这段逻辑重新生成：
import numpy as np H = 1.2 # 波高，单位m h = 32.0 # 静水深度 g = 9.81 c = np.sqrt(g * (h + H)) # 孤立波波速近似 t = np.linspace(0, 120, 12001) # 孤立波波面解析解 eta = H / np.cosh(np.sqrt(3.0 * H / (4.0 * h**3)) * (c * t - 320.0))**2 # 写边界文件 with open("solitary_series.bnd", "w") as f: for i in range(len(t)): f.write(f"{t[i]:.4f} {eta[i]:.6f}\n")这里c * t - 320.0中的320是让波峰在初始时刻位于x=320m附近，也就是锥形岛前一小段距离，给波一个短暂的传播预热期，避免初始时刻波峰直接压在岛坡上导致瞬时冲击。你可以根据想要的入射位置改这个数字。时间序列的采样间隔我习惯取0.01s，既不会让文件太大，也能保证SWASH在读取边界时有足够分辨率。采样太疏会让入射波形变得“棱角分明”，产生非物理的高频振荡。
4. 复现运行：编译、并行与后处理
4.1 编译SWASH前需要确认的三件事
这个算例包只包含输入文件，不包含SWASH源码本身，所以你拿到手的第一步是先确保本机computing环境到位。编译SWASH之前，我建议按这个顺序检查：
一是Fortran编译器。SWASH源码是用Fortran写的，本机至少要有一个可用的gfortran或者Intel编译器。Linux下检查命令是gfortran --version，macOS下如果没装，用Homebrew补一个。这一步没搞定，后面一切免谈。
二是MPI库。算例虽然不大，但SWASH的并行能力是标配，建议装OpenMPI或MPICH。检查命令是mpirun --version。如果你只想串行跑，可以跳过MPI，但后处理时你会发现算例稍微改细一点串行就很吃力，所以还是建议一次装到位。
三是NetCDF库。SWASH后处理输出经常需要NetCDF格式，编译时如果缺NetCDF，某些输出选项会被禁用。我的经验是，先装上libnetcdf-dev这类开发包再configure，省得后面为了一个输出格式重新编译一遍。
三件事确认完，按标准流程编译：
cd swash-src export FCFLAGS="-O2" ./configure --enable-mpi make sudo make install这里的-O2是优化级别，对绝大多数算例足够。--enable-mpi会启用MPI并行支持。如果你不需要NetCDF，可以在configure时加--disable-netcdf，但我不建议这么做，因为算例里的输出请求可能依赖NetCDF格式。
4.2 运行算例与输出控制：别让输出文件撑爆磁盘
编译完成后，进入算例目录，直接执行：
cd /path/to/runup_conical_island mpirun -np 4 swash runup_conical.island.sws-np 4表示用4个进程并行。先解释一下这里的进程数与网格的对应关系：算例网格是160×112左右，4个进程在二维区域分割时每个进程大约承担40×112个网格点，负载基本平衡。如果你用的是笔记本电脑，2个进程也能跑，只是速度慢一些。核心数超过8时我反而不建议盲目加进程，因为网格总量不大，通信开销占比会明显上升。
算例运行时会实时输出时间步和最大水位变化。看到STOP字样不代表出错，那可能是正常终止；真正需要警惕的是NaN或Inf出现在水位输出里，那说明计算已经发散。
输出控制方面，这个算例做了比较克制的配置：只在若干测点记录波面时间序列，同时用块输出保存每间隔一段时间的全场水位。你要特别注意块输出的时间间隔，如果设置成每0.1s保存一次全场，跑完120s会生成上千个全场文件，磁盘很快被塞满。我一般先看测点结果，确认波形合理再开全场输出。算例主控文件里的输出段落大致是：
OUTPUT 0.5 360.0 BLOCK 'T1' 'WATM' XP 0.0 640.0 YP 0.0 448.0 POINTS 'P1' 320.0 224.0 ... PTS 'P1' 'T1'第一行表示每0.5s写一次结果，最大输出时长360s。BLOCK行请求全场水位快照，POINTS行定义测点坐标，PTS行把测点数据输出到文件。如果你只关注爬高，测点输出才是核心，全场块输出可以在排错时再打开。
4.3 从输出文件中提取爬高过程
SWASH的测点输出文件名一般是P1.tbl之类的文本表格，第一列时间，第二列水位。用Python可以快速提取最大爬高：
import numpy as np import glob for pt_file in ["P1.tbl", "P2.tbl", "P3.tbl"]: data = np.loadtxt(pt_file, skiprows=1) t = data[:, 0] eta = data[:, 1] max_eta = np.max(eta) max_t = t[np.argmax(eta)] print(f"{pt_file}: 最大水位 {max_eta:.3f} m, 出现在 {max_t:.2f} s")这段代码的skiprows=1是因为SWASH输出文件首行是说明文字，实际文件中可能有两行表头，你根据打开文件后看到的情况灵活调整。计算爬高时要注意一个物理坑：测点如果布置在静水面以下的近岸区域，最大水位就是单纯的水面波动；如果测点布置在岛坡上且大部分时间处于暴露状态，测点输出的水位会在干湿边界切换时出现尖刺，这种尖刺不能算爬高。真正合理的爬高提取方法是让测点覆盖整个坡面，取每个测点在波到达后的最高水位，再沿坡面做最大值的包络线。
算例里测点布置在岛前、岛顶、背浪侧等关键位置，你第一次跑完直接看这几个点的最大水位，就能判断孤立波是否在背浪侧出现了明显的“二次爬高”特征。如果背浪侧测点的最大水位明显高于入射波高，说明绕射波在背浪侧发生了聚焦，这是锥形岛算例最有趣的现象。
5. 避坑与常见问题排查：这个算例最容易翻车的五个地方
5.1 编译与运行阶段的崩溃
现象1：编译时找不到NetCDF头文件，configure直接失败。
原因：NetCDF开发包装了但路径没暴露给configure。解决：先查NetCDF装在哪，通常是/usr/lib/x86_64-linux-gnu，然后设置环境变量export NETCDF_ROOT=/usr，再重新configure。如果你确认自己用不到NetCDF输出，也可以干脆--disable-netcdf，但算例里的某些输出配置可能失效。
现象2：运行到前几个时间步就崩溃，提示网格读取失败。
原因：地形文件的网格点数与主控制文件里的计算域、网格间距对不上。解决：检查地形文件的行数和列数，确保它等于Lx/dx + 1和Ly/dy + 1。我遇到过最笨的错是地形文件里第一行是说明文字，SWASH把它当成了水深数据，导致整个地形偏移一个网格。方案是用Python读文件时先看前五行，有说明文字就删掉。
现象3：算到40秒左右水位出现NaN，计算还能继续但结果全是乱码。
原因：干湿边界处网格被反复淹没问题。SWASH处理干湿边界时，如果某网格水深接近零且时间步长偏大，可能出现负水深振荡，最终爆炸。解决：把CFL从0.5降到0.3，或者把最小时间步再调小一档。这属于“杀敌一千自损八百”的招，计算时间会变长，但对干湿边界的稳定性很有效。
5.2 地形、边界与后处理阶段的问题
现象4：孤立波的波峰走到岛中央就停下来，波形严重变形。
原因：入射波时间序列的波速近似公式用得太粗糙。孤立波波速不是简单的sqrt(gh)，而是与波高有关，算例里用了sqrt(g(h+H))的近似。如果你把波高改得很大，这个近似误差会变大，波形在传播时位置会漂移。解决：用更完整的孤立波理论公式重新生成时间序列，或者让时间序列里的波峰初始位置再往左移一段距离，给波更多的调整空间。
现象5：背浪侧测点最大爬高比入射波高还大两倍，怎么看都不对。
原因：测点所处网格在退水过程中被输出了虚假水位。干湿边界处理的经典问题是网格从“干”转“湿”的瞬间会有一个数值冲击。解决：不要只看单点最大水位，把测点时间序列画出来，观察水位在波峰过后是否立刻回落；如果在回落之后还有一个尖锐峰值，那大概率是干湿边界尖刺。把测点往岸边方向移半个网格到1个网格，重新提取，通常能过滤掉这个假信号。
6. 进阶验证：网格敏感性测试与“后悔药”式工作流
第5章解决了跑通的问题，但跑通不等于算准。我拿到这个算例后做的第一件进阶工作是网格敏感性测试。做法很简单：在基准网格基础上生成两套衍生地形，一套网格加密一倍到2m，另一套放粗到8m，三套分别跑完同样的120s。对比三个结果中背浪侧最大爬高的差异：如果2m与4m的差异在5%以内，说明结果已基本收敛；如果差异还在15%以上，这4m网格就不够，需要继续加密。这个对比应当画成一张“网格分辨率—爬高值”的收敛曲线，而不是只比一个数字，爬高出现位置同样会随网格变化。
第二件常做的是时间步“后悔药”式验证。孤立波算例里，网格收敛不代表时间步收敛。实际做法是把CFL从0.5调到0.25，重跑一次，对比两次波峰到达测点的时间差。时间差如果超过0.5s，说明当前时间步对波峰追踪不够精确。对爬高研究来说，波峰到达时间偏差0.5s可能意味着与实测数据错开整整一个相位，直接影响结论可信度。
从那以后，我每次改完波高或地形参数，都强制走一遍“备份原始输入文件、生成新工况、粗网格快速核验、细网格正式计算”的流程。这套脚本化的流程不是为了让别人觉得我严谨，而是因为孤立波爬高这种算例太容易出现“看起来没问题但实际差一截”的情况。如果你也打算拿这个算例做工程预测或论文验证，建议从一开始就把这个习惯养成。希望帮到你。
本文还有配套的精品资源，点击获取
单调栈解法)
