简介：针对SWASH模型中的经典算例Runup of solitary waves on a conical island，提供一套自制输入文件与配套脚本，适合海岸工程、近岸水动力学方向的研究生和工程师用于复现孤立波在锥形岛上的爬高过程。压缩包共5个文件，包含3个.sws格式的SWASH配置文件，分别对应不同数值实验设置，以及2个Matlab脚本，一个用于生成.bot地形文件和.bnd边界文件，另一个用于结果可视化与后处理。资源包仅5KB，轻量紧凑，便于直接下载与修改。目前已有301人浏览学习。通过这套文件，读者可免去手工编写复杂输入文件的麻烦，快速搭建算例，同时结合脚本理解SWASH模型的配置逻辑、孤立波生成方式及边界文件构造方法，是学习SWASH模型和波浪爬高模拟的高效起点。
1. 锥形岛孤立波爬高的 SWASH 算例：为什么很多人卡在“没有输入文件”
SWASH 模型做近岸波浪爬高时，孤立波遇锥形岛是一个经典到不能再经典的验证场景：轴对称岛体、复杂的干湿边界、波前聚焦和绕射同时发生，数据敏感度比普通斜坡算例高一个量级。可真正动手时你会发现，官方示例里水槽、斜坡、恒定流一应俱全，唯独这种三维岛体算例的输入文件找不到现成的，只能自己造。这也是“算例输入文件（自制）”这个标题的由来。我当初也是从零开始拼出一套可用的输入文件，包括底部地形、孤立波初始场、开放边界和测站输出。这篇笔记把整套流程拆开讲，适合正在用 SWASH 复算海啸爬高、风暴潮越浪或锥形岛实验的工程师和研究生。
2. 算例输入文件的底牌：SWASH 到底需要哪几个文件
2.1 四类文件各管一摊
SWASH 的算例从来不是一个文件就能跑起来的，它由主控文件加若干外部数据文件组成。自制输入文件的第一件事，就是把这几张牌认全。
| 文件 | 内容 | 在这个算例中的角色 |
|---|---|---|
| 主控 INPUT 文件 | 网格参数、模式、摩擦、边界、输出设置 | 指挥中心，决定模型怎么跑 |
| 底部 BOTTOM 文件 | 每个网格点的平面坐标与底高程 | 告诉模型地形是锥形岛而不是平底 |
| 初始孤立波场 | 水面高度与流速的空间分布 | 把孤立波放进计算域，让它自己向岛传播 |
| 测站输出设置 | 观测点坐标与输出变量 | 提取爬高过程线，后续分析全靠它 |
常见的误区是把注意力全放在 INPUT 主文件上，结果跑起来才发现地形文件不匹配，或者初始场没赋速度，一启动模型就发散。其实自制算例的关键顺序应该是：先定几何，再定网格，然后做地形文件，最后才是边界和初始条件。
2.2 几何参数怎么落到数值模型里
“锥形岛”在数值模型里并不是一个尖锐的圆锥，而是一个顶部切平的锥台，否则平台上的干湿交替会带来不必要的数值噪声。控制几何的参数有五个：
- 底部半径 r_base
- 顶部平台半径 r_top
- 岛顶高程 z_top
- 海床高程 z_bed
- 静水深 h0
外加一个运动参数：孤立波波高 A。这些值可以直接取物理实验原型尺寸，也可以按自己的场景改。以我常用的典型工况为例：计算域 80m × 80m，网格间距 0.25m，岛底半径 30m，平台半径 12m，岛顶高程 0.5m，海床高程 0.0m，静水位 0.4m，孤立波高 0.1m。
这个取值里有几个隐藏约束。岛顶高程必须高于静水位，岛屿顶部才能露出来，爬高过程才有意义；波高又必须小于岛顶出露高度，否则整个岛都被淹没，测到的是漫顶水流而不是爬高。第一次做的时候建议按这个比例来，后面再根据兴趣调整。
2.3 目录组织与命名习惯
SWASH 输出文件很多，算例目录我一般这样组织：
conical_island/ make_bottom.py make_input.py run_swash.sh bottom.txt island.inp work/文件名全部用小写英文，不要带中文路径。某个开发者曾经把算例放在带空格的路径底下，SWASH 读取地形文件时路径解析出错，排查了一个小时才发现是空格的问题。work 目录单独放 SWASH 的输出，和输入文件隔离，避免中间文件覆盖掉原始配置。这套习惯看起来琐碎，但当你同时跑三组网格收敛性对比时，能省下的时间非常可观。
3. 用 Python 生成锥形岛地形文件：公式、脚本与检查
3.1 截锥地形的几何定义
地形生成的实质，是给计算域里每一个网格点算一个底高程 z。对锥形岛，我按径向距离判断三个区域：
- 当半径小于等于 r_top，底高程等于 z_top
- 当半径在 r_top 和 r_base 之间，底高程沿径向线性下降
- 当半径大于 r_base，底高程回到海床值 z_bed
这种分段函数用 numpy 的 where 组合写起来最顺手。网格用 meshgrid 生成，注意坐标系要让岛屿中心正好落在计算域中心，否则边界条件不对称，后续分析会很痛苦。
3.2 底部文件生成脚本
下面的脚本会输出 SWASH 可读的 bottom.txt，一行代表一个网格点，三列分别是 x、y、底高程。
import numpy as np # 计算域与网格 Lx, Ly = 80.0, 80.0 dx = dy = 0.25 nx, ny = int(Lx / dx) + 1, int(Ly / dy) + 1 # 几何参数 x0, y0 = Lx / 2, Ly / 2 r_base = 30.0 r_top = 12.0 z_top = 0.5 z_bed = 0.0 # 生成平面坐标网格 x = np.linspace(0, Lx, nx) y = np.linspace(0, Ly, ny) X, Y = np.meshgrid(x, y) R = np.sqrt((X - x0) ** 2 + (Y - y0) ** 2) # 分三段给底高程 Z = np.full_like(R, z_bed) slope = (z_top - z_bed) / (r_base - r_top) Z = np.where(R <= r_top, z_top, Z) Z = np.where((R > r_top) & (R <= r_base), z_bed + slope * (r_base - R), Z) # 写出 SWASH 底部文件 with open("bottom.txt", "w") as f: for iy in range(ny): for ix in range(nx): f.write(f"{X[iy, ix]:10.3f} {Y[iy, ix]:10.3f} {Z[iy, ix]:10.3f}\n")这里有两个细节值得展开。
第一个细节是写入顺序。SWASH 读取底部文件时，通常按行扫描平面坐标，先固定一个 y，x 从左到右变化，然后 y 再往上推进。上面代码外层循环是 iy，内层是 ix，正好匹配这个约定。如果两层循环写反，地形文件里就会出现“先固定 x 再变 y”的排列，模拟出来的岛屿形状会被转置，看起来像一条山脊。
第二个细节是浮点格式。用{:10.3f}保证每个数字占至少 10 个字符宽，列对齐之后用文本编辑器检查文件时，扫一眼就能看出坐标是否有跳变。文件大小也不用担心，321×321 也就是十万行左右，文本文件几 MB 级别，SWASH 读取没有压力。
3.3 运行前的地形“自检”
地形文件写出来之后，先别急着丢给 SWASH，花几十秒做一个静态检查。下面这组命令能快速找出高程异常：
# 查看高程最小值和最大值 sort -k3 -n bottom.txt | head -1 sort -k3 -n bottom.txt | tail -1 # 统计初始时刻的干区网格数（底高程高于静水位） awk '$3 > 0.4' bottom.txt | wc -l最大高程应该是 0.5，对应平台顶；最小高程应该是 0.0，对应外围海床。干区网格数可以手工估算：平台面积约 452m²，加上周边锥面上高于 0.4m 的环形带面积，除以单个网格面积 0.0625m²，数量级应该和 awk 统计结果对得上。对不上就说明几何公式或写入顺序有问题，这时候返工最划算，等进了 SWASH 再查就费时了。
# 若地形明显不对，先打印前 5 行确认行列顺序 head -5 bottom.txt前五行的 x 应该依次是 0、0.25、0.5、0.75、1.0，而 y 保持为 0。如果看到 x 不变、y 在变，行列顺序就是反的，改循环外层即可。
4. 组装 INPUT 文件：非静压、孤立波初始场与必调参数
4.1 INPUT 文件的主干结构
SWASH 的 INPUT 文件是按块驱动的文本，每个块承担一个职责。我的锥形岛算例主控件结构如下：
#!/usr/bin/env python3 # 生成锥形岛算例的 SWASH 输入文件框架 def write_input(): lines = [ "PROJECT 'island_runup'", "MODE 'NH'", "SETIONS", " CGRID S", " REGULAR 321 321 0.25 0.25", " DXY 0.25 0.25", " OFFSET 0.0 0.0", "BOTM BOTTOM", "FRIC 0.013", "DRY 0.005", "BOU T 1", " WLEVEL 0.40", "BOU T 2", " WLEVEL 0.40", "INI H 0.40", "OUT", " # 按你的版本手册填入测站坐标和输出变量", "TIM", " TUNIT 'SEC'", " STOP 60.0", "STOP", ] with open("island.inp", "w") as f: f.write("\n".join(lines) + "\n") if __name__ == "__main__": write_input()逻辑上要抓住三个关键点。第一是 MODE 'NH'，非静压模式。孤立波在传播过程中垂向加速度不可忽略，静压模式会把波形逐步抹平，爬高结果明显偏低。第二是 BOTM BOTTOM，它告诉 SWASH 去读当前目录下的 bottom.txt。第三是 DRY 0.005，这个值控制干湿判据，后面第五章会单独展开。
需要提醒一句：SWASH 各版本对关键字的大小写和缩写习惯不完全相同，上面给出的是一种常见写法。你拿到手应该做的是对照自己版本的参考手册，逐个块确认，而不是直接复制运行。自制算例本来就是“以你的模型版本为准”。
4.2 孤立波初始场：只抬水位会让模型“散架”
本算例推荐用初始场方式注入孤立波，而不是在边界上输入时间序列。原因是锥形岛周围有绕射和反射，边界输入会引入不必要的数值噪声；把孤立波整体放进计算域，让它自己朝岛传播，物理过程更干净。
孤立波的水面可以用解析解描述：
η(x) = A * sech^2( sqrt(3A / (4h^3)) * (x - x0) )其中 A 是波高，h 是静水深，x0 是初始波峰位置。在典型工况里，A=0.1m，h=0.4m，x0 取在岛中心上游 20m 处，这样孤立波有足够的空间完成自我整形，到达岛前保持稳定波形。
这里最常见的问题是只给水位，不给速度场。孤立波不是一坨静止的水包，它内部有特定的水平速度分布，近似关系为：
u(x) = c * η / (h + η) c ≈ sqrt(g * (h + A))如果初始时刻只抬水位，流速为零，重力会立刻让水包向下坍塌，一个孤立波会分裂成两个波，还会伴随寄生振荡。这是模型发散的经典元凶。把 η 和 u 同时写进初始场，波才能在传播过程中保持形状。
4.3 必调参数初值与边界处理
几个参数在第一次跑的时候建议直接用表里的初值：
| 参数 | 初值 | 调整方向 |
|---|---|---|
| 网格间距 | 0.25m | 岛附近加密，远场可放宽 |
| 最大时间步 | 0.05s | CFL 超 0.5 时减小 |
| 干湿阈值 DRY | 0.005m | 数值振荡时增大，爬高偏低时减小 |
| 底摩擦系数 | 0.013 | 模拟真实海床时调大 |
| 模拟时长 | 60s | 等爬高峰值过后 5s 再停 |
计算域的东西两侧边界都设成H水面为 0.4m 的开边界，南北两侧我一般会加吸收层。侧向边界如果没有吸收，波浪碰到边界会反射回来，二次经过岛后方，干扰测点数据。各版本对吸收层的指令命名有差异，但原理一致，照着手册设置即可。
5. SWASH 算例跑不动、跑不对：常见问题排查与避坑
5.1 地形行列顺序反了，锥形岛变成山脊
现象：模型跑起来后，水位云图里的岛屿形态沿某个方向被拉长，原本的圆形干区变成了条状，看起来像一条南北向或东西向的山脊。
原因：生成 bottom.txt 时，两层循环的顺序和 SWASH 读取顺序不一致。SWASH 默认先沿 x 方向扫描，再推进 y 方向；如果你写成了先固定 x、内层循环 y，地形就相当于转置了一次。
解决：修改生成脚本，让外层循环遍历 y、内层循环遍历 x。如果不想改代码，也可以重新排列文本文件的行顺序，用 Python 按 y 分组重写一遍。检查方法很简单：用head -5 bottom.txt看前几行坐标，x 变化而 y 不变就是正确的。
5.2 孤立波提前破碎，波前出现寄生振荡
现象：孤立波还没碰到岛，波形就开始不对称，波峰前面出现短周期小波浪，甚至像海浪一样翻卷。
原因：两个可能，一是波高与水深之比取得过大，强非线性导致孤立波本身不稳定；二是初始场只给了水位没给速度，水包在重力作用下分裂。
解决：把 A/h 控制在 0.25 左右，不要轻易超过 0.4。同时确认初始场文件里有一个速度分量，而不是只抬水面。判断方法是运行前单独画一个初始时刻的沿 x 轴水位与流速剖面，如果流速剖面不是平滑的山峰状，说明速度场有问题。
5.3 测站输出长时间为零，找不到波
现象：某个测站的输出文件是一条水平零线，连初始静水位都读不到，或者恒定在一个值不变。
原因：测站坐标落到了干网格上，或者坐标超出计算域范围。还有一个隐蔽原因：测站底高程高于静水位，模型认为这个网格初始是干的，自然不会更新水位。
解决：测站尽量布置在初始水深至少三倍波高的位置。输出时同时申请水位和水深两个变量，如果水位是零但水深为负，说明测站在干区。检查测站坐标的时候，注意单位要和计算域一致，有人把物理实验的缩比尺寸直接填了进去，结果坐标落在几十米之外。
5.4 计算时间不可接受，一个晚上还没跑完
现象：网格加密到 0.1m 以后，模型一个时间步叠代变得极慢，模拟 60 秒物理时间需要几十个小时。
原因：总网格数过大，加上非静压模式本身每步要解压力泊松方程，计算量随网格数超线性增长。频繁的输出也拖慢运行。
解决：先用 1m 网格跑通整个流程，确认边界和初始场没问题后，再全局加密到 0.5m，最后只在岛周围 1.5 倍底半径范围内加密到 0.25m。控制输出频率，不要每步都输出全场水位，只保留测站时间序列和每 1 秒一帧的全场水位就够了。垂直方向先保持单层，不要一上来就开三层。
6. 让锥形岛算例真正可信：收敛性、敏感性分析与进阶技巧
6.1 没有实验数据时怎么做验证
如果手头没有可对比的物模数据，网格收敛性就是最硬的验证手段。分别用 0.5m、0.25m、0.125m 三组网格跑同一工况，提取岛背风侧的最大爬高值，如果三组结果之间的偏差小于 5%，基本可以认为网格已经不是主要误差来源。同时检查质量守恒，SWASH 输出里的总水量漂移应控制在 1% 以内，超过这个量就要回头看边界条件和干湿阈值。
6.2 参数调优顺序
把参数从“能跑”调到“像样”，我习惯按这个顺序动：先加底摩擦，观察最大爬高下降幅度；再调干湿阈值，从 0.001 到 0.01 各跑一组，看爬高过程线是变得更平滑还是被削平；最后调吸收层强度，直到侧边界反射波不再影响岛后测点。每次都只动一个参数，否则出问题了根本分不清是哪个参数导致的。
6.3 从单层到多层，从过程线到动画
非静压模式的垂直分层从 1 层增加到 3 层，对最大爬高的影响通常不明显，但岛后绕射波形的相位会有所改善。如果只是关心爬高极值，单层够用；如果要做波前沿动画和细节对比，再考虑分层。最后把测站水位数据画成动画，观察波前沿在锥面上的推进方式，这一步能帮你判断结果是否物理合理。每次拿到新算例，我都会先画一遍地形、跑一次十几分钟的粗网格测试再决定全量计算，这个习惯帮我省掉了大量返工时间。希望帮到你。
本文还有配套的精品资源，点击获取
