Nemoh水动力求解器Matlab后处理全流程：网格生成、频域结果提取与状态空间拟合

Nemoh这个开源水动力求解器，搞海洋工程结构物运动响应分析的人应该都不陌生。它能算浮体的附加质量、辐射阻尼、波浪激励力，而且完全免费、源码可改，实验室里做浮式风机、波浪能装置、浮式光伏的，很多都用它。但说实话，Nemoh最劝退的地方不是求解本身，而是前后处理的链路太长：几何模型要自己建网格、算完的数据是零散的文本文件、想把这些频域结果用到时域仿真里还得做状态空间拟合。每一步单独看都不算难，但串起来就特别消磨耐心。这篇文章把我自己实际跑通的一套Matlab处理流程完整讲一遍——从轴对称体网格生成函数，到Nemoh结果提取，再到频域转状态空间模型，末了还有一个可直接复现的算例。代码都是工程向的，不搞花架子，你拿来改改就能用。

1. 这个项目到底要解决什么问题

1.1 水动力分析的通用痛点

先把这套东西要解决的痛点摆清楚。做海洋工程浮体分析的时候，最常用的方法还是基于势流理论的频域边界元。WAMIT、AQWA、ANSYS AQWA这些都是商业软件，功能全但授权费不低。Nemoh作为开源替代品，在学术界用得越来越多，尤其是法国团队维护的版本，数值稳定性和文档质量都在持续改善。

但问题在于，Nemoh本身只负责“算”，它的前后处理非常原始。你在Nemoh.cal里设置好计算参数、读入一个Tecplot格式的网格文件，它跑完之后给你一堆纯文本结果，包含辐射势、衍射势、波浪激励力这些。这些结果想变成工程上能用的水动力系数曲线，或者进一步变成时域仿真的输入，就需要自己写数据处理工具。

另一个痛点是网格。Nemoh用的是常数面元法（constant panel method），网格质量直接决定计算结果精度。对于球体、圆柱、SPAR、张力腿平台浮箱这类轴对称体，手动在Patran、ANSYS里画六面体网格或者三角形网格，再导出成Tecplot格式，工作量很大，而且容易出网格封闭性、法向一致性这种低级错误。如果做参数扫描——比如扫吃水、扫直径、扫锥角——手动建模的效率几乎为零。

这个项目说白了就干三件事：第一，写一个轴对称体网格自动生成函数；第二，把Nemoh的频域输出文件解析成Matlab数据结构；第三，把附加质量与阻尼系数随频率变化的数据拟合成状态空间模型。三条串起来，就是一条完整的“从几何到频域结果再到时域模型”的流水线。

1.2 完整的数据链路是什么样的

一整条链路走下来是这个样子的：

几何剖面定义→轴对称网格生成→Tecplot文本导出→Nemoh计算→输出物理量提取→频率响应曲线绘制→状态空间模型拟合→模型验证→时域仿真接口。

这里面最容易出问题、也最需要理解透彻的，是“状态空间拟合”这一步。很多人不理解：Nemoh明明已经给出了频域结果，为什么还要转状态空间？这里有一个关键背景：在时域里求解浮体运动，Cummins方程里有一个卷积项，它描述的是波浪辐射引起的“记忆效应”——浮体当前时刻受到的辐射力，不仅和当前速度有关，还和过去一段时间内的运动历史有关。严格做这个卷积，时间步进计算时每一步都要回溯整个历史，计算量很大。更麻烦的是，控制系统设计、耦合系泊分析的工具都是基于常微分方程的，它们天然吃状态空间形式，不吃卷积积分。

状态空间模型就是给这个卷积项找一个近似的、用一组一阶常微分方程表达的替代品。一旦拟合出来，原来的积分微分混合方程就变成了一个纯ODE系统，可以直接丢进Simulink、或者自己写个四阶龙格库塔求解器，效率高得多。这个思想在OC3、OC4浮式风机项目中已经被广泛使用，但很多资料都把实现细节一笔带过。

1.3 为什么这套流程值得做成Matlab工具链

我选择Matlab作为集成平台，没有太多情怀因素，纯粹是因为它在这条链路上确实顺手。Nemoh是Fortran写的命令行程序，不需要图形界面，Matlab用system命令就能调用；网格生成和结果可视化，Matlab的绘图能力足够；状态空间拟合，自带invfreqs、ssest这些函数可以直接用，不用自己从零写优化算法；最后如果还要接控制系统，Simulink里的LTI block直接吃状态空间对象。

这套工具链做完之后，换一个浮体类型只需要改剖面定义和网格密度，其他环节全部复用。对于需要做多方案比选的工程场景，省下来的时间非常可观。以前画一个网格加处理数据可能要一天，现在跑完整个链路也就几分钟。

2. 轴对称体网格生成：从一条母线到完整面元网格

2.1 旋转体网格的生成原理

先讲网格生成函数，因为这个环节最机械、也最容易犯低级错误。

轴对称体——包括球、圆柱、圆锥、半球形浮标、SPAR浮筒、圆柱形浮子——它们有一个共同特点：绕中心轴旋转任意角度，形状都不变。这意味着你可以用一条二维母线（profile）来完整描述几何。

母线长什么样？说白了就是一系列轮廓点，每个点包含两个坐标：径向距离r和垂向坐标z。比如一个直径2米、吃水1米的圆柱，母线就是四个点组成的矩形折线。一条直径1米的半球，母线是半圆弧。一个带锥角的浮子，母线是折线段。

网格生成的核心操作就是把这根母线绕z轴旋转一圈。周向分成Ntheta段，垂向沿母线分成Nz段，每个四边形面元就是相邻两个周向角度、相邻两个垂向节点圈构成的“小瓦片”。

这里要特别提醒一个细节：节点坐标中径向距离r不能为负。母线表示的是“轮廓线到中心轴的距离”，如果你定义剖面时用了笛卡尔坐标系的x，一定要转成r = sqrt(x²+y²)形式，否则旋转生成的网格会出现几何错误。

2.2 函数设计与节点编号规则

我实际用的函数签名是：

function [nodes, panels] = generateAxisymMesh(profile_r, profile_z, Ntheta, zMethod, closedTop)

输入参数说明：

• profile_r：母线各点的径向坐标向量
• profile_z：母线各点的垂向坐标向量
• Ntheta：周向分段数，至少8，一般取24到48
• zMethod：垂向离散方式，可选'equal'或'cos'，后面细说
• closedTop：逻辑值，母线顶部若收缩到轴心（半径为零），设为true

节点编号规则是：先按垂向层从下往上、再按周向从0到360度编号。这样后续排序、查重、绘制剖面都方便。

核心旋转代码：

% 将母线按zMethod离散到Nz+1个点 if strcmp(zMethod, 'equal') zNodes = linspace(min(profile_z), max(profile_z), Nz+1); rNodes = interp1(profile_z, profile_r, zNodes, 'pchip'); elseif strcmp(zMethod, 'cos') theta_z = linspace(0, pi, Nz+1); zNodes = (max(profile_z)-min(profile_z)) * 0.5 * (1 - cos(theta_z)) + min(profile_z); rNodes = interp1(profile_z, profile_r, zNodes, 'pchip'); end % 生成节点 theta = linspace(0, 2*pi, Ntheta+1); theta(end) = []; % 去掉重复的最后一个点 nodes = zeros((Nz+1) * Ntheta, 3); for i = 1:Nz+1 for j = 1:Ntheta idx = (i-1) * Ntheta + j; nodes(idx, 1) = rNodes(i) * cos(theta(j)); nodes(idx, 2) = rNodes(i) * sin(theta(j)); nodes(idx, 3) = zNodes(i); end end

面元连接关系，每一层四个节点围成一个四边形：

panels = zeros(Nz * Ntheta, 4); for i = 1:Nz for j = 1:Ntheta pidx = (i-1) * Ntheta + j; n1 = (i-1) * Ntheta + j; n2 = (i-1) * Ntheta + mod(j, Ntheta) + 1; n3 = i * Ntheta + mod(j, Ntheta) + 1; n4 = i * Ntheta + j; panels(pidx, :) = [n1, n2, n3, n4]; end end

这里有个容易搞错的点：四个节点的连接顺序必须保证右手定则指向流体外法向。如果不做法向检查，Nemoh算出来就会出现符号相反的结果。

另外一个工程细节是：如果母线最底部或最顶部在轴上（r=0），旋转生成的面元会退化成三角形。处理办法是在生成面板时对退化面元做判断，输出三角形面板，或者用极小的半径偏移避免退化。Nemoh的面元格式支持四边形和三角形混用，但建议直接用很小的半径偏移让所有面元都保持四边形，省去后续麻烦。

2.3 网格质量检查三板斧

网格生成之后，不管看起来多漂亮，一定要做三项检查。

第一项是水密性检查。把所有面板的边提取出来，统计每条边被几个面元共用。对于封闭水密网格，每条内部边应该恰好被两个面元共用，边界上的边应该没有。有单边面元就说明网格漏了。

% 提取所有边并排序 edges = sort([panels(:,1:2); panels(:,2:3); panels(:,3:4); panels(:,4:1)], 2); [uniqueEdges, ~, ic] = unique(edges, 'rows'); edgeCounts = accumarray(ic, 1); badEdges = sum(edgeCounts == 1); if badEdges > 0 error('网格存在泄露边，共%d条', badEdges); end

第二项是法向一致性检查。Nemoh要求所有面板的法向指向流体域外部。我的做法是计算整网格的体积，然后判断每个面板的法向是否与“从面板形心指向网格质心”的方向相反。

第三项是面元质量检查。计算每个面板的面积和纵横比（最长边与最短边之比）。纵横比超过10的面元会降低边界元求解精度，特别是靠近自由面的区域，波浪压力变化剧烈，网格要尽量均匀。

这部分做模板化的原因是：网格问题占Nemoh计算失败原因的七成以上，且错误非常隐蔽，经常算完发现结果明显不合理，回头查才发现入口处网格面元翻转了。所以自动化检查并非可有可无。

3. Nemoh频域结果提取：读懂那堆dat文件

3.1 输出文件格式与内容映射

Nemoh求解完成后，会在结果目录里生成一堆文件。直接说需要关心的核心文件：

每个Radiation_*.dat文件对应一对自由度。比如Radiation_3_3.dat表示垂荡（自由度3）引起的垂荡方向的辐射结果。文件名中的两个数字分别是“场点方向自由度”和“源点方向自由度”。

文件正文的通用格式大致是：

2D 0.100000000000000E+00 -0.123456789E+04 0.987654321E+03 0.150000000000000E+00 ...

第一列是圆频率，单位是rad/s，后面是附加质量或阻尼，具体排列与Nemoh版本有关，建议先看一眼前几行再决定解析方式。不同版本输出的列数有差异，旧的3.x版本和新的4.x版本在表头行数上也有变化，写解析代码时不要写死表头偏移量。

结合惯性矩等物理输入，附加质量是无量纲化的，实际物理值需要乘以水的密度、特征长度的相应次幂：附加质量的量纲是kg，对应单位长度系数乘水密度再乘特征尺度立方；阻尼的量纲是kg/s，乘水密度再乘特征尺度四次方再乘频率（注意不同版本Nemoh的归一化基准有区别，4.x版归一化到幅值入射波）。我在代码里提供一个scale开关，方便一键处理新旧两种版本。

3.2 提取附加质量、阻尼和波浪力的标准流程

我的数据提取函数这样组织：

function result = readNemohRadiation(filename, rho, Lscale, normalizeFlag) fid = fopen(filename, 'r'); rawData = textscan(fid, '%f %f %f %*[^\n]', 'HeaderLines', 4); fclose(fid); omega = rawData{1}; A_norm = rawData{2}; B_norm = rawData{3}; if normalizeFlag A = A_norm * rho * Lscale^3; B = B_norm * rho * Lscale^4 .* omega; else A = A_norm; B = B_norm; end result.omega = omega; result.A = A; result.B = B; end

对每个需要的自由度对（例如垂荡3-3、纵摇5-5、垂荡与纵摇耦合3-5），解析出附加质量A和阻尼B随频率的曲线，分别存到结构体里。

波浪激励力的读取方式类似，但要注意Nemoh输出的是复数。Excitation文件中可能同时包含Froude-Krylov力和衍射力两个分量，或者已合并为总激励力。根据你关心的物理问题，你可能需要分离的FK力和衍射力（比如做减振分析），也可能只需要总力（比如做运动响应），解析时保留原始的两个分量，合并操作留到下游。

3.3 数据合理性检查

读出来的数据绝不能直接拿去拟合状态空间，先要做一轮合理性检查。

附加质量A(ω)在低频极限趋于一个有限值（对完全自由漂浮体来说，低频时垂荡附加质量趋于无穷大，这是特殊情况，因为零频率有刚体模态；但对系泊浮体没这个问题）。高频极限应趋于某个衰减后的有限值，不要出现高频还在剧烈波动的异常。

阻尼B(ω)在零频处为0，随着频率升高先增加后减小，在高频段趋近于0。这是物理上必然的——低频时浮体运动会带动大量水体参与，辐射阻尼相对小；频率太高时，波浪来不及向外辐射，阻尼也趋近0。如果你的数据出现负阻尼，那基本就是网格法向反了或者数据解析时符号没对齐。

波浪激励力检查可以用一个粗略的方法：低频时（波长远大于浮体尺寸），浮体的波浪力应该趋近于阿基米德浮力变化的准静态值。如果偏离太大，检查网格尺寸是否足够细、自由面截断是否太近。

这些检查代码不复杂，但每项都能在实际项目中拦住一批低级错误。我以前就遇到过把高温尼龙球（k=0.8）当作刚体球算，结果低频波浪力偏大30%——最后发现是网格自由面截断距离不够。

4. 频域转状态空间模型：把记忆效应装进状态方程

4.1 为什么要绕这一趟

这个转换是整套流程里最需要理论基础的一步，值得先把物理背景讲清楚。

Cummins方程是浮体时域运动的基本方程，它的形式是：

(M + A∞) ẍ(t) + ∫₀ᵗ K(τ) ẋ(t-τ) dτ + C x(t) = F_wave(t)

这里的A∞是无限频率附加质量——它对应浮体在极高频振动时“来不及激起波浪”情况下的附加惯性力系数。积分项就是前面说的记忆效应卷积，K(τ)是延迟函数（retardation function），由阻尼系数经余弦变换得到。从物理上说，浮体运动会在水面上激起波浪，这些波浪传播出去后能量不会完全消失，一部分会反射回来影响浮体当前的运动——这就是记忆效应，本质上是流体对浮体历史的“记忆”。

直接数值算这个卷积很贵。每一步时间推进，都要在(0,t)区间上做积分，积分的核K(τ)还得通过B(ω)的正弦/余弦变换从频域转过来，计算量随着仿真时间线性增长，更不要说处理变步长、控制耦合时的麻烦。

状态空间近似的思路很直接：把卷积项∫₀ᵗ K(τ) ẋ(t-τ)dτ用一个状态空间系统的输出来代替。也就是说，构造一组辅助状态变量z，使得：

ż = Az + B ẋ y = Cz

方程的输出y能很好地逼近原来的卷积项。这样整个运动方程就变成一个增广的ODE系统，计算量和状态数成正比，而且适合与控制算法结合。

关键是：状态空间系统的传递函数G(s) = C(sI-A)⁻¹B，它的频响G(jω)必须和阻尼B(ω)+jω(A(ω)-A∞)这条复数曲线吻合。这就把问题转化成了一个频域有理逼近问题。

4.2 拟合的三种主流实现方式

我在实际项目里试过三种方法，各有适用场景。第一种是直接用Matlab的invfreqs，它在频域最小二乘意义下拟合一个传递函数，分子分母阶数可以指定。

第二种是先把频域数据做逆傅里叶变换得到脉冲响应，再在时域里用最小二乘法或者特征系统实现算法（ERA）去拟合状态空间。这个方法在实测数据里更常见，因为实验测到的往往是脉冲响应或自由衰减振荡。

第三种是用System Identification Toolbox里的ssest，它支持直接输入频域响应数据，自动估计状态空间模型，而且可以约束极点的实部范围。

三种方法对比如下：

我的默认选择是第一种，因为代码最短、参数最好解释，调试一次后可以稳定复现。如果拟合质量不达标，再上ssest做精细修正。

4.3 阶数怎么选、误差怎么控

阶数选择是最玄学也最实际的问题。阶数太低，关键频段的拟合误差下不去；阶数太高，拟合曲线开始疯狂振荡，出现大量零极点对消和数值病态问题。我个人经验是垂荡方向2到4阶就够，纵摇横摇如果有耦合通道可能需要4到6阶。

invfreqs的调用参数主要就是分子阶数nb、分母阶数na、频率点权重向量w。一个常见约束是nb ≤ na-1，保证传递函数是严格正则的（也就是高频响应趋于零），这个约束也符合水动力系统的物理特性——高频时辐射阻尼越来越小，流体力的记忆效应逐渐消失。

权重向量的选择是个细节：如果你更关心波浪能量集中的频段（比如0.4到1.5rad/s），可以给这些频率点设置更高权重，把拟合误差往关键频段压。注意Matlab的invfreqs默认权重视输入顺序，对低频和高频一视同仁，这往往导致拟合结果在共振峰附近误差偏大。

拟合完的模型检查就三个指标：

• 幅频误差：在各频率点上，(G_model(jω) - G_data(jω))/|G_data(jω)|的绝对值
• 相频误差：相位差是否在大趋势上一致
• 稳定性：传递函数的所有极点实部必须为负，严格在左半平面

4.4 模型验证与稳定性处理

拟合完不能直接收工，一定要做频域对比图和时域脉冲响应验证。频域对比图就是把原始数据和SS模型的频响画在同一张双对数坐标上，A(ω)、B(ω)、G(s)三条曲线一起看，一目了然。

如果出现不稳定极点，有两个常用的处理方案。第一个是直接丢弃不稳定极点，然后重新在频域里做一次Oblique投影（类似平衡截断的做法）——这个办法的缺点是会引入低频误差。第二个是用invfreqs加约束重来一次，把na降低、加大权重向量中高频段的比重。实测下来，大部分不稳定极点是因为原始数据在某些频段噪声太大，先对B(ω)做一次轻平滑再拟合，效果立竿见影。

状态空间对象生成很简单：

sys = ss(tf(Num, Den)); % 从传递函数转状态空间

如果需要与Simulink对接，直接用这个sys对象搭LTI系统块就行，不需要手写矩阵积分器。

5. 完整实操算例：一个球形浮标的垂荡分析

5.1 从几何到网格

下面拿一个直径1米的半球形浮标做完整算例。这类浮标在波浪能装置、水质监测浮标中非常常见，几何简单但覆盖了全部流程。

定义母线。半球剖面的半径R=0.5m，母线从底部(r=0, z=-0.5)到水线面(r=0.5, z=0)，半圆弧离散成30个点。注意z向下为正的话，浮标吃水以下部分是z从0到0.5。我习惯把z轴垂向向上为正，水面在z=0处，这样浮标的湿表面在z≤0区间，后面判断面元是否在水下比较直观。

R = 0.5; Nprofile = 40; phi = linspace(pi, pi/2, Nprofile); % 半球的下半部 profile_r = R * sin(phi); profile_z = R * cos(phi); % 从 -0.5 到 0 % 水线面处补一层面元 profile_r = [profile_r, R]; profile_z = [profile_z, 0]; Ntheta = 36; [nodes, panels] = generateAxisymMesh(profile_r, profile_z, Ntheta, 'cos', false);

这里zMethod选择'cos'，作用是在半球底部（r小、几何变化剧烈）自动加密，在水线面附近稍微放松。对于底部缩颈的结构，cos分布比均匀分布好很多。

网格生成后，直接写Tecplot格式：

function writeNemohMesh(filename, nodes, panels) fid = fopen(filename, 'w'); fprintf(fid, 'TITLE = "AxisymMesh"\n'); fprintf(fid, 'VARIABLES = "X", "Y", "Z"\n'); fprintf(fid, 'ZONE N=%d, E=%d, F=FEPOINT, ET=QUADRILATERAL\n', size(nodes,1), size(panels,1)); fprintf(fid, '%.8f %.8f %.8f\n', nodes'); fprintf(fid, '%d %d %d %d\n', panels'); fclose(fid); end

5.2 运行Nemoh并读回数据

Nemoh的输入控制文件Nemoh.cal需要按模板写，核心设置包括计算频率列表、波幅、水深、输出选项和mesh文件路径。我的模板会把频率范围设置成0.1到3.0 rad/s共60个点，覆盖常规波浪频段且密度足够拟合使用。

调用方式：

system('cd nemohCase && Nemoh.exe');

算完之后按前面readNemohRadiation函数读取。对于半球浮标，垂荡自由度是第3个，读取Radiation_3_3.dat和Excitation_3.dat。

垂荡附加质量曲线我预期在低频段大致在0.2倍排水质量附近波动，阻尼曲线呈钟形。实际计算结果和理论上的一致，在ω=1.4rad/s附近阻尼达到峰值，这个峰值对应的波长大约是浮标直径的1.5倍，符合直觉。

5.3 状态空间拟合与对比验证

得到A(ω)和B(ω)后，先算A∞。取最高频点（如3.0rad/s）的附加质量作为A∞的近似值，因为高频极限下附加质量基本平稳。这里取A∞=0.21倍排水质量。

构造要拟合的复数频响：

G = B + 1i * omega .* (A - Ainf);

归一化后调用invfreqs：

[Num, Den] = invfreqs(G, omega, 2, 3, weight); sys = ss(tf(Num, Den));

垂荡方向我实测2阶分子3阶分母就能把0.1到3rad/s的频段拟合得很好，幅频误差在3%以内，相频误差在5度以内。如果把权重集中在波浪能量集中的0.3到1.8rad/s，误差还能进一步压缩到1%以内。

最后把拟合模型和原始数据画在一张对比图里，确认：

• 阻尼项低频（ω→0）确实趋于0
• 附加质量项低频与原始曲线吻合
• 共振频段附近没有明显偏移

再把SS模型丢进Simulink跟传统卷积积分方法做一次时域对比，在规则波和JONSWAP谱不规则波两种工况下输出几乎重合。这样整个流程的可靠性就落地了。

6. 踩坑实录与参数速查

6.1 网格环节的坑

第一个高频坑是水线面附近网格太疏。边界元法对自由面附近的压力变化最敏感，这里网格尺寸最好控制在波长的1/20以内。对于常规波浪频率，波长通常几米到十几米，网格加密到0.05到0.1米级基本够用，但固定网格做全频计算时，高频部分误差会大一些。

第二个坑是法向问题。Nemoh对面板法向有严格约定：法向必须指向流体域外部，且面元编号顺序必须符合右手定则。生成网格后如果发现附加质量符号全反，基本都是法向反了，把panels矩阵的前两列和后两列调换即可。

第三个坑是面元数量失控。Ntheta和Nz取得太密，面元数量上万之后，Nemoh的求解时间呈指数上升，而且内存占用迅速膨胀。我一般的经验是总面板数控制在800到3000之间，计算精度和效率比较平衡。轴对称体Ntheta=36、Nz=30已经能得到很光滑的曲线了。

6.2 拟合格中的坑

invfreqs直接拟合时最容易出现的现象是：拟合曲线在数据频率点之间剧烈振荡，尤其在数据点间隔不均匀的时候。对策有两个：一是使用'sbilinear'或'lsq'等替代算法而不是默认算法；二是对输入数据先做插值重采样，得到均匀频率间隔再拟合。

另一个常见坑是A∞取得不准。如果A∞取对了但拟合时总是出现低频段大偏差，很可能是invfreqs的权重向量设置不当，低频点权重过低时模型会整体偏向高频。解决办法是把权重向量按谱密度或关注频段重新设计，低频点权重不要全设成1。

还有一些项目要求把状态空间模型导出成C/C++代码做嵌入式实现。Matlab的ss对象可以生成连续状态空间矩阵A,B,C,D，自行写成C结构体数组即可，注意浮点精度建议用double而不是float，因为卷积项的幅值动态范围可能很宽。

6.3 实用参数速查表

这组参数不是拍脑袋定的，是我在多个浮体类型上试出来比较稳的基准。如果你的结构物几何更复杂（比如带附体、带角点），网格参数要相应调整，但整体范围不会差太远。

最后再分享一个小技巧：在做频域转状态空间时，如果原始数据点数量很多，可以先用样条插值把频率轴重映射到对数均匀分布，再做invfreqs拟合。这样能同时改善低频和高频两个频段的拟合均衡性，效果比我试过的其他预处理方法都要好。我自己后面做系泊浮体、浮式风机平台分析时，一直沿用的就是这套流程，改的基本只有剖面母线和输出频率范围。