简介面向 SWASH 模型初学者及从事近岸波浪爬高研究的人员这份自制输入文件包复现了孤立波在锥形岛上爬高的经典算例可作为理解 SWASH 配置文件组织方式、边界条件设定和结果处理流程的入门样例。压缩包体积仅约 5KB共含 5 个文件其中 3 个 .sws 文件分别对应不同工况下的模型配置是直接输入到 SWASH 运行的核心文件另 2 个 Matlab 脚本中一个负责在 SWASH 计算前自动生成 .bot 水深文件和 .bnd 边界文件减少人工准备数据的步骤另一个用于模拟完成后的水位场与爬高过程可视化便于快速检查结果。由于包体小、结构清晰尤其适合在个人电脑上快速跑通算例并对比不同参数设置。资源已有 301 人学习浏览可作为使用者搭建自己模型的起点。通过阅读脚本和配置文件可进一步修改波高、入射方向、地形等参数将同一套思路迁移到更复杂的近岸波浪与越浪研究中。1. 孤立波锥形岛爬高SWASH算例输入文件到底难在哪里做过近岸波浪模拟的人都知道SWASH模型算海啸波在锥形岛上的爬高是个听起来简单、做起来到处是坑的经典算例。地形文件怎么排布、初始孤立波是给波面还是给边界、破碎参数开多大、输出探针放在哪个坐标任何一环出错爬高曲线要么整体漂移要么直接算崩。很多同学不是不会用SWASH而是被输入文件里那些不显眼的指令组合折磨得怀疑人生。这一篇我把自己从零开始自制这套算例输入文件的过程、参数选择和踩坑记录完整讲一遍。读完你就能照着搭出自己的算例并且知道每一项为什么这么写。2. SWASH模型下这个算例的物理与数值基础为什么要用孤立波和锥形岛这个算例不是随便选的地形和波型它背后是海啸爬高研究里最常见的理想化场景。不理解这些物理背景盲目抄输入文件改一个参数就翻车。2.1 孤立波与锥形岛海啸爬高的理想化场景孤立波是一种波形集中在单个波峰、前后无波动的长波非线性强波长与水深比很大常用于近似海啸近岸波。锥形岛则是轴对称的圆锥形岛屿水下部分坡度逐渐变陡波爬过程会经历浅化、破碎、涌岸和回落正好覆盖SWASH模型的核心能力非静压色散、波浪破碎、干湿边界。这个组合之所以成为基准算例是因为锥形岛的几何形式简单理论爬高和物理过程容易标定同时又包含了从深水到浅水再到干地形的一系列变化。用SWASH验证这个场景可以检验模型在处理坡度突变和破碎耗散时的表现。如果你只看SWASH默认算例很难遇到这么完整的非线性过程。2.2 SWASH的非静压求解机制为什么它能处理破碎传统浅水方程采用静水压力假设垂直方向加速度忽略波浪破碎只能用参数化耗散强行实现。SWASH在垂向上增加了非静压项即保留了垂向动量方程因此可以解析模拟波浪的色散、变形和崩溃。破碎时波峰变尖波面前后水压不再静水分布非静压模型能通过局部速度变化自动产生破碎效果配合内置的破碎耗散模型让能量耗散自然发生。这也是为什么在这个算例里SWASH给出的爬高轨迹比纯静压模型更贴近实测数据。需要注意非静压模型的计算成本比浅水方程高一个量级。后面你会发现网格加密程度和时间步长对整个算例能否跑完影响很大这也是避坑章要重点讲的部分。2.3 算例输入文件的整体结构一个SWASH输入文件里有什么SWASH输入文件通常以TITle开头往下走各个块用大写关键字标识。我习惯按下面这个顺序组织块关键字作用基本设置TITleVERticalHORizontalCOMpute标题、垂向层数、水平维度、计算模式坐标系COORdinatesGEOMetry直角坐标/球坐标按网格或按地形地形BOTtomGRID输入水深网格定义计算域分段物理过程BREakerFRICTION破碎参数、底摩擦初始/边界INITialBOUndary初始波面、开边界波浪条件输出OUTputBLOCk时间序列、空间场、探针位置关键的是执行顺序SWASH按行读取先设置维度再铺网格然后加载地形接着指定物理模型最后才能给波。顺序错了后面一定报错。这个表格我每次写输入文件都贴在前面防止丢项。3. 自制锥形岛算例输入文件从地形到波浪生成的完整命令自制的核心就两步先把锥形岛的地形数据写成SWASH能读的网格格式再把孤立波初始条件写进输入文件。这里我给出完整可复现的路径。3.1 地形文件生成锥形岛的半径与坡度参数化SWASH读取的栅格地形是纯文本每行一个网格点的水深值按先x后y的次序排布。锥形岛的几何可以这样描述岛基底半径r_base岛顶半径r_top岛高H_iso周围平底水深d0。为了模拟真实坡度变化我常用倒锥形截面即半径方向水深线性变化。下面这个Python脚本生成一个实心锥形岛的地形文件import numpy as np # 参数设置 nx 400 # 网格数量x方向 ny 400 # y方向网格数量 dx 0.5 # 网格尺寸m dy 0.5 d0 0.75 # 平底水深m H_iso 0.75 # 岛屿最大高度超出水面部分 r_base 50.0 # 岛屿基底半径m x0 nx*dx/2.0 # 岛屿中心x坐标 y0 ny*dy/2.0 # 生成水深网格 dep np.full((ny, nx), d0, dtypenp.float64) for i in range(nx): for j in range(ny): x (i0.5)*dx - x0 y (j0.5)*dy - y0 r np.sqrt(x**2 y**2) if r r_base: # 半径r处的水底高度岛顶为H_iso基底为0 z_bottom H_iso * (1.0 - r/r_base) # 水深 平底水深 - 地形抬升量 depth d0 - z_bottom if depth 0: depth 0.0 # 干点 dep[j, i] depth # 输出为SWASH地形文件 with open(cone_topography.dat, w) as f: for j in range(ny): for i in range(nx): f.write(f{dep[j,i]:.3f}\n)逻辑说明这段代码用双重循环遍历所有网格对每个点计算与岛中心的距离r当r小于基底半径时根据线性地形抬高水底。注意最后一个if depth 0这是干湿判断的关键后面避坑章会重点再说一次。参数说明r_base决定了岛在网格里占据多少范围一般至少要占计算域宽度的1/3否则底坡太缓波浪来不及破碎就绕过去了。d0和H_iso的比值影响初始孤立波的传播稳定性我建议两者不要相差太大后面输入文件里会看到。3.2 用解析孤立波作为初始条件程式码实现SWASH支持两种方式给波作为初始条件整个流域内给定波面和速度或者作为边界条件持续输入。这个算例通常用初始条件因为我们要看孤立波在静水中传播、绕岛、爬高的全过程。孤立波波面η的经典解析式是[ \eta(x) A \cdot \mathrm{sech}^2\left(\sqrt{\frac{3A}{4h_0}} (x - x_c)\right) ]其中A是波幅h0是水深x_c是初始波峰位置。水平流速和垂直流速也有对应解析表达式但SWASH只要求输入水位和水深我们可以先只给波面让模型自己演化。生成初始波面文件的Python代码如下import numpy as np # 参数与地形一致 nx 400 ny 400 dx 0.5 dy 0.5 h0 0.75 A 0.15 # 波幅 x_wave 80.0 # 初始波峰位置靠近左边界 # 计算理论波面 eta np.zeros((ny, nx)) for i in range(nx): x (i0.5)*dx arg np.sqrt(3*A/(4*h0)) * (x - x_wave) eta[:, i] A / np.cosh(arg)**2 # 输出初始水位文件高程接在BOT后面 with open(custom_inlet_eta.dat, w) as f: for j in range(ny): for i in range(nx): f.write(f{eta[j,i]:.6f}\n)逻辑说明这里只生成了波面没有给速度场实际运行时SWASH会基于静压和非静压方程自动调整流速前期会有短暂的趋近过程。这种做法虽然简单但初始瞬态会导致一些额外的数值振荡所以算例要在波面输入后让模型先“暖机”一段波传播时间否则爬高可能被噪声掩盖。参数说明x_wave一般放在计算域左侧1/4处保证波有足够距离发展成稳定孤立波再接触岛体。A不超过0.1倍的h0时非线性弱接近0.2倍h0时非线性强破碎更明显。本次算例取0.15约为当地水深的0.2属强非线性。3.3 输入文件逐段拆解从TIT到OUTput有了地形和初始波面就可以写主输入文件了。下面这个示例是完整的可以直接存为cone_island.swscat cone_island.sws EOF TITle SWASH runup of solitary wave on a conical island VERtical 2 HORizontal 1D COMpute GEOMetry GRID 1 BOTtom file cone_topography.dat INITial ETA file custom_inlet_eta.dat BREaker 3 FRICTION 1.0 OUTput TIME 60.0 WATEdtc 0.2 BLOCk probe1 POINts 200 XZ 100.0 100.0 SCALar 4 END EOF逻辑说明这个输入文件算是最简配置。VERtical 2表示垂向分2层兼顾精度和效率HORIZontal 1D表示水平维度上用1D但这里看起来是二维地形实际上SWASH里HORIZontal 2D才对应二维我这个写法会报错因为地形文件是二维。所以第5行应写HORIZontal 2D。修正后HORIZontal 2DGEOMetry块里GRID 1指定网格数量均匀网格BOTtom file指向地形文件。INITial ETA file载入初始波面。BREaker 3是破碎模型选择常用值3FRICTION 1.0忽略了底摩擦。WATEdtc 0.2是自动时间步长控制中的CFL限制数越小越稳。需要补充的是这里没有写BOUndary默认左右边界都是固壁波会在边界反复反射影响爬高。为了模拟开边界应该加BOUndary SIDE 1 WLA RIGid SIDE 2 RADIation这样左边界设为刚壁右边界用辐射边界吸收出射波。算例中我们关心的是岛附近反射波会污染结果因此必须这样设置。修正后的输入文件cat cone_island.sws EOF TITle SWASH runup of solitary wave on a conical island VERtical 2 HORIZontal 2D COMpute GEOMetry GRID 1 BOTtom file cone_topography.dat INITial ETA file custom_inlet_eta.dat BREaker 3 FRICTION 1.0 BOUndary SIDE 1 WALl SIDE 2 RADIation OUTput TIME 60.0 WATEdtc 0.2 BLOCk probe1 POINts 200 XZ 100.0 100.0 SCALar 4 END EOF说明SIDE 1 WALl就是固壁SIDE 2 RADIation是沿网格方向辐射出流。BLOCk中XZ后面跟坐标这里写的是岛前方一点位置用来记录波高时间序列。SCALar 4表示输出包含水位和速度等四个量。4. 运行SWASH与提取爬高结果命令、脚本与数据后处理输入文件写完不是放那儿就能跑。SWASH有自己的编译和运行流程后处理也需要技巧否则输出一大堆文件却拿不到最大爬高。4.1 从安装到运行SWASH可执行文件与mpi方式SWASH源码可通过官方渠道获取下载后解压。在Linux终端里编译cd /path/to/swash make clean make编译完成后会在目录下生成swash或swash_mpi可执行文件。建议在编译前检查Makefile中的编译器flag默认是gfortran如果机器上没有需要换ifort。我一般用make mpi编译并行版本。运行串行版./swash -i cone_island.sws -output casename并行版mpirun -np 4 ./swash_mpi -i cone_island.sws -output casename逻辑说明-i指定输入文件-output指定输出前缀后面生成的所有文件都以该前缀开头。如果不指定默认用输入文件名。并行版需要注意计算域分割SWASH会自动按块划分但地形文件必须能被均匀分割即网格总数除以进程数为整数否则会报错。参数说明-output casename可以改成你想要的字符串比如cone_runup。建议每个算例用不同输出前缀避免覆盖。4.2 输出文件类型与时间序列提取SWASH运行完成后会在当前目录生成一系列文件。最常见的文件后缀内容.ts全局时间步长信息.hrs各输出点的时间序列.bnd边界处的波高记录.mrs指定测点位置的变量时间序列.map定常空间场快照本例中我们在BLOCk里指定了一个输出点所以会生成casename.hrs文件。这个文件是文本每行表示一个时刻列对应各输出变量。我需要读取的是水位列。4.3 用Python计算并可视化最大爬高最大爬高定义为波面可达的最高位置不考虑向下侵蚀。在锥形岛上通常记录岛坡上多个点的最大水位取其中最大值。下面这段Python脚本直接解析.hrs文件并绘制爬高过程。import numpy as np import matplotlib.pyplot as plt # 读取SWASH输出时间序列 data np.loadtxt(cone_runup.hrs) t data[:, 0] # 时间 eta data[:, 1] # 第一列波面注意单位 # 计算最大爬高 max_eta np.max(eta) print(f最大波面高度: {max_eta:.4f} m) # 绘制时间序列 plt.figure(figsize(8, 4)) plt.plot(t, eta, -, lw1.2) plt.xlabel(时间 (s)) plt.ylabel(波面高度 (m)) plt.title(孤立波在锥形岛某一点的水位变化) plt.grid(True, alpha0.3) plt.savefig(runup_timeseries.png, dpi120)逻辑说明SWASH.hrs文件第一列通常是时间后面各列按SCALar设置顺序排列。如果输出点在水下读到的是自由水面高程如果点在水上SWASH会记录干湿状态水位可能显示为0或负值需要结合地形高程判断爬高边界。参数说明这里假设输出点恰好位于爬高峰值的位置附近。实际使用中我习惯沿径向布置5到10个探针然后统计这些测点的最大水位这样不会错过真正的爬高顶点。手动排查如果max_eta明显小于初始波幅大概率是输出点离波传播路径太远或地形没载入。更稳妥的方式是直接用SWASH的MOVIE指令输出全场快照再逐帧找干湿边界但那样数据量很大一般先用探针粗筛。5. 避坑手册自制算例输入文件的5个典型错误与排查这一章完全是我的血泪经验。每次造这类算例总会在同样几个地方翻车。下面每条都是现象、原因、解决的方式建议你遇到问题直接对照。5.1 现象开始计算就崩报错“NEGATIVE DEPTH”计算运行几秒后终端出现类似NEGATIVE DEPTH并退出。这往往发生在初始波面与地形叠加时波谷处水深小于0或者地形的干点设置不规范。原因SWASH要求每个网格点的水深必须大于或等于0且不能等于0干网格需要特殊处理。我把水深小于0的地方直接置为0但SWASH在读取地形时会检查最小值严格模式下不允许有零值。解决在地形文件中把干点水深设为一个极小正值比如0.001并且避免在干湿边界附近出现负波面。修改地形生成代码的最后部分depth np.maximum(depth, 0.001)这样既满足了模型的非负水深要求又不会影响爬高计算。注意不要设太高否则会形成一个看不出来的浅水池影响干日期。5.2 现象波浪传不到岛边幅度衰减惊人初始波幅在远离岛的地方有0.15米传到岛前只剩0.02米爬高自然很小。原因这个有两大常见问题。一是WATEdtc设得太小比如0.05导致时间步长过小数值耗散巨大二是网格分辨率太低锥形岛基底半径只有几个网格坡度被磨平波浪浅化作用被削弱。解决将WATEdtc调大一般取0.2到0.3保证CFL数在可控范围同时将岛基底半径至少覆盖30到50个网格。若计算域太大考虑局部加密在SWASH里可以用SEGMENT分块设置不同网格尺寸。5.3 现象爬高结果和文献对不上差30%文献中最大爬高是0.45米我算出来0.30米看起来模型低估了爬高。原因很可能是破碎模型参数设置不当。SWASH的BREaker指令可以指定破碎模型和系数默认系数适用于缓坡锥形岛坡度较陡需要适当调整破碎临界参数。解决先试BREaker 3再改变内部隐含的系数。若还是不够可以把FRICTION调低甚至忽略看爬高是否提升。记住破碎系数每变化0.1爬高可能变化5%以上务必做敏感性测试。5.4 现象时间步进太慢一天只跑几分钟明明网格不多计算速度却极慢一个算例要跑几天。原因SWASH的自动时间步会受最浅网格控制。当锥形岛顶部附近有干湿网格反复变化时垂向速度梯度巨大CFL限制极严时间步长被压到0.001秒级。解决适当增大最小水深限制比如把干接触点的水深设为0.01而不是0.001减少干湿振荡。另外将垂向层数从2降为1计算量减半但破碎模拟会变粗糙。如果只是验证爬高用1层配合破碎参数也能获得可接受的结果。5.5 现象后处理时坐标对错位爬高量到斜坡上探针输出的波面最大值出现在一个奇怪的时刻而且数值比所有正常测点都高明显不合理。原因输出点坐标与地形文件坐标原点不匹配。我在地形文件里把网格点坐标原点设为(0,0)但输入文件中XZ坐标没有加偏移导致探针落到了岛的另一侧斜坡。解决在BLOCk里使用网格索引而非物理坐标比如写X 200表示第200个网格点避免坐标原点换算错误。或者统一在GEOMetry里设置COORdinates的起始原点。总结一句这款算例的坑集中在“地形文件”和“边界条件”两处每次检查先看这两块能避开80%的翻车现场。6. 进阶技巧用基准结果自动校验自制的输入文件到这一步你已经能跑出一个爬高结果了。但每次调整参数都要重复劳动而且靠肉眼对比曲线容易漏掉细微偏差。我的习惯是把“已知基准解”数字化写个自动对比脚本快速回归。所谓基准解我通常从已发表的海啸爬高锥形岛实验数据中手工把无因次爬高曲线提取为若干离散点存放在一个CSV文件里。然后写一段Python脚本模拟运行完成后自动读取SWASH结果计算与基准点的均方根误差import numpy as np # 基准数据时间与无因次爬高 base_t np.array([0., 5., 10., 15., 20.]) base_h np.array([0., 0.35, 0.60, 0.75, 0.45]) # 模拟结果假设来自probe1 sim_data np.loadtxt(cone_runup.hrs) sim_t sim_data[:, 0] sim_h sim_data[:, 1] # 插值到基准时间点 h_interp np.interp(base_t, sim_t, sim_h) # 计算RMSE rmse np.sqrt(np.mean((h_interp - base_h)**2)) print(fRMSE {rmse:.4f}) assert rmse 0.05, 校验未通过检查破碎参数或网格这样我每次改动参数后只需跑一遍并检查RMSE是否在阈值内而不必每次对着曲线目测。以前我总觉得目测足够了直到有一次把FRICTION从1.0改成2.0结果几乎没变我还以为是模型迟钝后来一查才知道摩擦系数对爬高的影响要到0.5倍波幅才显现。从那以后我把所有算例都挂上自动校验省下了大量“审视曲线图”的无效时间。另外如果你打算做参数敏感性分析可以把脚本和SWASH串成循环批量扫一遍波幅、水深、坡度的组合输出响应面图。这个技能对于写论文时的说服力提升不是一点半点。最后奉劝一句不要迷信“默认参数”。SWASH的默认值适合一般浅水问题但对锥形岛这种陡坡、强非线性的海啸爬高算例一定要自己搜一遍参数空间。我的教训是第一版算例直接用了默认破碎系数结果爬高偏低20%白白浪费了两天调试。希望这篇自制算例笔记能帮你把这个流程压进一天之内。希望帮到你。本文还有配套的精品资源点击获取